%% Step 5 of CNN Threshold Process:
%Classify Annotated ROIs.
drive=['F:' filesep 'data' filesep];%Change to fit appropriate VP input
cruise='TN368_VPRII';
path=[drive cruise filesep 'threshold'];
classifier='C:\plgui\ifcb-run\models\VPR_2020-01-10b\model.pt';%Change later
[~,classifiername,~]=fileparts(fileparts(classifier));%Separate out the name of the classifier
refdir='PA_diatom_test';%Reference ROI directory
inputpath=[path filesep refdir];
outputpath=[path filesep 'outcsv'];
outcsv=[classifiername '_' refdir];
cmdtxt=['! conda activate ifcb && cd C:\plgui\ifcb-run\ && python neuston_run.py ' inputpath ' --model ' classifier ' --outdir ' outputpath ' --outfile ' outcsv '{dir}.csv' ' --loaders 4' ' --batch-size 108'];
eval(cmdtxt)
%% Step 6 of CNN Threshold Process:
%Create plots in line with Brownlee et al, as well as general classification results.
qc=0:0.05:1;

results=dir([outputpath filesep outcsv '*']);
for i=1:size(results,1)
    temp=results(i).name;
    rownames(i,1)={temp(length([outcsv])+1:end-4)};
    result(i,1)={temp};clear temp;
end
tempfolders=dir(inputpath);
for i=3:size(tempfolders,1),folders(i-2)=str2double(tempfolders(i).name);end
clear tempfolders results;
for qcind=1:length(qc) %Quality check increments
    for i=1:size(result,1)
        data=importdata([outputpath filesep char(result(i,1))]);
        [num,maxind]=max((data.data),[],2);
        data=data.textdata;
        if qcind==1
            if i==1 %Initialize column names
                colnames=data(1,3:end)'; %List of all taxa names
                colnames(end+1,:)={'unknown'};%For those that fail quality control
                matrix=zeros(size(rownames,1),size(colnames,1),length(folders),length(qc));
            end
            taxasum(i,1)=size(data,1)-1;%Overall number of manually annotated ROIs for each taxon
        end
        for j=1:size(num,1)
            [temp,~,~]=fileparts(fileparts(char(data(j+1,1))));
            [~,temp,~]=fileparts(temp);
            foldind=find(folders==str2double(temp));
            if num(j)>qc(qcind) %Quality Control Check
                matrix(i,maxind(j),foldind,qcind)=matrix(i,maxind(j),foldind,qcind)+1;
            else %for unknown
                matrix(i,size(colnames,1),foldind,qcind)=matrix(i,size(colnames,1),foldind,qcind)+1;
            end
        end
    end
end
%Below calculates number of manually annotated ROIs in each taxa and folder.
for i=1:length(folders),rowsum(:,i)=sum(matrix(:,:,i,1),2);end
save([path filesep outcsv '.mat'],'rownames','colnames','folders','qc','taxasum','rowsum','matrix');
clear temp i j qcind foldind num maxind;
%% Make Plots
if ~exist('matrix','var'),load([path filesep outcsv '.mat']);end
set(groot,'defaultFigureColor',[1 1 1]);%Because I am sick of cleaning the graphs.
g=0;
for i=1:size(rownames,1)
    for j=1:size(colnames,1)
        if strcmp(char(rownames(i,1)),char(colnames(j,1)))
            taxon(:,:)=matrix(i,j,:,:);%Rows=Folders, Columns=qc
            if sum(sum(taxon))>2*length(folders)
                %Linear Regression Calculation
                for k=1:size(taxon,2)%Manual=x, CNN=y
                    x=[ones(length(rowsum(i,:)),1) rowsum(i,:)'];
                    b(1:2,k)=x\taxon(:,k);%Row 1 = Y-int, Row 2 = Slope
                    yCalc(:,k)=b(1,k)+b(2,k)*rowsum(i,:);
                    Rsq(k)=1-sum((taxon(:,k)-yCalc(:,k)).^2)/sum((taxon(:,k)-mean(taxon(:,k))).^2);
                end
                g=g+1;
                figure(g)
                subplot(221)
                [~,ind]=sort(rowsum(i,:));
                plot(rowsum(i,ind),taxon(ind,1:4:end))
                xlabel(['Manually Annotated ' char(rownames(i,1))])
                ylabel(['CNN Annotated ' char(rownames(i,1))])
                legend(string(qc(1:4:end)')','Location','eastoutside')
                title(['Manual vs. CNN Annotated ' char(rownames(i,1))])
                axis([0 inf 0 inf])
                subplot(222)
                plot(qc,Rsq)
                xlabel('CNN Annotation Threshold')
                ylabel('Linear Regression R^2')
                title(['R^2 ' char(rownames(i,1))])
                subplot(223)
                plot(qc,b(1,:))
                xlabel('CNN Annotation Threshold')
                ylabel('Linear Regression Y-Intercept')
                title(['Y-Intercept ' char(rownames(i,1))])
                subplot(224)
                plot(qc,b(2,:))
                xlabel('CNN Annotation Threshold')
                ylabel('Linear Regression Slope')
                title(['Slope ' char(rownames(i,1))])
            end
        end
    end
end