Updating Inverse Transform Sampling class

[mohsensadr]

Allow updating bounds and local number of particles.

[mohsensadr]

Do we have a test for this function?

I added a test where we can check min/max of random variables before and after changing the boundaries. test/random/TestInverseTransformSamplingUpdateBounds.cpp

I use the same density and ITS to sample position and velocity. For position inside [-2,2]^2, and velocity inside [0,2]^2. The program writes

dim    min(pos)    max(pos)    min(vel)    max(vel)
0 -1.9796144197e+00 1.9890737993e+00 1.7195022872e-05 1.9954672006e+00
1 -1.9997797847e+00 1.9999229952e+00 3.6506341035e-06 1.9999535412e+00

[mohsensadr]

@matt-frey The current way of computing min/max of cdf for a bounded domain uses layout to find the bounds on each rank.

const typename RegionLayout::host_mirror_type regions = rlayout_r.gethLocalRegions();

                umin_m[d] = dist_m.getCdf(regions(rank)[d].min(),d);
                umax_m[d] = dist_m.getCdf(regions(rank)[d].max(),d);

https://github.com/IPPL-framework/ippl/blob/4510eedee0c32f9ed3467fbda2f3116d8f36ea6d/src/Random/InverseTransformSampling.h#L70C9-L97C11

What if one wants to use inverse transform sampling without any domain decomposition? Shouldn't we have another implementation where the user only specifies the bounds of the domain, and not rlayout...

[matt-frey]

@mohsensadr Instead of passing the RegionLayout, we could change to just pass the domain bounds. I don't know if RegionLayout is explicitly needed. @srikrrish may know more what is possible because Sri implemented the inverse sampling.

[srikrrish]

@mohsensadr Instead of passing the RegionLayout, we could change to just pass the domain bounds. I don't know if RegionLayout is explicitly needed. @srikrrish may know more what is possible because Sri implemented the inverse sampling.

I agree.. I was thinking the same during yesterday's meeting. Apart from the bounds I don't think RegionLayout is needed anywhere else.

[mohsensadr]

@matt-frey @srikrrish We essentially need bounds (min/max) of sub-domain (called local regions) corresponding to each rank. So, we can have another constructor that takes in rmin/rmax (bounds of the whole domain) as well as local_min/local_max (bounds of local regions on each rank).

By the way, is there a way to get the local bounds without rlayout?

           const typename RegionLayout::host_mirror_type regions = rlayout_r.gethLocalRegions();
            
            int rank = ippl::Comm->rank();
            Vector<T, Dim> nr_m, dr_m;
            for (unsigned d = 0; d < Dim; ++d) {
                nr_m[d] = dist_m.getCdf(regions(rank)[d].max(), d) - dist_m.getCdf(regions(rank)[d].min(),d);
                dr_m[d]   = dist_m.getCdf(rmax_m[d],d) - dist_m.getCdf(rmin_m[d],d);
                umin_m[d] = dist_m.getCdf(regions(rank)[d].min(),d);
                umax_m[d] = dist_m.getCdf(regions(rank)[d].max(),d);
            }

[mohsensadr]

@srikrrish I added another constructor for inverse transform sampling that takes local bounds for each rank, instead of rlayout. Please let me know if you have any comments.