function L = general_least_squares_error(theta) % USAGE L = general_least_squares_error(theta) % This program is used with fminsearch to find the optimum weights vector. global A prorate degrate active inactive tspan d wts b alpha lse_out penalty_out counter parms parmnames counter = counter + 1; % % n = length(theta); % % wts = theta(1:n/2); % % b = theta(n/2+1:n); b_or_tau = parms(1); alpha = parms(2); kk = parms(3); fix_P = parms(8); fix_b = parms(9); nedges = sum(A(:)); Ar = sum(A,2); Ai = Ar>0; no_inputs = find(Ai==0); i_forced = find(Ai == 1); n_forced = sum(Ai); n_active_genes = length(active); wts = theta(1:nedges); if b_or_tau == 0 b = theta(nedges+1:2*nedges); if fixP == 0 prorate = theta(2*nedges+1:2*nedges+n_active_genes); end end if b_or_tau == 1 if fix_b == 0 b = theta(nedges+1:n_forced+nedges); end if fix_P==0 prorate = theta(n_forced+nedges+1:nedges+n_forced+n_active_genes); end end % % % % % % if b_or_tau == 0 % % % % % % b = theta(nedges+1:2*nedges); % % % % % % prorate = theta(2*nedges+1:2*nedges+n_active_genes); % % % % % % end % % % % % % if b_or_tau == 1 % % % % % % b = theta(nedges+1:n_forced+nedges); % % % % % % prorate = theta(n_forced+nedges+1:nedges+n_forced+n_active_genes); % % % % % % end % Initial concentrations x0 = ones(n_active_genes,1); d1 = d(:,active); if tspan(1) > 1e-6 tspan1 = [0;tspan(:)]; addzero = 1; else tspan1 = tspan; addzero = 0; end % Matlab uses the o.d.e. solver function to obtain the data from our model [t,x] = ode45('general_network_dynamics',tspan1,x0); %disp('afterode45') % This is the square of the errors (the difference bwtn our data and % Schade's data) if addzero == 1; x1 = x(2:end,:); else x1 = x; end error = ((log2(x1) - d1).^2); % This is the sum of all the errors L = sum(error(:)); % Store the value of L because we're going to change it lse_out = L; % Add the penalty term to get the bowl shape for better optimization bp = 1; if b_or_tau == 1 bp = 0; end if fix_b == 1 bp = b; end proratep = prorate; if fix_P == 1 proratep = 0*prorate; end L = L + alpha*wts'*wts + alpha*(b-bp)'*(b-bp)*4 + 0.01*alpha*sum(proratep(:))^2; % Store the penalty value penalty_out = wts'*wts + (b-bp)'*(b-bp)*4 + 0.01*sum(proratep(:))^2;; if(100*round(counter/100)==counter) figure(1),subplot(211),plot(theta,'d'), title(['counter = ' num2str(counter) ' lseout = ' num2str(lse_out)]) subplot(212),plot(d1,'*'),hold on,plot(log2(x1)), hold off,pause(.1) end return