clear all;close all;

realvals=[0.1174;0.1467;0.2515;0.3427];
shuffvals1=[];shuffvals2=[];shuffvals3=[];shuffvals4=[];


% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\635\day9'; 



% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];


spike_data=data;
for num=1:1000;
% Get dimensions of the matrix
[numRows, numCols] = size(spike_data);

% Generate a vector of random shift amounts for each column
shiftAmounts = randi([0, numRows - 1], 1, numCols);

B_cell = arrayfun(@(colIdx) circshift(spike_data(:, colIdx), shiftAmounts(colIdx)), 1:numCols, 'UniformOutput', false);

% Convert the cell array back to a matrix
B_cell = cell2mat(B_cell);

% Assumed original matrix: `spike_data`
num_rows = size(B_cell, 1);
num_neurons = size(B_cell, 2);

% Calculate the number of rows to pad
rows_to_pad = mod(-num_rows, 3);

% Pad the matrix with zeros at the end
padded_data = [B_cell; zeros(rows_to_pad, num_neurons)];


% Reshape the padded data into a 3D array
num_new_rows = size(padded_data, 1) / 3;
binned_3d = reshape(padded_data, 3, num_new_rows, num_neurons);

% Sum the spikes within each bin
summed_binned = squeeze(sum(binned_3d, 1));

% Convert the summed data back to binary (0 or 1)
final_binned_data = double(summed_binned > 0);

sync_val=(sum(final_binned_data,2));
sync_val=sync_val/size(data,2);
shuffvals1=[shuffvals1;mean(nonzeros(sync_val))];
end;




% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\681\day9'; 



% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];


spike_data=data;
for num=1:1000;
% Get dimensions of the matrix
[numRows, numCols] = size(spike_data);

% Generate a vector of random shift amounts for each column
shiftAmounts = randi([0, numRows - 1], 1, numCols);

B_cell = arrayfun(@(colIdx) circshift(spike_data(:, colIdx), shiftAmounts(colIdx)), 1:numCols, 'UniformOutput', false);

% Convert the cell array back to a matrix
B_cell = cell2mat(B_cell);

% Assumed original matrix: `spike_data`
num_rows = size(B_cell, 1);
num_neurons = size(B_cell, 2);

% Calculate the number of rows to pad
rows_to_pad = mod(-num_rows, 3);

% Pad the matrix with zeros at the end
padded_data = [B_cell; zeros(rows_to_pad, num_neurons)];


% Reshape the padded data into a 3D array
num_new_rows = size(padded_data, 1) / 3;
binned_3d = reshape(padded_data, 3, num_new_rows, num_neurons);

% Sum the spikes within each bin
summed_binned = squeeze(sum(binned_3d, 1));

% Convert the summed data back to binary (0 or 1)
final_binned_data = double(summed_binned > 0);

sync_val=(sum(final_binned_data,2));
sync_val=sync_val/size(data,2);
shuffvals2=[shuffvals2;mean(nonzeros(sync_val))];
end;



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\691\day9'; 



% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];


spike_data=data;
for num=1:1000;
% Get dimensions of the matrix
[numRows, numCols] = size(spike_data);

% Generate a vector of random shift amounts for each column
shiftAmounts = randi([0, numRows - 1], 1, numCols);

B_cell = arrayfun(@(colIdx) circshift(spike_data(:, colIdx), shiftAmounts(colIdx)), 1:numCols, 'UniformOutput', false);

% Convert the cell array back to a matrix
B_cell = cell2mat(B_cell);

% Assumed original matrix: `spike_data`
num_rows = size(B_cell, 1);
num_neurons = size(B_cell, 2);

% Calculate the number of rows to pad
rows_to_pad = mod(-num_rows, 3);

% Pad the matrix with zeros at the end
padded_data = [B_cell; zeros(rows_to_pad, num_neurons)];


% Reshape the padded data into a 3D array
num_new_rows = size(padded_data, 1) / 3;
binned_3d = reshape(padded_data, 3, num_new_rows, num_neurons);

% Sum the spikes within each bin
summed_binned = squeeze(sum(binned_3d, 1));

% Convert the summed data back to binary (0 or 1)
final_binned_data = double(summed_binned > 0);

sync_val=(sum(final_binned_data,2));
sync_val=sync_val/size(data,2);
shuffvals3=[shuffvals3;mean(nonzeros(sync_val))];
end;



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\747\day9'; 



% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];


spike_data=data;
for num=1:1000;
% Get dimensions of the matrix
[numRows, numCols] = size(spike_data);

% Generate a vector of random shift amounts for each column
shiftAmounts = randi([0, numRows - 1], 1, numCols);

B_cell = arrayfun(@(colIdx) circshift(spike_data(:, colIdx), shiftAmounts(colIdx)), 1:numCols, 'UniformOutput', false);

% Convert the cell array back to a matrix
B_cell = cell2mat(B_cell);

% Assumed original matrix: `spike_data`
num_rows = size(B_cell, 1);
num_neurons = size(B_cell, 2);

% Calculate the number of rows to pad
rows_to_pad = mod(-num_rows, 3);

% Pad the matrix with zeros at the end
padded_data = [B_cell; zeros(rows_to_pad, num_neurons)];


% Reshape the padded data into a 3D array
num_new_rows = size(padded_data, 1) / 3;
binned_3d = reshape(padded_data, 3, num_new_rows, num_neurons);

% Sum the spikes within each bin
summed_binned = squeeze(sum(binned_3d, 1));

% Convert the summed data back to binary (0 or 1)
final_binned_data = double(summed_binned > 0);

sync_val=(sum(final_binned_data,2));
sync_val=sync_val/size(data,2);
shuffvals4=[shuffvals4;mean(nonzeros(sync_val))];
end;


count = sum(shuffvals1 >= realvals(1));
p_one_sided = (count + 1) / (length(shuffvals1) + 1)


count = sum(shuffvals2 >= realvals(2));
p_one_sided = (count + 1) / (length(shuffvals2) + 1)


count = sum(shuffvals3 >= realvals(3));
p_one_sided = (count + 1) / (length(shuffvals3) + 1)


count = sum(shuffvals4 >= realvals(4));
p_one_sided = (count + 1) / (length(shuffvals4) + 1)

figure;histogram(shuffvals1);hold on; vline(realvals(1),'r')
figure;histogram(shuffvals2);hold on; vline(realvals(2),'r')
figure;histogram(shuffvals3);hold on; vline(realvals(3),'r')
figure;histogram(shuffvals4);hold on; vline(realvals(4),'r')
