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