From b6175ee4ce0051a1eab3a1a44aea62bd093f7d75 Mon Sep 17 00:00:00 2001 From: chapuisk Date: Fri, 1 Nov 2024 09:36:03 +0100 Subject: [PATCH 1/2] [MORRIS] add matrix defined trajectories --- .../exploration/sampling/MorrisSampling.java | 255 +++++++++++++----- 1 file changed, 188 insertions(+), 67 deletions(-) diff --git a/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java b/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java index 0cffe0c8f4..0f771e970b 100644 --- a/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java +++ b/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java @@ -1,15 +1,20 @@ package gama.core.kernel.batch.exploration.sampling; import java.util.ArrayList; +import java.util.Arrays; import java.util.Collections; import java.util.LinkedHashMap; import java.util.List; import java.util.Map; import java.util.Random; +import java.util.stream.Collectors; import java.util.stream.IntStream; -import gama.core.kernel.experiment.ParametersSet; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; + import gama.core.kernel.experiment.IParameter.Batch; +import gama.core.kernel.experiment.ParametersSet; import gama.core.runtime.IScope; import gama.core.runtime.exceptions.GamaRuntimeException; @@ -40,6 +45,188 @@ public Trajectory(final List> points) { this.points = points; } } + + /* Retro-engine SAlib Morris sampling */ + + + /** + * Main method for Morris samples, give the samples with List and List> + * + * @param nb_levels + * the number of levels (Should be even, frequently 4) + * @param nb_sample + * the number of sample for each parameter. Add the end, the number of simulation is nb_sample * + * nb_parameters + * @param parameters + * Parameters of the model. + * @param scope + * @return samples for simulations. Size: nb_sample * inputs.size() + */ + public static List makeMorrisSampling(final int nb_levels, final int nb_sample, + final List parameters, final IScope scope) { + if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); + int nb_attributes = parameters.size(); + + + + List trajectories = + morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); + List nameInputs = new ArrayList<>(); + for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } + List> MorrisSamples = new ArrayList<>(); + trajectories.forEach(t -> { + t.points.forEach(p -> { + Map tmpMap = new LinkedHashMap<>(); + IntStream.range(0, nb_attributes).forEach(i -> { tmpMap.put(nameInputs.get(i), p.get(i)); }); + MorrisSamples.add(tmpMap); + }); + + }); + List result=new ArrayList<>(); + result.add(MorrisSamples); + result.add(buildParametersSetfromSample(scope, parameters, MorrisSamples)); + return result; + } + + /** + * Same as above but only give the sampling with a list of parameter set + * @param nb_levels + * the number of levels (Should be even, frequently 4) + * @param nb_sample + * the number of sample for each parameter. Add the end, the number of simulation is nb_sample * + * nb_parameters + * @param parameters + * Parameters of the model. + * @param scope + * @return + */ + public static List makeMorrisSamplingOnly(final int nb_levels, final int nb_sample, + final List parameters, final IScope scope) { + if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); + int nb_attributes = parameters.size(); + List trajectories = + morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); + List nameInputs = new ArrayList<>(); + for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } + List> MorrisSamples = new ArrayList<>(); + trajectories.forEach(t -> { + t.points.forEach(p -> { + Map tmpMap = new LinkedHashMap<>(); + IntStream.range(0, nb_attributes).forEach(i -> { tmpMap.put(nameInputs.get(i), p.get(i)); }); + MorrisSamples.add(tmpMap); + }); + }); + return buildParametersSetfromSample(scope, parameters, MorrisSamples); + } + + /* + * Create one random trajectory following classical random strategy: + * + * https://gsa-module.readthedocs.io/en/stable/implementation/morris_screening_method.html + * + * b∗ = (x∗ + Δ2 × ((2×b−j_k) × d∗ + j_k)) × p∗ + * + */ + private static Trajectory generateTraj(final int k, final int l, final Random rng) { + double delta = l / (2.0 * (l - 1)); + + RealMatrix g = MatrixUtils.createRealIdentityMatrix(k); + int numParams = g.getColumnDimension(); + int numGroups = g.getRowDimension(); + + RealMatrix B = MatrixUtils.createRealMatrix(numGroups+1, numGroups); + tril(B); + + RealMatrix P_star = MatrixUtils.createRealMatrix(numGroups, numGroups); + shuffleRow(P_star); + + RealMatrix J = ones(numGroups+1,numGroups); + + double[][] _Ds = new double[numParams][numParams]; + for (int i = 0; i < numParams; i++) { _Ds[i][i] = rng.nextBoolean() ? -1 : 1; } + RealMatrix D_star = MatrixUtils.createRealMatrix(_Ds); + + double[] x_star = seed(numParams, delta, l, rng); + + // b + RealMatrix b = g.multiply(P_star).transpose(); + // c + RealMatrix c = B.scalarMultiply(2d).multiply(b); + // d + RealMatrix d = c.subtract(J).multiply(D_star); + // EE + RealMatrix EE = d.add(J).scalarMultiply(delta/2); + + // Trajectory + RealMatrix B_star = rowise_add(x_star, EE); + + // Change to Gama Morris contract + List> trajectories = new ArrayList<>(); + for (int i = 0; i < B_star.getRowDimension(); i++) { + trajectories.add(Arrays.stream(B_star.getRow(i)).boxed().collect(Collectors.toList())); + } + + return new Trajectory(trajectories); + } + + // ----------------- INNER UTILITIES ----------------- // + + // Row wise addition to a 2D matrix : must be of same column dimension + private static RealMatrix rowise_add(double[] row, RealMatrix m) { + for (int i = 0; i < m.getRowDimension(); i++) { + final int fi = i; + m.setRow(i, + IntStream.range(0, m.getColumnDimension()) + .mapToDouble(j -> m.getEntry(fi, j)+row[j]).toArray()); + } + return m; + } + + // Original random point in the trajectory + private static double[] seed(int size, double delta, int levels, final Random rng) { + double[] seed = new double[size]; + double bound = 1 - delta; + double[] grid = new double[levels / 2]; + for (int i = 0; i < grid.length; i++) { + grid[i] = i * (bound / (grid.length - 1)); + } + + for (int i = 0; i < size; i++) { + seed[i] = grid[rng.nextInt(grid.length)]; + } + return seed; + } + + // 2D Matrix of ones + private static RealMatrix ones(int rows, int columns) { + double[][] ones = new double[rows][columns]; + for (int i = 0; i < rows; i++) { Arrays.fill(ones[i], 1.0); } + return MatrixUtils.createRealMatrix(ones); + } + + // Lower triangulation + private static RealMatrix tril(RealMatrix m) { + for (int i = 1; i < m.getRowDimension(); i++) { + for (int j = 0; j < i; j++) { + m.setEntry(i, j, 1); + } + } + return m; + } + + // Shuffles rows of the matrix + private static RealMatrix shuffleRow(RealMatrix m) { + List clm = IntStream.range(0, m.getColumnDimension()).boxed().collect(Collectors.toList()); + Collections.shuffle(clm); + for (int j = 0; j < m.getColumnDimension(); j++) { + m.setEntry(clm.remove(0), j, 1.0); + } + return m; + } + + ////////////////////////// + ////////////////////////// + /* OLD DIRTY TOM THINGS */ /** * For a given number of parameters k, a number of levels p, generation an initial seed for this parameters @@ -141,70 +328,4 @@ private static List morrisTrajectories(final int k, final int p, fin return addTrajectories(k, p, r - 1, rng, acc); } - /** - * Main method for Morris samples, give the samples with List and List> - * - * @param nb_levels - * the number of levels (Should be even, frequently 4) - * @param nb_sample - * the number of sample for each parameter. Add the end, the number of simulation is nb_sample * - * nb_parameters - * @param parameters - * Parameters of the model. - * @param scope - * @return samples for simulations. Size: nb_sample * inputs.size() - */ - public static List makeMorrisSampling(final int nb_levels, final int nb_sample, - final List parameters, final IScope scope) { - if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); - int nb_attributes = parameters.size(); - List trajectories = - morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); - List nameInputs = new ArrayList<>(); - for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } - List> MorrisSamples = new ArrayList<>(); - trajectories.forEach(t -> { - t.points.forEach(p -> { - Map tmpMap = new LinkedHashMap<>(); - IntStream.range(0, nb_attributes).forEach(i -> { tmpMap.put(nameInputs.get(i), p.get(i)); }); - MorrisSamples.add(tmpMap); - }); - - }); - List result=new ArrayList<>(); - result.add(MorrisSamples); - result.add(buildParametersSetfromSample(scope, parameters, MorrisSamples)); - return result; - } - - /** - * Same as above but only give the sampling with a list of parameter set - * @param nb_levels - * the number of levels (Should be even, frequently 4) - * @param nb_sample - * the number of sample for each parameter. Add the end, the number of simulation is nb_sample * - * nb_parameters - * @param parameters - * Parameters of the model. - * @param scope - * @return - */ - public static List makeMorrisSamplingOnly(final int nb_levels, final int nb_sample, - final List parameters, final IScope scope) { - if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); - int nb_attributes = parameters.size(); - List trajectories = - morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); - List nameInputs = new ArrayList<>(); - for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } - List> MorrisSamples = new ArrayList<>(); - trajectories.forEach(t -> { - t.points.forEach(p -> { - Map tmpMap = new LinkedHashMap<>(); - IntStream.range(0, nb_attributes).forEach(i -> { tmpMap.put(nameInputs.get(i), p.get(i)); }); - MorrisSamples.add(tmpMap); - }); - }); - return buildParametersSetfromSample(scope, parameters, MorrisSamples); - } } From 8e1e5ad92348b95db54f9abe795ac079cf5ba113 Mon Sep 17 00:00:00 2001 From: chapuisk Date: Fri, 1 Nov 2024 09:45:23 +0100 Subject: [PATCH 2/2] [MORRIS] fix #367 sampling issues out of bounds --- .../exploration/morris/MorrisExploration.java | 1 - .../exploration/sampling/MorrisSampling.java | 125 +++--------------- 2 files changed, 16 insertions(+), 110 deletions(-) diff --git a/gama.core/src/gama/core/kernel/batch/exploration/morris/MorrisExploration.java b/gama.core/src/gama/core/kernel/batch/exploration/morris/MorrisExploration.java index e67dce2bce..fd28e4775b 100644 --- a/gama.core/src/gama/core/kernel/batch/exploration/morris/MorrisExploration.java +++ b/gama.core/src/gama/core/kernel/batch/exploration/morris/MorrisExploration.java @@ -161,7 +161,6 @@ public void setChildren(final Iterable children) {} /** ######################### EVALUATE MORRIS INDEXES ######################### */ - @SuppressWarnings ("unchecked") @Override public void explore(final IScope scope) { this.sample = Cast.asInt(scope, getFacet(Exploration.SAMPLE_SIZE).value(scope)); diff --git a/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java b/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java index 0f771e970b..e1617f7095 100644 --- a/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java +++ b/gama.core/src/gama/core/kernel/batch/exploration/sampling/MorrisSampling.java @@ -67,10 +67,9 @@ public static List makeMorrisSampling(final int nb_levels, final int nb_ if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); int nb_attributes = parameters.size(); + List trajectories = new ArrayList<>(); + trajectories = morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator(), trajectories); - - List trajectories = - morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); List nameInputs = new ArrayList<>(); for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } List> MorrisSamples = new ArrayList<>(); @@ -104,8 +103,8 @@ public static List makeMorrisSamplingOnly(final int nb_levels, fi final List parameters, final IScope scope) { if (nb_levels % 2 != 0) throw GamaRuntimeException.error("The number of value should be even", scope); int nb_attributes = parameters.size(); - List trajectories = - morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator()); + List trajectories = new ArrayList<>(); + trajectories = morrisTrajectories(nb_attributes, nb_levels, nb_sample, scope.getRandom().getGenerator(), trajectories); List nameInputs = new ArrayList<>(); for (int i = 0; i < parameters.size(); i++) { nameInputs.add(parameters.get(i).getName()); } List> MorrisSamples = new ArrayList<>(); @@ -119,6 +118,18 @@ public static List makeMorrisSamplingOnly(final int nb_levels, fi return buildParametersSetfromSample(scope, parameters, MorrisSamples); } + /** + * Recursively generates r independent trajectories for k variables sampled with p levels. + * + * No specific strategy is used to sample trajectories - uniformely random + */ + private static List morrisTrajectories(final int k, final int p, final int r, final Random rng, final List acc) { + if (r == 0) return acc; + acc.add(generateTraj(k, p, rng)); + return morrisTrajectories(k, p, r - 1, rng, acc); + } + + /* * Create one random trajectory following classical random strategy: * @@ -223,109 +234,5 @@ private static RealMatrix shuffleRow(RealMatrix m) { } return m; } - - ////////////////////////// - ////////////////////////// - /* OLD DIRTY TOM THINGS */ - - /** - * For a given number of parameters k, a number of levels p, generation an initial seed for this parameters - */ - private static List seed(final int k, final int p, final Random rng) { - List seed = new ArrayList<>(); - double delta = 1 / (2 * ((double) p - 1)); - IntStream.range(0, k).forEach(i -> seed.add((rng.nextInt(p * 2 - 2) + 1) * delta)); - return seed; - } - - /** - * Build a trajectory (2nd function) - */ - private static List trajectoryBuilder(final double delta, final List order, final List direction, - final List seed, final List> accPoints, final List accdelta, final int index) { - if (order.isEmpty()) { - List trajectory = new ArrayList<>(); - trajectory.add(accPoints); - trajectory.add(accdelta); - return trajectory; - } - int idx = order.get(0); - double deltaOriented = delta * direction.get(0); - double valTemp = seed.get(idx) + deltaOriented; - List new_seed = new ArrayList<>(seed); - new_seed.set(idx, valTemp); - order.remove(0); - direction.remove(0); - accPoints.add(new_seed); - accdelta.add(deltaOriented); - return trajectoryBuilder(delta, order, direction, new_seed, accPoints, accdelta, index + 1); - } - - /** - * Build a trajectory (1st function) - */ - private static List trajectoryBuilder(final double delta, final List order, final List direction, - final List seed) { - List> accPoints = new ArrayList<>(); - List accDelta = new ArrayList<>(); - if (order.isEmpty()) { - // This is probably never used - List trajectory = new ArrayList<>(); - trajectory.add(accPoints); - trajectory.add(accDelta); - return trajectory; - } - int idx = order.get(0); - double deltaOriented = delta * direction.get(0); - double valTemp = seed.get(idx) + deltaOriented; - List new_seed = new ArrayList<>(seed); - new_seed.set(idx, valTemp); - order.remove(0); - direction.remove(0); - accPoints.add(new_seed); - accDelta.add(deltaOriented); - return trajectoryBuilder(delta, order, direction, new_seed, accPoints, accDelta, 1); - - } - - /** - * Create data for making trajectory k: Number of variable p: Number of levels (Should be even) return: new - * Trajectory composed of several points to visit - */ - @SuppressWarnings("unchecked") - private static Trajectory makeTrajectory(final int k, final int p, final Random rng) { - double delta = 1 / (2 * ((double) p - 1)); - List seed = seed(k, p, rng); - List orderVariables = new ArrayList<>(); - IntStream.range(0, k).forEach(orderVariables::add); - Collections.shuffle(orderVariables); - List directionVariables = new ArrayList<>(); - IntStream.range(0, k).forEach(s -> directionVariables.add(rng.nextInt(2) * 2 - 1)); - List List_p_d = trajectoryBuilder(delta, orderVariables, directionVariables, seed); - List> points = (List>) List_p_d.get(0); - return new Trajectory(points); - } - - /** - * Recursive function that add trajectories - */ - private static List addTrajectories(final int k, final int p, final int r, final Random rng, - final List acc) { - if (r == 0) return acc; - acc.add(makeTrajectory(k, p, rng)); - return addTrajectories(k, p, r - 1, rng, acc); - } - - /** - * Generates r independent trajectories for k variables sampled with p levels. - */ - private static List morrisTrajectories(final int k, final int p, final int r, final Random rng) { - List acc = new ArrayList<>(); - if (r == 0) - // Probably never used - return acc; - acc.add(makeTrajectory(k, p, rng)); - return addTrajectories(k, p, r - 1, rng, acc); - } }