clear all 
close all
clc

T2check = input('1 = Validate Curve or 2 = Predict Unknown ');
SC = input('Use a previously acquired standard curve? (Y/N) ', 's');

if SC == 'Y'
    filename = uigetfile('*.xlsx', 'Choose standard curve');
    expfit = [xlsread(filename, 'a2:g2')];
    exp_fit.a = expfit(1);
    exp_fit.b = expfit(2);
    posTHRESH = expfit(3);
    HThresh = expfit(4);
    FAMthresh = expfit(5);
elseif SC == 'y'
    filename = uigetfile('*.xlsx', 'Choose standard curve');
    expfit = [xlsread(filename, 'a2:g2')];
    exp_fit.a = expfit(1);
    exp_fit.b = expfit(2);
    posTHRESH = expfit(3);
    HThresh = expfit(4);
    FAMthresh = expfit(5);
else
    tf = input('Build a new standard curve. Enter the number of training files? ');
    ts = input('Number of samples in training files ');
    tr = input('Number of replicates per sample in each experiment? ');
    f = zeros((ts-1),1);
    f(1)=input('What is the lowest template concentration (in log10 copies)? ');
    step = input('How many orders of magnitude are there between each standard? '); 
    HThresh = input('Slope threshold? ');  % HThresh = slope threshold
    z = input('Number (z) standad deviations above background for positive threshold? ');          
    JPGcheck = input('Is datatype from 1 = thermal cycler, 2 = microscope TIF, 3 = microscope JPG? ');


    if JPGcheck == 1
        FTchange = input('input new IPC threshold? y/n? ', 's' );
        if FTchange == 'y'
            FAMthresh = input('set IPC FAMthreshold: '); 
        else
            FAMthresh = 327.0745;
        end

    elseif JPGcheck == 2
       FTchange = input('input new IPC threshold? y/n? ', 's' );
        if FTchange == 'y'
            FAMthresh = input('set IPC FAMthreshold: '); 
        else
            FAMthresh = 327.0745;
        end
    elseif JPGcheck == 3
       FTchange = input('input new IPC threshold? y/n? ', 's' );
        if FTchange == 'y'
            FAMthresh = input('set IPC FAMthreshold: '); 
        else
            FAMthresh = 9;
        end
    end

    % READ IN TRAINING FILES

    HEX = zeros(60,ts*tr*tf);
    FAM = zeros(60,ts*tr*tf);

    for i =1:tf
       [tfilename, tpathname] = uigetfile('*.xlsx', 'Choose training file to build standard curve');
        cd(tpathname);
        HEX(1:60,(ts*tr*(i-1)+1):(ts*tr*i))= [xlsread(tfilename, 'HEX', 'c2:af61')];
        FAM(1:60,(ts*tr*(i-1)+1):(ts*tr*i)) = [xlsread(tfilename, 'FAM', 'c2:af61')];
    end

    % CHECK HEX NTCs TO SET THRESHOLD FOR SUCCESSFULLY RUN SAMPLES
    HEXtest = max(HEX)-min(HEX);
    NTCmean = mean(HEXtest(1:ts:end));
    NTCstd = std(HEXtest(1:ts:end));
    posTHRESH = NTCmean + z*NTCstd;

    %% TRAINING SECTION

    % check for failed samples (no HEX or FAM amplification
    FAMtest = max(FAM)-min(FAM);
    HEXclean = zeros(60, ts*tr*tf);

    HEXfail = zeros(1, ts*tr*tf);  %TESTCHECK tells indicies of negatives and invalid
    for i=1:ts*tr*tf
        if HEXtest(i)>posTHRESH;
           HEXfail(i)=HEXtest(i);
           HEXclean(:,i)=HEX(:,i);
        else HEXfail(:,i)=NaN;
            HEXclean(:,i)=NaN;
        end
    end

    FAMfail = zeros(1, ts*tr*tf);
    for i=1:ts*tr*tf
        if FAMtest(i)>JPGcheck;      %!!!!!NEED TO SET THIS VALUE APPROPRIATELY
           FAMfail(i)=FAMtest(i);
        else FAMfail(:,i)=NaN;       
        end
    end

    BOTHfail = isnan(HEXfail).*isnan(FAMfail);
    HEXpos = ~isnan(HEXfail);
    FAMpos = ~isnan(FAMfail);
    valid = ~BOTHfail;

    % CHANGE SAMPLES THAT DIDN'T REACH THRESHOLD TO NANs IN HEXclean
    HEXclean = zeros(60, ts*tr*tf);
    for i=1:ts*tr*tf    
        if HEXtest(i)>posTHRESH;
           HEXclean(:,i)=HEX(:,i);
        else HEXclean(:,i)=NaN;   
        end
    end

    % smooth HEXclean
    HEXSMOOTH = {zeros(60,ts*tr*tf)};
    for i=1:ts*tr*tf                        
        HEXCOL = HEXclean(:,i);
        HEXSMOOTH{i} = smooth(HEXCOL,5);
    end
    HEXSMOOTH = cell2mat(HEXSMOOTH);

    % claculates rough first deriv
    HDIFF = {zeros(59,ts*tr*tf)};
    for i=1:ts*tr*tf                        
        HEXDIFFCOL = HEXSMOOTH(:,i);
        HDIFF{i} = diff(HEXDIFFCOL);
    end
    HDIFF = cell2mat(HDIFF);
    % cleans matrix back to 60 rows
    firstrow = zeros(1, ts*tr*tf);          
    HDIFF = [firstrow; HDIFF];

    % find threshold cycle
    HTHRESHcycle = {zeros(1,ts*tr*tf)};
    for i = 1:ts*tr*tf
        HTHRESHcycleCOL = HDIFF(:,i);
        HTHRESHcycle{i} = find(HTHRESHcycleCOL > HThresh, 1, 'first');
    end
    HTHRESHcycleCLEAN = cellfun('isempty', HTHRESHcycle);
    % for unamplified samples, HTHRESHcycle turned to NaN.
    HTHRESHcycle(HTHRESHcycleCLEAN)={nan};      
    HTHRESHcycle=cell2mat(HTHRESHcycle);

    % take every 6th value starting with 1, 2, etc to get a matrix with all
    % threshold values corresponding to a particular input log copy #
    HTHRESH0  = HTHRESHcycle(1:ts:end); 
    HY = zeros(1, tr*tf*(ts-1));
    for i = 1:(ts-1)
        HY((tr*tf*(i-1)+1):tr*tf*i) = HTHRESHcycle((i+1):ts:end);
    end

    % build an array with X values for each of the HY values

    for i = 2:(ts-1)
        f(i) =f(i-1)+step;
    end

    x=zeros((ts-1), tr*tf);
    for i=1:tr*tf                 
        x(:,i) = f;
    end

    xarray = zeros(1, (ts-1)*tr*tf);
    for i = 1:(ts-1)
        xarray(((i-1)*tr*tf+1):(i*tr*tf)) = x(i,:);
    end

    % clean the vector to get rid of x spots corresponding to HY NaNs
    xclean = ~isnan(HY).*xarray;
    xclean(xclean == 0) = [];
    xfit = [0:0.1:5];
    % get rid of NaN in HY
    HY(find(isnan(HY)))=[];

    % BEST FIT: EXPONENTIAL
    exp_fit = fit(xclean', HY', 'exp1');
    % use fit to calculate expected values
    ypred = exp_fit.a*exp(exp_fit.b*xclean);
    % use expected values to calculate R-squared
    SSE = sum((HY-ypred).^2);
    SST = sum((HY-mean(HY)).^2); 
    R2_EXPONENTIAL = 1-SSE/SST;
end


%% VALIDATION SECTION
vf = input('Number of files with validation/unknown samples? ');
vs = input('Number of validation/unknown samples per file? ');
vr= input('Number of replicates per sample in each file? ');
vs0 = zeros(vs, 1);
vs0(2) = expfit(6);
vstep = expfit(7);
    for i = 3:vs;
        vs0(i) = vs0(i-1)+vstep;
    end

vHEX = zeros(60,vs*vr*vf);
vFAM = zeros(60,vs*vr*vf);
sample_name = cell(1, vs*vr*vf);

for i =1:vf
   [vfilename, vpathname] = uigetfile('*.xlsx', 'Choose files with unknown sample data');
    cd(vpathname);
    vHEX(1:60,(vs*vr*(i-1)+1):(vs*vr*i))= [xlsread(vfilename, 'HEX', 'c2:af61')]; 
    vFAM(1:60,(vs*vr*(i-1)+1):(vs*vr*i))= [xlsread(vfilename, 'FAM', 'c2:af61')];
    [NUMERIC,TXT,RAW]=XLSREAD(vfilename, 'HEX', 'c1:af1');
    sample_name((vs*vr*(i-1)+1):vs*vr*i) = RAW(1:vs*vr);
end

% check for failed samples (no HEX or FAM amplification
vHEXtest = max(vHEX)-min(vHEX);
vFAMtest = max(vFAM)-min(vFAM);
vHEXclean = zeros(60, vs*vr*vf);

vHEXfail = zeros(1, vs*vr*vf);  %TESTCHECK tells indicies of negatives and invalid
for i=1:vs*vr*vf
    if vHEXtest(i)>posTHRESH;
       vHEXfail(i)=vHEXtest(i);
       vHEXclean(:,i)=vHEX(:,i);
    else vHEXfail(:,i)=NaN;
        vHEXclean(:,i)=NaN;
    end
end

vFAMfail = zeros(1, vs*vr*vf);
for i=1:vs*vr*vf
    if vFAMtest(i)>327.0745;         %from file FAM_Background_Threshold
       vFAMfail(i)=vFAMtest(i);
    else vFAMfail(:,i)=NaN;       
    end
end

vBOTHfail = isnan(vHEXfail).*isnan(vFAMfail);

% smooth hex values, calc 1st deriv, calculate threshold cycle
vHEXSMOOTH = {zeros(60,vs*vr*vf)};
for i=1:vs*vr*vf
    vHEXCOL = vHEXclean(:,i);
    vHEXSMOOTH{i} = smooth(vHEXCOL,5);
end
vHEXSMOOTH = cell2mat(vHEXSMOOTH);

vHDIFF = {zeros(59,vs*vr*vf)};
for i=1:vs*vr*vf
    vHEXDIFFCOL = vHEXSMOOTH(:,i);
    vHDIFF{i} = diff(vHEXDIFFCOL);
end
vHDIFF = cell2mat(vHDIFF);
firstrow = zeros(1, vs*vr*vf);
vHDIFF = [firstrow; vHDIFF];

vTHRESHcycle = {zeros(1,vs*vr*vf)};
for i = 1:vs*vr*vf
    vTHRESHcycleCOL = vHDIFF(:,i);
    vTHRESHcycle{i} = find(vTHRESHcycleCOL > HThresh, 1, 'first');
end
vTHRESHcycleCLEAN = cellfun('isempty', vTHRESHcycle);
vTHRESHcycle(vTHRESHcycleCLEAN)={NaN};
vTHRESHcycle=cell2mat(vTHRESHcycle);  

vTHRESH = zeros(vs, vf*vr*(vs-1));
vTHRESH0  = vTHRESHcycle(1:vs:end);
for i=1:(vs-1)
    vTHRESH(i, (vf*vr*(i-1)+1):vf*vr*i) = vTHRESHcycle((i+1):vs:end);
end

VPM = zeros((vs-1), vr*vf);
for i = 1:(vs-1)
    VPM(i, 1:vr*vf) = vTHRESHcycle((i+1):vs:end);
end


% use fit to back caculate x values for the validation threshold (y) values
vpx = (log(VPM/exp_fit.a))/exp_fit.b;
NANvpx = vpx;
vpx(isnan(vpx))=0;
vpzeros = (log(vTHRESH0/exp_fit.a))/exp_fit.b;

c = {zeros(1,vs*vr*vf)};
for i=1:vs*vr*vf
    if (isnan(vTHRESHcycle(i)))
        if(vBOTHfail(i)==1)
               c{i} = 'Invalid Run';
        else
            c{i} = 'Negative for HIV';
        end
    else
        c{i} = log(vTHRESHcycle(i)/exp_fit.a)/exp_fit.b;
    end             
end
Predicted_log_copies = c';
sample_name = sample_name';
Table1 = dataset(sample_name, Predicted_log_copies)

NANLCNO = bsxfun(@minus, NANvpx, vs0(2:vs)); %log copy number off expected value
d=[vpzeros; vpx];
d(isnan(d))=0;

avgPredicted = nanmean(d, 2);

if vr>1,
    stdOffExpected = (nanstd(d'))';
elseif vf>1,
    stdOffExpected = (nanstd(d'))';
else
    stdOffExpected = zeros(vs, 1);
end

VPP = zeros(vs, 1); % Validation Percent Positive
VPP(1) = sum(~isnan(vpzeros))/(vr*vf);
for i = 1:(vs-1)
    VPP(i+1) = sum(~isnan(NANLCNO((i), :)))/(vr*vf);
end

TrueConc = vs0;              
Percent_Pos = VPP;
Avg_Pred_Conc = avgPredicted;
StDev = stdOffExpected;
Pred_Conc = d;

if T2check == 1
    Table2 = dataset(TrueConc, Percent_Pos, Avg_Pred_Conc, StDev)
else 
end