% sweepFitDecayRate.m
% Jessica Parker, October 18, 2024
%
% This Matlab script takes the tail current data measured in measureTailCurrent.m and fits it with double exponential equation, and
% finally writes the parameters of best fit to an output file. It also plots the tail current data along with the best fit equation.
%
% NOTE! This code may take a few minutes to run. It is not the most efficient way to get the best fit, but it is effective.

close all;
clear all;

vrsn = 'A';
msrRealVrsn = 'A';

run1 = 4;
run2 = 0;
run3 = 0;
run4 = 2;
run5i = 1;
run5f = 3;
dataDir = 'data/';
figDir = 'plots/';

% Array of starting values for parameters to fit, with each element corresponding to a single run5 instance
strtingTau1s = [-0.8, -0.5, -0.25]; % Decay rate constant of single exponential or first part of double exponential
strtingTau2s = [-4.5, -3.5, -0.8]; % This parameter is only used if rtio is set less than 1. 2nd decay rate constant in double exponential
strtingRtios = [0.6, 0.65, 0.5]; % Proportion of double exponential belonging to tauA. Set to 1 if you want to fit with single 
                                 % exponential, set to < 1 for double exponential 

boltzType = 'h'; %Set to m or h depending on whether it is an activation or inactivation
percRange0 = 0.3; % Sets the range of parameter values you are testing as a proportion of the starting fit parameter value
nPercSteps = 200; % Sets the resolution (percRange/nPercSteps) you use to find the optimal parameter values
chngParThresh = 0.0005; %When the greatest proportion change in the parameters of best fit is less than this, the code stops searching
minRuns = 3; %The minimum number of searches that the sweep has to perform before settling on a parameter set

PlotResidualError = 1; % Set to 1 to plot the residual error of the last while loop run, set to 0 otherwise.
errVsblty = 'off'; %  % Set to 'off' to make residual error figures not pop up or set to 'on' otherwise. Final fit figures are always visible.
xmin = 0;
xmax = 15;
ymin = 90;
ymax = 250;
VariableYrange = 1;
PlotFit = 1;
CloseFigures = 0;
fntnm = 'Lucida Sans';
fntsz = 16;
fntwght = 'bold';
textYpostn = 0.6;
textXpostn = 0.2;
DullGreen = [155, 196, 114]/255;

dir2 = [num2str(run1) '_' num2str(run2) '_' num2str(run3)];
nruns = nPercSteps + 1;

tic;

for run5 = run5i:run5f
    tauA = strtingTau1s(run5);
    tauB = strtingTau2s(run5); 
    rtio = strtingRtios(run5); 

    dir3 = [dir2 '_' num2str(run4) '_' num2str(run5)];
    dlmwrite([dataDir 'startingValues' dir3 msrRealVrsn '_' vrsn '.txt'],[tauA, tauB, rtio]); % Save the starting fit parameter values used here

    current0 = load([dataDir 'tailCurrent' dir3 msrRealVrsn '.txt']); % Loading tail current
    tm = current0(2,:);
    dcay = -1*current0(1,:); % Flipping it upside down for simplicity's sake, since INaP is a negatively signed current.

    dcayStrt = dcay(1);
    dcayEnd = dcay(end);

    dcayFit = (dcayStrt-dcayEnd)*(rtio*exp(tm/tauA)+(1-rtio)*exp(tm/tauB))+dcayEnd; % Fit equation
    edif = (dcayFit - dcay).*(dcayFit - dcay); % Used to calculate residual error
    reserr = sum(edif); % (Residual error)^2

    chngPar = 1; % Just set to any number larger than any actual parameter proportion change calculated in the loop.
    bb = 1;
    sfty = 200; % Used for a safety break to prevent an infinite while loop
    while chngPar > chngParThresh && bb < sfty % Keeps looking for best fit until parameters stop changing or are barely changing

        if bb <= minRuns
            percRange = percRange0;
        elseif lastChngPar > chngPar && percRange > 0.1
            percRange = 0.97*percRange; % Range of parameter values tested narrows as it gets closer to best fit
        elseif percRange < 0.1 % Range width of parameters tested stops decreasing once it reaches 0.1*(current parameter value)
            percRange = 0.1;
        end
        percStep = percRange/nPercSteps; % Resolution of parameters values tested as a proportion of the starting parameter value

        %% Fitting 1rst exponential time constant, tauA
        tauAi = tauA-0.5*percRange*abs(tauA); % Start of tauA values to sweep
        parstep = percStep*abs(tauA); % Step size of tauA sweep
        lastPar = tauA; % tauA of best fit found in previous run of the loop
        tauAstep = parstep;

        for aa = 1:nruns
            tauA0 = tauAi + (aa-1)*parstep; % Sweeping tauA
            dcayFit0 = (dcayStrt-dcayEnd)*(rtio*exp(tm/tauA0)+(1-rtio)*exp(tm/tauB))+dcayEnd;  % Fit equation
            edif0 = (dcayFit0 - dcay).*(dcayFit0 - dcay); % used to calculate residual error
            reserrs(aa) = sum(edif0); % (residual error)^2
            aa = aa+1;
        end
        tauAReserrs = reserrs;
        [bf, bfIndx] = min(reserrs); % Finding best fit
        tauA = tauAi + (bfIndx-1)*parstep; % tauA of best fit
        tauAbf = bf;

        chngPar1 = abs(tauA - lastPar)/lastPar; % Proportion change of tauA from last run of the loop

        %% Fitting 2nd exponential time constant, tauB       
        if rtio < 1
            tauBi = tauB-0.5*percRange*abs(tauB); % Start of tauB values to sweep
            parstep = percStep*abs(tauB); % Step size of tauB sweep
            tauBstep = parstep;
            lastPar = tauB; % tauB of best fit found in previous run of the loop

            for aa = 1:nruns
                tauB0 = tauBi + (aa-1)*parstep; % Sweeping tauB
                dcayFit0 = (dcayStrt-dcayEnd)*(rtio*exp(tm/tauA)+(1-rtio)*exp(tm/tauB0))+dcayEnd;
                edif0 = (dcayFit0 - dcay).*(dcayFit0 - dcay);
                reserrs(aa) = sum(edif0);
                aa = aa+1;
            end
            tauBReserrs = reserrs;
            [bf, bfIndx] = min(reserrs);
            tauB = tauBi + (bfIndx-1)*parstep; % tauB of best fit
            tauBbf = bf;
            chngPar2 = abs(tauB - lastPar)/lastPar; % Proportion change of tauB from last run of the loop
        else
            chngPar2 = 0;
        end        

        %% Fitting double exponential ratio, rtio
        if rtio < 1
            rtioI = rtio-0.5*percRange*abs(rtio); % Start of rtio values to sweep
            parstep = percStep*abs(rtio); % Step size of rtio sweep
            rtiostep = parstep;
            lastPar = rtio; % rtio of best fit found in previous run of the loop

            for aa = 1:nruns
                rtio0 = rtioI + (aa-1)*parstep; % Sweeping rtio
                dcayFit0 = (dcayStrt-dcayEnd)*(rtio0*exp(tm/tauA)+(1-rtio0)*exp(tm/tauB))+dcayEnd;
                edif0 = (dcayFit0 - dcay).*(dcayFit0 - dcay);
                reserrs(aa) = sum(edif0);
                aa = aa+1;
            end
            rtioReserrs = reserrs;
            [bf, bfIndx] = min(reserrs);
            rtio = rtioI + (bfIndx-1)*parstep; % rtio of best fit
            rtiobf = bf;
            chngPar3 = abs(rtio - lastPar)/lastPar; % Proportion change of rtio from last run of the loop
        else
            rtio = 1;
            chngPar3 = 0;
        end        

        lastChngPar = chngPar;
        if bb >= minRuns % Forces while loop to run at least a minimum of minRuns times
            chngPar = max([chngPar1,chngPar2, chngPar3]); % Maximum proportion change of all fit parameters
        end

        bb = bb+1;

        %% Plot residual error of last fit parameter sweeps
        if chngPar < chngParThresh
            if PlotResidualError
                if rtio < 1
                    allpars = tauBi:tauBstep:tauBi+(nruns-1)*tauBstep;
                    f = figure('visible',errVsblty);
                    f.PaperUnits = 'inches';
                    f.Units = 'inches';
                    f.OuterPosition = [1 1 6.0 5.0];
                    f.InnerPosition = [0.25 0.25 5.5 4.5];
                    f.PaperPosition = [0.25 0.25 5.0 4.0];
                    f.RendererMode = 'manual';
                    axes('position',[0.15, 0.15, 0.8, 0.8]);
                    axes('position',[0.15, 0.15, 0.8, 0.8]);
                    plot(allpars,tauBReserrs);
                    hold on;
                    plot(tauB,tauBbf,'rx');
                    xlabel(['\tau_{' boltzType 'NaP2} (s)']);
                    ylabel('(Residual Error)^2');
                    print(f,[figDir 'tau2ResErr' dir3 msrRealVrsn '_' vrsn '.eps'],'-depsc','-r0');
                    if CloseFigures
                        close(f);
                    end

                    allpars = rtioI:rtiostep:rtioI+(nruns-1)*rtiostep;
                    f = figure('visible',errVsblty);
                    f.PaperUnits = 'inches';
                    f.Units = 'inches';
                    f.OuterPosition = [1 1 6.0 5.0];
                    f.InnerPosition = [0.25 0.25 5.5 4.5];
                    f.PaperPosition = [0.25 0.25 5.0 4.0];
                    f.RendererMode = 'manual';
                    axes('position',[0.15, 0.15, 0.8, 0.8]);
                    plot(allpars,rtioReserrs);
                    hold on;
                    plot(rtio,rtiobf,'rx');
                    xlabel(['Ratio']);
                    ylabel('(Residual Error)^2');
                    print(f,[figDir 'RatioResErr' dir3 msrRealVrsn '_' vrsn '.eps'],'-depsc','-r0');
                    if CloseFigures
                        close(f);
                    end
                end

                allpars = tauAi:tauAstep:tauAi+(nruns-1)*tauAstep;
                f = figure('visible',errVsblty);
                f.PaperPositionMode = 'manual';
                f.PaperUnits = 'inches';
                f.Units = 'inches';
                f.OuterPosition = [1 1 6.0 5.0];
                f.InnerPosition = [0.25 0.25 5.5 4.5];
                f.PaperPosition = [0.25 0.25 5.0 4.0];
                f.RendererMode = 'manual';
                axes('position',[0.15, 0.15, 0.8, 0.8]);
                plot(allpars,tauAReserrs);
                hold on;
                plot(tauA,tauAbf,'rx');
                xlabel(['\tau_{' boltzType 'NaP1} (s)']);
                ylabel('(Residual Error)^2');
                print(f,[figDir 'tau1ResErr' dir3 msrRealVrsn '_' vrsn '.eps'],'-depsc','-r0');
                if CloseFigures
                    close(f);
                end
            end
        end
        disp(['bb = ' num2str(bb) ', max par change = ' num2str(chngPar) ', tau1 = ' num2str(tauA) ', tau2 = ' ... 
            num2str(tauB) ', ratio = ' num2str(rtio) ', res. error = ' num2str(sqrt(bf))]);
    end
    disp(' ');

    %% Plot final fit
    if bb < sfty

        tmFit = xmin:0.01:xmax;
        dcayFit = (dcayStrt-dcayEnd)*(rtio*exp(tmFit/tauA)+(1-rtio)*exp(tmFit/tauB))+dcayEnd;

        tmPlt = tm;
        dcayPlt = dcay;

        if VariableYrange
            ymin0 = min(dcayFit);
            ymax0 = max(dcayFit);
            yl0 = ymax0 - ymin0;
            ymin = ymin0 - 0.1*yl0;
            ymax = ymax0 + 0.1*yl0;
        end
        yl = ymax - ymin;

        ymin0 = ymin;
        ymax0 = ymax;
        ymax = -1*ymin0;
        ymin = -1*ymax0;

        f = figure();
        f.PaperPositionMode = 'manual';
        f.PaperUnits = 'inches';
        f.Units = 'inches';
        f.OuterPosition = [1 1 6.0 5.4];
        f.InnerPosition = [0.25 0.25 5.5 4.9];
        f.PaperPosition = [0.25 0.25 5.0 4.4];
        f.RendererMode = 'manual';

        axes('position',[0.22 0.17 0.73 0.7]);
        hold on;
        plot(tmPlt,-1*dcayPlt,'k','linewidth',2);
        if PlotFit
            plot(tmFit,-1*dcayFit,'--','color',DullGreen,'linewidth',2);
        end
        text(xmin+textXpostn*(xmax-xmin),ymin+textYpostn*yl,['\tau_{' boltzType 'NaP} = ' num2str(round(abs(tauA),2)) ' s'],'fontsize',fntsz,'fontname',fntnm,'fontweight',fntwght);
        if rtio < 1
            text(xmin+(textXpostn+0.55)*(xmax-xmin),ymin+textYpostn*yl,[num2str(round(100*rtio)) ' %'],'fontsize',fntsz,'fontname',fntnm,'fontweight',fntwght);
            text(xmin+textXpostn*(xmax-xmin),ymin+(textYpostn-0.17)*yl,['\tau_{' boltzType 'NaP2} = ' num2str(round(abs(tauB),2)) ' s'],'fontsize',fntsz,'fontname',fntnm,'fontweight',fntwght);
            text(xmin+(textXpostn+0.55)*(xmax-xmin),ymin+(textYpostn-0.17)*yl,[num2str(round(100*(1-rtio))) ' %'],'fontsize',fntsz,'fontname',fntnm,'fontweight',fntwght);
        end
        ylabel(['I_{NaP} (pA)']);
        xlabel('Time (s)');
        xlim([xmin xmax]);
        ylim([ymin ymax]);
        box off;
        ax = gca;
        ax.FontName = fntnm;
        ax.FontSize = fntsz;
        ax.FontWeight = fntwght;
        ax.LineWidth = 3.0;

        print(f,[figDir boltzType 'DecayRateFit' dir3 msrRealVrsn '_' vrsn '.eps'],'-depsc','-r0');

        if CloseFigures
            close(f);
        end

        dlmwrite([dataDir boltzType 'DecayRate' dir3 msrRealVrsn '_' vrsn '.txt'], [tauA, tauB, rtio, sqrt(bf)]);

    else
        disp('Error! While loop reached safety break. Check starting parameter values and percRange value.');
    end
end

toc;