Added optional user function gridDensity

[allegroLeiden]

This is intended to specify the probability distribution function for the internal (i.e.
non-sink) grid points. If the user does not supply their own version of this function,
it defaults to the algorithm which was previously used.

If this is accepted, I'll re-write the treegrid PR #80 (now #93) to incorporate it, which will avoid the necessity for a lot more user parameters.

This PR is relevant to issues discussed in #62 and will also overlap with or affect #71.

[smaret]

Could you please document the gridDensity function in the the user manual?

[smaret]

@AttilaJuhasz : Your PR #71 has a similar purpose than this one. Could you please have a look at what Ian proposes here and let us know what do you think?

[smaret]

Looks good, but LIME seems to hang if fracDensity is set to values strictly lower than 1:

void
gridDensity(inputPars par, double x, double y, double z, double *fracDensity){
  *fracDensity = .5; // this hangs
  /* fracDensity = 1.0; // this works */
}

[smaret]

Also has par->samplingany effect if one sets fracDensity ?
IMS: yes (same as with standard Lime). See below.

[allegroLeiden]

I'll look into these. Jeez you guys, it is either nothing or a shitload with you. ;)

[allegroLeiden]

Ok, I had a look.

The root cause of the 'hang' is the slightly weird algorithm used for selecting grid points. Now that I look closer at it, it is not a straight rejection algorithm as I had assumed (albeit to a flattened density function), but a modification of it with some slightly strange results. In the iteration loop from lines 541 to 574 of grid.c, random point locations are chosen, with a weighting controlled by par->sampling. The weighting alone takes it away from pure rejection, but the real kicker is that the random number which controls whether a point is accepted or rejected is chosen outside this loop. A pure implementation of the rejection method would have the random chosen inside the loop and would not weight the pre-selection of points.

Consider an example where the random value (the variable temp, line 538) is 0.8. If there is no point in the model where fracDensity is greater than this, the loop will never terminate.

Possible solutions I can think of:

1. Change things such that the returned value of fracDensity is (as is evilly done in standard Lime) quasi-normalized by dividing it by some guessed maximum value.
2. Move to the treegrid algorithm, which gives a much closer approximation to true rejection method (with the accompanying fidelity to the desired distribution function) while avoiding its slowness.
3. Move the choice of temp inside the loop, with unknown effects on the point distribution and algorithm speed.
4. Add a warning to the documentation advising the user to choose a gridDensity() such that there is some place in the model where fracDensity>=1.

I favour 4 immediately and 2 eventually. I could also add another layer of loop outside the choice of temp, plus a maximum trial number on the inner loop. Pseudo code as follows:

for each point{
  while(1){
    temp=gsl_rng_uniform(ran);
    converged=0;
    i=0;
    while(!converged && i<max_number_tries){
      choose point location r;
      converged=(temp<fracDensity(r));
      i++;
    }
  }
}

[smaret]

Thanks for looking into this.

Add a warning to the documentation advising the user to choose a gridDensity() such that there is some place in the model where fracDensity>=1.

This sounds good to me. In fact I thought that fracDensity was some kind of probability (hence comprised between 0 and 1) for a grid point to be chosen at a specific place. I think that a more detailed description of the function and its return value in the documentation would be useful.

[allegroLeiden]

In fact I thought that fracDensity was some kind of probability (hence comprised between 0 and 1) for a grid point to be chosen at a specific place.

That is exactly what I intended, and that is exactly the effect it will have, as soon as we can alter the other parts of the algorithm to become a more pure rejection method. Without doing the latter, there is no present way to achieve this.