aguirreg

function script_equivalence_to_exponential_models()
%
%   This simulation aims to investigate the equivalence between the
%   drifting norm models and exponential decay models.
%   A number of time-series (data) is generated according to an exponential
%   integration model, and the drift covariate that best explains the data
%   is obtained for each generated dataset.
%   Finally, a curve displaying the match is plotted with the results
%   _______________________________________________
%   by Marcelo G Mattar (10/21/2013)
%   mattar@sas.upenn.edu

%% PARAMETERS
% Free-parameters of the model
adaptFactor = 1; % Magnitude of adaptation effects relative to direct effects. For large adaptFactor, adaptation will have large influence relative to main-effects, and vice-versa
convolve_data_and_covariates = false;
saveBool = false;
outputDir = './Figures/'; % Directory where output images will be saved

% Simulation and plotting parameters - do not affect the results
decayRange = [-1 1]; % Range of decay factors to test (10^decayFactor)
elasticity_coefficients_to_test = 0:0.01:1;
number_of_decay_factors_to_test = 1000;

%% INITIALIZE VARIABLES
baseDir = '/jet/aguirre/FaceSpace/';
REFsdir = [baseDir 'MDS_RGI/PrepScripts/REFs/'];
seqFile = [REFsdir 'MDS_RGI_Seq_Numerical.ref'];
stimSeq = [2,23,6,12,19,4,20,26,3,25,1,5,18,10,24,27,16,9,22,17,13,15,7,21,0,0,14,8,11,11,24,15,16,7,25,12,17,21,4,8,1,14,19,20,6,27,9,5,13,18,22,0,0,2,3,23,26,10,10,9,8,17,16,14,12,15,24,22,2,6,11,18,7,0,0,13,4,23,27,21,25,20,3,19,5,26,1,1,12,27,5,19,10,20,15,0,0,9,7,18,21,26,17,4,13,25,22,24,3,14,16,8,6,2,11,23,23,25,17,2,14,9,11,6,8,0,0,10,12,4,19,13,20,1,7,26,5,22,18,3,24,16,27,15,21,21,23,14,11,10,13,12,2,5,17,9,24,0,0,15,22,7,16,6,19,3,4,25,26,20,8,27,1,18,18,1,6,23,10,27,22,3,12,16,24,19,0,0,7,17,15,9,21,14,13,5,11,4,2,26,8,20,25,25,13,22,8,12,7,14,1,21,6,4,27,19,17,20,5,9,2,18,23,3,26,24,11,15,10,16,0,0,0,0,21,27,7,19,18,24,13,6,20,11,14,22,25,23,1,4,12,8,26,16,10,3,9,15,2,17,5,5,15,17,7,22,11,2,8,21,3,16,23,18,19,9,12,6,0,0,4,26,25,14,24,20,10,1,13,27,27,0,0,8,14,17,26,11,20,23,12,3,5,2,25,18,16,13,7,4,10,6,24,1,9,19,22,21,15,15,20,18,25,5,4,3,10,17,22,27,26,9,23,16,1,11,13,24,21,2,12,0,0,19,6,14,7,8,8,2,24,26,6,3,0,0,12,13,23,19,11,25,15,14,4,5,1,20,27,18,9,17,10,22,16,21,7,7,11,1,0,0,22,23,13,26,12,20,9,4,6,16,17,3,2,19,21,18,5,14,10,25,8,15,27,24,24,10,5,7,9,13,16,3,8,25,21,22,15,11,12,14,18,0,0,17,1,27,6,26,2,20,19,23,4,4,15,13,17,25,7,1,8,3,20,2,16,11,5,10,0,0,23,24,14,6,18,26,22,19,27,12,21,9,9,18,2,15,1,10,4,14,21,5,0,0,11,27,20,17,19,12,22,13,8,24,23,7,3,6,25,16,26,26,18,27,3,1,24,9,20,7,5,16,12,11,8,22,14,2,21,17,6,13,10,23,15,4,0,0,25,19,19,25,4,11,26,13,2,10,14,27,17,18,15,5,21,1,3,22,9,0,0,16,20,12,23,8,7,24,6,6,1,25,3,13,9,10,26,15,8,4,18,17,14,23,21,12,24,2,7,20,0,0,5,27,11,19,16,22,22,1,26,27,8,19,24,5,20,16,4,9,14,15,23,17,11,3,21,10,2,0,0,18,12,25,6,7,13,13,14,5,8,23,2,4,24,25,27,10,19,26,7,6,9,3,18,11,17,0,0,1,22,20,21,16,15,12,12,26,21,20,24,17,8,18,4,1,15,19,7,2,13,3,27,14,25,0,0,6,5,23,22,10,11,9,16,16,19,2,1,23,9,26,0,0,27,25,10,7,15,6,17,24,12,18,14,20,4,21,8,13,11,22,5,3,3,15,25,2,27,23,5,6,22,26,4,16,18,20,13,19,14,0,0,24,8,10,21,11,7,12,9,1,17,17,27,4,22,12,5,24,18,13,0,0,3,7,10,8,16,25,11,21,19,1,2,9,6,15,26,23,20,14,14,26,19,15,3,17,12,1,16,2,22,6,10,18,8,5,25,9,27,13,21,24,4,7,23,11,0,0,20,20,22,4,17,23,0,0,26,14,3,11,16,5,12,10,15,18,6,21,13,1,19,8,9,25,24,7,27,2,2,23,6,12,19,4,20,26,3,25,1,5,18,10,24,27,16,9,22,17,13,15,7,21,0,0,14,8,11,11,24,15,16,7,25,12,17,21,4,8,1,14,19,20,6,27,9,5,13,18,22,0,0,2,3,23,26,10,10,9,8,17,16,14,12,15,24,22,2,6,11,18,7,0,0,13,4,23,27,21,25,20,3,19,5,26,1,1,12,27,5,19,10,20,15,0,0,9,7,18,21,26,17,4,13,25,22,24,3,14,16,8,6,2,11,23,23,25,17,2,14,9,11,6,8,0,0,10,12,4,19,13,20,1,7,26,5,22,18,3,24,16,27,15,21,21,23,14,11,10,13,12,2,5,17,9,24,0,0,15,22,7,16,6,19,3,4,25,26,20,8,27,1,18,18,1,6,23,10,27,22,3,12,16,24,19,0,0,7,17,15,9,21,14,13,5,11,4,2,26,8,20,25,25,13,22,8,12,7,14,1,21,6,4,27,19,17,20,5,9,2,18,23,3,26,24,11,15,10,16,0,0,0,0,21,27,7,19,18,24,13,6,20,11,14,22,25,23,1,4,12,8,26,16,10,3,9,15,2,17,5,5,15,17,7,22,11,2,8,21,3,16,23,18,19,9,12,6,0,0,4,26,25,14,24,20,10,1,13,27,27,0,0,8,14,17,26,11,20,23,12,3,5,2,25,18,16,13,7,4,10,6,24,1,9,19,22,21,15,15,20,18,25,5,4,3,10,17,22,27,26,9,23,16,1,11,13,24,21,2,12,0,0,19,6,14,7,8,8,2,24,26,6,3,0,0,12,13,23,19,11,25,15,14,4,5,1,20,27,18,9,17,10,22,16,21,7,7,11,1,0,0,22,23,13,26,12,20,9,4,6,16,17,3,2,19,21,18,5,14,10,25,8,15,27,24,24,10,5,7,9,13,16,3,8,25,21,22,15,11,12,14,18,0,0,17,1,27,6,26,2,20,19,23,4,4,15,13,17,25,7,1,8,3,20,2,16,11,5,10,0,0,23,24,14,6,18,26,22,19,27,12,21,9,9,18,2,15,1,10,4,14,21,5,0,0,11,27,20,17,19,12,22,13,8,24,23,7,3,6,25,16,26,26,18,27,3,1,24,9,20,7,5,16,12,11,8,22,14,2,21,17,6,13,10,23,15,4,0,0,25,19,19,25,4,11,26,13,2,10,14,27,17,18,15,5,21,1,3,22,9,0,0,16,20,12,23,8,7,24,6,6,1,25,3,13,9,10,26,15,8,4,18,17,14,23,21,12,24,2,7,20,0,0,5,27,11,19,16,22,22,1,26,27,8,19,24,5,20,16,4,9,14,15,23,17,11,3,21,10,2,0,0,18,12,25,6,7,13,13,14,5,8,23,2,4,24,25,27,10,19,26,7,6,9,3,18,11,17,0,0,1,22,20,21,16,15,12,12,26,21,20,24,17,8,18,4,1,15,19,7,2,13,3,27,14,25,0,0,6,5,23,22,10,11,9,16,16,19,2,1,23,9,26,0,0,27,25,10,7,15,6,17,24,12,18,14,20,4,21,8,13,11,22,5,3,3,15,25,2,27,23,5,6,22,26,4,16,18,20,13,19,14,0,0,24,8,10,21,11,7,12,9,1,17,17,27,4,22,12,5,24,18,13,0,0,3,7,10,8,16,25,11,21,19,1,2,9,6,15,26,23,20,14,14,26,19,15,3,17,12,1,16,2,22,6,10,18,8,5,25,9,27,13,21,24,4,7,23,11,0,0,20,20,22,4,17,23,0,0,26,14,3,11,16,5,12,10,15,18,6,21,13,1,19,8,9,25,24,7,27,2]';
% these were taken directly from 'Eigenvectors and values for one perceptual matrix.xls'
coordinatelist = [NaN NaN NaN
    -0.43815233	0.043993216	0.000988356
    -0.16477091	-0.048723129	-0.090894096
    -0.24489513	-0.16840255	0.18525238
    -0.45287893	0.11428962	-0.009152772
    -0.1665505	-0.004077156	-0.10680729
    -0.1814501	-0.15005818	0.15631847
    -0.21176483	0.44336182	0.13605036
    0.10204564	0.37570476	0.03404298
    0.04136971	0.18141115	0.23324493
    -0.2850246	-0.020309502	-0.10849088
    0.032296558	-0.18105101	-0.14802057
    0.024659614	-0.30407585	0.13620073
    -0.29424431	0.051400923	-0.088352834
    0.12463116	-0.1310324	-0.10569481
    0.027527177	-0.29119301	0.18201885
    -0.061689789	0.39513857	-0.012518304
    0.30989956	0.22948687	-0.020096972
    0.31145212	0.10152802	0.1321107
    -0.089943402	-0.14891391	-0.15344147
    0.23075579	-0.21014674	-0.16836819
    0.14466016	-0.31300308	0.15095023
    -0.07980242	-0.053859279	-0.19630964
    0.23728748	-0.15372509	-0.14647582
    0.14153114	-0.31762302	0.058781603
    0.20348456	0.26825272	-0.067086623
    0.40232804	0.23268949	-0.11921149
    0.33723857	0.058936748	0.13496217];

decayFactor_list = logspace(decayRange(1),decayRange(2),number_of_decay_factors_to_test);
best_elasticity = zeros(size(decayFactor_list));
directEffect = 1; % Magnitude of a stimulus independent of adaptation effects

for d = 1:length(decayFactor_list)
    
    decayFactor = decayFactor_list(d);
    
    %% GENERATE DATA (TIME-SERIES)
    data = zeros(size(stimSeq));
    seqLen = length(stimSeq);
    
    % Create a sequence of 3-dimensional coordinates
    coordSeq = [coordinatelist(stimSeq+1,1) coordinatelist(stimSeq+1,2) coordinatelist(stimSeq+1,3)];
    
    % Construct the data, trial by trial
    for t=1:seqLen
        cummulativeAdaptation = 0;
        for i=1:(t-1)
            % temporalDecay is a multiplicative factor that attenuates
            % effects further back into the past
            temporalDecay = exp(decayFactor*(1+i-t));
            
            % distance is nan if it involves a blank
            distance = sqrt((coordSeq(t,1)-coordSeq(i,1))^2 + (coordSeq(t,2)-coordSeq(i,2))^2 + (coordSeq(t,3)-coordSeq(i,3))^2);
            % adaptEffect relates perceptual distance to adaptation
            adaptEffect = distance;
            
            %cummAdapt shows the cummulative adaptation over time
            if ~isnan(adaptEffect)
                cummulativeAdaptation = cummulativeAdaptation + temporalDecay*adaptEffect;
                % Notice that trials following a blank will have near-zero
                % cummulative adaptation, leading to small signal
            end
        end
        

        if stimSeq(t)==0
            data(t) = adaptFactor*cummulativeAdaptation;
            % However, in this case, cummulativeAdaptation=0, so data=0
        else
            data(t) = directEffect + adaptFactor*cummulativeAdaptation;
            % Notice the plus sign: large distances -> large recovery from
            % adaptation
        end
    end
    
    % Convolve simulated data with HRF
    if convolve_data_and_covariates
        data = downsample(hrfconv(data,1500),(length(hrfconv(data,1500)))/(length(data)));
    end
    
    
    %% FIND OPTIMAL ELASTICITY COEFFICIENT
    elasticity_coefficients = elasticity_coefficients_to_test;
    correlations_found = zeros(size(elasticity_coefficients));
    
    for i=1:length(elasticity_coefficients)
        drift_covariate = driftCov_simple(round(100*elasticity_coefficients(i)));
        
        % Convolve simulated data with HRF
        if convolve_data_and_covariates
            drift_covariate = downsample(hrfconv(drift_covariate,1500),(length(hrfconv(drift_covariate,1500)))/(length(drift_covariate)));
        end
        
        % Find correlation between data and covariate
        correlations_found(i) = corr2(drift_covariate,data);
    end
    
    [~, I] = max(correlations_found);
    best_elasticity(d) = elasticity_coefficients(I);
    
end

%% FINDS THE EXPONENTIAL CURVE CORRESPONDING TO ELASTICITY=0.40
[~,X] = min(abs(best_elasticity-0.4));
decayMatch = decayFactor_list(X);


%% PLOT THE RESULTS
figure(1);
hold on;
plot(1:0.01:11,adaptFactor*exp(10^decayRange(1)*(0:-0.01:-10)))
plot(1:0.01:11,adaptFactor*exp(10^(mean([decayRange(1),decayRange(2)]))*(0:-0.01:-10)),'g')
plot(1:0.01:11,adaptFactor*exp(10^decayRange(2)*(0:-0.01:-10)),'r')
plot(1:0.01:11,adaptFactor*exp(decayMatch*(0:-0.01:-10)),'k-')
title('Adaptation curves');
ylabel('Relative magnitude of attenuation');
xlabel('Number of trials back');
legend({['10^' num2str(decayRange(1))],['10^' num2str(mean([decayRange(1),decayRange(2)]))],['10^' num2str(decayRange(2))],'Matched adaptation to ec=0.4'});
ylim([0 1]);
xlim([1 10]);
if saveBool; print('-depsc',[outputDir 'Sample_adaptation_curves.eps']); end;
if saveBool; print('-dpng',[outputDir 'Sample_adaptation_curves.png']); end;

figure(2);
semilogx(decayFactor_list,best_elasticity);
title('Simulating the experiment sequence');
ylabel('Elasticity Coefficients');
xlabel('Adaptation curve decay factor');
ylim([0 1]);
h1=line([decayMatch decayMatch],get(gca,'YLim'),'Color','r');%h is the *handle*
set(h1,'LineStyle',':');%dash style
h2=line(get(gca,'XLim'),[0.4,0.4],'Color','r');%h is the *handle*
set(h2,'LineStyle',':');%dash style
if saveBool; print('-depsc',[outputDir 'Elasticity-decay_correspondence.eps']); end;
if saveBool; print('-dpng',[outputDir 'Elasticity-decay_correspondence.png']); end;

%{
figure(3);
plot(1:0.01:11,adaptFactor*exp(decayMatch*(0:-0.01:-10)))
title('Adaptation curve matching FFA optimal elasticity');
ylabel('Relative magnitude of attenuation');
xlabel('Number of trials back');
if saveBool; print('-depsc',[outputDir 'Matched_adaptation_curve.eps']); end;
if saveBool; print('-dpng',[outputDir 'Matched_adaptation_curve.png']); end;
ylim([0 1]);
xlim([1 10]);
%}
end












function [evokedSignal stimContextSeq coordSeq stimSeq] = driftCov_simple(adaptFactor)
%
%	adaptFactor should be between 0 and 100
%   _______________________________________________
%   by Marcelo G Mattar (05/24/2013)
%   mattar@sas.upenn.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INITIALIZE ALL VARIABLES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

baseDir = '/jet/aguirre/FaceSpace/';

stimSeq = [2,23,6,12,19,4,20,26,3,25,1,5,18,10,24,27,16,9,22,17,13,15,7,21,0,0,14,8,11,11,24,15,16,7,25,12,17,21,4,8,1,14,19,20,6,27,9,5,13,18,22,0,0,2,3,23,26,10,10,9,8,17,16,14,12,15,24,22,2,6,11,18,7,0,0,13,4,23,27,21,25,20,3,19,5,26,1,1,12,27,5,19,10,20,15,0,0,9,7,18,21,26,17,4,13,25,22,24,3,14,16,8,6,2,11,23,23,25,17,2,14,9,11,6,8,0,0,10,12,4,19,13,20,1,7,26,5,22,18,3,24,16,27,15,21,21,23,14,11,10,13,12,2,5,17,9,24,0,0,15,22,7,16,6,19,3,4,25,26,20,8,27,1,18,18,1,6,23,10,27,22,3,12,16,24,19,0,0,7,17,15,9,21,14,13,5,11,4,2,26,8,20,25,25,13,22,8,12,7,14,1,21,6,4,27,19,17,20,5,9,2,18,23,3,26,24,11,15,10,16,0,0,0,0,21,27,7,19,18,24,13,6,20,11,14,22,25,23,1,4,12,8,26,16,10,3,9,15,2,17,5,5,15,17,7,22,11,2,8,21,3,16,23,18,19,9,12,6,0,0,4,26,25,14,24,20,10,1,13,27,27,0,0,8,14,17,26,11,20,23,12,3,5,2,25,18,16,13,7,4,10,6,24,1,9,19,22,21,15,15,20,18,25,5,4,3,10,17,22,27,26,9,23,16,1,11,13,24,21,2,12,0,0,19,6,14,7,8,8,2,24,26,6,3,0,0,12,13,23,19,11,25,15,14,4,5,1,20,27,18,9,17,10,22,16,21,7,7,11,1,0,0,22,23,13,26,12,20,9,4,6,16,17,3,2,19,21,18,5,14,10,25,8,15,27,24,24,10,5,7,9,13,16,3,8,25,21,22,15,11,12,14,18,0,0,17,1,27,6,26,2,20,19,23,4,4,15,13,17,25,7,1,8,3,20,2,16,11,5,10,0,0,23,24,14,6,18,26,22,19,27,12,21,9,9,18,2,15,1,10,4,14,21,5,0,0,11,27,20,17,19,12,22,13,8,24,23,7,3,6,25,16,26,26,18,27,3,1,24,9,20,7,5,16,12,11,8,22,14,2,21,17,6,13,10,23,15,4,0,0,25,19,19,25,4,11,26,13,2,10,14,27,17,18,15,5,21,1,3,22,9,0,0,16,20,12,23,8,7,24,6,6,1,25,3,13,9,10,26,15,8,4,18,17,14,23,21,12,24,2,7,20,0,0,5,27,11,19,16,22,22,1,26,27,8,19,24,5,20,16,4,9,14,15,23,17,11,3,21,10,2,0,0,18,12,25,6,7,13,13,14,5,8,23,2,4,24,25,27,10,19,26,7,6,9,3,18,11,17,0,0,1,22,20,21,16,15,12,12,26,21,20,24,17,8,18,4,1,15,19,7,2,13,3,27,14,25,0,0,6,5,23,22,10,11,9,16,16,19,2,1,23,9,26,0,0,27,25,10,7,15,6,17,24,12,18,14,20,4,21,8,13,11,22,5,3,3,15,25,2,27,23,5,6,22,26,4,16,18,20,13,19,14,0,0,24,8,10,21,11,7,12,9,1,17,17,27,4,22,12,5,24,18,13,0,0,3,7,10,8,16,25,11,21,19,1,2,9,6,15,26,23,20,14,14,26,19,15,3,17,12,1,16,2,22,6,10,18,8,5,25,9,27,13,21,24,4,7,23,11,0,0,20,20,22,4,17,23,0,0,26,14,3,11,16,5,12,10,15,18,6,21,13,1,19,8,9,25,24,7,27,2,2,23,6,12,19,4,20,26,3,25,1,5,18,10,24,27,16,9,22,17,13,15,7,21,0,0,14,8,11,11,24,15,16,7,25,12,17,21,4,8,1,14,19,20,6,27,9,5,13,18,22,0,0,2,3,23,26,10,10,9,8,17,16,14,12,15,24,22,2,6,11,18,7,0,0,13,4,23,27,21,25,20,3,19,5,26,1,1,12,27,5,19,10,20,15,0,0,9,7,18,21,26,17,4,13,25,22,24,3,14,16,8,6,2,11,23,23,25,17,2,14,9,11,6,8,0,0,10,12,4,19,13,20,1,7,26,5,22,18,3,24,16,27,15,21,21,23,14,11,10,13,12,2,5,17,9,24,0,0,15,22,7,16,6,19,3,4,25,26,20,8,27,1,18,18,1,6,23,10,27,22,3,12,16,24,19,0,0,7,17,15,9,21,14,13,5,11,4,2,26,8,20,25,25,13,22,8,12,7,14,1,21,6,4,27,19,17,20,5,9,2,18,23,3,26,24,11,15,10,16,0,0,0,0,21,27,7,19,18,24,13,6,20,11,14,22,25,23,1,4,12,8,26,16,10,3,9,15,2,17,5,5,15,17,7,22,11,2,8,21,3,16,23,18,19,9,12,6,0,0,4,26,25,14,24,20,10,1,13,27,27,0,0,8,14,17,26,11,20,23,12,3,5,2,25,18,16,13,7,4,10,6,24,1,9,19,22,21,15,15,20,18,25,5,4,3,10,17,22,27,26,9,23,16,1,11,13,24,21,2,12,0,0,19,6,14,7,8,8,2,24,26,6,3,0,0,12,13,23,19,11,25,15,14,4,5,1,20,27,18,9,17,10,22,16,21,7,7,11,1,0,0,22,23,13,26,12,20,9,4,6,16,17,3,2,19,21,18,5,14,10,25,8,15,27,24,24,10,5,7,9,13,16,3,8,25,21,22,15,11,12,14,18,0,0,17,1,27,6,26,2,20,19,23,4,4,15,13,17,25,7,1,8,3,20,2,16,11,5,10,0,0,23,24,14,6,18,26,22,19,27,12,21,9,9,18,2,15,1,10,4,14,21,5,0,0,11,27,20,17,19,12,22,13,8,24,23,7,3,6,25,16,26,26,18,27,3,1,24,9,20,7,5,16,12,11,8,22,14,2,21,17,6,13,10,23,15,4,0,0,25,19,19,25,4,11,26,13,2,10,14,27,17,18,15,5,21,1,3,22,9,0,0,16,20,12,23,8,7,24,6,6,1,25,3,13,9,10,26,15,8,4,18,17,14,23,21,12,24,2,7,20,0,0,5,27,11,19,16,22,22,1,26,27,8,19,24,5,20,16,4,9,14,15,23,17,11,3,21,10,2,0,0,18,12,25,6,7,13,13,14,5,8,23,2,4,24,25,27,10,19,26,7,6,9,3,18,11,17,0,0,1,22,20,21,16,15,12,12,26,21,20,24,17,8,18,4,1,15,19,7,2,13,3,27,14,25,0,0,6,5,23,22,10,11,9,16,16,19,2,1,23,9,26,0,0,27,25,10,7,15,6,17,24,12,18,14,20,4,21,8,13,11,22,5,3,3,15,25,2,27,23,5,6,22,26,4,16,18,20,13,19,14,0,0,24,8,10,21,11,7,12,9,1,17,17,27,4,22,12,5,24,18,13,0,0,3,7,10,8,16,25,11,21,19,1,2,9,6,15,26,23,20,14,14,26,19,15,3,17,12,1,16,2,22,6,10,18,8,5,25,9,27,13,21,24,4,7,23,11,0,0,20,20,22,4,17,23,0,0,26,14,3,11,16,5,12,10,15,18,6,21,13,1,19,8,9,25,24,7,27,2]';
% these were taken directly from 'Eigenvectors and values for one perceptual
% matrix.xls'
coordinatelist = [NaN NaN NaN
    -0.43815233	0.043993216	0.000988356
    -0.16477091	-0.048723129	-0.090894096
    -0.24489513	-0.16840255	0.18525238
    -0.45287893	0.11428962	-0.009152772
    -0.1665505	-0.004077156	-0.10680729
    -0.1814501	-0.15005818	0.15631847
    -0.21176483	0.44336182	0.13605036
    0.10204564	0.37570476	0.03404298
    0.04136971	0.18141115	0.23324493
    -0.2850246	-0.020309502	-0.10849088
    0.032296558	-0.18105101	-0.14802057
    0.024659614	-0.30407585	0.13620073
    -0.29424431	0.051400923	-0.088352834
    0.12463116	-0.1310324	-0.10569481
    0.027527177	-0.29119301	0.18201885
    -0.061689789	0.39513857	-0.012518304
    0.30989956	0.22948687	-0.020096972
    0.31145212	0.10152802	0.1321107
    -0.089943402	-0.14891391	-0.15344147
    0.23075579	-0.21014674	-0.16836819
    0.14466016	-0.31300308	0.15095023
    -0.07980242	-0.053859279	-0.19630964
    0.23728748	-0.15372509	-0.14647582
    0.14153114	-0.31762302	0.058781603
    0.20348456	0.26825272	-0.067086623
    0.40232804	0.23268949	-0.11921149
    0.33723857	0.058936748	0.13496217];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CREATE CONTEXT SEQUENCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = adaptFactor/100; % k is the drifting factor, now scaled between 0 and 1

% Create a sequence of 3-dimensional coordinates
coordSeq = [coordinatelist(stimSeq+1,1) coordinatelist(stimSeq+1,2) coordinatelist(stimSeq+1,3)];
seqLen = size(coordSeq,1);
stimContextSeq = zeros(size(coordSeq));

% Start in the state [0 0 0]
stimContextSeq(1,:) = [0 0 0];

% Now, drift in the perceptual space depending on presented stimulus
for i=2:seqLen
    % When the current trial is a NaN (blank) trial
    if any(isnan(coordSeq(i-1,:)))
        % Model in which the context drifts back to the center (prototype)
        %stimContextSeq(i,:) = stimContextSeq(i-1,:) + k*([0 0 0] - stimContextSeq(i-1,:));
        % Model in which the contexts stays where it was
        stimContextSeq(i,:) = stimContextSeq(i-1,:);
    else % When the current trial is a regular (stimulus) trial
        stimContextSeq(i,:) = stimContextSeq(i-1,:) + k*(coordSeq(i-1,:) - stimContextSeq(i-1,:));
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CREATE EVOKED SIGNAL SEQUENCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

coordSeq_noNaNs = coordSeq;
coordSeq_noNaNs(isnan(coordSeq)) = 0;
evokedSignal3D = abs(coordSeq_noNaNs - stimContextSeq);
evokedSignal = sqrt(evokedSignal3D(:,1).^2 + evokedSignal3D(:,2).^2 + evokedSignal3D(:,3).^2);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MEAN-CENTER AND NORMALIZE TO HAVE VARIANCE 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mean center
evokedSignal = evokedSignal - mean(evokedSignal);
% Correct blank trials to have value of zero
evokedSignal(isnan(coordSeq(:,1))) = 0;

% Normalize evokedSignal
evokedSignal = evokedSignal / std(evokedSignal);

end