cba_fitness_function

 [f,s,v,fplus,fminus,zx, s, v] = cba_fitness_function(u,c_init,network,constraints)

 Fitness function used in CBA. Parameters are read from the fields of constraints

0001 function [f, s, v, fplus, fminus, zx] = cba_fitness_function(u, c_init, network, constraints)
0002 
0003 % [f,s,v,fplus,fminus,zx, s, v] = cba_fitness_function(u,c_init,network,constraints)
0004 %
0005 % Fitness function used in CBA. Parameters are read from the fields of constraints
0006 
0007 global my_c
0008 global my_v
0009 
0010 c_init = my_c;
0011 
0012 network.kinetics.u = u;
0013 [s, v, sdot]       = network_steady_state(network, c_init);
0014 fplus              = constraints.zv' * v;
0015 if isfield(constraints,'fitness_options'),
0016   fminus = feval(constraints.fitness_options.cost_function,u,constraints.fitness_options);
0017 else
0018   fminus = constraints.fx' * sum(u) + 0.5 * constraints.fxx * sum(u.^2);
0019 end
0020 f                  = fplus - fminus;
0021 
0022 % provide gradient?
0023 if 0,
0024   R   = basic_control_analysis(network,s);
0025   RJu = R.RJ(:,1:length(u));
0026   f_plus_grad  = [constraints.zv' * RJu]';
0027   f_minus_grad = constraints.fx * ones(size(u)) + constraints.fxx * u;
0028   f_grad       = f_plus_grad - f_minus_grad;
0029 end
0030 
0031 if nargout > 5,
0032   nr = length(u);
0033   zx = constraints.fx' * ones(nr,1) + constraints.fxx * u;
0034 end
0035 
0036 % punish solutions violating stationarity
0037 epsilon = 10^-3;
0038 error   = max(abs(sdot ./ [s(find(1-network.external)) + epsilon]));
0039 
0040 % v
0041 % sdot
0042 % s(find(1-network.external))
0043 
0044 if error > 0.0001, 
0045   f = -10^8*error; 
0046 else
0047   my_c = s;
0048   my_v = v;
0049 end
0050 
0051 display(sprintf('f: %f f_bene: %f f_cost: %f',f,fplus, fminus))
0052 if rand < 0.2,
0053   gp = struct('arrowvalues',v,'actstyle','fixed','flag_edges',1);
0054   netgraph_concentrations(network,s,u,1,gp); drawnow
0055 end