%% Obtain Spine density and dynamics values
% For in vivo spine data analyzed using MapManager V_20171028 
% Author: Richard Roth (rhroth@stanford.edu)
% Date: Sept 2023

%     Requires the following functions:
%     getmyMapListValuesCond.m
%     segmentCall.m
%     removeNanRows.m


%% Get Values
clearvars mySpinePool 

Condition = 'Gr1';
% Condition = 'Gr2';
SessionConditions = {'day-3', 'day-2', 'day-1', 'day0', 'day1', 'day2', 'day3', 'day4'};
StatName = 'pDist';

channel = 1;
PersCond = 1; %Minimum number of sessions that spine needs to be present in
MinSpines = 1; %Minimum number of spines per dendrite for dendrite to be included
myPool = getmyMapListValuesCond(myMapList, mySegmentTable, Condition, SessionConditions, StatName, channel, PersCond, MinSpines);

StatName = 'nntbutsmdtdmch1'; %Spine dsRed / Dendrite dsRed mean normalized

channel = 1;
PersCond = 1; %Minimum number of sessions that spine needs to be present in
MinSpines = 1; %Minimum number of spines per dendrite for dendrite to be included
myPoolInt = getmyMapListValuesCond(myMapList, mySegmentTable, Condition, SessionConditions, StatName, channel, PersCond, MinSpines);


%%
fieldNames = fieldnames(myPool.values.SpinePerSegmentNanRemoved);
for i = 1:length(fieldNames)
currentfieldName = fieldNames{i};
currentmatrix = myPool.values.SpinePerSegmentNanRemoved.(currentfieldName);

currentspinepDist = mean(currentmatrix, 2, 'omitnan'); %calculate average spine position

currentdendritelength = max(currentspinepDist)-min(currentspinepDist);

spinenum = sum(~isnan(currentmatrix)); %count number of spines
spinenum(spinenum == 0) = NaN;

[numRows, numCols] = size(currentmatrix);
Idx_add = NaN(numRows, 1);
Idx_sub = NaN(numRows, 1);


for row = 1:numRows
    %find the first column that is non NaN (session when spine is added)
    nonNaNCol_first = find(~isnan(currentmatrix(row, :)), 1, 'first');  
    Idx_add(row) = nonNaNCol_first;
    
    %find the last column that is non NaN (session when spine will be eliminated)
    nonNaNCol_last = find(~isnan(currentmatrix(row, :)));
    Idx_sub(row) = nonNaNCol_last(end);
end


if nnz(~isnan(spinenum)) == 4
spinenum_for_add = [NaN, NaN, NaN, spinenum(2), NaN, spinenum(4), NaN, spinenum(6)];
spinenum_for_sub = spinenum;
uniqueValues_first = [2,4,6,8];
uniqueValues_last = [2,4,6,8];

% Count occurrences using histcounts
occurrences_add = histcounts(Idx_add, [uniqueValues_first, max(uniqueValues_first) + 1]);
occurrences_add(1)=NaN;
occurrences_add = [NaN occurrences_add(1), NaN, occurrences_add(2), NaN, occurrences_add(3), NaN, occurrences_add(4)];

occurrences_sub = histcounts(Idx_sub, [uniqueValues_last, max(uniqueValues_last) + 1]);  %this is going by session in which spine is last seen (spine will be eliminated after this session)
occurrences_sub(end)=NaN;
occurrences_sub = [NaN occurrences_sub(1), NaN, occurrences_sub(2), NaN, occurrences_sub(3), NaN, occurrences_sub(4)];
occurrences_sub_2 = [NaN, NaN, occurrences_sub(1:end-2)]; %this is going by session after session in which spine is last seen (spine has been eliminated before this session)
else
    
spinenum_for_add = [NaN, spinenum([1,2,3,4,5,6,7])];
spinenum_for_sub = spinenum;
uniqueValues_first =[1,2,3,4,5,6,7,8];
uniqueValues_last =[1,2,3,4,5,6,7,8];

% Count occurrences using histcounts
occurrences_add = histcounts(Idx_add, [uniqueValues_first, max(uniqueValues_first) + 1]);
occurrences_add(1)=NaN;

occurrences_sub = histcounts(Idx_sub, [uniqueValues_last, max(uniqueValues_last) + 1]); %this is going by session in which spine is last seen (spine will be eliminated after this session)
occurrences_sub(end)=NaN;
occurrences_sub_2 = [NaN, occurrences_sub(1:end-1)]; %this is going by session after session in which spine is last seen (spine has been eliminated before this session)
end

percentages_add = occurrences_add./spinenum_for_add;
percentages_sub_2 = occurrences_sub_2./spinenum_for_add;

mySpinePool.spinenum.(currentfieldName) = spinenum;
mySpinePool.spinedensity.(currentfieldName) = spinenum./currentdendritelength;
mySpinePool.pAdd.(currentfieldName) = percentages_add;
mySpinePool.pSub_2.(currentfieldName) = percentages_sub_2;
end

spineNum_out = cell2mat(struct2cell(mySpinePool.spinenum));
spineDensity_out = cell2mat(struct2cell(mySpinePool.spinedensity));
spineAdd_out = cell2mat(struct2cell(mySpinePool.pAdd));
spineSub2_out = cell2mat(struct2cell(mySpinePool.pSub_2));

spineTOR_out = spineAdd_out + spineSub2_out;

%% make averages of neighboring 2 days

output_spineDensity_2davg = squeeze(mean(reshape(spineDensity_out(:,1:end-mod(size(spineDensity_out,2),2)),size(spineDensity_out,1),2,[]),2,'omitnan')); % SPINE DENSITY VALUES
output_spineAdd_2davg = squeeze(mean(reshape(spineAdd_out(:,1:end-mod(size(spineAdd_out,2),2)),size(spineAdd_out,1),2,[]),2,'omitnan')); % SPINE ADDITION VALUES
output_spineSub2_2davg = squeeze(mean(reshape(spineSub2_out(:,1:end-mod(size(spineSub2_out,2),2)),size(spineSub2_out,1),2,[]),2,'omitnan')); % SPINE SUBTRACTION VALUES
output_spineTOR_2davg = squeeze(mean(reshape(spineTOR_out(:,1:end-mod(size(spineTOR_out,2),2)),size(spineTOR_out,1),2,[]),2,'omitnan')); % SPINE TURNOVER RATIO

 