1 package net.bmahe.genetics4j.samples.mixturemodel;
2
3 import java.io.File;
4 import java.io.IOException;
5 import java.util.Comparator;
6 import java.util.List;
7 import java.util.Objects;
8 import java.util.random.RandomGenerator;
9
10 import org.apache.commons.cli.CommandLine;
11 import org.apache.commons.cli.CommandLineParser;
12 import org.apache.commons.cli.DefaultParser;
13 import org.apache.commons.cli.Options;
14 import org.apache.commons.cli.ParseException;
15 import org.apache.commons.io.FileUtils;
16 import org.apache.commons.lang3.Validate;
17 import org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution;
18 import org.apache.commons.math3.distribution.MultivariateNormalDistribution;
19 import org.apache.commons.math3.distribution.fitting.MultivariateNormalMixtureExpectationMaximization;
20 import org.apache.commons.math3.linear.RealMatrix;
21 import org.apache.commons.math3.util.Pair;
22 import org.apache.logging.log4j.LogManager;
23 import org.apache.logging.log4j.Logger;
24
25 import net.bmahe.genetics4j.core.Genotype;
26 import net.bmahe.genetics4j.core.chromosomes.FloatChromosome;
27 import net.bmahe.genetics4j.samples.CLIUtils;
28
29 public class Main {
30 public static final Logger logger = LogManager.getLogger(Main.class);
31
32 public static final String OPT_BASE_DIR = "base-directory";
33
34 public static final int NUM_POINTS_PER_DISTRIBUTION = 2_000;
35 public static final int MIN_NUM_DISTRIBUTIONS = 3;
36 public static final int MAX_NUM_DISTRIBUTIONS = 6;
37
38 public static final float MIN_MEAN = 0F;
39 public static final float MAX_MEAN = 30F;
40
41 public static final float MIN_STD_DEV = 0.5F;
42 public static final float MAX_STD_DEV = 5F;
43
44 public static final int DISTRIBUTION_NUM_PARAMETERS = 6;
45
46 public static final int MAX_GENERATIONS = 1_000;
47
48 public static final double EPSILON = 0.0001;
49
50 public static final Comparator<Genotype> deduplicator = (a, b) -> {
51 final var aChromosome = a.getChromosome(0, FloatChromosome.class);
52 final var bChromosome = b.getChromosome(0, FloatChromosome.class);
53
54 if (aChromosome.getSize() != bChromosome.getSize()) {
55 return Integer.compare(a.getSize(), b.getSize());
56 }
57
58 for (int i = 0; i < aChromosome.getSize(); i++) {
59 final var alleleA = aChromosome.getAllele(i);
60 final var alleleB = bChromosome.getAllele(i);
61
62 final var diff = alleleB - alleleA;
63 if (diff > EPSILON) {
64 return 1;
65 } else if (diff < -EPSILON) {
66 return -1;
67 }
68 }
69
70 return 0;
71 };
72
73 public static Genotype toGenotype(final int maxPossibleDistributions,
74 final MixtureMultivariateNormalDistribution mixtureMultivariateNormalDistribution) {
75 Objects.requireNonNull(mixtureMultivariateNormalDistribution);
76
77 final float[] values = new float[maxPossibleDistributions * DISTRIBUTION_NUM_PARAMETERS];
78
79 int i = 0;
80 final var components = mixtureMultivariateNormalDistribution.getComponents();
81 Validate.isTrue(components.size() <= maxPossibleDistributions);
82
83 for (final Pair<Double, MultivariateNormalDistribution> component : components) {
84
85 final var alpha = component.getFirst();
86 final var multivariateNormalDistribution = component.getSecond();
87
88 final double[] means = multivariateNormalDistribution.getMeans();
89 final RealMatrix covariances = multivariateNormalDistribution.getCovariances();
90
91 values[i + 0] = alpha.floatValue();
92 values[i + 1] = (float) means[0];
93 values[i + 2] = (float) means[1];
94 values[i + 3] = (float) covariances.getEntry(0, 0) + 15;
95 values[i + 4] = (float) covariances.getEntry(0, 1) + 15;
96 values[i + 5] = (float) covariances.getEntry(1, 1) + 15;
97
98 i += DISTRIBUTION_NUM_PARAMETERS;
99 }
100
101 return new Genotype(
102 new FloatChromosome(maxPossibleDistributions * DISTRIBUTION_NUM_PARAMETERS, 0, MAX_MEAN, values));
103 }
104
105 public static void main(String[] args) throws IOException {
106 logger.info("Starting");
107
108
109
110
111 logger.info("Parsing command line arguments");
112 final CommandLineParser parser = new DefaultParser();
113
114 final Options options = new Options();
115 options.addOption("h", "help", false, "print help");
116 options.addOption(null, "num-clusters", true, "number of clusters");
117 options.addOption(null, OPT_BASE_DIR, true, "base directory for output files");
118
119 String baseDir = "./";
120 try {
121 final CommandLine line = parser.parse(options, args);
122 if (line.hasOption("h")) {
123 CLIUtils.cliHelpAndExit(logger, Main.class, options, null);
124 }
125
126 if (line.hasOption(OPT_BASE_DIR)) {
127 baseDir = line.getOptionValue(OPT_BASE_DIR);
128 if (baseDir.endsWith("/") == false) {
129 baseDir = baseDir + "/";
130 }
131 }
132 } catch (ParseException exp) {
133 CLIUtils.cliHelpAndExit(logger, Main.class, options, "Unexpected exception:" + exp.getMessage());
134
135
136 throw new RuntimeException();
137
138 }
139
140 logger.info("Using base directory for files: {}", baseDir);
141 FileUtils.forceMkdir(new File(baseDir));
142
143 logger.info("Starting process");
144 final RandomGenerator randomGenerator = RandomGenerator.getDefault();
145 for (int i = 0; i < 100; i++) {
146 randomGenerator.nextInt();
147 }
148
149
150 final int maxPossibleDistributions = MAX_NUM_DISTRIBUTIONS - 1;
151 logger.info("Generating dataset");
152 final int numDistributions = randomGenerator.nextInt(MIN_NUM_DISTRIBUTIONS, MAX_NUM_DISTRIBUTIONS);
153 logger.info("Generating {} distributions:", numDistributions);
154
155 final double[][] means = new double[numDistributions][2];
156 final double[][][] covariances = new double[numDistributions][2][2];
157
158 final int[] cluster = new int[NUM_POINTS_PER_DISTRIBUTION * numDistributions];
159 final float[][] samples = new float[NUM_POINTS_PER_DISTRIBUTION * numDistributions][2];
160 final double[][] samplesDouble = new double[NUM_POINTS_PER_DISTRIBUTION * numDistributions][2];
161 for (int i = 0; i < numDistributions; i++) {
162 means[i] = new double[] { randomGenerator.nextDouble(MIN_MEAN, MAX_MEAN),
163 randomGenerator.nextDouble(MIN_MEAN, MAX_MEAN) };
164
165
166
167
168
169 final double[][] randomMatrix = new double[][] {
170 { randomGenerator.nextDouble(0, 4) - 2, randomGenerator.nextDouble(0, 4) - 2 },
171 { randomGenerator.nextDouble(0, 4) - 2, randomGenerator.nextDouble(0, 4) - 2 } };
172 final double[][] covariance = new double[][] {
173 { randomMatrix[0][0] * randomMatrix[0][0] + randomMatrix[0][1] * randomMatrix[0][1],
174 randomMatrix[0][0] * randomMatrix[1][0] + randomMatrix[0][1] * randomMatrix[1][1] },
175 { randomMatrix[0][0] * randomMatrix[1][0] + randomMatrix[0][1] * randomMatrix[1][1],
176 randomMatrix[1][0] * randomMatrix[1][0] + randomMatrix[1][1] * randomMatrix[1][1] } };
177 covariances[i] = covariance;
178
179 logger.info("\t index: {} - mean: {} - covariance: {}", i, means[i], covariance);
180 final var multivariateNormalDistribution = new MultivariateNormalDistribution(means[i], covariance);
181 for (int j = 0; j < 5; j++) {
182 final var sample = multivariateNormalDistribution.sample();
183 logger.info("\t\t{} - {}", sample, multivariateNormalDistribution.density(sample));
184 }
185 for (int j = 0; j < NUM_POINTS_PER_DISTRIBUTION; j++) {
186 final var sample = multivariateNormalDistribution.sample();
187 samples[i * NUM_POINTS_PER_DISTRIBUTION + j] = new float[] { (float) sample[0], (float) sample[1] };
188 samplesDouble[i * NUM_POINTS_PER_DISTRIBUTION + j] = sample;
189 cluster[i * NUM_POINTS_PER_DISTRIBUTION + j] = i;
190 }
191 }
192
193 final float[] x = new float[NUM_POINTS_PER_DISTRIBUTION * numDistributions];
194 final float[] y = new float[NUM_POINTS_PER_DISTRIBUTION * numDistributions];
195 for (int i = 0; i < NUM_POINTS_PER_DISTRIBUTION * numDistributions; i++) {
196 x[i] = samples[i][0];
197 y[i] = samples[i][1];
198 }
199
200 ClusteringUtils.persistClusters(x, y, cluster, baseDir + "original.csv");
201
202
203 final var singleObjectiveMethod = new SingleObjectiveMethod(DISTRIBUTION_NUM_PARAMETERS,
204 baseDir,
205 MAX_GENERATIONS);
206 singleObjectiveMethod.run(maxPossibleDistributions, samplesDouble, x, y, baseDir, List.of());
207
208 final var mooCPU = new MooCPU(DISTRIBUTION_NUM_PARAMETERS, deduplicator, baseDir, MAX_GENERATIONS);
209 final long startCPU = System.currentTimeMillis();
210 final var bestCPUResult = mooCPU.run(maxPossibleDistributions, samplesDouble, x, y, "GA", List.of());
211 final long endCPU = System.currentTimeMillis();
212 final long durationCPUMs = endCPU - startCPU;
213
214 final var mooGPU = new MooGPU(DISTRIBUTION_NUM_PARAMETERS, deduplicator, baseDir, MAX_GENERATIONS);
215 final long startGPU = System.currentTimeMillis();
216 mooGPU.run(
217 maxPossibleDistributions,
218 numDistributions,
219 samplesDouble,
220 samples,
221 x,
222 y,
223 "mixturemodel-moo-gpu",
224 List.of(),
225 bestCPUResult);
226 final long endGPU = System.currentTimeMillis();
227 final long durationGPUMs = endGPU - startGPU;
228
229 logger.info("Duration CPU: {} s - Duration GPU: {} s", durationCPUMs / 1000, durationGPUMs / 1000);
230
231
232 logger.info("Evaluating with commons math");
233 final var mmmmm = new MultivariateNormalMixtureExpectationMaximization(samplesDouble);
234
235 final var initialEstimatedDistribution = MultivariateNormalMixtureExpectationMaximization
236 .estimate(samplesDouble, numDistributions);
237 mmmmm.fit(initialEstimatedDistribution);
238 final var mixtureMultivariateNormalDistribution = mmmmm.getFittedModel();
239 final var components = mixtureMultivariateNormalDistribution.getComponents();
240 final var fittedCommonsMathGenotype = toGenotype(maxPossibleDistributions, mixtureMultivariateNormalDistribution);
241
242
243 logger.info("GPU run with best seeds from commons math results");
244 mooGPU.run(
245 maxPossibleDistributions,
246 numDistributions,
247 samplesDouble,
248 samples,
249 x,
250 y,
251 "mixturemodel-moo-gpu-seed",
252 List.of(fittedCommonsMathGenotype),
253 bestCPUResult);
254
255 logger.info("Best results from apache commons math:");
256 final int[] mnmeAssigned = new int[samples.length];
257 for (int i = 0; i < samples.length; i++) {
258 double bestScore = Double.NEGATIVE_INFINITY;
259
260 for (int c = 0; c < components.size(); c++) {
261 final var component = components.get(c);
262 final var alpha = component.getFirst();
263 final var distribution = component.getSecond();
264
265 final var score = alpha * distribution.density(samplesDouble[i]);
266 if (score > bestScore) {
267 bestScore = score;
268 mnmeAssigned[i] = c;
269 }
270 }
271 }
272
273 logger.info("Components:");
274 for (int c = 0; c < components.size(); c++) {
275 final var component = components.get(c);
276 final var alpha = component.getFirst();
277 final var distribution = component.getSecond();
278
279 logger.debug(
280 "\talpha: {} - means: {} - covariances: {}",
281 alpha,
282 distribution.getMeans(),
283 distribution.getCovariances());
284 }
285
286 logger.info("Commons-math: log likelyhood: {}", mmmmm.getLogLikelihood());
287 ClusteringUtils.persistClusters(x, y, mnmeAssigned, baseDir + "assigned-commons-math.csv");
288
289 final var fitnessCPUFunc = mooCPU.fitnessCPU(maxPossibleDistributions, samplesDouble);
290 final var mnmeFitness = fitnessCPUFunc.compute(fittedCommonsMathGenotype);
291
292 ClusteringUtils.writeCSVReferenceValue(baseDir + "commons-math.csv", MAX_GENERATIONS, mnmeFitness.get(0));
293 logger.info("Score of the Commons Math results: {}, comparing to {}", mnmeFitness, mmmmm.getLogLikelihood());
294 logger.info("Best CPU evolution: {}", bestCPUResult.bestFitness());
295 }
296 }