1 package net.bmahe.genetics4j.samples.symbolicregression;
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.Optional;
8 import java.util.Random;
9 import java.util.stream.IntStream;
10
11 import org.apache.commons.cli.CommandLine;
12 import org.apache.commons.cli.CommandLineParser;
13 import org.apache.commons.cli.DefaultParser;
14 import org.apache.commons.cli.HelpFormatter;
15 import org.apache.commons.cli.Options;
16 import org.apache.commons.cli.ParseException;
17 import org.apache.commons.io.FileUtils;
18 import org.apache.commons.lang3.StringUtils;
19 import org.apache.commons.lang3.Validate;
20 import org.apache.commons.lang3.time.DurationFormatUtils;
21 import org.apache.logging.log4j.LogManager;
22 import org.apache.logging.log4j.Logger;
23
24 import net.bmahe.genetics4j.core.EASystem;
25 import net.bmahe.genetics4j.core.EASystemFactory;
26 import net.bmahe.genetics4j.core.Fitness;
27 import net.bmahe.genetics4j.core.Genotype;
28 import net.bmahe.genetics4j.core.chromosomes.TreeChromosome;
29 import net.bmahe.genetics4j.core.evolutionlisteners.EvolutionListener;
30 import net.bmahe.genetics4j.core.evolutionlisteners.EvolutionListeners;
31 import net.bmahe.genetics4j.core.spec.EAConfiguration;
32 import net.bmahe.genetics4j.core.spec.EAExecutionContext;
33 import net.bmahe.genetics4j.core.spec.EvolutionResult;
34 import net.bmahe.genetics4j.core.spec.Optimization;
35 import net.bmahe.genetics4j.core.spec.mutation.MultiMutations;
36 import net.bmahe.genetics4j.core.termination.Terminations;
37 import net.bmahe.genetics4j.gp.Operation;
38 import net.bmahe.genetics4j.gp.math.SimplificationRules;
39 import net.bmahe.genetics4j.gp.program.Program;
40 import net.bmahe.genetics4j.gp.spec.GPEAExecutionContexts;
41 import net.bmahe.genetics4j.gp.spec.chromosome.ProgramTreeChromosomeSpec;
42 import net.bmahe.genetics4j.gp.spec.combination.ProgramRandomCombine;
43 import net.bmahe.genetics4j.gp.spec.mutation.NodeReplacement;
44 import net.bmahe.genetics4j.gp.spec.mutation.ProgramApplyRules;
45 import net.bmahe.genetics4j.gp.spec.mutation.ProgramRandomMutate;
46 import net.bmahe.genetics4j.gp.spec.mutation.ProgramRandomPrune;
47 import net.bmahe.genetics4j.gp.utils.ProgramUtils;
48 import net.bmahe.genetics4j.gp.utils.TreeNodeUtils;
49 import net.bmahe.genetics4j.moo.FitnessVector;
50 import net.bmahe.genetics4j.moo.nsga2.spec.TournamentNSGA2Selection;
51 import net.bmahe.genetics4j.moo.spea2.spec.replacement.SPEA2Replacement;
52
53 public class SymbolicRegressionWithMOOSPEA2 {
54 final static public Logger logger = LogManager.getLogger(SymbolicRegressionWithMOOSPEA2.class);
55
56 final static public String PARAM_DEST_CSV = "d";
57 final static public String LONG_PARAM_DEST_CSV = "csv-dest";
58
59 final static public String PARAM_POPULATION_SIZE = "p";
60 final static public String LONG_PARAM_POPULATION_SIZE = "population-size";
61
62 final static public String DEFAULT_DEST_CSV = SymbolicRegressionWithMOOSPEA2.class.getSimpleName() + ".csv";
63
64 final static public int DEFAULT_POPULATION_SIZE = 500;
65
66 public static void cliError(final Options options, final String errorMessage) {
67 final HelpFormatter formatter = new HelpFormatter();
68 logger.error(errorMessage);
69 formatter.printHelp(SymbolicRegressionWithMOOSPEA2.class.getSimpleName(), options);
70 System.exit(-1);
71 }
72
73 @SuppressWarnings("unchecked")
74 public void run(String csvFilename, int populationSize) {
75 Validate.isTrue(StringUtils.isNotBlank(csvFilename));
76 Validate.isTrue(populationSize > 0);
77
78 final Random random = new Random();
79
80 final Program program = SymbolicRegressionUtils.buildProgram(random);
81
82 final Comparator<Genotype> deduplicator = (a, b) -> TreeNodeUtils.compare(a, b, 0);
83
84
85 final Fitness<FitnessVector<Double>> computeFitness = (genoType) -> {
86 final TreeChromosome<Operation<?>> chromosome = (TreeChromosome<Operation<?>>) genoType.getChromosome(0);
87 final Double[][] inputs = new Double[100][1];
88 for (int i = 0; i < 100; i++) {
89 inputs[i][0] = (i - 50) * 1.2;
90 }
91
92 double mse = 0;
93 for (final Double[] input : inputs) {
94
95 final double x = input[0];
96 final double expected = SymbolicRegressionUtils.evaluate(x);
97 final Object result = ProgramUtils.execute(chromosome, input);
98
99 if (Double.isFinite(expected)) {
100 final Double resultDouble = (Double) result;
101 if (Double.isFinite(resultDouble)) {
102 mse += (expected - resultDouble) * (expected - resultDouble);
103 } else {
104 mse += 1_000_000_000;
105 }
106 }
107 }
108
109 return Double.isFinite(mse) ? new FitnessVector<Double>(mse / 100.0, (double) chromosome.getRoot().getSize())
110 : new FitnessVector<Double>(Double.MAX_VALUE, Double.MAX_VALUE);
111 };
112
113
114
115 final var eaConfigurationBuilder = new EAConfiguration.Builder<FitnessVector<Double>>();
116 eaConfigurationBuilder.chromosomeSpecs(ProgramTreeChromosomeSpec.of(program))
117 .parentSelectionPolicy(TournamentNSGA2Selection.ofFitnessVector(2, 3, deduplicator))
118 .replacementStrategy(SPEA2Replacement.ofFitnessVector(deduplicator))
119 .combinationPolicy(ProgramRandomCombine.build())
120 .mutationPolicies(
121 MultiMutations.of(
122 ProgramRandomMutate.of(0.15 * 3),
123 ProgramRandomPrune.of(0.15 * 3),
124 NodeReplacement.of(0.15 * 3)),
125 ProgramApplyRules.of(SimplificationRules.SIMPLIFY_RULES))
126 .optimization(Optimization.MINIMIZE)
127 .termination(
128 Terminations.or(
129 Terminations.<FitnessVector<Double>>ofMaxGeneration(200),
130 (eaConfiguration, generation, population, fitness) -> fitness.stream()
131 .anyMatch(fv -> fv.get(0) <= 0.00001 && fv.get(1) <= 20)))
132 .fitness(computeFitness);
133 final EAConfiguration<FitnessVector<Double>> eaConfiguration = eaConfigurationBuilder.build();
134
135
136 final var eaExecutionContextBuilder = GPEAExecutionContexts.<FitnessVector<Double>>forGP(random);
137 eaExecutionContextBuilder.populationSize(populationSize);
138 eaExecutionContextBuilder.numberOfPartitions(Math.max(1, Runtime.getRuntime().availableProcessors() - 3));
139
140 eaExecutionContextBuilder.addEvolutionListeners(
141 EvolutionListeners.ofLogTopN(
142 logger,
143 5,
144 Comparator.<FitnessVector<Double>, Double>comparing(fv -> fv.get(0)).reversed(),
145 (genotype) -> TreeNodeUtils.toStringTreeNode(genotype, 0)),
146 SymbolicRegressionUtils.csvLogger(
147 csvFilename,
148 evolutionStep -> evolutionStep.fitness().get(0),
149 evolutionStep -> evolutionStep.fitness().get(1)),
150 new EvolutionListener<FitnessVector<Double>>() {
151
152 long startTime = -1;
153 long previousTime = -1;
154
155 @Override
156 public void onEvolution(final long generation, final List<Genotype> population,
157 final List<FitnessVector<Double>> fitness, final boolean isDone) {
158 final long now = System.currentTimeMillis();
159
160 if (startTime < 0) {
161 startTime = now;
162 }
163
164 if (previousTime < 0) {
165 previousTime = now;
166 }
167
168 if (generation > 1) {
169 logger.info("Execution time:");
170 logger.info(
171 "\tCurrent generation duration: {}",
172 DurationFormatUtils.formatDurationHMS(now - previousTime));
173 logger.info(
174 "\tAverage duration: {}",
175 DurationFormatUtils.formatDurationHMS((now - startTime) / (generation - 1)));
176 }
177
178 previousTime = now;
179 }
180 });
181
182 final EAExecutionContext<FitnessVector<Double>> eaExecutionContext = eaExecutionContextBuilder.build();
183 final EASystem<FitnessVector<Double>> eaSystem = EASystemFactory.from(eaConfiguration, eaExecutionContext);
184
185 final EvolutionResult<FitnessVector<Double>> evolutionResult = eaSystem.evolve();
186 final Genotype bestGenotype = evolutionResult.bestGenotype();
187 final TreeChromosome<Operation<?>> bestChromosome = (TreeChromosome<Operation<?>>) bestGenotype.getChromosome(0);
188 logger.info("Best genotype: {}", bestChromosome.getRoot());
189 logger.info("Best genotype - pretty print: {}", TreeNodeUtils.toStringTreeNode(bestChromosome.getRoot()));
190
191 final int depthIdx = 1;
192 for (int i = 0; i < 15; i++) {
193 final int depth = i;
194 final Optional<Integer> optIdx = IntStream.range(0, evolutionResult.fitness().size())
195 .boxed()
196 .filter((idx) -> evolutionResult.fitness().get(idx).get(depthIdx) == depth)
197 .sorted(
198 (a, b) -> Double
199 .compare(evolutionResult.fitness().get(a).get(0), evolutionResult.fitness().get(b).get(0)))
200 .findFirst();
201
202 optIdx.stream().forEach((idx) -> {
203 final TreeChromosome<Operation<?>> treeChromosome = (TreeChromosome<Operation<?>>) evolutionResult
204 .population()
205 .get(idx)
206 .getChromosome(0);
207
208 logger.info(
209 "Best genotype for depth {} - score {} -> {}",
210 depth,
211 evolutionResult.fitness().get(idx).get(0),
212 TreeNodeUtils.toStringTreeNode(treeChromosome.getRoot()));
213 });
214 }
215 }
216
217 public static void main(String[] args) throws IOException {
218
219
220
221
222
223 final CommandLineParser parser = new DefaultParser();
224
225 final Options options = new Options();
226 options.addOption(PARAM_DEST_CSV, LONG_PARAM_DEST_CSV, true, "destination csv file");
227
228 options.addOption(PARAM_POPULATION_SIZE, LONG_PARAM_POPULATION_SIZE, true, "Population size");
229
230 String csvFilename = DEFAULT_DEST_CSV;
231 int populationSize = DEFAULT_POPULATION_SIZE;
232 try {
233 final CommandLine line = parser.parse(options, args);
234
235 if (line.hasOption(PARAM_DEST_CSV)) {
236 csvFilename = line.getOptionValue(PARAM_DEST_CSV);
237 }
238
239 if (line.hasOption(PARAM_POPULATION_SIZE)) {
240 populationSize = Integer.parseInt(line.getOptionValue(PARAM_POPULATION_SIZE));
241 }
242
243 } catch (ParseException exp) {
244 cliError(options, "Unexpected exception:" + exp.getMessage());
245 }
246
247 logger.info("Population size: {}", populationSize);
248
249 logger.info("CSV output located at {}", csvFilename);
250 FileUtils.forceMkdirParent(new File(csvFilename));
251
252 final var symbolicRegression = new SymbolicRegressionWithMOOSPEA2();
253 symbolicRegression.run(csvFilename, populationSize);
254 }
255 }