diff --git a/README.md b/README.md index 0f3a9bb..db7d1d5 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,7 @@ Stands for parallel magnetic field simulations. This repository offers a user fr - The magnetic flux density **B**, the field strength **H**, the polarization **J** and the magnetization **M** - Parametrised arrangements of them, built from the parameter file rather than written out one by one - Parallel magnetic force and torque simulation between magnets +- A genetic optimizer for the homogeneity of a Halbach cylinder, which writes its answer as an input file the simulator runs - A Qt viewer for the result, which is a separate program and needs neither Kokkos nor the input file - Easy and seamless workflow using a JSON parameter file - [Modern CMake practices](https://pabloariasal.github.io/2018/02/19/its-time-to-do-cmake-right/) @@ -63,16 +64,16 @@ The `/full_dir_path_to_data_output` is the directory in which the simulation res ### Build and run the standalone target -The compile script builds the simulator, the test suite and, when Qt 6 is -installed, the viewer: +The compile script builds the simulator, the Halbach optimizer, the test suite +and, when Qt 6 is installed, the viewer: ```bash ./compile.sh ``` It takes a few options: `--clean` starts over, `--debug` builds Debug, -`--no-tests` and `--no-gui` leave parts out, `-j N` sets the number of build -jobs. `./compile.sh --help` lists them all. +`--no-tests`, `--no-optimizer` and `--no-gui` leave parts out, `-j N` sets the +number of build jobs. `./compile.sh --help` lists them all. Then run an example: @@ -90,6 +91,17 @@ cmake --build build/standalone -j ./build/standalone/Greeter --help ``` +### Build and run the Halbach optimizer + +```bash +cmake -S optimizer -B build/optimizer -DCMAKE_BUILD_TYPE=Release -DKokkos_ENABLE_OPENMP=On -DCMAKE_CXX_COMPILER=g++ +cmake --build build/optimizer -j +./build/optimizer/halbach-optimizer --help +``` + +See [Optimizing the homogeneity of a Halbach cylinder](#optimizing-the-homogeneity-of-a-halbach-cylinder) +for what it does and what it writes. + ### Build and run test suite Use the following commands from the project's root directory to run the test suite. @@ -416,6 +428,155 @@ adding an arrangement never renames a magnet that was already there. They can be individually by those ids, or all at once by the arrangement they belong to, see below. See `arrangements.json` for a complete example. +### Optimizing the homogeneity of a Halbach cylinder + +A Halbach cylinder for imaging is built as a stack of rings, and how uniform the field in the +middle comes out depends on which ring diameter is used at which position along the stack. Picking +those diameters is a search, and `halbach-optimizer` runs it: + +```bash +./build/optimizer/halbach-optimizer -o output/halbach_optimized.json +``` + +With no configuration file the defaults are the design of `HalbachOptimisation/`, the genetic +algorithm this is a port of: 23 rings 22 mm apart, each position choosing one of 19 +(radius, magnet count) pairs, in two layers 21 mm apart, measured over a 200 mm sphere. + +What comes out is **an input file of this project**, not a report. It has `arrangements` and a +`field_of_view`, so the simulator runs it as it stands: + +```bash +./build/standalone/Greeter -i output/halbach_optimized.json -s output/snapshot.json +``` + +What the run found travels beside them in a `halbach_optimization` object, which the readers +ignore. That object is also a configuration file, so the output of one run is the input of the +next: + +```bash +./build/optimizer/halbach-optimizer -c output/halbach_optimized.json --seed 7 +``` + +#### How it works + +The search never simulates a magnet. It rests on superposition: the field of a set of permanent +magnets is the sum of their fields, and the rings at one position do not move when the choice at +another position changes. So the run has three stages, and only the middle one is a search. + +1. **Every candidate is sampled once.** Each (ring position, ring candidate) pair — 12 × 19 of + them for the defaults — is laid out through the same `halbach_ring` arrangement an input file + goes through, and its field is worked out at every observation point. The magnets of all 228 + configurations are packed end to end into one parameter array first, so this is a *single* + Kokkos parallel region a million iterations wide rather than 228 simulations. + +2. **A genetic algorithm picks one candidate per position.** A genome is 12 small integers, and + its field is the sum of 12 precomputed arrays, so scoring an individual is a walk over memory. + One team per individual, spread over the observation points, taking the peak, the trough and + the total in one pass. Selection, crossover and mutation are one further parallel region a + generation. + +3. **The answer is measured again, over the whole sphere, from the magnets themselves**, through + the same simulator everything else here uses. The search is reduced to an octant by default, + which a ring of a finite number of magnets is only nearly symmetric enough for, so the figure + that is optimized is not allowed to be the figure that is reported. + +#### Configuration + +`halbach_optimization.json` in the project root is the full set of defaults, written out. Every +key is optional: + +| Key | Meaning | +| --- | --- | +| `ring_count`, `ring_separation` | the stack of rings, centred on the origin | +| `candidates` | the `{"radius", "count"}` pairs a ring position chooses between | +| `outer_radius_offset`, `outer_count_offset` | the second layer, relative to the first | +| `order` | order of the rings, as `halbach_ring` means it; 1 is the dipolar ring | +| `element` | `size` of the cube [m] and `remanence` [T] | +| `dsv`, `resolution` | the sampled sphere and the grid over it [m] | +| `symmetry` | `octant`, `hemisphere` or `full` — how much of the sphere the search measures | +| `objective` | `ppm_bx` (the x component) or `ppm_bmag` (`|B|`) | +| `field_model` | `cuboid` for the analytic kernel, `dipole` for the far field | +| `genetic` | `population`, `generations`, `cx_prob`, `mut_prob`, `gene_mut_prob`, `tournament`, `elitism`, `seed` | + +The command line overrides the file: `--generations`, `--population`, `--seed`, `--symmetry`, +`--objective`, `--field-model`. `--emit magnets` writes every magnet out one by one instead of the +46 rings. `--history convergence.csv` writes the convergence curve. `--no-verify` skips stage +three, for a smoke test. `halbach-optimizer --help` lists them all. + +`--field-model dipole` replaces each cube with the point dipole of the same moment, +`m = V·Br/µ0`, which is the model `HalbachOptimisation/halbachFields.py` uses. It is there to +compare against that script's numbers; `cuboid`, the default, is the field the magnets actually +make. + +#### Differences from the Python it ports + +- `deap.tools.mutFlipBit` applies `not` to whatever it is given, which turns a gene naming one of + 19 candidates into a 0 or a 1 and collapses most of the search space. A mutated gene is redrawn + from the alphabet here. +- The best few individuals are carried into the next generation. The Python keeps its best in a + variable on the side and lets the population lose it. +- The homogeneity is peak-to-peak over the *magnitude* of the mean. Dividing by the signed mean, + as the Python does, rewards a genome that flips the sign of the field rather than measuring it. +- `innerRingRadii` has 19 entries there and `innerNumMagnets` has 20; the zip silently drops the + last. The pairing is written out here, and lists of different lengths are refused. +- The basis keeps all three field components rather than only `Bx`, which costs nothing and is + what makes `objective: ppm_bmag` possible. + +#### What to expect + +Two things are worth watching in the output. The first is that the octant figure and the +whole-sphere figure agree — on the default design they land within about 0.01% of each other, +which is the evidence that the symmetry reduction is sound. The second is that the field the +written file produces, run back through `Greeter`, gives the homogeneity the optimizer reported; +a test asserts exactly that. + +On the default design the run settles at about **805 ppm over the 200 mm sphere**, a mean field of +49.1 mT and a peak-to-peak spread of 0.04 mT, from a magnet of 2 971 cubes. It gets most of the +way there inside 40 generations. `--field-model dipole`, the model the Python uses, follows the +same trajectory generation for generation and lands at 651 ppm — close enough to say the two +models agree about the shape of the problem, far enough apart to say which one to trust. + +#### What it costs + +The default design, on four cores of an Intel i5-8365U with `Kokkos_ENABLE_OPENMP=On`: + +| Stage | What it does | Time | +| --- | --- | ---: | +| Basis | 228 configurations, 54 625 magnets, sampled at 4 662 points | 71 s | +| Evolution | 100 generations of 10 000 individuals | 24 s | +| Verification | 2 971 magnets at 33 401 points | 27 s | +| | **total** | **2 min 3 s** | + +The two stages that touch magnets dominate, and the search — a hundred generations of ten +thousand individuals, some five billion field samples summed — is a fifth of the run. That is the +whole point of precomputing the basis: the genetic algorithm is no longer where the time goes. + +`--field-model dipole` replaces the analytic cuboid kernel, with its two dozen transcendentals per +evaluation, by a handful of multiplies, and the two magnet stages fall to 2.8 s and 1.5 s. The +search itself is unchanged, since it never looks at a magnet either way. + +Across thread counts, on the same machine (`--kokkos-num-threads=N`, 20 generations): + +| Threads | Basis | Evolution | +| ---: | ---: | ---: | +| 1 | 199 s | 9.8 s | +| 2 | 133 s | 11.4 s | +| 4 | 131 s | 6.6 s | +| 8 | 75 s | 5.3 s | + +Read those with the machine in mind: it has four physical cores, and a U-series part running one +thread at nearly 3.8 GHz drops to about 2.2 GHz with all of them busy, so a good part of what four +cores gain they give back in clock. The measure that is not confounded by that is how busy the +cores are kept — the full run above spends 13 min 8 s of CPU time in 2 min 3 s of wall clock, +6.4 of the 8 hardware threads. On a machine that holds its clocks the table will look better than +this one does. + +Note that a run is repeatable for a given seed **and a given thread count**, but not across thread +counts: the genetic operators draw from a thread-partitioned random pool, so eight threads and +four do not walk the same sequence. Neither does a floating point sum of a hundred thousand +samples come out bit for bit the same when it is split up differently. The answers agree on the +design; they need not agree on the last digit. + ### Force and torque simulation Add an optional `force` section to the input file to compute the magnetic force (in Newton) and the diff --git a/all/CMakeLists.txt b/all/CMakeLists.txt index 8092382..7960dc6 100644 --- a/all/CMakeLists.txt +++ b/all/CMakeLists.txt @@ -11,5 +11,6 @@ include(../cmake/tools.cmake) enable_testing() add_subdirectory(${CMAKE_CURRENT_LIST_DIR}/../standalone ${CMAKE_BINARY_DIR}/standalone) +add_subdirectory(${CMAKE_CURRENT_LIST_DIR}/../optimizer ${CMAKE_BINARY_DIR}/optimizer) add_subdirectory(${CMAKE_CURRENT_LIST_DIR}/../test ${CMAKE_BINARY_DIR}/test) add_subdirectory(${CMAKE_CURRENT_LIST_DIR}/../documentation ${CMAKE_BINARY_DIR}/documentation) diff --git a/compile.sh b/compile.sh index 4a772b2..a1964c0 100755 --- a/compile.sh +++ b/compile.sh @@ -1,12 +1,13 @@ #!/usr/bin/env bash -# Builds ParaMagneticS: the standalone simulator, the test suite and, when -# Qt 6 is present, the viewer. Everything lands in build/. +# Builds ParaMagneticS: the standalone simulator, the Halbach optimizer, the +# test suite and, when Qt 6 is present, the viewer. Everything lands in build/. # -# ./compile.sh simulator + tests, and the viewer if Qt 6 is there +# ./compile.sh simulator, optimizer + tests, and the viewer if Qt 6 is there # ./compile.sh --clean wipe build/ and start over # ./compile.sh --debug Debug build instead of Release # ./compile.sh --no-tests skip the test suite +# ./compile.sh --no-optimizer skip the Halbach optimizer # ./compile.sh --gui insist on the viewer, fail without Qt 6 # ./compile.sh --no-gui skip the viewer even with Qt 6 installed # ./compile.sh -j N build with N jobs, all cores by default @@ -18,16 +19,18 @@ root="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" build_type=Release jobs="$(nproc 2>/dev/null || echo 4)" with_tests=1 +with_optimizer=1 with_gui=auto clean=0 -usage() { sed -n '3,12p' "$0" | sed 's/^# \{0,1\}//'; } +usage() { sed -n '3,14p' "$0" | sed 's/^# \{0,1\}//'; } while [[ $# -gt 0 ]]; do case "$1" in --clean) clean=1 ;; --debug) build_type=Debug ;; --no-tests) with_tests=0 ;; + --no-optimizer) with_optimizer=0 ;; --gui) with_gui=1 ;; --no-gui) with_gui=0 ;; -j|--jobs) jobs="$2"; shift ;; @@ -59,6 +62,12 @@ echo "== Simulator ==" cmake -S "$root/standalone" -B "$root/build/standalone" "${flags[@]}" cmake --build "$root/build/standalone" -j "$jobs" +if [[ "$with_optimizer" == 1 ]]; then + echo "== Halbach optimizer ==" + cmake -S "$root/optimizer" -B "$root/build/optimizer" "${flags[@]}" + cmake --build "$root/build/optimizer" -j "$jobs" +fi + if [[ "$with_tests" == 1 ]]; then echo "== Tests ==" # doctest 2.4.9 asks for a CMake older than CMake 4 will agree to, hence the @@ -76,6 +85,7 @@ fi echo echo "Built:" echo " $root/build/standalone/Greeter" +[[ "$with_optimizer" == 1 ]] && echo " $root/build/optimizer/halbach-optimizer" [[ "$with_tests" == 1 ]] && echo " $root/build/test/GreeterTests" [[ "$with_gui" == 1 ]] && echo " $root/build/gui/paramagnetics-viewer" echo diff --git a/halbach_optimization.json b/halbach_optimization.json new file mode 100644 index 0000000..c809f53 --- /dev/null +++ b/halbach_optimization.json @@ -0,0 +1,59 @@ +{ + "halbach_optimization": { + + "_comment": "Every value here is the default, so an optimizer run with no configuration file at all is this file. It is the design of HalbachOptimisation/homogeneityOptimisation.py: a cylinder of 23 rings 22 mm apart, each position choosing one of 19 (radius, magnet count) pairs, in two layers 21 mm apart.", + + "ring_count": 23, + "ring_separation": 0.022, + + "candidates": [ + {"radius": 0.148, "count": 50}, + {"radius": 0.151, "count": 51}, + {"radius": 0.154, "count": 52}, + {"radius": 0.156, "count": 53}, + {"radius": 0.159, "count": 54}, + {"radius": 0.162, "count": 55}, + {"radius": 0.165, "count": 56}, + {"radius": 0.168, "count": 57}, + {"radius": 0.171, "count": 58}, + {"radius": 0.174, "count": 59}, + {"radius": 0.177, "count": 60}, + {"radius": 0.180, "count": 61}, + {"radius": 0.183, "count": 62}, + {"radius": 0.186, "count": 63}, + {"radius": 0.189, "count": 64}, + {"radius": 0.192, "count": 65}, + {"radius": 0.195, "count": 66}, + {"radius": 0.198, "count": 67}, + {"radius": 0.201, "count": 68} + ], + + "outer_radius_offset": 0.021, + "outer_count_offset": 7, + + "order": 1, + + "element": { + "size": 0.012, + "remanence": 1.3 + }, + + "dsv": 0.2, + "resolution": 0.005, + + "symmetry": "octant", + "objective": "ppm_bx", + "field_model": "cuboid", + + "genetic": { + "population": 10000, + "generations": 100, + "cx_prob": 0.55, + "mut_prob": 0.4, + "gene_mut_prob": 0.05, + "tournament": 3, + "elitism": 1, + "seed": 42 + } + } +} diff --git a/include/greeter/io/HalbachOptimizationIO.h b/include/greeter/io/HalbachOptimizationIO.h new file mode 100644 index 0000000..cb532c7 --- /dev/null +++ b/include/greeter/io/HalbachOptimizationIO.h @@ -0,0 +1,81 @@ +#ifndef HALBACH_OPTIMIZATION_IO_H +#define HALBACH_OPTIMIZATION_IO_H + +#include +#include +#include +#include + + +namespace greeter { + +/* Whether the written file lists the rings or the magnets they stand for. */ +enum class HalbachEmit { + + /* + One "halbach_ring" per ring: forty six objects for a magnet of two + thousand nine hundred, and the form the optimizer actually decides in, + so a reader can see which candidate went where. + */ + Arrangements, + + /* Every magnet written out one by one, for a reader that wants them. */ + Magnets +}; + +/* + Reads what an optimization is asked to do, and writes what it decided. + + What comes out is an input file of this library, not a report: it has + "arrangements" and a "field_of_view", it passes MagnetIO::validateJSON, and + the standalone simulator runs it as it stands. What the optimization found + travels in a "halbach_optimization" object beside them, which the readers + ignore, so the one file is both the answer and the record of how it was + arrived at. +*/ +class HalbachOptimizationIO { + + public: + + HalbachOptimizationIO(); + ~HalbachOptimizationIO(); + + /* + Reads a specification. + + Either the whole file is the configuration, or it carries it under + "halbach_optimization", which is what the written file does. That is + what lets the output of one run be the input of the next. + */ + static greeter::optimization::HalbachSpec readSpec(const nlohmann::json& data); + + static greeter::optimization::HalbachSpec readSpecFile(const std::string& path); + + /* The solution as an input file that describes the optimized magnet. */ + static nlohmann::json write( + const greeter::optimization::HalbachSolution& solution, + const HalbachEmit& emit); + + static void writeFile( + const greeter::optimization::HalbachSolution& solution, + const HalbachEmit& emit, const std::string& path); + + /* The convergence curve, one row a generation. */ + static void writeHistoryCSV( + const greeter::optimization::HalbachSolution& solution, + const std::string& path); + + /* + One magnet in the schema of an input file. + + Only the shapes an optimized Halbach magnet is built out of, which is + the cuboid and the dipole. Anything else is refused by name rather than + written out wrongly. + */ + static nlohmann::json writeMagnet(const greeter::MagnetCollection& collection, + const size_t& index, const int64_t& id); +}; + +} // namespace greeter + +#endif // HALBACH_OPTIMIZATION_IO_H diff --git a/include/greeter/optimization/GeneticOptimizer.h b/include/greeter/optimization/GeneticOptimizer.h new file mode 100644 index 0000000..06e6900 --- /dev/null +++ b/include/greeter/optimization/GeneticOptimizer.h @@ -0,0 +1,105 @@ +#ifndef GENETIC_OPTIMIZER_H +#define GENETIC_OPTIMIZER_H + +#include +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* What one generation came to, for a convergence curve. */ +struct GenerationRecord { + uint32_t generation = 0; + float best = 0.0f; // best of this generation, in ppm + float best_ever = 0.0f; // best seen since the run began + float mean = 0.0f; // mean of the generation, a spread indicator + float seconds = 0.0f; +}; + +/* Told how a run is going, and asked whether to carry on. */ +class GenerationSink { + + public: + + virtual ~GenerationSink() = default; + + /* Returns false to stop the run at the end of this generation. */ + virtual bool onGeneration(const GenerationRecord& record) = 0; +}; + +/* + A genetic algorithm over genomes of small integers, run on the device. + + A genome is one index per gene, drawn from a common alphabet: for the + Halbach problem, which ring candidate each ring position uses. This is the + port of the DEAP loop of homogeneityOptimisation.py, with the same operators + at the same rates, but with the population living in a Kokkos view and a + whole generation produced by a single parallel region. + + The one operator that is not a port is the mutation. DEAP's mutFlipBit + applies "not" to whatever it is given, which turns a gene of the nineteen + candidates into a zero or a one and quietly collapses most of the search + space. A mutated gene is redrawn from the alphabet here, which is what a + mutation on an integer genome is meant to do. +*/ +class GeneticOptimizer { + + public: + + GeneticOptimizer(const GeneticSettings& settings, + const size_t& num_genes, const size_t& num_candidates); + + ~GeneticOptimizer(); + + /* + Runs the evolution and returns the best genome it ever saw. + + "Ever" rather than "at the end": a generation can lose its best + individual to a crossover, so the best of the last generation is not + the best of the run. The elite carried over makes the two agree in + practice, and this makes them agree by construction. + */ + std::vector run(const HomogeneityObjective& objective, + GenerationSink* sink, const bool& verbose); + + float getBestFitness() const; + + const std::vector& getHistory() const; + + /* Whether a sink asked the run to stop before the last generation. */ + bool wasStopped() const; + + private: + + /* + Selection, crossover and mutation of the whole population, in one + parallel region over pairs of offspring. They fuse because selection + only reads the fitness of the generation that is being replaced, so + nothing here has to wait for anything else here. + */ + void reproduce(const PopulationView& from, const PopulationView& into, + const FitnessView& fitness); + + GeneticSettings settings; + + size_t num_genes; + size_t num_candidates; + + RandPoolType random_pool; + + std::vector best_genome; + float best_fitness; + + std::vector history; + + bool stopped = false; +}; + +} // namespace optimization +} // namespace greeter + +#endif // GENETIC_OPTIMIZER_H diff --git a/include/greeter/optimization/HalbachBasis.h b/include/greeter/optimization/HalbachBasis.h new file mode 100644 index 0000000..fc9176e --- /dev/null +++ b/include/greeter/optimization/HalbachBasis.h @@ -0,0 +1,100 @@ +#ifndef HALBACH_BASIS_H +#define HALBACH_BASIS_H + +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* + The field every choice the optimizer can make would produce, worked out once. + + The whole optimization rests on one fact: the field of a set of permanent + magnets is the sum of their fields. A genome picks one ring candidate per + ring position, and the rings of one position do not move when the choice at + another position changes. So the field of a genome is the sum of as many + precomputed fields as there are genes, and the genetic algorithm never has + to evaluate a magnet again after this class has run. + + What is stored is + + fields(configuration, component, point) + + where a configuration is a (gene, candidate) pair flattened as + + configuration = gene * num_candidates + candidate + + and holds every ring that pair stands for: both layers, at +z and -z. The + point index is last because it is the one that is walked, so LayoutRight + makes each component of each configuration a contiguous run: a vector loop + on a CPU, a coalesced read on a GPU. + + All three components are kept rather than the x component the Python script + keeps. Superposition is linear per component, so summing the three and + taking the norm afterwards is exact, and that is what makes an objective on + |B| available at no cost beyond the memory. +*/ +class HalbachBasis { + + public: + + HalbachBasis(); + ~HalbachBasis(); + + /* + Samples the field of every candidate configuration. + + One Kokkos kernel over (configuration, observation point), not one + simulation per configuration: the magnets of every configuration are + packed end to end first, exactly as MagnetCollection::fillMagnetParameters + packs one collection, so that the whole precomputation is a single + parallel region a million iterations wide. + */ + static HalbachBasis build(const HalbachSpec& spec, const bool& verbose); + + /* The observation points that were sampled, three floats each. */ + const std::vector>& getPoints() const; + + size_t getNumPoints() const; + size_t getNumGenes() const; + size_t getNumCandidates() const; + + /* How many magnets went into the precomputation, for a report. */ + size_t getNumMagnets() const; + + Kokkos::View getFields() const; + + /* + The field of one genome at one point, on the host, summed the plain way. + The device objective is checked against this. + */ + void evaluateAt(const std::vector& genome, const size_t& point, + float* b) const; + + private: + + Kokkos::View fields; + + std::vector> points; + + size_t num_genes = 0; + size_t num_candidates = 0; + size_t num_magnets = 0; +}; + +/* + The points of the sampled volume, in the order the basis stores them. + + A regular grid over the box of the DSV, thinned to the points that lie + inside the sphere and on the side of it the symmetry asks for. The Python + script builds the same set with a numpy mask. +*/ +std::vector> makeSamplePoints(const HalbachSpec& spec); + +} // namespace optimization +} // namespace greeter + +#endif // HALBACH_BASIS_H diff --git a/include/greeter/optimization/HalbachOptimizer.h b/include/greeter/optimization/HalbachOptimizer.h new file mode 100644 index 0000000..e405569 --- /dev/null +++ b/include/greeter/optimization/HalbachOptimizer.h @@ -0,0 +1,69 @@ +#ifndef HALBACH_OPTIMIZER_H +#define HALBACH_OPTIMIZER_H + +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* How a run is carried out, as against what it is looking for. */ +struct RunOptions { + + /* Whether every step announces itself on the standard output. */ + bool verbose = true; + + /* + Whether the answer is measured again over the whole sphere, from the + magnets rather than from the basis. + + On by default. It costs one more field simulation, and it is the only + part of a run that does not take the symmetry reduction on trust, so + turning it off is for a smoke test and not for a result. + */ + bool verify = true; +}; + +/* + The whole optimization, from a specification to a magnet. + + Three stages, and the middle one is the only one that is a search: + + 1. Every ring candidate at every ring position is simulated once, into + HalbachBasis. This is where the magnets are, and it is a single + parallel region a million iterations wide. + + 2. A genetic algorithm picks one candidate per position. It never touches + a magnet: superposition means the field of a choice is the sum of the + precomputed fields, so a generation of ten thousand is a walk over + memory rather than a simulation. + + 3. The choice is turned back into "halbach_ring" arrangements and + simulated over the whole sphere, which is the number the run reports. +*/ +class HalbachOptimizer { + + public: + + static HalbachSolution run(const HalbachSpec& spec, GenerationSink* sink, + const RunOptions& options); + + /* + Measures a magnet that already exists, the way stage three does. + + Given to the solution's own collection it is the verification, and given + to any other it is the comparison that says whether the optimization was + worth doing. + */ + static FieldMetrics measure(const greeter::MagnetCollection& collection, + const HalbachSpec& spec, + const Objective& objective, + size_t& num_points); +}; + +} // namespace optimization +} // namespace greeter + +#endif // HALBACH_OPTIMIZER_H diff --git a/include/greeter/optimization/HalbachSolution.h b/include/greeter/optimization/HalbachSolution.h new file mode 100644 index 0000000..1e9dd1d --- /dev/null +++ b/include/greeter/optimization/HalbachSolution.h @@ -0,0 +1,69 @@ +#ifndef HALBACH_SOLUTION_H +#define HALBACH_SOLUTION_H + +#include +#include +#include +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* + What the optimizer settled on, and how good it turned out to be. + + There are two homogeneity figures here, and they answer different questions. + `optimized` is what the genetic algorithm was steered by: the objective over + the reduced volume, summed out of the precomputed basis. `verified` is the + same figure taken again over the whole sphere, from the magnets themselves, + through the simulator every other part of this library uses. + + They are kept apart on purpose. The reduction to an octant assumes a + symmetry that a ring of a finite number of magnets only nearly has, and a + figure that was optimized cannot also be the evidence that the optimization + worked. The second number is the one to quote. +*/ +struct HalbachSolution { + + HalbachSpec spec; + + /* One ring candidate index per ring position of the symmetric half. */ + std::vector genome; + + FieldMetrics optimized; // over the reduced volume, out of the basis + size_t optimized_points = 0; + + FieldMetrics verified; // over the whole sphere, out of the simulator + size_t verified_points = 0; + bool was_verified = false; + + size_t num_magnets = 0; + + /* The size of stage one: every candidate at every position, sampled. */ + size_t basis_configurations = 0; + size_t basis_magnets = 0; + + std::vector history; + + double basis_seconds = 0.0; + double evolution_seconds = 0.0; + double verification_seconds = 0.0; + + /* The rings the genome chose, in the schema an input file writes. */ + std::vector buildArrangements() const; + + /* Those rings expanded into magnets, which is what a simulation runs. */ + greeter::MagnetCollection buildCollection() const; + + /* The candidates the genome chose, for a report. */ + std::vector getChosenCandidates() const; +}; + +} // namespace optimization +} // namespace greeter + +#endif // HALBACH_SOLUTION_H diff --git a/include/greeter/optimization/HalbachSpec.h b/include/greeter/optimization/HalbachSpec.h new file mode 100644 index 0000000..f3a25e4 --- /dev/null +++ b/include/greeter/optimization/HalbachSpec.h @@ -0,0 +1,214 @@ +#ifndef HALBACH_SPEC_H +#define HALBACH_SPEC_H + +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* + Which quantity the homogeneity is measured on. + + The basis fields are stored per component and summed per component, so both + of these are linear in the ring contributions right up to the last step. + Only the norm of BMagnitude is not, and it is taken once, after the sum. +*/ +enum class Objective { + Bx, // the x component alone, which is what the Python script uses + BMagnitude // |B|, the quantity an MRI homogeneity figure usually means +}; + +/* Which part of the sampled sphere the objective is evaluated over. */ +enum class Symmetry { + Octant, // x, y, z all >= 0, the eighth the Python script uses + Hemisphere, // z >= 0 + Full +}; + +/* + Which field model the basis is built with. + + Cuboid is the analytic kernel of this library, which is what the magnets + actually are. Dipole is the far field approximation the Python script uses, + kept because it is the only way to compare against the numbers that script + printed. +*/ +enum class FieldModel { + Cuboid, + Dipole +}; + +/* + One entry of the list the optimizer chooses from. + + A gene of the genome is an index into this list, and picking entry s for + ring position p means: an inner ring of `count` magnets at `radius`, and an + outer ring of `count + outer_count_offset` magnets at + `radius + outer_radius_offset`. The two layers move together, which is what + makes one gene enough per position. +*/ +struct RingCandidate { + float radius = 0.0f; // [m], inner layer + uint32_t count = 0; // magnets in the inner ring +}; + +/* + The magnet the rings are built out of. + + Only the side length is asked for, because a Halbach ring of cubes is what + both the Python script and every arrangement in this repository build, and + because the dipole model has to derive a moment from a volume. Everything + else about the element follows from the field model. +*/ +struct Element { + float size = 0.012f; // [m], side of the cube + float remanence = 1.3f; // [T], Br of the material +}; + +/* How the genetic algorithm is run. */ +struct GeneticSettings { + + uint32_t population = 10000; + uint32_t generations = 100; + + float cx_prob = 0.55f; // chance a pair is crossed + float mut_prob = 0.4f; // chance an individual is mutated at all + float gene_mut_prob = 0.05f; // chance each gene of that individual moves + + uint32_t tournament = 3; + + /* + How many of the best are carried into the next generation untouched. + + The Python script keeps none, and tracks the best it ever saw by hand + while letting the population lose it. One elite costs a single genome + and makes the reported best a member of the population that produced it. + */ + uint32_t elitism = 1; + + uint64_t seed = 42; +}; + +/* + Everything the optimizer needs, and the whole of what a configuration file + can say. + + Constructed with the defaults of homogeneityOptimisation.py, so that an + optimizer run with no configuration file at all is the port of that script. +*/ +class HalbachSpec { + + public: + + HalbachSpec(); + + /* Reads what a configuration object overrides, leaving the rest alone. */ + void readJSON(const nlohmann::json& data); + + /* Throws with a message naming the field that is wrong. */ + void validate() const; + + /* + Where the rings sit along the axis, mirrored about the origin. + + A cylinder of `ring_count` rings spaced `ring_separation` apart, centred + on the origin. An even count has no ring at zero, an odd one does. + */ + std::vector getRingPositions() const; + + /* + The non negative half of the above, which is what the genome indexes. + + The magnet is built mirror symmetric in z, so a gene at a position other + than zero stands for two rings per layer rather than one. + */ + std::vector getSymmetricRingPositions() const; + + size_t getNumGenes() const; + size_t getNumCandidates() const; + + /* The moment [A m^2] a cube of this element is worth, m = V * Br / mu0. */ + float getDipoleMoment() const; + + /* + The element as an arrangement writes it: a cuboid carrying a + polarization along its local x, or a dipole carrying the equivalent + moment along the same axis. The Halbach ring turns it from there. + */ + nlohmann::json getElementJSON() const; + + /* + One "halbach_ring" arrangement, written exactly as an input file of this + library writes one. Ring layout is decided here and nowhere else, so the + magnet the optimizer measures and the magnet it writes out cannot drift + apart. + */ + nlohmann::json makeRingJSON(const float& radius, const uint32_t& count, + const float& z) const; + + /* + Every ring one gene stands for: the inner and the outer layer, at +z and, + unless the gene sits at zero, at -z as well. Two rings or four. + */ + std::vector makeRingsForGene(const size_t& gene, + const size_t& candidate) const; + + /* The whole spec back out again, for the record written with a solution. */ + nlohmann::json toJSON() const; + + uint32_t ring_count = 23; + float ring_separation = 0.022f; // [m] + + std::vector candidates; + + float outer_radius_offset = 0.021f; // [m] + uint32_t outer_count_offset = 7; + + /* + Order of the Halbach rings, as HalbachRingArrangement means it: a member + is turned by (order + 1) times the angle at which it sits. Order one is + the dipolar ring, whose field inside is uniform and transverse, and is + the kValue = 2 of the Python script. + */ + int64_t order = 1; + + Element element; + + float dsv = 0.2f; // [m], diameter of the sampled sphere + float resolution = 0.005f; // [m], spacing of the sample grid + + Symmetry symmetry = Symmetry::Octant; + Objective objective = Objective::Bx; + FieldModel field_model = FieldModel::Cuboid; + + GeneticSettings genetic; +}; + +/* + A float as a JSON number a person can read. + + Every length here is a float, and a float widened to the double a JSON + number is made of prints as seventeen digits of which ten are the noise of + the conversion: a twelve millimetre cube comes out as 0.012000000104308128. + Rounded to the seven significant digits a float is good for, it comes out + as 0.012. + + A float needs nine digits to be certain of surviving the trip, so this + gives up the last bit or two of a number that was not written as a short + decimal to begin with. That is a fraction of a nanometre on the radius of a + ring, weighed against a file nobody can read. +*/ +double tidy(const float& value); + +const char* toString(const Objective& objective); +const char* toString(const Symmetry& symmetry); +const char* toString(const FieldModel& model); + +} // namespace optimization +} // namespace greeter + +#endif // HALBACH_SPEC_H diff --git a/include/greeter/optimization/HomogeneityObjective.h b/include/greeter/optimization/HomogeneityObjective.h new file mode 100644 index 0000000..14d3fd7 --- /dev/null +++ b/include/greeter/optimization/HomogeneityObjective.h @@ -0,0 +1,151 @@ +#ifndef HOMOGENEITY_OBJECTIVE_H +#define HOMOGENEITY_OBJECTIVE_H + +#include +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +/* One integer index per gene. Sixteen bits is more candidates than anyone + would list, and keeps a population of ten thousand under a megabyte. */ +typedef Kokkos::View PopulationView; +typedef Kokkos::View FitnessView; + +/* + The peak, the trough and the total of a field, gathered together. + + Kokkos combines several reducers in one pass only for a top level policy, + not for a reduction nested inside a team, and the whole point of the + evaluation below is that it walks the volume once. So the three quantities + the homogeneity needs travel as one value with one join. + + The total is kept in double. A hundred thousand single precision samples of + nearly the same size lose, in their sum, digits that the answer is measured + in. +*/ +struct FieldExtent { + float min = 0.0f; + float max = 0.0f; + double sum = 0.0; +}; + +template +struct FieldExtentReducer { + + public: + + typedef FieldExtentReducer reducer; + typedef FieldExtent value_type; + typedef Kokkos::View result_view_type; + + KOKKOS_INLINE_FUNCTION + explicit FieldExtentReducer(value_type& _value): value(_value) {} + + KOKKOS_INLINE_FUNCTION + void join(value_type& into, const value_type& from) const { + if (from.min < into.min) { into.min = from.min; } + if (from.max > into.max) { into.max = from.max; } + into.sum += from.sum; + } + + KOKKOS_INLINE_FUNCTION + void init(value_type& into) const { + into.min = Kokkos::reduction_identity::min(); + into.max = Kokkos::reduction_identity::max(); + into.sum = Kokkos::reduction_identity::sum(); + } + + KOKKOS_INLINE_FUNCTION + value_type& reference() const { return value; } + + KOKKOS_INLINE_FUNCTION + result_view_type view() const { return result_view_type(&value); } + + KOKKOS_INLINE_FUNCTION + bool references_scalar() const { return true; } + + private: + + value_type& value; +}; + +/* What a field over the sampled volume looks like, in one line. */ +struct FieldMetrics { + + float min = 0.0f; + float max = 0.0f; + float mean = 0.0f; + + /* + Peak to peak over the mean, in parts per million, which is how the + homogeneity of a magnet is quoted. + + The mean is taken in magnitude. The Python script divides by the signed + mean, which is the same number whenever the field points the way it is + expected to, and is a trap when it does not: a genome that flipped the + sign of the mean would be rewarded rather than measured. + */ + float ppm = 0.0f; +}; + +/* + How good a genome is: how flat the field it produces is. + + The genome is not turned back into magnets to answer that. Every gene picks + one precomputed field out of HalbachBasis, and the field of the genome is + their sum, so the whole evaluation is a walk over as many arrays as there + are genes. That is the reason the basis exists and the reason a generation + of ten thousand individuals is seconds rather than hours. +*/ +class HomogeneityObjective { + + public: + + HomogeneityObjective(); + HomogeneityObjective(const HalbachBasis& basis, const Objective& objective); + + /* + Scores a whole population in one parallel region. + + A team per individual, the team spread over the observation points, and + the peak, the trough and the total taken in a single pass with three + reducers rather than three passes. The configuration each gene points at + is worked out once per team into scratch, since it is otherwise + recomputed at every one of the thousands of points. + */ + void evaluate(const PopulationView& population, const FitnessView& fitness) const; + + /* One genome, on the host, the plain way. What a report quotes. */ + FieldMetrics evaluateGenome(const std::vector& genome) const; + + /* + The same figure for a field that was sampled some other way, so that a + run can be checked against the simulator rather than against itself. + The values are the ones the objective measures: Bx, or |B|. + */ + static FieldMetrics summarise(const std::vector& values); + + size_t getNumGenes() const; + size_t getNumCandidates() const; + size_t getNumPoints() const; + + private: + + Kokkos::View fields; + + size_t num_genes = 0; + size_t num_candidates = 0; + size_t num_points = 0; + + Objective objective = Objective::Bx; +}; + +} // namespace optimization +} // namespace greeter + +#endif // HOMOGENEITY_OBJECTIVE_H diff --git a/optimized_halbach.json b/optimized_halbach.json new file mode 100644 index 0000000..435e88b --- /dev/null +++ b/optimized_halbach.json @@ -0,0 +1,1875 @@ +{ + "arrangements": [ + { + "id": 1, + "parameters": { + "count": 61, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.0 + ], + "radius": 0.18 + }, + "type": "halbach_ring" + }, + { + "id": 2, + "parameters": { + "count": 68, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.0 + ], + "radius": 0.201 + }, + "type": "halbach_ring" + }, + { + "id": 3, + "parameters": { + "count": 63, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.022 + ], + "radius": 0.186 + }, + "type": "halbach_ring" + }, + { + "id": 4, + "parameters": { + "count": 70, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.022 + ], + "radius": 0.207 + }, + "type": "halbach_ring" + }, + { + "id": 5, + "parameters": { + "count": 63, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.022 + ], + "radius": 0.186 + }, + "type": "halbach_ring" + }, + { + "id": 6, + "parameters": { + "count": 70, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.022 + ], + "radius": 0.207 + }, + "type": "halbach_ring" + }, + { + "id": 7, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.04400001 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 8, + "parameters": { + "count": 69, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.04400001 + ], + "radius": 0.204 + }, + "type": "halbach_ring" + }, + { + "id": 9, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.04400001 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 10, + "parameters": { + "count": 69, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.04400001 + ], + "radius": 0.204 + }, + "type": "halbach_ring" + }, + { + "id": 11, + "parameters": { + "count": 61, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.066 + ], + "radius": 0.18 + }, + "type": "halbach_ring" + }, + { + "id": 12, + "parameters": { + "count": 68, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.066 + ], + "radius": 0.201 + }, + "type": "halbach_ring" + }, + { + "id": 13, + "parameters": { + "count": 61, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.066 + ], + "radius": 0.18 + }, + "type": "halbach_ring" + }, + { + "id": 14, + "parameters": { + "count": 68, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.066 + ], + "radius": 0.201 + }, + "type": "halbach_ring" + }, + { + "id": 15, + "parameters": { + "count": 60, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.08799998 + ], + "radius": 0.177 + }, + "type": "halbach_ring" + }, + { + "id": 16, + "parameters": { + "count": 67, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.08799998 + ], + "radius": 0.198 + }, + "type": "halbach_ring" + }, + { + "id": 17, + "parameters": { + "count": 60, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.08799998 + ], + "radius": 0.177 + }, + "type": "halbach_ring" + }, + { + "id": 18, + "parameters": { + "count": 67, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.08799998 + ], + "radius": 0.198 + }, + "type": "halbach_ring" + }, + { + "id": 19, + "parameters": { + "count": 59, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.11 + ], + "radius": 0.174 + }, + "type": "halbach_ring" + }, + { + "id": 20, + "parameters": { + "count": 66, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.11 + ], + "radius": 0.195 + }, + "type": "halbach_ring" + }, + { + "id": 21, + "parameters": { + "count": 59, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.11 + ], + "radius": 0.174 + }, + "type": "halbach_ring" + }, + { + "id": 22, + "parameters": { + "count": 66, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.11 + ], + "radius": 0.195 + }, + "type": "halbach_ring" + }, + { + "id": 23, + "parameters": { + "count": 68, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.132 + ], + "radius": 0.201 + }, + "type": "halbach_ring" + }, + { + "id": 24, + "parameters": { + "count": 75, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.132 + ], + "radius": 0.222 + }, + "type": "halbach_ring" + }, + { + "id": 25, + "parameters": { + "count": 68, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.132 + ], + "radius": 0.201 + }, + "type": "halbach_ring" + }, + { + "id": 26, + "parameters": { + "count": 75, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.132 + ], + "radius": 0.222 + }, + "type": "halbach_ring" + }, + { + "id": 27, + "parameters": { + "count": 55, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.154 + ], + "radius": 0.162 + }, + "type": "halbach_ring" + }, + { + "id": 28, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.154 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 29, + "parameters": { + "count": 55, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.154 + ], + "radius": 0.162 + }, + "type": "halbach_ring" + }, + { + "id": 30, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.154 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 31, + "parameters": { + "count": 50, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.176 + ], + "radius": 0.148 + }, + "type": "halbach_ring" + }, + { + "id": 32, + "parameters": { + "count": 57, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.176 + ], + "radius": 0.169 + }, + "type": "halbach_ring" + }, + { + "id": 33, + "parameters": { + "count": 50, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.176 + ], + "radius": 0.148 + }, + "type": "halbach_ring" + }, + { + "id": 34, + "parameters": { + "count": 57, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.176 + ], + "radius": 0.169 + }, + "type": "halbach_ring" + }, + { + "id": 35, + "parameters": { + "count": 55, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.198 + ], + "radius": 0.162 + }, + "type": "halbach_ring" + }, + { + "id": 36, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.198 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 37, + "parameters": { + "count": 55, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.198 + ], + "radius": 0.162 + }, + "type": "halbach_ring" + }, + { + "id": 38, + "parameters": { + "count": 62, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.198 + ], + "radius": 0.183 + }, + "type": "halbach_ring" + }, + { + "id": 39, + "parameters": { + "count": 51, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.22 + ], + "radius": 0.151 + }, + "type": "halbach_ring" + }, + { + "id": 40, + "parameters": { + "count": 58, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.22 + ], + "radius": 0.172 + }, + "type": "halbach_ring" + }, + { + "id": 41, + "parameters": { + "count": 51, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.22 + ], + "radius": 0.151 + }, + "type": "halbach_ring" + }, + { + "id": 42, + "parameters": { + "count": 58, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.22 + ], + "radius": 0.172 + }, + "type": "halbach_ring" + }, + { + "id": 43, + "parameters": { + "count": 50, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.242 + ], + "radius": 0.148 + }, + "type": "halbach_ring" + }, + { + "id": 44, + "parameters": { + "count": 57, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + 0.242 + ], + "radius": 0.169 + }, + "type": "halbach_ring" + }, + { + "id": 45, + "parameters": { + "count": 50, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.242 + ], + "radius": 0.148 + }, + "type": "halbach_ring" + }, + { + "id": 46, + "parameters": { + "count": 57, + "element": { + "parameters": { + "dimensions": [ + 0.012, + 0.012, + 0.012 + ], + "magnetization": [ + 1.3, + 0.0, + 0.0 + ] + }, + "type": "cuboid" + }, + "order": 1, + "orientation": [ + 1.0, + 0.0, + 0.0, + 0.0 + ], + "position": [ + 0.0, + 0.0, + -0.242 + ], + "radius": 0.169 + }, + "type": "halbach_ring" + } + ], + "field_of_view": { + "x": { + "max": 0.1, + "min": -0.1, + "n": 41 + }, + "y": { + "max": 0.1, + "min": -0.1, + "n": 41 + }, + "z": { + "max": 0.1, + "min": -0.1, + "n": 41 + } + }, + "halbach_optimization": { + "basis": { + "configurations": 228, + "magnets": 54625, + "points": 4662 + }, + "candidates": [ + { + "count": 50, + "radius": 0.148 + }, + { + "count": 51, + "radius": 0.151 + }, + { + "count": 52, + "radius": 0.154 + }, + { + "count": 53, + "radius": 0.156 + }, + { + "count": 54, + "radius": 0.159 + }, + { + "count": 55, + "radius": 0.162 + }, + { + "count": 56, + "radius": 0.165 + }, + { + "count": 57, + "radius": 0.168 + }, + { + "count": 58, + "radius": 0.171 + }, + { + "count": 59, + "radius": 0.174 + }, + { + "count": 60, + "radius": 0.177 + }, + { + "count": 61, + "radius": 0.18 + }, + { + "count": 62, + "radius": 0.183 + }, + { + "count": 63, + "radius": 0.186 + }, + { + "count": 64, + "radius": 0.189 + }, + { + "count": 65, + "radius": 0.192 + }, + { + "count": 66, + "radius": 0.195 + }, + { + "count": 67, + "radius": 0.198 + }, + { + "count": 68, + "radius": 0.201 + } + ], + "chosen": [ + { + "candidate": 11, + "inner_count": 61, + "inner_radius": 0.18, + "outer_count": 68, + "outer_radius": 0.201, + "ring_position": 0.0 + }, + { + "candidate": 13, + "inner_count": 63, + "inner_radius": 0.186, + "outer_count": 70, + "outer_radius": 0.207, + "ring_position": 0.022 + }, + { + "candidate": 12, + "inner_count": 62, + "inner_radius": 0.183, + "outer_count": 69, + "outer_radius": 0.204, + "ring_position": 0.04400001 + }, + { + "candidate": 11, + "inner_count": 61, + "inner_radius": 0.18, + "outer_count": 68, + "outer_radius": 0.201, + "ring_position": 0.066 + }, + { + "candidate": 10, + "inner_count": 60, + "inner_radius": 0.177, + "outer_count": 67, + "outer_radius": 0.198, + "ring_position": 0.08799998 + }, + { + "candidate": 9, + "inner_count": 59, + "inner_radius": 0.174, + "outer_count": 66, + "outer_radius": 0.195, + "ring_position": 0.11 + }, + { + "candidate": 18, + "inner_count": 68, + "inner_radius": 0.201, + "outer_count": 75, + "outer_radius": 0.222, + "ring_position": 0.132 + }, + { + "candidate": 5, + "inner_count": 55, + "inner_radius": 0.162, + "outer_count": 62, + "outer_radius": 0.183, + "ring_position": 0.154 + }, + { + "candidate": 0, + "inner_count": 50, + "inner_radius": 0.148, + "outer_count": 57, + "outer_radius": 0.169, + "ring_position": 0.176 + }, + { + "candidate": 5, + "inner_count": 55, + "inner_radius": 0.162, + "outer_count": 62, + "outer_radius": 0.183, + "ring_position": 0.198 + }, + { + "candidate": 1, + "inner_count": 51, + "inner_radius": 0.151, + "outer_count": 58, + "outer_radius": 0.172, + "ring_position": 0.22 + }, + { + "candidate": 0, + "inner_count": 50, + "inner_radius": 0.148, + "outer_count": 57, + "outer_radius": 0.169, + "ring_position": 0.242 + } + ], + "dsv": 0.2, + "element": { + "remanence": 1.3, + "size": 0.012 + }, + "field_model": "cuboid", + "generations_run": 100, + "genetic": { + "cx_prob": 0.55, + "elitism": 1, + "gene_mut_prob": 0.05, + "generations": 100, + "mut_prob": 0.4, + "population": 10000, + "seed": 7, + "tournament": 3 + }, + "genome": [ + 11, + 13, + 12, + 11, + 10, + 9, + 18, + 5, + 0, + 5, + 1, + 0 + ], + "magnet_count": 2819, + "objective": "ppm_bx", + "optimized": { + "homogeneity_ppm": 1503.87, + "max_T": 0.05321901, + "mean_T": 0.05317418, + "min_T": 0.05313905, + "peak_to_peak_T": 7.996708e-05, + "points": 4662 + }, + "order": 1, + "outer_count_offset": 7, + "outer_radius_offset": 0.021, + "resolution": 0.005, + "ring_count": 23, + "ring_separation": 0.022, + "symmetry": "octant", + "timing_s": { + "basis": 73.046293053, + "evolution": 23.015295853, + "verification": 26.125628915 + }, + "verified": { + "homogeneity_ppm": 1534.972, + "max_T": 0.05322012, + "mean_T": 0.05317432, + "min_T": 0.0531385, + "peak_to_peak_T": 8.162111e-05, + "points": 33401 + } + } +} diff --git a/optimizer/CMakeLists.txt b/optimizer/CMakeLists.txt new file mode 100644 index 0000000..3c5470a --- /dev/null +++ b/optimizer/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required(VERSION 3.14...3.22) + +project(HalbachOptimizer LANGUAGES CXX) + +# --- Import tools ---- + +include(../cmake/tools.cmake) + +# ---- Dependencies ---- + +include(../cmake/CPM.cmake) + +CPMAddPackage( + GITHUB_REPOSITORY jarro2783/cxxopts + VERSION 3.0.0 + OPTIONS "CXXOPTS_BUILD_EXAMPLES NO" "CXXOPTS_BUILD_TESTS NO" "CXXOPTS_ENABLE_INSTALL YES" +) + +CPMAddPackage(NAME Greeter SOURCE_DIR ${CMAKE_CURRENT_LIST_DIR}/..) + +# ---- Create the optimizer executable ---- + +file(GLOB sources CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/source/*.cpp) + +add_executable(${PROJECT_NAME} ${sources}) + +set_target_properties( + ${PROJECT_NAME} PROPERTIES CXX_STANDARD 17 OUTPUT_NAME "halbach-optimizer" +) + +target_link_libraries(${PROJECT_NAME} Greeter::Greeter cxxopts) +target_link_libraries(${PROJECT_NAME} Kokkos::kokkos) +target_link_libraries(${PROJECT_NAME} nlohmann_json::nlohmann_json) diff --git a/optimizer/source/main.cpp b/optimizer/source/main.cpp new file mode 100644 index 0000000..66e7583 --- /dev/null +++ b/optimizer/source/main.cpp @@ -0,0 +1,272 @@ +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + + +namespace { + +/* + Says how the evolution is going, on one line a generation that rewrites + itself. A hundred generations of ten thousand individuals is a wait, and a + wait with nothing on the screen looks like a hang. +*/ +class ConsoleGenerations: public greeter::optimization::GenerationSink { + + public: + + explicit ConsoleGenerations(const uint32_t& _total): total(_total) {} + + bool onGeneration( + const greeter::optimization::GenerationRecord& record) override { + + std::cout << "\r generation " << record.generation << " of " << total + << ", best " << std::fixed << std::setprecision(1) + << record.best_ever << " ppm " << std::flush; + + if (record.generation == total) { + std::cout << std::endl; + } + + return true; + } + + private: + + uint32_t total; +}; + + +void report(const greeter::optimization::HalbachSolution& solution) { + + const std::vector positions = + solution.spec.getSymmetricRingPositions(); + + std::cout << std::endl << "Chosen rings" << std::endl; + + for (size_t gene = 0; gene < solution.genome.size(); gene++) { + + const greeter::optimization::RingCandidate& candidate = + solution.spec.candidates[solution.genome[gene]]; + + std::cout << " z = " << std::fixed << std::setprecision(4) + << positions[gene] << " m candidate " + << solution.genome[gene] << " inner " + << std::setprecision(3) << candidate.radius << " m of " + << candidate.count << " outer " + << candidate.radius + solution.spec.outer_radius_offset + << " m of " + << candidate.count + solution.spec.outer_count_offset + << std::endl; + } + + std::cout << std::endl + << "Basis: " << solution.basis_configurations + << " configurations, " << solution.basis_magnets + << " magnets sampled at " << solution.optimized_points + << " points" << std::endl + << "Magnets: " << solution.num_magnets << std::endl; + + std::cout << std::setprecision(2); + + std::cout << "Optimized over " << solution.optimized_points + << " points: " << std::fixed << std::setprecision(1) + << solution.optimized.ppm << " ppm, mean " + << 1e3f * solution.optimized.mean << " mT" << std::endl; + + if (solution.was_verified) { + std::cout << "Verified over " << solution.verified_points + << " points: " << solution.verified.ppm << " ppm, mean " + << 1e3f * solution.verified.mean << " mT, peak to peak " + << std::setprecision(4) + << 1e3f * (solution.verified.max - solution.verified.min) + << " mT" << std::endl; + } + + std::cout << std::setprecision(2) + << "Time: basis " << solution.basis_seconds << " s, evolution " + << solution.evolution_seconds << " s"; + + if (solution.was_verified) { + std::cout << ", verification " << solution.verification_seconds << " s"; + } + + std::cout << std::endl; +} + +} // namespace + + +auto main(int argc, char** argv) -> int { + + cxxopts::Options options( + "halbach-optimizer", + "ParaMagneticS, homogeneity optimization of a Halbach cylinder"); + + options.add_options() + ("c,config", "Configuration file. Without one the defaults are those " + "of HalbachOptimisation/homogeneityOptimisation.py", + cxxopts::value()) + ("o,output", "Where to write the optimized magnet, as an input file " + "this project's simulator runs", + cxxopts::value()->default_value("output/halbach_optimized.json")) + ("history", "Also write the convergence curve, one row a generation", + cxxopts::value()) + ("emit", "Write the result as \"arrangements\" or as \"magnets\"", + cxxopts::value()->default_value("arrangements")) + ("field-model", "\"cuboid\" for the analytic kernel of this library, " + "\"dipole\" for the far field the Python script uses", + cxxopts::value()) + ("objective", "\"ppm_bx\" or \"ppm_bmag\"", + cxxopts::value()) + ("symmetry", "Part of the sphere the search is measured over: " + "\"octant\", \"hemisphere\" or \"full\"", + cxxopts::value()) + ("g,generations", "Generations to evolve", cxxopts::value()) + ("p,population", "Individuals in a generation", cxxopts::value()) + ("s,seed", "Seed of the random number pool", cxxopts::value()) + ("no-verify", "Skip measuring the answer over the whole sphere") + ("q,quiet", "Do not announce every step") + ("h,help", "Print this"); + + /* + Kokkos reads its own "--kokkos-..." arguments out of the command line, + and this parser refuses anything it does not know about, which is worth + keeping: a mistyped option is a run that quietly does something else. + So the two are separated rather than the parser being made permissive, + and Kokkos still sees the whole of the original line. + */ + std::vector arguments; + arguments.reserve(argc); + + for (int i = 0; i < argc; i++) { + if (std::string(argv[i]).rfind("--kokkos-", 0) != 0) { + arguments.push_back(argv[i]); + } + } + + int filtered_argc = (int) arguments.size(); + char** filtered_argv = arguments.data(); + + cxxopts::ParseResult parsed; + + try { + parsed = options.parse(filtered_argc, filtered_argv); + } catch (const std::exception& error) { + std::cerr << "Error: " << error.what() << std::endl; + return 1; + } + + if (parsed.count("help")) { + std::cout << options.help() << std::endl; + return 0; + } + + Kokkos::initialize(argc, argv); + + int status = 0; + + // In its own scope: everything holding a Kokkos view has to be gone before + // Kokkos::finalize, and a solution that outlives it takes the process down. + { + try { + + greeter::optimization::HalbachSpec spec; + + if (parsed.count("config")) { + spec = greeter::HalbachOptimizationIO::readSpecFile( + parsed["config"].as()); + } + + // The command line has the last word over the file, so that a + // shorter run of the same configuration does not need a second one. + nlohmann::json overrides = nlohmann::json::object(); + + if (parsed.count("field-model")) { + overrides["field_model"] = parsed["field-model"].as(); + } + if (parsed.count("objective")) { + overrides["objective"] = parsed["objective"].as(); + } + if (parsed.count("symmetry")) { + overrides["symmetry"] = parsed["symmetry"].as(); + } + if (parsed.count("generations")) { + overrides["genetic"]["generations"] = + parsed["generations"].as(); + } + if (parsed.count("population")) { + overrides["genetic"]["population"] = + parsed["population"].as(); + } + if (parsed.count("seed")) { + overrides["genetic"]["seed"] = parsed["seed"].as(); + } + + if (!overrides.empty()) { + spec.readJSON(overrides); + } + + const std::string emit_name = parsed["emit"].as(); + + if (emit_name != "arrangements" && emit_name != "magnets") { + throw std::invalid_argument( + "--emit is \"arrangements\" or \"magnets\", not \"" + + emit_name + "\""); + } + + const greeter::HalbachEmit emit = + emit_name == "magnets" ? greeter::HalbachEmit::Magnets + : greeter::HalbachEmit::Arrangements; + + greeter::optimization::RunOptions run_options; + run_options.verbose = false; + run_options.verify = !parsed.count("no-verify"); + + const bool quiet = parsed.count("quiet") > 0; + + ConsoleGenerations generations(spec.genetic.generations); + + const greeter::optimization::HalbachSolution solution = + greeter::optimization::HalbachOptimizer::run( + spec, quiet ? nullptr : &generations, run_options); + + if (!quiet) { + report(solution); + } + + const std::string output_path = parsed["output"].as(); + + greeter::HalbachOptimizationIO::writeFile(solution, emit, output_path); + + std::cout << "Written to " << output_path << std::endl; + + if (parsed.count("history")) { + + const std::string history_path = + parsed["history"].as(); + + greeter::HalbachOptimizationIO::writeHistoryCSV( + solution, history_path); + + std::cout << "Convergence written to " << history_path + << std::endl; + } + + } catch (const std::exception& error) { + std::cerr << "Error: " << error.what() << std::endl; + status = 1; + } + } + + Kokkos::finalize(); + + return status; +} diff --git a/source/io/CubicMagnetIO.cpp b/source/io/CubicMagnetIO.cpp index 3706986..6e435d6 100644 --- a/source/io/CubicMagnetIO.cpp +++ b/source/io/CubicMagnetIO.cpp @@ -31,14 +31,9 @@ std::vector greeter::CubicMagnetIO::readMagnetization(const nlohmann::jso } std::unique_ptr greeter::CubicMagnetIO::createMagnet(const nlohmann::json& magnet) { - std::cout << "Reading position" << std::endl; std::vector position = readPosition(magnet); - std::cout << "position read" << std::endl; std::vector dimensions = readDimensions(magnet); - std::cout << "dimensions read" << std::endl; std::vector orientation = readOrientation(magnet); - std::cout << "orientation read" << std::endl; std::vector magnetization = readMagnetization(magnet); - std::cout << "extracted all data from JSON" << std::endl; return std::make_unique(position, dimensions, orientation, magnetization); } diff --git a/source/io/HalbachOptimizationIO.cpp b/source/io/HalbachOptimizationIO.cpp new file mode 100644 index 0000000..c023610 --- /dev/null +++ b/source/io/HalbachOptimizationIO.cpp @@ -0,0 +1,268 @@ +#include + +#include +#include + +#include +#include +#include + + +namespace greeter { + +namespace { + +using greeter::optimization::HalbachSolution; +using greeter::optimization::HalbachSpec; +using greeter::optimization::FieldMetrics; +using greeter::optimization::tidy; + +/* The three numbers a vector of a JSON file is written as. */ +nlohmann::json vector3(const float* values) { + return nlohmann::json::array({tidy(values[0]), tidy(values[1]), tidy(values[2])}); +} + +/* + The box the optimization measured in, as a "field_of_view". + + The written file is meant to be runnable, and a file with no field of view + and no force section is not one this library reads at all. The grid is the + one HalbachBasis samples, so a run of the written file lands on the points + the answer was measured at. +*/ +nlohmann::json fieldOfView(const HalbachSpec& spec) { + + const int64_t steps = std::lround(spec.dsv / spec.resolution); + + const float half = 0.5f * spec.dsv; + + const nlohmann::json axis = { + {"min", tidy(-half)}, {"max", tidy(half)}, {"n", steps + 1} + }; + + return nlohmann::json{{"x", axis}, {"y", axis}, {"z", axis}}; +} + +nlohmann::json metrics(const FieldMetrics& measured, const size_t& points) { + return nlohmann::json{ + {"homogeneity_ppm", tidy(measured.ppm)}, + {"mean_T", tidy(measured.mean)}, + {"min_T", tidy(measured.min)}, + {"max_T", tidy(measured.max)}, + {"peak_to_peak_T", tidy(measured.max - measured.min)}, + {"points", points} + }; +} + +} // namespace + + +HalbachOptimizationIO::HalbachOptimizationIO() {} + +HalbachOptimizationIO::~HalbachOptimizationIO() {} + + +greeter::optimization::HalbachSpec HalbachOptimizationIO::readSpec( + const nlohmann::json& data) { + + greeter::optimization::HalbachSpec spec; + + if (data.contains("halbach_optimization")) { + spec.readJSON(data["halbach_optimization"]); + } else { + spec.readJSON(data); + } + + return spec; +} + + +greeter::optimization::HalbachSpec HalbachOptimizationIO::readSpecFile( + const std::string& path) { + + std::ifstream file(path); + + if (!file.is_open()) { + throw std::invalid_argument("Could not open " + path); + } + + return readSpec(nlohmann::json::parse(file)); +} + + +nlohmann::json HalbachOptimizationIO::writeMagnet( + const greeter::MagnetCollection& collection, const size_t& index, + const int64_t& id) { + + const std::vector parameters = collection.getMagnetParameters(index); + const uint16_t type = collection.getMagnetTypeID(index); + + // position (3), orientation (4), the geometry of the shape, magnetization + // (3). See MagnetParameters.h, which is the one place that layout is set. + const float* position = parameters.data(); + const float* orientation = parameters.data() + 3; + const float* dimensions = parameters.data() + 7; + + if (type == greeter::CuboidMagnet::getStaticTypeID()) { + + return nlohmann::json{ + {"id", id}, + {"type", greeter::CuboidMagnet::getStaticTypeName()}, + {"parameters", { + {"dimensions", vector3(dimensions)}, + {"magnetization", vector3(dimensions + 3)}, + {"position", vector3(position)}, + {"orientation", nlohmann::json::array({ + tidy(orientation[0]), tidy(orientation[1]), + tidy(orientation[2]), tidy(orientation[3]) + })} + }} + }; + } + + if (type == greeter::DipoleMagnet::getStaticTypeID()) { + + // A dipole has no geometry at all, so its moment sits where the + // dimensions of a shaped magnet would start. + return nlohmann::json{ + {"id", id}, + {"type", greeter::DipoleMagnet::getStaticTypeName()}, + {"parameters", { + {"moment", vector3(dimensions)}, + {"position", vector3(position)}, + {"orientation", nlohmann::json::array({ + tidy(orientation[0]), tidy(orientation[1]), + tidy(orientation[2]), tidy(orientation[3]) + })} + }} + }; + } + + throw std::invalid_argument( + "A Halbach optimization writes out cuboids and dipoles, not the magnet " + "type " + std::to_string(type)); +} + + +nlohmann::json HalbachOptimizationIO::write( + const HalbachSolution& solution, const HalbachEmit& emit) { + + nlohmann::json record = solution.spec.toJSON(); + + nlohmann::json genome = nlohmann::json::array(); + for (const auto& gene : solution.genome) { + genome.push_back(gene); + } + + nlohmann::json chosen = nlohmann::json::array(); + + const std::vector positions = solution.spec.getSymmetricRingPositions(); + const std::vector candidates = + solution.getChosenCandidates(); + + for (size_t gene = 0; gene < candidates.size(); gene++) { + chosen.push_back({ + {"ring_position", tidy(positions[gene])}, + {"candidate", solution.genome[gene]}, + {"inner_radius", tidy(candidates[gene].radius)}, + {"inner_count", candidates[gene].count}, + {"outer_radius", + tidy(candidates[gene].radius + solution.spec.outer_radius_offset)}, + {"outer_count", + candidates[gene].count + solution.spec.outer_count_offset} + }); + } + + record["genome"] = genome; + record["chosen"] = chosen; + record["magnet_count"] = solution.num_magnets; + + record["basis"] = { + {"configurations", solution.basis_configurations}, + {"magnets", solution.basis_magnets}, + {"points", solution.optimized_points} + }; + + record["optimized"] = metrics(solution.optimized, solution.optimized_points); + + if (solution.was_verified) { + record["verified"] = metrics(solution.verified, solution.verified_points); + } + + record["timing_s"] = { + {"basis", solution.basis_seconds}, + {"evolution", solution.evolution_seconds}, + {"verification", solution.verification_seconds} + }; + + record["generations_run"] = + solution.history.empty() ? 0 : solution.history.back().generation; + + nlohmann::json data; + + // First, so that a person opening the file reads what was found before + // the two thousand magnets it was found in. + data["halbach_optimization"] = record; + + if (emit == HalbachEmit::Magnets) { + + const greeter::MagnetCollection collection = solution.buildCollection(); + + nlohmann::json magnets = nlohmann::json::array(); + + for (size_t i = 0; i < collection.get_num_magnets(); i++) { + magnets.push_back(writeMagnet(collection, i, (int64_t) i + 1)); + } + + data["magnets"] = magnets; + + } else { + + nlohmann::json arrangements = nlohmann::json::array(); + + for (const auto& ring : solution.buildArrangements()) { + arrangements.push_back(ring); + } + + data["arrangements"] = arrangements; + } + + data["field_of_view"] = fieldOfView(solution.spec); + + return data; +} + + +void HalbachOptimizationIO::writeFile( + const HalbachSolution& solution, const HalbachEmit& emit, + const std::string& path) { + + std::ofstream file(path); + + if (!file.is_open()) { + throw std::invalid_argument("Could not open " + path + " for writing"); + } + + file << write(solution, emit).dump(1, ' ') << std::endl; +} + + +void HalbachOptimizationIO::writeHistoryCSV( + const HalbachSolution& solution, const std::string& path) { + + std::ofstream file(path); + + if (!file.is_open()) { + throw std::invalid_argument("Could not open " + path + " for writing"); + } + + file << "generation,best_ppm,best_ever_ppm,mean_ppm,seconds" << std::endl; + + for (const auto& record : solution.history) { + file << record.generation << "," << record.best << "," + << record.best_ever << "," << record.mean << "," + << record.seconds << std::endl; + } +} + +} // namespace greeter diff --git a/source/io/SphericalMagnetIO.cpp b/source/io/SphericalMagnetIO.cpp index 9a79706..253cd1f 100644 --- a/source/io/SphericalMagnetIO.cpp +++ b/source/io/SphericalMagnetIO.cpp @@ -83,15 +83,10 @@ std::vector greeter::SphericalMagnetIO::readMagnetization(const nlohmann: } std::unique_ptr greeter::SphericalMagnetIO::createMagnet(const nlohmann::json& magnet) { - std::cout << "Reading position" << std::endl; std::vector position = readPosition(magnet); - std::cout << "position read" << std::endl; std::vector dimensions = readDimensions(magnet); - std::cout << "dimensions read" << std::endl; std::vector orientation = readOrientation(magnet); - std::cout << "orientation read" << std::endl; std::vector magnetization = readMagnetization(magnet); - std::cout << "extracted all data from JSON" << std::endl; return std::make_unique( position, orientation, dimensions[0], magnetization[2]); } diff --git a/source/io/TetrahedronMagnetIO.cpp b/source/io/TetrahedronMagnetIO.cpp index 1549a71..b3efee1 100644 --- a/source/io/TetrahedronMagnetIO.cpp +++ b/source/io/TetrahedronMagnetIO.cpp @@ -85,15 +85,10 @@ std::vector greeter::TetrahedronMagnetIO::readMagnetization(const nlohman std::unique_ptr greeter::TetrahedronMagnetIO::createMagnet( const nlohmann::json& magnet) { - std::cout << "Reading position" << std::endl; std::vector position = readPosition(magnet); - std::cout << "position read" << std::endl; std::vector vertices = readVertices(magnet); - std::cout << "vertices read" << std::endl; std::vector orientation = readOrientation(magnet); - std::cout << "orientation read" << std::endl; std::vector magnetization = readMagnetization(magnet); - std::cout << "extracted all data from JSON" << std::endl; return std::make_unique( position, vertices, orientation, magnetization); diff --git a/source/optimization/GeneticOptimizer.cpp b/source/optimization/GeneticOptimizer.cpp new file mode 100644 index 0000000..0b4d20d --- /dev/null +++ b/source/optimization/GeneticOptimizer.cpp @@ -0,0 +1,394 @@ +#include + +#include +#include +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +namespace { + +typedef RandPoolType::generator_type Generator; + +/* + The best of a few individuals drawn at random. + + Tournament selection, at the size the Python script uses. It is the whole + of the selection pressure here: a larger tournament converges sooner and + explores less. +*/ +KOKKOS_INLINE_FUNCTION +uint32_t tournament(Generator& generator, const FitnessView& fitness, + const uint32_t& population_size, const uint32_t& size) { + + uint32_t winner = Kokkos::rand::draw( + generator, population_size); + + for (uint32_t i = 1; i < size; i++) { + + const uint32_t challenger = Kokkos::rand::draw( + generator, population_size); + + if (fitness(challenger) < fitness(winner)) { + winner = challenger; + } + } + + return winner; +} + +} // namespace + + +GeneticOptimizer::GeneticOptimizer(const GeneticSettings& _settings, + const size_t& _num_genes, + const size_t& _num_candidates): + settings(_settings), + num_genes(_num_genes), + num_candidates(_num_candidates), + random_pool(_settings.seed), + best_fitness(std::numeric_limits::max()) { + + if (num_genes == 0) { + throw std::invalid_argument("A genome needs at least one gene"); + } + + if (num_candidates < 1) { + throw std::invalid_argument("A gene needs at least one candidate value"); + } + + if (settings.population < 2) { + throw std::invalid_argument("A population needs at least two individuals"); + } + + if (settings.elitism >= settings.population) { + throw std::invalid_argument( + "The elitism carries over the whole population"); + } +} + + +GeneticOptimizer::~GeneticOptimizer() {} + + +void GeneticOptimizer::reproduce(const PopulationView& from, + const PopulationView& into, + const FitnessView& fitness) { + + const uint32_t population_size = settings.population; + const uint32_t the_num_genes = (uint32_t) num_genes; + const uint32_t alphabet = (uint32_t) num_candidates; + + const uint32_t tournament_size = settings.tournament; + + const float cx_prob = settings.cx_prob; + const float mut_prob = settings.mut_prob; + const float gene_mut_prob = settings.gene_mut_prob; + + RandPoolType pool = random_pool; + + // Offspring are made two at a time because a crossover needs two of them, + // which is also what lets the whole generation be one parallel region: a + // pair is the unit nothing outside of it touches. + const uint32_t pairs = (population_size + 1) / 2; + + Kokkos::parallel_for( + "reproduce", Kokkos::RangePolicy(0, pairs), + KOKKOS_LAMBDA(const uint32_t pair) { + + Generator generator = pool.get_state(); + + const uint32_t first = 2 * pair; + const uint32_t second = first + 1; + const bool has_second = second < population_size; + + const uint32_t parent_a = + tournament(generator, fitness, population_size, tournament_size); + const uint32_t parent_b = + tournament(generator, fitness, population_size, tournament_size); + + for (uint32_t gene = 0; gene < the_num_genes; gene++) { + into(first, gene) = from(parent_a, gene); + if (has_second) { + into(second, gene) = from(parent_b, gene); + } + } + + // Two point crossover: the stretch between two cut points changes + // hands. With one gene there is nothing to cut. + if (has_second && the_num_genes > 1 + && generator.frand() < cx_prob) { + + uint32_t cut_a = Kokkos::rand::draw( + generator, the_num_genes); + uint32_t cut_b = Kokkos::rand::draw( + generator, the_num_genes); + + if (cut_a > cut_b) { + const uint32_t swap = cut_a; + cut_a = cut_b; + cut_b = swap; + } + + for (uint32_t gene = cut_a; gene <= cut_b; gene++) { + const uint16_t swap = into(first, gene); + into(first, gene) = into(second, gene); + into(second, gene) = swap; + } + } + + /* + Mutation, in two stages as DEAP applies it: an individual is + picked at all with one probability, and then each of its genes + moves with another. A gene that moves is redrawn from the + alphabet, see the note on mutFlipBit in the header. + */ + for (uint32_t child = 0; child < 2; child++) { + + if (child == 1 && !has_second) { + break; + } + + if (generator.frand() >= mut_prob) { + continue; + } + + const uint32_t row = child == 0 ? first : second; + + for (uint32_t gene = 0; gene < the_num_genes; gene++) { + if (generator.frand() < gene_mut_prob) { + into(row, gene) = (uint16_t) + Kokkos::rand::draw( + generator, alphabet); + } + } + } + + pool.free_state(generator); + }); + + Kokkos::fence(); +} + + +std::vector GeneticOptimizer::run( + const HomogeneityObjective& objective, GenerationSink* sink, + const bool& verbose) { + + if (objective.getNumGenes() != num_genes + || objective.getNumCandidates() != num_candidates) { + throw std::invalid_argument( + "The objective and the optimizer disagree about the shape of a genome"); + } + + history.clear(); + stopped = false; + best_fitness = std::numeric_limits::max(); + best_genome.assign(num_genes, 0); + + PopulationView population("population", settings.population, num_genes); + PopulationView offspring("offspring", settings.population, num_genes); + + FitnessView fitness("fitness", settings.population); + + // A first generation drawn uniformly from the alphabet, which is what + // toolbox.attr_bool does in the script this ports. + { + RandPoolType pool = random_pool; + const uint32_t alphabet = (uint32_t) num_candidates; + const uint32_t the_num_genes = (uint32_t) num_genes; + + Kokkos::parallel_for( + "seed_population", + Kokkos::RangePolicy(0, settings.population), + KOKKOS_LAMBDA(const uint32_t individual) { + + Generator generator = pool.get_state(); + + for (uint32_t gene = 0; gene < the_num_genes; gene++) { + population(individual, gene) = (uint16_t) + Kokkos::rand::draw(generator, alphabet); + } + + pool.free_state(generator); + }); + + Kokkos::fence(); + } + + objective.evaluate(population, fitness); + + typename FitnessView::HostMirror host_fitness = + Kokkos::create_mirror_view(fitness); + + /* + A host buffer of its own rather than a mirror of `population`. + + The two population buffers are swapped every generation, and a mirror + would go on naming whichever one it was made from: on a host build it + *is* that buffer, so after a swap the copy below would write the new + offspring over the buffer they are being made in. Re-making the mirror + each generation would work and would allocate once a generation. + */ + Kokkos::View host_population( + "host_population", settings.population, num_genes); + + // Where the elite of a generation are, kept on the device so that copying + // them over is a parallel region rather than a row at a time. + Kokkos::View elite_rows( + "elite_rows", settings.elitism > 0 ? settings.elitism : 1); + + typename Kokkos::View::HostMirror + host_elite_rows = Kokkos::create_mirror_view(elite_rows); + + std::vector order(settings.population); + + Kokkos::Timer timer; + + // The clock is read at the top of each pass, so what a record holds is the + // time from the start of the previous pass: the making and the scoring of + // the generation it describes. Timing only the bookkeeping below it would + // report tens of microseconds for work that took a second. + double previous = 0.0; + + for (uint32_t generation = 0; generation <= settings.generations; + generation++) { + + const double now = timer.seconds(); + const double elapsed = now - previous; + previous = now; + + Kokkos::deep_copy(host_fitness, fitness); + + /* + The generation is summed up on the host. It is a few thousand floats + against a parallel region that has just read tens of millions, so + the copy costs nothing, and it buys the mean, the best and the + ranking the elite are taken from without three more reductions. + */ + uint32_t best_index = 0; + double total = 0.0; + + for (uint32_t i = 0; i < settings.population; i++) { + total += (double) host_fitness(i); + if (host_fitness(i) < host_fitness(best_index)) { + best_index = i; + } + } + + const float generation_best = host_fitness(best_index); + + if (generation_best < best_fitness) { + + best_fitness = generation_best; + + Kokkos::deep_copy(host_population, population); + + for (size_t gene = 0; gene < num_genes; gene++) { + best_genome[gene] = host_population(best_index, gene); + } + } + + GenerationRecord record; + record.generation = generation; + record.best = generation_best; + record.best_ever = best_fitness; + record.mean = (float) (total / (double) settings.population); + record.seconds = (float) elapsed; + + history.push_back(record); + + if (verbose) { + std::cout << " generation " << generation << ": best " + << generation_best << " ppm, best ever " << best_fitness + << " ppm" << std::endl; + } + + if (sink != nullptr && !sink->onGeneration(record)) { + stopped = true; + break; + } + + // The last pass through is a report on the final population, not a + // generation of its own. + if (generation == settings.generations) { + break; + } + + reproduce(population, offspring, fitness); + + /* + The best few of the outgoing generation are put back untouched. + Without this a crossover can throw away the best genome found so + far, and the population is free to drift back uphill; the Python + script keeps its best in a variable on the side and lets exactly + that happen. + */ + if (settings.elitism > 0) { + + std::iota(order.begin(), order.end(), 0u); + + const size_t keep = + std::min((size_t) settings.elitism, order.size()); + + std::partial_sort( + order.begin(), order.begin() + keep, order.end(), + [&](const uint32_t& a, const uint32_t& b) { + return host_fitness(a) < host_fitness(b); + }); + + for (size_t i = 0; i < keep; i++) { + host_elite_rows(i) = order[i]; + } + + Kokkos::deep_copy(elite_rows, host_elite_rows); + + const uint32_t the_num_genes = (uint32_t) num_genes; + const PopulationView from = population; + const PopulationView into = offspring; + const Kokkos::View rows = elite_rows; + + Kokkos::parallel_for( + "carry_elite", + Kokkos::RangePolicy(0, (uint32_t) keep), + KOKKOS_LAMBDA(const uint32_t i) { + for (uint32_t gene = 0; gene < the_num_genes; gene++) { + into(i, gene) = from(rows(i), gene); + } + }); + + Kokkos::fence(); + } + + // The offspring become the population. Two buffers swapped rather + // than one copied, so a generation costs no traffic beyond making it. + PopulationView swap = population; + population = offspring; + offspring = swap; + + objective.evaluate(population, fitness); + } + + return best_genome; +} + + +float GeneticOptimizer::getBestFitness() const { + return best_fitness; +} + +const std::vector& GeneticOptimizer::getHistory() const { + return history; +} + +bool GeneticOptimizer::wasStopped() const { + return stopped; +} + +} // namespace optimization +} // namespace greeter diff --git a/source/optimization/HalbachBasis.cpp b/source/optimization/HalbachBasis.cpp new file mode 100644 index 0000000..9af2a16 --- /dev/null +++ b/source/optimization/HalbachBasis.cpp @@ -0,0 +1,304 @@ +#include + +#include +#include +#include +#include + +#include +#include +#include + + +namespace greeter { +namespace optimization { + +HalbachBasis::HalbachBasis() {} + +HalbachBasis::~HalbachBasis() {} + + +std::vector> makeSamplePoints(const HalbachSpec& spec) { + + // The same grid numpy.linspace gives for the box of the DSV: the number of + // steps that fit in it, and one more for the far edge. + const int32_t steps = (int32_t) std::lround(spec.dsv / spec.resolution); + + if (steps < 1) { + throw std::invalid_argument( + "The \"resolution\" does not divide the \"dsv\" into any step"); + } + + const int32_t n = steps + 1; + + const float radius = 0.5f * spec.dsv; + const float spacing = spec.dsv / (float) steps; + + // A point exactly on the surface of the sphere is inside it, and one a + // rounding error outside would otherwise drop out of an octant along its + // own axis. The tolerance is a thousandth of a step. + const float tolerance = 1e-3f * spacing; + const float radius_squared = radius * radius + tolerance * radius; + + std::vector> points; + + for (int32_t i = 0; i < n; i++) { + + const float x = -radius + (float) i * spacing; + + if (spec.symmetry == Symmetry::Octant && x < -tolerance) { + continue; + } + + for (int32_t j = 0; j < n; j++) { + + const float y = -radius + (float) j * spacing; + + if (spec.symmetry == Symmetry::Octant && y < -tolerance) { + continue; + } + + for (int32_t k = 0; k < n; k++) { + + const float z = -radius + (float) k * spacing; + + if (spec.symmetry != Symmetry::Full && z < -tolerance) { + continue; + } + + if (x * x + y * y + z * z > radius_squared) { + continue; + } + + points.push_back({x, y, z}); + } + } + } + + if (points.empty()) { + throw std::invalid_argument( + "The sampled volume came out empty, which means the \"resolution\" " + "is too coarse for the \"dsv\""); + } + + return points; +} + + +HalbachBasis HalbachBasis::build(const HalbachSpec& spec, const bool& verbose) { + + spec.validate(); + + HalbachBasis basis; + + basis.points = makeSamplePoints(spec); + basis.num_genes = spec.getNumGenes(); + basis.num_candidates = spec.getNumCandidates(); + + const size_t num_points = basis.points.size(); + const size_t num_configurations = basis.num_genes * basis.num_candidates; + + /* + Every configuration's magnets, one after the other in a single + collection, with a note of where each configuration's run begins. The + rings themselves are laid out by HalbachRingArrangement, which is what + an input file asking for a "halbach_ring" goes through, so the magnet + measured here is the magnet the solution writes out. + */ + greeter::MagnetCollection all; + + std::vector configuration_first(num_configurations, 0); + std::vector configuration_count(num_configurations, 0); + + for (size_t gene = 0; gene < basis.num_genes; gene++) { + for (size_t candidate = 0; candidate < basis.num_candidates; candidate++) { + + const size_t configuration = gene * basis.num_candidates + candidate; + + configuration_first[configuration] = all.get_num_magnets(); + + for (const auto& ring : spec.makeRingsForGene(gene, candidate)) { + for (auto& member : + greeter::HalbachRingArrangement::expand(ring)) { + all.addMagnet(std::move(member)); + } + } + + configuration_count[configuration] = + all.get_num_magnets() - configuration_first[configuration]; + } + } + + basis.num_magnets = all.get_num_magnets(); + + if (verbose) { + std::cout << " " << num_configurations << " configurations, " + << basis.num_magnets << " magnets, " << num_points + << " observation points" << std::endl; + } + + // The parameters of every magnet end to end, in the order the field + // kernels read them. One array for all of them, as the simulators use. + UInt32VectorView magnet_types("magnet_types", basis.num_magnets); + UInt32VectorView parameter_offsets("parameter_offsets", basis.num_magnets); + FloatVectorView magnet_parameters( + "magnet_parameters", all.getTotalNumOfParameters()); + + all.fillMagnetParameters(magnet_parameters, parameter_offsets, magnet_types); + + /* + The kernel of every magnet, looked up once. Doing it inside the loop + would put a hash lookup and a branch on the type in front of every one + of the hundreds of millions of field evaluations below. This is what + MagneticFieldSimulator::resolveMagnetTypes does, for the same reason. + */ + MagnetKernelView kernels("magnet_kernels", basis.num_magnets); + + const greeter::MagneticFieldMethodFactory& factory = + greeter::MagneticFieldMethodFactory::getInstance(); + + for (size_t i = 0; i < basis.num_magnets; i++) { + kernels(i).kernel = + factory.getComputeMagneticField((u_int16_t) magnet_types(i)); + kernels(i).polarization = nullptr; + kernels(i).parameter_offset = parameter_offsets(i); + } + + Float3VectorView observation_points("observation_points", num_points); + + for (size_t i = 0; i < num_points; i++) { + observation_points(i, 0) = basis.points[i][0]; + observation_points(i, 1) = basis.points[i][1]; + observation_points(i, 2) = basis.points[i][2]; + } + + UInt32VectorView first("configuration_first", num_configurations); + UInt32VectorView count("configuration_count", num_configurations); + + for (size_t i = 0; i < num_configurations; i++) { + first(i) = configuration_first[i]; + count(i) = configuration_count[i]; + } + + basis.fields = Kokkos::View( + "basis_fields", num_configurations, 3, num_points); + + Kokkos::View fields = basis.fields; + + Kokkos::Timer timer; + + /* + One parallel region for the whole precomputation, a configuration and an + observation point wide. Every iteration is independent and writes three + floats nobody else writes, so there is no reduction and no atomic here. + */ + Kokkos::parallel_for( + "halbach_basis", + Kokkos::MDRangePolicy>( + {0, 0}, {(int64_t) num_configurations, (int64_t) num_points}), + KOKKOS_LAMBDA(const int64_t configuration, const int64_t point) { + + const float observation_point[3] = { + observation_points(point, 0), + observation_points(point, 1), + observation_points(point, 2) + }; + + float sum_x = 0.0f; + float sum_y = 0.0f; + float sum_z = 0.0f; + + const uint32_t begin = first(configuration); + const uint32_t end = begin + count(configuration); + + for (uint32_t m = begin; m < end; m++) { + + const MagnetKernel magnet = kernels(m); + + const float* parameters = greeter::magnetParameters( + magnet_parameters, magnet.parameter_offset); + + float b_x = 0.0f; + float b_y = 0.0f; + float b_z = 0.0f; + + magnet.kernel(parameters, observation_point, b_x, b_y, b_z); + + sum_x += b_x; + sum_y += b_y; + sum_z += b_z; + } + + fields(configuration, 0, point) = sum_x; + fields(configuration, 1, point) = sum_y; + fields(configuration, 2, point) = sum_z; + }); + + Kokkos::fence(); + + if (verbose) { + std::cout << " basis built in " << timer.seconds() << " s" << std::endl; + } + + return basis; +} + + +const std::vector>& HalbachBasis::getPoints() const { + return points; +} + +size_t HalbachBasis::getNumPoints() const { + return points.size(); +} + +size_t HalbachBasis::getNumGenes() const { + return num_genes; +} + +size_t HalbachBasis::getNumCandidates() const { + return num_candidates; +} + +size_t HalbachBasis::getNumMagnets() const { + return num_magnets; +} + +Kokkos::View HalbachBasis::getFields() const { + return fields; +} + + +void HalbachBasis::evaluateAt(const std::vector& genome, + const size_t& point, float* b) const { + + if (genome.size() != num_genes) { + throw std::invalid_argument( + "The genome has " + std::to_string(genome.size()) + + " genes, the basis was built for " + std::to_string(num_genes)); + } + + if (point >= points.size()) { + throw std::out_of_range("Observation point index out of range"); + } + + b[0] = 0.0f; + b[1] = 0.0f; + b[2] = 0.0f; + + for (size_t gene = 0; gene < num_genes; gene++) { + + if (genome[gene] >= num_candidates) { + throw std::out_of_range("Gene value out of range"); + } + + const size_t configuration = gene * num_candidates + genome[gene]; + + for (size_t component = 0; component < 3; component++) { + b[component] += fields(configuration, component, point); + } + } +} + +} // namespace optimization +} // namespace greeter diff --git a/source/optimization/HalbachOptimizer.cpp b/source/optimization/HalbachOptimizer.cpp new file mode 100644 index 0000000..4a1ce71 --- /dev/null +++ b/source/optimization/HalbachOptimizer.cpp @@ -0,0 +1,130 @@ +#include + +#include + +#include + +#include +#include + + +namespace greeter { +namespace optimization { + +FieldMetrics HalbachOptimizer::measure( + const greeter::MagnetCollection& collection, const HalbachSpec& spec, + const Objective& objective, size_t& num_points) { + + /* + Over the whole sphere, whatever symmetry the search was reduced to. An + eighth of a sphere is a fair place to search, because the choice that is + best there is very nearly the choice that is best everywhere, but it is + not a fair place to report from: the ring has a finite number of magnets + and so is only nearly symmetric under the rotations the reduction + assumes. + */ + HalbachSpec full = spec; + full.symmetry = Symmetry::Full; + + const std::vector> points = makeSamplePoints(full); + + num_points = points.size(); + + // The simulator is driven directly rather than through + // MagnetCollection::simulate, which announces every run it makes. This is + // one step of several and has its own line to print. + std::unique_ptr simulator = + collection.createSimulator(); + + simulator->fillObservationPoints(points); + simulator->simulate(false); + + const std::vector fields = simulator->getMagneticFieldsFlat(); + + std::vector values; + values.reserve(num_points); + + for (size_t i = 0; i < num_points; i++) { + + const float b_x = fields[3 * i + 0]; + const float b_y = fields[3 * i + 1]; + const float b_z = fields[3 * i + 2]; + + values.push_back( + objective == Objective::BMagnitude + ? std::sqrt(b_x * b_x + b_y * b_y + b_z * b_z) + : b_x); + } + + return HomogeneityObjective::summarise(values); +} + + +HalbachSolution HalbachOptimizer::run(const HalbachSpec& spec, + GenerationSink* sink, + const RunOptions& options) { + + spec.validate(); + + HalbachSolution solution; + solution.spec = spec; + + Kokkos::Timer timer; + + if (options.verbose) { + std::cout << "Sampling every ring candidate at every ring position" + << std::endl; + } + + const HalbachBasis basis = HalbachBasis::build(spec, options.verbose); + + solution.basis_seconds = timer.seconds(); + solution.optimized_points = basis.getNumPoints(); + solution.basis_magnets = basis.getNumMagnets(); + solution.basis_configurations = basis.getNumGenes() * basis.getNumCandidates(); + + const HomogeneityObjective objective(basis, spec.objective); + + GeneticOptimizer optimizer( + spec.genetic, basis.getNumGenes(), basis.getNumCandidates()); + + if (options.verbose) { + std::cout << "Evolving " << spec.genetic.population << " individuals of " + << basis.getNumGenes() << " genes over " + << spec.genetic.generations << " generations" << std::endl; + } + + timer.reset(); + + solution.genome = optimizer.run(objective, sink, options.verbose); + + solution.evolution_seconds = timer.seconds(); + solution.history = optimizer.getHistory(); + + solution.optimized = objective.evaluateGenome(solution.genome); + + const greeter::MagnetCollection collection = solution.buildCollection(); + + solution.num_magnets = collection.get_num_magnets(); + + if (options.verify) { + + if (options.verbose) { + std::cout << "Measuring the chosen magnet over the whole sphere" + << std::endl; + } + + timer.reset(); + + solution.verified = measure( + collection, spec, spec.objective, solution.verified_points); + + solution.verification_seconds = timer.seconds(); + solution.was_verified = true; + } + + return solution; +} + +} // namespace optimization +} // namespace greeter diff --git a/source/optimization/HalbachSolution.cpp b/source/optimization/HalbachSolution.cpp new file mode 100644 index 0000000..39e53f0 --- /dev/null +++ b/source/optimization/HalbachSolution.cpp @@ -0,0 +1,70 @@ +#include + +#include + +#include + + +namespace greeter { +namespace optimization { + +std::vector HalbachSolution::buildArrangements() const { + + if (genome.size() != spec.getNumGenes()) { + throw std::invalid_argument( + "The genome has " + std::to_string(genome.size()) + + " genes, the specification asks for " + + std::to_string(spec.getNumGenes())); + } + + std::vector arrangements; + + int64_t id = 1; + + for (size_t gene = 0; gene < genome.size(); gene++) { + + for (auto& ring : spec.makeRingsForGene(gene, genome[gene])) { + ring["id"] = id++; + arrangements.push_back(std::move(ring)); + } + } + + return arrangements; +} + + +greeter::MagnetCollection HalbachSolution::buildCollection() const { + + greeter::MagnetCollection collection; + + // Through the same arrangement an input file goes through, so that the + // magnet checked here is the magnet the written file describes. + for (const auto& ring : buildArrangements()) { + for (auto& member : greeter::HalbachRingArrangement::expand(ring)) { + collection.addMagnet(std::move(member)); + } + } + + return collection; +} + + +std::vector HalbachSolution::getChosenCandidates() const { + + std::vector chosen; + chosen.reserve(genome.size()); + + for (const auto& gene : genome) { + + if (gene >= spec.candidates.size()) { + throw std::out_of_range("Gene value out of range"); + } + + chosen.push_back(spec.candidates[gene]); + } + + return chosen; +} + +} // namespace optimization +} // namespace greeter diff --git a/source/optimization/HalbachSpec.cpp b/source/optimization/HalbachSpec.cpp new file mode 100644 index 0000000..98e8eee --- /dev/null +++ b/source/optimization/HalbachSpec.cpp @@ -0,0 +1,557 @@ +#include + +#include +#include +#include +#include +#include +#include + + +namespace greeter { +namespace optimization { + +namespace { + +constexpr double MU_0 = 4.0e-7 * M_PI; + +/* + A number that has to be there and has to be positive. Reading a length as + zero and carrying on gives a ring of no radius rather than a message. +*/ +float readPositive(const nlohmann::json& data, const std::string& name, + const float& fallback) { + + if (!data.contains(name)) { + return fallback; + } + + if (!data[name].is_number()) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization must be a number"); + } + + const float value = data[name].get(); + + if (!(value > 0.0f)) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization must be positive"); + } + + return value; +} + +float readNumber(const nlohmann::json& data, const std::string& name, + const float& fallback) { + + if (!data.contains(name)) { + return fallback; + } + + if (!data[name].is_number()) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization must be a number"); + } + + return data[name].get(); +} + +uint32_t readCount(const nlohmann::json& data, const std::string& name, + const uint32_t& fallback, const uint32_t& least) { + + if (!data.contains(name)) { + return fallback; + } + + if (!data[name].is_number_integer() && !data[name].is_number_unsigned()) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization must be a whole number"); + } + + const int64_t value = data[name].get(); + + if (value < (int64_t) least) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization must be at least " + + std::to_string(least)); + } + + return (uint32_t) value; +} + +/* A probability, which is a number between zero and one inclusive. */ +float readProbability(const nlohmann::json& data, const std::string& name, + const float& fallback) { + + const float value = readNumber(data, name, fallback); + + if (value < 0.0f || value > 1.0f) { + throw std::invalid_argument( + "\"" + name + "\" of a Halbach optimization is a probability, " + "which lies between zero and one"); + } + + return value; +} + +/* + One of a fixed set of names. Anything else is refused with the list, rather + than quietly falling back to the default and optimizing something the file + did not ask for. +*/ +std::string readChoice(const nlohmann::json& data, const std::string& name, + const std::vector& allowed, + const std::string& fallback) { + + if (!data.contains(name)) { + return fallback; + } + + if (!data[name].is_string()) { + throw std::invalid_argument("\"" + name + "\" must be a string"); + } + + const std::string value = data[name].get(); + + for (const auto& candidate : allowed) { + if (value == candidate) { + return value; + } + } + + std::string message = "\"" + name + "\" is one of"; + for (const auto& candidate : allowed) { + message += " \"" + candidate + "\""; + } + message += ", not \"" + value + "\""; + + throw std::invalid_argument(message); +} + +} // namespace + + +double tidy(const float& value) { + + if (!std::isfinite(value)) { + return (double) value; + } + + char buffer[32]; + std::snprintf(buffer, sizeof(buffer), "%.7g", (double) value); + + return std::strtod(buffer, nullptr); +} + + +const char* toString(const Objective& objective) { + return objective == Objective::Bx ? "ppm_bx" : "ppm_bmag"; +} + +const char* toString(const Symmetry& symmetry) { + switch (symmetry) { + case Symmetry::Octant: return "octant"; + case Symmetry::Hemisphere: return "hemisphere"; + default: return "full"; + } +} + +const char* toString(const FieldModel& model) { + return model == FieldModel::Cuboid ? "cuboid" : "dipole"; +} + + +HalbachSpec::HalbachSpec() { + + /* + The candidates of homogeneityOptimisation.py. + + That script lists nineteen radii and twenty magnet counts, and pairs + them by index, so its last count is never used. The pairing is written + out here rather than left to two lists of different lengths that happen + to be zipped. + */ + const float radii[] = { + 0.148f, 0.151f, 0.154f, 0.156f, 0.159f, 0.162f, 0.165f, + 0.168f, 0.171f, 0.174f, 0.177f, 0.180f, 0.183f, 0.186f, + 0.189f, 0.192f, 0.195f, 0.198f, 0.201f + }; + + const size_t count = sizeof(radii) / sizeof(radii[0]); + + candidates.reserve(count); + + for (size_t i = 0; i < count; i++) { + RingCandidate candidate; + candidate.radius = radii[i]; + candidate.count = (uint32_t) (50 + i); + candidates.push_back(candidate); + } +} + + +void HalbachSpec::readJSON(const nlohmann::json& data) { + + if (!data.is_object()) { + throw std::invalid_argument( + "A Halbach optimization configuration is an object"); + } + + ring_count = readCount(data, "ring_count", ring_count, 1); + ring_separation = readPositive(data, "ring_separation", ring_separation); + + outer_radius_offset = + readNumber(data, "outer_radius_offset", outer_radius_offset); + + outer_count_offset = + readCount(data, "outer_count_offset", outer_count_offset, 0); + + if (data.contains("order")) { + if (!data["order"].is_number_integer()) { + throw std::invalid_argument( + "The \"order\" of a Halbach ring is a whole number"); + } + order = data["order"].get(); + } + + if (data.contains("candidates")) { + + const nlohmann::json& given = data["candidates"]; + + if (!given.is_array() || given.empty()) { + throw std::invalid_argument( + "\"candidates\" is a non empty array of {\"radius\", \"count\"}"); + } + + candidates.clear(); + candidates.reserve(given.size()); + + for (const auto& entry : given) { + + if (!entry.is_object()) { + throw std::invalid_argument( + "A ring candidate is an object with \"radius\" and \"count\""); + } + + RingCandidate candidate; + candidate.radius = readPositive(entry, "radius", 0.0f); + candidate.count = readCount(entry, "count", 0, 1); + + if (candidate.radius <= 0.0f) { + throw std::invalid_argument( + "A ring candidate needs a \"radius\""); + } + + candidates.push_back(candidate); + } + } + + if (data.contains("element")) { + + const nlohmann::json& given = data["element"]; + + element.size = readPositive(given, "size", element.size); + element.remanence = readNumber(given, "remanence", element.remanence); + } + + dsv = readPositive(data, "dsv", dsv); + resolution = readPositive(data, "resolution", resolution); + + const std::string symmetry_name = readChoice( + data, "symmetry", {"octant", "hemisphere", "full"}, toString(symmetry)); + + symmetry = symmetry_name == "octant" ? Symmetry::Octant + : symmetry_name == "hemisphere" ? Symmetry::Hemisphere + : Symmetry::Full; + + const std::string objective_name = readChoice( + data, "objective", {"ppm_bx", "ppm_bmag"}, toString(objective)); + + objective = objective_name == "ppm_bx" ? Objective::Bx : Objective::BMagnitude; + + const std::string model_name = readChoice( + data, "field_model", {"cuboid", "dipole"}, toString(field_model)); + + field_model = model_name == "cuboid" ? FieldModel::Cuboid : FieldModel::Dipole; + + if (data.contains("genetic")) { + + const nlohmann::json& given = data["genetic"]; + + genetic.population = readCount(given, "population", genetic.population, 2); + genetic.generations = readCount(given, "generations", genetic.generations, 1); + genetic.tournament = readCount(given, "tournament", genetic.tournament, 1); + genetic.elitism = readCount(given, "elitism", genetic.elitism, 0); + + genetic.cx_prob = readProbability(given, "cx_prob", genetic.cx_prob); + genetic.mut_prob = readProbability(given, "mut_prob", genetic.mut_prob); + genetic.gene_mut_prob = + readProbability(given, "gene_mut_prob", genetic.gene_mut_prob); + + if (given.contains("seed")) { + if (!given["seed"].is_number_integer() && + !given["seed"].is_number_unsigned()) { + throw std::invalid_argument("The \"seed\" is a whole number"); + } + genetic.seed = given["seed"].get(); + } + } + + validate(); +} + + +void HalbachSpec::validate() const { + + if (candidates.empty()) { + throw std::invalid_argument( + "A Halbach optimization needs at least one ring candidate"); + } + + // A gene is stored as a 16 bit index, which is what the population view + // holds, and there is no sensible run with more candidates than that. + if (candidates.size() > 65535) { + throw std::invalid_argument( + "A Halbach optimization takes at most 65535 ring candidates"); + } + + for (const auto& candidate : candidates) { + + if (!(candidate.radius > 0.0f)) { + throw std::invalid_argument("A ring candidate needs a positive radius"); + } + + if (candidate.count < 1) { + throw std::invalid_argument("A ring candidate needs at least one magnet"); + } + + if (candidate.radius + outer_radius_offset <= 0.0f) { + throw std::invalid_argument( + "The \"outer_radius_offset\" leaves the outer ring at a radius " + "of zero or less"); + } + } + + // The sampled sphere has to fit inside the bore, or the observation points + // land inside the magnets and the field there is neither uniform nor, + // for a dipole, finite. + float smallest = candidates[0].radius; + for (const auto& candidate : candidates) { + smallest = std::min(smallest, candidate.radius); + smallest = std::min(smallest, candidate.radius + outer_radius_offset); + } + + if (0.5f * dsv >= smallest - 0.5f * element.size) { + throw std::invalid_argument( + "The sampled sphere reaches the magnets: half the \"dsv\" is at " + "least the smallest candidate radius less half the element size"); + } + + if (resolution > dsv) { + throw std::invalid_argument( + "The \"resolution\" is coarser than the whole \"dsv\""); + } + + if (genetic.elitism >= genetic.population) { + throw std::invalid_argument( + "The \"elitism\" carries over the whole population, which leaves " + "no room for a next generation"); + } + + if (element.remanence == 0.0f) { + throw std::invalid_argument( + "An element of zero remanence produces no field at all"); + } +} + + +std::vector HalbachSpec::getRingPositions() const { + + std::vector positions; + positions.reserve(ring_count); + + const float length = (float) (ring_count - 1) * ring_separation; + + for (uint32_t i = 0; i < ring_count; i++) { + positions.push_back(-0.5f * length + (float) i * ring_separation); + } + + return positions; +} + + +std::vector HalbachSpec::getSymmetricRingPositions() const { + + std::vector positions; + + /* + A ring the arithmetic puts at the centre may land a rounding error to + either side of it, and the ring at the centre is the one that is not + mirrored. A residue read as a position of its own would put two rings a + nanometre apart where there should be one, so the centre is snapped + shut. The tolerance is relative to the spacing, since that is the only + length here that says what "close to zero" means. + */ + const float tolerance = 1e-4f * ring_separation; + + for (const auto& position : getRingPositions()) { + + if (std::fabs(position) < tolerance) { + positions.push_back(0.0f); + } else if (position > 0.0f) { + positions.push_back(position); + } + } + + return positions; +} + + +size_t HalbachSpec::getNumGenes() const { + return getSymmetricRingPositions().size(); +} + + +size_t HalbachSpec::getNumCandidates() const { + return candidates.size(); +} + + +float HalbachSpec::getDipoleMoment() const { + + // m = V * J / mu0, with J the remanence of the material and V the volume + // of the cube. This is the same number as the dip_mom of halbachFields.py, + // which writes it as Br * a^3 / (4 pi * 1e-7). + const double volume = (double) element.size * (double) element.size + * (double) element.size; + + return (float) (volume * (double) element.remanence / MU_0); +} + + +nlohmann::json HalbachSpec::getElementJSON() const { + + // The polarization, or the moment, points along the local x of the + // element. A Halbach ring turns its members about its own axis, carrying + // that direction round with them, so a member at angle theta ends up + // polarized along (cos((order + 1) theta), sin((order + 1) theta), 0). + // With order one that is the 2 theta of the Python script. + if (field_model == FieldModel::Dipole) { + + return nlohmann::json{ + {"type", "dipole"}, + {"parameters", { + {"moment", {tidy(getDipoleMoment()), 0.0, 0.0}} + }} + }; + } + + return nlohmann::json{ + {"type", "cuboid"}, + {"parameters", { + {"dimensions", {tidy(element.size), tidy(element.size), + tidy(element.size)}}, + {"magnetization", {tidy(element.remanence), 0.0, 0.0}} + }} + }; +} + + +nlohmann::json HalbachSpec::makeRingJSON( + const float& radius, const uint32_t& count, const float& z) const { + + return nlohmann::json{ + {"type", "halbach_ring"}, + {"parameters", { + {"radius", tidy(radius)}, + {"count", count}, + {"order", order}, + {"position", {0.0, 0.0, tidy(z)}}, + {"orientation", {1.0, 0.0, 0.0, 0.0}}, + {"element", getElementJSON()} + }} + }; +} + + +std::vector HalbachSpec::makeRingsForGene( + const size_t& gene, const size_t& candidate) const { + + const std::vector positions = getSymmetricRingPositions(); + + if (gene >= positions.size()) { + throw std::out_of_range("Gene index out of range"); + } + + if (candidate >= candidates.size()) { + throw std::out_of_range("Ring candidate index out of range"); + } + + const RingCandidate& chosen = candidates[candidate]; + + const float z = positions[gene]; + + std::vector rings; + rings.reserve(4); + + // The two layers move together, which is what makes one gene enough for a + // position: an inner ring, and an outer ring further out with more magnets + // in it. + rings.push_back(makeRingJSON(chosen.radius, chosen.count, z)); + rings.push_back(makeRingJSON(chosen.radius + outer_radius_offset, + chosen.count + outer_count_offset, z)); + + // The magnet is built mirror symmetric about the middle, so every gene but + // the one at the centre stands for a pair of positions. + if (z != 0.0f) { + rings.push_back(makeRingJSON(chosen.radius, chosen.count, -z)); + rings.push_back(makeRingJSON(chosen.radius + outer_radius_offset, + chosen.count + outer_count_offset, -z)); + } + + return rings; +} + + +nlohmann::json HalbachSpec::toJSON() const { + + nlohmann::json candidate_list = nlohmann::json::array(); + + for (const auto& candidate : candidates) { + candidate_list.push_back({ + {"radius", tidy(candidate.radius)}, + {"count", candidate.count} + }); + } + + return nlohmann::json{ + {"ring_count", ring_count}, + {"ring_separation", tidy(ring_separation)}, + {"candidates", candidate_list}, + {"outer_radius_offset", tidy(outer_radius_offset)}, + {"outer_count_offset", outer_count_offset}, + {"order", order}, + {"element", { + {"size", tidy(element.size)}, + {"remanence", tidy(element.remanence)} + }}, + {"dsv", tidy(dsv)}, + {"resolution", tidy(resolution)}, + {"symmetry", toString(symmetry)}, + {"objective", toString(objective)}, + {"field_model", toString(field_model)}, + {"genetic", { + {"population", genetic.population}, + {"generations", genetic.generations}, + {"cx_prob", tidy(genetic.cx_prob)}, + {"mut_prob", tidy(genetic.mut_prob)}, + {"gene_mut_prob", tidy(genetic.gene_mut_prob)}, + {"tournament", genetic.tournament}, + {"elitism", genetic.elitism}, + {"seed", genetic.seed} + }} + }; +} + +} // namespace optimization +} // namespace greeter diff --git a/source/optimization/HomogeneityObjective.cpp b/source/optimization/HomogeneityObjective.cpp new file mode 100644 index 0000000..30fa840 --- /dev/null +++ b/source/optimization/HomogeneityObjective.cpp @@ -0,0 +1,243 @@ +#include + +#include +#include +#include + + +namespace greeter { +namespace optimization { + +namespace { + +typedef Kokkos::TeamPolicy TeamPolicy; +typedef TeamPolicy::member_type TeamMember; + +typedef Kokkos::View ScratchIndices; + +/* + A mean of zero would divide the homogeneity by nothing. It means the field + averages out over the volume, which is a genome worth rejecting rather than + one worth an infinity, so it is scored as the worst there is. +*/ +KOKKOS_INLINE_FUNCTION +float partsPerMillion(const FieldExtent& extent, const double& count) { + + const double mean = extent.sum / count; + const double magnitude = mean < 0.0 ? -mean : mean; + + if (!(magnitude > 0.0)) { + return Kokkos::reduction_identity::min(); + } + + return (float) (1e6 * (double) (extent.max - extent.min) / magnitude); +} + +} // namespace + + +HomogeneityObjective::HomogeneityObjective() {} + +HomogeneityObjective::HomogeneityObjective(const HalbachBasis& basis, + const Objective& _objective): + fields(basis.getFields()), + num_genes(basis.getNumGenes()), + num_candidates(basis.getNumCandidates()), + num_points(basis.getNumPoints()), + objective(_objective) {} + + +void HomogeneityObjective::evaluate(const PopulationView& population, + const FitnessView& fitness) const { + + const size_t num_individuals = population.extent(0); + + if (population.extent(1) != num_genes) { + throw std::invalid_argument( + "The population has " + std::to_string(population.extent(1)) + + " genes an individual, the basis was built for " + + std::to_string(num_genes)); + } + + if (fitness.extent(0) != num_individuals) { + throw std::invalid_argument( + "There is one fitness per individual"); + } + + // Copied into locals: a lambda that captured `this` would dereference a + // host pointer inside the parallel region. + const Kokkos::View the_fields = fields; + const size_t the_num_genes = num_genes; + const size_t the_num_candidates = num_candidates; + const size_t the_num_points = num_points; + + const double count = (double) num_points; + + TeamPolicy policy((int) num_individuals, Kokkos::AUTO); + + policy.set_scratch_size( + 0, Kokkos::PerTeam((int) ScratchIndices::shmem_size(the_num_genes))); + + const bool magnitude = objective == Objective::BMagnitude; + + Kokkos::parallel_for( + "homogeneity", policy, + KOKKOS_LAMBDA(const TeamMember& team) { + + const int individual = team.league_rank(); + + /* + Where in the basis each gene of this individual points. Worked + out once for the team rather than once per observation point, + which is the difference between one multiply-add per gene and + one per gene per point. + */ + ScratchIndices configuration(team.team_scratch(0), the_num_genes); + + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team, the_num_genes), + [&](const size_t& gene) { + configuration(gene) = (uint32_t) + (gene * the_num_candidates + population(individual, gene)); + }); + + team.team_barrier(); + + FieldExtent extent; + + /* + One pass over the volume. The peak, the trough and the total are + all the homogeneity needs, and reading the basis three times to + get them separately would cost three times the memory traffic of + the only expensive step here. + */ + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange(team, the_num_points), + [&](const size_t& point, FieldExtent& local) { + + float value = 0.0f; + + if (magnitude) { + + float b_x = 0.0f; + float b_y = 0.0f; + float b_z = 0.0f; + + for (size_t gene = 0; gene < the_num_genes; gene++) { + const uint32_t c = configuration(gene); + b_x += the_fields(c, 0, point); + b_y += the_fields(c, 1, point); + b_z += the_fields(c, 2, point); + } + + value = sqrt(b_x * b_x + b_y * b_y + b_z * b_z); + + } else { + + for (size_t gene = 0; gene < the_num_genes; gene++) { + value += the_fields(configuration(gene), 0, point); + } + } + + if (value < local.min) { local.min = value; } + if (value > local.max) { local.max = value; } + local.sum += (double) value; + }, + FieldExtentReducer(extent)); + + Kokkos::single(Kokkos::PerTeam(team), [&]() { + fitness(individual) = partsPerMillion(extent, count); + }); + }); + + Kokkos::fence(); +} + + +FieldMetrics HomogeneityObjective::evaluateGenome( + const std::vector& genome) const { + + if (genome.size() != num_genes) { + throw std::invalid_argument( + "The genome has " + std::to_string(genome.size()) + + " genes, the objective was built for " + std::to_string(num_genes)); + } + + std::vector values; + values.reserve(num_points); + + for (size_t point = 0; point < num_points; point++) { + + float b[3] = {0.0f, 0.0f, 0.0f}; + + for (size_t gene = 0; gene < num_genes; gene++) { + + if (genome[gene] >= num_candidates) { + throw std::out_of_range("Gene value out of range"); + } + + const size_t configuration = gene * num_candidates + genome[gene]; + + b[0] += fields(configuration, 0, point); + b[1] += fields(configuration, 1, point); + b[2] += fields(configuration, 2, point); + } + + values.push_back( + objective == Objective::BMagnitude + ? std::sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]) + : b[0]); + } + + return summarise(values); +} + + +FieldMetrics HomogeneityObjective::summarise(const std::vector& values) { + + FieldMetrics metrics; + + if (values.empty()) { + return metrics; + } + + metrics.min = values[0]; + metrics.max = values[0]; + + // In double, because a sum of a hundred thousand single precision samples + // of nearly the same size loses digits the homogeneity is measured in. + double total = 0.0; + + for (const auto& value : values) { + metrics.min = std::min(metrics.min, value); + metrics.max = std::max(metrics.max, value); + total += (double) value; + } + + metrics.mean = (float) (total / (double) values.size()); + + const float magnitude = std::fabs(metrics.mean); + + metrics.ppm = magnitude > 0.0f + ? 1e6f * (metrics.max - metrics.min) / magnitude + : std::numeric_limits::max(); + + return metrics; +} + + +size_t HomogeneityObjective::getNumGenes() const { + return num_genes; +} + +size_t HomogeneityObjective::getNumCandidates() const { + return num_candidates; +} + +size_t HomogeneityObjective::getNumPoints() const { + return num_points; +} + +} // namespace optimization +} // namespace greeter diff --git a/test/source/halbachOptimization.cpp b/test/source/halbachOptimization.cpp new file mode 100644 index 0000000..d54af78 --- /dev/null +++ b/test/source/halbachOptimization.cpp @@ -0,0 +1,580 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +/* + The optimizer of a Halbach cylinder. + + Everything here rests on one claim: that the field of a choice of rings is + the sum of the fields of those rings, worked out once each. So the claim + worth checking hardest is that claim itself, against the simulator the rest + of this library uses. The genetic algorithm on top of it is checked for the + properties an algorithm of that kind has to have rather than for a number, + because the number it lands on is not repeatable across machines and is not + supposed to be. +*/ + +namespace { + +using namespace greeter::optimization; + +/* + A specification small enough to run inside a test: three ring positions, + two candidates, and a volume of a few dozen points. The physics is the same + as the full one, only less of it. +*/ +HalbachSpec smallSpec() { + + HalbachSpec spec; + + spec.ring_count = 3; + spec.ring_separation = 0.03f; + + spec.candidates.clear(); + spec.candidates.push_back({0.10f, 8}); + spec.candidates.push_back({0.12f, 10}); + + spec.outer_radius_offset = 0.02f; + spec.outer_count_offset = 2; + + spec.element.size = 0.01f; + spec.element.remanence = 1.3f; + + spec.dsv = 0.04f; + spec.resolution = 0.01f; + + spec.genetic.population = 64; + spec.genetic.generations = 6; + spec.genetic.seed = 7; + + return spec; +} + +} // namespace + + +TEST_CASE("the specification lays the cylinder out the way the Python script does") { + + const HalbachSpec spec; // the defaults, which are that script's + + SUBCASE("the rings are centred and evenly spaced") { + + const std::vector positions = spec.getRingPositions(); + + REQUIRE(positions.size() == 23); + + CHECK(positions.front() == doctest::Approx(-0.242f)); + CHECK(positions.back() == doctest::Approx(0.242f)); + + for (size_t i = 1; i < positions.size(); i++) { + CHECK(positions[i] - positions[i - 1] == doctest::Approx(0.022f)); + } + } + + SUBCASE("a genome covers the non negative half, mirrored") { + + const std::vector half = spec.getSymmetricRingPositions(); + + CHECK(half.size() == 12); + CHECK(spec.getNumGenes() == 12); + CHECK(half.front() == doctest::Approx(0.0f)); + CHECK(half.back() == doctest::Approx(0.242f)); + } + + SUBCASE("the gene at the centre stands for one pair of rings, the others two") { + + CHECK(spec.makeRingsForGene(0, 0).size() == 2); + CHECK(spec.makeRingsForGene(1, 0).size() == 4); + } + + SUBCASE("the two layers of a gene differ by the offsets") { + + const std::vector rings = spec.makeRingsForGene(3, 5); + + const float inner = rings[0]["parameters"]["radius"].get(); + const float outer = rings[1]["parameters"]["radius"].get(); + + CHECK(outer - inner == doctest::Approx(spec.outer_radius_offset)); + + CHECK(rings[1]["parameters"]["count"].get() + - rings[0]["parameters"]["count"].get() + == spec.outer_count_offset); + } + + SUBCASE("the mirrored rings sit opposite each other") { + + const std::vector rings = spec.makeRingsForGene(4, 0); + + REQUIRE(rings.size() == 4); + + CHECK(rings[0]["parameters"]["position"][2].get() + == doctest::Approx(-rings[2]["parameters"]["position"][2].get())); + } + + SUBCASE("the dipole model carries the moment of the cube it replaces") { + + HalbachSpec dipole_spec; + dipole_spec.field_model = FieldModel::Dipole; + + // m = V * Br / mu0, which is the dip_mom of halbachFields.py. + const double mu_0 = 4.0e-7 * M_PI; + const double expected = std::pow((double) dipole_spec.element.size, 3) + * (double) dipole_spec.element.remanence / mu_0; + + CHECK(dipole_spec.getDipoleMoment() == doctest::Approx((float) expected)); + + CHECK(dipole_spec.getElementJSON()["type"].get() == "dipole"); + } + + SUBCASE("a specification that reaches into the magnets is refused") { + + HalbachSpec too_wide; + too_wide.dsv = 0.6f; // wider than the smallest candidate radius + + CHECK_THROWS_AS(too_wide.validate(), std::invalid_argument); + } +} + + +TEST_CASE("the centre ring is the one that is not mirrored, at any ring count") { + + /* + The position of the middle ring is arrived at by arithmetic that need not + land on exactly zero, and a residue read as a position of its own would + put two rings a nanometre apart where there should be one. So this is + checked over a range of counts and spacings rather than the one pair the + defaults happen to use. + */ + for (const uint32_t count : {1u, 2u, 3u, 4u, 7u, 22u, 23u, 24u, 41u}) { + for (const float separation : {0.022f, 0.007f, 0.13f, 0.001f}) { + + HalbachSpec spec = smallSpec(); + spec.ring_count = count; + spec.ring_separation = separation; + + const std::vector half = spec.getSymmetricRingPositions(); + + // An odd count has a ring at the centre, an even one does not, and + // either way the halves account for every ring exactly once. + const bool odd = count % 2 == 1; + + CHECK(half.size() == (odd ? (count + 1) / 2 : count / 2)); + + size_t rings = 0; + + for (size_t gene = 0; gene < half.size(); gene++) { + + const size_t made = spec.makeRingsForGene(gene, 0).size(); + + // Two layers at one position, or two layers at each of two. + CHECK((made == 2 || made == 4)); + CHECK((made == 2) == (half[gene] == 0.0f)); + + rings += made; + } + + CHECK(rings == 2 * count); + + if (odd) { + CHECK(half.front() == 0.0f); + } else { + CHECK(half.front() == doctest::Approx(0.5f * separation)); + } + } + } +} + + +TEST_CASE("the sampled volume is the part of the sphere the symmetry asks for") { + + HalbachSpec spec = smallSpec(); + + spec.symmetry = Symmetry::Full; + const std::vector> full = makeSamplePoints(spec); + + spec.symmetry = Symmetry::Hemisphere; + const std::vector> half = makeSamplePoints(spec); + + spec.symmetry = Symmetry::Octant; + const std::vector> octant = makeSamplePoints(spec); + + CHECK(octant.size() < half.size()); + CHECK(half.size() < full.size()); + + const float radius = 0.5f * spec.dsv; + + for (const auto& point : full) { + CHECK(point[0] * point[0] + point[1] * point[1] + point[2] * point[2] + <= doctest::Approx(radius * radius).epsilon(1e-3)); + } + + for (const auto& point : octant) { + CHECK(point[0] >= -1e-6f); + CHECK(point[1] >= -1e-6f); + CHECK(point[2] >= -1e-6f); + } +} + + +TEST_CASE("a basis field is the field of the rings it stands for") { + + const HalbachSpec spec = smallSpec(); + + const HalbachBasis basis = HalbachBasis::build(spec, false); + + REQUIRE(basis.getNumGenes() == spec.getNumGenes()); + REQUIRE(basis.getNumCandidates() == spec.getNumCandidates()); + REQUIRE(basis.getNumPoints() > 0); + + /* + The claim the whole optimizer rests on, checked against the simulator + every other part of this library uses: what the packed kernel wrote into + the basis for one configuration is what a collection of those same rings + simulates to. + */ + for (size_t gene = 0; gene < basis.getNumGenes(); gene++) { + for (size_t candidate = 0; candidate < basis.getNumCandidates(); candidate++) { + + greeter::MagnetCollection collection; + + for (const auto& ring : spec.makeRingsForGene(gene, candidate)) { + for (auto& member : greeter::HalbachRingArrangement::expand(ring)) { + collection.addMagnet(std::move(member)); + } + } + + const std::vector> expected = + collection.simulate(basis.getPoints()); + + const size_t configuration = gene * basis.getNumCandidates() + candidate; + + for (size_t point = 0; point < basis.getNumPoints(); point++) { + for (size_t component = 0; component < 3; component++) { + CHECK(basis.getFields()(configuration, component, point) + == doctest::Approx(expected[point][component]).epsilon(1e-5)); + } + } + } + } +} + + +TEST_CASE("summing the basis is simulating the magnet it describes") { + + const HalbachSpec spec = smallSpec(); + + const HalbachBasis basis = HalbachBasis::build(spec, false); + + // One candidate per gene, chosen so that not every gene picks the same one. + std::vector genome; + for (size_t gene = 0; gene < basis.getNumGenes(); gene++) { + genome.push_back((uint16_t) (gene % basis.getNumCandidates())); + } + + HalbachSolution solution; + solution.spec = spec; + solution.genome = genome; + + const greeter::MagnetCollection collection = solution.buildCollection(); + + const std::vector> expected = + collection.simulate(basis.getPoints()); + + for (size_t point = 0; point < basis.getNumPoints(); point++) { + + float b[3]; + basis.evaluateAt(genome, point, b); + + for (size_t component = 0; component < 3; component++) { + CHECK(b[component] == doctest::Approx(expected[point][component]).epsilon(1e-4)); + } + } +} + + +TEST_CASE("the objective on the device agrees with the objective on the host") { + + HalbachSpec spec = smallSpec(); + + for (const auto& kind : {Objective::Bx, Objective::BMagnitude}) { + + spec.objective = kind; + + const HalbachBasis basis = HalbachBasis::build(spec, false); + const HomogeneityObjective objective(basis, kind); + + const size_t genes = basis.getNumGenes(); + const size_t candidates = basis.getNumCandidates(); + + const size_t individuals = 8; + + PopulationView population("population", individuals, genes); + FitnessView fitness("fitness", individuals); + + std::vector> genomes; + + for (size_t i = 0; i < individuals; i++) { + + std::vector genome; + + for (size_t gene = 0; gene < genes; gene++) { + const uint16_t value = (uint16_t) ((i + gene) % candidates); + genome.push_back(value); + population(i, gene) = value; + } + + genomes.push_back(genome); + } + + objective.evaluate(population, fitness); + + for (size_t i = 0; i < individuals; i++) { + CHECK(fitness(i) + == doctest::Approx(objective.evaluateGenome(genomes[i]).ppm) + .epsilon(1e-3)); + } + } +} + + +TEST_CASE("the evolution improves, keeps its best, and stays inside the alphabet") { + + const HalbachSpec spec = smallSpec(); + + const HalbachBasis basis = HalbachBasis::build(spec, false); + const HomogeneityObjective objective(basis, spec.objective); + + GeneticOptimizer optimizer( + spec.genetic, basis.getNumGenes(), basis.getNumCandidates()); + + const std::vector best = optimizer.run(objective, nullptr, false); + + REQUIRE(best.size() == basis.getNumGenes()); + + SUBCASE("every gene names a candidate that exists") { + for (const auto& gene : best) { + CHECK(gene < basis.getNumCandidates()); + } + } + + SUBCASE("the best ever never gets worse") { + + const std::vector& history = optimizer.getHistory(); + + REQUIRE(history.size() == spec.genetic.generations + 1); + + for (size_t i = 1; i < history.size(); i++) { + CHECK(history[i].best_ever <= history[i - 1].best_ever); + } + } + + SUBCASE("the reported best is the score of the reported genome") { + CHECK(optimizer.getBestFitness() + == doctest::Approx(objective.evaluateGenome(best).ppm).epsilon(1e-3)); + } + + SUBCASE("the answer is at least as good as the first generation was") { + CHECK(optimizer.getBestFitness() <= optimizer.getHistory().front().best); + } +} + + +TEST_CASE("the same seed gives the same answer") { + + const HalbachSpec spec = smallSpec(); + + const HalbachBasis basis = HalbachBasis::build(spec, false); + const HomogeneityObjective objective(basis, spec.objective); + + GeneticOptimizer first( + spec.genetic, basis.getNumGenes(), basis.getNumCandidates()); + GeneticOptimizer again( + spec.genetic, basis.getNumGenes(), basis.getNumCandidates()); + + CHECK(first.run(objective, nullptr, false) + == again.run(objective, nullptr, false)); + + GeneticSettings other = spec.genetic; + other.seed = spec.genetic.seed + 1; + + GeneticOptimizer different( + other, basis.getNumGenes(), basis.getNumCandidates()); + + different.run(objective, nullptr, false); + + // Not an equality check the other way round: two seeds may land on the same + // answer, and on a problem this small they often do. What has to hold is + // that the second run is a run of the same problem. + CHECK(different.getBestFitness() > 0.0f); +} + + +TEST_CASE("a run asked to stop stops") { + + const HalbachSpec spec = smallSpec(); + + const HalbachBasis basis = HalbachBasis::build(spec, false); + const HomogeneityObjective objective(basis, spec.objective); + + struct StopAtTwo: public GenerationSink { + bool onGeneration(const GenerationRecord& record) override { + return record.generation < 2; + } + }; + + StopAtTwo sink; + + GeneticOptimizer optimizer( + spec.genetic, basis.getNumGenes(), basis.getNumCandidates()); + + optimizer.run(objective, &sink, false); + + CHECK(optimizer.wasStopped()); + CHECK(optimizer.getHistory().size() == 3); +} + + +TEST_CASE("what the optimizer writes is what the simulator reads") { + + HalbachSpec spec = smallSpec(); + spec.genetic.generations = 3; + + greeter::optimization::RunOptions options; + options.verbose = false; + options.verify = true; + + const HalbachSolution solution = + HalbachOptimizer::run(spec, nullptr, options); + + REQUIRE(solution.was_verified); + REQUIRE(solution.num_magnets > 0); + + for (const auto& emit : {greeter::HalbachEmit::Arrangements, + greeter::HalbachEmit::Magnets}) { + + const nlohmann::json written = + greeter::HalbachOptimizationIO::write(solution, emit); + + // The file this library reads, not a report beside one. + CHECK(greeter::MagnetIO::validateJSON(written)); + + const greeter::Scene scene = greeter::SceneIO::read(written); + + CHECK(scene.collection.get_num_magnets() == solution.num_magnets); + + /* + And the same magnet, not merely the same number of them: read back and + simulated, it gives the homogeneity the run reported. + */ + size_t points = 0; + + const FieldMetrics measured = HalbachOptimizer::measure( + scene.collection, solution.spec, solution.spec.objective, points); + + CHECK(points == solution.verified_points); + CHECK(measured.mean == doctest::Approx(solution.verified.mean).epsilon(1e-4)); + CHECK(measured.ppm == doctest::Approx(solution.verified.ppm).epsilon(1e-3)); + } +} + + +TEST_CASE("the written record reads back as the specification that made it") { + + HalbachSpec spec = smallSpec(); + spec.field_model = FieldModel::Dipole; + spec.objective = Objective::BMagnitude; + spec.symmetry = Symmetry::Hemisphere; + spec.genetic.generations = 2; + + greeter::optimization::RunOptions options; + options.verbose = false; + options.verify = false; + + const HalbachSolution solution = + HalbachOptimizer::run(spec, nullptr, options); + + const nlohmann::json written = greeter::HalbachOptimizationIO::write( + solution, greeter::HalbachEmit::Arrangements); + + const HalbachSpec read = greeter::HalbachOptimizationIO::readSpec(written); + + CHECK(read.ring_count == spec.ring_count); + CHECK(read.getNumCandidates() == spec.getNumCandidates()); + // Through their names: doctest cannot print a scoped enumeration, and the + // name is what a failure would want to show anyway. + CHECK(std::string(toString(read.field_model)) == toString(spec.field_model)); + CHECK(std::string(toString(read.objective)) == toString(spec.objective)); + CHECK(std::string(toString(read.symmetry)) == toString(spec.symmetry)); + CHECK(read.genetic.seed == spec.genetic.seed); + CHECK(read.dsv == doctest::Approx(spec.dsv)); + CHECK(read.element.size == doctest::Approx(spec.element.size)); + + SUBCASE("a dipole element is written as a dipole") { + CHECK(written["arrangements"][0]["parameters"]["element"]["type"] + .get() == "dipole"); + } +} + + +TEST_CASE("a configuration is read, and a wrong one is refused by name") { + + SUBCASE("what a file overrides is overridden, and the rest is left alone") { + + const nlohmann::json data = nlohmann::json::parse(R"({ + "halbach_optimization": { + "ring_count": 5, + "dsv": 0.1, + "field_model": "dipole", + "genetic": { "population": 32, "generations": 4 } + } + })"); + + const HalbachSpec spec = greeter::HalbachOptimizationIO::readSpec(data); + + CHECK(spec.ring_count == 5); + CHECK(spec.dsv == doctest::Approx(0.1f)); + CHECK(std::string(toString(spec.field_model)) == "dipole"); + CHECK(spec.genetic.population == 32); + CHECK(spec.genetic.generations == 4); + + // Untouched by the file, so still the default of the Python script. + CHECK(spec.ring_separation == doctest::Approx(0.022f)); + CHECK(spec.getNumCandidates() == 19); + CHECK(spec.genetic.tournament == 3); + } + + SUBCASE("an unknown choice is refused rather than quietly ignored") { + + HalbachSpec spec; + + CHECK_THROWS_AS( + spec.readJSON(nlohmann::json::parse(R"({"field_model": "monopole"})")), + std::invalid_argument); + + CHECK_THROWS_AS( + spec.readJSON(nlohmann::json::parse(R"({"symmetry": "quadrant"})")), + std::invalid_argument); + + CHECK_THROWS_AS( + spec.readJSON(nlohmann::json::parse(R"({"dsv": -1})")), + std::invalid_argument); + + CHECK_THROWS_AS( + spec.readJSON(nlohmann::json::parse(R"({"genetic": {"cx_prob": 2}})")), + std::invalid_argument); + } +}