%%%%Brain-guided network training code, set up for parallel computing
%%%%Author: Bruce C Hansen, Colgate University, Neuroscience Program, Department of Psychological & Brain Sciences

clear all

NumStim = 80;
ImDim = 128;

%%%%%%Get Network Input%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('Path to StimImages.mat');
ImagePath = 'Path to stimuli';

ValidStim = [];
count = 0;
for loop = 1:NumStim
    if(StimImages(loop,2) == 1)
        count = count+1;
        ValidStim(count,1) = loop;
        ValidStim(count,2) = StimImages(loop,1);
    end
end
TotalTargets = size(ValidStim,1);
        

TAREGT_RMS = 0.2;
for loop = 1:TotalTargets
    I = imread(strcat(ImagePath,num2str(ValidStim(loop,2)),'.png'));
    I = imresize(I,[ImDim ImDim],'bicubic');
    I = double(I);
    I_Mu = mean2(I);
    I_SD = std2(I);
    I_RMS = I_SD/I_Mu;
    I_NORM = I-I_Mu;
    I_NORM_Scaled = I_NORM./max(max(abs(I_NORM)))./2;
    RMS_altered = std2(I_NORM_Scaled+.5)./mean2(I_NORM_Scaled+.5);
    rms_alt_scale = (TAREGT_RMS*2)/RMS_altered;
    I_ADJST = I_NORM_Scaled.*rms_alt_scale;
    I_ADJST = (I_ADJST*127);
    NetworkInput(:,loop) = reshape(I_ADJST,ImDim*ImDim,1);
end
NetworkInput = NetworkInput.*0.01;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%Get Brain Supervision%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

JunkEs = [127,17,126,125,128,48,43,44,39,49,45,50,57,56,119,120,114,113,115,108,101,100,107];
FullNet = 1:128;
ValidEs = setdiff(FullNet,JunkEs)';
load('Path to TaskCorrs_ValidStim_ERPs.mat');
TaskR2s = 1-(TaskCorrs.^2);
R2_Thresh = 0.70;
SelectedEs = {};

for TimeStep = 1:494
    SelectedEs{TimeStep,1} = ValidEs(find(TaskR2s(ValidEs,TimeStep) > R2_Thresh));
end

DataFile = 'pAVG_n24'; %<---training using the participant averaged data


%%%%Time 1%%%%%%%%%%%%%%%%%%%
TimeStep = 164;
HalfWin = 5;
StartTime = TimeStep-HalfWin;
StopTime = TimeStep+HalfWin;

Electrodes = SelectedEs{TimeStep-100,1};
DataSelectAll = [];
load(strcat('Path to data',DataFile,'_Otask_TrialAvgData.mat'));
DataSelect1 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect1,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
ObjDists(:,1) = DataSelectVect;
ObjDat1(:,:) = mean(VEPdat,3);

load(strcat('Path to data',DataFile,'_Ftask_TrialAvgData.mat'));
DataSelect2 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect2,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
FuncDists(:,1) = DataSelectVect;
FuncDat1(:,:) = mean(VEPdat,3);

%%%%Time 2%%%%%%%%%%%%%%%%%%%
TimeStep = 259;
HalfWin = 5;
StartTime = TimeStep-HalfWin;
StopTime = TimeStep+HalfWin;

Electrodes = SelectedEs{TimeStep-100,1};
DataSelectAll = [];
load(strcat('Path to data',DataFile,'_Otask_TrialAvgData.mat'));
DataSelect1 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect1,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
ObjDists(:,2) = DataSelectVect;
ObjDat2(:,:) = mean(VEPdat,3);

load(strcat('Path to data',DataFile,'_Ftask_TrialAvgData.mat'));
DataSelect2 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect2,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
FuncDists(:,2) = DataSelectVect;
FuncDat2(:,:) = mean(VEPdat,3);

%%%%Time 3%%%%%%%%%%%%%%%%%%%
TimeStep = 344;
HalfWin = 5;
StartTime = TimeStep-HalfWin;
StopTime = TimeStep+HalfWin;

Electrodes = SelectedEs{TimeStep-100,1};
DataSelectAll = [];
load(strcat('Path to data',DataFile,'_Otask_TrialAvgData.mat'));
DataSelect1 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect1,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
ObjDists(:,3) = DataSelectVect;
ObjDat3(:,:) = mean(VEPdat,3);

load(strcat('Path to data',DataFile,'_Ftask_TrialAvgData.mat'));
DataSelect2 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect2,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
FuncDists(:,3) = DataSelectVect;
FuncDat3(:,:) = mean(VEPdat,3);

%%%%Time 4%%%%%%%%%%%%%%%%%%%
TimeStep = 502;
HalfWin = 5;
StartTime = TimeStep-HalfWin;
StopTime = TimeStep+HalfWin;

Electrodes = SelectedEs{TimeStep-100,1};
DataSelectAll = [];
load(strcat('Path to data',DataFile,'_Otask_TrialAvgData.mat'));
DataSelect1 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect1,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
ObjDists(:,4) = DataSelectVect;
ObjDat3(:,:) = mean(VEPdat,3);

load(strcat('Path to data',DataFile,'_Ftask_TrialAvgData.mat'));
DataSelect2 = TrialAvgData(Electrodes,StartTime:StopTime,:); %<----set up for trial averaged data
holdVEPs = TrialAvgData(Electrodes,:,:);
VEPdat(:,:) = mean(holdVEPs,1);
DataSelectAll(:,:) = mean(DataSelect2,1);
DataSelectVect = mean(DataSelectAll,1)';
DataSelectVect = DataSelectVect./(max(abs(DataSelectVect)));
FuncDists(:,4) = DataSelectVect;
FuncDat3(:,:) = mean(VEPdat,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%Define Target Responses%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NumSims = 1000;

IterNum = 64.*2;
IterHlf = IterNum./2;

SimTargets = [];

for loop = 1:TotalTargets
    
    SimTargets(loop,1) = loop;
    SimTargets(loop+IterHlf,1) = loop;
    
    SimTargets(loop,2) = ValidStim(loop,1);
    SimTargets(loop+IterHlf,2) = ValidStim(loop,1);
    
    StimTargets(loop,3) = 0;
    SimTargets(loop+IterHlf,3) = 1;
    
    SimTargets(loop,4) = ObjDists(ValidStim(loop,1),1);
    SimTargets(loop+IterHlf,4) = FuncDists(ValidStim(loop,1),1);
    
    SimTargets(loop,5) = ObjDists(ValidStim(loop,1),2);
    SimTargets(loop+IterHlf,5) = FuncDists(ValidStim(loop,1),2);
    
    SimTargets(loop,6) = ObjDists(ValidStim(loop,1),3);
    SimTargets(loop+IterHlf,6) = FuncDists(ValidStim(loop,1),3);
    
    SimTargets(loop,7) = ObjDists(ValidStim(loop,1),4);
    SimTargets(loop+IterHlf,7) = FuncDists(ValidStim(loop,1),4);

end

SimTargets(:,8) = 1:IterNum;
    

BrainRespMat(:,1:4) = SimTargets(1:IterHlf,4:7);%%%Obj t1 t2 t3 t4
BrainRespMat(:,5:8) = SimTargets(IterHlf+1:IterNum,4:7);%%%func t1 t2 t3 t4
BrainCorrSamp = corrcoef(BrainRespMat)


%Whiten neural data (comment out if raw data desired)
epsilon = 0.0001;
ZCAdat = ZCAwhiten(BrainRespMat,epsilon);
ZCACorrSamp = corrcoef(ZCAdat)

for col = 1:4
    for row = 1:4
        ZCACorrSamp(col,row) = corr(BrainRespMat(:,col),ZCAdat(:,row));
    end
end

ZCACorrSamp

SimTargets(1:IterHlf,4:7) = ZCAdat(:,1:4);
SimTargets(IterHlf+1:IterNum,4:7) = ZCAdat(:,5:8);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Sig_b = 1;
Sig_Hmax = 0;
LreLu = .1;

%%%%%%Initialize Network Weights%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MXpoolDim = 3;
MXpoolHlf = floor(MXpoolDim./2);

H1_Nodes = 100;
H1stride = 2;
H1mxStride = 2;
H1dim = 7;
H1dimhlf = floor(H1dim./2);
H1Scalar = sqrt(2./H1_Nodes);
H1w = (rand(H1dim,H1dim,H1_Nodes).*(H1Scalar.*2))-H1Scalar;
H1RespDim = ImDim./H1stride;
H1mxDim = H1RespDim./H1mxStride;

SubSampDim = length(1:H1stride:ImDim);
mxpDim1 = floor(SubSampDim./H1mxStride);

H2_Nodes = 170;
H2stride = 1;   
H2mxStride = 1;
H2dim = 5;
H2dimhlf = floor(H2dim./2);
H2Scalar = sqrt(2./H2_Nodes);
H2w = (rand(H2dim,H2dim,H1_Nodes,H2_Nodes).*(H2Scalar.*2))-H2Scalar;
mxpDim2 = floor(mxpDim1./H2mxStride);
H2RespDim = H1RespDim./H1mxStride./H2stride;

H3_Nodes = 170;
H3stride = 1;   
H3mxStride = 1;
H3dim = 3;
H3dimhlf = floor(H3dim./2);
H3Scalar = sqrt(2./H3_Nodes);
H3w = (rand(H3dim,H3dim,H2_Nodes,H3_Nodes).*(H3Scalar.*2))-H3Scalar;
mxpDim3 = floor(mxpDim2./H3mxStride);
H3RespDim = H2RespDim./H2mxStride./H3stride;

H4_Nodes = 150;
H4stride = 1;   
H4mxStride = 1;
H4dim = 3;
H4dimhlf = floor(H4dim./2);
H4Scalar = sqrt(2./H4_Nodes);
H4w = (rand(H4dim,H4dim,H3_Nodes,H4_Nodes).*(H4Scalar.*2))-H4Scalar;
mxpDim4 = floor(mxpDim3./H4mxStride);
H4RespDim = H3RespDim./H3mxStride./H4stride;

FC1_Nodes = 280;
FC1Scalar = sqrt(1./FC1_Nodes);
FC1w = (rand(FC1_Nodes,mxpDim4*mxpDim4*H4_Nodes).*(FC1Scalar.*2))-FC1Scalar;
OutW = (rand(FC1_Nodes,1).*(FC1Scalar.*2))-FC1Scalar;


Lrate = 0.0065; 
Lrate_H4 = 0.0015;
Lrate_H3 = 0.0015;
Lrate_H2 = 0.0055;
Lrate_H1 = 0.0035;


f = figure;
set(gcf,'color','w');
f.Position = [200 500 1800 1200];

%Begin training

for SimNum = 1:NumSims
    
    
    ESUM = 0;
    RESSUM1 = 0;
    RESSUM2 = 0;
    RESSUM3 = 0;
    RESSUM4 = 0;
    
    IterInfo = SimTargets;
    R = randperm(IterNum)';
    IterInfo(:,9) = R;
    IterInfo = sortrows(IterInfo,9);
    IterInfo = IterInfo(:,1:8);
    
    TrackCorr = [];
    
    for IterLoop = 1:IterNum
        
        Stim = IterInfo(IterLoop,1);
        TargetResp = IterInfo(IterLoop,3);
        BrainResp_t1 = IterInfo(IterLoop,4);
        BrainResp_t2 = IterInfo(IterLoop,5);
        BrainResp_t3 = IterInfo(IterLoop,6);
        BrainResp_t4 = IterInfo(IterLoop,7);
        
        %%%Layer 1 (passive)
        L1 = reshape(NetworkInput(:,Stim),ImDim,ImDim);
        
        %%%Tracking matrices
        RowColOut = cell(H1_Nodes,1);
        RowColIm = cell(H1_Nodes,1);
        RowColOut2 = cell(H2_Nodes,1);
        RowColIm2 = cell(H2_Nodes,1);
        RowColOut3 = cell(H3_Nodes,1);
        RowColIm3 = cell(H3_Nodes,1);
        RowColOut4 = cell(H4_Nodes,1);
        RowColIm4 = cell(H4_Nodes,1);
        
        RespResMap_h1 = cell(H1_Nodes,1);
        RespResMap_h2 = cell(H2_Nodes,1);
        RespResMap_h2(:,:) = {zeros(H2RespDim,H2RespDim)};
        RespResMap_h3 = cell(H3_Nodes,1);
        RespResMap_h3(:,:) = {zeros(H3RespDim,H3RespDim)};
        RespResMap_h4 = cell(H4_Nodes,1);
        RespResMap_h4(:,:) = {zeros(H4RespDim,H4RespDim)};
        
        H1poolsize = length(1:H1mxStride:length(1:H1stride:ImDim));
        H1pool = cell(H1_Nodes,1);
        H1RESpool = cell(H1_Nodes,1);
        
        
        %%%Hidden Layer 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        parfor node = 1:H1_Nodes
            
            RespMap = [];
            
            H1resp = conv2(L1,H1w(:,:,node),'same');
            NormResp = H1resp-mean2(H1resp);
            NormResp = NormResp./std2(H1resp);
            
            for y = 1:H1stride:ImDim
                for x = 1:H1stride:ImDim
                    Resp = H1resp(y,x);
                    RespMap(y,x) = Resp;
                    RespN = NormResp(y,x);
                    Residual = BrainResp_t1-RespN;
                    RespResMap_h1{node}(y,x) = Residual;
                end
            end

            
            tmpdrop = nonzeros(RespMap);
            RespMap = reshape(tmpdrop,H1RespDim,H1RespDim);
            
            tmpdrop = nonzeros(RespResMap_h1{node}(:,:));
            RespResMap_h1{node} = reshape(tmpdrop,H1RespDim,H1RespDim);
            
            
            H1size = size(RespMap,1);
            tmp1 = padarray(RespMap,[MXpoolHlf MXpoolHlf],0,'both'); %%%pad with nan to avoid zero pad being minimum
            tmp2 = padarray(RespResMap_h1{node},[MXpoolHlf MXpoolHlf],0./0,'both'); %%%pad with nan to avoid zero pad being minimum
            
            
            for y = 1:H1mxStride:H1size
                for x = 1:H1mxStride:H1size
                    kern1 = tmp1(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter response
                    kern2 = tmp2(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter residual
                    
                    [Val,IDX] = nanmin(abs(kern2(:)));%%%get abs min error
                    [Y X] = ind2sub(size(kern2),IDX);
                    H1pool{node}(y,x) = kern1(Y,X);%%%get filter response that corresponds to min error
                    H1RESpool{node}(y,x) = kern2(Y,X);

                end
            end
            
            
            tmpdrop = nonzeros(H1pool{node});
            H1pool{node} = reshape(tmpdrop,H1mxDim,H1mxDim);
            tmpdrop = nonzeros(H1RESpool{node});
            H1RESpool{node} = reshape(tmpdrop,H1mxDim,H1mxDim);
            
            H1Out(:,:,node) = max(H1pool{node}.*LreLu,H1pool{node}); %%%leaky ReLu
            [RowColOut{node,1}(:,1) RowColOut{node,1}(:,2)] = find(H1Out(:,:,node)~=0);
            
            count = 0;
            for y = 1:H1mxDim
                for x = 1:H1mxDim
                    count = count+1;
                    tmp = H1pool{node}(y,x);
                    [RowColIm{node,1}(count,1) RowColIm{node,1}(count,2)] = find(H1resp == tmp,1);
                end
            end
            
        end

        
        %%%%Hidden Layer 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        H1Outsize = size(H1Out(:,:,:),1);
        tmp2 = padarray(H1Out,[H2dimhlf H2dimhlf],0,'both');
        
        H2poolsize = length(1:H2mxStride:H1Outsize);
        H2pool = cell(H2_Nodes,1);
        H2pool(:,:) = {zeros(H2poolsize,H2poolsize)};
        H2RESpool = cell(H2_Nodes,1);
        H2RESpool(:,:) = {zeros(H2poolsize,H2poolsize)};
        

        parfor node = 1:H2_Nodes
            
            H2resp = [];
            
            Vfilt = reshape(H2w(:,:,:,node),H2dim*H2dim*H1_Nodes,1);
            tmpresp = [];
            for y = 1:H2stride:H1Outsize
                for x = 1:H2stride:H1Outsize
                    Vpatch = reshape(tmp2(y:y+H2dim-1,x:x+H2dim-1,:),H2dim*H2dim*H1_Nodes,1);
                    H2resp(y,x) = sum(Vfilt.*Vpatch);
                end
            end
            
            %%%%normalize for brainresp comparison
            H2resp = H2resp-mean2(H2resp);
            H2resp = H2resp./std2(H2resp);
            H2Outsize = size(H2resp,1);
            
            
            %%%%Get layer 2 residuals
            for y = 1:H2Outsize
                for x = 1:H2Outsize
                    RespN = H2resp(y,x);
                    Residual = BrainResp_t2-RespN;
                    RespResMap_h2{node}(y,x) = Residual;
                end
            end
            
            
            tmp21 = padarray(H2resp,[MXpoolHlf MXpoolHlf],0,'both');
            tmp22 = padarray(RespResMap_h2{node},[MXpoolHlf MXpoolHlf],0./0,'both');
            
            for y = 1:H2mxStride:H2Outsize
                for x = 1:H2mxStride:H2Outsize
                    kern1 = tmp21(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter response
                    kern2 = tmp22(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter residual
                    
                    [Val,IDX] = nanmin(abs(kern2(:)));%%%get abs min error
                    [Y X] = ind2sub(size(kern2),IDX);
                    H2pool{node}(y,x) = kern1(Y,X);
                    H2RESpool{node}(y,x) = kern2(Y,X);
                end
            end
            H2Out(:,:,node) = max(H2pool{node}.*LreLu,H2pool{node});
            [RowColOut2{node,1}(:,1) RowColOut2{node,1}(:,2)] = find(H2Out(:,:,node)~=0);
            
            count = 0;
            for y = 1:H2RespDim
                for x = 1:H2RespDim
                    count = count+1;
                    tmp = H2pool{node}(y,x);
                    [RowColIm2{node,1}(count,1) RowColIm2{node,1}(count,2)] = find(H2resp == tmp,1);
                end
            end

        end
    
        
        %%%%Hidden Layer 3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        H2Outsize = size(H2Out(:,:,:),1);
        tmp2 = padarray(H2Out,[H2dimhlf H2dimhlf],0,'both');
        
        H3poolsize = length(1:H3mxStride:H2Outsize);
        H3pool = cell(H3_Nodes,1);
        H3pool(:,:) = {zeros(H3poolsize,H3poolsize)};
        H3RESpool = cell(H3_Nodes,1);
        H3RESpool(:,:) = {zeros(H3poolsize,H3poolsize)};
        

        parfor node = 1:H3_Nodes
            
            H3resp = [];
            
            Vfilt = reshape(H3w(:,:,:,node),H3dim*H3dim*H2_Nodes,1);
            tmpresp = [];
            for y = 1:H3stride:H2Outsize
                for x = 1:H3stride:H2Outsize
                    Vpatch = reshape(tmp2(y:y+H3dim-1,x:x+H3dim-1,:),H3dim*H3dim*H2_Nodes,1);
                    H3resp(y,x) = sum(Vfilt.*Vpatch);
                end
            end
            
            %%%%normalize for brainresp comparison
            H3resp = H3resp-mean2(H3resp);
            H3resp = H3resp./std2(H3resp);
            H3Outsize = size(H3resp,1);
            
            
            %%%%Get layer 3 residuals
            for y = 1:H3Outsize
                for x = 1:H3Outsize
                    RespN = H3resp(y,x);
                    Residual = BrainResp_t3-RespN;
                    RespResMap_h3{node}(y,x) = Residual;
                end
            end
            
            
            tmp21 = padarray(H3resp,[MXpoolHlf MXpoolHlf],0,'both');
            tmp22 = padarray(RespResMap_h3{node},[MXpoolHlf MXpoolHlf],0./0,'both');
            
            for y = 1:H3mxStride:H3Outsize
                for x = 1:H3mxStride:H3Outsize
                    kern1 = tmp21(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter response
                    kern2 = tmp22(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter residual
                    
                    [Val,IDX] = nanmin(abs(kern2(:)));%%%get abs min error
                    [Y X] = ind2sub(size(kern2),IDX);
                    H3pool{node}(y,x) = kern1(Y,X);
                    H3RESpool{node}(y,x) = kern2(Y,X);
                end
            end
            H3Out(:,:,node) = max(H3pool{node}.*LreLu,H3pool{node});
            [RowColOut3{node,1}(:,1) RowColOut3{node,1}(:,2)] = find(H3Out(:,:,node)~=0);
            
            count = 0;
            for y = 1:H3RespDim
                for x = 1:H3RespDim
                    count = count+1;
                    tmp = H3pool{node}(y,x);
                    [RowColIm3{node,1}(count,1) RowColIm3{node,1}(count,2)] = find(H3resp == tmp,1);
                end
            end

        end
        
        
        %%%%Hidden Layer 4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        H3Outsize = size(H3Out(:,:,:),1);
        tmp2 = padarray(H3Out,[H3dimhlf H3dimhlf],0,'both');
        
        H4poolsize = length(1:H4mxStride:H3Outsize);
        H4pool = cell(H4_Nodes,1);
        H4pool(:,:) = {zeros(H4poolsize,H4poolsize)};
        H4RESpool = cell(H4_Nodes,1);
        H4RESpool(:,:) = {zeros(H4poolsize,H4poolsize)};
        
        parfor node = 1:H4_Nodes
            
            H4resp = [];
            
            Vfilt = reshape(H4w(:,:,:,node),H4dim*H4dim*H3_Nodes,1);
            tmpresp = [];
            for y = 1:H4stride:H3Outsize
                for x = 1:H4stride:H3Outsize
                    Vpatch = reshape(tmp2(y:y+H4dim-1,x:x+H4dim-1,:),H4dim*H4dim*H3_Nodes,1);
                    H4resp(y,x) = sum(Vfilt.*Vpatch);
                end
            end
            
            %%%%normalize for brainresp comparison
            H4resp = H4resp-mean2(H4resp);
            H4resp = H4resp./std2(H4resp);
            H4Outsize = size(H4resp,1);
            
            
            %%%%Get layer 4 residuals
            for y = 1:H4Outsize
                for x = 1:H4Outsize
                    RespN = H4resp(y,x);
                    Residual = BrainResp_t4-RespN;
                    RespResMap_h4{node}(y,x) = Residual;
                end
            end
            
            
            tmp21 = padarray(H4resp,[MXpoolHlf MXpoolHlf],0,'both');
            tmp22 = padarray(RespResMap_h4{node},[MXpoolHlf MXpoolHlf],0./0,'both');
            
            for y = 1:H4mxStride:H4Outsize
                for x = 1:H4mxStride:H4Outsize
                    kern1 = tmp21(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter response
                    kern2 = tmp22(y:y+MXpoolDim-1,x:x+MXpoolDim-1);%%%filter residual
                    
                    [Val,IDX] = nanmin(abs(kern2(:)));%%%get abs min error
                    [Y X] = ind2sub(size(kern2),IDX);
                    H4pool{node}(y,x) = kern1(Y,X);
                    H4RESpool{node}(y,x) = kern2(Y,X);
                end
            end
            H4Out(:,:,node) = max(H4pool{node}.*LreLu,H4pool{node});
            [RowColOut4{node,1}(:,1) RowColOut4{node,1}(:,2)] = find(H4Out(:,:,node)~=0);
            
            count = 0;
            for y = 1:H4RespDim
                for x = 1:H4RespDim
                    count = count+1;
                    tmp = H4pool{node}(y,x);
                    [RowColIm4{node,1}(count,1) RowColIm4{node,1}(count,2)] = find(H4resp == tmp,1);
                end
            end

        end
        
        %%%%Fully Connected Layer 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        [h w d] = size(H1Out);
        H1OutV = reshape(H1Out,h.*w.*d,1);
        
        [h w d] = size(H2Out);
        H2OutV = reshape(H2Out,h.*w.*d,1);
        
        [h w d] = size(H3Out);
        H3OutV = reshape(H3Out,h.*w.*d,1);
        
        [h w d] = size(H4Out);
        H4OutV = reshape(H4Out,h.*w.*d,1);
        
        FC1dot = FC1w*H4OutV;
        FC1_Eval = 1./(1+exp(-Sig_b.*(FC1dot-Sig_Hmax)));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
        %%%%Output Layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        Wout = sum(FC1_Eval.*OutW);
        Netw_Eval = 1./(1+exp(-Sig_b.*(Wout-Sig_Hmax)));
        
        TrackPatternClassification(IterLoop,1) = Netw_Eval;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        ERR = (TargetResp-Netw_Eval);
        ESUM = ESUM+abs(ERR);
        
        %%%Binary cross entropy/log loss
        %ERR = -((TargetResp.*log(Netw_Eval))+(1-TargetResp).*log(1-Netw_Eval));
        
        
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%Backpropagation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        OutSlope = Netw_Eval.*(1-Netw_Eval);
        FC1Slope = FC1_Eval.*(1-FC1_Eval);

        
        %%%H4 Layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        LossGradW = ((FC1w'*FC1Slope.*OutSlope).*ERR);
        LossGradW = reshape(LossGradW,H4poolsize,H4poolsize,H4_Nodes);

        OrigH4W = H4W;
        
        H3OutPad = padarray(H3Out,[H4dimhlf H4dimhlf],0,'both');
        
        dW = zeros(H4dim,H4dim,H3_Nodes,H4_Nodes);
        iterCount = 0;
        ResTally4 = 0;

        parfor node = 1:H4_Nodes
            
            Coords = RowColIm4{node};
            gCoords = RowColOut4{node};
            numCoord = size(Coords,1);
            
            for Cloop = 1:numCoord
                iterCount = iterCount+1;
                Y = Coords(Cloop,1)+H4dimhlf;%%%shift to account for padding
                X = Coords(Cloop,2)+H4dimhlf;
                gY = gCoords(Cloop,1);
                gX = gCoords(Cloop,2);
                
                GradW = LossGradW(gY,gX,node);
                ResW = H4RESpool{node}(gY,gX);
                ResTally4 = ResTally4+abs(ResW);
                Patch = H3OutPad(Y-H4dimhlf:Y+H4dimhlf,X-H4dimhlf:X+H4dimhlf,:);
                Patch = Patch.*(max(sign(Patch).*LreLu,sign(Patch)).*(sign(Patch))); %%%%Calculate and apply leaky reLu here
                
                %%%Regular model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                dW(:,:,:,node) = dW(:,:,:,node)+((H4w(:,:,:,node).*Patch.*GradW.*ResW).*mean2(abs(H4RESpool{node})));   %%%%adjusting by residual mean so that filters that were doing well aren't overly corrected by network error
                %dW(:,:,:,node) = dW(:,:,:,node)+((H4w(:,:,:,node).*Patch.*GradW));
                
            end
            
            H4w(:,:,:,node) = H4w(:,:,:,node)+((Lrate_H4.*dW(:,:,:,node))); 
            
            %%%%%%L2 norm
            tmp = reshape(H4w(:,:,:,node),H4dim*H4dim*H3_Nodes,1);
            tmp = tmp./norm(tmp);
            H4w(:,:,:,node) =reshape(tmp,H4dim,H4dim,H3_Nodes);
            
        end
        
        H4lossGrad = LossGradW;
        
        RESSUM4 = RESSUM4+(ResTally4./iterCount);
        
        
        %%%H3 Layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        rH4w = rot90(rot90(OrigH4W));
        LossGradW = convn(H3Out,rH4w,'same');

        OrigH3W = H3w;
        
        %%norm to control smashing H3 weights
        tmp = reshape(LossGradW,H3poolsize*H3poolsize*H3_Nodes,1);
        tmp = tmp./norm(tmp);
        LossGradW = reshape(tmp,H3poolsize,H3poolsize,H3_Nodes).*OutSlope.*ERR; %%%%<----error and outslope applied here because of normalization
        
        H2OutPad = padarray(H2Out,[H3dimhlf H3dimhlf],0,'both');
        
        dW = zeros(H3dim,H3dim,H2_Nodes,H3_Nodes);
        iterCount = 0;
        ResTally3 = 0;

        parfor node = 1:H3_Nodes
            
            Coords = RowColIm3{node};
            gCoords = RowColOut3{node};
            numCoord = size(Coords,1);
            
            for Cloop = 1:numCoord
                iterCount = iterCount+1;
                Y = Coords(Cloop,1)+H3dimhlf;%%%shift to account for padding
                X = Coords(Cloop,2)+H3dimhlf;
                gY = gCoords(Cloop,1);
                gX = gCoords(Cloop,2);
                
                GradW = LossGradW(gY,gX,node);
                ResW = H3RESpool{node}(gY,gX);
                ResTally3 = ResTally3+abs(ResW);
                Patch = H2OutPad(Y-H3dimhlf:Y+H3dimhlf,X-H3dimhlf:X+H3dimhlf,:);
                Patch = Patch.*(max(sign(Patch).*LreLu,sign(Patch)).*(sign(Patch))); %%%%Calculate and apply leaky reLu here
                
                %%%Regular model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                dW(:,:,:,node) = dW(:,:,:,node)+((H3w(:,:,:,node).*Patch.*GradW.*ResW).*mean2(abs(H3RESpool{node})));   %%%%adjusting by residual mean so that filters that were doing well aren't overly corrected by network error
                %dW(:,:,:,node) = dW(:,:,:,node)+((H3w(:,:,:,node).*Patch.*GradW));
                
            end
            
            H3w(:,:,:,node) = H3w(:,:,:,node)+((Lrate_H3.*dW(:,:,:,node))); 
            
            %%%%%%L2 norm
            tmp = reshape(H3w(:,:,:,node),H3dim*H3dim*H2_Nodes,1);
            tmp = tmp./norm(tmp);
            H3w(:,:,:,node) =reshape(tmp,H3dim,H3dim,H2_Nodes);
            
        end
        
        H3lossGrad = LossGradW;
        
        RESSUM3 = RESSUM3+(ResTally3./iterCount);
        
        
        %%%H2 Layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        rH3w = rot90(rot90(OrigH3W));
        LossGradW = convn(H2Out,rH3w,'same');

        OrigH2W = H2W;
        
        %%norm to control smashing H2 weights
        tmp = reshape(LossGradW,H2poolsize*H2poolsize*H2_Nodes,1);
        tmp = tmp./norm(tmp);
        LossGradW = reshape(tmp,H2poolsize,H2poolsize,H2_Nodes).*OutSlope.*ERR; %%%%<----error and outslope applied here because of normalization
        
        H1OutPad = padarray(H1Out,[H2dimhlf H2dimhlf],0,'both');
        
        dW = zeros(H2dim,H2dim,H1_Nodes,H2_Nodes);
        iterCount = 0;
        ResTally2 = 0;

        parfor node = 1:H2_Nodes
            
            Coords = RowColIm2{node};
            gCoords = RowColOut2{node};
            numCoord = size(Coords,1);
            
            for Cloop = 1:numCoord
                iterCount = iterCount+1;
                Y = Coords(Cloop,1)+H2dimhlf;%%%shift to account for padding
                X = Coords(Cloop,2)+H2dimhlf;
                gY = gCoords(Cloop,1);
                gX = gCoords(Cloop,2);
                
                GradW = LossGradW(gY,gX,node);
                ResW = H2RESpool{node}(gY,gX);
                ResTally2 = ResTally2+abs(ResW);
                Patch = H1OutPad(Y-H2dimhlf:Y+H2dimhlf,X-H2dimhlf:X+H2dimhlf,:);
                Patch = Patch.*(max(sign(Patch).*LreLu,sign(Patch)).*(sign(Patch))); %%%%Calculate and apply leaky reLu here
                
                %%%Regular model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                dW(:,:,:,node) = dW(:,:,:,node)+((H2w(:,:,:,node).*Patch.*GradW.*ResW).*mean2(abs(H2RESpool{node})));   %%%%adjusting by residual mean so that filters that were doing well aren't overly corrected by network error
                %dW(:,:,:,node) = dW(:,:,:,node)+((H2w(:,:,:,node).*Patch.*GradW));
                
            end
            
            H2w(:,:,:,node) = H2w(:,:,:,node)+((Lrate_H2.*dW(:,:,:,node))); 
            
            %%%%%%L2 norm
            tmp = reshape(H2w(:,:,:,node),H2dim*H2dim*H1_Nodes,1);
            tmp = tmp./norm(tmp);
            H2w(:,:,:,node) =reshape(tmp,H2dim,H2dim,H1_Nodes);
            
        end
        
        H2lossGrad = LossGradW;
        
        RESSUM2 = RESSUM2+(ResTally2./iterCount);
        
        %%%H1 Layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        rH2w = rot90(rot90(OrigH2W));
        LossGradW = convn(H1Out,rH2w,'same');
        
        %%norm to control smashing H1 weights
        tmp = reshape(LossGradW,H1poolsize*H1poolsize*H1_Nodes,1);
        tmp = tmp./norm(tmp);
        LossGradW = reshape(tmp,H1poolsize,H1poolsize,H1_Nodes).*OutSlope.*ERR; %%%%<----error and outslope applied here because of normalization
        

        L1pad = padarray(L1,[H1dimhlf H1dimhlf],0,'both');
        dW = zeros(H1dim,H1dim,H1_Nodes);
        iterCount = 0;
        ResTally1 = 0;


        parfor node = 1:H1_Nodes
            
            Coords = RowColIm{node};
            gCoords = RowColOut{node};
            numCoord = size(Coords,1);
            
            for Cloop = 1:numCoord
                iterCount = iterCount+1;
                Y = Coords(Cloop,1)+H1dimhlf;%%%shift to account for padding
                X = Coords(Cloop,2)+H1dimhlf;
                gY = gCoords(Cloop,1);
                gX = gCoords(Cloop,2);
                
                GradW = LossGradW(gY,gX,node);
                ResW = H1RESpool{node}(gY,gX);
                ResTally1 = ResTally1+abs(ResW);
                Patch = L1pad(Y-H1dimhlf:Y+H1dimhlf,X-H1dimhlf:X+H1dimhlf);
                
                %%%Regular model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                a = rot90(rot90(H1w(:,:,node))).*Patch; %%%account for rotation for conv2 rotation in layer 1
                HW1slope = a.*(max(sign(a).*LreLu,sign(a)).*(sign(a))); %%%%Get leaky reLu slope here
                dW(:,:,node) = dW(:,:,node)+((HW1slope.*rot90(rot90(H1w(:,:,node))).*Patch.*GradW.*ResW).*mean2(abs(H1RESpool{node}))); %%%adjusting by residual mean so that filters that were doing well aren't overly corrected by network error
                %dW(:,:,node) = dW(:,:,node)+((HW1slope.*H1w(:,:,node).*Patch.*GradW));
            end
            
            rot90(rot90(H1w(:,:,node))) = rot90(rot90(H1w(:,:,node)))+(Lrate_H1.*dW(:,:,node));
            
            %%%%%%L2 norm
            tmp = reshape(H1w(:,:,node),H1dim*H1dim,1);
            tmp = tmp./norm(tmp);
            H1w(:,:,node) =reshape(tmp,H1dim,H1dim);
            
        end
        
        H1lossGrad = LossGradW;
        
        RESSUM1 = RESSUM1+(ResTally1./iterCount);

        
        %%%%%FC layer 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        dOutdW = ((FC1w.*FC1Slope).*H4OutV').*OutSlope;
        FC1w = FC1w+(Lrate.*(dOutdW.*ERR));

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%%%%Output layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        dOutdW = FC1_Eval.*OutSlope;
        OutW = OutW+(Lrate.*(dOutdW.*ERR));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
    end
    
    TRACK_ERROR(SimNum,1) = ESUM;
    TRACK_resERROR1(SimNum,1) = RESSUM1;
    TRACK_resERROR2(SimNum,1) = RESSUM2;
    TRACK_resERROR3(SimNum,1) = RESSUM3;
    TRACK_resERROR4(SimNum,1) = RESSUM4;
    
    tmpPL = [];
    tmpPL(:,1) = TrackPatternClassification;
    tmpPL(:,2) = IterInfo(:,8);
    tmpPL = sortrows(tmpPL,2);
    tmpPL = tmpPL(:,1);
    
    ACC = (sum(tmpPL(1:IterHlf)<=.40)+sum(tmpPL(IterHlf+1:IterNum)>=.60))./IterNum;
    
    figure(f)
    subplot(3,5,1)
    plot(TRACK_ERROR,'-o');
    title('Network Error');
    subplot(3,5,2)
    bar(tmpPL);
    title(strcat('Accuracy:',{' '},num2str(ACC)));
    xline(IterHlf+1,'-','linewidth',3);
    yline(.4,'--','linewidth',2);
    yline(.5,'-','linewidth',2);
    yline(.6,'--','linewidth',2);
    subplot(3,5,3)
    plot(TRACK_resERROR1,'-o');
    title('Residual Sum hL1');
    subplot(3,5,4)
    plot(TRACK_resERROR2,'-o');
    title('Residual Sum hL2');
    subplot(3,5,5)
    plot(TRACK_resERROR3,'-o');
    title('Residual Sum hL3');
    subplot(3,5,6)
    plot(TRACK_resERROR4,'-o');
    title('Residual Sum hL4');    
    subplot(3,5,7)
    histogram(H1w);
    title('HL1 Weights');
    subplot(3,5,8)
    histogram(H1lossGrad);
    title('HL1 Loss Grad');
    subplot(3,5,9)
    histogram(H2w);
    title('HL2 Weights');
    subplot(3,5,10)
    histogram(H2lossGrad);
    title('HL2 Loss Grad'); 
    subplot(3,5,11)
    histogram(H3w);
    title('HL3 Weights');
    subplot(3,5,12)
    histogram(H3lossGrad);
    title('HL3 Loss Grad');
    subplot(3,5,13)
    histogram(H4w);
    title('HL4 Weights');
    subplot(3,5,14)
    histogram(H4lossGrad);
    title('HL4 Loss Grad');        
    
    pause(.01)
    
   
    
    save(strcat('Path to where you want to save Conv 1 weights /BrainGuidedCNN_Conv1.mat'),'H1w');
    save(strcat('Path to where you want to save Conv 2 weights /BrainGuidedCNN_Conv2.mat'),'H2w');
    save(strcat('Path to where you want to save Conv 3 weights /BrainGuidedCNN_Conv3.mat'),'H3w');
    save(strcat('Path to where you want to save Conv 4 weights /BrainGuidedCNN_Conv4.mat'),'H4w');
    save(strcat('Path to where you want to save fully connected layer weights weights /BrainGuidedCNN_FC1.mat'),'FC1w');
    save(strcat('Path to where you want to save decision node weights /BrainGuidedCNN_OutLayer.mat'),'OutW');
    
    INFO.Conv1Nodes = H1_Nodes;
    INFO.Conv2Nodes = H2_Nodes;
    INFO.Conv3Nodes = H3_Nodes;
    INFO.Conv4Nodes = H4_Nodes;
    INFO.FC1Nodes = FC1_Nodes;
    INFO.SIMNUM = SimNum;
    INFO.Lrate_FC = Lrate;
    INFO.Lrate_H1 = Lrate_H1;
    INFO.Lrate_H2 = Lrate_H2;
    INFO.Lrate_H3 = Lrate_H3;
    INFO.Lrate_H4 = Lrate_H4;
    INFO.TRACK_ERROR = TRACK_ERROR;
    INFO.TRACK_resERROR1 = TRACK_resERROR1;
    INFO.TRACK_resERROR2 = TRACK_resERROR2;
    INFO.TRACK_resERROR3 = TRACK_resERROR3;
    INFO.TRACK_resERROR4 = TRACK_resERROR4;
    INFO.Performance = tmpPL;
    save(strcat('Pat to where you want to save parameters /BrainGuidedCNN_Info.mat'),'INFO');
    
    
end








