function [output] = general_parameter_estimation_function(input_file_name,igraph,iestimate);
% USAGE: [output] = general_parameter_estimation_function(input_file_name,igraph,iestimate);
%
% LSE driver for optimizing the weights 
% The global statement is a list of all the variables that are shared by
% each program in the folder
global A prorate degrate active inactive tspan d wts b alpha lse_out penalty_out counter parms parmnames

if nargin<2
    igraph = 1;
    iestimate = 1;
end
if nargin <3
    iestimate = 1;
end

input_file  = [input_file_name '.xls'];

% These variables call data from Excel files
[degrate,TX0] = xlsread(input_file,'degradation_rates');
[wtmat,TX2]   = xlsread(input_file,'network_weights');
[A,TX2]       = xlsread(input_file,'network');
[dd,TX1]      = xlsread(input_file,'log2_concentrations');
[simtime,TX9] = xlsread(input_file,'simulation_times');

nn = length(dd(:,1));

for ii = 1:nn
    TX11{ii,1} = TX1{ii,2};
end

[type, sheets] = xlsfinfo(input_file);

% default optimization parameters
alpha       = 0.5;
kk_max      = 5;
cssigflg    = 0;
b_or_tau    = 1;

parms(1) = b_or_tau;
parms(2) = alpha;
parms(3) = kk_max;
parms(4) = 1e6;
parms(5) = 1e-8;
parms(6) = 1e6;
parms(7) = 1e-8;
parms(8) = 0;
parms(9) = 0;
parms    = parms(:);

parmnames{1,1} = 'optimization_parameter';
parmnames{2,1} = 'b_or_tau';
parmnames{3,1} = 'alpha';
parmnames{4,1} = 'kk';
parmnames{5,1} = 'MaxIter';
parmnames{6,1} = 'TolFun';
parmnames{7,1} = 'MaxFunEval';
parmnames{8,1} = 'TolX';
parmnames{9,1} = 'fix_P';
parmnames{10,1} = 'fix_b';

prorate  = 2*degrate;

csigma      = 'concentration_sigmas';
coptprms    = 'optimization_parameters';
cprorates   = 'production_rates';
ctau        = 'network_thresholds';
cb          = 'network_b';

bflag = 0;

for ii = 1:length(sheets)
    if length(sheets{ii}) == length(csigma)
        if sum(sheets{ii} == csigma) == length(csigma)
            cssigflg = 1;
            [sigmam,TX4] = xlsread(input_file,'concentration_sigmas');
            sigmas = sigmam(2:end,:)';
        end
    end
    if length(sheets{ii}) == length(coptprms)
        if sum(sheets{ii} == coptprms) == length(coptprms)
            [parms0,parmnames0]  = xlsread(input_file,'optimization_parameters');
            for jj = 1:length(parms0);
                parms(jj) = parms0(jj);
                parmnames{jj+1,1} = parmnames0{jj+1,1};
            end
        end
    end
    if length(sheets{ii}) == length(cprorates)
        if sum(sheets{ii} == cprorates) == length(cprorates)
           [prorate,TX5]  = xlsread(input_file,'production_rates');
        end
    end
    if length(sheets{ii}) == length(ctau)
        if sum(sheets{ii} == ctau) == length(ctau)
           [taumat,TX6]  = xlsread(input_file,'network_thresholds');
        end
    end
    if length(sheets{ii}) == length(cb)
        if sum(sheets{ii} == cb) == length(cb)
           [bvec,TX6]  = xlsread(input_file,'network_b');
           bflag = 1;
           b = bvec;
        end
    end

end

b_or_tau = parms(1);
alpha = parms(2);
kk = parms(3);

fix_P = parms(8);
fix_b = parms(9);


output_file = [input_file_name '_estimation_output_' num2str(b_or_tau) '.xls'];
output_mat  = [input_file_name '_estimation_output_' num2str(b_or_tau) '.mat'];


% To be read by the program correctly, prorate and degrate need to be row
% vectors instead of column vectors, so we use the transpose of the array
% obtained from the datafile
prorate         = prorate';
degrate         = degrate';

% We use variables wherever we would need numbers that might change if we
% are using a different network.
n_genes         = length(dd(:,1))-1;
n_active_genes  = length(A(1,:));
n_times         = length(dd(1,:));

tspan           = dd(1,:);
d               = dd(2:end,:)';

nedges          = sum(A(:));
% % % % nparms          = 2*nedges;

% % [active,inactive] = detect_active_genes(TX11,TX2);

active      = 1:n_active_genes;
inactive    = [];

Ar = sum(A,2);
Ai = Ar>0;

no_inputs = find(Ai==0);
i_forced  = find(Ai == 1);
n_forced  = sum(Ai);

[ig,jg] = find(A == 1);

if b_or_tau == 0
    for ii = 1:nedges
        w0(ii,1) = wtmat(ig(ii),jg(ii));
        w0(ii+nedges,1)   = taumat(ig(ii),jg(ii));
    end

    for ii = 1:n_active_genes
        w0(ii+2*nedges) = prorate(ii);
    end
    
    lb = zeros(size(w0));
    ub = 100*ones(size(w0));
    lb(1:2*nedges) = -100*ones(2*nedges,1);
end
if b_or_tau == 1
    for ii = 1:nedges
        w0(ii,1) = wtmat(ig(ii),jg(ii));
    end
    
    if fix_b == 0
        if bflag == 0
            for ii = i_forced
                w0(ii+nedges,1)   = 0;
            end
        else
            for ii = i_forced
                w0(ii+nedges,1)   = bvec(ii);
            end
        end
    end
    if fix_P == 0
        for ii = 1:n_active_genes
            w0(ii+n_forced+nedges)   = prorate(ii);
        end
    end
    lb = zeros(size(w0));
    ub = 10*ones(size(w0));
    lb(1:nedges) = -10*ones(nedges,1);
    if fix_b == 0
        lb(nedges+1:n_forced+nedges) = -10*ones(n_forced,1);
    end
end

x0              = ones(n_active_genes,1);

% Call the least squares error program, store the sum of the squares of the
% errors in lse_0
counter = 0;
L       = general_least_squares_error(w0);
lse_0   = lse_out;
w1      = w0;

if iestimate == 1
    % This performs the optimization
    for kk = 1:kk_max
        w1      = fmincon(@general_least_squares_error,w1,[],[],[],[],lb,ub,[],optimset(parmnames{5,1},parms(4),parmnames{6,1},parms(5),parmnames{7,1},parms(6),parmnames{8,1},parms(7)));
        L       = general_least_squares_error(w1);
        lse_1   = lse_out;
    end

end
 
time        = simtime;
[t,model]   = ode45(@general_network_dynamics,time,x0);

n_model_times = length(t);

tmin = t(1);
tmax = t(end);

 if igraph == 1

     if cssigflg == 1
         concd = 2.^d(:,active);
         error_up = (d(:,active)+1.96*sigmas);
         error_dn = (d(:,active)-1.96*sigmas);
         for ii=1:n_active_genes
             figure(ii+2),hold on
             plot(time,log2(model(:,ii)),tspan,d(:,active(ii)),'o',tspan,error_up(:,ii),'r+',tspan,error_dn(:,ii),'r+','LineWidth',3),axis([tmin tmax -3 3])
             title(TX1{1+(ii),2},'FontSize',16)
             xlabel('Time (minutes)','FontSize',16)
             ylabel('Expression (log2 fold change)','FontSize',16)
         end
     else
         for ii=1:n_active_genes
             figure(ii+2),hold on
             plot(time,log2(model(:,ii)),tspan,d(:,active(ii)),'o','LineWidth',3),axis([tmin tmax -3 3])
             title(TX1{1+(ii),2},'FontSize',16)
             xlabel('Time (minutes)','FontSize',16)
             ylabel('Expression (log2 fold change)','FontSize',16)
         end
     end
 end

outputtimes{1,1} = 'time';

for ii = 1:n_genes+1
    
    outputnet{1,ii} = TX2{1,ii};
    outputnet{ii,1} = TX2{1,ii};
    outputcells{ii,1} = TX0{ii,1};
    outputcells{ii,2} = TX0{ii,2};
    outputdata{ii,1} = TX0{ii,1};
    outputdata{ii,2} = TX0{ii,2};
    outputdeg{ii,1} = TX0{ii,1};
    outputdeg{ii,2} = TX0{ii,2};
    outputpro{ii,1} = TX0{ii,1};
    outputpro{ii,2} = TX0{ii,2};

    if ii>=2
        for jj = 2:n_model_times+1
            outputcells{ii,jj+1} = log2(model(jj-1,ii-1));
        end
        for jj = 2:n_times+1
            outputdata{ii,jj+1} = dd(ii,jj-1);
        end
        for jj = 2:n_genes+1
            outputnet{jj,ii} = A(jj-1,ii-1);
        end
        outputpro{ii,3} = prorate(ii-1);
        outputdeg{ii,3} = degrate(ii-1);
    else
        outputdeg{ii,3} = TX0{ii,3};
        outputpro{ii,3} = cprorates;
        
        for jj = 2:n_model_times+1
            outputcells{ii,jj+1} = time(jj-1);
        end
        for jj = 2:n_times+1
            outputdata{ii,jj+1} = tspan(jj-1);
            outputtimes{1,jj} = tspan(jj-1);
        end        
    end
end
xlswrite(output_file,outputdata,'log2_concentrations');

xlswrite(output_file,outputcells,'log2_optimized_concentrations');
xlswrite(output_file,outputdeg,'degradation_rates');
xlswrite(output_file,outputpro,'production_rates');
xlswrite(output_file,outputtimes,'measurement_times');
xlswrite(output_file,outputnet,'network');

for ii = 1:nedges
        outputnet{ig(ii)+1,jg(ii)+1} = w0(ii);
end
xlswrite(output_file,outputnet,'network_weights');

if b_or_tau == 0
    for ii = 1:nedges
        outputnet{ig(ii)+1,jg(ii)+1} = w0(ii+nedges);
    end
    xlswrite(output_file,outputnet,'network_thresholds');
    
    for ii = 1:nedges
        outputnet{ig(ii)+1,jg(ii)+1} = w1(ii+nedges);
    end
    xlswrite(output_file,outputnet,'network_optimized_thresholds');
end
if b_or_tau == 1
    outputpro{1,3} = 'b';
    if fix_b == 0
        for ii = 1:n_forced
            outputpro{i_forced(ii)+1,3} = w1(ii+nedges);
        end
        for ii = 1:length(no_inputs)
            outputpro{no_inputs(ii)+1,3} = 0;
        end
        xlswrite(output_file,outputpro,'network_optimized_b');
    end
    if fix_b == 1
        for ii = 1:n_forced
            outputpro{i_forced(ii)+1,3} = b(ii);
        end
        for ii = 1:length(no_inputs)
            outputpro{no_inputs(ii)+1,3} = 0;
        end
        xlswrite(output_file,outputpro,'network_b');
    end
end

for ii = 1:nedges
        outputnet{ig(ii)+1,jg(ii)+1} = w1(ii);
end
xlswrite(output_file,outputnet,'network_optimized_weights');

my_string = ['save(''' output_mat ''')'];
eval(my_string);

output.t        = t;
output.model    = model;
output.name     = input_file;
output.prorate = prorate;
output.degrate = degrate;
output.wts     = wts;
output.b       = b;
output.A       = A;
output.active  = active;
output.tspan   = tspan;
output.d       = d;
output.lse_out = lse_out;
output.reg_out = penalty_out;
output.alpha   = alpha;
output.b_or_tau= b_or_tau;

if cssigflg == 1
    
    outsigmas = outputdata;
    for ii = 1:n_genes+1
        if ii>=2
            for jj = 2:n_times+1
                outsigmas{ii,jj+1} = sigmam(ii,jj-1);
            end
        end
    end
    
    xlswrite(output_file,outsigmas,'concentration_sigmas');
    output.sigmas = sigmas;
end
return