writeccocov
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
genccocov
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;
gentestmat
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);
mncnleavezeros
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;