test_levmar

(tmp_for_tar/levmar/test_levmar.m)

Function Synopsis

[v,dv,d2v] = d2gg(u,x)

Help text

 Test whether levmar_min() functions correctly, 

 first with a simple quadratic programming problem, (minimize 'ff'
 declared below).

 Then with some tests from ftp://ftp.mathworks.com/pub/contrib/, but
 the list of tests is commented out in the present script. 

 Last modified: March 2000

Listing of function file tmp_for_tar/levmar/test_levmar.m

## Test whether levmar_min() functions correctly, 
##
## first with a simple quadratic programming problem, (minimize 'ff'
## declared below).
##
## Then with some tests from ftp://ftp.mathworks.com/pub/contrib/, but
## the list of tests is commented out in the present script. 
##
##

## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>
## Last modified: March 2000

1 ;

P = 10+floor(30*rand(1)) ;	# Nparams
R = P+floor(30*rand(1)) ;	# Nobses
noise = 0 ;
global obsmat ;
obsmat = randn(R,P) ;
global truep ;
truep = randn(P,1) ;
xinit = randn(P,1) ;

global obses ;
obses = obsmat*truep ;
if noise, obses = adnois(obses,noise); end


function v = ff(u,x)
  global obsmat ;
  v = msq( u - obsmat*x ) ;
endfunction


function [v,dv,d2v] = d2ff(u,x)
  global obsmat ;
  er = -u + obsmat*x ;
  dv = er'*obsmat ;
  v = msq( er ) ;
  d2v = pinv( obsmat'*obsmat ) ;
endfunction



printf("gonna do : levmar_min\n");
mytic() ;
[xlev,vlev,nlev] = levmar_min(obses,"ff","d2ff",xinit) ;
tlev = mytic() ;

printf("Levmar should find in one iteration + one more to check\n");
printf(["LEVMAR : niter=%4d  nobs=%4d,nparams=%4d\n",...
	"  time=%8.3g errx=%8.3g   minv=%8.3g\n"],...
       nlev,R,P,tlev,max(abs( xlev-truep )),vlev);


## error("Work ahead")

global probm;
global probn;
global probnum ;
				# What follows requires some test
				# functions

for probnum = [],		# [1,11,17,19],
  printf("Gonna try to solve : \n");
  problems(probnum)

  [probn,probm,x0] = initf(probnum) ;
  ## x0 = [-1;1];
  u = zeros(probm,1) ;
  
  function v=gg(u,x)
    global probm;
    global probn;
    global probnum ;
    vv = vecfcn(probn,probm,x,probnum) ;
    
    v = msq( u(:) - vv(:) ) ;
  end

  function [v,dv,d2v] = d2gg(u,x)
    global probm;
    global probn;
    global probnum ;
    [vv,jj] = func(probn,probm,x,probnum,3);
    er = -u(:) + vv(:) ;
    dv = er'*jj ;
    v = msq( er ) ;
    d2v = jj'*jj ;
    d2v = pinv( d2v ) ;
    ## d2v = pinv( d2v, max(eig(d2v))/10000 ) ;
  endfunction
  
  printf("gonna do : levmar_min\n");
  mytic() ;
  [xlev,vlev,nlev] = levmar_min(u,"gg","d2gg",x0) ;
  tlev = mytic() ;
  
  
  printf(["LEVMAR : niter=%4d  nobs=%4d,nparams=%4d\n",...
	  "  time=%8.3g  minv=%8.3g\n"],...
	 nlev,probm,probn,tlev,vlev);
  xlev'
  printf("Check for yourself\n");
  
  eval(['help ',prob2name(probnum)]);
end

if 0,
  prob2num = name2probnum("helix");
  [probn,probm,x0] = initf(probnum) ;
  yy=xx=-4:0.3:4;n=size(xx,2);zz=zeros(n);
  for i=1:n,for j=1:n,zz(i,j)=sum(helix(probn,probm,[xx(i),yy(j)],1).^2);end;end;
  mesh(xx,yy,zz)
end