function [f,p,var,iter,cor,std]=leasqr1(x,y,pin,func,stol,niter,wt)
% [f,p,var,iter,cor,std]=leasqr1(x,y,pin,func,{stol,niter,wt})
%  {}= optional parameters,can use [] to hold place at default value
% Levenberg-Marquardt nonlinear regression of y=f(x,p) , where:
% ALGORITHM FROM OPTIMISATION TOOLBOX (LEASTSQF)
% x=vec or matrix of indep variables, 1 row/observation: x(i)=[x1(i)..xn(i)]
% y=vec of observed values, same no. of rows as x.
% e.g. y=p0 +p1*x1 +p2*x1^2 +p3*x2 +p4/x2 ==> y=f([x1 x2],[p0..p4])
% wt=vec(dim=1 or length(y))of stat. wts:(sqrts of var(y))^-1; default=1
% pin=vector of initial parameters to be adjusted by leasqr.
% func=name of function in quotes,of the form y=f(x,p) y MUST RETURN AS COL
% stol=scalar tolerances on fractional improvement in ss,p,default stol=.0001
% niter=scalar max no. of iterations=Function evaluations.default=200*length(p)
%	    OUTPUT VARIABLES
% f=vec function values computed in function func.
% p=vec trial or final parameters. i.e, the solution.
% iter=scalar no. of iterations(Function evaluations) used.
% var = variance of residuals =ss/(m-n)
% cor= correlation matrix for parameters
% std=standard deviation error for each parameter


%	      A.Jutan (519)-679-2111	 (apr 91/feb94)

%Arthur Jutan 
%Dept of Chemical and Biochemical Enginering
%Univ of Western Ontario
%London Ontario Canada, N6A 5B9

%     uses modified LEASTSQ-->LEASTSQF from OPTIM toolbox
%set default if args not given or if arg given as []
% set remaining params to []
if nargin ==4,stol=[];niter=[];wt=[]; end
if nargin ==5,niter=[];wt=[]; end
if nargin ==6,wt=[]; end
%
if (nargin <=6|isempty(wt)) ,wt=1.0;end; %LS WEIGHT =1.0 DEFAULT
if (nargin <=5|isempty(niter)) ,niter=0 ;end; %default
if (nargin <=4|isempty(stol)) ,stol=0;end;   %default
%
y=y(:);wt=wt(:);pin=pin(:); %change all vectors to columns
m=length(y); n=length(pin); p=pin;[m1,m2]=size(x);
% check data vectors- same length?
if m1~=m ,error('input(x)/output(y) data must have same # of rows ') ,end;
f=eval([func,'(x,p)']);
[lf,wf]=size(f);
if lf~=m,error('output from func must be column same length as y data'),end;
%
% USE LATEST LEVENBERG -MARQUARDT METHOD   (OPT TOOLBOX)
options(2)=stol; options(3)=stol;
options(14)=niter; % Max # of function evals
options(1)=1; % DISPLAY PROGRESS FROM LEASTSQF **************

% Weighted LS (default wt=1.0)

[f,p,jac,options]=leastsqf(x,y,pin,wt,func,options);
mean(jac');
iter=options(10); % function eval count
e=f-y; %error
ss=e'*e; % best sum of squares
if (m~=n),var=ss./(m-n);else, var=NaN;end % variance of Residuals
% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
% re-evaluate the Jacobian at optimal values
     %jac=feval(dfdp,x,f,p,dp,func);
% Note Jacobian Scaled by wt  in leastsqf
%	That is true Jacobian(no weighting)=jac(from leastsqf)/wt
	jac=jac*diag(wt)^(-1); % calc true (unweighted)Jacobian
	    xtx=jac*diag(wt)^(-1)*jac';% to allow for non equal wts(D&Sp109)
%	Error Variance = diag(wt)*sigma^2
      % del rows and cols corresponding to zero variance
      % occurs when some elemnts of dp set to zero (i.e. dont adjust)
	% does not apply in this version since dfdp.m not used(aj '94)
% calc reduced size of xtxinv
      i=find(diag(xtx)==0);
      xtx(:,i)=[];
      xtx(i,:)=[];
      xtxinv=inv(xtx);
      
      %calc correlation matrix cor
	d1 = diag(xtxinv);
	cor = xtxinv./sqrt(d1*d1');

% Plot the fit vs data
	
      t=1:m;
      plot(t,y,'o',t,f,'g')
      title(' y(data) and f(model) versus observation number')
      xlabel(' observation number');
      ylabel(' data and model')
      grid;
      
      disp(' Least Squares Estimates of Parameters')
      disp(p')
      disp(' correlation matrix for parameters ')
      disp(cor)
      vcv=xtxinv.*var; % mult by var of residuals~=pure error
%Formulae for vcv=(x'.vo.x)^-1 *sigma^2 where meas error Var, v=[vo]*sigma^2	
      std=sqrt(diag(vcv)); % calc std error for each param
      disp('2*standard deviation (95%CL) for each parameter')
      disp(2*std')
%Calculate R^2 (Ref Draper & Smith p.46)
      r=corrcoef(y,f);
      r2=r(1,2).^2;
      disp('Variance of Residuals  ' )
      disp(  var )
      disp( 'Correlation Coefficient R^2')
      disp(r2)
