%%------------------------------------------------
% Main script for generating a population of artifical spike trains
%%------------------------------------------------
% Input (in parameters cell, below): 
%   -> Population statistics for biological data that you are trying to
%   reproduce in an artifical spike train population (matrix with columns: ??)
%   -> At least one PC spike train that does not have behavioral
%   modifications 
%   -> Peri-stimulus time histogram for behavioral response 
% Output: 
%   -> A collection of files is saved to the current directory, one for
%   each AST. Filenames are gammaspikes_{cellnum}_{ASTnum}.txt and content is spiketimes in seconds 
%   -> Population statistics are compared between biological data (input)
%   and the produced AST population 
%Updated on 6/2/2017 by Samira Abbasi
%% MAIN Parameters
clear all; close all, clc

BehModFrac = 1;
BehModStren = 1;
ShiftFrac = [0];

for BMF=1:length(BehModFrac)
 % portion of ASTs with behavioral modulation
    ASTbehmodfraction=BehModFrac(BMF);
for BMS=1:length(BehModStren)
    BehModStrength = BehModStren(BMS);
for ShF=1:length(ShiftFrac)
% portion of final signal based on shifted signal 
shiftfraction=ShiftFrac(ShF);
% number of ASTs to produce
numAST = 50;
savedir = strcat('AST','\','Shift',num2str(shiftfraction),'_BehMod',num2str(ASTbehmodfraction),'_BehModStrength',num2str(BehModStrength));
mkdir(savedir);
addpath(savedir);
% bio population stats (numASTs, FR, LV, CVS, CVF, behmodfraction)
% popstat = makepoptargets(numAST, 66.82, 0.47, 0.31, 0.17, ASTbehmodfraction); %based on "representative" bio data #27
% popstat = makepoptargets(numAST, 80, 0.47, 0.24, 0.17, ASTbehmodfraction); %based on "representative" bio data #28
% popstat = makepoptargets(numAST, 84.76, 0.65, 0.34, 0.18, ASTbehmodfraction); %based on "representative" bio data #29
popstat = makepoptargets(numAST, 64.86, 0.5, 0.2, 0.19, ASTbehmodfraction); %based on "representative" bio data #30
% popstat = makepoptargets(numAST, 'dist', 0.5, 0.2, 0.19, ASTbehmodfraction);
% popstat = makepoptargets(numAST, 37.82, 0.86, 0.46, 0.34, ASTbehmodfraction); %based on "representative" bio data #31
% popstat = makepoptargets(numAST, 59.92, 0.28, 0.16, 0.08, ASTbehmodfraction); %based on "representative" bio data #32
% popstat = makepoptargets(numAST, 46.59, 0.51, 0.26, 0.15, ASTbehmodfraction); %based on "representative" bio data #34
% popstat = makepoptargets(numAST, 71.34, 0.94, 0.26, 0.21, ASTbehmodfraction); %based on "representative" bio data #37
% popstat = makepoptargets(numAST, 77.48, 0.43, 0.19, 0.21, ASTbehmodfraction); %based on "representative" bio data #42
% popstat = makepoptargets(numAST, 136.1, 0.46, 0.25, 0.24, ASTbehmodfraction); %based on "representative" bio data #43
% popstat = makepoptargets(numAST, 102, 0.31, 0.18, 0.08, ASTbehmodfraction); %based on "representative" bio data #44
% popstat = makepoptargets(numAST, 134.46, 0.36, 0.24, 0.08, ASTbehmodfraction); %based on "representative" bio data #45
% popstat = makepoptargets(numAST, 64.5, 0.44, 0.22, 0.13, ASTbehmodfraction);%bio data 46
% popstat = makepoptargets(numAST, 167.95, 0.48, 0.18, 0.12, ASTbehmodfraction);%bio data 47
% popstat = makepoptargets(numAST, 100.44, 0.56, 0.27, 0.23, ASTbehmodfraction);%bio data 48
% popstat = makepoptargets(numAST, 89.39, 0.74, 0.3, 0.26, ASTbehmodfraction); %based on "representative" bio data #49
% popstat = makepoptargets(numAST, 83.1, 0.49, 0.23, 0.13, ASTbehmodfraction); %based on "representative" bio data #51
% popstat = makepoptargets(numAST, 143.15, 0.67, 0.26, 0.31, ASTbehmodfraction); %based on "representative" bio data #52
% popstat = makepoptargets(numAST, 105.23, 0.89, 0.24, 0.12, ASTbehmodfraction); %based on "representative" bio data #82
% popstat = makepoptargets(numAST, 88.94, 0.73, 0.24, 0.12, ASTbehmodfraction); %based on "representative" bio data #83
% popstat = makepoptargets(numAST, 128.91, 0.67, 0.2, 0.1, ASTbehmodfraction); %based on "representative" bio data #86
%popstat = makepoptargets(numAST, 'dist', 'dist', 'dist','dist', ASTbehmodfraction); %bio based ranges
%popstat = makepoptargets(numAST, [50 120], 0.47, 0.24, 0.17, 0); %FR sweep
%popstat = load('aststats_HTsweep.txt'); % for hand-tuning population stats
save('astno_rate_lv_cvsl_cvft_psth.txt','popstat','-ascii');

% time units 
bio_t=0.0001; % conversion factor for input bio data - into seconds
desired_t=0.001; % desired units of ASTs (aka conversion factor for AST spiketimes to seconds)
win = 450;
order = 30;
% refractory period in seconds
refper = 0.003; 

% Which algorithm to use? 1:Selva (make sure filter=4) 2:UngammaForwardLooking (Filter=1) 3:GammaForwardLooking (filter=4)
algorithmflag=2;
% Smoothing of adaptive Gaussian rate estimate
filter = 1;

% plotting? y/n 1/0
plotflag = 0; 

%% Make ASTs 
% i=recorded PC number
for i= [30]
    %% Load bio data 
    
    % example bio spike times for initial rate estimate 
    fname=data_spt(i);
    sptimes=load(strcat(fname,'_refp3.txt'));   
        
    % corresponding event times 
    bevtfname=data_insp(i); %or data_exp(i);
    bevt=load(bevtfname)*bio_t; % convert to seconds  
  
    %% Get rate template from bio data 
    
    % Get fixed Gaussian estimate
    [fixedwinRate] = fixedGauss_FRest(sptimes,bio_t,plotflag,'.b'); 
       % Get adaptive Gaussian estimate
    [adapwinRate] = adaptGauss_FRest(sptimes,fixedwinRate,bio_t,filter,plotflag,'.r');
    
    % Loop for AST creation, given number of ASTs per base rate template     
    for j=1:numAST
%         Grab statistical targets for AST
        outputmean = popstat(j,2);
        reg = lv_kappa(popstat(j,3));
        gainsl = popstat(j,4);
        gainfa = popstat(j,5);   
        if popstat(j,6)==0 || BehModStrength==0
            psth = load('PSTH_flat.txt');
        else
            psth = load('PSTH_movavgorder_30_cellnum_43.txt');
            Mp = mean(psth);
            psth = psth - Mp;
            psth = psth * BehModStrength;
            psth = psth + Mp;
             
        end
        
%         add "noise" (shifted template) to rate template
            noisyfixrate=fixedwinRate; % these two lines were added by Samira
            noisyadaprate=adapwinRate;
        if ~(shiftfraction==0)
%             Get shifted version of rate function
            [shiftedf,shift]=templateshift(fixedwinRate(:,2)); %random shift
            [shifteda,~]=templateshift(adapwinRate(:,2),shift); %matches shift from fixedwinrate to make later processing easier 
            
%             Add in shifted signal
            noisyfixrate=fixedwinRate; %these two lines became comment by
%             Samira, because when shiftfraction=0 this is an error
            noisyadaprate=adapwinRate;
            noisyfixrate(:,2)=((1-shiftfraction)*fixedwinRate(:,2))+(shiftfraction*shiftedf);
            noisyadaprate(:,2)=((1-shiftfraction)*adapwinRate(:,2))+(shiftfraction*shifteda);
        end
%         measure CV for adjusted rate template (split by frequency)
        [adaptive,slow,templateCV_s, templateCV_f] = splitsignal_adjust(noisyfixrate,noisyadaprate,1,1,0);
    
%         adjust final templates with gain modulation 
        [FRadj,sl_FRadj,~,~] = splitsignal_adjust(noisyfixrate,noisyadaprate,gainsl/templateCV_s,gainfa/templateCV_f,shiftfraction,plotflag);


%       Create AST
        singleASTgen(savedir,ASTbehmodfraction,FRadj,bevt,outputmean,reg,refper,psth,algorithmflag,strcat(num2str(i),'_',num2str(j)));
    
    end
      
    % Compare AST pop stats to bio pop stats
%     compareASTs
    [bioStats, astStats] = bioVSast(num2cell(i),i, bio_t, numAST,savedir,1);
    savefig(strcat(savedir,'\','bar.fig'))
     close all
% 
end
end
end
end
