Merge pull request #15 from samthiriot/test_functions1

add test functions for NSGA2

README.md

@@ -15,6 +15,8 @@ Still they should be easily reusable and any material required by the workflow (
 
 All the accepted submissions will be part of the market place and will be compiled against the last version of OpenMOLE. The available entries are also available from the OpenMOLE website: [http://next.openmole.org/Documentation_Market%20Place_All.html](http://next.openmole.org/Documentation_Market%20Place_All.html)
 
+Note: if you are asked for a password for "gitlab.openmole.org" when pushing your changes into your fork of the repository, you might simply add the `--no-verify` flag to your git push.
+
 ## Available workflows ##
 
   - [Advanced methods](https://github.com/guillaumecherel/TutorialEAForModelling): advanced methods for calibrating, validating and analyzing complex systems models.
@@ -30,5 +32,7 @@ All the accepted submissions will be part of the market place and will be compil
   - [R hello world](R-hello): an example of how to embed R code in OpenMole.  This workflow executes an R program with 100 different inputs, makes a computation, and saves to a file.
   - [Random Forest classifier](randomforest): This workflow explores the parameters of a random forest image classifier written in Python using scikit-learn.
   - [SimPopLocal model](simpoplocal): a geographical model calibrated using genetic algorithms.
-  - [Sensitivity-Screening analysis](sensitivity/morris): a method to quickly analyze which inputs are influential on large spaces of parameters.
-  - [Global Sensitivity Analysis](sensitivity/saltelli): a variance based sensitivity analysis of model output.
+  - [Sensitivity-Screening analysis](sensitivity/morris): a method to quickly analyze which inputs are influential on large spaces of parameters.
+  - [Global Sensitivity Analysis](sensitivity/saltelli): a variance based sensitivity analysis of model output
+  - [Test Functions for NSGA2](nsga2-test-functions): reference functions to double check the correctedness of the NSGA2 algorithm, and also view examples of usage of NSGA2
+  

nsga2-test-functions/README.md

+# Test Functions for NSGA2 
+
+## Why Test? 
+
+Metaheuristics for optimization such as genetic algorithms leave a huge freedom to developers for their implementation.
+Genetic Algorithms such as NSGA2 are even less trivial. 
+Their implementation in the [mgo](https://github.com/openmole/mgo) used in [OpenMole](https://openmole.org/) was adapted so they are efficient for computation on clusters or grids. 
+This leaves a lot of uncertainty to the user on what really happens.
+Why should you, as a user, trust this implementation?
+The answer is simple: you should not.
+Any implementation in scientific computing should be verified to be trusted.
+
+
+## Test Functions
+
+Many functions named "Test Functions" were proposed over time to test and compare multi-objective optimization algorithms.
+The [wikipedia page](https://en.wikipedia.org/wiki/Test_functions_for_optimization) lists several of them.
+
+A test function is an function  you can use as the problem to optimize, for which the expected results are known.
+Test functions might test:
+  - problems with multiple parameters and multiple objectives
+  - problems with constraints: a part of the space of solutions can not be computed, so the algorithm should achieve to deal with it
+  - problems with parts of the space of solution that are tricky to detect, either because they are statistically unlikely to find, or because they are in a part of the space of solutions which is heavily constrained, etc.
+
+The results of a test function are usually compared in both terms of:
+  - coverage of the Pareto front,
+  - speed of convergence
+
+
+## Test and Learn 
+
+Test functions have another role: they constitute a simple example of optimization problem.
+As a user, you can also learn to tune the parameters of the genetic algorithm as check when convergence 
+occurs.
+
+
+## Test Workflows for OpenMole
+
+We provide here a few examples of test functions which you can open with OpenMole. 
+
+  - ConstrEx
+  - CP1
+  - Schaffer N1
+  - Schaffer N2
+
+
+To run a test 
+  - Choose one of the workflows starting with "test function"
+  - Run the workflow
+  - Update the view on the left using the "refresh" button
+  - Download the graph of the last Pareto front, and compare it with the literature (you might use the [wikipedia page](https://en.wikipedia.org/wiki/Test_functions_for_optimization) as a start)
+
+You might then tune the parameters of the optimization algorithm and analyze the results, to understand how to better use the optimization method.
+

nsga2-test-functions/execEnvironment.oms

+
+// about the local computer 
+val threads = 4
+
+val envMultiThread = LocalEnvironment(threads, name="multithread") 
\ No newline at end of file

nsga2-test-functions/plottingFunctions.oms

+
+val last_pareto = Val[File]
+val directoryWithResults = Val[File]
+val filesHaveHeaders = Val[Int]
+val countInputs = Val[Int]
+
+val taskPlotLastParetoFront = RTask("""
+    library(ggplot2)
+    
+    colnames <- if (countInputs == 2) c("iteration", "x", "y", "f1", "f2") else c("iteration", "x", "f1", "f2")
+    coltypes <- if (countInputs == 2) c("integer", "numeric", "numeric", "numeric", "numeric") else c("integer", "numeric", "numeric", "numeric") 
+    names(coltypes) <- colnames
+
+    directoryWithResultsName <- "mydirectory"
+    
+    # ensure check the directory exists
+    if (!file.exists(directoryWithResultsName)) { stop(paste("ERROR: the directory", directoryWithResultsName, "does not exists!")) }
+    
+    # get the most recent file (will be the last result)
+    allfiles <- file.info(list.files(directoryWithResultsName, full.names = T))
+    lastfilename <- rownames(allfiles)[which.max(allfiles$mtime)]
+    if (!file.exists(lastfilename)) { stop(paste("ERROR: no file found in", directoryWithResultsName) ) }
+    print(lastfilename)
+    
+    # read the last file 
+    pop <- read.csv(header = filesHaveHeaders>0, col.names=colnames, colClasses=coltypes, file=lastfilename)
+    print(paste("there are", nrow(pop), "points on the last Pareto front"))
+    
+    # plot
+    g <- ggplot(pop, aes(x=f1,y=f2)) + geom_point()
+    ggsave(filename="/tmp/last_pareto.png", plot=g) 
+    
+    """,
+    install = Seq(
+        // install the libs required for the compilation of R packages
+        "fakeroot apt-get install -y libssl-dev libcurl4-openssl-dev libudunits2-dev",
+        // install required R packages in their binary version (quicker, much stable!)
+        "fakeroot apt-get install -y r-cran-ggplot2",
+        ),
+    libraries = Seq()
+) set (
+    inputFiles += (directoryWithResults, "mydirectory"),
+    outputFiles += ("/tmp/last_pareto.png", last_pareto),
+    inputs += filesHaveHeaders.mapped,
+    inputs += countInputs.mapped,
+    (inputs, outputs) += directoryWithResults,
+    filesHaveHeaders := 1,
+    countInputs := 2
+)
+  

nsga2-test-functions/test function Binh Korn.oms

+import _file_.plottingFunctions._
+import _file_.execEnvironment._
+
+//model inputs
+val x = Val[Double]
+val y = Val[Double]
+
+//model outputs
+val f1 = Val[Double]
+val f2 = Val[Double]
+
+// about the current experiment
+val relativePath = "results/BinhKorn"
+
+// the test function 
+val testFunctionBinhKorn =
+    ScalaTask("""
+        val g1 = Math.pow(x - 5, 2) + Math.pow(y, 2)
+        val g2 = Math.pow(x - 8, 2) + Math.pow(y + 2, 3)
+        val constrained = g1 > 25 || g2 < 7.7
+        
+        val f1 = if (constrained) Double.NaN else 4*Math.pow(x,2) + 4*Math.pow(y,2)
+        val f2 = if (constrained) Double.NaN else Math.pow(x-5,2) + Math.pow(y-5,2)
+        """) set (
+    inputs += x,
+    inputs += y,
+    outputs += f1,
+    outputs += f2
+    )
+    
+// the optimisation algorithm under test 
+val evolutionBinhKorn =
+  NSGA2Evolution(
+    genome = Seq(
+      x in (0.0, 5.0),
+      y in (0.0, 3.0)
+    ),
+    objective = Seq(f1, f2),
+    evaluation = testFunctionBinhKorn,
+    parallelism = threads * 2,
+    termination = 1000
+  )
+
+// compute evolution on the test Function
+(evolutionBinhKorn on envMultiThread hook SavePopulationHook(evolutionBinhKorn, workDirectory/relativePath)
+// then plot the last Pareto front
+) -- (taskPlotLastParetoFront set ( directoryWithResults := workDirectory/relativePath) hook CopyFileHook(last_pareto, workDirectory/"last Pareto front BinhKorn.png" ) )
+

nsga2-test-functions/test function ConstrEx.oms

+import _file_.plottingFunctions._
+import _file_.execEnvironment._
+
+//model inputs
+val x = Val[Double]
+val y = Val[Double]
+
+//model outputs
+val f1 = Val[Double]
+val f2 = Val[Double]
+
+// about the current experiment
+val relativePath = "results/ConstrEx"
+
+// the test function 
+val testFunctionConstrEx =
+    ScalaTask("""
+        val g1 = y + 9*x 
+        val g2 = -y + 9*x
+        val constrained = g1 < 6 || g2 < 1
+
+        val f1 = if (constrained) Double.NaN else x;
+        val f2 = if (constrained) Double.NaN else (1 + y) / x;        
+        """) set (
+    inputs += x, 
+    inputs += y,
+    outputs += f1, 
+    outputs += f2
+    )
+
+// the optimisation algorithm under test 
+val evolutionConstrEx =
+  NSGA2Evolution(
+    genome = Seq(
+      x in (0.1, 1.0),
+      y in (0.0, 5.0)
+    ),
+    objective = Seq(f1, f2),
+    evaluation = testFunctionConstrEx, 
+    parallelism = threads * 2,
+    termination = 5000
+  )
+
+// compute evolution on the test Function
+(evolutionConstrEx on envMultiThread hook SavePopulationHook(evolutionConstrEx, workDirectory/relativePath)
+// then plot the last Pareto front
+) -- (taskPlotLastParetoFront set ( directoryWithResults := workDirectory/relativePath)  hook CopyFileHook(last_pareto, workDirectory/"last Pareto front ConstrEx.png" ) )

nsga2-test-functions/test function SchafferN1.oms

+import _file_.plottingFunctions._
+import _file_.execEnvironment._
+
+//model inputs
+val x = Val[Double]
+
+//model outputs
+val f1 = Val[Double]
+val f2 = Val[Double]
+
+// explored space
+val A = 1000.0 // up to 10^5
+
+// about the current experiment
+val relativePath = "results/SchafferN1"
+
+// the test function 
+val testFunctionSchafferN1 =
+    ScalaTask("""
+        val f1 = Math.pow(x, 2);
+        val f2 = Math.pow(x - 2, 2);        
+        """) set (
+    inputs += x,
+    outputs += f1,
+    outputs += f2
+    )
+    
+// the optimisation algorithm under test 
+val evolutionSchafferN1 =
+  NSGA2Evolution(
+    genome = Seq(
+      x in (-A, A)
+    ),
+    objective = Seq(f1, f2),
+    evaluation = testFunctionSchafferN1, 
+    parallelism = threads * 2,
+    termination = 1000
+  )
+  
+// compute evolution on the test Function
+(evolutionSchafferN1 on envMultiThread hook SavePopulationHook(evolutionSchafferN1, workDirectory/relativePath)
+// then plot the last Pareto front
+) -- (taskPlotLastParetoFront set ( directoryWithResults := workDirectory/relativePath, countInputs := 1) hook CopyFileHook(last_pareto, workDirectory/"last Pareto front SchafferN1.png" ) )

nsga2-test-functions/test function SchafferN2.oms

+import _file_.plottingFunctions._
+import _file_.execEnvironment._
+
+//model inputs
+val x = Val[Double]
+
+//model outputs
+val f1 = Val[Double]
+val f2 = Val[Double]
+
+// about the current experiment
+val relativePath = "results/SchafferN2"
+
+// the test function 
+val testFunctionSchafferN2=
+    ScalaTask("""
+        val f1 = if (x <= 1) {
+            -x*1.0
+        } else if (x <= 3) {
+            x - 2.0
+        } else if (x <= 4) {
+            4.0 - x
+        } else {
+            x - 4.0
+        }
+        val f2 = Math.pow(x - 5, 2);        
+        """) set (
+    inputs += x,
+    outputs += f1,
+    outputs += f2
+    )
+    
+// the optimisation algorithm under test 
+val evolutionSchafferN2 =
+  NSGA2Evolution(
+    genome = Seq(
+      x in (-5.0, 10.0)
+    ),
+    objective = Seq(f1, f2),
+    evaluation = testFunctionSchafferN2, 
+    parallelism = threads * 2,
+    termination = 1000
+  )
+  
+// compute evolution on the test Function
+(evolutionSchafferN2 on envMultiThread hook SavePopulationHook(evolutionSchafferN2, workDirectory/relativePath)
+// then plot the last Pareto front
+) -- (taskPlotLastParetoFront set ( directoryWithResults := workDirectory/relativePath, countInputs := 1) hook CopyFileHook(last_pareto, workDirectory/"last Pareto front SchafferN2.png" ) )