[SciPy-User] scipy.optimize.leastsq question

[Ralph Kube]

Hello people,
I am having a problem using the leastsq routine. My goal is to
determine three parameters r_i, r_s and ppw so that the residuals
to a model function a(r_i, r_s, ppw) to a measurement are minimal.
When I call the leastsq routine with a good guess of starting values, it
iterates 6 times without changing the vales of the initial parameters 
and then exits without an error.
The function a is very complicated and expensive to evaluate. Some
evaluation is done by using the subprocess module of python. Can this 
pose a problem for the leastsq routine?


This is in the main routine:

import numpy as N

for t_idx, t in enumerate(time_var):

	r_i = 300.
	r_s = 1.0
	ppw=1e-6
	sza = 70.
	wl = N.arange(300., 3001., 1.)		

	albedo_true = compute_albedo(r_i, r_s, ppw, sza, wl)
	# This emulates the measurement data
	albedo_meas = albedo_true + 0.01*N.random.randn(len(wl))
		
	print 'Optimizing albedo'
	p0 = [2.*r_i, 1.4*r_s, 4.*ppw]
	plsq2 = leastsq(albedo_residual, p0, args=(albedo_meas, sza,
wl))
	print '... done: ', plsq2[0][0], plsq2[0][1], plsq2[0][2]
	albedo_model = compute_albedo(plsq2[0][0], plsq2[0][1], plsq2[0][2], 
sza, wl)

The residual function:
def albedo_residual(p, y, sza, wvl):
	r_i, r_s, ppw = p
	albedo = compute_albedo(r_i, r_s, ppw, sza, wvl)
	err = albedo - y
	print 'Albedo for  r_i = %4.0f, r_s = %4.2f, ppw = %3.2e \
		Residual squared: %5f' % (r_i, r_s, ppw, N.sum(err**2))

	return err

The definition of the function a(r_i, r_s, ppw)
def compute_albedo(radius_ice, radius_soot, ppw, sza, wvl):

The output is:
Optimizing albedo
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared: 
0.973819
... done:  600.0 1.4 4e-06

To check for errors, I implemented the example code from 
http://www.tau.ac.il/~kineret/amit/scipy_tutorial/ in my code and it 
runs successfully.

I would be glad for any suggestion.


Cheers, Ralph