SwitchingFCalongAtask/hsmm/gammaParamsCondshD.m

function [params events temp] = gammaParamsCondshD(eventprobs,lens,shape,conds,map,maxDur)
% @ lens: mean trial duration per condition [N Conditions]
offset1=3;
spacing=5;
nconds=max(conds);
map=cat(2,map,ones(nconds,1)); % add one column maps[conditions, N stages]
nstates=size(eventprobs,3);  % number of bumps
temp=zeros(maxDur,nstates,nconds);
params=zeros(nstates+1,2,nconds);
for i = 1:nconds
    temp(:,:,i)=reshape(mean(eventprobs(:,conds==i,:),2),maxDur,nstates);
    % temp[samples,bumps,condition] mean accros trials per sample, like an ERP
end
maxobs=(maxDur+1)-(2*(offset1-1)+(nstates-1)*spacing); % max duration of trial without bumps' length and offset
events=zeros(nstates+2,nconds);
for j = 1:nconds
    offset=-2;
    for i = 1:nstates
        if map(j,i)~=0 % if condition/current stage estimation requested
            offset=offset+spacing;
            temp(1:maxobs,i,j)=temp(offset:offset+maxobs-1,i,j);
            % global expected location of a bump
            events(i+1,j)=[1:maxobs]*temp(1:maxobs,i,j);
        else % if not requested: assign to next stage 'events' of current stage/state
            events(i+1,j)=events(i,j);
        end
    end
    events(i+2,j)=lens(j)-offset-1; % mean trial duration per condition - offest - 1
end
temp(maxobs+1:end,:,:)=[];
for j = 1:nstates+1 % iterate stages
    % get the index of requested stages to compute the duration of the flats.
    uniques=unique(map(:,j)); % take uinque requested map/condition values in bump
    uniques=uniques(uniques>0); % discard zeros (not requested)
    for i = 1:length(uniques) % iteratee unique values
        subs=find(map(:,j)==i); % index of unque value in a stage
        % Gamma scale : duration of a flat (event(j+1) - event(j)) divided by shape
        params(j,:,subs)=repmat([shape,mean(events(j+1,subs)-events(j,subs))/shape],[1,1,length(subs)]);
    end
end
params(1,2,:)=params(1,2,:)-.5/shape;
params(2:nstates,2,:)=params(2:nstates,2,:)+.5/shape;