Dear Patrick, an interesting mail!
Oliver and me both pondered your questions, and now had a brief discussion. We
think this kind of change is possible in principle.
Suggestions:
1. Hard to do this automatically, and perhaps dangerous. Perhaps we could add
the threshhold for considering a new point to be exactly on an old grid point
as a new argument to gridpos? With a default of zero.
2. The awkward and potentially inefficient “tr.nonfrac?1:2” expressions in the
loop indices could be avoided by not using a bool like nonfrac, but instead
simply the loop dimension, 1 or 2. The standard value of 2 would give the old
behaviour, the value 1 the new “nonfrac” case. Your piece of code would be:
tia = 0;
for ( Index r=0; r<tr.loop_dim; ++r )
for ( Index c=0; c<tc.loop_dim; ++c )
{
tia += a.get(tr.idx+r,tc.idx+c) * itw.get(ir,ic,r*2+c);
}
}
If you need/want a boolean for some external use cases, one could simply add a
member function nonfrac() that returns true for loop_dim==1.
3. Yes, the “r*2+c” expression in the itw.get command is awkward, and will
waste time for the high-dimensional cases. Perhaps a better solution, that you
already hint at, is to simply modify the interpweights functions consistently.
They have similar nested for loops as in the example above. If the new loop_dim
variables are used also there, then the “ tia +=“ expression in the kernel of
the nested loops can stay exactly as it is now. Only the needed weights will be
generated in the first place. (But this may screw up code that does things with
the weights directly, I don’t know if you do this anywhere. One should do a
systematic search for places where the weights are accessed outside the
interpolation routines themselves.)
As far as we can see, with these modifications there will be no performance
loss. As you point out, it would even be faster for cases where old and new
grid are exactly the same. (Although I am not so optimistic that this will be a
big gain in practice for most users.)
Best wishes,
Stefan
> On 15. Mar 2019, at 08:26, Patrick Eriksson <patrick.eriks...@chalmers.se>
> wrote:
>
> 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
smime.p7s
Description: S/MIME cryptographic signature
_______________________________________________ arts_dev.mi mailing list arts_dev.mi@lists.uni-hamburg.de https://mailman.rrz.uni-hamburg.de/mailman/listinfo/arts_dev.mi
