Hi Stefan, Oliver and other interested in this long email,
I just made a commit where I try to fix things to make retrievals using
grids having length 1 possible. This is very handy in some cases, but
causes problems with the interpolation inside arts. For the moment I
have solved this by special functions, and this generated quite a bit of
code.
But it could in fact be solved in a general manner. And then also fix
other things.
For running OEM the issue of concern appears when mapping the x vector
back to arts' variables. If we use 3D and a length-1 pressure retrieval
grid as example, then we end up with a situation that we want to
interpolate a Tensor3 having size X(1,nlat,nlon)
X can be a scaling factor for H2O as indicated in my ChangeLog message
In the calculation of the Jacobian, I use the grid_pos system to map
data from line-of-sights to the retrieval grids, and I once made some
function that gives a correct grid position for length-1 retrieval
grids. That is, the grid position is set to
0 0 1
It works in fact to calculate interpolation weights for this case.
On the other hand, you can not apply the interp functions, as these
functions always include idx+1, even if fd[0]=0. That is, there is a
blind summation for the values inside the grid range, even if some
values will get weight 0. My new functions could be removed if the
interp functions noticed when it is unnecessary to involve idx+1.
This would also solve:
Now it is not possible to set the grid position to the end grid point.
That is,
n-1 0 1
as idx+1 is then an invalid index. This generated extra code in the
ppath part. For basically all observations there is a ppath point
exactly at the uppermost pressure level.
Quite a lot of our interpolations could be avoided. As retrieval grids
normally are sub-sets of the forward model grids (and is the recommended
choice for numerical accuracy), there is a lot of interpolation that in
fact is not interpolation (i.e. points of new and old grid are the
same). And is this not also the case for interpolation of absorption
lookup tables? The normal case should be that the pressure grid of the
lookup table equals p_grid. As we here apply polynomial interpolation, a
lot of calculations could be saved by identifying when fd[0]=0. (or is
there special code for abs_lookup?)
The question is how to add this without making the general interpolation
slower. I assume that a boolean in the grid_pos structure could help,
that is true if fd[0] can be treated to be exactly zero. If we call this
new field nonfrac, the green 2D interpolation could like this
tia = 0;
for ( Index r=0; r<tr.nonfrac?1:2; ++r )
for ( Index c=0; c<tc.nonfrac?1:2; ++c )
{
tia += a.get(tr.idx+r,tc.idx+c) * itw.get(ir,ic,r*2+c);
}
}
(cf. interpolation.cc line 2580. Note that I had to add some temporary
asserts to make sure that my special code is OK.)
This has the drawback that r*2+c is more costly than ++iti (as in the
present solution), but we could still gain on average. Another option
would be that a similar thing is done when calculating the interpolation
weights, i.e. the number of weights is smaller when nonfrac is true. The
weights should be ordered in such way that ++iti will work. And we
should then really gain on average!?
And now realize that this bool would have saved me a lot of headache in
the ppath_step part. There are a lot of checks if fd[0] is effectively
zero (which it is for the end points of each step). Introducing the flag
would allow me to clean up this.
What do you say? Comments? Or an even better solution? If you don't
follow, ask or let's discuss next Friday.
Bye,
Patrick
_______________________________________________
arts_dev.mi mailing list
arts_dev.mi@lists.uni-hamburg.de
https://mailman.rrz.uni-hamburg.de/mailman/listinfo/arts_dev.mi
