[A,info]=greyboxbuild(condfn,resfn,data,vndata,optn)

function [A,info]=greyboxbuild(condfn,resfn,data,vndata,optn)
% greyboxbuild  Greybox build an improved model
%  2012-01-24  Matlab8  Copyright (c) 2013, W J Whiten  BSD License
%
% [pr,rmse1,rmse2,info]=greyboxbuild(resfn,condfn,data,vndata,optn)
%  condfn  Condition function  
%           [c,pscale,pbase,pvalue]=condfn(data,idata)
%            data  User supplied data in any form
%            idata Data set number
%            c      Experimental condition column vector: data set idata
%            pscale Vector of scale (ie typical nonzero size) values 
%                    for parameters (optional if optn.Amat given)
%            pbase  Vector of base values for parameters (optional:
%                    default zero)
%  resfn   Residue calculation  [r,dr]=resfn(data,idata,p)
%            data  User supplied data in any form
%            idata Data set number
%            p     Parameter values
%            r     Column vector of residues for data set nbr idata
%            dr    Derivatives of r wrt p (ith column wrt p(i)) (optional)
%  data    Data for use in resfn & condfn (possibly a struct)
%  vndata  Vector of selected data set indices (typically 1:ndata)
%  optn    Optional additional information:
%           .Amat  Initial value for A matrix
%           .Acode Code value for use of each A matrix element
%                   0 do not use, 1 always use, 2 select
%           .mres  Estimate of total number of residuals
%           .reg   Ratio of largest eigenvalue for reject in linfitreg
%           .maxitr1 Maximum iterations for outer loop
%           .maxitr2 Maximum iterations for inner loop
%           .print   0 not printing, 1 final result, 2 partial results
%           .conv    Test level for convergence
%           .linfit  Function to use for linear fitting (default
%           .linoptn Struct options for linear fitting
%              .reg   Ratio of largest eigenvalue for reject in linfit
%              .tval  t ratio for rejection in linfitreg
%
%  pr      Estimated probability of reduction using random data
%  rmse0   Root mean squared error of unadjusted model
%  rmse2   Root mean squared error of adjusted model
%  info    Struct giving:
%           .Amat  Calculated A matrix value 
%           .Acode Matrix code for terms in A matrix (>0 use term)
%           .reg   Regularisation fraction for rejection of eigenvalues
%           .nreg  Number of eigenvectors used in regression
%           .ncond Number of condition terms
%           .npar  Number of parameters
%           .nres  Toal number of residuals
%           .conv  0 did not converge, 1 normal convergence,
%                   2 unable to improve current estimate.

% default option values
if(nargin<5)
    optn=struct();
end
optn0=struct('Amat',[],'Acode',[],'nres',1000,  ...
    'conv',1e-10,'maxitr1',50,'maxitr2',20,'print',2,  ...
    'linfit',@linfitregvar,'linoptn',struct());
optn=optndfts(optn,optn0);
A=optn.Amat;
Acode=optn.Acode;
linoptn=optn.linoptn;

optncore=selectfields(optn,'Amat','Acode','nres');

if(optn.print>0)
    fprintf('\nIter #   RMS error   Inner loop #\n')
end
indconv=0;
for i=1:optn.maxitr1
    optncore.Amat=A;
    [X,y,r1,info]=greyboxcore(condfn,resfn,data,vndata,optncore);
    
    if(optn.print>1 && i==1)
        fprintf('%4i %#15.9g\n',0,rms(r1))
    end
    
    if(i==1)
        ssq0=r1'*r1;
        if(isempty(optn.Acode))
            Acode=2*ones(info.npar,info.ncond);
        end
        
        linoptn.sel=Acode(Acode>0);

    end
    
    % solve for result matrix, default linfitregvar
    [a2,vm,sel,infolf]=optn.linfit(X,y,linoptn);
    a1=zeros(size(X,2),1);
    a1(sel)=a2;

    % form new A matrix
    if(isempty(A))
        A=zeros(info.npar,info.ncond);
    end
    A(Acode>0)=a1;

    if(optn.maxitr1>1)
        
        for j=1:optn.maxitr2
            r2=zeros(info.nres,1);
            k1=0;
            for kk=1:length(vndata)
                k=vndata(kk);
                
                par=parcalc(condfn,data,k,A);
                r2a=resfn(data,k,par);
                k1a=k1+1;
                k1=k1+numel(r2a);
                r2(k1a:k1)=r2a(:);
            end
            
            r12=(r1+r2)'*(r1-r2);
            if(abs(r12)/sum((abs(r1)+abs(r2)).^2)<optn.conv)
                indconv=1;
                break
            end
            
            if(optn.maxitr1==1)
                break
            end
            
            if(r12<0 && ~isempty(optn.Amat))
                A=(optn.Amat+A)/2;
            else
                break
            end
            
        end
        
        if(j==optn.maxitr2)
            indconv=2;
            break
        end
        
        if(indconv>0)
            break
        end
    
        if(optn.print>1)
            fprintf('%4i %#15.9g %4i\n',i,rms(r2),j)
        end
        
    elseif(optn.maxitr1==1)
        % if not iterating use predicted new residues
        r2=X*a1-y;
        j=0;
    end
    
    optn.Amat=A;
end

ssq2=r2'*r2;

% calculate output values
rx=1-ssq2/ssq0;
rx(rx<0)=0;
nreg=infolf.nreg;
pr=betainc(rx,nreg/2,(info.nres-nreg)/2,'upper');

rmse0=sqrt(ssq0/info.nres);
rmse2=sqrt(ssq2/info.nres);

info.indconv=indconv;
t=find(Acode>0);
info.sel=t(sel);
info.Acode=Acode;
info.vm=vm;
info.itr1=i;
info.nreg=nreg;
info.rmse0=rmse0;
info.rmse2=rmse2;

if(optn.print>0)
    fprintf('%4i %#15.9g %4i\n',i,rms(r2),j)
    fprintf('Test percent  %#9.3g\n',pr)
end
if(indconv~=1 && optn.maxitr1~=1)
    warning('greyboxeval  Did not converge')
end

return
end