John, nice work.
I encountered a minor glitch.
When specifying the modelterms as an array of exponents, the variable 'varlist' was not assigned a value.
I fixed this by inserting the following at line 160:
else varlist = repmat({''},1,p);

I had a similar issue. I group multiple test cases in a file and place all of the test files in one directory. Then I call runtests with the directoryname as the argument.
I also noticed that if one of the test cases in one file has a syntax error, the whole file will be skipped silently.

To mend this, I made a small change in 'fromName' method in 'TestSuite.m'.

Original:
>>>
else

try
if nargout(name) == 0
<<<

modified:
>>>
else
try
nout = nargout(name);
catch
warning('''%s'' does not appear to be a valid function.',...
name)
end
try
if nout == 0
<<<

If the file has a syntax error, the nargout test will fail and a warning is printed before the actual testing starts.

Whatever value of z0 that seems to work should be fine there.
You can think of it as a truncated Taylor expansion if you like, which it essentially is. I simply chose a straight line that matched the slope and value of the log function at the break point, to create a continuous, differentiable function that will extrapolate linearly. So sort of a spline too.

Thank you, John; your suggestion worked. I chose Z0 to be 0.5% of max(Z) of 2E12 in the data I'm fitting. Now, using gridfit to fit the conditionally log transformed data yields a positive surface which does not overshoot.
A question: how did you arrive at the second part of the transform, i.e. ((z - z0)/z0 + log(z0)).*(z>z0)? It works but it's not immediately obvious to me why. (Some sort of Taylor expansion of the log?)
Thanks,
Neil

Neil,
Of course, my response should be - picky, picky, picky!
This is a consequence of a log transform I guess. You could use the 'springs' method as a regularizer to avoid that, but it might be too extreme.
So an alternative is to use a piecewise transformation. Thus, choose some breakpoint, above which the transformation will be linear, and below which, it will behave as a log transform. I'll call that value z0. So we can write the transform function W(z) as
W = @(z,z0) log(z).*(z<=z0) + ((z - z0)/z0 + log(z0)).*(z>z0);
You can plot it, picking perhaps z0=3 as a breakpoint for an example.
ezplot(@(z) W(z,3))
Of course, the transformation is invertible. This should work:
Winv = @(w,z0) exp(w).*(w<=log(z0)) + z0*(w - log(z0) + 1).*(w>log(z0));
Again, we can plot it, here for z0 = 3 again, as:
ezplot(@(w) Winv(w,3))
Something like that should give you the best of all worlds, and all you need to do is choose some value for z0 that makes sense to you. I might start with the mean of z as a good choice for z0.

Hi John. Thanks for the tip regarding using a log transform. It seems to be a sensible approach and I tried it. The fitted surface is now all positive, as expected. However, it also overshoots: the highest peak (Z) in the fitted surface is now 3 times the largest data point. Previously the fit closely followed the data.
I should add that the data did have 0's in it and I didn't want to completely remove them because they represent boundary conditions for the somewhat sparse data I am fitting to; for example, the spectrum vanishes at a certain energy (e.g. E=0) or at some angles. Instead, as you suggested, I replaced the 0's with a small number; I tried 1E-10, 1E-5, 1E-1, and 1. The overshooting appeared in all cases and is probably unrelated to the presence of 0's or small numbers in the data.
I tried to attach pictures of the fitted surface before and after the log transform, but there doesn't seem to be a way to do that on this forum.
Thanks,
Neil

Comment only