aguirreg

function [evokedSignal, stimContextSeq, coordSeq, stimSeq] = driftCov(dataset, subjectID, elasticityCoefficient, saveBool)
%
%	Create drift covariate with specified elasticityCoefficient, using
%	perceptual space sequences specific to FaceSpace project
%
%   PS: elasticityCoefficient should be between 0 and 100
%   _______________________________________________
%   by Marcelo G Mattar (05/24/2013)
%   mattar@sas.upenn.edu

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    saveBool = false;
end

if nargin < 3
    elasticityCoefficient = 50;
end

% Check whether adaptFactor is integer and within acceptable values
if mod(elasticityCoefficient,1) || elasticityCoefficient<0 || elasticityCoefficient>100
    error('adaptFactor should be an integer between 0 and 100');
end


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

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

% Load the numerical sequence
switch dataset
    case 'MDS_RGI'
        REFsdir = [baseDir dataset '/PrepScripts/REFs/'];
        seqFile = [REFsdir 'MDS_RGI_Seq_Numerical.ref'];
        stimSeq = readref(seqFile,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];
    case 'MDS_Multi'
        REFsdir = [baseDir dataset '/PrepScripts/REFs/' subjectID '/'];
        seqFile = [REFsdir 'sequence_numerical.ref'];
        fid = fopen(seqFile);
        temp = textscan(fid,'%s'); % reading in as strings as the # of headerlines is variable
        temp = temp{1}; %textscan outputs a cell
        temp = temp(length(temp)-1623:end);
        for m = 1:length(temp)
            if isempty(strfind(temp{m},';'))
                sequence(m,1) = str2num(temp{m});
            end
        end
        fclose('all');
        stimSeq = sequence;
        % these were taken directly from 'Eigenvectors and values for one perceptual
        % matrix.xls'
        coordinatelist = [NaN NaN NaN
            -0.429506025	0.14859489	0.057221819
            -0.446558075	-0.072662325	0.082842149
            -0.39130837	-0.247629855	0.099019244
            -0.306944082	0.269967421	0.042639309
            -0.433483984	-0.028763763	0.066148744
            -0.093075036	-0.246233099	0.113647528
            -0.082223497	0.304004574	0.19419913
            0.184970441	0.114906441	0.338685284
            0.190919212	-0.130203479	0.237115284
            -0.370009793	0.128845431	-0.044180839
            -0.337408573	-0.147048324	-0.05693726
            -0.234586561	-0.268410842	-0.033561283
            -0.232464348	0.280904757	-0.012625924
            -0.102973475	-0.102022412	-0.065337964
            0.177664453	-0.274895784	-0.022533763
            0.241218979	0.328820781	0.062569699
            0.437924802	0.035934599	0.151015462
            0.34263309	-0.143838716	0.092513187
            -0.17387477	0.207994257	-0.244505531
            -0.181519001	-0.135689451	-0.252204393
            0.080987665	-0.281117325	-0.161772301
            0.074844511	0.302484149	-0.260652118
            0.446598907	-0.030520112	-0.129587393
            0.260651704	-0.243830992	-0.079453789
            0.388182313	0.277024364	-0.103839934
            0.538898617	0.065782352	-0.072075491
            0.450440895	-0.112397537	0.001651143];
end



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

k = elasticityCoefficient/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);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAVE COVARIATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if saveBool
    driftCovDir = [REFsdir 'driftCovs/'];
    % Create necessary directories
    if ~exist(driftCovDir, 'dir')
        mkdir(driftCovDir);
    end
    
    % Create file ID
    this_filename = [driftCovDir 'driftCov_' sprintf('%03d',elasticityCoefficient) '.ref'];
    driftCov_fileID = fopen(this_filename,'w');
    
    % Populate file
    for trial=1:(seqLen-1)
        fprintf(driftCov_fileID,'%f\n',evokedSignal(trial));
    end
    fprintf(driftCov_fileID,'%f',evokedSignal(seqLen));
    
    fclose('all');
end