|
GRALE
|
Apart from being able to simulate gravitational lensing scenarios, the real purpose of GRALE is lens inversion, i.e. trying to find a gravitational lens model which explains the configuration and shapes of observed images. To make this possible, a genetic algorithm is used, the basics of which are explained in these two articles:
The first article explains the genetic algorithm in case only a single criterion is optimized (a 'single-objective' genetic algorithm). The second one explains the case in which multiple criteria are optimized at the same time (a 'multi-objective' genetic algorithm).
The core of the genetic algorithm is not a part of GRALE, but was moved to a separate library so that it can be used for other problems as well. This library is called MOGAL (Multi-Objective Genetic Algorithm Library). The lensing specific part is a part of GRALE of course.
A genetic algorithm is an optimization method inspired by evolution as it happens in the world around us. It mimics the principle of 'survival of the fittest' to 'breed' solutions to a particular problem. As you will notice, many terms used are actually borrowed from the field of biology.
To look for a solution, a genetic algorithm (GA) typically works with a fixed size set of trial solutions, which are updated iteratively. Such sets are typically referred to as 'generations', in the sense that the GA tries to improve the trial solutions generation after generation. Starting from a first generation that contains randomly initialized solutions, a new generation is created by combining and mutating the trial solutions in the current generation. Often, such a trial solution is referred to as a 'genome', a 'chromosome' or something similar. This same procedure is then repeated to create another generation and another, etc.
The key ingredient to make the generations evolve towards better solutions is called 'selection pressure'. To mimic natural selection and the principle of survival of the fittest, we would like to see that better solutions generate more offspring. There are many different ways to do this, but the lens inversion method uses the way described in the mogal::GeneticAlgorithmParams documentation. Basically, the genomes in a generation are ranked according to how well they're able to explain the observed lens data, and genomes with a higher rank in general produce more offspring.
To be able to rank the genomes in a generation, we need to assign a number (or several numbers in case of a multi-objective GA) to the genome which describes how well it performs, which in this case means how well it explains the observed lensing scenario. Such a number (or set of numbers), is called the 'fitness measure', or just 'fitness'.
To be able to implement your own inversion module, you do not need to know much about the inner workings of the GA used. But there is something that's absolutely essential to understand. A trial solution, or a genome, describes the mass distribution by specifying the weights of basis functions arranged on a specific grid. It is these weights that the GA tries to optimize. But the weights on their own do not fully describe the mass distribution. The weights are always rescaled so that they lie between 0 and 1, so they only describe the shape of the mass distribution. The scale factor to really fix the mass distribution still needs to be fixed, and this is not a part of the genome.
Instead, the GA performs two steps: first it looks for the best possible mass scale weight for the genome, and then, using this weight, all fitness values (remember that it can be multi-objective) are evaluated. What the best weight is, is up to you: you need to specify a function calculateMassScaleFitness that evaluates the situation for a specific mass scale weight, and returns a value indicating how well this mass scale weight performs.
After looking for the lowest value of this mass scale weight, this then fixes the entire mass distribution, and using this mass distribution the current situation is evaluated in full. To do this, a call to calculateOverallFitness is made, which should fill in all the fitness values. But we shall encounter these functions below, the important thing to remember is that it's a two-step process, and especially for multi-objective algorithms, a call to calculateOverallFitness takes much more time than a call to calculateMassScaleFitness.
To let the GA optimize a lens model according to your own fitness criteria, you basically need to do two things. The first is creating a class derived from grale::LensFitnessObject. It is this object that will evaluate a specific lens model and calculate one or more associated fitness measures. The second part is creating a so-called module for the genetic algorithm infrastructure. We'll get to the latter in the next section, first we'll focus on the LensFitnessObject part.
If you want to create your own lens inversion module with your own LensFitnessObject derived code, it is best to work in the grale_gamodules package. The documentation is a part of GRALE since all the base classes are a part of it, but to actually program a new lens inversion module you should not modify anything from the GRALE code, which can just be compiled and installed as usual. The actual programming can all be done like the other modules in the grale_gamodules package.
The first thing you'll need to do is derive a class from grale::LensFitnessObject. This is an abstract class, so you'll need to provide an implementation of all the member functions. The first function to consider is the grale::LensFitnessObject::init function, which takes two parameters: images and massScaleImages.
In GRALESHELL, you can create a list of images data files, and when starting the actual lens inversion GA, this list is passed as the images parameter of the grale::LensFitnessObject::init function in your derived class. The second list, massScaleImages, is empty when the function is called, and it is up to you to store some members from images in this list, if necessary.
The init function can remove some elements of the images list, and for all remaining entries some properties will be calculated. Usually you won't need to modify the list, but I mention it for completeness. The properties that are calculated can consist of the backprojected positions of image points, derivatives of the deflection angles at those points etc. All this information will be presented to your fitness calculation code later in the form of a grale::ProjectedImagesInterface instance.
As mentioned before, the fitness calculations happen in two steps: first a weight for the basis functions are calculated, and then, using the fully determined mass model, the full fitness calculation occurs. In the massScaleImages list, you should store those entries of the images list that you need in the mass scale fitness calculation. If you leave the list empty, it is assumed to be the same list as the images list.
In all the GA modules for lens inversion that currently exist, the mass scale fitness is based on the backprojected images alone. The fitness is the positional fitness described in the articles, so that the mass scale is fixed to yield the model that corresponds to the best positional fitness measure. Using this model, the other fitness measures (null space, time delay, ...) are then evaluated.
After processing the input data in the init function and providing a massScaleImages list if necessary, you'll need to create the two fitness-related functions. The grale::LensFitnessObject::calculateMassScaleFitness in your implementation should evalutate the backprojected image data etc in the grale::ProjectedImagesInterface instance that was passed as an argument. Based on this information it should return a value which will be used to optimize the scale of the basis functions (a lower value is considered better).
The second function you'll need to implement is the grale::LensFitnessObject::calculateOverallFitness function. The first argument is again a grale::ProjectedImagesInterface instance and it is these data you should evalutate. All the fitness measures you calculate should be stored in the pFitnessValues array that is the second argument. Possibly one of the fitness measures will correspond to the calculation you did in the calculateMassScaleFitness function, so you can call this function again and avoid writing the code twice.
The length of the array pFitnessValues that is passed to the function will correspond to the value that you return in grale::LensFitnessObject::getNumberOfFitnessComponents. As the name suggests, this function should return the number of fitness components in your procedure. This is one for a single-objective algorithm or more for a multi-objective GA.
The grale::ProjectedImagesInterface abstract class has two full implementations: grale::ImagesBackProjector and grale::BackProjectMatrixNew. The first one is used by the GRALESHELL command 'imgdata/fitness', to calculate the fitness values for the currently active lens model and the current images data list. This command performs the following actions:
images list from the init function As you can see, for this call only the full fitness measure is calculated. This is only natural of course, since the mass model is already fully fixed and no scale factor needs to be calculated. The grale::ImagesBackProjector instance also makes all possible calculations in the grale::ProjectedImagesInterface available, so this is a good way to test your fitness calculation code. Inside the GA, not all of those functions may give valid results, it depends on the values you return in other member functions of grale::LensFitnessObject.
As mentioned above, inside the GA a grale::BackProjectMatrixNew instance is used to describe the situation created by the current basis functions and mass scale. This instance is much more optimized than the grale::ImagesBackProjector class, and will only calculate/store those things that you specify using other member functions of the grale::LensFitnessObject class.
First, there's a number of functions which decide what should be calculated or stored for the images data list images from the init function. The function grale::LensFitnessObject::getTotalCalcFlags is called with three empty lists of booleans. It is up to you to resize these to have the same length as the images list from the init function. For each entry, they should indicate if that source needs deflection calculations for the image points, the derivatives of the deflection angle at those points or even the lens potential. The grale::LensFitnessObject::getTotalStoreFlags indicate what properties should be stored from the original images list. This can same some memory if necessary, but if this isn't an issue you can just set all of these flags to true. It won't affect the speed at which the GA is executed.
You can also indicate if some extra properties that depend on the derivatives of the deflection angle should be calculated. These can be the inverse magnification, the shear components and the convergence. If a function
returns false, it is assumed that the corresponding calculations can be skipped for all sources. If true is returned, a further call to the functions
is made to determine for which precise sources the calculations need to be done. The list returned should again be of the same length as the images list from the init function.
There is a very similar set of functions for the massScaleImages list from the init function. They have slightly different names and refer to this other list of images data, but otherwise their meaning is the same.
This is the file that will be loaded when specify a specific algorithm in the 'invert/grid/local' command in GRALESHELL for example. Such a module will not only contain the LensFitnessObject you created, but also some other code. The easiest way to get started is to take a look at the templates that are a part of the grale_gamodules package.
There's a template for both a single and multi-objective algorithm. First and foremost, the class defined in such a file derives from grale::GridLensInversionGAFactoryBase. This is a genetic algorithm 'factory' that's derived from a class in MOGAL and which contains most of the GA specific code. The templates contain some documentation which tells you what each function needs to do, which usually isn't that much.
One thing you'll definitely need is some kind of stopping criterion that signal the GA when the optimization is done. For a multi-objective algorithm you'll also need to provide a function that picks a genome from a 'non-dominated' set. When optimizing multiple criteria, each generation no longer has a single best genome. Instead, there's a set of genomes called the non-dominated set. At some point, the GA will have to come up with a single solution though, so you'll have to provide a function that picks a single genome from a non-dominated set.
Apart from this class, your module will also need to contain two functions, within an extern "C" block. The first function is called mogal::GAFactory *CreateFactoryInstance() This function should create a new instance of the class you defined, and return this object. The second function is grale::LensFitnessObject *CreateFitnessObject() and should allocate an instance of the grale::LensFitnessObject based class you wrote to evaluate the fitness.
To make sure that the new module is created, take a look at the CMakeLists.txt file to see how it's done for other modules, and just add a similar line for your own module. Note that the module name should always start with galens_ as this is assumed by GRALESHELL. Just typing 'make' again in the build directory should then create your module. To test it, it's usually easiest to create some lens model in GRALESHELL, create an images data list that you wish to use as input, and call the 'imgdata/fitness' command with your modulename (the part after galens_ that you used as the file name). If that seems to work, you can then try the generic algorithm using your module, but remember that a lot more functions are used here, so it's usually much more difficult to debug).
In case I forgot to mention it above: with the fitness values it is always assumed that less is better.
About writing multi-objective algorithms. This is usually tricky. Furthermore, currently it will only work if all the criteria are non-conficting. For example, when using the null-space fitness measure in combination with the overlap fitness, this works since the two criteria are quite compatible: you can have a solution that both has a good overlap and a good null space value, and, as the overlap improves, the null space fitness will also usually improve (a better source estimate will cause fewer additional images).
If you're using criteria that are not compatible, it will probably not well. For example, you could combine the overlap fitness criterion with some criterion testing the smoothness of the solution. At the one end you'll then have a completely flat mass distribution that's optimal with respect to one criterion, and at the other end a solution that causes a good overlap but definitely isn't flat. To handle such cases well, you'd need to implement additional techniques that guarantee some diversity in the non-dominated set yourself, since these currently are not present in MOGAL. But note that trying to do such an optimization will make it inherently difficult to select a single best solution anyway, so in my view it's usually a better idea to put some extra thought into it and try to create criteria which are easier to combine.
1.8.1.1