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.


arts_dev.mi mailing list

Reply via email to