clc 
clear all
close all
%% Choose which output to plot
helix_2d = 0;
helix_down_layer1 = 0;
helix_up_layer1 = 0;
% pegboard_2d_layer1 = 0;
pegboard_3d = 1;

%% load the files
% if helix_2d == 1
%     out = readtable('encoded_helix_t20p3_b7_2d_train.csv');
%     traj_data = readtable('traj_wo_20k.csv');
% end
% if helix_down_layer1 == 1
%     out = readtable('encoded_helix_t20p3_b7_3d_test_down.csv');
%     traj_data = readtable('Trajectory_fivecoil_down.csv');
% end
% if helix_up_layer1 == 1
%     out = readtable('encoded_helix_t20p3_b7_3d_test_up.csv');
%     traj_data = readtable('Trajectory_fivecoil_up.csv');
% end
% if pegboard_2d_layer1 == 1
%     out = readtable('encoded_Pegboard_2d_l1.csv');
%     traj_data = readtable('traj_peg2d_50k.csv');
% end
if pegboard_3d == 1
    out = readtable('encoded_peg_t20p3_b0.4pi_3d_testA.csv');
    traj_data = readtable('Trajectory_interpolated_pegboard_new.csv');
end

out = table2array(out);
traj_data = traj_data(:,1:3);
% if pegboard_2d_layer1
%     traj_data = traj_data(:,1:2);
% end
% if pegboard_3d_layer1
%     traj_data = traj_data(:,2:3);
% end
pos = table2array(traj_data);
%% plotting output
if helix_2d || helix_down_layer1 || helix_up_layer1 
    grid_cells = [6, 10, 13, 14, 15, 18, 23, 26, 34, 37, 38, 40];
    place_cells = [ 2, 4, 7, 8, 9, 12, 16, 17, 19, 20, 21, 22, 24, 25, 27, 29, 31, 35, 36, 39, 42, 43, 44, 46, 47, 48, 49, 50];
 else
%     grid_cells = [ 4  8 13 24 28 30 34 35 37 42 48];
%     place_cells = [ 1 2 3 5 9 10 11 12 14 15 16 19 21 22 26 27 29 33 36 39 40 43 44 46 47 48 49 50];
end   
for n = place_cells
    n
    ot = out(:,n);
    ot_mean = mean(ot);
    ot_std = std(ot);
    thresh = ot_mean + 1.5*ot_std;
    %ot = out(n,:);  
    %thresh=max(ot)*lim;  
    firr=find((ot)>thresh); 
    firposgrid=pos(firr,:); %Firing positions on the trajectory
    if helix_down_layer1 || helix_up_layer1 || helix_2d
        res = 35;
        x_mesh = -1:1/res:1;
        y_mesh = -1:1/res:1; 
    end
    if pegboard_2d_layer1 || pegboard_3d_layer1
        res = 4;
        x_mesh = 0:1/res:10;
        y_mesh = 0:1/res:10;
    end
    [fx,fy] = meshgrid(x_mesh,y_mesh);

    firingmap = zeros(length(fx));
    gridpoint = [reshape(fx,prod(size(fx)),1) reshape(fy,prod(size(fy)),1)];
    roundinggridpoint = round(gridpoint);
    firposround = round(firposgrid); 
    firingvalue = abs(ot(firr));
    for ii = 1:length(firposgrid(:,1))
        [~,q1]=min(abs(firposgrid(ii,1)-fx(1,:)));
        [~,q2]=min(abs(firposgrid(ii,2)-fx(1,:)));
        firingmap(q1,q2) = firingvalue(ii);
    end
    firingmap = firingmap/max(max(firingmap));
    gaussian = fspecial('gaussian',[10 10],3); % Design the gaussian filter for smoothing the map
    spikes_smooth=conv2(gaussian,firingmap);
    spikes_smooth = imrotate(spikes_smooth/max(max(spikes_smooth)),90);
    spikes_smooth = spikes_smooth(end-length(y_mesh)-5:end,:);

clf;
figure(1);
subplot(2,1,1);plot(pos(:,1),pos(:,2),'b');
axis equal;
title(['Neuron ',num2str(n)], 'Fontsize', 20);
hold on;
plot(firposgrid(:,1),firposgrid(:,2),'.r', 'markersize', 10);axis off
subplot(2,1,2);imagesc(spikes_smooth); axis off
colormap(jet)
axis equal;
savefig(['place_2D', num2str(n)]);
saveas(gcf, ['place_2D', num2str(n), '.png'])
% saveas(gcf,sprintf('place_helixD_%d.png',n))
% saveas(gcf,sprintf('place_helixD_%d.fig',n))
pause;
close(gcf);
end
% %% firing rate maps 
% if pegboard_2d_layer1
% hex = 0;
% squ = 0;
% for n = 1:50
%     ot = out(:,n); 
%     %ot = out(:,n);
%     ot_mean = mean(ot);
%     ot_std = std(ot);
%     thresh = ot_mean + 1.5*ot_std;
%     %ot = out(n,:);  
%     %thresh=max(ot)*lim;  
%     firr=find((ot)>thresh); 
%     %thresh=max(ot)*.80; 
%     %firr=find((ot)>thresh);
%     firposgrid=pos(firr,:); %Firing positions on the trajectory
% 
%     res = 10; %resolution of the rate map
%     fx = meshgrid(0:1/res:10);
%     fy = meshgrid(0:1/res:10);
%     firingmap = zeros(length(fx));
%     gridpoint = [reshape(fx,prod(size(fx)),1) reshape(fy,prod(size(fx)),1)];
%     roundinggridpoint = round(gridpoint);
%     firposround = round(firposgrid);
%     firingvalue = abs(ot(firr));
%     for ii = 1:length(firposgrid)
%         [~,q1]=min(abs(firposgrid(ii,1)-fx(1,:)));
%         [~,q2]=min(abs(firposgrid(ii,2)-fy(1,:)));
%         firingmap(q1,q2) = firingvalue(ii);
%     end
%     firingmap = firingmap/max(max(firingmap));
%     gaussian = fspecial('gaussian',[10 10],3); % Design the gaussian filter for smoothing the map
%     spikes_smooth=conv2(gaussian,firingmap);
%     figure(1); imagesc(imrotate(spikes_smooth/max(max(spikes_smooth)),90)); axis off
%     colormap(jet)
%     title(['Rate map of neuron ',num2str(n)]);
%  
% % Autocorrelation
% Rxy = correlation_map(spikes_smooth,spikes_smooth);
% figure(2); imagesc(Rxy); axis off; colormap(jet)
% title(['Autocorrelation map of neuron ',num2str(n)]);
% % Grid scale computation
% % Computing the central peak index of the autocorrelation map
% c = Rxy;
% [s1,s2] = size(c);
% centralpeakindex = 0.5*[s1+1 s2+1]; cpx = centralpeakindex(1); cpy = centralpeakindex(2);
% cc=imregionalmax(c);
% [peakx,peaky] = find(cc ==1); peakmatrix = [peakx peaky];
% centralpeakrepmat = repmat(centralpeakindex,length(peakx),1);
% % Minimum and median grid scale computation
% distfromcentralpeak = diag(pdist2(peakmatrix,centralpeakrepmat));
% if distfromcentralpeak == 0
%     continue;
% end
% [distascending,index] = sort(distfromcentralpeak);
% gridscale = distascending(2); %Minimum grid scale
% mediangridscale = median(distascending(2:min(6,length(distascending)))); %Medial grid scale
% % Grid score computation 
% c=Rxy;
% radius1 = max(distascending(2:min(6,length(distascending)))) + 3;
% [p1 q1] = meshgrid(1:length(c));
% circlemask1 = sqrt((p1-centralpeakindex(1)).^2 + (q1-centralpeakindex(2)).^2) <= radius1;
% c2 = c; c2 = c2.*circlemask1;
% figure(3); imagesc(c2)
% colormap(jet)
% title('Masking the outer peaks of the autocorrelation map')
% axis off
% radius2 = 12;
% [p2 q2] = meshgrid(1:length(c));
% circlemask2 = sqrt((p2-centralpeakindex(1)).^2 + (q2-centralpeakindex(2)).^2) <= radius2;
% circlemask2 = not(circlemask2);
% c3 = c2.*circlemask2;
% figure(4); imagesc(c3)
% colormap(jet)
% title('Masking the central peak of the autocorrelation map')
% axis off
% hgs = gridscore(c3); %Hexagonal Gridness Score of the neuron
% sgs = squaregridscore(c3); %Square Gridness Score of the neuron
% if hgs > 0 && sgs < 0
%     hex = hex + 1;
%     n
% end
% if sgs > 0 && hgs < 0
%     squ = squ + 1;
% end
% figure(5);plot(-1:0.1:1,zeros(21),'k',zeros(21),-1:0.1:1,'k',hgs,sgs,'*r','markersize', 10);%hold on;
% title('HGS and SGS scores','FontSize',20);
% hold on;
% %pause;
% end
% end