aguirreg

function writeccocov(varargin)
% WRITECCOCOV(DEGREES,POINTSLIST,USETARGET,CUTOFFFIRST,RUNFILES,SEQUENCE, MINKEXPONENT)
% Generates covariates CCO design given a set of points and a degrees of assumed axis rotation. 
% sequence is a list of stimulus ids, including the first CUTOFFFIRST which will be discarded
% DEGREES         default = no rotation (0 deg)
% POINTSLIST      default = 4x4 linearly scaled grid
% USETARGET       default = no target
% CUTOFFFIRST     default = 10
% SEQUENCE        if specified, use this sequence instead of reading 
%
% reads sequence file, calls genccocov to generate the covs, writes ref files
% 2007 ddrucker@psych.upenn.edu
 
defaultValues={0,points(0,4),false,10,5,[],2};
nonemptyIdx = ~cellfun('isempty',varargin);
defaultValues(nonemptyIdx) = varargin(nonemptyIdx);
[degrees pointslist usetarget cutofffirst runfiles sequence minkexponent] = deal(defaultValues{:});
useseq = isempty(sequence);
for i=1:runfiles
    if useseq
        sequence = readref(['dioct-seq-' num2str(i) '-optCE.txt']);
%        sequence = readref(['dioct18-target-' num2str(i) '-optCE.txt']);
%        sequence = readref(['out' num2str(i) '.txt']);
    end
    c(i) = genccocov(sequence,degrees,pointslist,usetarget,cutofffirst, minkexponent);
end
 
 
writeref('main.ref',            cat(1, c(:).main        ));
writeref('new.ref',             cat(1, c(:).new         ));
writeref('rept.ref',            cat(1, c(:).rept        ));
writeref('directP.ref',         cat(1, c(:).directP     ));
writeref('directQ.ref',         cat(1, c(:).directQ     ));
writeref('adaptorig.ref',       cat(1, c(:).adaptorig    ));
% the following is unix-dependent. I'd really like to do this a more portable way, but Matlab, strings, and I don't get along.
!sed -i -e 's/-9/new/' adaptorig.ref
!sed -i -e 's/-5/rept/' adaptorig.ref
delete adaptorig.ref-e  % this is retarded I know; it's a osx vs linux sed difference bug thing.
writeref('adapteuorig.ref',     cat(1, c(:).adapteuorig  ));
!sed -i -e 's/-9/new/' adapteuorig.ref
!sed -i -e 's/-5/rept/' adapteuorig.ref
delete adapteuorig.ref-e
%writeref('adaptN.ref',          cat(1, c(:).adaptN       ));
%writeref('adaptNP.ref',         cat(1, c(:).adaptNP      ));
%writeref('adaptNQ.ref',         cat(1, c(:).adaptNQ      ));
writeref('adaptP.ref',          rangeone(cat(1, c(:).adaptP      )));
writeref('adaptQ.ref',          rangeone(cat(1, c(:).adaptQ      )));
writeref('adapteu.ref',         rangeone(cat(1, c(:).adapteu     )));
writeref('adapt.ref',           rangeone(cat(1, c(:).adapt       )));
writeref('adaptOrthoEuclid.ref', rangeone(cat(1,c(:).adaptOrthoEuclid)));
 
dlmwrite('sequence.ref',cat(1,c(:).sequence),'precision','%02d');
PrependHeader('sequence.ref',';VB98\n;REF1\n');
 
fid=fopen('pairings.ref','w');
for r=1:runfiles
    for i=1:length(cat(1,c(r).apd))
        fprintf(fid,'%s\n',cat(1,c(r).apd{i}));
    end
end
fclose(fid);
PrependHeader('pairings.ref',';VB98\n;REF\n');
 
fid=fopen('saulpairings.ref','w');
for r=1:runfiles
    for i=1:length(cat(1,c(r).apd))
        fprintf(fid,'%s\n',cat(1,c(r).saulapd{i}));
    end
end
fclose(fid);
PrependHeader('saulpairings.ref',';VB98\n;REF\n');
 
if usetarget  % alison wants these, i don't
    writeref('target.ref',      cat(1, c(:).target      ));
    writeref('aftertarget.ref', cat(1, c(:).aftertarget ));
    i=0;
    ddists = cat(1,c(:).ddists);
%    for angle=c(1).anglelist'
%        i=i+1;
%        writeref(sprintf('angle%d.ref',angle), mncnleavezeros(ddists(:,i)));
%    end
end

function covs = genccocov(v,varargin)
% COVS = GENCCOCOV(V,DEGREES,POINTSLIST,USETARGET,CUTOFFFIRST,MINKEXPONENT)
% Generates covariates for a continuous carryover fMRI design (Aguirre 2007) given a set of points and a degrees of assumed axis rotation. 
% V is a list of stimulus ids, including the first CUTOFFFIRST which will be discarded
% DEGREES         default = no rotation (0 deg)
% POINTSLIST      default = 4x4 linearly scaled grid
% USETARGET       default = no target
% CUTOFFFIRST     default = 10
% MINKEXPONENT    default = 2 (applies to adaptMink only)
%
% 2007-9 Daniel M. Drucker ddrucker@psych.upenn.edu
 
%% initialize default values
defaultValues={0,points(0,4),false,10,2};
nonemptyIdx = ~cellfun('isempty',varargin);
defaultValues(nonemptyIdx) = varargin(nonemptyIdx);
[degrees pointslist usetarget cutofffirst minkexponent] = deal(defaultValues{:});
 
%% set up template variables
template = zeros(length(v),1);
notzero = find(v~=0);
iszero  = find(v==0);
 
%% MAIN 
% + where a stimulus is on the screen, - for blanks --> mean-centered
main = template;
main(:) = 1;                       % main is normally +
main(iszero)=-1;                   % main is - if v is 0
 
%% NEW 
new = template;
new(notzero)=-1;                   % new is - if v is +
new(iszero+1)=1;                   % new is + following v being 0
new(iszero)=0;                     % new is 0 if v is 0, even if the previous element is 0
new = new(1:length(v));            % if last element is 0, don't lengthen
 
%% REPT
rept = template;
rept(notzero)=-1;                  % - for most stims
rept(find(diff(v)==0)+1)=1;        % 1 for repeats
rept(iszero)=0;                    % 0 for blanks
rept = rept(1:length(v));          % if last element is 0, don't lengthen
 
%% TARGET
t=0; % how many extra rows/columns do we need, besides the one for blank?
if usetarget
    target = template;
    target(v==max(v)) = 1;         % target is + if v is max
    target(v~=max(v)) = -1;        % target is - otherwise
    target(iszero) = 0;            % target is 0 if v is 0
 
    % AFTERTARGET
    aftertarget = template;
    aftertarget(notzero) = -1;                % normally -
    aftertarget(find(v==max(v))+1) = 1;       % aftertarget is + AFTER v is max
    aftertarget(iszero) = 0;                  % aftertarget is 0 if v is 0
    aftertarget = aftertarget(1:length(v));   % if last element is 0, don't lengthen
 
    % general target-related housekeeping 
    notzero = intersect(find(v~=0),find(v~=max(v)));
    t=1;  % we need an extra row/column for target
end
 
%% direct
% what are the actual direct and adapt effects?
% we need to get this information from GenTestMat
 
% note that pointslist should NOT contain the target if there is one
d = GenTestMat(pointslist,degrees,minkexponent);
 
% v is the vector containing the id number of the row/col in GenTestMat's output
 
directP = template;
directQ = template;
% rotate coords of pointslist
pointslist = rotm(pointslist,degrees)+1;
directP(notzero) = pointslist(v(notzero),1);
directQ(notzero) = pointslist(v(notzero),2);
 
%% adapt
% i need a (length(v),2) of all the sequence pairs
pairs = [];
pairs(:,1) = v(1:end-1);
pairs(:,2) = v(2:end);
pairs = [0 0 ; pairs];
% and pairs will be indexed by 1, keeping in mind that 1 actually is 0
pairs = pairs + 1;
 
% matlab, being retarded, can't deal with 0s as indices so we need to
% generate a sim matrix that's off by one
 
% mean center
d.CityOrig = d.City;
d.City  = mncnsim(d.City);
d.CityP = mncnsim(d.CityP);
d.CityQ = mncnsim(d.CityQ);
 
 
% initialize all the matrices we're going to populate
[City CityOrig CityP CityQ CityN CityNP CityNQ Mink Euclid EuclidOrig OrthoEuclid anglemat] = deal(zeros(length(d.pointslist)+t+1));
 
City(2:end-t,2:end-t)      = d.City;
CityOrig(2:end-t,2:end-t)  = d.CityOrig;
CityP(2:end-t,2:end-t)     = d.CityP;
CityQ(2:end-t,2:end-t)     = d.CityQ;
 
adapt      = diag(City(pairs(:,1),pairs(:,2)));
adaptOrig  = diag(CityOrig(pairs(:,1),pairs(:,2)));
adapt(new==1)=-9;
adaptOrig(new==1)=-9;
adapt(rept==1)=-5;
adaptOrig(rept==1)=-5;
adaptP     = diag(CityP(pairs(:,1),pairs(:,2)));
adaptQ     = diag(CityQ(pairs(:,1),pairs(:,2)));
 
CityN(2:end-t,2:end-t)  = d.CityN;
CityNP(2:end-t,2:end-t) = d.CityNP;
CityNQ(2:end-t,2:end-t) = d.CityNQ;
 
adaptN  = diag( CityN(pairs(:,1),pairs(:,2)));
adaptNP = diag(CityNP(pairs(:,1),pairs(:,2)));
adaptNQ = diag(CityNQ(pairs(:,1),pairs(:,2)));
 
%% euclidean
OrthoEuclid(2:end-t,2:end-t) = vec2sim(d.OrthoEuclidVec);
adaptOrthoEuclid = diag(OrthoEuclid(pairs(:,1),pairs(:,2)));
 
%% adaptpairs
% encodes adapt IDs for all 120 pairs of stimuli
apdraw = sort(pairs-1,2);
for i=1:length(apdraw)
    apd(i)={sprintf('%02d_%02d',apdraw(i,1),apdraw(i,2))};
end
saulapd = apd;
apd(iszero)  ={'0'};
%apd(rept==1) ={'0'};
apd(iszero+1)={'0'};
if usetarget
    apd(target==1)  ={'0'};
    apd(aftertarget==1) = {'0'};
end
apd = apd(1:length(v)); % if last element is 0, don't lengthen
 
%% direction distances
% there are length(unique(round(d.Angle90Vec))) many of these
% their value is the Euclidean distance between the pair if the pair
% is of the matching angle for that cov, 0 elsewhere, and 0 for blanks
 
% set up all of them
d.EuclidOrig     = d.Euclid;
d.Euclid = mncnsim(d.Euclid);
Euclid(2:end-t,2:end-t)     = d.Euclid;
EuclidOrig(2:end-t,2:end-t) = d.EuclidOrig;
adaptEuclid     = diag(Euclid(pairs(:,1),pairs(:,2)));
adaptEuclidOrig = diag(EuclidOrig(pairs(:,1),pairs(:,2)));
adaptEuclid(new==1)=-9;
adaptEuclidOrig(new==1)=-9;
adaptEuclid(rept==1)=-5;
adaptEuclidOrig(rept==1)=-5;
 
d.Mink = mncnsim(d.Mink);
Mink(2:end-t,2:end-t) = d.Mink;
adaptMink = diag(Mink(pairs(:,1),pairs(:,2)));
 
angles = rounddig(d.Angle90,2);                   % interpoint angle matrix
anglemat = anglemat-1;                              % will be -1 for blanks
 
anglemat(2:end-t,2:end-t) = angles;                % insert the matrix
anglelist = unique(angles);
ddists = zeros(length(v),length(anglelist));   % one cov per unique angle
adaptangle = diag(anglemat(pairs(:,1),pairs(:,2)));  % angle covariate
i=0;
for angle=anglelist'
    i=i+1;
    ddists(:,i) = adaptEuclid;                 % it's the Euclidean distance
    ddists(adaptangle ~= angle,i) = 0;         % but clear it out unless we match the angle.
end
 
%% prepare output - cut off the first set and mean-center w/o zeros
 
covs.adapteuorig       =                  adaptEuclidOrig(cutofffirst+1:end) ;
covs.adapteu           = mncnleavezeros(  adaptEuclid(cutofffirst+1:end));
covs.adaptMink         = mncnleavezeros(    adaptMink(cutofffirst+1:end));
covs.main              = mncnleavezeros(         main(cutofffirst+1:end));
covs.new               = mncnleavezeros(          new(cutofffirst+1:end));
covs.rept              = mncnleavezeros(         rept(cutofffirst+1:end));
covs.directP           = mncnleavezeros(      directP(cutofffirst+1:end));
covs.directQ           = mncnleavezeros(      directQ(cutofffirst+1:end));
covs.adaptorig         =                        adaptOrig(cutofffirst+1:end) ;
covs.adapt             = mncnleavezeros(        adapt(cutofffirst+1:end));
covs.adaptP            = mncnleavezeros(       adaptP(cutofffirst+1:end));
covs.adaptQ            = mncnleavezeros(       adaptQ(cutofffirst+1:end));
covs.adaptN            = mncnleavezeros(       adaptN(cutofffirst+1:end));
covs.adaptNP           = mncnleavezeros(      adaptNP(cutofffirst+1:end));
covs.adaptNQ           = mncnleavezeros(      adaptNQ(cutofffirst+1:end));
if usetarget
    covs.target        = mncnleavezeros(       target(cutofffirst+1:end));
    covs.aftertarget   = mncnleavezeros(  aftertarget(cutofffirst+1:end));
end
covs.ddists            =                      ddists(cutofffirst+1:end,:);
covs.adaptOrthoEuclid  = mncnleavezeros(adaptOrthoEuclid(cutofffirst+1:end));
covs.sequence          =                            v(cutofffirst+1:end) ;
covs.apd               =                          apd(cutofffirst+1:end) ;
covs.saulapd           =                      saulapd(cutofffirst+1:end) ;
covs.anglelist         = anglelist;

function d = GenTestMat(varargin)
% D = GENTESTMAT(POINTSLIST,DEGREES,MINKEXPONENT)
% degrees is the assumed rotation of axes
% size(pointslist,2) must = 2
% pointslist defaults to a 4x4 square linear grid
% degrees defaults to 0
% minkexponent defaults to 2, specifies the mink power for MinkVec
%
% 2008-08 added minkexponent
% 2007 ddrucker@psych.upenn.edu
 
 
defaultValues={points(0,4),0,2};
nonemptyIdx = ~cellfun('isempty',varargin);
defaultValues(nonemptyIdx) = varargin(nonemptyIdx);
[d.pointslist degrees minkexponent] = deal(defaultValues{:});
 
%% generate matrices
p=mncn(d.pointslist(:,1));
q=mncn(d.pointslist(:,2));
d.OrderG=nchoosek(length(p),2);
 
ext=length(p);
 
for i=1:ext
    for j=1:ext
        % what is the x,y coordinates of this pair?
        pqnosym = rotm([p(i) q(i)] - [p(j) q(j)],degrees);
        pqdiff = abs(pqnosym);
        d.CityP(i,j)  = pqdiff(1);
        d.CityQ(i,j)  = pqdiff(2);
        d.CityNP(i,j) = pqnosym(1);
        d.CityNQ(i,j) = pqnosym(2);
    end
end
 
d.City = d.CityP + d.CityQ;
d.CityN = d.CityNP + d.CityNQ;
 
d.Euclid = vec2sim(getalldistances(d.pointslist,degrees,2)); % degrees doesn't really matter since it's Euclidean
d.Mink = vec2sim(getalldistances(d.pointslist,degrees,minkexponent));
d.Angle  = vec2sim(getallangles(d.pointslist,degrees,45));
d.Angle90 = vec2sim(getallangles(d.pointslist,degrees,90));
d.AngleVec         = sim2vec(d.Angle);
d.Angle90Vec       = sim2vec(d.Angle90);
 
d.CityVec       = mncn(   sim2vec(d.City)        );
d.CityPVec      = mncn(   sim2vec(d.CityP)       );
d.CityQVec      = mncn(   sim2vec(d.CityQ)       );
d.EuclidVec     = mncn(   sim2vec(d.Euclid)      );
d.MinkVec       = mncn(   sim2vec(d.Mink)        );
 
%% orthogonalize EuclidVec wrt CityVec to make OrthoEuclidVec
G=zeros(d.OrderG,2);
G(:,1)=1; % intercept
G(:,2)=d.CityVec;
d.OrthoEuclidVec=orthogonalize(G,d.EuclidVec);

function meancentered = mncnleavezeros(vec)
% MEANCENTERED = MNCNLEAVEZEROS(VEC)
% mean-centers a vector (or matrix), ignoring 0 values (0s stay 0 and are not included
% as part of the calculation of the mean)
%
% 2007 ddrucker@psych.upenn.edu
 
 
nonzeros = find(vec~=0);
nz = numel(nonzeros);
 
sig = sum(vec(nonzeros));
mean = sig/nz;
 
meancentered = vec;
meancentered(nonzeros) = meancentered(nonzeros)-mean;