% Script 5
% goals of this script:
% - compare PCA results across subjects
% if you have questions, feel free to contact Briana R. McRae: brm2144@columbia.edu

% NOTE: this script also requires specialization to your experiment,
% specifically carrying out the group comparisons of interest

%%  Section 1: read in PCA results

clc; clear all; close all;

% inputs: edit!
outputFolder = '/path/to/supplementalMaterials/exampleData/outputs/';

% analysis options: edit!
% NOTE: for example data, group 1 is mothers, group 2 is virgins, 
% stimulus 1 is pup calls, and stimulus 2 is pure tones
stimnum = 1; % which stimulus event type to analyze
dsfactor = 1; % spatial downsample factor (default = 1, no downsampling)
fnfactor = 1; % frame number downsample factor (default = 1, no downsampling)
animals = {'Br011', 'Br014'};
group1idx = [2];
group2idx = [1];
PCn = 1; % which PC to look at

% set up folders
keyInfo = ['stimnum', num2str(stimnum), 'dsfactor', num2str(dsfactor), 'fnfactor', num2str(fnfactor)];
pcaOutputPath = [outputFolder, 'pcaOutputs_', keyInfo, '/']; % makes a new folder for each combination of stimnum, dsfactor, and fnfactor
cd(outputFolder)

% ROI input
inputCoords = load('ROIcoords.mat'); % edit!
coords = inputCoords.coords;
ROIwidth = round(coords(1)):round(coords(1)+coords(3));
xd = ROIwidth(end)-ROIwidth(1)+1;
ROIheight = round(coords(2)):round(coords(2)+coords(4));
yd = ROIheight(end)-ROIheight(1)+1;

pcaStruct = struct;

% edit the code block below!
for anm=1:length(animals)
    pcaStruct.(animals{anm}) = [];
    runs = dir([pcaOutputPath '*' animals{anm} '*']);
    for rr=1:length(runs)
        load([runs(rr).folder '/' runs(rr).name])
        pcaStruct.(animals{anm}) = [pcaStruct.(animals{anm}) pcs(:,PCn)];
    end
end

disp('done with script 5, section 1!')

%% Section 2: define groups and average PCA results per group

% all animals
all_mat = [];
for anm=1:length(animals)
    all_mat = [all_mat; mean(pcaStruct.(animals{anm}),2)'];
end

% group1
group1_mat = [];
for anm=group1idx
    group1_mat = [group1_mat; mean(pcaStruct.(animals{anm}),2)'];
end

% group2
group2_mat = [];
for anm=group2idx
    group2_mat = [group2_mat; mean(pcaStruct.(animals{anm}),2)'];
end

disp('done with script 5, section 2!')

%% Section 3: visualize PC loadings per group

% inputs: edit!
group = (all_mat); % options: all_mat, group1_mat, or group2_mat
myTitle = ['All animals: stimulus ', int2str(stimnum)];
threshT = 0; % edit! see publicly-available T-statistic to p-value tables

% calculate magnitude of PC loadings
mag = mean(group);

% calculate T-statistics
Tscore = sqrt(size(group,1))*mean(group)./std(group);

% visualize magnitude map thresholded at certain T (determined by p-value of choice and df)
figure()
exhibit = mag.*(Tscore>threshT);
exhibit(exhibit==0) = NaN;
toshow = reshape(exhibit, [yd xd]);
img = imrotate(toshow,-90); % edit! rotate according to orientation of videos
imalpha = ones(size(img));
imalpha(isnan(img)) = 0;
imagesc(img,'AlphaData',imalpha)
colormap hot
axis equal
xlim([1,yd])
ylim([1,xd])
caxis([0,.01]) % edit!
title(myTitle)
colorbar

disp('done with script 5, section 3!')

%% Section 4: group comparison (group1 > group2)

% inputs: edit based on analysis needs!
mag = (mean(group1_mat,1)-mean(group2_mat,1));
myTitle = ['Group 1 > group 2: stimulus ', int2str(stimnum)]; % edit according to group comparison chosen
threshT = 0; % edit! see publicly-available T-statistic to p-value tables

% calculate T-statistics
Tscore = (mean(group1_mat,1)-mean(group2_mat,1))./sqrt( ((std(group1_mat).^2)/(size(group1_mat,1))) + ((std(group2_mat).^2)/(size(group2_mat,1)))); % also here?!!?!??!

% visualize magnitude map thresholded at certain T (determined by p-value of choice and df)
figure()
exhibit = mag.*(Tscore>threshT);
exhibit(exhibit==0) = NaN;
toshow = reshape(exhibit, [yd xd]);
img = imrotate(toshow,-90); % edit! rotate according to orientation of videos
imalpha = ones(size(img));
imalpha(isnan(img)) = 0;
imagesc(img,'AlphaData',imalpha)
colormap hot
caxis([0,.01]) % edit!
axis equal
xlim([1,yd])
ylim([1,xd])
title(myTitle)
colorbar

disp('done with script 5, section 4!')