function [outfit] = grtnewtrap(freq)

% function [outfit] = grtnewtrap(freq);

% single 4x4 matrix; assumes bounds are coordinate axes
% assumes row/col order: aa, ab, ba, bb

delta = 1/10000; % adjust as needed

% cols: 1 x-sep; 2 y-sep; 3-6 indep. for aa, ab, ba, bb, respectively
% to be extended to allow (failure of) indep in each dist.
% and single correlation across dists

mods = [0 0 0 0 0 0 0; 1 0 0 0 0 0 0; 0 1 0 0 0 0 0; 1 1 0 0 0 0 0; ...
    0 0 0 0 0 0 1; 1 0 0 0 0 0 1; 0 1 0 0 0 0 1; 1 1 0 0 0 0 1; ...
    0 0 1 1 1 1 0; 1 0 1 1 1 1 0; 0 1 1 1 1 1 0; 1 1 1 1 1 1 0];

%mods = [0 0 0 0 0 0 0; 0 1 0 0 0 0 0; 1 0 0 0 0 0 0];

[nmod,nrst] = size(mods);

% comment w.r.t. restricted models:
%  currently maintains length 12 d vector and 12x12 info matrix
%  with repeated elements, as necessary; may need to use shorter vectors
%  and a smaller info matrix instead of repeated elements

% protection against zeros, which seem to cause the algorithm to explode
freq = freq+1;

% parameters in order:
% mu_x_** mu_y_** rho_**
% ** = aa, ab, ba, bb

% initialization

% notes:
% should initialize correlation parameters in more informative way
% need to figure right out criterion for continued iterations

[nrow,ncol] = size(freq);

for m = 1:nmod
    
    w = 1;
    
    prob = zeros(nrow,ncol);
    
    for i = 1:nrow
        prob(i,:) = freq(i,:)./sum(freq(i,:));
    end
    
    rows = [];
    xpar = [];
    ypar = [];
    rpar = [];
    nx = 0;
    ny = 0;
    nr = 0;
    xc_aa = 0;
    xc_ab = 0;
    xc_ba = 0;
    xc_bb = 0;
    yc_aa = 0;
    yc_ab = 0;
    yc_ba = 0;
    yc_bb = 0;
    
    if mods(m,1)==1
        xpar = [1:2];
        nx = 2;
        rows(xpar,:) = [1 1 0 0; 0 0 1 1];
    else
        xpar = [1:4];
        nx = 4;
        rows(xpar,:) = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    end
    if mods(m,2)==1
        ypar = nx + [1:2];
        ny = 2;
        rows(ypar,:) = [1 0 1 0; 0 1 0 1];
    else
        ypar = nx + [1:4];
        ny = 4;
        rows(ypar,:) = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    end
    if mods(m,7)==1
        rpar = nx + ny + 1;
        nr = 1;
        rows(nx + ny + 1,:) = [1 1 1 1];
    else
        if mods(m,3)==0
            rpar(end+1) = nx + ny + 1;
            rows(nx + ny + nr + 1,:) = [1 0 0 0];
            nr = nr + 1;
        end
        if mods(m,4)==0
            rpar(end+1) = nx + ny + nr + 1;
            rows(nx + ny + nr + 1,:) = [0 1 0 0];
            nr = nr + 1;
        end
        if mods(m,5)==0
            rpar(end+1) = nx + ny + nr + 1;
            rows(nx + ny + nr + 1,:) = [0 0 1 0];
            nr = nr + 1;
        end
        if mods(m,6)==0
            rpar(end+1) = nx + ny + nr + 1;
            rows(nx + ny + nr + 1,:) = [0 0 0 1];
            nr = nr + 1;
        end
    end
    npar = nx + ny + nr;
    d = zeros(npar,1);
    ps_old = zeros(npar,1);
    E = zeros(npar);
    
    rows = logical(rows);
    
    % stim aa means
    
    if mods(m,1)==1
        xc_aa = ...
            norminv(.5*(prob(1,1)+prob(1,2)) + .5*(prob(2,1)+prob(2,2)));
    else
        xc_aa = norminv(prob(1,1) + prob(1,2));
    end
    
    if mods(m,2)==1
        yc_aa = ...
            norminv(.5*(prob(1,1)+prob(1,3)) + .5*(prob(3,1)+prob(3,3)));
    else
        yc_aa = norminv(prob(1,1) + prob(1,3));
    end
    
    ps_old(xpar(1)) = -xc_aa;
    ps_old(ypar(1)) = -yc_aa;
    
    
    % stim ab means
    
    if mods(m,1)==0
        xc_ab = norminv(prob(2,1) + prob(2,2));
        ps_old(xpar(2)) = -xc_ab;
    end
    
    if mods(m,2)==1
        yc_ab = ...
            norminv(.5*(prob(2,1)+prob(2,3)) + .5*(prob(4,1)+prob(4,3)));
    else
        yc_ab = norminv(prob(2,1) + prob(2,3));
    end
    ps_old(ypar(2)) = -yc_ab;
    
    
    % stim ba means
    
    if mods(m,1)==1
        xc_ba = ...
            norminv(.5*(prob(3,1)+prob(3,2)) + .5*(prob(4,1)+prob(4,2)));
        ps_old(xpar(2)) = -xc_ba;
    else
        xc_ba = norminv(prob(3,1) + prob(3,2));
        ps_old(xpar(3)) = -xc_ba;
    end
    
    if mods(m,2)==0
        yc_ba = norminv(prob(3,1) + prob(3,3));
        ps_old(ypar(3)) = -yc_ba;
    end
    
    
    % stim bb means
    
    if mods(m,1)==0
        xc_bb = norminv(prob(4,1) + prob(4,2));
        ps_old(xpar(4)) = -xc_bb;
    end
    
    if mods(m,2)==0
        yc_bb = norminv(prob(4,1) + prob(4,3));
        ps_old(ypar(4)) = -yc_bb;
    end
    
    
    % correlations parameters
    
    % single correlation parameter across distributions
    if mods(m,7)==1
        a = .25*sum(prob(:,4));
        b = .25*sum(prob(:,3));
        c = .25*sum(prob(:,2));
        e = .25*sum(prob(:,1));
        r = cos(pi/(1+sqrt((a*e)/(b*c))));
        if r <= -1
            r = -.95;
        elseif r >= 1
            r = .95;
        end
        ps_old(rpar) = r;
    else
        % aa
        if mods(m,3)==0
            a = prob(1,4);
            b = prob(1,3);
            c = prob(1,2);
            e = prob(1,1);
            r = cos(pi/(1+sqrt((a*e)/(b*c))));
            if r <= -1
                r = -.95;
            elseif r >= 1
                r = .95;
            end
            rho_aa = r;
        else
            rho_aa = [];
        end
        % ab
        if mods(m,4)==0
            a = prob(2,4);
            b = prob(2,3);
            c = prob(2,2);
            e = prob(2,1);
            r = cos(pi/(1+sqrt((a*e)/(b*c))));
            if r <= -1
                r = -.95;
            elseif r >= 1
                r = .95;
            end
            rho_ab = r;
        else
            rho_ab = [];
        end
        % ba
        if mods(m,5)==0
            a = prob(3,4);
            b = prob(3,3);
            c = prob(3,2);
            e = prob(3,1);
            r = cos(pi/(1+sqrt((a*e)/(b*c))));
            if r <= -1
                r = -.95;
            elseif r >= 1
                r = .95;
            end
            rho_ba = r;
        else
            rho_ba = [];
        end
        % bb
        if mods(m,6)==0
            a = prob(4,4);
            b = prob(4,3);
            c = prob(4,2);
            e = prob(4,1);
            r = cos(pi/(1+sqrt((a*e)/(b*c))));
            if r <= -1
                r = -.95;
            elseif r >= 1
                r = .95;
            end
            rho_bb = r;
        else
            rho_bb = [];
        end
        ps_old(rpar) = [rho_aa; rho_ab; rho_ba; rho_bb];
    end
    
    % stimulus aa
    
    alpha = ps_old(xpar(1)); % x mean, stim aa
    kappa = ps_old(ypar(1)); % y mean, stim aa
    if mods(m,3)==1
        rho = 0;
    else % mods(m,3)==0 or mods(m,7)==1
        rho = ps_old(rpar(1));
    end
    
    prob(1,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    v(1,:,:) = vcalc(alpha,kappa,rho);
    
    % stimulus ab
    
    if mods(m,1)==1
        alpha = ps_old(xpar(1)); % x mean, stim ab
    else
        alpha = ps_old(xpar(2));
    end
    kappa = ps_old(ypar(2)); % y mean, stim ab
    if mods(m,7)==1
        rho = ps_old(rpar(1));
    elseif mods(m,4)==0
        offset = 0+[mods(m,3)==0];
        rho = ps_old(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(2,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    v(2,:,:) = vcalc(alpha,kappa,rho);
    
    % stimulus ba
    
    if mods(m,1)==1
        alpha = ps_old(xpar(2));
    else
        alpha = ps_old(xpar(3)); % x mean, stim ba
    end
    if mods(m,2)==1
        kappa = ps_old(ypar(1));
    else
        kappa = ps_old(ypar(3)); % y mean, stim ba
    end
    if mods(m,7)==1
        rho = ps_old(rpar(1));
    elseif mods(m,5)==0
        offset = sum(mods(m,3:4)==0);
        rho = ps_old(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(3,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    v(3,:,:) = vcalc(alpha,kappa,rho);
    
    % stimulus bb
    
    if mods(m,1)==1
        alpha = ps_old(xpar(2));
    else
        alpha = ps_old(xpar(4)); % x mean, stim bb
    end
    if mods(m,2)==1
        kappa = ps_old(ypar(2));
    else
        kappa = ps_old(ypar(4)); % y mean, stim bb
    end
    if mods(m,7)==1
        rho = ps_old(rpar(1));
    elseif mods(m,6)==0
        offset = sum(mods(m,3:5)==0);
        rho = ps_old(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(4,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    v(4,:,:) = vcalc(alpha,kappa,rho);
    
    % log-likelihood gradient, information matrix
    for i = 1:npar
        if sum(i==xpar)
            g = 1;
        elseif sum(i==ypar)
            g = 2;
        elseif sum(i==rpar)
            g = 3;
        end
        ir = rows(i,:);
        d(i) = sum(sum( (freq(ir,:)./prob(ir,:)).*v(ir,:,g) ));
        for j = 1:npar
            if sum(j==xpar)
                h = 1;
            elseif sum(j==ypar)
                h = 2;
            elseif sum(j==rpar)
                h = 3;
            end
            jr = rows(j,:);
            
            ijr = ir==1 & jr==1;
            
            if sum(ijr) > 0
                
                K = sum( freq(ijr,:) , 2 )./sum( prob(ijr,:) , 2 );
                
                L = sum( v(ijr,:,g).*v(ijr,:,h)./prob(ijr,:) , 2 ) - ...
                    sum(v(ijr,:,g),2).*sum(v(ijr,:,h),2)./sum(prob(ijr,:),2);
                
                E(i,j) = -K'*L;
            end
        end
    end
    
    Ei = E\eye(npar);
    
    ps_new = ps_old - Ei*d;
    
    % iterate!
    
    it = 1;
    
    df = abs( ps_new-ps_old )./abs(ps_new);
    
    dfp_new = df'*df;
    
    dfp_old = dfp_new;
    
    while dfp_new > delta
        
        disp([m it])
        
        ps_old = ps_new;
        
        % stimulus aa
        
        % stimulus aa
        
        alpha = ps_old(xpar(1)); % x mean, stim aa
        kappa = ps_old(ypar(1)); % y mean, stim aa
        if mods(m,3)==1
            rho = 0;
        else % mods(m,3)==0 or mods(m,7)==1
            rho = ps_old(rpar(1));
        end
        
        prob(1,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
        
        v(1,:,:) = vcalc(alpha,kappa,rho);
        
        % stimulus ab
        
        if mods(m,1)==1
            alpha = ps_old(xpar(1)); % x mean, stim ab
        else
            alpha = ps_old(xpar(2));
        end
        kappa = ps_old(ypar(2)); % y mean, stim ab
        if mods(m,7)==1
            rho = ps_old(rpar(1));
        elseif mods(m,4)==0
            offset = 0+[mods(m,3)==0];
            rho = ps_old(rpar(offset+1));
        else
            rho = 0;
        end
        
        prob(2,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
        
        v(2,:,:) = vcalc(alpha,kappa,rho);
        
        % stimulus ba
        
        if mods(m,1)==1
            alpha = ps_old(xpar(2));
        else
            alpha = ps_old(xpar(3)); % x mean, stim ba
        end
        if mods(m,2)==1
            kappa = ps_old(ypar(1));
        else
            kappa = ps_old(ypar(3)); % y mean, stim ba
        end
        if mods(m,7)==1
            rho = ps_old(rpar(1));
        elseif mods(m,5)==0
            offset = sum(mods(m,3:4)==0);
            rho = ps_old(rpar(offset+1));
        else
            rho = 0;
        end
        
        prob(3,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
        
        v(3,:,:) = vcalc(alpha,kappa,rho);
        
        % stimulus bb
        
        if mods(m,1)==1
            alpha = ps_old(xpar(2));
        else
            alpha = ps_old(xpar(4)); % x mean, stim bb
        end
        if mods(m,2)==1
            kappa = ps_old(ypar(2));
        else
            kappa = ps_old(ypar(4)); % y mean, stim bb
        end
        if mods(m,7)==1
            rho = ps_old(rpar(1));
        elseif mods(m,6)==0
            offset = sum(mods(m,3:5)==0);
            rho = ps_old(rpar(offset+1));
        else
            rho = 0;
        end
        
        prob(4,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
        
        v(4,:,:) = vcalc(alpha,kappa,rho);
        
        % log-likelihood gradient, information matrix
        for i = 1:npar
            if sum(i==xpar)
                g = 1;
            elseif sum(i==ypar)
                g = 2;
            elseif sum(i==rpar)
                g = 3;
            end
            ir = rows(i,:);
            d(i) = sum(sum( v(ir,:,g).*freq(rows(i,:),:)./prob(rows(i,:),:) ));
            for j = 1:npar
                if sum(j==xpar)
                    h = 1;
                elseif sum(j==ypar)
                    h = 2;
                elseif sum(j==rpar)
                    h = 3;
                end
                jr = rows(j,:);
                if (ir+0)*(jr+0)' > 0
                    
                    ijr = ir==1 & jr==1;
                    
                    K = -1*(sum( freq(ijr,:) , 2 )./sum( prob(ijr,:) , 2 ));
                    
                    L = sum( v(ijr,:,g).*v(ijr,:,h)./prob(ijr,:) , 2 ) - ...
                        sum(v(ijr,:,g),2).*sum(v(ijr,:,h),2)./sum(prob(ijr,:),2);
                    
                    E(i,j) = K'*L;
                end
            end
        end
        
        Ei = E\eye(npar);
        
        if dfp_new > dfp_old %w halving procedure
            w = .5*w;
        end
        
        ps_new = ps_old - w*Ei*d;
        
        df = abs(ps_new-ps_old)./abs(ps_new);
        
        %        disp(df'*df)
        
        dfp_old = dfp_new;
        
        dfp_new = df'*df;
        
        it = it + 1;
        
    end
    
    %outfit.d(:,m) = zeros(12,1);
    %outfit.d(1:npar,m) = d;
    if mods(m,1)==1
        outfit.p(1,m) = ps_new(1);
        outfit.p(2,m) = ps_new(1);
        outfit.p(3,m) = ps_new(2);
        outfit.p(4,m) = ps_new(2);
    else
        outfit.p(1,m) = ps_new(1);
        outfit.p(2,m) = ps_new(2);
        outfit.p(3,m) = ps_new(3);
        outfit.p(4,m) = ps_new(4);
    end
    if mods(m,2)==1
        outfit.p(5,m) = ps_new(nx+1);
        outfit.p(6,m) = ps_new(nx+2);
        outfit.p(7,m) = ps_new(nx+1);
        outfit.p(8,m) = ps_new(nx+2);
    else
        outfit.p(5,m) = ps_new(nx+1);
        outfit.p(6,m) = ps_new(nx+2);
        outfit.p(7,m) = ps_new(nx+3);
        outfit.p(8,m) = ps_new(nx+4);
    end
    if mods(m,7)==1
        outfit.p(9:12,m) = ps_new(nx+ny+1)*ones(4,1);
    else
        if mods(m,3)==1
            outfit.p(9,m) = 0;
        else
            outfit.p(9,m) = ps_new(nx+ny+1);
        end
        if mods(m,4)==1
            outfit.p(10,m) = 0;
        else
            offset = 0+mods(m,3)==0;
            outfit.p(10,m) = ps_new(nx+ny+offset+1);
        end
        if mods(m,5)==1
            outfit.p(11,m) = 0;
        else
            offset = sum(mods(m,3:4)==0);
            outfit.p(11,m) = ps_new(nx+ny+offset+1);
        end
        if mods(m,6)==1
            outfit.p(12,m) = 0;
        else
            offset = sum(mods(m,3:5)==0);
            outfit.p(12,m) = ps_new(nx+ny+offset+1);
        end
    end
    
    outfit.M(:,:,m) = zeros(12);
    infoM = (-E)\eye(npar);
    outfit.M(1:npar,1:npar,m) = infoM;
    
    % for calculating log-likelihood, AIC, BIC
    
    alpha = ps_new(xpar(1)); % x mean, stim aa
    kappa = ps_new(ypar(1)); % y mean, stim aa
    if mods(m,3)==1
        rho = 0;
    else % mods(m,3)==0 or mods(m,7)==1
        rho = ps_new(rpar(1));
    end
    
    prob(1,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    if mods(m,1)==1
        alpha = ps_new(xpar(1)); % x mean, stim ab
    else
        alpha = ps_new(xpar(2));
    end
    kappa = ps_new(ypar(2)); % y mean, stim ab
    if mods(m,7)==1
        rho = ps_new(rpar(1));
    elseif mods(m,4)==0
        offset = 0+[mods(m,3)==0];
        rho = ps_new(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(2,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    if mods(m,1)==1
        alpha = ps_new(xpar(2));
    else
        alpha = ps_new(xpar(3)); % x mean, stim ba
    end
    if mods(m,2)==1
        kappa = ps_new(ypar(1));
    else
        kappa = ps_new(ypar(3)); % y mean, stim ba
    end
    if mods(m,7)==1
        rho = ps_new(rpar(1));
    elseif mods(m,5)==0
        offset = sum(mods(m,3:4)==0);
        rho = ps_new(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(3,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    if mods(m,1)==1
        alpha = ps_new(xpar(2));
    else
        alpha = ps_new(xpar(4)); % x mean, stim bb
    end
    if mods(m,2)==1
        kappa = ps_new(ypar(2));
    else
        kappa = ps_new(ypar(4)); % y mean, stim bb
    end
    if mods(m,7)==1
        rho = ps_new(rpar(1));
    elseif mods(m,6)==0
        offset = sum(mods(m,3:5)==0);
        rho = ps_new(rpar(offset+1));
    else
        rho = 0;
    end
    
    prob(4,:) = prcalc([alpha kappa],eye(2) + [0 rho; rho 0]);
    
    outfit.pr = prob;
    
    loglike = sum(sum(freq.*log(prob)));
    outfit.nll(m) = -loglike;
    outfit.aic(m) = 2*npar - 2*loglike;
    outfit.bic(m) = npar*log(sum(sum(freq))) - 2*loglike;
    outfit.icomp(m) = -loglike + (npar/2)*log(trace(infoM)/npar)-.5*log(det(infoM));
end

% calculate predicted probabilities
function pr = prcalc(mv,SM);

pr(1,1) = mvncdf([-inf -inf],[0 0],mv,SM);
pr(1,2) = mvncdf([-inf 0],[0 inf],mv,SM);
pr(1,3) = mvncdf([0 -inf],[inf 0],mv,SM);
pr(1,4) = mvncdf([0 0],[inf inf],mv,SM);

% calculate v-matrix elements
function ve = vcalc(ap,kp,rh)

d_mx_arg = (rh*ap-kp)/sqrt(1-rh^2);
d_mx(1,1) = -normpdf(-ap)*normcdf( d_mx_arg );
d_mx(1,2) = -normpdf(-ap);
d_mx(2,1) = 0;
d_mx(2,2) = 0;

ve(1,1,1) =  d_mx(1,1);
ve(1,2,1) =  d_mx(1,2) - d_mx(1,1);
ve(1,3,1) =  d_mx(2,1) - d_mx(1,1);
ve(1,4,1) =  d_mx(2,2) - d_mx(2,1) - d_mx(1,2) + d_mx(1,1);

d_my_arg = (rh*kp-ap)/sqrt(1-rh^2);
d_my(1,1) = -normpdf(-kp)*normcdf( d_my_arg );
d_my(1,2) = 0;
d_my(2,1) = -normpdf(-kp);
d_my(2,2) = 0;

ve(1,1,2) =  d_my(1,1);
ve(1,2,2) =  d_my(1,2) - d_my(1,1);
ve(1,3,2) =  d_my(2,1) - d_my(1,1);
ve(1,4,2) =  d_my(2,2) - d_my(1,2) - d_my(2,1) + d_my(1,1);

S_aa = [1 rh; rh 1];
d_rho = mvnpdf([-ap -kp],[0 0],S_aa);

ve(1,1,3) =  d_rho;
ve(1,2,3) = -d_rho;
ve(1,3,3) = -d_rho;
ve(1,4,3) =  d_rho;