function varargout = pcr_sim(varargin)
% PCR_SIM Application M-file for pcr_sim.fig
%    FIG = PCR_SIM launch pcr_sim GUI.
%    PCR_SIM('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 16-Feb-2004 16:40:36


if nargin == 0  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
	guidata(fig, handles);

	if nargout > 0
		varargout{1} = fig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end



% --------------------------------------------------------------------
%This function runs the PCR simulator
function varargout = pcr_button_Callback(h, eventdata, handles, varargin)

    %randomize the seed
rand('state',sum(100*clock))

tic

    %get the values of the inputs
[numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles);

    %get the expected number of calculations
exptdNumCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep);
set(handles.num_calc, 'String', ['Expected number of calculations: ' num2str(round(exptdNumCalc))]);

    %if we don't want to pause, eliminate divide by zero error by setting
    %puaseStep higher than numTrial
if pauseStep==0
    pauseStep = numTrial+1;
end
    
    
for curTrial = 1:numTrial
    
        %start with the initial number of copies
    totalCopy(curTrial) = numCopy;

    %for each iteration, count the number that are successfully copied, 
    %and add to total number of copies... repeat for each iteration
    for curIter = 1:numIter
        
        successfulCopy = 0;

        for i=1:totalCopy(curTrial)
            if rand<=probRep
                successfulCopy=successfulCopy+1;
            end
        end
        
        
        totalCopy(curTrial) = totalCopy(curTrial) + successfulCopy;
    end
    
        %pause if appropriate, and plot histogram
    if (curTrial/pauseStep) == round(curTrial/pauseStep)
        hist(totalCopy, round(sqrt(curTrial)));
        pause(pauseLength);
    end
    
end

toc

hist(totalCopy, round(sqrt(numTrial)));




% --------------------------------------------------------------------
%This function calculates the expected number of calculations to be
%performed.  
%numCalc = #trials*#copies*SUM(1+probRep)^(n-1), for n=1:numIter
function numCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep)
temp=0;
for i = 1:numIter
    temp = temp + (1 + probRep)^(i-1);
end
numCalc = temp * numTrial * numCopy;




% --------------------------------------------------------------------
%This function gets the inputs from appropriate fields on the form, checks
%the input for problems, and rounds the numbers that should be integers
function [numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles)

numIter = str2num(get(handles.num_iter, 'String'));
numCopy = str2num(get(handles.num_copy, 'String'));
probRep = str2num(get(handles.prob_rep, 'String'));
numTrial = str2num(get(handles.num_trial, 'String'));
pauseStep = str2num(get(handles.pause_step, 'String'));
pauseLength = str2num(get(handles.pause_length, 'String'));


if (isempty(numIter) & isempty(numCopy) & isempty(probRep) & isempty(numTrial) & isempty(pauseStep) & isempty(pauseLength))
    numIter = 0;
    numTrial = 0;
    msgbox('One of the entries is not a number.');
end
    
numIter = round(numIter);
numCopy = round(numCopy);
numTrial = round(numTrial);
pauseStep = round(pauseStep);
if probRep > 1 
    probRep = 1;
    set(handles.prob_rep, 'String', '1');
end





% --------------------------------------------------------------------
%When the value in the num_copy box is changed, update the expected 
%number of calculations
function num_copy_Callback(hObject, eventdata, handles)
[numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles);
exptdNumCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep);
set(handles.num_calc, 'String', ['Expected number of calculations: ' num2str(round(exptdNumCalc))]);
if exptdNumCalc > 1E8
    msgbox(['Warning, this could take up to ' num2str(round(exptdNumCalc/1E7)) 's to simulate']);
end

% --------------------------------------------------------------------
%When the value in the prob_rep box is changed, update the expected 
%number of calculations
function prob_rep_Callback(hObject, eventdata, handles)
[numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles);
exptdNumCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep);
set(handles.num_calc, 'String', ['Expected number of calculations: ' num2str(round(exptdNumCalc))]);
if exptdNumCalc > 1E8
    msgbox(['Warning, this could take up to ' num2str(round(exptdNumCalc/1E7)) 's to simulate']);
end

% --------------------------------------------------------------------
%When the value in the num_iter box is changed, update the expected 
%number of calculations
function num_iter_Callback(hObject, eventdata, handles)
[numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles);
exptdNumCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep);
set(handles.num_calc, 'String', ['Expected number of calculations: ' num2str(round(exptdNumCalc))]);
if exptdNumCalc > 1E8
    msgbox(['Warning, this could take up to ' num2str(round(exptdNumCalc/1E7)) 's to simulate']);
end

% --------------------------------------------------------------------
%When the value in the num_trial box is changed, update the expected 
%number of calculations
function num_trial_Callback(hObject, eventdata, handles)
[numIter, numCopy, probRep, numTrial, pauseStep, pauseLength] = get_input(handles);
exptdNumCalc = get_num_calc(numIter, numCopy, probRep, numTrial, pauseStep);
set(handles.num_calc, 'String', ['Expected number of calculations: ' num2str(round(exptdNumCalc))]);
if exptdNumCalc > 1E8
    msgbox(['Warning, this could take up to ' num2str(round(exptdNumCalc/1E7)) 's to simulate']);
end


    
% --------------------------------------------------------------------
function pause_step_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function pause_length_Callback(hObject, eventdata, handles)



% --------------------------------------------------------------------
function num_copy_CreateFcn(hObject, eventdata, handles)

% --------------------------------------------------------------------
function prob_rep_CreateFcn(hObject, eventdata, handles)

% --------------------------------------------------------------------
function num_iter_CreateFcn(hObject, eventdata, handles)

% --------------------------------------------------------------------
function num_trial_CreateFcn(hObject, eventdata, handles)

% --------------------------------------------------------------------
function pause_step_CreateFcn(hObject, eventdata, handles)

% --------------------------------------------------------------------
function pause_length_CreateFcn(hObject, eventdata, handles)