We Helped With This MATLAB Programming Homework: Have A Similar One?

clear all
close all

dat=[];

cd results

obsList=[1:10,12:22];

for f=1:length(obsList);

observer=num2str(obsList(f));

dir=['Obs' observer];
cd (dir)

fileName=['output' observer '.txt'];

res=load(fileName);
dat(end+1,1:3)=res(1:3);
cd ..

end;

cd ..

results=[transpose(obsList),dat];
dlmwrite('globalMotionResults.txt', results, '	');

close all; clear all
observer=eval(input('Type the observer number: ','s'));

%--------------------------------------------------------------------------
%                       Stimulus setup
%--------------------------------------------------------------------------

%this will want to update later
%signalDots=200;
%noiseDots=50;
    nframes     = 120; % number of animation frames in loop
    mon_width   = 40.5;   % horizontal dimension of viewable screen (cm)
    v_dist      = 50;   % viewing distance (cm)
    dot_speed   = 1;  % dot speed (deg/sec)
    ndots       = 100; % total number of dots
    max_d       = 10;   % maximum radius of  annulus (degrees)
    min_d       = 0.2;    % minumum (deg)
    dot_w       = 0.1;  % width of dot (deg)
    fix_r       = 0.15; % radius of fixation point (deg)

%--------------------------------------------------------------------------
%                       Staircase input
%--------------------------------------------------------------------------

% stair input
probeset    = -15:0.5:15;       % set of possible probe values
meanset     = -10:0.2:10;       % sampling of pses, doesn't have to be the same as probe set
slopeset    = [.5:.1:5].^2;     % set of slopes, quad scale
lapse       = 0.05;             % lapse/mistake rate
guess       = 0.50;             % guess rate / minimum correct response rate (for detection expt: 1/num_alternative, 0 for discrimination expt)

% probeset    = 0:round(ndots./50):ndots; 
% meanset     = 0:round(ndots./50):round(ndots./2); 
% slopeset    = [.5:.1:5].^2;
% lapse       = 0.05;  
% guess       = 0.5; 
% general settings
ntrial  = 15;
qpause  = true;    % pause after every iteration? (press any key to continue)
qplot   = false;    % plot information about each trial? (this pauses as well, regardless of whether you specified qpause as true)

% model observer parameters 
%qusemodel = true;   % use model observer to get responses? Or, if false, input responses by hand (0/1)
%truepse = 0;        % inflection point (50% if guess rate is 0)
%truedl  = 4;        % (75%-25% point)/2 if guess rate is 0. In general, take the position of the halfway points between the inflection point and the upper and lower asymptotes, then its the distance between them

% model observer definition. uses a cum normal for psychometric function,
% the formula for which is equivalent to what is used by the staircase
% internally (if it is set up to use a cumulative Gaussian)
if guess==0
    resp    = @(probe)  lapse/2 + (1-lapse)      *normcdf(probe,truepse,truedl/sqrt(2)/erfinv(0.5)) > rand;
else
    resp    = @(probe)   guess  + (1-lapse-guess)*normcdf(probe,truepse,truedl/sqrt(2)/erfinv(0.5)) > rand;
end
% Create staircase instance. If you want to interleave staircases, or
% otherwise run multiple conditions, create one staircase per condition.
% You can store these in a cell-array and of course use different settings
% for each as needed.
stair = MinExpEntStair('v2');
% use lookup table to cache pvalues and avoid unnecessary evaluations of
% psychometric function? Can require lots of memory, especially when
% stepsize of probeset and meanset is not equal. Call before calling
% stair.init.
%stair.set_use_lookup_table(true);
% option: use logistic instead of default cumulative normal. best to call
% before stair.init
% stair('set_psychometric_func','logistic');
% init stair
stair.init(probeset,meanset,slopeset,lapse,guess);
% option: use a subset of all data for choosing the next probe, use
% proportion of available data (good idea for robustness - see docs)
stair.toggle_use_resp_subset_prop(10,.9);
% option: instead of choosing a random value for the first probe,
% you can set which value is to be tested first.
stair.set_first_value(3);

%--------------------------------------------------------------------------
%                   Screen initialisation
%--------------------------------------------------------------------------

% Make sure that the computer is running the OpenGL psych toolbox
AssertOpenGL;
Screen('Preference', 'SkipSyncTests', 1 );

% Find the screen to use for display
screenid = max(Screen('Screens'));

% Set the black and white index
black = BlackIndex(screenid);
white = WhiteIndex(screenid);

% Color codes
red = [255 0 0];
blue = [0 0 255];
grey=[128 128 128];

grey=127.5;

	%calibration stuff 
	a=transpose([0:5.1:255]);
	lum=load('cali1.txt'); %use 1 if using monitor 1
	p=polyfit(lum,a,4);
	maxLum=max(lum);
	minLum=min(lum);
	lumrange=(maxLum-minLum);
	midgrey=lumrange./2+minLum;
	midgrey=p(5)+(p(4).*midgrey)+(p(3).*midgrey.^2)+(p(2).*midgrey.^3)+(p(1).*midgrey.^4);
	%black=p(5)+(p(4).*black)+(p(3).*black.^2)+(p(2).*black.^3)+(p(1).*black.^4);
	%white=p(5)+(p(4).*white)+(p(3).*white.^2)+(p(2).*white.^3)+(p(1).*white.^4);

% Open a double buffered fullscreen window and draw a gray background
% to front and back buffers enable 32 floating point
[window winRect] = Screen('OpenWindow', screenid, black);%, [], 32, 2, [], []);

% Hide the cursor
HideCursor;

% Enable alpha blending for antialiasing
Screen('BlendFunction', window, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

% Get the width and height of the window in pixels
[screenXpix, screenYpix] = Screen('WindowSize', window);
[center(1), center(2)] = RectCenter(winRect);

 
% Viewing distance in cm    
distCm = 50 ;
% Screen information (cm)

screenXcm = 40.5;
screenYcm = 30.0;

%atand(screenXcm./distCm);
%=39.01 degrees
%1024 pixels = 39.01 deg
%26.25 pixPerDeg

	% Calculate the pixels per centimetre
	pixPerCm = screenXpix / screenXcm;

% Measure the vertical refresh rate of the monitor
ifi = Screen('GetFlipInterval', window);

% Get the maximum priority level
topPriorityLevel = MaxPriority(window);
%--------------------------------------------------------------------------
%                       Keyboard information
%--------------------------------------------------------------------------

% Unify the keyboard names for mac and windows computers
 KbName('UnifyKeyNames');
spaceKey = KbName('space');
rightKey = KbName('rightArrow');
leftKey = KbName('leftArrow');
escapeKey=KbName('ESCAPE');
% No keys are down
keyIsDown = 0;
   fps=Screen('FrameRate',window);

    ppd = pi*(winRect(3)-winRect(1))/atan(mon_width/v_dist/2)/360;% pix per deg
    pfs = dot_speed * ppd / fps;                       % dot speed (pixels/frame)
    s = dot_w * ppd;                                   % dot size (pixels)
    fix_cord = [center-fix_r*ppd center+fix_r*ppd];
    rmax = max_d * ppd;	% maximum radius of annulus (pixels from center)
    rmin = min_d * ppd; % minimum
    rad = rmax * sqrt(rand(ndots,1));	               % r
    rad(rad<rmin) = rmin;
    t = 2*pi*rand(ndots,1);                            % theta polar coordinate
    cs = [cos(t), sin(t)];
    xy = [rad rad] .* cs;   % dot positions in Cartesian coordinates (pix from center)

    mdir = 2 * floor(rand(ndots,1)+0.5) - 1;    % motion direction 
    dr = pfs * mdir;                            % change in r / frame (pixels)
    dxdy = [dr dr] .* cs;                       % change in x and y / frame(pixels)
		xy=transpose(xy); %the original column vector		
%--------------------------------------------------------------------------
%                       Start trials
%--------------------------------------------------------------------------

allResponsesStair1=[];

for ktrial = 1:ntrial
    if nframes<12; nframes=12; end;
    
    % trial
    [p,entexp,ind]  = stair.get_next_probe();      % get next probe to test
    %fprintf('%d, new sample point: %f
expect ent: %f
', ...
    %    ktrial,p,entexp(ind));
    
    noiseDots=round(ndots-(p));%noise chosen by staircase
    signalDots=round(p);%signal dots are a ratio
    		%chose initial signal dots
		xymatrix=xy(:,1:signalDots); %the first vector of signal dots

	        %chose noise dots
		noiseindex=round(rand(1,noiseDots).*(length(xy)-1));
		noise=xy(:,noiseindex+1);
        
        
        %--------------------------------------------------------------------------
%                       Stimulus presentation
%--------------------------------------------------------------------------

try

direction=round(rand(1));

flipCount=0;
dotLifetime=8; %number of frames a dot ilves for
dotLivesNoise=round(rand(1,length(noise)).*dotLifetime); %each dot is assigned lifetime
dotLivesSignal=round(rand(1,length(xymatrix)).*dotLifetime); %each dot is assigned lifetime

        Screen('FillRect', window, midgrey);
        Screen('Flip', window);

time2=GetSecs;

    for numFrames=1:nframes
          Screen('DrawDots',window,xymatrix,5,black,center,1);
	  Screen('DrawDots',window,noise,5,black,center,1);

             Screen('Flip', window);
		flipCount=flipCount+1;
		dotLivesNoise=dotLivesNoise+1;
		dotLivesSignal=dotLivesSignal+1;

	        %update the matrix
		%signal moves
		    if direction == 1; %direction 1 = right;
                xymatrix=[xymatrix(1,:)+pfs;xymatrix(2,:)];
		    elseif direction == 0; %0 is left
                xymatrix=[xymatrix(1,:)-pfs;xymatrix(2,:)];
		    end;
		%noise jitters
		noise=noise+randn(2,noiseDots).*pfs;
		
		%update a certain number of signal dots
		for f=1:length(dotLivesSignal);
		if dotLivesSignal(f)>dotLifetime;
		t1=2*pi*rand(1,1);
		xymatrix(:,f)=[cos(t1);sin(t1)].*rmax.*sqrt(rand(1));
		dotLivesSignal(f)=0;
		end;
		end;

		%update a certain number of noise dots
		for f=1:length(dotLivesNoise);
		if dotLivesNoise(f)>dotLifetime;
		t1=2*pi*rand(1,1);
		noise(:,f)=[cos(t1);sin(t1)].*rmax.*sqrt(rand(1));
		dotLivesNoise(f)=0;
		end;
		end;

		%kill dots that go out
		if max(sqrt(xymatrix(1,:).^2+xymatrix(2,:).^2))>rmax; g=1; end;
		for f=1:length(xymatrix);
		if sqrt(xymatrix(1,f).^2+xymatrix(2,f).^2)>rmax; 
		t1=2*pi*rand(1,1);
		xymatrix(:,f)=[cos(t1);sin(t1)].*rmax.*sqrt(rand(1));
        end;
		end;

		if max(sqrt(noise(1,:).^2+noise(2,:).^2))>rmax; g=1; end;
		for f=1:length(noise);
		if sqrt(noise(1,f).^2+noise(2,f).^2)>rmax; 
		%noise(:,f)=round(randn(2,1).*125); end;
		t1=2*pi*rand(1,1);
		noise(:,f)=[cos(t1);sin(t1)].*rmax.*sqrt(rand(1));
        end;
		end;
    end %end presentation for loop
	%put up a blank screen;
        Screen('FillRect', window, midgrey);
        Screen('Flip', window);

time1=GetSecs;

%response code
respToBeMade=true;
while respToBeMade;
        [keyIsDown,secs, keyCode] = KbCheck;
        if keyCode(escapeKey)
            respToBeMade=false;
		Screen('CloseAll');
          elseif keyCode(rightKey)
            respToBeMade=false;
            if direction == 1;%going right
            response=1;
            WaitSecs(0.1);
            else; response=0;
            WaitSecs(0.1);
            end;
	      elseif keyCode(leftKey);
            respToBeMade=false;
            if direction == 0;%going left
            response=1;
            WaitSecs(0.1);
            else; response=0;
            WaitSecs(0.1);
            end;
        end;
end; %end while loop

catch
Screen('CloseAll');
end;        %end try catch loop
        
    
   % if qusemodel % set whether model creates response or you do by manual input
   %     % get observer response from model observer
   %     r = resp(p);
   %     fprintf('response: %d
',r);
   % else
   %     % make the response yourself, provide either 0 or 1 (actually,
   %     % anything below and including 0 or anything above 0)
   %     r = input(sprintf('r(%d): ',ktrial));
   %     qpause = false;
   % end
r=response; qpause=false;
    stair.process_resp(r);                        % store response in staircase
    % end trial
    
    if ktrial == ntrial || qplot
        
        [m,s,loglik]    = stair.get_fit();
        [ps,rs]         = stair.get_history();
        
        figure(1);
        subplot(1,3,1);
        imagesc(exp(.5*loglik));
        set(gca,'YTick',1:4:length(slopeset));
        set(gca,'YTickLabel',slopeset(1:4:end));
        set(gca,'XTick',1:5:length(meanset));
        set(gca,'XTickLabel',meanset(1:5:end));
        title('estimated likelihood function');
        xlabel('mean (a)')
        ylabel('slope (b)')
        
        subplot(1,3,2);
        hold off;
        if ~isempty(entexp)
            plot(probeset,entexp,'k-o');
            hold on;
            plot(ps(ktrial)*[1,1],[min(entexp),max(entexp)],'r-');
        else
            plot(ps(ktrial)*[1,1],[0,1],'r-');
        end
        title('expected entropy vs probe');
        xlabel('possible probe values')
        xlim([min(probeset) max(probeset)]);
        
        subplot(1,3,3);
        pind = find(rs>0);
        nind = setdiff(1:length(ps),pind);
        plot(1:length(ps),ps,'k-',pind,ps(pind),'bo',nind,ps(nind),'ro');
        ylim([min(probeset) max(probeset)]);
        title('probe-resp history');
     end
    
    % pause if needed
    if (ktrial ~= ntrial) && (qpause || qplot)
        pause;
    end
    
end % loop over trials
Screen('CloseAll');

%%% show final results
%  [mu,sigma] = stair('get_fit');    % get fitted parameters of cumulative Gaussian
%  DL = sigma*erfinv(.5)*sqrt(2)     % convert sigma (std) to DL (75%-25% point)/2
% get DL from staircase directly, NB: the space of the outputted
% loglikelihood is the mean/slope space as defined atop this script, its
% not a PSE/DL space
[PSEfinal,DLfinal,loglikfinal]  = stair.get_PSE_DL();
finalent                        = sum(-exp(loglikfinal(:)).*loglikfinal(:));
fprintf('final estimates:
PSE: %f
DL: %f
ent: %f
',PSEfinal,DLfinal,finalent);
% for actual offline fitting of your data, you would probably want to use a
% dedicated toolbox such as Prins, N & Kingdom, F. A. A. (2009) Palamedes:
% Matlab routines for analyzing psychophysical data.
% http://www.palamedestoolbox.org.
% Also note that while the staircase runs far more rebust when a small
% lapse rate is assumed, it is common to either fit the psychometric
% function without a lapse rate, or otherwise with the lapse rate as a free
% parameter (possibily varying only over subjects, but not over conditions
% within each subject).

figure(2);
imagesc(exp(.5*loglikfinal));
set(gca,'YTick',1:4:length(slopeset));
set(gca,'YTickLabel',slopeset(1:4:end));
set(gca,'XTick',1:5:length(meanset));
set(gca,'XTickLabel',meanset(1:5:end));
xlabel('$mu$','interpreter','latex')
switch stair.get_psychometric_func()
    case 'cumGauss'
        title('estimated likelihood function - cumulative Gaussian')
        ylabel('$sigma$','interpreter','latex')
    case 'logistic'
        title('estimated likelihood function - logistic')
        ylabel('$s$','interpreter','latex')
end