=====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;