classdef WriFish < matlab.mixin.SetGet & dynamicprops
    % The software models the change in frequency of the B allele subject to initial f(B) and genotype fitness rules.
    
    %-----------------------------------------------
    % Author:           Adam Gogacz $ Steven M. Carr
    % Version:          1.0
    % Last Revision:    2020-10-30
    % Copyright:        Adam Gogacz
    %-----------------------------------------------
    
    %% --- parameters and constants
    properties (SetAccess=private,SetObservable=true,Hidden=true)
        %---- modelling parameters and variables
        q_init      =0.5;       % f(a) = q;
        w0          =1;      % fitness of AA genotype; must be between 0 and Inf
        w1          =1;         % fitness of AB genotype; must be between 0 and Inf
        w2          =1;      % fitness of BB genotype; must be between 0 and Inf
        Nindiv      =100;       % number of individuals in the population
        Ngen        =200;       % number of generations
        Npop        =100;        % number of populations
        N_BB        double;
        N_AB        double;
        qT          double;
        %--- GUI controls
        axes        matlab.ui.control.UIAxes;
    end
    
    %--- GUI constants
    properties (Constant=true,Hidden=true)
        uiHt=20;
        uiGap=10;
        uiHelp=15;
        Name='Wright-Fisher Model';
    end
    
    methods
        function this=WriFish(varargin)
            matlabScreenrSize=get(0,'screensize'); scW=matlabScreenrSize(1,3); scH=matlabScreenrSize(1,4); figTag='wrightfisher';
            delete(findall(0,'Type','figure','-and','Tag',figTag));
            hf=uifigure('Name',this.Name,'Tag',figTag,'Position',[1,0.25*scH,scW,0.75*scH-36],'Visible',0);
            mainGridLayout=matlab.ui.container.GridLayout('Parent',hf,'RowHeight',{'1x','1x'},'ColumnWidth',{220,'1x','1x'},...
                'RowSpacing',48,'ColumnSpacing',10,'Scrollable','off');
            this.axes=matlab.ui.control.UIAxes('Parent',mainGridLayout); this.axes.Layout.Row=[1,2]; this.axes.Layout.Column=[2,3];
            %--- labels
            this.axes.XAxis.Label.String='Generations';
            this.axes.XAxis.Label.FontWeight='bold';
            this.axes.YAxis.Label.String = 'f(B)';
            this.axes.YAxis.Label.FontWeight='bold';
            this.axes.XLim=[1,this.Ngen];
            this.axes.YLim=[-0.01,1.01];
            this.axes.XGrid='on';
            this.axes.YGrid='on';
            
            
            % --- add input controls
            glTP=matlab.ui.container.GridLayout('Parent',mainGridLayout,'RowSpacing',6,'ColumnSpacing',6,...
                'ColumnWidth',{'2x','1x'},'Scrollable','on','Padding',[0,0,0,0],...
                'RowHeight',{...
                this.uiHt,...                           % EVOLUTION CONTROLS header
                this.uiHt,this.uiHt,this.uiHt,...       % w0, w1, and w2 uicontrols
                this.uiHt,...                           % initial f(B) frequency
                this.uiHt,this.uiHt,this.uiHt...        % population size, number of generations, number of populations
                this.uiHt,...                           % spacer
                this.uiHt,...                           % STATS header
                2.7*this.uiHt,...                       % stats table
                '1x',...
                1.5*this.uiHt,...                       % RUN button
                });glTP.Layout.Row=[1,2];glTP.Layout.Column=1;
            
            % --- modelling parameters
            ir=1;
            hdummy=matlab.ui.container.Panel('Parent',glTP,'BorderType','none','Title','Model Parameters',...
                'FontWeight','bold','TitlePosition','centertop','ForegroundColor','#0072BD');
            hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=[1,2];
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','AA fitness W0','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w0,'Limits',[0,Inf],'UpperLimitInclusive',0,...
                'ValueChangedFcn',@(h,~)set(this,'W0',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','AB fitness W1','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w1,'Limits',[0,Inf],'UpperLimitInclusive',0,...
                'ValueChangedFcn',@(h,~)set(this,'W1',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','BB fitness W2','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w2,'Limits',[0,Inf],'UpperLimitInclusive',0,...
                'ValueChangedFcn',@(h,~)set(this,'W2',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','Initial q=f(B)','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.q_init,'Limits',[0,Inf],'UpperLimitInclusive',0,...
                'ValueChangedFcn',@(h,~)set(this,'q_init',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','Population size','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Nindiv,'Limits',[1,Inf],'LowerLimitInclusive',0,...
                'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Nindiv',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','No. generations','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Ngen,'Limits',[1,Inf],'LowerLimitInclusive',0,...
                'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Ngen',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            hdummy=matlab.ui.control.Label('Parent',glTP,'Text','No. populations','FontWeight','bold',...
                'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1;
            hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Npop,'Limits',[1,Inf],'LowerLimitInclusive',0,...
                'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Npop',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2;
            
            % --- output stats
            ir=ir+1;
            hdummy=matlab.ui.container.Panel('Parent',glTP,'BorderType','none','Title','Results',...
                'FontWeight','bold','TitlePosition','centertop','ForegroundColor','#0072BD');
            hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=[1,2];
            
            htable=matlab.ui.control.Table('Parent',glTP,'ColumnEditable',[false,false],'ColumnFormat',{'char','char'},'ColumnWidth',{'2x','1x'},'Data',cell.empty,...
                'RowName',{},'ColumnName',{}); addStyle(htable,uistyle('FontWeight','bold'),'column',1);
            htable.Layout.Row=ir; htable.Layout.Column=[1,2];
            
            hdummy=matlab.ui.control.Button('Parent',glTP,'Text','Run Evolution','FontWeight','bold','BackgroundColor','#EDB120',...
                'FontSize',this.uiHt,'Interruptible',false,'ButtonPushedFcn',@(~,~)run_buttonpushedfcn());
            hdummy.Layout.Row=numel(glTP.RowHeight);hdummy.Layout.Column=[1,2];
            function run_buttonpushedfcn()
                try
                    run(); plotData();
                catch ex
                    hdummy=uiconfirm(hf,{'CRITICAL ERROR,';' ';['Identifier: ',ex.identifier];['Message: ',ex.message]},...
                        this.Name,'Options',{'OK'},'Icon','error');
                end
            end
            
            %--- turn figure visibility on
            hf.Visible=1;
            
            %% --- supporting functions
            function run()
                try
                    hbar=matlab.ui.dialog.ProgressDialog(hf,'Indeterminate',0,'ShowPercentage',1,'Title',this.Name,'Message','Evolution in progress...');
                    wmax=max([this.w0,this.w1,this.w2]);
                    if(0==wmax); uialert(hf,'At least one fitness value must be non-zero.',this.Name); return; end
                    %--- calculated values
                    this.N_BB = zeros(this.Npop,this.Ngen);
                    this.N_AB = zeros(this.Npop,this.Ngen);
                    w0loc=this.w0/wmax; w1loc=this.w1/wmax; w2loc=this.w2/wmax;
                    Nindiv=this.Nindiv; Ngen=this.Ngen; Npop=this.Npop;
                    
                    for s = 1:Npop
                        hbar.Value=s/Npop;
                        q = this.q_init;
                        %--- initialize variables
                        this.N_BB(s,1) = round(Nindiv*q^2);
                        this.N_AB(s,1) = round(Nindiv*2*q*(1-q));
                        
                        %--- proceed with looping over generations
                        for m = 2:Ngen                            
                            %--- reset variables
                            g=1; ct=1;
                            genGenotype=nan(1,Nindiv);
                            
                            %--- begin looping until maximum capacity is reached (Npop)
                            while(g <= Nindiv)
                                %--- generate a genotype
                                donAlpha = rand(1) >= q;            % first allele; 1 represents presence of "A" allele, and 0 for "B" allele
                                donBeta = rand(1) >= q;             % second allele; 1 represents presence of "A" allele, and 0 for "B" allele
                                genotype = donAlpha+donBeta;        % genotype of the individual before fitness application
                                if(2==genotype)         % AA genotype
                                    if(1==w0loc || w0loc >= rand(1)); genGenotype(g)=2; g=g+1; end
                                elseif(1==genotype)     % AB genotype
                                    if(1==w1loc || w1loc >= rand(1)); genGenotype(g)=1; g=g+1; end
                                else                    % BB genotype
                                    if(1==w2loc || w2loc >= rand(1)); genGenotype(g)=0; g=g+1; end
                                end
                                %--- update loop counter
                                ct=ct+1;
                            end
                            this.N_BB(s,m)=nnz(0==genGenotype);
                            this.N_AB(s,m)=nnz(1==genGenotype);
                            %--- update q
                            q = (this.N_BB(s,m) + this.N_AB(s,m)/2)/Nindiv;        % frequency of "B" allele
                            %--- test if q went to 0 or 1, if so, break out of the loop and update the N_AB and N_BB matrices
                            if(0==q || 1==q)
                                this.N_BB(s,m+1:end)=this.N_BB(s,m);
                                this.N_AB(s,m+1:end)=this.N_AB(s,m);
                                break;
                            end
                        end
                    end
                    %--- update variables
                    qT=transpose((this.N_BB + this.N_AB/2)/Nindiv);
                    
                    %--- update table
                    htable.Data={...
                        'f(B)=1 fixation',nnz(1==max(qT,[],1));...
                        'f(B)=0 loss',nnz(0==min(qT,[],1))};
                    
                    %--- clean up
                    delete(hbar);
                catch ex
                    delete(hbar); rethrow(ex);
                end
            end
            
            function plotData()
                g = 1:this.Ngen;
                plot(this.axes,g,transpose((this.N_BB + this.N_AB/2)/this.Nindiv));
                this.axes.XLim=[1,this.Ngen];
            end
        end
    end
    
    %--- ICONS
    properties(Constant=true,Hidden=true)
        helpIcon=[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255;
            255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255;
            255,255,255,255,255,255,255,239,180,146,121,155,154,129,146,180,239,255,255,255,255,255,255,255;
            255,255,255,255,255,255,178,126,186,238,255,254,250,246,227,186,132,178,255,255,255,255,255,255;
            255,255,255,255,247,140,158,246,255,255,255,254,250,246,242,239,227,158,140,247,255,255,255,255;
            255,255,255,255,137,164,255,255,255,255,242,205,202,234,242,238,234,230,183,137,255,255,255,255;
            255,255,255,171,143,255,255,255,255,155,57,60,64,78,177,237,233,230,226,148,171,255,255,255;
            255,255,236,102,235,255,255,255,155,57,60,87,90,70,74,186,232,229,225,208,108,236,255,255;
            255,255,167,154,251,253,254,254,70,60,144,247,244,126,79,121,231,228,222,211,143,167,255,255;
            255,255,124,199,248,249,250,250,107,109,223,244,242,190,85,125,229,226,218,206,176,124,255,255;
            255,255,89,214,244,245,246,246,246,245,243,241,239,124,90,145,227,222,212,201,189,89,255,255;
            255,255,86,239,241,242,242,242,242,241,240,120,104,90,104,203,224,216,205,195,182,86,255,255;
            255,255,81,235,237,238,238,238,238,237,236,85,90,153,211,224,218,208,198,186,175,81,255,255;
            255,255,78,193,233,234,234,234,234,233,232,90,96,227,224,218,209,200,190,179,162,78,255,255;
            255,255,108,161,229,230,231,231,230,230,229,127,132,222,216,208,200,191,180,170,133,108,255,255;
            255,255,152,100,225,226,227,227,227,226,224,163,165,200,205,198,190,180,171,160,95,152,255,255;
            255,255,233,69,171,219,220,220,219,217,214,107,112,178,194,186,179,170,159,126,69,233,255,255;
            255,255,255,149,75,189,207,207,207,205,202,132,135,175,182,175,168,157,137,79,149,255,255,255;
            255,255,255,255,101,89,178,194,194,192,189,185,181,176,170,163,151,132,76,101,255,255,255,255;
            255,255,255,255,244,98,69,137,179,178,175,172,168,163,154,142,110,69,98,244,255,255,255,255;
            255,255,255,255,255,255,145,59,73,106,137,134,129,121,97,68,59,145,255,255,255,255,255,255;
            255,255,255,255,255,255,255,231,144,94,57,57,57,57,94,144,231,255,255,255,255,255,255,255;
            255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255;
            255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]
    end
end