%Make multiple runs of ReactionSim and plot results.

fclose('all');
fignum = 0;
clear;

%Creates a new file if the old one doesn't exist
filename = input('Project name: ','s');

A = exist([filename '.m'],'file');
if A == 2 %if the file is found it does nothing
else  %else if the file is not found, a default file is made
    make_file(filename); %make_file creates a .m file for the user to enter their reactions
end
disp(' ');
disp('The file you just named should pop up.');
disp('Make changes to this example to create your model.');
disp('Save and close the file once you''re done.');

st = ['!' filename '.m '];  %The .m file pops up

eval(st);

%The .m file is called as a function
s = sprintf('[CV reactions rev_reactions init unSortedSpecies start stop increment] = %s();',filename);
eval(s);
[rWells cWells] = size(init(:,:,1));
nWells = rWells*cWells;

if nWells > 1
    disp(' ');
    disp('Mulitple wells were detected.') %because reactionsim organizes the outputs in alphabetical order
    disp('How many iterations of simple diffusion')
    iterations = input('would you like to run? ');

    disp(' ');
    disp('How many trials (# of simulations for multiple wells');
    trials = input('and all diffusions)? ');

    disp(' ');
    disp('Choose your dispersion pattern:')
    disp('1 - Input file name')
    disp('2 - Simple diffusion')
    pattern = input('Pattern choice: ');

    nRuns = 1;
    if pattern==1
        patternFilename = input('Input filename for pattern file: ','s');
        excel=fopen(patternFilename,'r');       %open file for reading
        test=fscanf(excel,'%s %s %s %s %s',5);     %get rid of header
        C = textscan(excel, '%*s%c%f%*s%f%*s%c%f', 'delimiter', ',','treatAsEmpty', {'NA', 'na'});
        fclose(excel);

    elseif pattern==2
        wellToWell = input('How much to keep in same well (uL): ');
        wellToAdj = input('How much to put in adjacent wells (uL): ');

    end
else
    iterations = 1;
    nRuns = input('How many simulations per well? ');
    trials = 1;
    pattern = 3;
end

tic



%This block sorts the species, entered by the user, into alphabetical order
%and sorts the initial conditions to match the ordered species
species = sort(unSortedSpecies);
nvars = length(unSortedSpecies);

%preallocating DiffMatrix
DiffMatrix = cell(iterations+1);
for SDiff = 1:iterations+1
    DiffMatrix{SDiff} = zeros(rWells,cWells,nvars);
end


for wr = 1:rWells
    for wc = 1:cWells
        for ind = 1:nvars
            for rand_ind = 1:nvars
                if species{ind} == unSortedSpecies{rand_ind}
                    DiffMatrix{1}(wr,wc,ind) = init(wr,wc,rand_ind); %matches the orderedspecies index to the init index
                end
            end
        end
    end
end

%This block saves the constant initial values

%preallocating const
const = zeros(rWells,cWells,nvars);

init = DiffMatrix{1};
count = 0;
for num = 1:length(CV)
    for ind = 1:nvars
        if species{ind} == CV{num}
            const(:,:,ind) = init(:,:,ind);
            count = count + 1;
            match(count) = ind;
        end
    end
end

ninst = (stop - start)/increment + 1;


%preallocating Sdata
Sdata = zeros(nRuns,ninst,nvars);

%preallocating Wells
Wells = zeros(1,1,nvars);

%preallocating Out
Out = zeros(rWells,cWells,nvars);

%preallocating ReducedWells
ReducedWells = zeros(rWells,cWells,nvars);

%SDiff means Simple Diffusions
%If the user wanted 3 diffusions, we would have 3 plates for each
%diffusion; therefore, you can think of SDiff as the plate number that we're
%currently on
for trial = 1:trials;
    init = DiffMatrix{1};
    for SDiff = 1:iterations
        plotcount = 1;
        for wr = 1:rWells
            for wc = 1:cWells

                disp(sprintf('%s %d %s %d %s %d %s','Currently processing iteration ', SDiff, ' Well(', wr,',',wc,')'));
                %Resets the initial conditions for the constant value species
                if nWells > 1
                for count = 1:length(match)
                    init(:,:,match(count)) = const(:,:,match(count));
                end
                end

                % write_file writes the information entered by the user into a
                % text file that hybridsim can read

                write_file(wr,wc,species,init,reactions,rev_reactions,start,stop,increment);
                for runs=1:nRuns
                    randnum = round(10000*rand);
                    s=sprintf('%s %d %s','!..\bin\ReactionSim cvode temporary.txt -runs 1 -batches temp -seed ', randnum, '> garbage.txt');
                    eval(s);
                    %Read batches output file produced by ReactionSim
                    fileID=fopen('temp','r');       %open file for reading
                    test=fscanf(fileID,'%s',1);     %get rid of header
                    A=fscanf(fileID,'%*d,%*d,%*s%e',[nvars,ninst]);

                    for ind=1:nvars
                        for time=1:ninst
                            Sdata(runs,time,ind) = A(ind,time); %A 3d matrix called Sdata is created to store the run(runs), time, and species
                        end
                    end
                    fclose(fileID);
                end

                %condensed averages loops
                if nRuns == 1
                    for ind = 1:nvars
                        for time = 1:ninst
                            averages(1,time,ind) = A(ind,time);
                            AllInst{trial}{ind,time,SDiff}(wr,wc) = A(ind,time);
                        end
                    end
                else  %nRuns might always be 1 for multiple iterations but it doesn't have to be
                    averages = mean(Sdata);
                    for ind = 1:nvars
                        for time = 1:ninst
                            AllInst{trial}{ind,time,SDiff}(wr,wc) = averages(1,time,ind);
                        end
                    end
                end

                %readerplots organizes the plotting if there is only 1 Well
                %It also does the histograms and error bars
                if nWells == 1
                    toc
                    readerplots(Sdata,species,ninst,nvars,nRuns,averages)
                else



                    %If the user enters only 1 run, the values
                    %sent from hybridsim will be put into the averages matrix

                    %Stores the percentages of each well for mixing that will occur later

                    for ind = 1:nvars
                        Wells(wr,wc,ind) = averages(1,ninst,ind);
                        if pattern == 2;
                            Out(wr,wc,ind) = (wellToAdj/200).*Wells(wr,wc,ind);
                            ReducedWells(wr,wc,ind) = (wellToWell/200).*Wells(wr,wc,ind);
                        end
                    end
                end
            end
        end
        DiffMatrix{SDiff+1} = Wells;
        %At this point, a matrix called Wells is created with all the final
        %values of the individual wells stored.
        %Now we are going to mix them

        %mixes based on user input above and only mixes for multiple wells
        if pattern==1 && nWells > 1
            [init] = patternFileReader(Wells,C,rWells,cWells,nvars);

        elseif pattern == 2 && nWells > 1
            [init] = mix_wells(nvars,rWells,cWells,ReducedWells,Out);
        end
    end
end




if nWells > 1

    disp(' ');
    disp('The amount of time the simulation took to run: ');
    toc

    %setting up the range of values that were generated for plotting
    %purposes.  They're stored into clims, which is used for the color bar range.

    for ind = 1:nvars
        for plate = 1:(iterations)
            for time = 1:ninst
                large(ind,time,plate) = max(max(AllInst{1}{ind,time,plate,1}));
            end
        end

        largest(ind) = max(max(max(large(ind,:,:))));

        for plate = 1:(iterations+1)

            if largest(ind) == 0
                largest(ind) = 1;
            else
            end
            clims{ind} = [0 largest(ind)];

        end
    end


    disp(' ');
    disp('Copy and paste ');
    disp('''options_menu(" ");''');
    disp('from readerworks.m ');
    disp('into the command window if you encounter any ');
    disp('errors while trying to use the menu');
    disp(' ');
    disp('Menu: ');



    %Copy and paste "options_menu(" ");"
    %into the command window
    %if you encounter any errors while trying to view the plots
    %and movies

    options_menu(nvars,species,iterations,ninst,AllInst,clims,DiffMatrix,trials)
end

