feat(loading): add nelder mead method launcher

alchemist-loading/build.gradle.kts

@@ -38,8 +38,11 @@ dependencies {
     implementation(libs.kasechange)
     implementation(libs.kotlin.reflect)
     implementation(libs.kotlin.coroutines.core)
+    implementation(libs.kotlinx.serialization.json)
+    implementation(libs.rrmxmx)
     implementation(libs.mongodb)
     implementation(libs.snakeyaml)
+    implementation(libs.caffeine)
 
     runtimeOnly(libs.groovy.jsr223)
     runtimeOnly(kotlin("scripting-jsr223"))
@@ -50,6 +53,7 @@ dependencies {
     testImplementation(libs.appdirs)
     testImplementation(libs.caffeine)
     testImplementation(libs.embedmongo)
+    testImplementation(libs.kotlin.test)
     testRuntimeOnly(incarnation("sapere"))
     testRuntimeOnly(incarnation("protelis"))
 }

alchemist-loading/src/main/kotlin/it/unibo/alchemist/boundary/NelderMeadMethod.kt

@@ -0,0 +1,181 @@
+package it.unibo.alchemist.boundary
+
+import com.github.benmanes.caffeine.cache.Caffeine
+import java.util.concurrent.ExecutorService
+import java.util.concurrent.Executors
+import java.util.concurrent.Future
+
+/**
+ * Nelder-Mead optimization method.
+ * Given an initial [simplex], this method iteratively refines the simplex to minimize a given [objective] function.
+ * The method is suitable for optimizing functions that are continuous but not differentiable.
+ * Other parameters are:
+ * - [alpha]: Reflection coefficient (standard value is 1.0);
+ * - [gamma]: Expansion coefficient (standard value is 2.0);
+ * - [rho]: Contraction coefficient (standard value is 0.5);
+ * - [sigma]: Shrink coefficient (standard value is 0.5);
+ * - [maxIterations]: Maximum number of iterations;
+ * - [tolerance]: Termination condition (small variation in function values).
+ */
+class NelderMeadMethod(
+    val simplex: List<Vertex>,
+    private val maxIterations: Int,
+    private val tolerance: Double,
+    private val alpha: Double, // Reflection coefficient
+    private val gamma: Double, // Expansion coefficient
+    private val rho: Double, // Contraction coefficient
+    private val sigma: Double, // Shrink coefficient
+    executorService: ExecutorService = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()),
+    private val objective: (List<Double>) -> Future<Double>,
+) {
+
+    private val cache = Caffeine.newBuilder().executor(executorService)
+        .build<List<Double>, Future<Double>> { coordinates ->
+            objective(coordinates)
+        }
+
+    /**
+     * Apply the Nelder-Mead optimization method to the given [simplex] and [objective] function.
+     */
+    fun optimize(): Vertex {
+        require(simplex.isNotEmpty()) { "The initial simplex must not be empty" }
+        val dimensions = simplex.first().size
+        require(dimensions > 0) { "The number of dimensions must be greater than 0" }
+        require(simplex.size == dimensions + 1) {
+            "The vertices of the initial simplex must be one more than the number of dimensions"
+        }
+        require(simplex.all { it.size == dimensions }) {
+            "All vertices of the initial simplex must have the same number of dimensions"
+        }
+        var symplexUpdated = simplex
+        repeat(maxIterations) {
+            // Sort simplex by function values
+            val sortedSimplex = symplexUpdated
+                .map { it to cache[it.valuesToList()] }
+                .sortedBy { it.second.get() }
+                .map { it.first }
+            val bestVertex: Vertex = sortedSimplex.first()
+            val worstVertex: Vertex = sortedSimplex.last()
+            val secondWorstVertex: Vertex = sortedSimplex[simplex.size - 2]
+            val bestValue = cache[bestVertex.valuesToList()]
+            val worstValues = worstVertex.valuesToList()
+            // Compute centroid (excluding worst point)
+            val centroid =
+                DoubleArray(dimensions) { index ->
+                    sortedSimplex.dropLast(1).sumOf { it[index] } / (sortedSimplex.size - 1)
+                }.toList()
+            // Reflections
+            val reflected: List<Double> = centroid.mapCentroid(alpha, worstValues)
+            val reflectedValue = cache[reflected]
+            require(reflectedValue.get().isFinite() && !reflectedValue.get().isNaN()) {
+                "Invalid objective function return value for reflection with $reflected = $reflectedValue.\n" +
+                    "Check the objective function implementation, the result should be a finite number."
+            }
+//            check (!reflectedValue.get().isFinite()) {
+//                error("Invalid objective function return value for reflection")
+//            }
+            val newSimplex = when {
+                reflectedValue < bestValue -> { // expansion
+                    val expanded: List<Double> = centroid.mapCentroid(gamma, reflected)
+                    when {
+                        cache[expanded] < reflectedValue -> sortedSimplex.updateLastVertex(expanded)
+                        else -> sortedSimplex.updateLastVertex(reflected)
+                    }
+                }
+                reflectedValue < cache[secondWorstVertex.valuesToList()] -> { // accept reflection
+                    sortedSimplex.updateLastVertex(reflected)
+                }
+                else -> { // contraction
+                    val contracted = when {
+                        reflectedValue < cache[worstValues] -> centroid.mapCentroid(rho, reflected)
+                        else -> centroid.mapCentroid(rho, worstValues)
+                    }
+                    when {
+                        cache[contracted] < cache[worstValues] -> sortedSimplex.updateLastVertex(contracted)
+                        else -> { // shrink simplex
+                            sortedSimplex.map { vertex ->
+                                Vertex(
+                                    vertex.keys().associateWith { key ->
+                                        // Find the index corresponding to the key
+                                        val index = bestVertex.keys().indexOf(key)
+                                        val oldValue = bestVertex[index]
+                                        oldValue + sigma * (vertex[index] - oldValue) // Apply shrink transformation
+                                    },
+                                )
+                            }
+                        }
+                    }
+                }
+            }
+            // Check termination condition (small variation in function values)
+            val functionValues = newSimplex.map { cache[it.valuesToList()].get() }
+            symplexUpdated = newSimplex
+            val maxDiff = functionValues.maxOrNull()!! - functionValues.minOrNull()!!
+            if (maxDiff < tolerance) return symplexUpdated.first()
+        }
+        return symplexUpdated.first()
+    }
+
+    private fun List<Double>.mapCentroid(coefficient: Double, values: List<Double>): List<Double> =
+        mapIndexed { index, value -> value + coefficient * (values[index] - value) }
+
+    private fun List<Vertex>.updateLastVertex(newVertex: List<Double>): List<Vertex> = mapIndexed { index, vertex ->
+        if (index == size - 1) {
+            Vertex(vertex.keys().associateWith { newVertex[vertex.keys().indexOf(it)] })
+        } else {
+            vertex
+        }
+    }
+
+    /**
+     * This companion object contains the comparison operators for [Future] of [Double] objects and [Double]s.
+     */
+    companion object {
+        /**
+         * Compares two [Future] of [Double] objects.
+         */
+        operator fun Future<Double>.compareTo(other: Future<Double>): Int = get().compareTo(other.get())
+
+        /**
+         * Compares a [Future] of [Double] object with a [Double].
+         */
+        operator fun Future<Double>.compareTo(other: Double): Int = get().compareTo(other)
+
+        /**
+         * Compares a [Double] with a [Future] of [Double] object.
+         */
+        operator fun Double.compareTo(other: Future<Double>): Int = compareTo(other.get())
+    }
+}
+
+/**
+ * A vertex of the simplex in the Nelder-Mead method.
+ */
+@JvmInline
+value class Vertex(private val vertex: Map<String, Double>) {
+    /**
+     * Returns amount of dimensions of the vertex.
+     */
+    val size: Int
+        get() = vertex.size
+
+    /**
+     * Returns the values of the vertex as a list.
+     */
+    fun valuesToList(): List<Double> = vertex.values.toList()
+
+    /**
+     * Returns the keys of the vertex as a set.
+     */
+    fun keys(): Set<String> = vertex.keys
+
+    /**
+     * Returns the value of the vertex at the given index.
+     */
+    operator fun get(index: Int) = valuesToList()[index]
+
+    /**
+     * Returns the key of the vertex at the given index.
+     */
+    fun keyAt(index: Int) = vertex.keys.elementAt(index)
+}

alchemist-loading/src/main/kotlin/it/unibo/alchemist/boundary/launchers/NelderMeadLauncher.kt

@@ -0,0 +1,166 @@
+package it.unibo.alchemist.boundary.launchers
+
+import it.unibo.alchemist.boundary.Launcher
+import it.unibo.alchemist.boundary.Loader
+import it.unibo.alchemist.boundary.NelderMeadMethod
+import it.unibo.alchemist.boundary.Variable
+import it.unibo.alchemist.boundary.Vertex
+import it.unibo.alchemist.model.Environment
+import kotlinx.serialization.ExperimentalSerializationApi
+import kotlinx.serialization.json.JsonPrimitive
+import kotlinx.serialization.json.buildJsonObject
+import org.danilopianini.rrmxmx.RrmxmxRandom
+import org.danilopianini.rrmxmx.RrmxmxRandom.Companion.DEFAULT_SEED
+import java.io.File
+import java.time.LocalDateTime
+import java.util.concurrent.Callable
+import java.util.concurrent.ConcurrentLinkedDeque
+import java.util.concurrent.ExecutorService
+import java.util.concurrent.Executors
+import java.util.concurrent.ForkJoinPool
+import java.util.concurrent.Future
+import java.util.concurrent.TimeUnit
+import java.util.concurrent.atomic.AtomicInteger
+import kotlin.Int.Companion.MAX_VALUE
+
+/**
+ * Alchemist launcher for the Nelder-Mead method optimization.
+ * This launcher optimize the [variables] of the simulation by a given [objectiveFunction].
+ * The optimization is done by the Nelder-Mead method, which is a derivative-free optimization algorithm.
+ * The optimization is done in parallel for a given number of [repetitions] of the [seedName] and [maxIterations].
+ * It will stop when the [tolerance] is reached.
+ * The [alpha], [gamma], [rho] and [sigma] parameters are the standard values for the Nelder-Mead method.
+ * The result of the optimization is written into `[outputPath]` folder in [outputFileName] as a JSON file.
+ */
+class NelderMeadLauncher
+@JvmOverloads
+constructor(
+    private val outputPath: String = "data",
+    private val outputFileName: String = "nelderMead",
+    private val objectiveFunction: Environment<*, *>.() -> Double,
+    private val variables: List<String> = emptyList(),
+    private val seedName: String,
+    private val repetitions: Int = 1,
+    private val maxIterations: Int = MAX_VALUE,
+    private val seed: ULong = DEFAULT_SEED,
+    private val tolerance: Double = 1e-6,
+    private val alpha: Double = 1.0, // standard value for the reflection in Nelder-Mead method
+    private val gamma: Double = 2.0, // standard value for the expansion in Nelder-Mead method
+    private val rho: Double = 0.5, // standard value for the contraction in Nelder-Mead method
+    private val sigma: Double = 0.5, // standard value for the shrink in Nelder-Mead method
+) : Launcher {
+    @OptIn(ExperimentalSerializationApi::class)
+    @Synchronized
+    override fun launch(loader: Loader) {
+        require(loader.variables.isNotEmpty() || variables.isNotEmpty()) {
+            "No variables found, can not optimize anything."
+        }
+        val simplexVertices: List<Map<String, Double>> = generateSymplexVertices(loader.variables)
+        val seeds: List<Int> =
+            loader.variables[seedName]
+                ?.stream()
+                ?.map {
+                    check(it is Number) { "Seed must be a number. $it is not." }
+                    it.toInt()
+                }?.toList()
+                ?.take(repetitions) ?: listOf(repetitions)
+        val executorID = AtomicInteger(0)
+        val executor =
+            Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()) {
+                Thread(it, "Alchemist Nelder Mead worker #${executorID.getAndIncrement()}")
+            }
+        val errorQueue = ConcurrentLinkedDeque<Throwable>()
+        loader.executeWithNelderMead(simplexVertices, executor) { vertex ->
+            val futureValues = seeds.map<Int, Future<Double>> { currentSeed ->
+                // associate keys to vertex values
+                val simulationParameters = variables
+                    .associateWith { vertex[variables.indexOf(it)] } + (seedName to currentSeed)
+                check(loader.variables.keys == simulationParameters.keys) {
+                    "Variables do not match: ${loader.variables.keys} != ${simulationParameters.keys}"
+                }
+                executor.submit<Double>(
+                    Callable {
+                        val simulation = loader.getWith<Any?, Nothing>(simulationParameters)
+                        simulation.play()
+                        simulation.run()
+                        if (simulation.error.isPresent) {
+                            errorQueue.add(simulation.error.get())
+                        }
+                        objectiveFunction(simulation.environment)
+                    },
+                )
+            }
+            ForkJoinPool.commonPool().submit(
+                Callable {
+                    futureValues.map { it.get() }.average()
+                },
+            )
+        }.also { result ->
+            // if not exists create the directory
+            File(outputPath).mkdirs()
+            val outputFile = File(
+                "$outputPath${File.separator}${outputFileName}" +
+                    "_${LocalDateTime.now().toString().replace(":", "-")}.json",
+            )
+            val outputContent = buildJsonObject {
+                put(seedName, JsonPrimitive(repetitions))
+                put("variables", buildJsonObject {
+                    variables.forEach { variable ->
+                        put(variable, JsonPrimitive(result[variable]))
+                    }
+                })
+            }.toString()
+            outputFile.writeText(outputContent)
+        }
+        executor.shutdown()
+        executor.awaitTermination(Long.MAX_VALUE, TimeUnit.HOURS)
+        if (errorQueue.isNotEmpty()) {
+            throw errorQueue.reduce { previous, other ->
+                previous.addSuppressed(other)
+                previous
+            }
+        }
+    }
+
+    private fun Loader.executeWithNelderMead(
+        simplexVertices: List<Map<String, Double>>,
+        executorService: ExecutorService,
+        executeFunction: (List<Double>) -> Future<Double>,
+    ): Map<String, Double> = NelderMeadMethod(
+        simplex = simplexVertices.map { Vertex(it) },
+        maxIterations = maxIterations,
+        tolerance = tolerance,
+        alpha = alpha,
+        gamma = gamma,
+        rho = rho,
+        sigma = sigma,
+        executorService = executorService,
+        objective = executeFunction,
+    ).optimize()
+        .let { result ->
+            this@NelderMeadLauncher.variables.associateWith {
+                result[this@NelderMeadLauncher.variables.indexOf(it)]
+            }
+        }
+
+    private fun generateSymplexVertices(loaderVariables: Map<String, Variable<*>>): List<Map<String, Double>> {
+        val randomGenerator = RrmxmxRandom(seed)
+        val instances: Map<String, ClosedRange<Double>> =
+            variables.associateWith { varName ->
+                val variable = loaderVariables.getValue(varName)
+                val allValues =
+                    variable
+                        .stream()
+                        .map {
+                            check(it is Number) {
+                                "All variables to optimize must be Numbers. $varName has value $it."
+                            }
+                            it.toDouble()
+                        }.toList()
+                allValues.min()..allValues.max()
+            }
+        return (0..variables.size).map {
+            instances.mapValues { (_, range) -> randomGenerator.nextDouble(range.start, range.endInclusive) }
+        }
+    }
+}