Improved Grid_points and rand_p

[ryboselm]

This PR contains only the grid_points and rand_p related changes.

• rand_p function
• grid_points improvement for RealIndexMaps
• grid_points and rand_p tests

[ryanlevy]

Thanks Ryan!

I think the biggest change here is that rand_p changes the return value from Float to BigFloat. I have some concerns with this:

• I'm not sure we have tightened up passing Number everywhere (we could maybe just have a test that we can evaluate a basic function on a grid or something)
• grid_points doesn't match this behavior (basically everywhere just does base^L).
• I'm also not sure we want to enforce BigInt/BigFloat for say L=4. We could disbatch on having base be BigInt and force the user to deal with this, but I'm not 100% on this at the moment

For performance reasons I previously swapped index_value_to_scalar to return float(base)^L instead of base^L, but that won't work here due to InexactErrors. We could try to reconcile these choices though

@JoeyT1994 what do you think?

[JoeyT1994]

I think we should also add support for multidimensional grid points.

Like we have a function where we pass a vector of dimensions and then zip up the result of grid_points called on each of the dimensions

[ryboselm]

could you elaborate more on what your idea for this would look like?

for example, if you passed in a vector dimensions = [1,3], and suppose that grid_points(s,1) = [0, 1] and grid_points(s,3) = [0.3, 0.4, 0.5], would grid_points(s, [1, 3]) then be something like [[0, 0.3], [0, 0.4], [0, 0.5], [1, 0.3], [1, 0.4], [1, 0.5]]?

[JoeyT1994]

Yes so if for each dimension I you input you generate some vector of Ni points then it would return the iterative product which consists of every unique combination of points:

e.g. [0.25, 0.5], [0.625, 0.975, 0.9875] -> [(0.25, 0.625), (0.25, 0.975),(0.25, 0.9875), (0.5, 0.625), (0.5, 0.975), (0.5, 0.9875)]

[ryboselm]

would it be better to return an array of tuples or an array of arrays here?

I initially leaned towards having an array of arrays because evaluate takes in an array to specify the point it should evaluate the function at.

[JoeyT1994]

Thanks Ryan! I think we're getting there

[ryboselm]

For now, to deal with the BigFloat issue, I made it such that it only uses big() when L is too large to use a regular Float for. I also cleaned up a lot of the tests and simplified the code according to the suggestions.

[ryanlevy]

This assumes base 2, but we shouldnt right? However, for base 2 we can use a trick

round(ldexp(point,L))/2.0^L

That seems to work for any L that can be represented by point.

[ryboselm]

Good point, I generalized this function to take in the base as an argument

[ryanlevy]

Thanks Ryan, maybe hold off until Joey can comment, but could you specialize anything with min(63,L) etc to be a Float64 (double) instead of Number? (Or determine the max L based on the type and base). There's an assumption that the Number is Float64 and that base=2 for 'crossoverL = 63'

[JoeyT1994]

Thanks Ryan!

This looks nearly there.

I think the biggest change I would say we should do is to generate the random numbers in a way that avoids the Big issue entirely.

The easiest way to do that is to generate a random bit string of length L and then convert that to decimal. So the code would look like

function rand_point(L::Int, base::Int)
    bitstring = rand([j for j in 0:(base-1)], L)
    decimal = sum([b * (base ^ -i) for i in 1:L])
    return decimal
end

which now relies on negative powers which Julia is capable of calculating even for values of L in the 100s. This should circumnavigate any of the large L issues we are having.

[ryanlevy]

Quick comment, a better way to do the sum that reduces allocations is

decimal = sum(((i, b),) -> b * ((base * 1.0)^-i), enumerate(bitstring))"

[ryboselm]

thanks! @JoeyT1994 @ryanlevy made a quick fix