clear
addpath(genpath('../../cellLib/Scripts/'))
%% Build new components from the .csv
csvFile=[pwd filesep 'Vars.csv'];
libPath="C:/Users/wai484/Documents/ECC_MSK_Rios_1993/cellLib";
%% Reuse existing cellml components
cellmlFiles=[libPath+ "/Protocols/Patch_clamp_protocol.cellml",libPath+ "/Components/time.cellml",...
    libPath+ "/Components/MarkovS.cellml",libPath+ "/Components/MarkovS.cellml",...
    ];
whichComp={["dPulse_protocol_ms"],["time_ms"],["S2"],["S3"]};
% Specify the inputs of the existing components, otherwise the variable would
% be categorized as parameters if pub:in
labelRead={};
nm=1;
labelRead(nm).comp="S2";
labelRead(nm).vars=["S_1","S_2","norm","t","k_i1","k_i2","k_1i","k_2i","Si_init"];
nm=nm+1;
labelRead(nm).comp="S3";
labelRead(nm).vars=["S_1","S_2","S_3","norm","t","k_i1","k_i2","k_i3","k_1i","k_2i","k_3i","Si_init"];
%% Encapsuate new components
encap={};
m=1;
encap(m).name="C0_S2";
encap(m).chd.name=["S2",];
encap(m).chd.def=["S2"];
encap(m).chd.label=[];
m=m+1;
encap(m).name="CN_S2";
encap(m).chd.name=["S2",];
encap(m).chd.def=["S2"];
encap(m).chd.label=[];
m=m+1;
encap(m).name="C_S3";
encap(m).chd.name=["S3",];
encap(m).chd.def=["S3"];
encap(m).chd.label=[];
m=m+1;
encap(m).name="O0_S2";
encap(m).chd.name=["S2",];
encap(m).chd.def=["S2"];
encap(m).chd.label=[];
m=m+1;
encap(m).name="ON_S2";
encap(m).chd.name=["S2",];
encap(m).chd.def=["S2"];
encap(m).chd.label=[];
m=m+1;
encap(m).name="O_S3";
encap(m).chd.name=["S3",];
encap(m).chd.def=["S3"];
encap(m).chd.label=[];

m=m+1;
encap(m).name="MWC_10";
encap(m).chd.name=["C0","C1","C2","C3","C4","O0","O1","O2","O3","O4"];
encap(m).chd.def=["C0_S2","C_S3","C_S3","C_S3","CN_S2","O0_S2","O_S3","O_S3","O_S3","ON_S2"];
label={};
N=10;
for nm=1:N/2
    i=nm-1;
    label(nm).comp=sprintf("C%d",i);
    if i==0 || i==N/2-1
        label(nm).vars=["Si_init",sprintf("C%d_init",i);"S_i",sprintf("C_%d",i);"dS_i",sprintf("dC_%d",i);];
    else
        label(nm).vars=["Si_init",sprintf("C%d_init",i);"S_i",sprintf("C_%d",i);"dS_i",sprintf("dC_%d",i);"i",sprintf("C%d_i",i);];
    end
end
for nm=N/2+1:N
    i=nm-(N/2+1);
    label(nm).comp=sprintf("O%d",i);
    if i==0 || i==N/2-1
        label(nm).vars=["Si_init",sprintf("O%d_init",i);"S_i",sprintf("O_%d",i);"dS_i",sprintf("dO_%d",i);];
    else
        label(nm).vars=["Si_init",sprintf("O%d_init",i);"S_i",sprintf("O_%d",i);"dS_i",sprintf("dO_%d",i);"i",sprintf("O%d_i",i);];
    end
    
end

encap(m).chd.label=label;
m=m+1;
encap(m).name="MWC_18";
encap(m).chd.name=["C0","C1","C2","C3","C4","C5","C6","C7","C8","O0","O1","O2","O3","O4","O5","O6","O7","O8"];
encap(m).chd.def=["C0_S2","C_S3","C_S3","C_S3","C_S3","C_S3","C_S3","C_S3","CN_S2","O0_S2","O_S3","O_S3","O_S3","O_S3","O_S3","O_S3","O_S3","ON_S2"];
label={};
N=18;
for nm=1:N/2
    i=nm-1;
    label(nm).comp=sprintf("C%d",i);
    if i==0 || i==N/2-1
        label(nm).vars=["Si_init",sprintf("C%d_init",i);"S_i",sprintf("C_%d",i);"dS_i",sprintf("dC_%d",i);];
    else
        label(nm).vars=["Si_init",sprintf("C%d_init",i);"S_i",sprintf("C_%d",i);"dS_i",sprintf("dC_%d",i);"i",sprintf("C%d_i",i);];
    end
end
for nm=N/2+1:N
    i=nm-(N/2+1);
    label(nm).comp=sprintf("O%d",i);
    if i==0 || i==N/2-1
        label(nm).vars=["Si_init",sprintf("O%d_init",i);"S_i",sprintf("O_%d",i);"dS_i",sprintf("dO_%d",i);];
    else
        label(nm).vars=["Si_init",sprintf("O%d_init",i);"S_i",sprintf("O_%d",i);"dS_i",sprintf("dO_%d",i);"i",sprintf("O%d_i",i);];
    end
    
end
encap(m).chd.label=label;
%% Add equations to components
Eqs={};
m=1;
Eqs(m).comp="C0_S2";
Eqs(m).listEqs=["k_i1 = k_L;",...
    "k_i2 = N*k_C;",...
    "k_1i = k_Lminus;",...
    "k_2i = k_Cminus;"];
m=m+1;
Eqs(m).comp="CN_S2";
Eqs(m).listEqs=[" k_i1 = k_Cminus*N;",...
       " k_i2 = k_L/pow(f, N);",...
       " k_1i = k_C;",...
       " k_2i = k_Lminus*pow(f, N);",...
    ];
m=m+1;
Eqs(m).comp="C_S3";
Eqs(m).listEqs=["k_i1 = k_Cminus*i;",...
    "k_i2 = k_L/pow(f, i);",...
    "k_i3 = (N-i)*k_C;",...
    "k_1i = (N-i+1{dimensionless})*k_C;",...
    "k_2i = k_Lminus*pow(f, i);",...
    "k_3i = (i+1{dimensionless})*k_Cminus;",...
    ];
m=m+1;
Eqs(m).comp="O0_S2";
Eqs(m).listEqs=["k_i1 = k_Lminus;",...
    "k_i2 = N*k_C/f;",...
    "k_1i = k_L;",...
    "k_2i = f*k_Cminus;",...
    ];
m=m+1;
Eqs(m).comp="ON_S2";
Eqs(m).listEqs=["k_i1 = k_Cminus*N*f;",...
        "k_i2 = k_Lminus*pow(f, N);",...
        "k_1i = k_C/f;",...
        "k_2i = k_L/pow(f, N);",...
    ];
m=m+1;
Eqs(m).comp="O_S3";
Eqs(m).listEqs=["k_i1 = k_Cminus*i*f;",...
    "k_i2 = k_Lminus*pow(f, i);",...
    "k_i3 = (N-i)*k_C/f;",...
    "k_1i = (N-i+1{dimensionless})*k_C/f;",...
    "k_2i = k_L/pow(f, i);",...
    "k_3i = (i+1{dimensionless})*k_Cminus*f;",...
    ];
m=m+1;
Eqs(m).comp="MWC_10";
Eqs(m).listEqs=["k_C = 0.5{dimensionless}*alpha*exp((V-V0)/(8{dimensionless}*K));",...
    "k_Cminus = 0.5{dimensionless}*alpha*exp(-(V-V0)/(8{dimensionless}*K));",...
    "O = O_0+O_1+O_2+O_3+O_4;",...
    "C = C_0+C_1+C_2+C_3+C_4;",...
    "norm = C+O;",...
    "L=k_Lminus/k_L;",...
    "P_o = O/norm;",...
    "Q = (O_1+2{dimensionless}*O_2+3{dimensionless}*O_3+4{dimensionless}*O_4+C_1+2{dimensionless}*C_2+3{dimensionless}*C_3+4{dimensionless}*C_4)/4{dimensionless};",...
    "dQ = (dO_1+2{dimensionless}*dO_2+3{dimensionless}*dO_3+4{dimensionless}*dO_4+dC_1+2{dimensionless}*dC_2+3{dimensionless}*dC_3+4{dimensionless}*dC_4)/4{dimensionless};",...
    "P_ss=pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))*pow(f,-2{dimensionless})),4{dimensionless})/(pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))*pow(f,-2{dimensionless})),4{dimensionless})+L*pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))),4{dimensionless}));",...
    "Q_ss=(pow(f,-2{dimensionless})*exp((V-V0)/(4{dimensionless}*K))*pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))*pow(f,-2{dimensionless})),3{dimensionless})+L*exp((V-V0)/(4{dimensionless}*K))*pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))),3{dimensionless}))/(pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))*pow(f,-2{dimensionless})),4{dimensionless})+L*pow((1{dimensionless}+exp((V-V0)/(4{dimensionless}*K))),4{dimensionless}));",...
    ];
m=m+1;
Eqs(m).comp="MWC_18";
Eqs(m).listEqs=["k_C = 0.5{dimensionless}*alpha*exp((V-V0)/(8{dimensionless}*K));",...
    "k_Cminus = 0.5{dimensionless}*alpha*exp(-(V-V0)/(8{dimensionless}*K));",...
    "O = O_0+O_1+O_2+O_3+O_4+O_5+O_6+O_7+O_8;",...
    "C = C_0+C_1+C_2+C_3+C_4+C_5+C_6+C_7+C_8;",...
    "norm = C+O;",...
    "P_o = O/norm;",...
    "Q = (O_1+2{dimensionless}*O_2+3{dimensionless}*O_3+4{dimensionless}*O_4+5{dimensionless}*O_5+6{dimensionless}*O_6+7{dimensionless}*O_7+8{dimensionless}*O_8+C_1+2{dimensionless}*C_2+3{dimensionless}*C_3+4{dimensionless}*C_4+5{dimensionless}*C_5+6{dimensionless}*C_6+7{dimensionless}*C_7+8{dimensionless}*C_8)/8{dimensionless};",...
    "dQ = (dO_1+2{dimensionless}*dO_2+3{dimensionless}*dO_3+4{dimensionless}*dO_4+5{dimensionless}*dO_5+6{dimensionless}*dO_6+7{dimensionless}*dO_7+8{dimensionless}*dO_8+dC_1+2{dimensionless}*dC_2+3{dimensionless}*dC_3+4{dimensionless}*dC_4+5{dimensionless}*dC_5+6{dimensionless}*dC_6+7{dimensionless}*dC_7+8{dimensionless}*dC_8)/8{dimensionless};",...
    ];
%%
[comp,idx]=buildComponents(csvFile,cellmlFiles,whichComp,labelRead,encap,Eqs);
save("comp.mat","comp","idx")