function varargout = mixsir(varargin)

% MIXSIR M-file for mixsir.fig
%      MIXSIR, by itself, creates a new MIXSIR or raises the existing
%      singleton*.
%
%      H = MIXSIR returns the handle to a new MIXSIR or the handle to
%      the existing singleton*.
%
%      MIXSIR('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MIXSIR.M with the given input arguments.
%
%      MIXSIR('Property','Value',...) creates a new MIXSIR or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MixSIR_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to mixsir_openingfcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help mixsir

% Last Modified by GUIDE v2.5 25-Feb-2008 17:03:37

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @mixsir_OpeningFcn, ...
                   'gui_OutputFcn',  @mixsir_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before mixsir is made visible.
function mixsir_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to mixsir (see VARARGIN)


% Choose default command line output for mixsir
handles.output = hObject;
priors = 0; % initialize priors (this specifies the uniformative case as default)
save 'priors.mat' priors %saves prior placeholder to working directory
handles.printresult=0; %sets default to not write results files
handles.xscale=1; %by default scale histogram results on interval [0,1]

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes mixsir wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = mixsir_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in "About mixsir...".
function about_Callback(hObject, eventdata, handles)
% hObject    handle to about (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
open about.txt;

function iterations_Callback(hObject, eventdata, handles)
% hObject    handle to iterations (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of iterations as text
handles.iters =str2double(get(hObject,'String'));% returns contents of iterations as a double
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function iterations_CreateFcn(hObject, eventdata, handles)
% hObject    handle to iterations (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in runmodelbutton.
function runmodelbutton_Callback(hObject, eventdata, handles)
% hObject    handle to runmodelbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


%get the user defined number of model iterations to be used
iterations_Callback(hObject, eventdata, handles);
try
if isnan(handles.iters)
   msgbox('You have not entered a valid number into the iterations field.','Oops...','error')
   return;
elseif handles.iters<10
   msgbox('Less than 10,000 iterations not allowed.','Oops...','error')
   return;
elseif mod(handles.iters,1)~=0
    msgbox('iterations must be in increments of 1000. Iterations will be rounded up.','Oops...','warn')
    handles.iters = handles.iters+(1-mod(handles.iters,1)); %ceilings to the nearest 1000
end
catch
   msgbox('You have not entered a valid number into the iterations field.','Oops...','error')
   return;
end
%READ THESE IN...
%but first check to make sure they exist
%_________________________________________________________________________________________
if exist('mean_source.txt','file')==0
    msgbox('The file "mean_source.txt" is missing from the MixSIR directory!')
    return;
else
    load mean_source.txt;       %creates matrix called mean_source
end
if exist('SD_source.txt','file')==0
    msgbox('The file "SD_source.txt" is missing from the MixSIR directory!')
    return;
else
    load SD_source.txt; 		%creates matrix called SD_source
end
if exist('mix_data.txt','file')==0
    msgbox('The file "mix_data.txt" is missing from the MixSIR directory!','File Missing...','error')
    return;
else
    load mix_data.txt;	handles.mean_mix = mean(mix_data,1);%creates matrix called mean_mix
end
if exist('mean_frac.txt','file')==0
    msgbox('The file "mean_frac.txt" is missing from the MixSIR directory!','File Missing...','error')
    return;
else
    load mean_frac.txt;         %creates matrix called mean_frac 
end
if exist('SD_frac.txt','file')==0
    msgbox('The file "SD_frac.txt" is missing from the MixSIR directory!','File Missing...','error')
    return;
else
    load SD_frac.txt;           %creates matrix called SD_frac
end
%convert from SD to variance
var_source = SD_source.^2;
var_frac = SD_frac.^2;
%_________________________________________________________________________________________

% Check that data are in the right format (at least in terms of matrix sizes)
if (size(mean_source,1)-size(SD_source,1) + size(mean_source,2)-size(SD_source,2)) +...
        (size(mean_frac,1)-size(SD_frac,1) + size(mean_frac,2)-size(SD_frac,2))+...
        (size(mix_data,2)-size(SD_frac,2))+...
        (size(mean_source,2)-size(mean_frac,2) + (size(mean_source,1)-size(mean_frac,1))) ~= 0
    msgbox('Dimensions of source data matrices do not agree!','Oops...','error');
    return;
end

%MAKE THESE FOR WORK....
%_________________________________________________________________________________________
handles.source_num = size(mean_source,1);       % how many sources contribute to the mixator spp of concern?
handles.iso_num = size(mean_source,2);    	% how many isotopes are you using?
handles.datasize = size(mix_data,1);        %how many data points are you using?
%set the priors based on user specs
if ~exist('priors.mat','file') %initializes the first time only....
        priors = 0; % toggle for using priors (set to off)
        save 'priors.mat' priors %saves prior placeholder to working directory
end
load priors.mat;  %load prior placeholder file
   if priors == 0; 
      handles.A = ones(1,handles.source_num);%sets uninformative priors as default
      handles.B = ones(1,handles.source_num);%sets uninformative priors as default
      clear priors;
   else
        if exist('define.txt')==2
        load define.txt; %brings in defined priors A and B (columns) for each source (rows)
        handles.A = define(1,:);
        handles.B = define(2,:);
        save 'priors.mat' priors%saves prior placeholder to working directory
        else
            msgbox('ERROR: The file define.txt does not exist in the program directory','File Missing...','error')
        end
   end

source = zeros(handles.source_num, handles.iso_num); 			% this used to store sampled iso vals for source
mix = zeros(1,handles.iso_num); 						% this used to store sampled iso vals for mix
contribution = zeros(1,handles.source_num)';			% this used to store sampled contribs
handles.contrib_out = zeros(1,handles.source_num); 			% this used to store accepted contrib vectors (OUTPUT!)
out_ticker = 0;									% keeps track of the # of SIR accepts
fraction = zeros(handles.source_num,handles.iso_num);					% this used to store sampled fractionation for each iso
Total = 0;											% this used to keep tally of likes for SIR accept
MaxDup = 0;											% keeps track of duplicate accepts. Used to evaluate threshold.
likestor = 0;										% this used to store best threshold guess
ljump = 0;											% keeps track of the # likes accepted
resampold = 0;                                      %used to keep of duplicates
handles.L =zeros(1,1);										% this used to store all the accepted likes
progressbar;                                        %Initializes a GUI progress bar
%_________________________________________________________________________________________

%IDENTIFY AN APPROPRIATE THRESHOLD.....
%_________________________________________________________________________________________

for irep=1:(round(handles.iters*0.1))%uses ~10% run size to get threshold
    
    progressbar(irep/(round(handles.iters*0.1)+handles.iters)); %begins updating progress bar

    % Generate random proportions that sum to one (one for each source)
   contribution= rand_props2(handles.source_num,1000); %calls function randprops2 to make random proportions (1000 columns)
   
   %Below uses an analytic solution to the mean and var of pred (without frationation) based on combining source
   %means and vars.  Note that I have implemented a vectorized version of this calc below, so beware of sneaky matrix math!
    mproposal = (mean_source'*contribution)';%generates a row of proposed mean mixs (one for each iso)
    vproposal = ((var_source)'*(contribution.^2))'; %does same for var
    
   % Now use the same method as above to get proposed fractionation
   mfraction = (mean_frac'*contribution)';%generates a row of proposed fractionations (one for each iso)
   vfraction = ((var_frac)'*(contribution.^2))'; %does same for var
   
   %when adding distributions, both the mean and variance are additive. So
   %we simply add the fractionation mean and var to the mean and var of the
   %mix calculated above
   mproposal=mproposal+mfraction; vproposal=vproposal+vfraction;
      
	% EVALUATE THE NEGLOGLIKE FOR THE PROPOSED MIXES BASED ON DATA
   neglog = zeros(1000,1); % clear to zero for each run
   for dloop = 1:handles.datasize %loop through all mix data points, calc NNLs for each of the 1000 draws
   for iso_i = 1:handles.iso_num %loops over # of isotopes
       z = (mix_data(dloop,iso_i)-mproposal(:,iso_i))./vproposal(:,iso_i);
       neglog = neglog + (.5.*z.*z + log(sqrt(2.*pi).*vproposal(:,iso_i)));
   end
   end
   
   % EVALUATE THE NEGLOGLIKE FOR source CONTRIBS BASED ON PRIORS
   for ploop = 1:handles.source_num
    neglog = neglog - (log(betapdf(contribution(ploop,:)',handles.A(1,ploop),handles.B(1,ploop))+1e-100));
   end
   
    likelihood = exp(-(neglog)); %presto! The likelihoods appear like magic
   
   if max(likelihood) > likestor
      likestor = max(likelihood); %keeps the very best like found
   end
end

Threshold = likestor*2;			% Define SIR accept. criteria (2x max like for safety)
%_________________________________________________________________________________________

%RUN THE SIR ALGORITHM....
%_________________________________________________________________________________________

% DO for each replicate
for irep = 1:handles.iters
   
    progressbar((round(handles.iters*0.1)+irep)/(round(handles.iters*0.1)+handles.iters)); %update progress bar
    
    % Generate random proportions that sum to one (one for each source)
   contribution= rand_props2(handles.source_num,1000); %calls function randprops2 to make random proportions (1000 columns)
   
   %Below uses an analytic solution to the mean and var of pred (without frationation) based on combining source
   %means and vars.  Note that I have implemented a vectorized version of this calc below, so beware of sneaky matrix math!
    mproposal = (mean_source'*contribution)';%generates a row of proposed mean mixs (one for each iso)
    vproposal = ((var_source)'*(contribution.^2))'; %does same for var
    
   % Now use the same method as above to get proposed fractionation
   mfraction = (mean_frac'*contribution)';%generates a row of proposed fractionations (one for each iso)
   vfraction = ((var_frac)'*(contribution.^2))'; %does same for var
   
   %when adding distributions, both the mean and variance are additive. So
   %we simply add the fractionation mean and var to the mean and var of the
   %mix calculated above
   mproposal=mproposal+mfraction; vproposal=vproposal+vfraction;
      
	% EVALUATE THE NEGLOGLIKE FOR THE PROPOSED MIXES BASED ON DATA
   neglog = zeros(1000,1); % clear to zero for each run
   for dloop = 1:handles.datasize %loop through all mix data points, calc NNLs for each of the 1000 draws
   for iso_i = 1:handles.iso_num %loops over # of isotopes
       z = (mix_data(dloop,iso_i)-mproposal(:,iso_i))./vproposal(:,iso_i);
       neglog = neglog + (.5.*z.*z + log(sqrt(2.*pi).*vproposal(:,iso_i)));
   end
   end
   
   % EVALUATE THE NEGLOGLIKE FOR source CONTRIBS BASED ON PRIORS
   for ploop = 1:handles.source_num
    neglog = neglog - (log(betapdf(contribution(ploop,:)',handles.A(1,ploop),handles.B(1,ploop))+1e-100));
   end
   
    likelihood = exp(-(neglog)); %presto! The likelihoods appear like magic
    
   % HERE COMES THE Hilborn SIR... (mostly vectorized for speed)
   Ltemp = likelihood(1,1); %stores first like val
   likelihood(1,1)=likelihood(1,1)+Total; %carries the total over from the previous 1000 vector
   cumu_full= cumsum(likelihood)./Threshold; %gets the cumulative sum of the like vector
   likelihood(1,1) = Ltemp; %puts it back the way it was
   cumu = floor(cumu_full); %floor to integers for cell refs
   cumu(2:1000) = cumu(2:1000) - cumu(1:999); %sets references for the resample off of the like vector
   
   if max(cumu)>0 %if there are any resamples to be had in this 1000 vector...
    if max(cumu)>MaxDup; MaxDup = max(cumu); end %keeps track of the max duplicate draws
   
    exnd = nonzeros(cumu); %gets #  to be taken for each resample
    refvec = find(cumu>0); %gets the cell refences for the resample values
        for zz = 1:size(refvec)%loops over the resamples
            for qq = 1:exnd(zz)%loops over the number of redraws for each resample (in most cases it should be =1)
               out_ticker = out_ticker+1; %records the total number of resamples thus far...
               handles.contrib_out(out_ticker,:) = contribution(:,refvec(zz))'; %records all resampled contribs
               handles.L(out_ticker,:) = likelihood(refvec(zz))';%records all resamlped likelihoods
            end
        end
        
     Total = ((1+mod(cumu_full(max(find(cumu>0))),1))*Threshold) - Threshold+...
     sum(likelihood(refvec(end)+1:1000)); %sets up total for next go round
   else
       Total = Total+sum(likelihood);%sets up total for next go round
   end
   
      
   
end % ends iterations


%DEVELOP RESULTS OF MODEL RUN BELOW.....
if handles.xscale==1
    histbin = 0:.01:1;
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii,'Parent',handles.g_panel);
        [y,x] = hist(handles.contrib_out(:,ii),histbin); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        axis([0 1 -Inf Inf])
        end
else
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii,'Parent',handles.g_panel);
        [y,x] = hist(handles.contrib_out(:,ii),100); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        end
end
    
if handles.printresult==1
    % write out results
    dlmwrite('contrib_out.txt', handles.contrib_out, '\t')
    dlmwrite('likelihoods.txt', handles.L, '\t')
end
handles.runique = length(unique(handles.contrib_out,'rows'));% # of unique parameter vectors in resample
handles.Lmax = max(handles.L)/sum(handles.L); % best likelihood (% of total)
printout(irep,out_ticker,MaxDup,handles.contrib_out,handles.source_num, handles.runique,handles.Lmax);

guidata(hObject, handles);
% --- Executes on button press in hist_butt.

function hist_butt_Callback(hObject, eventdata, handles)
% hObject    handle to hist_butt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if handles.boxprint == 1 %if user requests figure in its own window...
    figure  %first make the figure window
    if handles.xscale==1
        histbin = 0:.01:1;
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii); 
        [y,x] = hist(handles.contrib_out(:,ii),histbin); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        axis([0 1 -Inf Inf])
        end
    else
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii); 
        [y,x] = hist(handles.contrib_out(:,ii),100); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        end
    end
else
    % plot pairwise contours of spp posterior proportions
    if handles.xscale==1
        histbin = 0:.01:1;
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii); 
        [y,x] = hist(handles.contrib_out(:,ii),histbin); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        axis([0 1 -Inf Inf])
        end
    else
        for ii = 1:handles.source_num
        a1=subplot(handles.source_num,1,ii); 
        [y,x] = hist(handles.contrib_out(:,ii),100); %this step does the division in Bayes theorem
        bar(a1,x,(y./sum(y)),1);
        if ii == handles.source_num; xlabel('Source Contributions to the Mix'); end %put in x axis label
        if ii == 1+floor(handles.source_num/2); ylabel('Posterior Probability');end %put in y axis label
        end
    end
end


% --- Executes on button press in likelihood_distrib.
function likelihood_distrib_Callback(hObject, eventdata, handles)
% hObject    handle to likelihood_distrib (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% all accepted likelihoods relative to max accepted like (MLE)
if handles.likeprint == 1 %if user requests figure in separate window...
    figure; %first make the window
    hist(handles.L/max(handles.L),20);	
    h = findobj(gca,'Type','patch');
    set(h,'FaceColor','y','EdgeColor','k')
    xlabel('Relative Likelihood')
    ylabel('Count')
    title('Likelihoods for SIR Resample')
else
    a1=subplot(1,1,1,'Parent',handles.g_panel); %use subplot because it automatically writes over old plot in uipanel
    hist(a1,handles.L/max(handles.L),20);	
    h = findobj(gca,'Type','patch');
    set(h,'FaceColor','y','EdgeColor','k')
    xlabel('Relative Likelihood')
    ylabel('Count')
    title('Likelihoods for SIR Resample')
end

% --- Executes on button press in results_check.
function results_check_Callback(hObject, eventdata, handles)
% hObject    handle to results_check (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.printresult =get(hObject,'Value');% returns toggle state of results_check
guidata(hObject,handles);

% --- Executes on button press in boxprint.
function boxprint_Callback(hObject, eventdata, handles)
% hObject    handle to boxprint (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.boxprint =get(hObject,'Value');% returns toggle state of results_check
guidata(hObject,handles);

% --- Executes on button press in hist_x_scale.
function hist_x_scale_Callback(hObject, eventdata, handles)
% hObject    handle to hist_x_scale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.xscale =get(hObject,'Value');% returns toggle state of results_check
guidata(hObject,handles);

% --- Executes on button press in likeprint.
function likeprint_Callback(hObject, eventdata, handles)
% hObject    handle to likeprint (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.likeprint =get(hObject,'Value');% returns toggle state of results_check
guidata(hObject,handles);

% --------------------------------------------------------------------
function priorpanel_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to priorpanel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

switch get(hObject,'Tag')   % Get Tag of selected object
    case 'uninformative'
        % code piece when radiobutton1 is selected goes here
        priors = 0; % 0 value set prior distrib to Beta(1,1) upon program execution
        save 'priors.mat' priors %saves prior placeholder to working directory
    case 'define'
        % code piece when radiobutton2 is selected goes here
        if exist('define.txt')==2
        load define.txt; %brings in defined priors A and B (columns) for each source (rows)
        priors = 1;  %1 value sets program to draw priors from prior file (define.txt) each model run
        save 'priors.mat' priors%saves prior placeholder to working directory
        else
            msgbox('ERROR: The file define.txt does not exist in the program directory','File Missing...','error')
        end
end

  
%--- STATUS BAR BELOW
function progressbar(fractiondone, position)
% Description:
%   progressbar(fractiondone,position) provides an indication of the progress of
% some task using graphics and text. Calling progressbar repeatedly will update
% the figure and automatically estimate the amount of time remaining.
%   This implementation of progressbar is intended to be extremely simple to use
% while providing a high quality user experience.
%
% Features:
%   - Can add progressbar to existing m-files with a single line of code.
%   - The figure closes automatically when the task is complete.
%   - Only one progressbar can exist so old figures don't clutter the desktop.
%   - Remaining time estimate is accurate even if the figure gets closed.
%   - Minimal execution time. Won't slow down code.
%   - Random color and position options. When a programmer gets bored....
%
% Usage:
%   fractiondone specifies what fraction (0.0 - 1.0) of the task is complete.
% Typically, the figure will be updated according to that value. However, if
% fractiondone == 0.0, a new figure is created (an existing figure would be
% closed first). If fractiondone == 1.0, the progressbar figure will close.
%   position determines where the progressbar figure appears on screen. This
% argument only has an effect when a progress bar is first created or is reset
% by calling with fractiondone = 0. The progress bar's position can be specifed
% as follows:
%       [x, y]  - Position of lower left corner in normalized units (0.0 - 1.0)
%           0   - Centered (Default)
%           1   - Upper right
%           2   - Upper left
%           3   - Lower left
%           4   - Lower right
%           5   - Random [x, y] position
%   The color of the progressbar is choosen randomly when it is created or
% reset. Clicking inside the figure will cause a random color change.
%   For best results, call progressbar(0) (or just progressbar) before starting
% a task. This sets the proper starting time to calculate time remaining.
%
% Example Function Calls:
%   progressbar(fractiondone,position)
%   progressbar               % Initialize/reset
%   progressbar(0)            % Initialize/reset
%   progressbar(0,4)          % Initialize/reset and specify position
%   progressbar(0,[0.2 0.7])  % Initialize/reset and specify position
%   progressbar(0.5)          % Update
%   progressbar(1)            % Close
%
% Demo:
%   n = 1000;
%   progressbar % Create figure and set starting time
%   for i = 1:n
%       pause(0.01) % Do something important
%       progressbar(i/n) % Update figure
%   end
%
% Author: Steve Hoelzer
%
% Revisions:
% 2002-Feb-27   Created function
% 2002-Mar-19   Updated title text order
% 2002-Apr-11   Use floor instead of round for percentdone
% 2002-Jun-06   Updated for speed using patch (Thanks to waitbar.m)
% 2002-Jun-19   Choose random patch color when a new figure is created
% 2002-Jun-24   Click on bar or axes to choose new random color
% 2002-Jun-27   Calc time left, reset progress bar when fractiondone == 0
% 2002-Jun-28   Remove extraText var, add position var
% 2002-Jul-18   fractiondone input is optional
% 2002-Jul-19   Allow position to specify screen coordinates
% 2002-Jul-22   Clear vars used in color change callback routine
% 2002-Jul-29   Position input is always specified in pixels
% 2002-Sep-09   Change order of title bar text
% 2003-Jun-13   Change 'min' to 'm' because of built in function 'min'
% 2003-Sep-08   Use callback for changing color instead of string
% 2003-Sep-10   Use persistent vars for speed, modify titlebarstr
% 2003-Sep-25   Correct titlebarstr for 0% case
% 2003-Nov-25   Clear all persistent vars when percentdone = 100
% 2004-Jan-22   Cleaner reset process, don't create figure if percentdone = 100
% 2004-Jan-27   Handle incorrect position input
% 2004-Feb-16   Minimum time interval between updates
% 2004-Apr-01   Cleaner process of enforcing minimum time interval
% 2004-Oct-08   Seperate function for timeleftstr, expand to include days
% 2004-Oct-20   Efficient if-else structure for sec2timestr
%

persistent progfig progpatch starttime lastupdate

% Set defaults for variables not passed in
if nargin < 1
    fractiondone = 0;
end
if nargin < 2
    position = 0;
end

try
    % Access progfig to see if it exists ('try' will fail if it doesn't)
    dummy = get(progfig,'UserData');
    % If progress bar needs to be reset, close figure and set handle to empty
    if fractiondone == 0
        delete(progfig) % Close progress bar
        progfig = []; % Set to empty so a new progress bar is created
    end
catch
    progfig = []; % Set to empty so a new progress bar is created
end

% If task completed, close figure and clear vars, then exit
percentdone = floor(100*fractiondone);
if percentdone == 100 % Task completed
    delete(progfig) % Close progress bar
    clear progfig progpatch starttime lastupdate % Clear persistent vars
    return
end

% Create new progress bar if needed
if isempty(progfig)
    
    % Calculate position of progress bar in normalized units
    scrsz = [0 0 1 1];
    width = scrsz(3)/4;
    height = scrsz(4)/50;
    if (length(position) == 1)
        hpad = scrsz(3)/64; % Padding from left or right edge of screen
        vpad = scrsz(4)/24; % Padding from top or bottom edge of screen
        left   = scrsz(3)/2 - width/2; % Default
        bottom = scrsz(4)/2 - height/2; % Default
        switch position
            case 0 % Center
                % Do nothing (default)
            case 1 % Top-right
                left   = scrsz(3) - width  - hpad;
                bottom = scrsz(4) - height - vpad;
            case 2 % Top-left
                left   = hpad;
                bottom = scrsz(4) - height - vpad;
            case 3 % Bottom-left
                left   = hpad;
                bottom = vpad;
            case 4 % Bottom-right
                left   = scrsz(3) - width  - hpad;
                bottom = vpad;
            case 5 % Random
                left   = rand * (scrsz(3)-width);
                bottom = rand * (scrsz(4)-height);
            otherwise
                warning('position must be (0-5). Reset to 0.')
        end
        position = [left bottom];
    elseif length(position) == 2
        % Error checking on position
        if (position(1) < 0) | (scrsz(3)-width < position(1))
            position(1) = max(min(position(1),scrsz(3)-width),0);
            warning('Horizontal position adjusted to fit on screen.')
        end
        if (position(2) < 0) | (scrsz(4)-height < position(2))
            position(2) = max(min(position(2),scrsz(4)-height),0);
            warning('Vertical position adjusted to fit on screen.')
        end
    else
        error('position is not formatted correctly')
    end
    
    % Initialize progress bar
    progfig = figure(...
        'Units',            'normalized',...
        'Position',         [position width height],...
        'NumberTitle',      'off',...
        'Resize',           'off',...
        'MenuBar',          'none',...
        'BackingStore',     'off' );
    progaxes = axes(...
        'Position',         [0.02 0.15 0.96 0.70],...
        'XLim',             [0 1],...
        'YLim',             [0 1],...
        'Box',              'on',...
        'ytick',            [],...
        'xtick',            [] );
    progpatch = patch(...
        'XData',            [0 0 0 0],...
        'YData',            [0 0 1 1],...
        'EraseMode',        'none' );
    set(progfig,  'ButtonDownFcn',{@changecolor,progpatch});
    set(progaxes, 'ButtonDownFcn',{@changecolor,progpatch});
    set(progpatch,'ButtonDownFcn',{@changecolor,progpatch});
    changecolor(0,0,progpatch)
    
    % Set time of last update to ensure a redraw
    lastupdate = clock - 1;
    
    % Task starting time reference
    if isempty(starttime) | (fractiondone == 0)
        starttime = clock;
    end
    
end

% Enforce a minimum time interval between updates
if etime(clock,lastupdate) < 0.01
    return
end

% Update progress patch
set(progpatch,'XData',[0 fractiondone fractiondone 0])

% Update progress figure title bar
if (fractiondone == 0)
    titlebarstr = ' 0%';
else
    runtime = etime(clock,starttime);
    timeleft = runtime/fractiondone - runtime;
    timeleftstr = sec2timestr(timeleft);
    titlebarstr = sprintf('%2d%%    %s remaining',percentdone,timeleftstr);
end
set(progfig,'Name',titlebarstr)

% Force redraw to show changes
drawnow

% Record time of this update
lastupdate = clock;


% ------------------------------------------------------------------------------
function changecolor(h,e,progpatch)
% Change the color of the progress bar patch

colorlim = 2.8; % Must be <= 3.0 - This keeps the color from being too light
thiscolor = rand(1,3);
while sum(thiscolor) > colorlim
    thiscolor = rand(1,3);
end
set(progpatch,'FaceColor',thiscolor);


% ------------------------------------------------------------------------------
function timestr = sec2timestr(sec)
% Convert a time measurement from seconds into a human readable string.

% Convert seconds to other units
d = floor(sec/86400); % Days
sec = sec - d*86400;
h = floor(sec/3600); % Hours
sec = sec - h*3600;
m = floor(sec/60); % Minutes
sec = sec - m*60;
s = floor(sec); % Seconds

% Create time string
if d > 0
    if d > 9
        timestr = sprintf('%d day',d);
    else
        timestr = sprintf('%d day, %d hr',d,h);
    end
elseif h > 0
    if h > 9
        timestr = sprintf('%d hr',h);
    else
        timestr = sprintf('%d hr, %d min',h,m);
    end
elseif m > 0
    if m > 9
        timestr = sprintf('%d min',m);
    else
        timestr = sprintf('%d min, %d sec',m,s);
    end
else
    timestr = sprintf('%d sec',s);
end