dhmunro
commented
11 years ago Actually, the bug is in your variance_kernel function.
Yorick vector functions, like the builtin exp(x), cos(x), etc., must return values with the same dimensions as their inputs. However, when presented with a scalar value, variance_kernel returns an array of length 1. For example variance_kernel(.001) returns [1.10637e-14], rather than simply 1.10637e-14. Furthermore, variance_kernel always returns only a 1D array, even if its argument is multi-dimensional, again breaking yorick's vector syntax rules. This is a little picky, I admit, but it will cause failures in a large number of places in yorick, so applying the band-aid you suggest in the romberg routine will not address the real substance of your complaint.
You also failed to vectorize your function - the for loop will slow yorick down by a factor of several hundred if you present variance_kernel with a big array. It is slightly tricky to fix all this without falling prey to the same problem of mismatching input and output dimensionality. The trick is to use yorick's ".." index function, which works properly for zero dimensions as well as for multiple dimensions. Here is a correct implementation of variance_kernel, for which the output has the same dimensionality as the input, and without the for loop so it will take full advantage of yorick's vector syntax:
func variance_kernel(x)
{
// unused??? z = [25.,275.,425.,1250.,4000.,8000.,13000.];
Cn2 = [0.4393, 0.3054, 0.2328, 1.2228, 0.7947, 0.2381, 0.2627]*1e-13 ;
R = 4. ; // telescope diameter
L0 = 60.; // outer scale (m)
x = x(-,..); // prepend leading index matching Cn2 index
fZ = (4.*bessj(2,2.*pi*R*x)/(2.*pi*R*x))^2 ;
fA = 1.;
phi = 0.0097 * Cn2 * (x^2 + 1./L0^2)^(-11./6.);
return (2.*pi * x*phi*fZ*fA)(sum,..); // sum over Cn2 index, leaving original x dimensions
}This has the advantage of making the formulas a little clearer, more like you would write them algebraically, although just as in algebra, you have to learn to read implicit multiple indices as you read implicit multiple function arguments in an abbreviated algebraic expression.
Incidentally, "float" is almost never what you want. Yorick's default real and integer types are "double" and "long", and you should stick with those unless there is a very good reason not to do so.
mvandam
The error message I get is: ERROR (romberg) non-scalar operand to BranchFalse LINE: 39 FILE: /usr/lib/yorick/i/romberg.i
I fixed it locally by changing line 39 to: if (abs(dss) <= (epsilon*abs(ss))(1)) return ss;
See code below...
func variance_kernel(x) { z = [25.,275.,425.,1250.,4000.,8000.,13000.]; Cn2 = [0.4393, 0.3054, 0.2328, 1.2228, 0.7947, 0.2381, 0.2627]*1e-13 ; R = 4. ; // telescope diameter L0 = 60.; // outer scale (m)
nx = numberof(x); output = array(float,nx);
fZ = (4._bessj(2,2._pi_R_x)/(2._pi_R*x))^2 ; fA = 1.;
for (c1=1;c1<=nx;c1++){
phi = 0.0097Cn2(x(c1)^2. + 1./L0^2)^(-11./6.); output(c1) = sum(2_pi_x(c1)_phi_fZ(c1)*fA); }
return output; }
tt_var = romberg(variance_kernel,1e-6,1e-2);