master_genecircuit()

% System Biology
% Direct transcription
% four nodes

function gene = master_genecircuit()
clc;clear
p.nt = 201;
p.t0 = 0;
p.tf = 100;
p.t1 = 40;  % time when input changest
p.ns = 4;
p.tstep = p.tf/(p.nt-1);
p.tspan = 0:p.tstep:p.tf;
p.h = p.tstep;
p.I1 = 0.5;
p.I2 = 0.6;
p.n_param = 64; 
p.xinit = [0;0;0;0];
p.peak = 5; 
p.start = 5;
p.end = 15;

p.pho = 1000;
p.lhs_num = 40;

xd0_lhs = lhsdesign(p.lhs_num,p.n_param)';
xd_lhs = zeros(p.n_param,p.lhs_num);
index_lhs = cell(p.lhs_num,1);
log_precision_lhs = zeros(1,p.lhs_num)';
sensitivity_lhs = zeros(1,p.lhs_num)';

exitflag_lhs = zeros(1,p.lhs_num)';
inform_lhs = zeros(1,p.lhs_num)';
iter_lhs = zeros(1,p.lhs_num)';

tic
parfor j = 1: p.lhs_num
xd0 = xd0_lhs(:,j);
[tout,z] = ode45(@(t,X)system(t,X,xd0,p),p.tspan,p.xinit);
x_state = reshape(z',p.nt*p.ns,1);

x0 = [xd0;x_state]; 

lb_xd = [];
ub_xd = [];
for i = 1: p.n_param/2
    lb_xd = [lb_xd 0 0];
    ub_xd = [ub_xd 10 100];
end

lb = [lb_xd -0.01*ones(1,p.ns*p.nt)];
ub = [ub_xd ones(1,p.ns*p.nt)];
sensitivity_den = (p.I2-p.I1)/p.I1;

A = zeros(6,p.n_param+p.ns*p.nt);
A(1,:) = [zeros(1,p.n_param+p.ns*(p.t1/p.h)+3) sensitivity_den+1 ...
     zeros(1,p.ns*(p.peak/p.h-1)+3) -1 ...
     zeros(1, p.ns*(p.nt-(p.t1+p.peak)/p.h-1))];
A(2,:) = [zeros(1,p.n_param+p.ns*(p.t1/p.h)+3) -1 ...
    zeros(1,p.ns*(2*p.peak/p.h-1)+3) 1 ...
    zeros(1, p.ns*(p.nt-(p.t1+2*p.peak)/p.h-1))];
A(3,p.n_param+1) = 1;
A(4,p.n_param+2) = 1;
A(5,p.n_param+3) = 1;
A(6,p.n_param+4) = 1;

b_U = [0;0.01;0;0;0;0];
b_L = [-inf*ones(1,1); zeros(5,1)];

fprintf('\n Direct Transcription');
tic;
% TOMLAB Syntax 
Name = 'DT';

x_L = lb;
x_U = ub;


c_L = zeros(1,p.nt*p.ns);
c_U = zeros(1,p.nt*p.ns);
   
Prob = conAssign(@(x)obj(x,p),[],[],[],x_L,x_U,Name,x0,[],[], ...
    A,b_L,b_U,@(x)non_con_tom(x,p),@(x)non_con_jacobian_tom(x,p),[],[],c_L,c_U);


Prob.PriLevOpt = 3;
Prob.FEASTOL = 1e-4;
Prob.OPTTOL = 1e-4;
Prob.XTOL = 1e-4;
Prob.KNITRO.options.MAXIT = 300;

Result = tomRun('knitro',Prob,1);

xopt = Result.x_k;
fopt = Result.f_k;

xd = xopt(1:p.n_param);
index = find(xd<0.001);
x_state = xopt(p.n_param+1:end);
x_state_frame = reshape(x_state, p.ns, p.nt);

error = sqrt(fopt);
precision = 1/error;
log_precision = log10(1/error);

sensitivity_num = (x_state_frame(4,(p.t1+p.peak)/p.h+1)-x_state_frame(4,p.t1/p.h+1))/...
        x_state_frame(4,p.t1/p.h+1);
sensitivity_den = (p.I2-p.I1)/p.I1;
sensitivity = sensitivity_num/sensitivity_den;

xd_lhs(:,j) = xd;
index_lhs{j} = index;
log_precision_lhs(j) = log_precision;
sensitivity_lhs(j) = sensitivity;
x_state_lhs{j} = x_state_frame';
tout_lhs{j}= tout;

exitflag_lhs(j) = Result.ExitFlag;
inform_lhs(j) = Result.Inform;
iter_lhs(j) = Result.Iter;
end

feasible_index = find(exitflag_lhs == 0); % actuall optimal
feasible_log_precision = log_precision_lhs(feasible_index);
feasible_tout_lhs = tout_lhs(feasible_index);
feasible_x_state_lhs = x_state_lhs(feasible_index);

for k = 1:length(feasible_index)
h = figure;
plot(feasible_tout_lhs{k},feasible_x_state_lhs{k}(:,1),'b:','LineWidth',2);
hold on
plot(feasible_tout_lhs{k},feasible_x_state_lhs{k}(:,2),'r--','LineWidth',2);
hold on 
plot(feasible_tout_lhs{k},feasible_x_state_lhs{k}(:,3),'c.','LineWidth',2);
xlabel('Time','FontSize',14);
hold on
plot(feasible_tout_lhs{k},feasible_x_state_lhs{k}(:,4),'g-','LineWidth',2);
xlabel('Time','FontSize',14);

ylabel('Concentration','FontSize',14);
legend('A','B','C','D');

save2pdf(['Adaptation_Fig_',num2str(k)],h,600);
end

toc
gene.xd0 = xd0_lhs;
gene.xd = xd_lhs;
gene.index = index_lhs;
gene.log_precision_lhs = log_precision_lhs;
gene.sensitivity_lhs = sensitivity_lhs;
gene.exitflag = exitflag_lhs;
gene.inform = inform_lhs;
gene.iter = iter_lhs;
gene.x_state = x_state_lhs;
gene.tout = tout_lhs;
save('four_node_DT','gene');
end