[X,y,r,info]=greyboxcore(condfn,resfn,data,vndata,optn)

function [X,y,r,info]=greyboxcore(condfn,resfn,data,vndata,optn)
% greyboxcore  Greybox evaluate possibility of improving model
%  2012-12-31  Matlab8  Copyright (c) 2013, W J Whiten  BSD License
%
% [X,y,info]=greyboxcore(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 conditions vector or matrix
%                     for data set idata
%             pscale Vector of scale (ie typical nonzero size) values 
%                     for parameters (optional if optn.Amat given)
%                     if optn.Amat not given, initial parameter values
%             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 (a matrix if c is a matrix)
%             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 A matrix elements
%                   0 do not use, 1 always use, 2 select
%           .mres  Estimate of total number of residuals
%
%  X       Regression matrix
%  y       Regresion right hand side
%  r       Model residues
%  info    Struct giving:
%           .Amat  A matrix value from input
%           .Acode Matrix code for terms in A matrix (>0 use term)
%           .ncond Number of condition terms
%           .npar  Number of parameters
%           .nres  Total number of residuals
%
% includes c & p as matrices

if(nargin<5)
    optn=struct();
end

nresfn=nargout(resfn);
ncondfn=nargout(condfn);

% default option values, including those not used
optn=optndfts(optn,'Amat',[],'Acode',[],'nres',1000);
A=optn.Amat;
Acode=optn.Acode;
mr=optn.nres;
mrr=mr;

if(ncondfn==1 && isempty(A))
    error('greybox:  condfn must give pscale if optn.Amat not given')
end

% calculate regression matrix
nr=0;
nrr=0;
ssq1=0;
sqreps=sqrt(eps);
for ii=1:length(vndata)
    i=vndata(ii);
    
    % get condition and parameter vakues
    [p,pscale,pbase,c]=parcalc(condfn,data,i,A);
    nc2=size(c,2);
    
    % initial values now sizes are known
    if(ii==1)
        nc=size(c,1);
        np=size(p,1);
        if(isempty(Acode))
            Acode=ones(np,nc);
        end
        nx=sum(Acode>0);
        mx=sum(nx);
        Acode1=Acode>0;
        if(isempty(A))
            Ainit=zeros(np,nc);
        else
            Ainit=A;
            Ainit(Acode1)=0;
        end
        r=zeros(mr,1);
        y=zeros(mr,1);
        X=zeros(mr,mx);
    end
    
    % calculate residual and its derivative
    if(nresfn==2)
        [rr,dr1]=resfn(data,i,p);
        rr=rr(:);
        ir=length(rr);
    else
        rr=resfn(data,i,p);
        rr=rr(:);
        ir=length(rr);
        dr=zeros(ir,np);
    end
    
    for k=1:nc2
        if(nresfn==2)
            dr=dr1(:,:,k);
        else
            for j=1:np;
                pj=p(j,k);
                if(abs(pj)>1e-150 || pj==0)
                    pscj=pscale(j,k);
                    dp=pscj*sqreps;
                    p(j,k)=pj+dp;
                    r1=resfn(data,i,p);
                    dr(:,j)=(r1(:)-rr)/dp;
                    p(j,k)=pj;
                else
                    dr(:,j)=0;
                end
            end
        end
        ssq1=ssq1+rr'*rr;

        % include into regression matrices
        nr1=nr+1;
        nr=nr+ir;
        X1=zeros(ir,mx);
        ix=0;
        for j=1:nc
            ix1=ix+1;
            ix=ix+nx(j);
            X1(:,ix1:ix)=c(j,k)*dr(:,Acode1(:,j));
        end

        if(nr>mr)
            mr=nr+round(1.5*nr);
            y(mr)=0;
            X(mr,1)=0;
        end
        y(nr1:nr,1)=dr*(p(:,k)-pbase-Ainit*c(:,k))-rr;
        X(nr1:nr,:)=X1; 
    end
    nrr1=nrr+1;
    nrr=nrr+ir;
    if(nrr>mrr)
        mrr=nrr+round(1.5*nrr);
        r(mrr)=0;
    end
    r(nrr1:nrr,1)=rr;
end

if(nr<mr)
    y=y(1:nr);
    X=X(1:nr,:);
end
if(nrr<mrr)
    r=r(1:nrr);
end

info=optn;
info.ncond=nc;
info.npar=np;
info.nres=nrr;

return
end