Conversation
Previously there was a singularity at r==0.
Code that was in grid.c and weights.c has been put in a new module, random_points.c.
This is in preparation for the introduction of a new grid-point-selection algorithm.
The module weights.c is no longer compiled, since it is now obsolete.
Additional changes connected with the above:
- Added a new (optional) user parameter 'samplingAlgorithm'.
- Added external variables modelRadiusSquared, densityNormalizer,
numCollisionPartners.
- Added macro DENSITY_POWER.
- New struct locusType.
The new algorithm calls the function randomsViaTree() recursively in a 2^D-tree divide-and-conquer approach. The aim of this is to allow more accurate representation of models for which the density varies strongly across the model, without the severe time penalties that straightforward use of the rejection method would entail. (The previous algorithm uses rejection method, but it 'flattens' the density function in order to avoid long run times, which has the effect that high-density regions are undersampled.) The algorithm assesses the density gradients over a sub-volume of the model space, and if those gradients are judged too large for quick results from the rejection method, the sub-volume is divided into 8 smaller sub-volumes and the same function is called for each one of these. The division is performed by cutting the volume approximately in halves (at the 0.5 point plus some random dither) along each of the (3) axes of the coordinate space. Additional changes connected with the above: - New macros: N_TREE_RANDOMS MAX_RECURSION MAX_N_TRIALS_TREE TREE_DITHER
The new grid-selection algorithm estimates the density gradients within a sub-volume. This estimate may be too low if the function has narrow peaks. In this case, the functioning of the algorithm is improved if the user pre-supplies a list of the maximal values of the density function. Additional changes connected with the above: - New (optional) user parameters numDensityMaxima, densityMaxLoc, densityMaxValue. Note that I have left the appropriate lines commented out for now in the example model.c file. - New macro: MAX_N_HIGH
In general it seems desirable to have the grid point number density follow the total density of the collision partners, as is now more accurately provided by the new algorithm. If the density function is heavily peaked, whereas the chosen number of grid points is kept relatively small, a situation could arise in which the distribution of points in the low-density regions of the model was very sparse. To prevent this occurring, the user can now set a parameter minPointNumDensity. The units of this should if possible save the user calculation, and to be scale-invariant: the chosen unit therefore is points per spherical model volume. Note that this feature is only implemented when par->samplingAlgorithm==1. Additional changes connected with the above: - A new extern variable minDensity is defined.
Previously there was a singularity at r==0.
Code that was in grid.c and weights.c has been put in a new module, random_points.c.
This is in preparation for the introduction of a new grid-point-selection algorithm.
The module weights.c is no longer compiled, since it is now obsolete.
Additional changes connected with the above:
- Added a new (optional) user parameter 'samplingAlgorithm'.
- Added external variables modelRadiusSquared, densityNormalizer,
numCollisionPartners.
- Added macro DENSITY_POWER.
- New struct locusType.
The new algorithm calls the function randomsViaTree() recursively in a 2^D-tree divide-and-conquer approach. The aim of this is to allow more accurate representation of models for which the density varies strongly across the model, without the severe time penalties that straightforward use of the rejection method would entail. (The previous algorithm uses rejection method, but it 'flattens' the density function in order to avoid long run times, which has the effect that high-density regions are undersampled.) The algorithm assesses the density gradients over a sub-volume of the model space, and if those gradients are judged too large for quick results from the rejection method, the sub-volume is divided into 8 smaller sub-volumes and the same function is called for each one of these. The division is performed by cutting the volume approximately in halves (at the 0.5 point plus some random dither) along each of the (3) axes of the coordinate space. Additional changes connected with the above: - New macros: N_TREE_RANDOMS MAX_RECURSION MAX_N_TRIALS_TREE TREE_DITHER
The new grid-selection algorithm estimates the density gradients within a sub-volume. This estimate may be too low if the function has narrow peaks. In this case, the functioning of the algorithm is improved if the user pre-supplies a list of the maximal values of the density function. Additional changes connected with the above: - New (optional) user parameters numDensityMaxima, densityMaxLoc, densityMaxValue. Note that I have left the appropriate lines commented out for now in the example model.c file. - New macro: MAX_N_HIGH
In general it seems desirable to have the grid point number density follow the total density of the collision partners, as is now more accurately provided by the new algorithm. If the density function is heavily peaked, whereas the chosen number of grid points is kept relatively small, a situation could arise in which the distribution of points in the low-density regions of the model was very sparse. To prevent this occurring, the user can now set a parameter minPointNumDensity. The units of this should if possible save the user calculation, and to be scale-invariant: the chosen unit therefore is points per spherical model volume. Note that this feature is only implemented when par->samplingAlgorithm==1. Additional changes connected with the above: - A new extern variable minDensity is defined.
The 'init_lte' parameter was gone.
- Module random_points.c has been deleted and replaced by tree_random.c with header file tree_random.h. The algorithm has been changed in three main ways: * Construction of the sub-cell tree is now done separately to filling the sub-cells with points. * The number of points assigned to each sub-cell is calculated via the multinomial distribution, instead of (incorrectly) via a kind of progressive binomial distribution. * The set of points which are used as input to all the sub-cell rejection sessions is now chosen via a quasi-random routine which spaces them much more evenly apart. This renders post-grid smoothing unnecessary and thus allows much sharper changes in grid point density distribution. - There is a new macro TREE_POWER by which the density function within gridDensity should be raised in order to approximately match the grid density achieved by the previous flattened-density function with par->sample==2 (logarithmic radius sampling). Why is this a good thing? Because the old routine only gave a good match for the example model, which has a single central condensation. A model with concentrations away from the central core would *not* be handled well by this corner-cutting, but *would* be handled well by the new routine.
|
I've merged the current master branch into this PR and, if there are no objections in the next day or so, I will merge it into master (i.e. accept the PR). A summary of the changes in this PR is as follows. The PR addresses two separate challenges:
A: The canonical algorithm for choosing randoms to follow some generic distribution is called the rejection method (RM). It is robust but can be very slow if the desired distribution function has large variations. LIME<=1.6 basically uses RM but applies a number of tweaks and biases designed to speed it up. These work well for the example model, but in fact only because they are tuned to suit the example model - for a different model they will be about as quick but will not necessarily produce anything close to an accurate sampling of the model. The present PR applies pure RM but essentially on a local level, which allows it to be fast. It does this by dividing the model space into 8 sub-cubes, then dividing each of these again, and so forth, until the variation of the distribution function within a sub-cube falls to an acceptable level. B: LIME<=1.6 achieves this by smoothing the points after the initial choice. There are two problems with this: it is an iterative process, and fairly slow; and it blurs out sharp boundaries in the point density. The present PR chooses points, not randomly, but sub-randomly, using a Halton-sequence algorithm, producing points which 'look' random but which avoid approaching too close to each other. Some careful attention was necessary to avoid artefacts at the sub-cube boundaries but test distributions of points from the total algorithm (8-tree recursive RM plus Halton-sequence choice of points in each sub-cube) ended up looking very convincing. I added a macro TREE_POWER to lime.h whose value I have somewhat arbitrarily set at 2.1. The reason for this is as follows. LIME<=1.6 (at least with |
|
Note that the Halton-sequence thing removes the need to smooth after grid point choice. |
3 participants
This is a rebased version of @imckstewart 's PR #80.