% Jon Newman
% Estimate STA (unfiltered)
% 2013 03 07

function [h, stm] = make_h_3(config,T)

% Cycle through drug conditions
clear h;
for r = 1:length(config.rec_map)
    
    % Load pre-parsed data struct
    load([config.parse_path config.expID '_' config.rec_map{r} '.mat']);
    load([config.here config.corr_path config.expID '_' config.rec_map{r} '_psth.mat']);
    
    if ~isfield(data,'rep') || isempty(data.rep)
        continue;
    end
    
    % Sample period
    dt = psth.dt;
    
    % Regularization parameter
    lam = 10;
    
    % Current drug condition, formatted as valied field name
    this_drug = strrep(config.rec_map{r}, '-','_');
    
    % Find the number of trials
    flds = fields(data.rep);
    
    % Calculate STA vectors bounds
    N = ceil(T/dt);
    
    % h structure
    h.(this_drug).dt = dt;
    h.(this_drug).T = T;
    h.(this_drug).N = N;
    h.(this_drug).h = zeros(...
        numel(psth.units),...     % Number of units
        N+1);                           % Max lag in h
    
    % We will gather down-sampled stimuli while we are here.
    stm.(this_drug).dt = dt;
    stm.(this_drug).stm = zeros(...
        length(flds),...                        % Number of trials
        ceil(config.trial_length_sec/dt + 1)); 	% Number of data points per stimulus trial
    
    % Get analytical version of autocorrelation matrix
    TT = dt*(-N:N);
    ACmat = zeros(N+1,N+1);
    if strcmp(config.type,'OU')
        AC = exp(-abs(TT)/config.tau_stim);
        AC = fftshift(AC);
        for i = 1:N+1
            AC = circshift(AC,[0 1]);
            ACmat(i,:) = AC(1:N+1);
        end
    end
    
    % Cycle through trials
    %     for i = 1:length(psth.trials)
    
    % Stimulus time/value on trial i, pad stimulus waveform with zeros
    % HACK
    if dt == 0.0005
        ss_ti = interp(data.rep.(flds{1}).stm.lightDensitymWmm2,2);
        ss_ti = (ss_ti - nanmean(ss_ti))/nanstd(ss_ti);
    else
        [~,uidx] = unique(dt*round(data.rep.(flds{1}).stm.time/dt));
        ss_ti = data.rep.(flds{1}).stm.lightDensitymWmm2(uidx);
        ss_ti = (ss_ti - nanmean(ss_ti))/nanstd(ss_ti);
    end
    
    stm.(this_drug).stm(i,1:length(ss_ti)) = ss_ti;
    ss_ti = padarray(ss_ti',N)';
    
    % Cycle through units
    for j = 1:length(psth.units)
        
        % Stim/response matracies
        X = zeros(size(psth.fr_rep(j,:),2),N+1); % Holds stim snipets that preceed each response bin
        Y = squeeze(psth.fr_rep(j,:))'; % Binned firing rate
        
        % You only need to change the stim shift matrix where the cell
        % fired
        iX = find(Y>0);
        
        % If the cell fired
        if sum(Y) > 0
            
            % Stimulus shift matrix
            for k = 1:length(iX);
                X(iX(k),:) = ss_ti(iX(k):iX(k)+N); % Account for padding with the +N
            end
            
            % The regularized, whitened, STA
            sta = (T*(ACmat + lam*eye(N+1)))\((1/sum(Y))*(X'*Y));
        else
            sta = zeros(N+1,1);
        end
        
        % The STA for the ith trial, jth neuron
        h.(this_drug).h(j,:) = sta;
        
    end
    
    % Calculate h for the average network response
    % Stim/response matracies
    X = zeros(size(psth.fr_rep(j,:),2),N+1); % Holds stim snipets that preceed each response bin
    Y = psth.fr_pop.fr'; % Binned firing rate
    
    % Stimulus time/value on trial i, pad stimulus waveform with NaNs
    ss_ti = mean(stm.(this_drug).stm,1);
    stm.(this_drug).avgstm = ss_ti;
    ss_ti = padarray(ss_ti',N)';
    
    % You only need to change the stim shift matrix where the cell
    % fired
    iX = find(Y>0);
    
    % If the cell fired
    if sum(Y) > 0
        
        % Stimulus shift matrix
        for k = 1:length(iX);
            X(iX(k),:) = ss_ti(iX(k):iX(k)+N); % Account for padding with the +N
        end
        
        % The regularized, whitened, STA
        sta = (T*(ACmat + lam*eye(N+1)))\((1/sum(Y))*(X'*Y));
    end
    
    h.(this_drug).poph = sta';
    
end


end