function [res_bestIndv, res_bestScore] = fitEP_Rap1_4param(n, numGeneration)
    % Rap1モデルのパラメータをEPで推定するモデル。
    % 引数のnは系列数。numGenerationは世代数。

    res_bestIndv = zeros(n, 4); % 結果(パラメータ)を格納する変数
    res_bestScore = zeros(n, 1); % 結果(RSS)を格納する変数
    
    for i = 1:n
        % 各系列で最小のRSSとなるパラメータを探し出し、変数に格納
        [res_bestIndv(i, :), res_bestScore(i, :)] = work_EP_Rap1(numGeneration);
    end
    save('result_fitEP_Rap1_4para.mat') % 結果をmatファイルとして保存    
end

function [bestIndv, bestScore, report] = work_EP_Rap1(numGeneration)
    
    y0 = [1, 0, 0.005, 0]; % [S, GEFa, GAPa, Rap1a]
    param = zeros(1, 4);
    load('data_Rap1_4para.mat'); % x: S, GEFa, GAPa, Rap1aのデータ、t:時間
    
    %%
    %%%%  EPパートここから  %%%%
    currentTime = clock; % 現在の時刻を取得
    second = currentTime(6); % 現在の時刻から’秒’を取得
    
    strSeed = '2018';
    seed = str2num(strSeed)*second*1000; % 乱数の種に‘秒’をかけることで、よりランダムに
    try RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', seed)); % 乱数の初期設定
    catch; RandStream.setDefaultStream(RandStream('mt19937ar', 'Seed', seed)); % ver. 2012以前
    end
    
    lb = 1*10^-3 * ones(size(param)); % パラメータ探索範囲の下限
    ub = 1*10^2 * ones(size(param)); % パラメータ探索範囲の上限
    
    numParents = 20; % 子を生成する親の数
    
    [bestIndv, bestScore, report] = copasiep(@(param) cal_RSS(x, t, param, y0), ...
        numParents, numGeneration, lb, ub); % copasiのEPでパラメータサーチ
    %bestIndv: パラメータの最適解
    %bestScore: RSSの最小値
    %%%%  EPパートここまで  %%%%
    
    %% 各系列で最小のRSSのときのパラメータを用いてシミュレーション
    time = [0, 100];
    [t_est, x_est] = ode15s(@(t, y) ODE(t, y, bestIndv), time, y0);

    % データ点とシミュレーション結果の重ね描き
    figure(1);
    for i = 1:2
        
        subplot(1, 2, i);
        hold on
        plot(t, x(:, i*2), 'bo'); % データ点
        hold on
        plot(　　　　　　　　　　　　, 'r:'); % シミュレーション結果
        
        if i == 1
            title('GEFa');
        else
            title('Rap1');
        end
    end

    % EPでRSSが減少して行く様子を描画
    figure(2);
    subplot(1, 2, 1);
    hold on
    plot(report(:, 1), 'b-');
    title('最も良い個体');
    xlabel('世代数');
    ylabel('残差2乗和');
    
    subplot(1, 2, 2);
    hold on
    plot(report(:, 2), 'b:');
    title('全個体の平均');
    xlabel('世代数');
    ylabel('残差2乗和');
    
end

function [RSS] = cal_RSS(x_data, SamplingTime, param, y0)
    % 残差2乗和を計算
    [~, tc] = ode15s(@(t, y) ODE(t, y, param), SamplingTime, y0);
    residual = 　　　　　　　　　　; % 残差
    RSS = 　　　　　　　　　　　　; % 残差2乗和
end


function dydt = ODE(t, y, param)
    % ここにモデルをコピー
    　　　　　　　　　　　　　　　　　
    
end

