function [fitps, best_fitps, best_pM] = grtfit(oM,inits,fp,nfits,js_fact,dbtype,mod_id);

% function [fitps, best_fitps, best_pM] = grtfit(oM,inits,fp,nfits,js_fact,dbtype,mod_id);
%
% fitps     = structure containing fitted parameters and fit statistic
% oM        = confusion matrices (4x4x5 array)
% inits     = structure containing initial parameter values
% fp        = structure of logicals indicating which parameters are free
% nfits     = number of accepted fits wanted
% js_fact   = scales size of jumps in parameter space
% dbtype    = 'linear' or 'piecewise' or 'maxpost'
% mod_id    = name of model being fit (e.g., ns_m, ns_PI, etc...)
%
% see file for more details

% fields in inits must be identical to those in fitps, namely
%
% for dbtype = 'linear' or 'piecewise'
%
% fitps = struct('maa',{},'mab',{},'mba',{},'mbb',{},'caa',{},'cab',{},'cba',{},'cbb',{},...
%     'ax',{},'bx',{},'cx1',{},'dx1',{},'cx3',{},'dx3',{},'cx5',{},'dx5',{},'sx',{},...
%     'ay',{},'by',{},'cy2',{},'dy2',{},'cy3',{},'dy3',{},'cy4',{},'dy4',{},'sy',{},...
%     'LL',{},'g2',{},'bic',{});
%
% for dbtype = 'maxpost'
%
% fitps = struct('maa',{},'mab',{},'mba',{},'mbb',{},'caa',{},'cab',{},'cba',{},'cbb',{},...
%     'c1aa',{},'c1ab',{},'c1ba',{},'c1bb',{},'c2aa',{},'c2ab',{},'c2ba',{},'c2bb',{},'c3aa',{},...
%     'c3ab',{},'c3ba',{},'c3bb',{},'c4aa',{},'c4ab',{},'c4ba',{},'c4bb',{},'c5aa',{},'c5ab',{},...
%     'c5ba',{},'c5bb',{},'g2',{},'LL',{},'bic',{});%
%
% fp fields are as follows
%
% fp = struct('ax',{},'bx',{},'cx',{},'cx_k',{},'dx',{},'dx_k',{},'sx',{},...
%                 'ay',{},'by',{},'cy',{},'cy_k',{},'dy',{},'dy_k',{},'sy',{},...
%                 'vxa',{},'vxb',{},'vya',{},'vyb',{},'raa',{},'rab',{},'rba',{},'rbb',{},...
%                 'mxa',{},'mxb',{},'mya',{},'myb',{},'mab',{},'mba',{},'mbb',{},...
%                 'vaa',{},'vab',{},'vba',{},'vbb',{});
%
% 1 = it can vary; 0 = it can't:
% bounds:   'ax',{},'bx',{},'cx',{},'dx',{},'sx',{},'ay',{},'by',{},'cy',{},'dy',{},'sy',{},...
% (co)var:  'raa',{},'rab',{},'rba',{},'rbb',{},'vaa',{},'vab',{},'vba',{},'vbb',{});
% means:    'mab',{},'mba',{},'mbb',{},...
%
% 1 = fixed across levels; 0 = not fixed across levels:
% variance: 'vxa',{},'vxb',{},'vya',{},'vyb',{},...
% means:    'mxa',{},'mxb',{},'mya',{},'myb',{},...
% bounds: 'cx_k',{},'dx_k',{}'cy_k',{},'dy_k',{}
%
% if strcmp(dbtype,'linear'), the c and d parameters are slope and
% intercept, respectively, when free
% if strcmp(dbtype,'piecewise'), the c and d parameters are intercepts,
% when free.
% if strcmp(dbtype,'maxpost'), the a, b, c, and d parameters are baserate
% 'scaling' parameters, when free.
%

timeseed = 1;

if timeseed
    rand('state',sum(100*clock));
else
    rand('state',5); % for debugging, checking variations in zpro (see below)
end
        
if strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')
    fitps = struct('maa',{},'mab',{},'mba',{},'mbb',{},'caa',{},'cab',{},'cba',{},'cbb',{},...
        'ax',{},'bx',{},'cx1',{},'dx1',{},'cx3',{},'dx3',{},'cx5',{},'dx5',{},'sx',{},...
        'ay',{},'by',{},'cy2',{},'dy2',{},'cy3',{},'dy3',{},'cy4',{},'dy4',{},'sy',{},'g2',{},'LL',{},'bic',{});
elseif strcmp(dbtype,'maxpost')
    fitps = struct('maa',{},'mab',{},'mba',{},'mbb',{},'caa',{},'cab',{},'cba',{},'cbb',{},...
        'c1aa',{},'c1ab',{},'c1ba',{},'c1bb',{},'c2aa',{},'c2ab',{},'c2ba',{},'c2bb',{},'c3aa',{},...
        'c3ab',{},'c3ba',{},'c3bb',{},'c4aa',{},'c4ab',{},'c4ba',{},'c4bb',{},'c5aa',{},'c5ab',{},...
        'c5ba',{},'c5bb',{},'g2',{},'LL',{},'bic',{});
end

% initial parameter settings

if strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')
    npar = 34; % 5 m, 11 c, 18 db
elseif strcmp(dbtype,'maxpost')
    npar = 31; % 5 m, 11 c, 15 'db'
end

maa = [0; 0]; % not free
m = [maa inits.mab inits.mba inits.mbb];
if fp.mxa
    m(1,2) = m(1,1);
    npar = npar - 1;
end
if fp.mxb
    m(1,4) = m(1,3);
    npar = npar - 1;
end
if fp.mya
    m(2,3) = m(2,1);
    npar = npar - 1;
end
if fp.myb
    m(2,4) = m(2,2);
    npar = npar - 1;
end
c = zeros(2,2,4);
c(:,:,1) = inits.caa;
c(1,1,1) = 1; % not free
c(:,:,2) = inits.cab;
c(:,:,3) = inits.cba;
c(:,:,4) = inits.cbb;
if fp.vxa
    c(1,1,2) = c(1,1,1);
    npar = npar - 1;
end
if fp.vxb
    c(1,1,4) = c(1,1,3);
    npar = npar - 1;
end
if fp.vya
    c(2,2,3) = c(2,2,1);
    npar = npar - 1;
end
if fp.vyb
    c(2,2,4) = c(2,2,2);
    npar = npar - 1;
end
if ~fp.vaa
    npar = npar - 1;
end
if ~fp.vab
    npar = npar - 2;
end
if ~fp.vba
    npar = npar - 2;
end
if ~fp.vbb
    npar = npar - 2;
end
if ~fp.mab
    npar = npar - 2;
end
if ~fp.mba
    npar = npar - 2;
end
if ~fp.mbb
    npar = npar - 1;
end
if ~fp.raa
    c(1,2,1) = 0;
    c(2,1,1) = 0;
    npar = npar - 1;
end
if ~fp.rab
    c(1,2,2) = 0;
    c(2,1,2) = 0;
    npar = npar - 1;
end
if ~fp.rba
    c(1,2,3) = 0;
    c(2,1,3) = 0;
    npar = npar - 1;
end
if ~fp.rbb
    c(1,2,4) = 0;
    c(2,1,4) = 0;
    npar = npar - 1;
end
if strcmp(dbtype,'maxpost')
    dbp(1,:) = [inits.c1aa inits.c1ab inits.c1ba inits.c1bb];
    dbp(2,:) = [inits.c2aa inits.c2ab inits.c2ba inits.c2bb];
    dbp(3,:) = [inits.c3aa inits.c3ab inits.c3ba inits.c3bb];
    dbp(4,:) = [inits.c4aa inits.c4ab inits.c4ba inits.c4bb];
    dbp(5,:) = [inits.c5aa inits.c5ab inits.c5ba inits.c5bb];    
else
    dbp(1,:) = [inits.ax inits.bx inits.cx1 inits.dx1 inits.sx];
    dbp(2,:) = [inits.ay inits.by inits.cy2 inits.dy2 inits.sy];
    dbp(3,:) = [inits.ax inits.bx inits.cx3 inits.dx3 inits.sx];
    dbp(4,:) = [inits.ay inits.by inits.cy4 inits.dy4 inits.sy];
    dbp(5,:) = [inits.ax inits.bx inits.cx5 inits.dx5 inits.sx];
    dbp(6,:) = [inits.ay inits.by inits.cy3 inits.dy3 inits.sy];
end
[nstim,nresp,nbr] = size(oM);

if strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')
    if ~fp.ax
        dbp(1,1) = 0;
        dbp(3,1) = 0;
        dbp(5,1) = 0;
        npar = npar - 1;
    end
    if ~fp.bx
        dbp(1,2) = 0;
        dbp(3,2) = 0;
        dbp(5,2) = 0;
        npar = npar - 1;
    end
    if ~fp.ay
        dbp(2,1) = 0;
        dbp(6,1) = 0;
        dbp(4,1) = 0;
        npar = npar - 1;
    end
    if ~fp.by
        dbp(2,2) = 0;
        dbp(6,2) = 0;
        dbp(4,2) = 0;
        npar = npar - 1;
    end
    if ~fp.sx
        dbp(1,5) = 0;
        dbp(3,5) = 0;
        dbp(5,5) = 0;
        npar = npar - 1;
    end
    if ~fp.sy
        dbp(2,5) = 0;
        dbp(6,5) = 0;
        dbp(4,5) = 0;
        npar = npar - 1;
    end
    if ~fp.cx
        dbp(1,3) = 0;
        dbp(3,3) = 0;
        dbp(5,3) = 0;
        npar = npar - 3;
    end
    if ~fp.cy
        dbp(2,3) = 0;
        dbp(6,3) = 0;
        dbp(4,3) = 0;
        npar = npar - 3;
    end
    if fp.cx_k
        dbp(1,3) = dbp(3,3);
        dbp(5,3) = dbp(3,3);
        if fp.cx
            npar = npar - 2;
        end
    end
    if fp.cy_k
        dbp(2,3) = dbp(6,3);
        dbp(4,3) = dbp(6,3);
        if fp.cy
            npar = npar - 2;
        end
    end
    if fp.dx_k
        dbp(1,4) = dbp(3,4);
        dbp(5,4) = dbp(3,4);
        if fp.dx
            npar = npar - 2;
        end
    end
    if fp.dy_k
        dbp(2,4) = dbp(6,4);
        dbp(4,4) = dbp(6,4);
        if fp.dy
            npar = npar - 2;
        end
    end
    if ~fp.dx
        dbp(1,4) = mean(m(1,:));
        dbp(3,4) = dbp(1,4);
        dbp(5,4) = dbp(1,4);
        npar = npar - 3;
    end
    if ~fp.dy
        dbp(2,4) = mean(m(2,:));
        dbp(6,4) = dbp(2,4);
        dbp(4,4) = dbp(2,4);
        npar = npar - 3;
    end
elseif strcmp(dbtype,'maxpost')
    mp_free = 0;
    mp_feat_g = 0;
    mp_feat_r = 0;
    if ~fp.ay & ~fp.by & ~fp.ax & ~fp.bx
        if fp.cy & fp.dy & fp.cx & fp.dx
            if ~fp.cy_k & ~fp.cx_k
                % feat = feature; g = general; 'upper middle' model
                % other dimension's 'neutral' parameters can vary from .5
                mp_feat_g = 1; 
                npar = npar - 9;
            elseif fp.cy_k & fp.cx_k
                % feat = feature; r = restricted; 'lower middle' model
                % other dimensions 'neutral' parameters are .5
                mp_feat_r = 1; 
                npar = npar - 11;
            end
            feat_params = zeros(2,6);
        else
            mp_br = 1; % br = base rate; restricted model
            npar = npar - 15;
        end
    else
        mp_free = 1; % each distribution gets its own weight
    end
end


% implementing BIC calculation; added 5/30/07
big_N = sum(sum(sum(oM)));
complexity = npar*log(big_N);

% bin widths
bin = .05;

cf = -1e15; % zero likelihood, used to avoid unnecessary calculations with obviously bad parameter settings

% initialization of fit statistics
g2_new = -cf;
g2_old = inits.g2;
g2_min = g2_old;
bic_new = -cf;
bic_old = inits.bic;
bic_min = bic_old;
LL_new = cf;
LL_old = inits.LL;
LL_max = LL_old;
best_fitps = inits;
%best_pM = predcon(best_fitps,dbtype,0,fp,oM);

% initialize # of fits
n = 0;

% initialize assorted other things
pv = zeros(4,1);
n_rejected = 0;

zerofix = '1in4';

if strcmp(zerofix,'loglinear')
    zpro = .5;
elseif strcmp(zerofix,'1in4')
    zpro = .25;
elseif strcmp(zerofix,'1in100')
    zpro = .01;
elseif strcmp(zerofix,'1/N')
    zpro = 1/sum(sum(sum(oM)));
else
    error('A zero-cell correction must be specified.')
end

oM = oM + zpro;

Nr_oM = zeros(nstim,nbr);
for ci = 1:nbr
    for i = 1:nstim
        Nr_oM(i,ci) = sum(oM(i,:,ci));
    end
end

% initialize dbtype maxpost baserate scaling parameters
if strcmp(dbtype,'maxpost')
    if sum(dbp(3,1:4))>1
    
        c1_N = sum(Nr_oM(:,1));
        c1_w = Nr_oM(:,1)/c1_N;
        c2_N = sum(Nr_oM(:,2));
        c2_w = Nr_oM(:,2)/c2_N;
        c3_N = sum(Nr_oM(:,3));
        c3_w = Nr_oM(:,3)/c3_N;
        c4_N = sum(Nr_oM(:,4));
        c4_w = Nr_oM(:,4)/c4_N;
        c5_N = sum(Nr_oM(:,5));
        c5_w = Nr_oM(:,5)/c5_N;
        
        dbp(1,1) = c1_w(1);
        dbp(1,2) = c1_w(2);
        dbp(1,3) = c1_w(3);
        dbp(1,4) = c1_w(4);
        
        dbp(2,1) = c2_w(1);
        dbp(2,2) = c2_w(2);
        dbp(2,3) = c2_w(3);
        dbp(2,4) = c2_w(4);
        
        dbp(3,1) = c3_w(1);
        dbp(3,2) = c3_w(2);
        dbp(3,3) = c3_w(3);
        dbp(3,4) = c3_w(4);
        
        dbp(4,1) = c4_w(1);
        dbp(4,2) = c4_w(2);
        dbp(4,3) = c4_w(3);
        dbp(4,4) = c4_w(4);
        
        dbp(5,1) = c5_w(1);
        dbp(5,2) = c5_w(2);
        dbp(5,3) = c5_w(3);
        dbp(5,4) = c5_w(4);
        
    end
    % feat_params structure
    % (1,1) = stim.1 condition 1 weight
    % (1,2) = stim1. condition 2 weight
    % (1,3) = stim.1 condition 3 ('neutral') weight
    % (1,4) = stim1. condition 4 weight
    % (1,5) = stim.1 condition 5 weight
    % (1,6) = stim1. condition 3 ('neutral') weight
    % (2,:) = 1-(1,:)
    
    if mp_feat_g || mp_feat_r
        feat_params(1,1) = dbp(1,1)/( dbp(1,1) + dbp(1,2) );
        feat_params(1,2) = dbp(2,1)/( dbp(2,1) + dbp(2,3) );
        feat_params(1,3) = dbp(3,1)/( dbp(3,1) + dbp(3,2) );
        feat_params(1,4) = dbp(4,1)/( dbp(4,1) + dbp(4,3) );
        feat_params(1,5) = dbp(5,1)/( dbp(5,1) + dbp(5,2) );
        feat_params(1,6) = dbp(3,1)/( dbp(3,1) + dbp(3,3) );
    end
    
end

while n <= nfits
    
    disp(mod_id)
    disp(['working on fit ' num2str(n) ' out of ' num2str(nfits)])
        
    %LL_new = [];
    
    % jumping distribution funtime party hour
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % want to make it adaptive in some %
    % sense, e.g., smaller jumps for a %
    % bit before going big when 'stuck'%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    big_jump = 3;
    wee_jump = .25;
    if n_rejected > 9
        if n_rejected < 20
            js_m = wee_jump*.01*js_fact; % scales the jumping distributions for means
            js_v = wee_jump*.01*js_fact; % scales the js for variances
            js_r = wee_jump*.01*js_fact; % scales the js for covariances
            if strcmp(dbtype,'maxpost')
                js_a = wee_jump*.01*js_fact; % scales the js for the cubic terms
                js_b = wee_jump*.01*js_fact; % scales the js for the quadratic terms
                js_c = wee_jump*.01*js_fact; % scales the js for the slope terms
                js_d = wee_jump*.01*js_fact; % scales the js for the intercepts
            else
                js_a = wee_jump*.2*js_fact; % scales the js for the cubic terms
                js_b = wee_jump*.2*js_fact; % scales the js for the quadratic terms
                js_c = wee_jump*.01*js_fact; % scales the js for the slope terms
                js_d = wee_jump*.01*js_fact; % scales the js for the intercepts
            end
            js_s = wee_jump*.15*js_fact; % scales the js for the offset terms
        else
            js_m = big_jump*.02*js_fact; % scales the jumping distributions for means
            js_v = big_jump*.02*js_fact; % scales the js for variances
            js_r = big_jump*.02*js_fact; % scales the js for covariances
            if strcmp(dbtype,'maxpost')
                js_a = big_jump*.01*js_fact; % scales the js for the cubic terms
                js_b = big_jump*.01*js_fact; % scales the js for the quadratic terms
                js_c = big_jump*.01*js_fact; % scales the js for the slope terms
                js_d = big_jump*.01*js_fact; % scales the js for the intercepts
            else
                js_a = big_jump*.2*js_fact; % scales the js for the cubic terms
                js_b = big_jump*.2*js_fact; % scales the js for the quadratic terms
                js_c = big_jump*.02*js_fact; % scales the js for the slope terms
                js_d = big_jump*.02*js_fact; % scales the js for the intercepts
            end
            js_s = big_jump*.15*js_fact; % scales the js for the offset terms
            accept = 1;
        end
    end
    if n_rejected == 0
        js_m = .01*js_fact; % scales the jumping distributions for means
        js_v = .01*js_fact; % scales the js for variances
        js_r = .01*js_fact; % scales the js for covariances (uses betainv)
        if strcmp(dbtype,'maxpost')
            js_a = .01*js_fact; % scales the js for the cubic terms
            js_b = .01*js_fact; % scales the js for the quadratic terms
            js_c = .01*js_fact; % scales the js for the slope terms
            js_d = .01*js_fact; % scales the js for the intercepts
        else
            js_a = .01*js_fact; % scales the js for the cubic terms
            js_b = .01*js_fact; % scales the js for the quadratic terms
            js_c = .01*js_fact; % scales the js for the slope terms
            js_d = .01*js_fact; % scales the js for the intercepts
        end
        js_s = .05*js_fact; % scales the js for the offset terms
        accept = 0; % logical that determines whether or not to accept the current fit
    end
    
    m_old = m;
    c_old = c;
    dbp_old = dbp;
    
    % means

    if fp.mab
        if ~fp.mxa
            m(1,2) = m(1,2) + js_m*randn(1); % m_xab free
        end
    end
    if fp.mba & fp.mbb
        if ~fp.mxb
            m(1,3) = m(1,3) + js_m*randn(1); % new m_xba
            m(1,4) = m(1,4) + js_m*randn(1); % new m_xbb
        else
            m(1,3) = m(1,3) + js_m*randn(1); % new m_xb
            m(1,4) = m(1,3); % only one m_xb
        end
    end
    if fp.mba
        if ~fp.mya
            m(2,3) = m(2,3) + js_m*randn(1); % m_yba free
        end
    end
    if fp.mab & fp.mbb
        if ~fp.myb
            m(2,2) = m(2,2) + js_m*randn(1); % new m_yab
            m(2,4) = m(2,4) + js_m*randn(1); % new m_ybb
        else
            m(2,2) = m(2,2) + js_m*randn(1); % new m_yb
            m(2,4) = m(2,2); % only one m_yb
        end
    end
    
    % covariance

    if fp.raa
        raa = c(1,2,1)/sqrt(c(1,1,1)*c(2,2,1));
    end
    if fp.rab
        rab = c(1,2,2)/sqrt(c(1,1,2)*c(2,2,2));
    end
    if fp.rba
        rba = c(1,2,3)/sqrt(c(1,1,3)*c(2,2,3));
    end
    if fp.rbb
        rbb = c(1,2,4)/sqrt(c(1,1,4)*c(2,2,4));
    end
    
    if fp.vab
        if ~fp.vxa
            c(1,1,2) = abs(c(1,1,2) + js_v*randn(1)); % new v_x for ab
        else
            c(1,1,2) = c(1,1,1);
        end
    end
    if fp.vba & fp.vbb
        if ~fp.vxb
            c(1,1,3) = abs(c(1,1,3) + js_v*randn(1)); % new v_x for ba
            c(1,1,4) = abs(c(1,1,4) + js_v*randn(1)); % new v_x for bb
        else
            c(1,1,3) = abs(c(1,1,3) + js_v*randn(1)); % new v_x for ba
            c(1,1,4) = c(1,1,3); % only one v_xb
        end
    end
    if fp.vaa
        c(2,2,1) = abs(c(2,2,1) + js_v*randn(1)); % new v_y for aa
    end
    if fp.vba
        if ~fp.vya
            c(2,2,3) = abs(c(2,2,3) + js_v*randn(1)); % new v_y for ba
        else
            c(2,2,3) = c(2,2,1);
        end
    end
    if fp.vab & fp.vbb
        if ~fp.vyb
            c(2,2,2) = abs(c(2,2,2) + js_v*randn(1)); % new v_y for ab
            c(2,2,4) = abs(c(2,2,4) + js_v*randn(1)); % new v_y for bb
        else
            c(2,2,2) = abs(c(2,2,2) + js_v*randn(1)); % new v_y for ab
            c(2,2,4) = c(2,2,2); % there can be only one
        end
    end

    if fp.raa
        %probe_r = 0;
        %while probe_r == 0
        raa_new = raa + js_r*randn(1);
        if raa_new <= -1
            raa_new = -.9999;
        elseif raa_new >= 1
            raa_new = .9999;
        end
        %if raa_new < 1 & raa_new > -1
        raa = raa_new;
        %probe_r = 1;
        %end
        %end            
        c(1,2,1) = raa*sqrt(c(1,1,1)*c(2,2,1));
        c(2,1,1) = c(1,2,1);
    end

    if fp.rab
        %probe_r = 0;
        %while probe_r == 0
        rab_new = rab + js_r*randn(1);
        if rab_new <= -1
            rab_new = -.9999;
        elseif rab_new >= 1
            rab_new = .9999;
        end
        %if raa_new < 1 & raa_new > -1
        rab = rab_new;
        %probe_r = 1;
        %end
        %end            
        c(1,2,2) = rab*sqrt(c(1,1,2)*c(2,2,2));
        c(2,1,2) = c(1,2,2);
    end

    if fp.rba
        %probe_r = 0;
        %while probe_r == 0
        rba_new = rba + js_r*randn(1);
        if rba_new <= -1
            rba_new = -.9999;
        elseif rba_new >= 1
            rba_new = .9999;
        end
        %if raa_new < 1 & raa_new > -1
        rba = rba_new;
        %probe_r = 1;
        %end    
        %end
        c(1,2,3) = rba*sqrt(c(1,1,3)*c(2,2,3));
        c(2,1,3) = c(1,2,3);
    end

    if fp.rbb
        %probe_r = 0;
        %while probe_r == 0
        rbb_new = rbb + js_r*randn(1);
        if rbb_new <= -1
            rbb_new = -.9999;
        elseif rbb_new >= 1
            rbb_new = .9999;
        end
        %if raa_new < 1 & raa_new > -1
        rbb = rbb_new;
        %probe_r = 1;
        %end
        %end            
        c(1,2,4) = rbb*sqrt(c(1,1,4)*c(2,2,4));
        c(2,1,4) = c(1,2,4);
    end
    
    % 'horizontal' decision bound parameters
    % some of this code depends explicitly on exactly four disributions
    if strcmp(dbtype,'maxpost')
        if mp_free
            dbp(1,:) = abs(dbp(1,:) + js_a*randn(1,4));
            dbp(1,:) = dbp(1,:)/sum(dbp(1,:));
            dbp(2,:) = abs(dbp(2,:) + js_a*randn(1,4));
            dbp(2,:) = dbp(2,:)/sum(dbp(2,:));
            dbp(3,:) = abs(dbp(3,:) + js_a*randn(1,4));
            dbp(3,:) = dbp(3,:)/sum(dbp(3,:));
            dbp(4,:) = abs(dbp(4,:) + js_a*randn(1,4));
            dbp(4,:) = dbp(4,:)/sum(dbp(4,:));
            dbp(5,:) = abs(dbp(5,:) + js_a*randn(1,4));
            dbp(5,:) = dbp(5,:)/sum(dbp(5,:));
        elseif mp_feat_g || mp_feat_r
            
            if mp_feat_g
                feat_params = abs(feat_params + js_a*randn(2,6));
            else
                feat_params(:,1:2) = abs(feat_params(:,1:2) + js_a*randn(2,2));
                feat_params(:,4:5) = abs(feat_params(:,4:5) + js_a*randn(2,2));
            end
            normalizer = sum(feat_params,1);
            for i = 1:2
               feat_params(i,:) = feat_params(i,:)./normalizer;
            end
            
            % feat_params structure
            % (1,1) = stim.1 condition 1 weight
            % (1,2) = stim1. condition 2 weight
            % (1,3) = stim.1 condition 3 ('neutral') weight
            % (1,4) = stim1. condition 4 weight
            % (1,5) = stim.1 condition 5 weight
            % (1,6) = stim1. condition 3 ('neutral') weight
            % (2,:) = 1-(1,:)
            
            dbp(1,1) = feat_params(1,1)*feat_params(1,6);
            dbp(1,2) = feat_params(2,1)*feat_params(1,6);
            dbp(1,3) = feat_params(1,1)*feat_params(2,6);
            dbp(1,4) = feat_params(2,1)*feat_params(2,6);
            dbp(1,:) = dbp(1,:)/sum(dbp(1,:));
            
            dbp(2,1) = feat_params(1,2)*feat_params(1,3);
            dbp(2,2) = feat_params(1,2)*feat_params(2,3);
            dbp(2,3) = feat_params(2,2)*feat_params(1,3);
            dbp(2,4) = feat_params(2,2)*feat_params(2,3);
            dbp(2,:) = dbp(2,:)/sum(dbp(2,:));
            
            dbp(3,1) = feat_params(1,3)*feat_params(1,6);
            dbp(3,2) = feat_params(2,3)*feat_params(1,6);
            dbp(3,3) = feat_params(1,3)*feat_params(2,6);
            dbp(3,4) = feat_params(2,3)*feat_params(2,6);
            dbp(3,:) = dbp(3,:)/sum(dbp(3,:));
            
            dbp(4,1) = feat_params(1,4)*feat_params(1,3);
            dbp(4,2) = feat_params(1,4)*feat_params(2,3);
            dbp(4,3) = feat_params(2,4)*feat_params(1,3);
            dbp(4,4) = feat_params(2,4)*feat_params(2,3);
            dbp(4,:) = dbp(4,:)/sum(dbp(4,:));

            dbp(5,1) = feat_params(1,5)*feat_params(1,6);
            dbp(5,2) = feat_params(2,5)*feat_params(1,6);
            dbp(5,3) = feat_params(1,5)*feat_params(2,6);
            dbp(5,4) = feat_params(2,5)*feat_params(2,6);
            dbp(5,:) = dbp(5,:)/sum(dbp(5,:));
                 
        elseif mp_br
            % nothing
        end
        %disp([dbp(1,1:4) sum(dbp(1,1:4))])
    elseif strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')
        if fp.ax
            dbp(1,1) = dbp(1,1) + js_a*randn(1);
            dbp(3,1) = dbp(1,1);
            dbp(5,1) = dbp(1,1);
        end
        if fp.bx
            dbp(1,2) = dbp(1,2) + js_b*randn(1);
            dbp(3,2) = dbp(1,2);
            dbp(5,2) = dbp(1,2);
        end
        if fp.cx
            if ~fp.cx_k
                dbp(1,3) = dbp(1,3) + js_c*randn(1);
                dbp(3,3) = dbp(3,3) + js_c*randn(1);
                dbp(5,3) = dbp(5,3) + js_c*randn(1);
            else
                dbp(3,3) = dbp(3,3) + js_c*randn(1);
                dbp(1,3) = dbp(3,3);
                dbp(5,3) = dbp(3,3);
            end            
        end
        if fp.dx
            if ~fp.dx_k
                dbp(1,4) = dbp(1,4) + js_d*randn(1);
                dbp(3,4) = dbp(3,4) + js_d*randn(1);
                dbp(5,4) = dbp(5,4) + js_d*randn(1);
            else
                dbp(3,4) = dbp(3,4) + js_d*randn(1);
                dbp(1,4) = dbp(3,4);
                dbp(5,4) = dbp(3,4);
            end
        end
        if fp.sx
            dbp(1,5) = dbp(1,5) - abs(js_s*randn(1));
            dbp(3,5) = dbp(1,5);
            dbp(5,5) = dbp(1,5);
        end
        
        % 'vertical' decision bound parameters
        if fp.ay
            dbp(2,1) = dbp(2,1) + js_a*randn(1);
            dbp(6,1) = dbp(2,1);
            dbp(4,1) = dbp(2,1);
        end
        if fp.by
            dbp(2,2) = dbp(2,2) + js_b*randn(1);
            dbp(6,2) = dbp(2,2);
            dbp(4,2) = dbp(2,2);
        end
        if fp.cy
            if ~fp.cy_k
                dbp(2,3) = dbp(2,3) + js_c*randn(1);
                dbp(6,3) = dbp(6,3) + js_c*randn(1);
                dbp(4,3) = dbp(4,3) + js_c*randn(1);
            else
                dbp(6,3) = dbp(6,3) + js_c*randn(1);
                dbp(2,3) = dbp(6,3);
                dbp(4,3) = dbp(6,3);
            end
        end
        if fp.dy
            if ~fp.dy_k
                dbp(2,4) = dbp(2,4) + js_d*randn(1);
                dbp(6,4) = dbp(6,4) + js_d*randn(1);
                dbp(4,4) = dbp(4,4) + js_d*randn(1);
            else
                dbp(6,4) = dbp(6,4) + js_d*randn(1);
                dbp(2,4) = dbp(6,4);
                dbp(4,4) = dbp(6,4);
            end
        end
        if fp.sy
            dbp(2,5) = dbp(2,5) - abs(js_s*randn(1));
            dbp(6,5) = dbp(2,5);
            dbp(4,5) = dbp(2,5);
        end
    end
    
    % a priori conditions placed on the parameters
    
    % mean vectors maa mab mba mbb
    if m(1,3) < m(1,1) || m(1,4) < m(1,2) || ...
            m(2,2) < m(2,1) || m(2,4) < m(2,3)
        pv(1) = 1;
    elseif m(1,2) < -3 || m(1,3) > 5 || m(1,4) > 5 || ...
        m(2,2) > 4 || m(2,3) < -3 || m(2,4) > 5
        pv(1) = 1;
    else
        pv(1) = 0;
    end
        
    % covariance, variance
    
    for i = 1:4
        if abs(c(1,2,i))>( sqrt(c(1,1,i)) * sqrt(c(2,2,i)) )
            pv(2) = 1;
            break
        elseif sqrt(c(1,1,i)) > 3 || sqrt(c(2,2,i)) > 3
                pv(2) = 1;
                break
        else
            pv(2) = 0;
        end
    end
    
    if ~strcmp(dbtype,'maxpost') & (dbp(1,5) > 0 || dbp(3,5) > 0 || dbp(5,5) > 0 ...
            || dbp(2,5) > 0 || dbp(6,5) > 0 || dbp(4,5) > 0)
        pv(4) = 1;
    else
        pv(4) = 0;
    end
    
    if pv(1)
        disp('parameter violation: means in disarray')
        LL_new = cf;
    elseif pv(2)
        disp('parameter violation: covariance too big')
        LL_new = cf;
    elseif pv(3)
        disp('parameter violation: variance too big')
        LL_new = cf;
    elseif pv(4)
        disp('parameter violation: bound offset > 0')        
        LL_new = cf;
    else
        
        % created bivar_norm function on 6/17/08
        [dist,x,y] = bivar_norm(m,c,bin);
        
        % created dec_reg function on 6/17/08
        dec = logical(dec_reg(dbtype,dist,x,y,dbp,nbr));
                
        pM = zeros(nstim,nresp,nbr);
        
        binsize = bin^2;
        
        %disp(['calculating predicted confusion probabilities'])
        
        for id = 1:nbr % base rate condition
            for s = 1:nstim % stimulus
                pdf = dist(:,:,s);
                for r = 1:nresp % response
                    pM(s,r,id) = max(binsize*sum( sum( pdf( dec(:,:,r,id) ) ) ),eps);
                end
            end
        end
        
        for id = 1:nbr
            for s = 1:nstim
                pM(s,:,id) = pM(s,:,id)/sum(pM(s,:,id));
            end
        end
    
        efM = zeros(nstim,nresp,nbr);
        for ci = 1:nbr
            for i = 1:nstim
                efM(i,:,ci) = pM(i,:,ci)*Nr_oM(i,ci);
            end
        end

        LL_M = zeros(nstim,nresp,nbr);
        LL_M = oM.*log(pM); 
        LL_new = sum(sum(sum(LL_M(:,:,:))));
        bic_new = -2*LL_new + complexity;
        g2_new = 2*sum(sum(sum( oM.*log(oM./efM) ))); % 2*(LL_sat - LL_new);
        
    end % if pv(i)==1, etc...

    if n==0
        g2_min = g2_new;
        bic_min = bic_new;
        LL_max = LL_new;
        bic_old = bic_new;
        g2_old = g2_new;
        LL_old = LL_new;
        n = n+1;
        best_pM = pM;
        best_fitps.maa = m(:,1);
        best_fitps.mab = m(:,2);
        best_fitps.mba = m(:,3);
        best_fitps.mbb = m(:,4);
        best_fitps.caa = c(:,:,1);
        best_fitps.cab = c(:,:,2);
        best_fitps.cba = c(:,:,3);
        best_fitps.cbb = c(:,:,4);
        if strcmp(dbtype,'maxpost')
            best_fitps.c1aa = dbp(1,1);
            best_fitps.c1ab = dbp(1,2);
            best_fitps.c1ba = dbp(1,3);
            best_fitps.c1bb = dbp(1,4);
            best_fitps.c2aa = dbp(2,1);
            best_fitps.c2ab = dbp(2,2);
            best_fitps.c2ba = dbp(2,3);
            best_fitps.c2bb = dbp(2,4);
            best_fitps.c3aa = dbp(3,1);
            best_fitps.c3ab = dbp(3,2);
            best_fitps.c3ba = dbp(3,3);
            best_fitps.c3bb = dbp(3,4);
            best_fitps.c4aa = dbp(4,1);
            best_fitps.c4ab = dbp(4,2);
            best_fitps.c4ba = dbp(4,3);
            best_fitps.c4bb = dbp(4,4);
            best_fitps.c5aa = dbp(5,1);
            best_fitps.c5ab = dbp(5,2);
            best_fitps.c5ba = dbp(5,3);
            best_fitps.c5bb = dbp(5,4); 
        else
            best_fitps.ax = dbp(1,1);
            best_fitps.bx = dbp(1,2);
            best_fitps.cx1 = dbp(1,3);
            best_fitps.cx3 = dbp(3,3);
            best_fitps.cx5 = dbp(5,3);
            best_fitps.dx1 = dbp(1,4);
            best_fitps.dx3 = dbp(3,4);
            best_fitps.dx5 = dbp(5,4);
            best_fitps.sx = dbp(1,5);
            best_fitps.ay = dbp(2,1);
            best_fitps.by = dbp(2,2);
            best_fitps.cy2 = dbp(2,3);
            best_fitps.cy3 = dbp(6,3);
            best_fitps.cy4 = dbp(4,3);
            best_fitps.dy2 = dbp(2,4);
            best_fitps.dy3 = dbp(6,4);
            best_fitps.dy4 = dbp(4,4);
            best_fitps.sy = dbp(2,5);
        end
        best_fitps.g2 = g2_new;
        best_fitps.bic = bic_new;
        best_fitps.LL = LL_new;
    else
        
        if bic_new < bic_old
            accept = 1;
        end
        
        if bic_new < bic_min
            g2_min = g2_new;
            bic_min = bic_new;
            LL_max = LL_new;
            best_pM = pM;
            best_fitps.maa = m(:,1);
            best_fitps.mab = m(:,2);
            best_fitps.mba = m(:,3);
            best_fitps.mbb = m(:,4);
            best_fitps.caa = c(:,:,1);
            best_fitps.cab = c(:,:,2);
            best_fitps.cba = c(:,:,3);
            best_fitps.cbb = c(:,:,4);
            if strcmp(dbtype,'maxpost')
                best_fitps.c1aa = dbp(1,1);
                best_fitps.c1ab = dbp(1,2);
                best_fitps.c1ba = dbp(1,3);
                best_fitps.c1bb = dbp(1,4);
                best_fitps.c2aa = dbp(2,1);
                best_fitps.c2ab = dbp(2,2);
                best_fitps.c2ba = dbp(2,3);
                best_fitps.c2bb = dbp(2,4);
                best_fitps.c3aa = dbp(3,1);
                best_fitps.c3ab = dbp(3,2);
                best_fitps.c3ba = dbp(3,3);
                best_fitps.c3bb = dbp(3,4);
                best_fitps.c4aa = dbp(4,1);
                best_fitps.c4ab = dbp(4,2);
                best_fitps.c4ba = dbp(4,3);
                best_fitps.c4bb = dbp(4,4);
                best_fitps.c5aa = dbp(5,1);
                best_fitps.c5ab = dbp(5,2);
                best_fitps.c5ba = dbp(5,3);
                best_fitps.c5bb = dbp(5,4);                
            elseif strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')
                best_fitps.ax = dbp(1,1);
                best_fitps.bx = dbp(1,2);
                best_fitps.cx1 = dbp(1,3);
                best_fitps.cx3 = dbp(3,3);
                best_fitps.cx5 = dbp(5,3);
                best_fitps.dx1 = dbp(1,4);
                best_fitps.dx3 = dbp(3,4);
                best_fitps.dx5 = dbp(5,4);
                best_fitps.sx = dbp(1,5);
                best_fitps.ay = dbp(2,1);
                best_fitps.by = dbp(2,2);
                best_fitps.cy2 = dbp(2,3);
                best_fitps.cy3 = dbp(6,3);
                best_fitps.cy4 = dbp(4,3);
                best_fitps.dy2 = dbp(2,4);
                best_fitps.dy3 = dbp(6,4);
                best_fitps.dy4 = dbp(4,4);
                best_fitps.sy = dbp(2,5);
            end
            best_fitps.g2 = g2_new;
            best_fitps.bic = bic_new;
            best_fitps.LL = LL_new;
        end
    end
    
    if accept
        
        % put the parameters in the fitps struct
        fitps(n).maa = m(:,1);
        fitps(n).mab = m(:,2);
        fitps(n).mba = m(:,3);
        fitps(n).mbb = m(:,4);
        fitps(n).caa = c(:,:,1);
        fitps(n).cab = c(:,:,2);
        fitps(n).cba = c(:,:,3);
        fitps(n).cbb = c(:,:,4);
        if strcmp(dbtype,'maxpost')
            fitps(n).c1aa = dbp(1,1);
            fitps(n).c1ab = dbp(1,2);
            fitps(n).c1ba = dbp(1,3);
            fitps(n).c1bb = dbp(1,4);
            fitps(n).c2aa = dbp(2,1);
            fitps(n).c2ab = dbp(2,2);
            fitps(n).c2ba = dbp(2,3);
            fitps(n).c2bb = dbp(2,4);
            fitps(n).c3aa = dbp(3,1);
            fitps(n).c3ab = dbp(3,2);
            fitps(n).c3ba = dbp(3,3);
            fitps(n).c3bb = dbp(3,4);
            fitps(n).c4aa = dbp(4,1);
            fitps(n).c4ab = dbp(4,2);
            fitps(n).c4ba = dbp(4,3);
            fitps(n).c4bb = dbp(4,4);
            fitps(n).c5aa = dbp(5,1);
            fitps(n).c5ab = dbp(5,2);
            fitps(n).c5ba = dbp(5,3);
            fitps(n).c5bb = dbp(5,4);                
        else
            fitps(n).ax = dbp(1,1);
            fitps(n).bx = dbp(1,2);
            fitps(n).cx1 = dbp(1,3);
            fitps(n).cx3 = dbp(3,3);
            fitps(n).cx5 = dbp(5,3);
            fitps(n).dx1 = dbp(1,4);
            fitps(n).dx3 = dbp(3,4);
            fitps(n).dx5 = dbp(5,4);
            fitps(n).sx = dbp(1,5);
            fitps(n).ay = dbp(2,1);
            fitps(n).by = dbp(2,2);
            fitps(n).cy2 = dbp(2,3);
            fitps(n).cy3 = dbp(6,3);
            fitps(n).cy4 = dbp(4,3);
            fitps(n).dy2 = dbp(2,4);
            fitps(n).dy3 = dbp(6,4);
            fitps(n).dy4 = dbp(4,4);
            fitps(n).sy = dbp(2,5);
        end
        fitps(n).g2 = g2_new;
        fitps(n).bic = bic_new;
        fitps(n).LL = LL_new;
        fitps(n).nrej = n_rejected;
        
        bic_old = bic_new;
        g2_old = g2_new;
        LL_old = LL_new;
        n = n+1; % increment the index for the current fit
        n_rejected = 0;
    else
        m = m_old;
        c = c_old;
        if strcmp(dbtype,'maxpost')
            dbp(1,:) = dbp_old(1,:);
            dbp(3,:) = dbp_old(3,:);
            dbp(5,:) = dbp_old(5,:);
            dbp(2,:) = dbp_old(2,:);
            dbp(4,:) = dbp_old(4,:);
        else
            dbp(1,:) = dbp_old(1,:);
            dbp(3,:) = dbp_old(3,:);
            dbp(5,:) = dbp_old(5,:);
            dbp(2,:) = dbp_old(2,:);
            dbp(6,:) = dbp_old(6,:);
            dbp(4,:) = dbp_old(4,:);
        end
        n_rejected = n_rejected + 1;
    end
    %disp(['LL_sat = ' num2str(LL_sat)])
    %disp(['LL_max = ' num2str(LL_max)])
    %disp(['LL_old = ' num2str(LL_old)])
    %disp(['LL_new = ' num2str(LL_new)])
    if n==1 || n==50 || n==100 || n==150
        disp(['complexity = ' num2str(complexity)])
    end
    disp(['bic_min = ' num2str(bic_min)])
    disp(['bic_old = ' num2str(bic_old)])
    disp(['bic_new = ' num2str(bic_new)])
    %disp(['LL_dif = ' num2str(LL_sat - LL_new)])
    
    disp(['n_rejected = ' num2str(n_rejected)])
    disp(' ')
end