1 package net.bmahe.genetics4j.samples.mixturemodel;
2
3 import java.io.IOException;
4 import java.util.Collection;
5 import java.util.Comparator;
6 import java.util.HashSet;
7 import java.util.List;
8 import java.util.Set;
9 import java.util.concurrent.ExecutionException;
10 import java.util.concurrent.Executors;
11 import java.util.stream.IntStream;
12
13 import org.apache.commons.collections4.CollectionUtils;
14 import org.apache.commons.math3.distribution.MultivariateNormalDistribution;
15 import org.apache.commons.math3.exception.MathUnsupportedOperationException;
16 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
17 import org.apache.commons.math3.linear.EigenDecomposition;
18 import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
19 import org.apache.commons.math3.linear.SingularMatrixException;
20 import org.apache.logging.log4j.LogManager;
21 import org.apache.logging.log4j.Logger;
22
23 import net.bmahe.genetics4j.core.Genotype;
24 import net.bmahe.genetics4j.core.chromosomes.FloatChromosome;
25 import net.bmahe.genetics4j.core.evolutionlisteners.EvolutionListeners;
26 import net.bmahe.genetics4j.core.spec.EvolutionResult;
27 import net.bmahe.genetics4j.core.spec.chromosome.FloatChromosomeSpec;
28 import net.bmahe.genetics4j.core.spec.combination.MultiCombinations;
29 import net.bmahe.genetics4j.core.spec.combination.MultiPointArithmetic;
30 import net.bmahe.genetics4j.core.spec.combination.MultiPointCrossover;
31 import net.bmahe.genetics4j.core.spec.combination.SinglePointArithmetic;
32 import net.bmahe.genetics4j.core.spec.combination.SinglePointCrossover;
33 import net.bmahe.genetics4j.core.spec.mutation.CreepMutation;
34 import net.bmahe.genetics4j.core.spec.mutation.MultiMutations;
35 import net.bmahe.genetics4j.core.spec.mutation.RandomMutation;
36 import net.bmahe.genetics4j.core.spec.mutation.SwapMutation;
37 import net.bmahe.genetics4j.core.spec.replacement.Elitism;
38 import net.bmahe.genetics4j.core.spec.statistics.distributions.NormalDistribution;
39 import net.bmahe.genetics4j.core.termination.Terminations;
40 import net.bmahe.genetics4j.extras.evolutionlisteners.CSVEvolutionListener;
41 import net.bmahe.genetics4j.extras.evolutionlisteners.ColumnExtractor;
42 import net.bmahe.genetics4j.gpu.GPUEASystemFactory;
43 import net.bmahe.genetics4j.gpu.spec.DeviceFilters;
44 import net.bmahe.genetics4j.gpu.spec.GPUEAConfiguration;
45 import net.bmahe.genetics4j.gpu.spec.GPUEAExecutionContext;
46 import net.bmahe.genetics4j.gpu.spec.Program;
47 import net.bmahe.genetics4j.gpu.spec.fitness.FitnessExtractor;
48 import net.bmahe.genetics4j.gpu.spec.fitness.SingleKernelFitness;
49 import net.bmahe.genetics4j.gpu.spec.fitness.SingleKernelFitnessDescriptor;
50 import net.bmahe.genetics4j.gpu.spec.fitness.cldata.DataLoaders;
51 import net.bmahe.genetics4j.gpu.spec.fitness.cldata.DataSupplier;
52 import net.bmahe.genetics4j.gpu.spec.fitness.cldata.ResultAllocators;
53 import net.bmahe.genetics4j.gpu.spec.fitness.cldata.StaticDataLoaders;
54 import net.bmahe.genetics4j.gpu.spec.fitness.kernelcontext.KernelExecutionContextComputers;
55 import net.bmahe.genetics4j.moo.FitnessVector;
56 import net.bmahe.genetics4j.moo.nsga2.impl.NSGA2SelectionPolicyHandler;
57 import net.bmahe.genetics4j.moo.nsga2.impl.TournamentNSGA2SelectionPolicyHandler;
58 import net.bmahe.genetics4j.moo.nsga2.spec.NSGA2Selection;
59 import net.bmahe.genetics4j.moo.nsga2.spec.TournamentNSGA2Selection;
60 import net.bmahe.genetics4j.moo.spea2.replacement.SPEA2ReplacementStrategyHandler;
61
62 public class MooGPU {
63 public static final Logger logger = LogManager.getLogger(MooGPU.class);
64
65 private final int distributionNumParameters;
66 private final Comparator<Genotype> deduplicator;
67 private final String baseDir;
68 private final int maxGenerations;
69
70 public MooGPU(final int _distributionNumParameters,
71 final Comparator<Genotype> _deduplicator,
72 final String _baseDir,
73 final int _maxGenerations) {
74
75 this.distributionNumParameters = _distributionNumParameters;
76 this.deduplicator = _deduplicator;
77 this.baseDir = _baseDir;
78 this.maxGenerations = _maxGenerations;
79 }
80
81
82 public final FitnessExtractor<FitnessVector<Float>> fitnessExtractor(final int maxPossibleDistributions,
83 final double[][] samples) {
84 return (openCLExecutionContext, kernelExecutionContext, executorService, generation, genotypes,
85 resultExtractor) -> {
86
87 final float[] results = resultExtractor.extractFloatArray(openCLExecutionContext, 5);
88
89 try {
90 final var extractionTask = executorService.submit(() -> {
91 return IntStream.range(0, genotypes.size()).boxed().parallel().map(genotypeIndex -> {
92 float tally = 0;
93
94 final int startIndex = genotypeIndex * samples.length;
95 final int endIndex = startIndex + samples.length;
96 int notFiniteCount = 0;
97 for (int i = startIndex; i < endIndex; i++) {
98 if (Float.isFinite(results[i])) {
99 tally += results[i];
100 } else {
101 notFiniteCount++;
102 tally -= 10_000f;
103 }
104 }
105 final int[] assignedClusters = ClusteringUtils.assignClustersFloatChromosome(
106 distributionNumParameters,
107 samples,
108 genotypes.get(genotypeIndex));
109 final Set<Integer> uniqueAssigned = new HashSet<>();
110 uniqueAssigned.addAll(IntStream.of(assignedClusters).boxed().toList());
111
112 if (notFiniteCount > 0 || notFiniteCount == samples.length || uniqueAssigned.isEmpty()) {
113 return new FitnessVector<>(-100_000_000f, -100f);
114 }
115
116 final float numUnusedCluster = (float) (maxPossibleDistributions - uniqueAssigned.size());
117 return new FitnessVector<>(tally / samples.length, numUnusedCluster);
118 }).toList();
119 });
120
121 return extractionTask.get();
122 } catch (InterruptedException | ExecutionException e) {
123 logger.error("Could not extract results", e);
124 throw new RuntimeException(e);
125 }
126 };
127 }
128
129
130 public DataSupplier<float[]> computeSideData(final int maxPossibleDistributions) {
131 return (openCLExecutionContext, generation, genotypes) -> {
132
133
134
135
136
137
138 final float[] densityDataHelper = new float[genotypes.size() * maxPossibleDistributions * (2 + 4)];
139 int genotypeIndex = 0;
140 for (final var genotype : genotypes) {
141 final var fChromosome = genotype.getChromosome(0, FloatChromosome.class);
142
143 for (int distributionIndex = 0; distributionIndex < maxPossibleDistributions; distributionIndex++) {
144 final int baseDistributionIndex = distributionIndex * distributionNumParameters;
145
146 final double[] mean = new double[] { fChromosome.getAllele(baseDistributionIndex + 1),
147 fChromosome.getAllele(baseDistributionIndex + 2) };
148 final double[][] covariances = new double[2][2];
149 covariances[0][0] = fChromosome.getAllele(baseDistributionIndex + 3) - 15;
150 covariances[0][1] = fChromosome.getAllele(baseDistributionIndex + 4) - 15;
151 covariances[1][0] = fChromosome.getAllele(baseDistributionIndex + 4) - 15;
152 covariances[1][1] = fChromosome.getAllele(baseDistributionIndex + 5) - 15;
153
154 try {
155
156
157
158 final var multivariateNormalDistribution = new MultivariateNormalDistribution(mean, covariances);
159
160 final var covarianceMatrix = new Array2DRowRealMatrix(covariances);
161 final var eigenDecomposition = new EigenDecomposition(covarianceMatrix);
162 final var eigenDecompositionSolver = eigenDecomposition.getSolver();
163 final var covarianceInverse = eigenDecompositionSolver.getInverse();
164 final var covarianceDeterminant = eigenDecomposition.getDeterminant();
165
166 final int baseDensityDataHelperIndex = genotypeIndex * maxPossibleDistributions * 6
167 + distributionIndex * 6;
168 densityDataHelper[baseDensityDataHelperIndex + 1] = (float) covarianceDeterminant;
169 densityDataHelper[baseDensityDataHelperIndex + 2] = (float) covarianceInverse.getEntry(0, 0);
170 densityDataHelper[baseDensityDataHelperIndex + 3] = (float) covarianceInverse.getEntry(0, 1);
171 densityDataHelper[baseDensityDataHelperIndex + 4] = (float) covarianceInverse.getEntry(1, 0);
172 densityDataHelper[baseDensityDataHelperIndex + 5] = (float) covarianceInverse.getEntry(1, 1);
173
174 densityDataHelper[baseDensityDataHelperIndex + 0] = 10;
175
176 } catch (NonPositiveDefiniteMatrixException | MathUnsupportedOperationException
177 | SingularMatrixException e) {
178
179 }
180 }
181
182 genotypeIndex++;
183 }
184 return densityDataHelper;
185 };
186 }
187
188 public void run(final int maxPossibleDistributions, final int numDistributions, final double[][] samplesDouble,
189 final float[][] samples, final float[] x, final float[] y, final String algorithmName,
190 final Collection<Genotype> seedPopulation, final EvolutionResult<FitnessVector<Float>> bestCPUResult)
191 throws IOException {
192
193 final var kernelName = "mixtureModelKernel";
194 final int chromosomeSize = distributionNumParameters * maxPossibleDistributions;
195
196
197 final var singleKernelFitnessDescriptorBuilder = SingleKernelFitnessDescriptor.builder();
198 singleKernelFitnessDescriptorBuilder.kernelName(kernelName)
199 .kernelExecutionContextComputer(KernelExecutionContextComputers.ofGlobalWorkSize1D(samples.length))
200 .putStaticDataLoaders(0, StaticDataLoaders.ofLinearize(samples))
201 .putStaticDataLoaders(
202 1,
203 StaticDataLoaders
204 .of(distributionNumParameters, maxPossibleDistributions, chromosomeSize, samples.length))
205 .putDataLoaders(2, DataLoaders.ofGenerationAndPopulationSize())
206 .putDataLoaders(3, DataLoaders.ofLinearizeFloatChromosome(0))
207 .putDataLoaders(4, DataLoaders.ofFloatSupplier(computeSideData(maxPossibleDistributions)))
208 .putResultAllocators(5, ResultAllocators.ofMultiplePopulationSizeFloat(samples.length));
209 final var singleKernelFitnessDescriptor = singleKernelFitnessDescriptorBuilder.build();
210
211
212
213 final var singleKernelFitness = SingleKernelFitness
214 .of(singleKernelFitnessDescriptor, fitnessExtractor(maxPossibleDistributions, samplesDouble));
215
216 final var eaConfigurationBuilder = GPUEAConfiguration.<FitnessVector<Float>>builder()
217 .chromosomeSpecs(FloatChromosomeSpec.of(chromosomeSize, 0, 30))
218 .parentSelectionPolicy(TournamentNSGA2Selection.ofFitnessVector(2, 5, deduplicator))
219 .replacementStrategy(
220 Elitism.builder()
221 .offspringRatio(0.950)
222 .offspringSelectionPolicy(TournamentNSGA2Selection.ofFitnessVector(2, 5, deduplicator))
223 .survivorSelectionPolicy(NSGA2Selection.ofFitnessVector(2, deduplicator))
224 .build())
225 .combinationPolicy(
226 MultiCombinations.of(
227 SinglePointArithmetic.of(0.9),
228 SinglePointCrossover.build(),
229 MultiPointCrossover.of(3),
230 MultiPointArithmetic.of(3, 0.9)))
231 .mutationPolicies(
232 MultiMutations.of(
233 RandomMutation.of(0.40),
234 CreepMutation.of(0.40, NormalDistribution.of(0.0, 1)),
235 SwapMutation.of(0.40, 5, false)))
236 .program(Program.ofResource("/opencl/mixturemodel/main.cl", kernelName, ""))
237 .fitness(singleKernelFitness)
238 .termination(Terminations.ofMaxGeneration(maxGenerations));
239
240 if (CollectionUtils.isNotEmpty(seedPopulation)) {
241 eaConfigurationBuilder.seedPopulation(seedPopulation);
242 }
243
244 final var eaConfiguration = eaConfigurationBuilder.build();
245
246
247
248 final var csvEvolutionListener = CSVEvolutionListener.<FitnessVector<Float>, Void>of(
249 baseDir + algorithmName + ".csv",
250 List.of(
251 ColumnExtractor.of("generation", e -> e.generation()),
252 ColumnExtractor.of("fitness", e -> e.fitness().get(0)),
253 ColumnExtractor.of("numUnusedCluster", e -> e.fitness().get(1)),
254 ColumnExtractor.of("numCluster", e -> maxPossibleDistributions - e.fitness().get(1))));
255
256 final var topNloggerEvolutionListener = EvolutionListeners.<FitnessVector<Float>>ofLogTopN(
257 logger,
258 5,
259 Comparator.<FitnessVector<Float>, Float>comparing(fv -> fv.get(0)));
260
261 final var eaExecutionContextBuilder = GPUEAExecutionContext.<FitnessVector<Float>>builder()
262 .deviceFilters(DeviceFilters.ofGPU())
263 .populationSize(250)
264 .addEvolutionListeners(csvEvolutionListener, topNloggerEvolutionListener);
265
266
267
268
269 eaExecutionContextBuilder.addSelectionPolicyHandlerFactories(
270 gsd -> new NSGA2SelectionPolicyHandler<>(),
271 gsd -> new TournamentNSGA2SelectionPolicyHandler<>(gsd.randomGenerator()));
272
273 eaExecutionContextBuilder.addReplacementStrategyHandlerFactories(gsd -> new SPEA2ReplacementStrategyHandler<>());
274
275 final var eaExecutionContext = eaExecutionContextBuilder.build();
276
277
278 final var executorService = Executors.newWorkStealingPool();
279 final var eaSystem = GPUEASystemFactory.from(eaConfiguration, eaExecutionContext, executorService);
280
281 final EvolutionResult<FitnessVector<Float>> evolutionResult = eaSystem.evolve();
282 logger.info("Best genotype: {}", evolutionResult.bestGenotype());
283 logger.info(" with fitness: {}", evolutionResult.bestFitness());
284 logger.info(" at generation: {}", evolutionResult.generation());
285
286 ClusteringUtils.categorizeByNumClusters(
287 distributionNumParameters,
288 maxPossibleDistributions,
289 x,
290 y,
291 samplesDouble,
292 evolutionResult,
293 baseDir,
294 algorithmName);
295 }
296 }