From bdc07e6229364f0f6d222378f84b4561803ceca9 Mon Sep 17 00:00:00 2001 From: = <=> Date: Thu, 4 Jun 2026 01:09:19 +0200 Subject: [PATCH 01/20] Add builtin function lhs that performs latin hypercube sampling. The function only supports 2 dimensional matrices and return an NxN matrix of 0s and 1s. Also add a rudimentary test of the function, that checks if the number of samples is correct. --- scripts/builtin/lhs.dml | 41 +++++++++++ .../org/apache/sysds/common/Builtins.java | 3 +- .../builtin/part1/BuiltinLHSTest.java | 73 +++++++++++++++++++ src/test/scripts/functions/builtin/lhs.dml | 24 ++++++ 4 files changed, 140 insertions(+), 1 deletion(-) create mode 100644 scripts/builtin/lhs.dml create mode 100644 src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java create mode 100644 src/test/scripts/functions/builtin/lhs.dml diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml new file mode 100644 index 00000000000..ef206fb07eb --- /dev/null +++ b/scripts/builtin/lhs.dml @@ -0,0 +1,41 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# This script runs latin hypercube sampling. +# The Breslow method is used for handling ties and the regression parameters +# are computed using trust region newton method with conjugate gradient +# +# INPUT: +# ------------------------------------------------------------------------------------------------- +# N The number of samples to be taken +# +# ------------------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------------------------ +# M A N x N matrix M +# Contains N*(N-1) 0s and N 1s +# ------------------------------------------------------------------------------------------------------ + +lhs = function(Integer N) + return (Matrix[Double] M) { + M = table(seq(1,N),sample(N,N)) + } \ No newline at end of file diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java index e21c539d6d8..40d48d54f70 100644 --- a/src/main/java/org/apache/sysds/common/Builtins.java +++ b/src/main/java/org/apache/sysds/common/Builtins.java @@ -409,7 +409,8 @@ public enum Builtins { UPPER_TRI("upper.tri", false, true), EINSUM("einsum", false, false), XDUMMY1("xdummy1", true), //error handling test - XDUMMY2("xdummy2", true); //error handling test + XDUMMY2("xdummy2", true), //error handling test + LHS("lhs",true); Builtins(String name, boolean script) { this(name, null, script, false, ReturnType.SINGLE_RETURN); diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java new file mode 100644 index 00000000000..e6751c744cc --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -0,0 +1,73 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part1; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; + +public class BuiltinLHSTest extends AutomatedTestBase +{ + private final static String TEST_NAME = "lhs"; + private final static String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinLHSTest.class.getSimpleName() + "/"; + + + @Override + public void setUp() { + addTestConfiguration(TEST_NAME, new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"B"})); + } + + + @Test + public void testMisc() { + runLhsTest(5); + } + private void runLhsTest(int N) + { + ExecMode platformOld = setExecMode(ExecMode.HYBRID); + + try + { + String NString = Integer.toString(N); + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + String HOME = SCRIPT_DIR + TEST_DIR; + fullDMLScriptName = HOME + TEST_NAME + ".dml"; + programArgs = new String[]{"-args", NString, output("C")}; + + //execute test + runTest(true, false, null, -1); + + //compare number of values in sample and requested N + double val = readDMLMatrixFromOutputDir("C") + .values() + .stream() + .mapToDouble(x -> x) + .sum(); + Assert.assertEquals(N, val,0); + } + finally { + rtplatform = platformOld; + } + } +} \ No newline at end of file diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml new file mode 100644 index 00000000000..15c343b298c --- /dev/null +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -0,0 +1,24 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +N = $1 +M = lhs(N) +write(M, $2) From be8ad7fb4798759a98951a094b551c7af3aceee9 Mon Sep 17 00:00:00 2001 From: = <=> Date: Sun, 21 Jun 2026 22:51:44 +0200 Subject: [PATCH 02/20] Change lhs function to return a matrix of points coordinates instead of a NxN binary matrix The function now handles the multidimensional case. Also update the unit tests to check if the returned matrix is a valid latin hypercube. --- scripts/builtin/lhs.dml | 12 +++-- .../builtin/part1/BuiltinLHSTest.java | 51 ++++++++++++++----- src/test/scripts/functions/builtin/lhs.dml | 5 +- 3 files changed, 48 insertions(+), 20 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index ef206fb07eb..f8b62137d12 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -26,16 +26,20 @@ # INPUT: # ------------------------------------------------------------------------------------------------- # N The number of samples to be taken +# d The number of dimension # # ------------------------------------------------------------------------------------------------- # # OUTPUT: # ------------------------------------------------------------------------------------------------------ -# M A N x N matrix M -# Contains N*(N-1) 0s and N 1s +# M A N x d matrix M +# Contains the coordinates of sampled points # ------------------------------------------------------------------------------------------------------ -lhs = function(Integer N) +lhs = function(Integer N, Integer d) return (Matrix[Double] M) { - M = table(seq(1,N),sample(N,N)) + M = sample(N,N) + for (i in 1:d-1) { + M = cbind(M,sample(N,N)); + } } \ No newline at end of file diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index e6751c744cc..a27d676e0e5 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -19,12 +19,17 @@ package org.apache.sysds.test.functions.builtin.part1; -import org.junit.Assert; import org.junit.Test; +import static org.junit.Assert.assertFalse; +import static org.junit.Assert.assertTrue; + +import java.util.HashMap; + import org.apache.sysds.common.Types.ExecMode; import org.apache.sysds.test.AutomatedTestBase; import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; public class BuiltinLHSTest extends AutomatedTestBase { @@ -32,7 +37,24 @@ public class BuiltinLHSTest extends AutomatedTestBase private final static String TEST_DIR = "functions/builtin/"; private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinLHSTest.class.getSimpleName() + "/"; - + public void checkLatinHypercubeValidity(HashMap m,int N, int d ) { + for (int col = 1; col <= d; col++) { + boolean[] seen = new boolean[N + 1]; + for (int row = 1; row <= N; row++) { + + double val = m.get(new CellIndex(row, col)); + int intVal = (int) val; + assertTrue(intVal >= 1 && intVal <= N); + + assertFalse(seen[intVal]); + seen[intVal] = true; + } + for (int i = 1; i <= N; i++) { + assertTrue(seen[i]); + } + } + } + @Override public void setUp() { addTestConfiguration(TEST_NAME, new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"B"})); @@ -40,31 +62,32 @@ public void setUp() { @Test - public void testMisc() { - runLhsTest(5); + public void testTwoDim() { + runLhsTest(5,2); } - private void runLhsTest(int N) + + @Test + public void testMultiDim() { + runLhsTest(10,4); + } + + private void runLhsTest(int N,int d) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); try { - String NString = Integer.toString(N); loadTestConfiguration(getTestConfiguration(TEST_NAME)); String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; - programArgs = new String[]{"-args", NString, output("C")}; + programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), output("C")}; //execute test runTest(true, false, null, -1); - //compare number of values in sample and requested N - double val = readDMLMatrixFromOutputDir("C") - .values() - .stream() - .mapToDouble(x -> x) - .sum(); - Assert.assertEquals(N, val,0); + + HashMap m = readDMLMatrixFromOutputDir("C"); + checkLatinHypercubeValidity(m,N,d); } finally { rtplatform = platformOld; diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml index 15c343b298c..f78e299c981 100644 --- a/src/test/scripts/functions/builtin/lhs.dml +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -20,5 +20,6 @@ #------------------------------------------------------------- N = $1 -M = lhs(N) -write(M, $2) +d = $2 +M = lhs(N,d) +write(M, $3) From 29c561d6350b56e3ae899186b56db6a060898068 Mon Sep 17 00:00:00 2001 From: = <=> Date: Wed, 24 Jun 2026 20:29:44 +0200 Subject: [PATCH 03/20] add two new methods for lhs generation: improved and maximin the methods generate candidates and pick the best one based on distance to already chosen points TODO: there is a lot of code duplication between the two methods, candidate generation can maybe be vectorized --- scripts/builtin/lhs.dml | 139 +++++++++++++++++- .../builtin/part1/BuiltinLHSTest.java | 21 ++- src/test/scripts/functions/builtin/lhs.dml | 6 +- 3 files changed, 153 insertions(+), 13 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index f8b62137d12..4ad0c671599 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -36,10 +36,139 @@ # Contains the coordinates of sampled points # ------------------------------------------------------------------------------------------------------ -lhs = function(Integer N, Integer d) +minLhsDistance = function(Matrix[Double] M) + return (Double minD) { + D = dist(M) + diag(matrix(1/0, rows = nrow(M), cols = 1)) + minD = min(D) + } + +lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer random_repetitions = 1) return (Matrix[Double] M) { - M = sample(N,N) - for (i in 1:d-1) { - M = cbind(M,sample(N,N)); + + if (method == "random"){ + + current_best = 0 + M = matrix(0,rows = N, cols = d) + for (i in 1:random_repetitions) { + + newM = matrix(0,rows = N, cols = d) + for (i in 1:d) { + newM[,i] = sample(N,N); + } + + current = minLhsDistance(newM) + if(current > current_best){ + M = newM + current_best = current + } + } + } + + else if (method == "improved"){ + + # optimal distance + opt2 = (N / (N^(1/d)))^2 + + # available ints per column + available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) + + # generate first point at random + first_sample = sample(N,d,TRUE) + M = t(first_sample) + + # replace used ints for each column, with last sample + # in subsequent runs we only look at the first N-1,N-2,... rows + replacement_mask = table(first_sample,seq(1,N),N,d) + available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask + + # for the next N -2 points: + # generate reps candidates and choose the candidate, + # for which the minimal distance to the already chosen points + # is closest to the optimum distance + # and update the available_coords matrix + + for (i in 2:N-1){ + coords_left = N - i + 1 + current_best = 1/0 + row = matrix(0,rows = 1, cols = d) + replacement_mask = matrix(0,rows = coords_left, cols = d) + + for(j in 1:random_repetitions){ + # random mask for coordinates generation + rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) + + # get coordinates and flatten to row vector + vals_matrix = rand * available_coords + row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix + + # calculate minimal squared euclidean distance between candidate and already picked points + minD = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) + + # check if current candidate is better then the best + if(abs(opt2-minD) < current_best){ + row = row_candidate + current_best = abs(opt2-minD) + replacement_mask = rand + } + } + M = rbind(M,row) + # replace chosen coords with last, not picked coords + available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask + } + print(available_coords) + # for the last point there is only one choice left + M = rbind(M,available_coords[1,]) + } + + else if (method == "maximin"){ + + # available ints per column + available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) + + # generate first point at random + first_sample = sample(N,d,TRUE) + M = t(first_sample) + + # replace used ints for each column, with last sample + # in subsequent runs we only look at the first N-1,N-2,... rows + replacement_mask = table(first_sample,seq(1,N),N,d) + available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask + + # for the next N -2 points: + # generate reps candidates and choose the candidate, + # for which the minimal distance to the already chosen points is maximal + # and update the available_coords matrix + + for (i in 2:N-1){ + coords_left = N - i + 1 + current_best = 0 + row = matrix(0,rows = 1, cols = d) + replacement_mask = matrix(0,rows = coords_left, cols = d) + + for(j in 1:random_repetitions){ + # random mask for coordinates generation + rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) + + # get coordinates and flatten to row vector + vals_matrix = rand * available_coords + row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix + + # calculate minimal squared euclidean distance between candidate and already picked points + minD = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) + + # check if current candidate is better then the best + if(minD > current_best){ + row = row_candidate + current_best = minD + replacement_mask = rand + } + } + M = rbind(M,row) + # replace chosen coords with last, not picked coords + available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask + } + # for the last point there is only one choice left + M = rbind(M,available_coords[1,]) } - } \ No newline at end of file + + } diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index a27d676e0e5..42e9a9a153b 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -62,16 +62,25 @@ public void setUp() { @Test - public void testTwoDim() { - runLhsTest(5,2); + public void testRandomTwoDim() { + runLhsTest(5,2,5,"random"); } @Test - public void testMultiDim() { - runLhsTest(10,4); + public void testRandomMultiDim() { + runLhsTest(10,4,5,"random"); } - private void runLhsTest(int N,int d) + @Test + public void testImproved() { + runLhsTest(10,4,5,"improved"); + } + @Test + public void testMaximin() { + runLhsTest(10,4,5,"maximin"); + } + + private void runLhsTest(int N,int d, int repetitions, String method) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); @@ -80,7 +89,7 @@ private void runLhsTest(int N,int d) loadTestConfiguration(getTestConfiguration(TEST_NAME)); String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; - programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), output("C")}; + programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method ,output("C")}; //execute test runTest(true, false, null, -1); diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml index f78e299c981..839f41feaf4 100644 --- a/src/test/scripts/functions/builtin/lhs.dml +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -21,5 +21,7 @@ N = $1 d = $2 -M = lhs(N,d) -write(M, $3) +reps = $3 +method = $4 +M = lhs(N = N,d = d,random_repetitions = reps, method = method) +write(M, $5) From 3819a6d721c093bbb20f52fbfeb352a54bebcb4c Mon Sep 17 00:00:00 2001 From: = <=> Date: Fri, 3 Jul 2026 20:52:27 +0200 Subject: [PATCH 04/20] consolidate improved and maximin method and add cp_sweep method the improved and maximin method differ only in how they evaluate candidates, therefore they are consolidated into one method "build" additional goal parameter is added, which determines how the candidates are determined also add new method "cp_sweep", which instead of creating potential candidates for new points, swaps values in columns to create the biggest improvement --- scripts/builtin/lhs.dml | 248 ++++++++++++------ .../builtin/part1/BuiltinLHSTest.java | 21 +- src/test/scripts/functions/builtin/lhs.dml | 5 +- 3 files changed, 182 insertions(+), 92 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 4ad0c671599..7dc4bdd999a 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -20,14 +20,25 @@ #------------------------------------------------------------- # This script runs latin hypercube sampling. -# The Breslow method is used for handling ties and the regression parameters -# are computed using trust region newton method with conjugate gradient +# # # INPUT: # ------------------------------------------------------------------------------------------------- # N The number of samples to be taken # d The number of dimension -# +# method The method to be used for sampling, possible values are: +# - "random" - random sampling +# - "build" - build the sample iteratively +# - "cp_sweep" - iterative columnwise pairwise swapping +# reps The number of repetitions for each method: +# - for "random" method, the best sample is chosen from reps random samples +# - for "build" method, reps candidates are generated for each point and the best candidate is chosen +# - for "cp_sweep" method, reps sweeps are performed for each column +# goal The goal of the sampling, possible values are: +# - "maximin" - maximize the minimum distance between points +# - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) +# - "sum_inv" - minimize the sum of the inverse distances between points +# eps The threshold for stopping the cp_sweep method # ------------------------------------------------------------------------------------------------- # # OUTPUT: @@ -37,23 +48,37 @@ # ------------------------------------------------------------------------------------------------------ minLhsDistance = function(Matrix[Double] M) - return (Double minD) { + return (Double min_dist) { D = dist(M) + diag(matrix(1/0, rows = nrow(M), cols = 1)) - minD = min(D) + min_dist = min(D) + } + +sumInvDistance = function(Matrix[Double] M) + return (Double sum_inv) { + inv_dist = 1/lower.tri(target = dist(M),diag=FALSE, values = TRUE) + # replace infinite values with 0 and sum the inverse distances + inv_dist = replace(target = inv_dist,pattern = 1/0,replacement = 0) + sum_inv = sum(inv_dist) } -lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer random_repetitions = 1) +swap = function(Matrix[Double] M, Integer col, Integer row1, Integer row2) return (Matrix[Double] M) { + temp = M[row1,col] + M[row1,col] = M[row2,col] + M[row2,col] = temp + } +lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer reps = 1, String goal = "maximin", Double eps = 0.001) + return (Matrix[Double] M) { if (method == "random"){ current_best = 0 M = matrix(0,rows = N, cols = d) - for (i in 1:random_repetitions) { + for (i in 1:reps) { newM = matrix(0,rows = N, cols = d) - for (i in 1:d) { - newM[,i] = sample(N,N); + for (j in 1:d) { + newM[,j] = sample(N,N); } current = minLhsDistance(newM) @@ -64,8 +89,8 @@ lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer } } - else if (method == "improved"){ - + + else if (method == "build"){ # optimal distance opt2 = (N / (N^(1/d)))^2 @@ -83,92 +108,147 @@ lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer # for the next N -2 points: # generate reps candidates and choose the candidate, - # for which the minimal distance to the already chosen points - # is closest to the optimum distance + # for which the minimal distance to the already chosen points is maximal # and update the available_coords matrix - for (i in 2:N-1){ - coords_left = N - i + 1 + for (i in 2:(N-1)){ + #number of coordinates left to choose from + coords_left = N - i + 1 + # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv + if(goal == "maximin"){ + current_best = 0 + } + else{ current_best = 1/0 - row = matrix(0,rows = 1, cols = d) - replacement_mask = matrix(0,rows = coords_left, cols = d) - - for(j in 1:random_repetitions){ - # random mask for coordinates generation - rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) - - # get coordinates and flatten to row vector - vals_matrix = rand * available_coords - row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix - - # calculate minimal squared euclidean distance between candidate and already picked points - minD = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) - - # check if current candidate is better then the best - if(abs(opt2-minD) < current_best){ - row = row_candidate - current_best = abs(opt2-minD) - replacement_mask = rand + } + + #initialize row vector and replacement mask for available coordinates + row = matrix(0,rows = 1, cols = d) + replacement_mask = matrix(0,rows = coords_left, cols = d) + + for(j in 1:reps){ + # random mask for coordinates generation + rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) + + # get coordinates and flatten to row vector + vals_matrix = rand * available_coords + row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix + + # check if the candidate is better then the current best + # accoring to the goal and update the row vector and replacement mask accordingly + if(goal == "maximin"){ + min_dist = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) + if(min_dist > current_best){ + row = row_candidate + current_best = min_dist + replacement_mask = rand } - } - M = rbind(M,row) - # replace chosen coords with last, not picked coords - available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask + } + else if(goal == "opt"){ + dist_from_opt = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2))) + if(abs(dist_from_opt - opt2) < current_best){ + row = row_candidate + current_best = dist_from_opt + replacement_mask = rand + } + } + else if(goal == "sum_inv"){ + #calculate new matrix and inverse distances + sum_inv = sumInvDistance(rbind(M,row_candidate)) + + if(sum_inv < current_best){ + row = row_candidate + current_best = sum_inv + replacement_mask = rand + } + } + } + # add the best candidate to the matrix of chosen points + M = rbind(M,row) + # replace chosen coordinates with last, not yet 1picked coordinates + available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask } - print(available_coords) # for the last point there is only one choice left M = rbind(M,available_coords[1,]) } - - else if (method == "maximin"){ - - # available ints per column - available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) - - # generate first point at random - first_sample = sample(N,d,TRUE) - M = t(first_sample) + # TODO: vectorize the inner loops, currently this method is very slow + else if (method == "cp_sweep"){ + # optimal distance + opt2 = (N / (N^(1/d)))^2 - # replace used ints for each column, with last sample - # in subsequent runs we only look at the first N-1,N-2,... rows - replacement_mask = table(first_sample,seq(1,N),N,d) - available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask + # generate initial random latin hypercube + M = matrix(0,rows = N, cols = d) + for (i in 1:d) { + M[,i] = sample(N,N); + } - # for the next N -2 points: - # generate reps candidates and choose the candidate, - # for which the minimal distance to the already chosen points is maximal - # and update the available_coords matrix + # initialize G according to the goal + if(goal == "maximin"){ + G = minLhsDistance(M) + } + else if(goal == "opt"){ + G = abs(opt2 - minLhsDistance(M)) + } + else if(goal == "sum_inv"){ + G = sumInvDistance(M) + } - for (i in 2:N-1){ - coords_left = N - i + 1 - current_best = 0 - row = matrix(0,rows = 1, cols = d) - replacement_mask = matrix(0,rows = coords_left, cols = d) - - for(j in 1:random_repetitions){ - # random mask for coordinates generation - rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) - - # get coordinates and flatten to row vector - vals_matrix = rand * available_coords - row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix - - # calculate minimal squared euclidean distance between candidate and already picked points - minD = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) - - # check if current candidate is better then the best - if(minD > current_best){ - row = row_candidate - current_best = minD - replacement_mask = rand + pred = TRUE + iter = 0 + while (pred){ + iter = iter + 1; + if(iter >= reps){ + pred = FALSE; + } + #iterate over all columns and swap all pairs of rows in each column + for(j in 1:d){ + # initialize improvement vector, for easier calculation we flatten the 2D matrix to a 1D vector + # where the index of the improvement vector corresponds to the row pair (k,l) with k < l + # so the 1D vector index = l + (k-1)*N + improvement_vector = matrix(0,rows = 1, cols = N*(N-1)) + for(k in 1:(N-1)){ + for(l in (k+1):(N)){ + # add improvement to the improvement vector according to the goal + if(goal == "maximin"){ + newM = swap(M,j,k,l) + improvement_vector[1, l + (k-1)*N] = minLhsDistance(newM) - G + } + else if(goal == "opt"){ + newM = swap(M,j,k,l) + improvement_vector[1, l + (k-1)*N] = G - (abs(opt2 - minLhsDistance(newM))) + } + else if(goal == "sum_inv"){ + newM = swap(M,j,k,l) + improvement_vector[1, l + (k-1)*N] = G - sumInvDistance(newM) } + } } - M = rbind(M,row) - # replace chosen coords with last, not picked coords - available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask + #if there was an improvement, swap the rows with the best improvement + if(max(improvement_vector) > 0){ + #calculate k and l from the index of the improvement vector + idx = as.scalar(rowIndexMax(improvement_vector)) + row1 = ((idx - 1) %/% N) + 1 + row2 = ((idx - 1) %% N) + 1 + M = swap(M,j,row1,row2) + } + } + if(goal == "maximin"){ + newG = minLhsDistance(M) + } + else if(goal == "opt"){ + newG = abs(opt2 - minLhsDistance(M)) + } + else if(goal == "sum_inv"){ + newG = sumInvDistance(M) + } + # stop if no significant improvement was made, otherwise update G + # we can use abs() here, because we only update M if there was an improvement + if(abs(G - newG) < eps){ + pred = FALSE + } + else{ + G = newG + } } - # for the last point there is only one choice left - M = rbind(M,available_coords[1,]) } - } diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index 42e9a9a153b..53472c1cc87 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -63,24 +63,33 @@ public void setUp() { @Test public void testRandomTwoDim() { - runLhsTest(5,2,5,"random"); + runLhsTest(5,2,5,"random","opt"); } @Test public void testRandomMultiDim() { - runLhsTest(10,4,5,"random"); + runLhsTest(10,4,5,"random","opt"); } @Test public void testImproved() { - runLhsTest(10,4,5,"improved"); + runLhsTest(10,4,5,"build","opt"); + } + @Test + public void testSumInv() { + runLhsTest(10,4,5,"build","sum_inv"); } @Test public void testMaximin() { - runLhsTest(10,4,5,"maximin"); + runLhsTest(10,4,5,"build","maximin"); + } + + @Test + public void testCPsweep() { + runLhsTest(10,4,5,"cp_sweep","sum_inv"); } - private void runLhsTest(int N,int d, int repetitions, String method) + private void runLhsTest(int N,int d, int repetitions, String method, String goal) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); @@ -89,7 +98,7 @@ private void runLhsTest(int N,int d, int repetitions, String method) loadTestConfiguration(getTestConfiguration(TEST_NAME)); String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; - programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method ,output("C")}; + programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method, goal, output("C")}; //execute test runTest(true, false, null, -1); diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml index 839f41feaf4..dae502ab0c2 100644 --- a/src/test/scripts/functions/builtin/lhs.dml +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -23,5 +23,6 @@ N = $1 d = $2 reps = $3 method = $4 -M = lhs(N = N,d = d,random_repetitions = reps, method = method) -write(M, $5) +goal = $5 +M = lhs(N = N,d = d,reps = reps, method = method, goal = goal) +write(M, $6) From d96b2c647f38bf4a2f44776e53364d5efa3c2fc0 Mon Sep 17 00:00:00 2001 From: Mikolaj <83309779+grzybkos@users.noreply.github.com> Date: Sun, 5 Jul 2026 04:02:51 +0200 Subject: [PATCH 05/20] add new method for lhs generation: genetic --- scripts/builtin/lhs.dml | 134 +++++++++++++++++- .../builtin/part1/BuiltinLHSTest.java | 5 +- 2 files changed, 134 insertions(+), 5 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 4ad0c671599..45854b5729c 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -25,9 +25,14 @@ # # INPUT: # ------------------------------------------------------------------------------------------------- -# N The number of samples to be taken -# d The number of dimension -# +# N Number of samples to be taken +# d Number of dimensions +# method Method to be used for sampling ("random", "improved", "maximin" or "genetic") +# random_repetitions Number of random repetitions +# genetic_generations Number of generations ("genetic" method) +# genetic_population_size Size of the population ("genetic" method) +# genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) +# # ------------------------------------------------------------------------------------------------- # # OUTPUT: @@ -42,7 +47,41 @@ minLhsDistance = function(Matrix[Double] M) minD = min(D) } -lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer random_repetitions = 1) +sOptimality = function(Matrix[Double] M) + return (Double score) { + D = dist(M) + score = sum(D) / (nrow(M) * (nrow(M) - 1)) + } + +designStart = function(Integer index, Integer N) + return (Integer start) { + start = (index - 1) * N + 1 + } +designEnd = function(Integer index, Integer N) + return (Integer end) { + end = index * N + } +getDesignByIndex = function(Matrix[Double] population, Integer index, Integer N) + return (Matrix[Double] design) { + start = designStart(index, N) + end = designEnd(index, N) + design = population[start:end,] + } +attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) + return (Matrix[Double] design) { + if (as.scalar(rand(rows=1, cols=1, min=0, max=1)) < mutation_rate) { + column_to_mutate = as.scalar(sample(ncol(design), 1)) + rows_to_mutate = sample(nrow(design), 2) + row1 = as.scalar(rows_to_mutate[1,1]) + row2 = as.scalar(rows_to_mutate[2,1]) + + tmp = design[row1, column_to_mutate] + design[row1, column_to_mutate] = design[row2, column_to_mutate] + design[row2, column_to_mutate] = tmp + } + } + +lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer random_repetitions = 1, Integer genetic_generations = 5, Integer genetic_population_size = 10, Double genetic_mutation_rate = 0.25) return (Matrix[Double] M) { if (method == "random"){ @@ -171,4 +210,91 @@ lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer M = rbind(M,available_coords[1,]) } + # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html + else if (method == "genetic") { + population = matrix(0, rows=genetic_population_size*N, cols=d) + scores = matrix(0, rows=genetic_population_size, cols=1) + + # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix + for (p in 1:genetic_population_size) { + design = lhs(N = N, d = d, method = "random", random_repetitions = random_repetitions) + population_start = designStart(p, N) + population_end = designEnd(p, N) + population[population_start:population_end,] = design + + # 2. Calculate the S-optimality score for each design + scores[p,1] = sOptimality(design) + } + + for (generation in 1:genetic_generations) { + # 3. Put the best design in the first position of survivors and throw away the worst half of the remaining population + survivors_cutoff = floor(genetic_population_size / 2) + 1 + survivors = matrix(0, rows=survivors_cutoff * N, cols=d) + for (survivor in 1:survivors_cutoff) { + current_best_index = as.scalar(rowIndexMax(t(scores))) + current_best_design = getDesignByIndex(population, current_best_index, N) + + # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again + survivors_start = designStart(survivor, N) + survivors_end = designEnd(survivor, N) + survivors[survivors_start:survivors_end,] = current_best_design + scores[current_best_index, 1] = -1/0 + } + + # 4. Create new designs from the survivors + new_population = matrix(0, rows=genetic_population_size * N, cols=d) + best_design = survivors[1:N,] + new_population[1:N,] = best_design + next_free_slot = 2 + + # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. + for (survivor_index in 2:survivors_cutoff) { + new_design = getDesignByIndex(survivors, survivor_index, N) + source_col = as.scalar(sample(d, 1)) + target_col = as.scalar(sample(d, 1)) + new_design[,target_col] = best_design[,source_col] + + # 5. Attempt to mutate the child design by switching two elements in a random column. + new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + + # Store the generated child design in the new population + new_population_start = designStart(next_free_slot, N) + new_population_end = designEnd(next_free_slot, N) + new_population[new_population_start:new_population_end,] = new_design + next_free_slot = next_free_slot + 1 + } + + # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. + for (survivor_index in 2:survivors_cutoff) { + if (next_free_slot <= genetic_population_size) { + source_design = getDesignByIndex(survivors, survivor_index, N) + new_design = best_design + source_col = as.scalar(sample(d, 1)) + target_col = as.scalar(sample(d, 1)) + new_design[,target_col] = source_design[,source_col] + + # 5. Attempt to mutate the child design by switching two elements in a random column. + new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + + # Store the generated child design in the new population + new_population_start = designStart(next_free_slot, N) + new_population_end = designEnd(next_free_slot, N) + new_population[new_population_start:new_population_end,] = new_design + next_free_slot = next_free_slot + 1 + } + } + + population = new_population + + # Score the new population + for (p in 1:genetic_population_size) { + design = getDesignByIndex(population, p, N) + scores[p,1] = sOptimality(design) + } + } + + # Return the best design from the final population + best_index = as.scalar(rowIndexMax(t(scores))) + M = getDesignByIndex(population, best_index, N) + } } diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index 42e9a9a153b..ef4ec544042 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -79,7 +79,10 @@ public void testImproved() { public void testMaximin() { runLhsTest(10,4,5,"maximin"); } - + @Test + public void testGenetic() { + runLhsTest(10,4,5,"genetic"); + } private void runLhsTest(int N,int d, int repetitions, String method) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); From abdbe644be9e5cadfed4ac406d1eda22a9f6a7cc Mon Sep 17 00:00:00 2001 From: Mikolaj <83309779+grzybkos@users.noreply.github.com> Date: Sun, 5 Jul 2026 04:05:10 +0200 Subject: [PATCH 06/20] remove leftover comment --- scripts/builtin/lhs.dml | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 45854b5729c..640b4285aea 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -154,7 +154,6 @@ lhs = function(Integer N = 10, Integer d = 2, String method = "random", Integer # replace chosen coords with last, not picked coords available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask } - print(available_coords) # for the last point there is only one choice left M = rbind(M,available_coords[1,]) } From 17b19e6d590f04e36b09c5f9e5923392eb40a23d Mon Sep 17 00:00:00 2001 From: = <=> Date: Sun, 5 Jul 2026 16:42:50 +0200 Subject: [PATCH 07/20] add functions for goal evaluation for full matrix and for a new row to simplify objective calculation between methods change maximin and sOptimality to minimize the negatie value objective to match other metods add return_type parameter, which controls whether matrix with integer coordianates or double uniform distribution 0-1 is returned --- scripts/builtin/lhs.dml | 660 ++++++++++++++++++++-------------------- 1 file changed, 330 insertions(+), 330 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 5bbaa36c99c..3690a8d7567 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -1,364 +1,340 @@ -#------------------------------------------------------------- -# -# Licensed to the Apache Software Foundation (ASF) under one -# or more contributor license agreements. See the NOTICE file -# distributed with this work for additional information -# regarding copyright ownership. The ASF licenses this file -# to you under the Apache License, Version 2.0 (the -# "License"); you may not use this file except in compliance -# with the License. You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, -# software distributed under the License is distributed on an -# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY -# KIND, either express or implied. See the License for the -# specific language governing permissions and limitations -# under the License. -# -#------------------------------------------------------------- - -# This script runs latin hypercube sampling. -# -# -# INPUT: -# ------------------------------------------------------------------------------------------------- -# N The number of samples to be taken -# d The number of dimension -# method The method to be used for sampling, possible values are: -# - "random" - random sampling -# - "build" - build the sample iteratively -# - "cp_sweep" - iterative columnwise pairwise swapping -# - "genetic" - genetic algorithm -# reps The number of repetitions for each method: -# - for "random" method, the best sample is chosen from reps random samples -# - for "build" method, reps candidates are generated for each point and the best candidate is chosen -# - for "cp_sweep" method, reps sweeps are performed for each column -# goal The goal of the sampling, possible values are: -# - "maximin" - maximize the minimum distance between points -# - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) -# - "sum_inv" - minimize the sum of the inverse distances between points -# eps The threshold for stopping the cp_sweep method -# genetic_generations Number of generations ("genetic" method) -# genetic_population_size Size of the population ("genetic" method) -# genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) -# ------------------------------------------------------------------------------------------------- -# -# OUTPUT: -# ------------------------------------------------------------------------------------------------------ -# M A N x d matrix M -# Contains the coordinates of sampled points -# ------------------------------------------------------------------------------------------------------ - -minLhsDistance = function(Matrix[Double] M) - return (Double min_dist) { - D = dist(M) + diag(matrix(1/0, rows = nrow(M), cols = 1)) - min_dist = min(D) - } -swap = function(Matrix[Double] M, Integer col, Integer row1, Integer row2) - return (Matrix[Double] M) { - temp = M[row1,col] - M[row1,col] = M[row2,col] - M[row2,col] = temp - } + #------------------------------------------------------------- + # + # Licensed to the Apache Software Foundation (ASF) under one + # or more contributor license agreements. See the NOTICE file + # distributed with this work for additional information + # regarding copyright ownership. The ASF licenses this file + # to you under the Apache License, Version 2.0 (the + # "License"); you may not use this file except in compliance + # with the License. You may obtain a copy of the License at + # + # http://www.apache.org/licenses/LICENSE-2.0 + # + # Unless required by applicable law or agreed to in writing, + # software distributed under the License is distributed on an + # "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + # KIND, either express or implied. See the License for the + # specific language governing permissions and limitations + # under the License. + # + #------------------------------------------------------------- + + # This script runs latin hypercube sampling. + # + # + # INPUT: + # ------------------------------------------------------------------------------------------------- + # N The number of samples to be taken + # d The number of dimension + # method The method to be used for sampling, possible values are: + # - "random" - random sampling + # - "build" - build the sample iteratively + # - "cp_sweep" - iterative columnwise pairwise swapping + # - "genetic" - genetic algorithm + # reps The number of repetitions for each method: + # - for "random" method, the best sample is chosen from reps random samples + # - for "build" method, reps candidates are generated for each point and the best candidate is chosen + # - for "cp_sweep" method, reps sweeps are performed for each column + # - for "genetic" method, reps candidates are generate in the initial design + # goal The goal to be minimized during the sampling procedure, possible values are: + # - "maximin" - maximize the minimum distance between points (in pracitce the negative minimum distance is minimized) + # - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) + # - "sum_inv" - minimize the sum of the inverse distances between points + # - "avg_dist" - maximize avg euclidean distance between points (in practice the negative is minimized) + # eps The threshold for stopping the cp_sweep method + # genetic_generations Number of generations ("genetic" method) + # genetic_population_size Size of the population ("genetic" method) + # genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) + # return_type How should the matrix be returned: + # - "coordinates": return the matrix of vectors with point integer point coordinates in the d-dimensional space + # - "distribution": return the matrix as vectors sampled from the d-dimensional unifom distribution in range 0-1 + # ------------------------------------------------------------------------------------------------- + # + # OUTPUT: + # ------------------------------------------------------------------------------------------------------ + # M A N x d matrix M + # Contains the coordinates of sampled points + # ------------------------------------------------------------------------------------------------------ + + minLhsDistance = function(Matrix[Double] M) + return (Double min_dist) { + D = dist(M) + diag(matrix(1/0, rows = nrow(M), cols = 1)) + min_dist = min(D) + } + swap = function(Matrix[Double] M, Integer col, Integer row1, Integer row2) + return (Matrix[Double] M) { + temp = M[row1,col] + M[row1,col] = M[row2,col] + M[row2,col] = temp + } -sOptimality = function(Matrix[Double] M) - return (Double score) { - D = dist(M) - score = sum(D) / (nrow(M) * (nrow(M) - 1)) - } + sOptimality = function(Matrix[Double] M) + return (Double score) { + D = dist(M) + score = sum(D) / (nrow(M) * (nrow(M) - 1)) + } -designStart = function(Integer index, Integer N) - return (Integer start) { - start = (index - 1) * N + 1 - } -designEnd = function(Integer index, Integer N) - return (Integer end) { - end = index * N + designStart = function(Integer index, Integer N) + return (Integer start) { + start = (index - 1) * N + 1 + } + designEnd = function(Integer index, Integer N) + return (Integer end) { + end = index * N + } + getDesignByIndex = function(Matrix[Double] population, Integer index, Integer N) + return (Matrix[Double] design) { + start = designStart(index, N) + end = designEnd(index, N) + design = population[start:end,] + } + attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) + return (Matrix[Double] design) { + if (as.scalar(rand(rows=1, cols=1, min=0, max=1)) < mutation_rate) { + column_to_mutate = as.scalar(sample(ncol(design), 1)) + rows_to_mutate = sample(nrow(design), 2) + row1 = as.scalar(rows_to_mutate[1,1]) + row2 = as.scalar(rows_to_mutate[2,1]) + design = swap(design,column_to_mutate,row1,row2) + } + } + + sumInvDistance = function(Matrix[Double] M) + return (Double sum_inv) { + inv_dist = 1/lower.tri(target = dist(M),diag=FALSE, values = TRUE) + # replace infinite values with 0 and sum the inverse distances + inv_dist = replace(target = inv_dist,pattern = 1/0,replacement = 0) + sum_inv = sum(inv_dist) } -getDesignByIndex = function(Matrix[Double] population, Integer index, Integer N) - return (Matrix[Double] design) { - start = designStart(index, N) - end = designEnd(index, N) - design = population[start:end,] + + evaluateGoalFullMatrix = function (Matrix[Double] M, String goal) + return(Double G){ + if(goal == "maximin"){ + G = -minLhsDistance(M) + } + else if(goal == "opt"){ + N = nrow(M) + d = ncol(M) + opt = N/(N^(1/d)) + G = abs(opt - minLhsDistance(M)) + } + else if(goal == "sum_inv"){ + G = sumInvDistance(M) + } + else if(goal == "avg_dist"){ + G = -sOptimality(M) + } } -attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) - return (Matrix[Double] design) { - if (as.scalar(rand(rows=1, cols=1, min=0, max=1)) < mutation_rate) { - column_to_mutate = as.scalar(sample(ncol(design), 1)) - rows_to_mutate = sample(nrow(design), 2) - row1 = as.scalar(rows_to_mutate[1,1]) - row2 = as.scalar(rows_to_mutate[2,1]) - design = swap(design,column_to_mutate,row1,row2) - } + + evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal) + return(Double G){ + if(goal == "maximin"){ + G = -min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)) + } + else if(goal == "opt"){ + N = nrow(M) + d = ncol(M) + opt2 = (N/(N^(1/d)))^2 + G = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2))) + } + else if(goal == "sum_inv"){ + G = sumInvDistance(rbind(M,row)) + } + else if(goal == "avg_dist"){ + G = -sOptimality(rbind(M,row)) + } } -sumInvDistance = function(Matrix[Double] M) -return (Double sum_inv) { - inv_dist = 1/lower.tri(target = dist(M),diag=FALSE, values = TRUE) - # replace infinite values with 0 and sum the inverse distances - inv_dist = replace(target = inv_dist,pattern = 1/0,replacement = 0) - sum_inv = sum(inv_dist) -} - - - -lhs = function( - Integer N = 10, - Integer d = 2, - String method = "random", - Integer reps = 1, - String goal = "maximin", - Double eps = 0.001, - Integer genetic_generations = 5, - Integer genetic_population_size = 10, - Double genetic_mutation_rate = 0.25 - ) - return (Matrix[Double] M) { - if (method == "random"){ - - current_best = 0 - M = matrix(0,rows = N, cols = d) - for (i in 1:reps) { - - newM = matrix(0,rows = N, cols = d) - for (j in 1:d) { - newM[,j] = sample(N,N); - } - current = minLhsDistance(newM) - if(current > current_best){ - M = newM - current_best = current + + lhs = function( + Integer N = 10, + Integer d = 2, + String method = "random", + Integer reps = 1, + String goal = "maximin", + Double eps = 0.001, + Integer genetic_generations = 5, + Integer genetic_population_size = 10, + Double genetic_mutation_rate = 0.25, + String return_type = "coordinates" + ) + return (Matrix[Double] M) { + if (method == "random"){ + + G = 1/0 + M = matrix(0,rows = N, cols = d) + for (i in 1:reps) { + + newM = matrix(0,rows = N, cols = d) + for (j in 1:d) { + newM[,j] = sample(N,N); + } + + newG = evaluateGoalFullMatrix(newM,goal) + if(newG < G){ + M = newM + G = newG + } } } - } - else if (method == "build"){ - # optimal distance - opt2 = (N / (N^(1/d)))^2 + else if (method == "build"){ + # optimal distance + opt2 = (N / (N^(1/d)))^2 - # available ints per column - available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) + # available ints per column + available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) - # generate first point at random - first_sample = sample(N,d,TRUE) - M = t(first_sample) + # generate first point at random + first_sample = sample(N,d,TRUE) + M = t(first_sample) - # replace used ints for each column, with last sample - # in subsequent runs we only look at the first N-1,N-2,... rows - replacement_mask = table(first_sample,seq(1,N),N,d) - available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask + # replace used ints for each column, with last sample + # in subsequent runs we only look at the first N-1,N-2,... rows + replacement_mask = table(first_sample,seq(1,N),N,d) + available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask - # for the next N -2 points: - # generate reps candidates and choose the candidate, - # for which the minimal distance to the already chosen points is maximal - # and update the available_coords matrix + # for the next N -2 points: + # generate reps candidates and choose the candidate, + # for which the minimal distance to the already chosen points is maximal + # and update the available_coords matrix - for (i in 2:(N-1)){ - #number of coordinates left to choose from - coords_left = N - i + 1 - # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv - if(goal == "maximin"){ - current_best = 0 - } - else{ - current_best = 1/0 - } + for (i in 2:(N-1)){ + #number of coordinates left to choose from + coords_left = N - i + 1 + # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv + G = 1/0 - #initialize row vector and replacement mask for available coordinates - row = matrix(0,rows = 1, cols = d) - replacement_mask = matrix(0,rows = coords_left, cols = d) + #initialize row vector and replacement mask for available coordinates + row = matrix(0,rows = 1, cols = d) + replacement_mask = matrix(0,rows = coords_left, cols = d) - for(j in 1:reps){ - # random mask for coordinates generation - rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) + for(j in 1:reps){ + # random mask for coordinates generation + rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) - # get coordinates and flatten to row vector - vals_matrix = rand * available_coords - row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix + # get coordinates and flatten to row vector + vals_matrix = rand * available_coords + row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix - # check if the candidate is better then the current best - # accoring to the goal and update the row vector and replacement mask accordingly - if(goal == "maximin"){ - min_dist = min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2)) - if(min_dist > current_best){ - row = row_candidate - current_best = min_dist - replacement_mask = rand - } - } - else if(goal == "opt"){ - dist_from_opt = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row_candidate))^2))) - if(abs(dist_from_opt - opt2) < current_best){ - row = row_candidate - current_best = dist_from_opt - replacement_mask = rand - } - } - else if(goal == "sum_inv"){ - #calculate new matrix and inverse distances - sum_inv = sumInvDistance(rbind(M,row_candidate)) + # check if the candidate is better then the current best + # accoring to the goal and update the row vector and replacement mask accordingly - if(sum_inv < current_best){ + newG = evaluateGoalNewRow(M,row_candidate,goal) + if(newG < G){ row = row_candidate - current_best = sum_inv + G = newG replacement_mask = rand } - } + } + # add the best candidate to the matrix of chosen points + M = rbind(M,row) + # replace chosen coordinates with last, not yet 1picked coordinates + available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask } - # add the best candidate to the matrix of chosen points - M = rbind(M,row) - # replace chosen coordinates with last, not yet 1picked coordinates - available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask - } - # for the last point there is only one choice left - M = rbind(M,available_coords[1,]) - } - # TODO: vectorize the inner loops, currently this method is very slow - else if (method == "cp_sweep"){ - # optimal distance - opt2 = (N / (N^(1/d)))^2 - - # generate initial random latin hypercube - M = matrix(0,rows = N, cols = d) - for (i in 1:d) { - M[,i] = sample(N,N); + # for the last point there is only one choice left + M = rbind(M,available_coords[1,]) } + # TODO: vectorize the inner loops, currently this method is very slow + else if (method == "cp_sweep"){ + # optimal distance + opt2 = (N / (N^(1/d)))^2 + + # generate initial random latin hypercube + M = matrix(0,rows = N, cols = d) + for (i in 1:d) { + M[,i] = sample(N,N); + } - # initialize G according to the goal - if(goal == "maximin"){ - G = minLhsDistance(M) - } - else if(goal == "opt"){ - G = abs(opt2 - minLhsDistance(M)) - } - else if(goal == "sum_inv"){ - G = sumInvDistance(M) - } + # initialize G according to the goal + G = evaluateGoalFullMatrix(M,goal) - pred = TRUE - iter = 0 - while (pred){ - iter = iter + 1; - if(iter >= reps){ - pred = FALSE; - } - #iterate over all columns and swap all pairs of rows in each column - for(j in 1:d){ - # initialize improvement vector, for easier calculation we flatten the 2D matrix to a 1D vector - # where the index of the improvement vector corresponds to the row pair (k,l) with k < l - # so the 1D vector index = l + (k-1)*N - improvement_vector = matrix(0,rows = 1, cols = N*(N-1)) - for(k in 1:(N-1)){ - for(l in (k+1):(N)){ - # add improvement to the improvement vector according to the goal - if(goal == "maximin"){ - newM = swap(M,j,k,l) - improvement_vector[1, l + (k-1)*N] = minLhsDistance(newM) - G - } - else if(goal == "opt"){ - newM = swap(M,j,k,l) - improvement_vector[1, l + (k-1)*N] = G - (abs(opt2 - minLhsDistance(newM))) - } - else if(goal == "sum_inv"){ + pred = TRUE + iter = 0 + while (pred){ + iter = iter + 1; + if(iter >= reps){ + pred = FALSE; + } + #iterate over all columns and swap all pairs of rows in each column + for(j in 1:d){ + # initialize new goal vector, for easier calculation we flatten the 2D matrix to a 1D vector + # where the index of the improvement vector corresponds to the row pair (k,l) with k < l + # so the 1D vector index = l + (k-1)*N + new_goal_vector = matrix(0,rows = 1, cols = N*(N-1)) + for(k in 1:(N-1)){ + for(l in (k+1):(N)){ + # add new goal to the goal vector according to the goal newM = swap(M,j,k,l) - improvement_vector[1, l + (k-1)*N] = G - sumInvDistance(newM) + new_goal_vector[1, l + (k-1)*N] = G - evaluateGoalFullMatrix(newM,goal) } } + #if there was an improvement, swap the rows with the best improvement + if(min(new_goal_vector) < G){ + #calculate k and l from the index of the goal vector + idx = as.scalar(rowIndexMax(new_goal_vector)) + row1 = ((idx - 1) %/% N) + 1 + row2 = ((idx - 1) %% N) + 1 + M = swap(M,j,row1,row2) + + } } - #if there was an improvement, swap the rows with the best improvement - if(max(improvement_vector) > 0){ - #calculate k and l from the index of the improvement vector - idx = as.scalar(rowIndexMax(improvement_vector)) - row1 = ((idx - 1) %/% N) + 1 - row2 = ((idx - 1) %% N) + 1 - M = swap(M,j,row1,row2) + newG = evaluateGoalFullMatrix(M,goal) + # stop if no significant improvement was made, otherwise update G + # we can use abs() here, because we only update M if there was an improvement + if(abs(G - newG) < eps){ + pred = FALSE + } + else{ + G = newG } } - if(goal == "maximin"){ - newG = minLhsDistance(M) - } - else if(goal == "opt"){ - newG = abs(opt2 - minLhsDistance(M)) - } - else if(goal == "sum_inv"){ - newG = sumInvDistance(M) - } - # stop if no significant improvement was made, otherwise update G - # we can use abs() here, because we only update M if there was an improvement - if(abs(G - newG) < eps){ - pred = FALSE - } - else{ - G = newG - } - } - } - # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html - else if (method == "genetic") { - population = matrix(0, rows=genetic_population_size*N, cols=d) - scores = matrix(0, rows=genetic_population_size, cols=1) - - # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix - for (p in 1:genetic_population_size) { - design = lhs(N = N, d = d, method = "random", reps = reps) - population_start = designStart(p, N) - population_end = designEnd(p, N) - population[population_start:population_end,] = design - - # 2. Calculate the S-optimality score for each design - scores[p,1] = sOptimality(design) } + # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html + else if (method == "genetic") { + population = matrix(0, rows=genetic_population_size*N, cols=d) + scores = matrix(0, rows=genetic_population_size, cols=1) - for (generation in 1:genetic_generations) { - # 3. Put the best design in the first position of survivors and throw away the worst half of the remaining population - survivors_cutoff = floor(genetic_population_size / 2) + 1 - survivors = matrix(0, rows=survivors_cutoff * N, cols=d) - for (survivor in 1:survivors_cutoff) { - current_best_index = as.scalar(rowIndexMax(t(scores))) - current_best_design = getDesignByIndex(population, current_best_index, N) - - # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again - survivors_start = designStart(survivor, N) - survivors_end = designEnd(survivor, N) - survivors[survivors_start:survivors_end,] = current_best_design - scores[current_best_index, 1] = -1/0 - } + # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix + for (p in 1:genetic_population_size) { + design = lhs(N = N, d = d, method = "random", reps = reps) + population_start = designStart(p, N) + population_end = designEnd(p, N) + population[population_start:population_end,] = design - # 4. Create new designs from the survivors - new_population = matrix(0, rows=genetic_population_size * N, cols=d) - best_design = survivors[1:N,] - new_population[1:N,] = best_design - next_free_slot = 2 - - # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. - for (survivor_index in 2:survivors_cutoff) { - new_design = getDesignByIndex(survivors, survivor_index, N) - source_col = as.scalar(sample(d, 1)) - target_col = as.scalar(sample(d, 1)) - new_design[,target_col] = best_design[,source_col] - - # 5. Attempt to mutate the child design by switching two elements in a random column. - new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) - - # Store the generated child design in the new population - new_population_start = designStart(next_free_slot, N) - new_population_end = designEnd(next_free_slot, N) - new_population[new_population_start:new_population_end,] = new_design - next_free_slot = next_free_slot + 1 + # 2. Calculate the S-optimality score for each design + scores[p,1] = sOptimality(design) } - # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. - for (survivor_index in 2:survivors_cutoff) { - if (next_free_slot <= genetic_population_size) { - source_design = getDesignByIndex(survivors, survivor_index, N) - new_design = best_design + for (generation in 1:genetic_generations) { + # 3. Put the best design in the first position of survivors and throw away the worst half of the remaining population + survivors_cutoff = floor(genetic_population_size / 2) + 1 + survivors = matrix(0, rows=survivors_cutoff * N, cols=d) + for (survivor in 1:survivors_cutoff) { + current_best_index = as.scalar(rowIndexMax(t(scores))) + current_best_design = getDesignByIndex(population, current_best_index, N) + + # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again + survivors_start = designStart(survivor, N) + survivors_end = designEnd(survivor, N) + survivors[survivors_start:survivors_end,] = current_best_design + scores[current_best_index, 1] = -1/0 + } + + # 4. Create new designs from the survivors + new_population = matrix(0, rows=genetic_population_size * N, cols=d) + best_design = survivors[1:N,] + new_population[1:N,] = best_design + next_free_slot = 2 + + # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. + for (survivor_index in 2:survivors_cutoff) { + new_design = getDesignByIndex(survivors, survivor_index, N) source_col = as.scalar(sample(d, 1)) target_col = as.scalar(sample(d, 1)) - new_design[,target_col] = source_design[,source_col] + new_design[,target_col] = best_design[,source_col] # 5. Attempt to mutate the child design by switching two elements in a random column. new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) @@ -369,19 +345,43 @@ lhs = function( new_population[new_population_start:new_population_end,] = new_design next_free_slot = next_free_slot + 1 } - } - population = new_population + # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. + for (survivor_index in 2:survivors_cutoff) { + if (next_free_slot <= genetic_population_size) { + source_design = getDesignByIndex(survivors, survivor_index, N) + new_design = best_design + source_col = as.scalar(sample(d, 1)) + target_col = as.scalar(sample(d, 1)) + new_design[,target_col] = source_design[,source_col] + + # 5. Attempt to mutate the child design by switching two elements in a random column. + new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + + # Store the generated child design in the new population + new_population_start = designStart(next_free_slot, N) + new_population_end = designEnd(next_free_slot, N) + new_population[new_population_start:new_population_end,] = new_design + next_free_slot = next_free_slot + 1 + } + } - # Score the new population - for (p in 1:genetic_population_size) { - design = getDesignByIndex(population, p, N) - scores[p,1] = sOptimality(design) + population = new_population + + # Score the new population + for (p in 1:genetic_population_size) { + design = getDesignByIndex(population, p, N) + scores[p,1] = sOptimality(design) + } } + + # Return the best design from the final population + best_index = as.scalar(rowIndexMax(t(scores))) + M = getDesignByIndex(population, best_index, N) } - # Return the best design from the final population - best_index = as.scalar(rowIndexMax(t(scores))) - M = getDesignByIndex(population, best_index, N) + if(return_type == "distribution"){ + r = rand(rows = N, cols = d, min = 0, max = 1/N) + M = M/N - r + } } - } From 61582cc10bc7ba7e9b4235113ebce097f0d72b3a Mon Sep 17 00:00:00 2001 From: = <=> Date: Wed, 8 Jul 2026 17:45:41 +0200 Subject: [PATCH 08/20] add checks for correct parameters and break up lhs function into seperate functions by method --- scripts/builtin/lhs.dml | 755 +++++++++++++++++++++------------------- 1 file changed, 392 insertions(+), 363 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 3690a8d7567..5eed04670ca 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -1,387 +1,416 @@ - #------------------------------------------------------------- - # - # Licensed to the Apache Software Foundation (ASF) under one - # or more contributor license agreements. See the NOTICE file - # distributed with this work for additional information - # regarding copyright ownership. The ASF licenses this file - # to you under the Apache License, Version 2.0 (the - # "License"); you may not use this file except in compliance - # with the License. You may obtain a copy of the License at - # - # http://www.apache.org/licenses/LICENSE-2.0 - # - # Unless required by applicable law or agreed to in writing, - # software distributed under the License is distributed on an - # "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY - # KIND, either express or implied. See the License for the - # specific language governing permissions and limitations - # under the License. - # - #------------------------------------------------------------- - - # This script runs latin hypercube sampling. - # - # - # INPUT: - # ------------------------------------------------------------------------------------------------- - # N The number of samples to be taken - # d The number of dimension - # method The method to be used for sampling, possible values are: - # - "random" - random sampling - # - "build" - build the sample iteratively - # - "cp_sweep" - iterative columnwise pairwise swapping - # - "genetic" - genetic algorithm - # reps The number of repetitions for each method: - # - for "random" method, the best sample is chosen from reps random samples - # - for "build" method, reps candidates are generated for each point and the best candidate is chosen - # - for "cp_sweep" method, reps sweeps are performed for each column - # - for "genetic" method, reps candidates are generate in the initial design - # goal The goal to be minimized during the sampling procedure, possible values are: - # - "maximin" - maximize the minimum distance between points (in pracitce the negative minimum distance is minimized) - # - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) - # - "sum_inv" - minimize the sum of the inverse distances between points - # - "avg_dist" - maximize avg euclidean distance between points (in practice the negative is minimized) - # eps The threshold for stopping the cp_sweep method - # genetic_generations Number of generations ("genetic" method) - # genetic_population_size Size of the population ("genetic" method) - # genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) - # return_type How should the matrix be returned: - # - "coordinates": return the matrix of vectors with point integer point coordinates in the d-dimensional space - # - "distribution": return the matrix as vectors sampled from the d-dimensional unifom distribution in range 0-1 - # ------------------------------------------------------------------------------------------------- - # - # OUTPUT: - # ------------------------------------------------------------------------------------------------------ - # M A N x d matrix M - # Contains the coordinates of sampled points - # ------------------------------------------------------------------------------------------------------ - - minLhsDistance = function(Matrix[Double] M) - return (Double min_dist) { - D = dist(M) + diag(matrix(1/0, rows = nrow(M), cols = 1)) - min_dist = min(D) - } - swap = function(Matrix[Double] M, Integer col, Integer row1, Integer row2) - return (Matrix[Double] M) { - temp = M[row1,col] - M[row1,col] = M[row2,col] - M[row2,col] = temp - } +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# This script runs latin hypercube sampling. +# +# INPUT: +# ------------------------------------------------------------------------------------------------- +# N The number of samples to be taken +# d The number of dimension +# method The method to be used for sampling, possible values are: +# - "random" - random sampling +# - "build" - build the sample iteratively +# - "cp_sweep" - iterative columnwise pairwise swapping +# - "genetic" - genetic algorithm +# reps The number of repetitions for each method: +# - for "random" method, the best sample is chosen from reps random samples +# - for "build" method, reps candidates are generated for each point and the best candidate is chosen +# - for "cp_sweep" method, reps sweeps are performed for each column +# - for "genetic" method, reps candidates are generate in the initial design +# goal The goal to be minimized during the sampling procedure, possible values are: +# - "maximin" - maximize the minimum distance between points (in practice the negative minimum distance is minimized) +# - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) +# - "sum_inv" - minimize the sum of the inverse distances between points +# - "avg_dist" - maximize avg euclidean distance between points (in practice the negative is minimized) +# eps The threshold for stopping the cp_sweep method +# genetic_generations Number of generations ("genetic" method) +# genetic_population_size Size of the population ("genetic" method) +# genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) +# return_type How should the matrix be returned: +# - "coordinates": return the matrix of vectors with point integer point coordinates in the d-dimensional space +# - "distribution": return the matrix as vectors sampled from the d-dimensional uniform distribution in range 0-1 +# ------------------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------------------------ +# M A N x d matrix M +# Contains the coordinates of sampled points +# ------------------------------------------------------------------------------------------------------ + +lhs = function( + Integer N = 10, + Integer d = 2, + String method = "random", + Integer reps = 1, + String goal = "maximin", + Double eps = 0.001, + Integer genetic_generations = 5, + Integer genetic_population_size = 10, + Double genetic_mutation_rate = 0.25, + String return_type = "buckets" +) +return (Matrix[Double] M) { + if (goal != "maximin" & goal != "opt" & goal != "sum_inv" & goal != "avg_dist") { + stop("Goal should be one of: maximin, opt, sum_inv, avg_dist") + } + if (return_type != "buckets" & return_type != "numerical") { + stop("Return type should be one of: buckets, numerical") + } - sOptimality = function(Matrix[Double] M) - return (Double score) { - D = dist(M) - score = sum(D) / (nrow(M) * (nrow(M) - 1)) - } + if (method == "random") { + M = randomLHS(N, d, reps, goal) + } + else if (method == "build") { + M = buildLHS(N, d, reps, goal) + } + else if (method == "cp_sweep") { + M = cpSweepLHS(N, d, reps, goal, eps) + } + # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html + else if (method == "genetic") { + M = geneticLHS(N, d, reps, goal, genetic_generations, genetic_population_size, genetic_mutation_rate) + } + else { + stop("Method should be one of: random, build, cp_sweep, genetic") + } - designStart = function(Integer index, Integer N) - return (Integer start) { - start = (index - 1) * N + 1 - } - designEnd = function(Integer index, Integer N) - return (Integer end) { - end = index * N + if (return_type == "numerical") { + r = rand(rows = N, cols = d, min = 0, max = 1 / N) + M = M / N - r + } +} + +randomLHS = function(Integer N, Integer d, Integer reps, String goal) +return (Matrix[Double] M) { + G = 1 / 0 + M = matrix(0, rows = N, cols = d) + for (i in 1:reps) { + newM = matrix(0, rows = N, cols = d) + for (j in 1:d) { + newM[, j] = sample(N, N) } - getDesignByIndex = function(Matrix[Double] population, Integer index, Integer N) - return (Matrix[Double] design) { - start = designStart(index, N) - end = designEnd(index, N) - design = population[start:end,] + + newG = evaluateGoalFullMatrix(newM, goal) + if (newG < G) { + M = newM + G = newG } - attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) - return (Matrix[Double] design) { - if (as.scalar(rand(rows=1, cols=1, min=0, max=1)) < mutation_rate) { - column_to_mutate = as.scalar(sample(ncol(design), 1)) - rows_to_mutate = sample(nrow(design), 2) - row1 = as.scalar(rows_to_mutate[1,1]) - row2 = as.scalar(rows_to_mutate[2,1]) - design = swap(design,column_to_mutate,row1,row2) + } +} + +buildLHS = function(Integer N, Integer d, Integer reps, String goal) +return (Matrix[Double] M) { + # optimal distance + opt2 = (N / (N^(1 / d)))^2 + + # available buckets per column + free_buckets = seq(1, N) %*% matrix(1, rows = 1, cols = d) + + # generate first point at random + first_sample = sample(N, d, TRUE) + M = t(first_sample) + + # replace used ints for each column, with last sample + # in subsequent runs we only look at the first N-1,N-2,... rows + replacement_mask = table(first_sample, seq(1, N), N, d) + free_buckets = (1 - replacement_mask) * free_buckets + (matrix(1, rows = N, cols = 1) %*% (free_buckets[N, ])) * replacement_mask + + # for the next N - 2 points: + # generate reps candidates and choose the candidate, + # for which the minimal distance to the already chosen points is maximal + # and update the free_buckets matrix + for (i in 2:(N - 1)) { + # number of buckets left to choose from + num_free_buckets = N - i + 1 + # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv + G = 1 / 0 + + # initialize row vector and replacement mask for free buckets + row = matrix(0, rows = 1, cols = d) + replacement_mask = matrix(0, rows = num_free_buckets, cols = d) + + for (j in 1:reps) { + # random mask for bucket generation + rand = table(sample(num_free_buckets, d, TRUE), seq(1, N), N, d) + + # get bucket and flatten to row vector + vals_matrix = rand * free_buckets + row_candidate = matrix(1, rows = 1, cols = N) %*% vals_matrix + + # check if the candidate is better then the current best + # according to the goal and update the row vector and replacement mask accordingly + newG = evaluateGoalNewRow(M, row_candidate, goal) + if (newG < G) { + row = row_candidate + G = newG + replacement_mask = rand } } - sumInvDistance = function(Matrix[Double] M) - return (Double sum_inv) { - inv_dist = 1/lower.tri(target = dist(M),diag=FALSE, values = TRUE) - # replace infinite values with 0 and sum the inverse distances - inv_dist = replace(target = inv_dist,pattern = 1/0,replacement = 0) - sum_inv = sum(inv_dist) + # add the best candidate to the matrix of chosen points + M = rbind(M, row) + # replace chosen buckets with last, not yet picked buckets + free_buckets = (1 - replacement_mask) * free_buckets + (matrix(1, rows = N, cols = 1) %*% (free_buckets[num_free_buckets, ])) * replacement_mask } - evaluateGoalFullMatrix = function (Matrix[Double] M, String goal) - return(Double G){ - if(goal == "maximin"){ - G = -minLhsDistance(M) - } - else if(goal == "opt"){ - N = nrow(M) - d = ncol(M) - opt = N/(N^(1/d)) - G = abs(opt - minLhsDistance(M)) - } - else if(goal == "sum_inv"){ - G = sumInvDistance(M) - } - else if(goal == "avg_dist"){ - G = -sOptimality(M) - } - } + # for the last point there is only one bucket left + M = rbind(M, free_buckets[1, ]) +} - evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal) - return(Double G){ - if(goal == "maximin"){ - G = -min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)) - } - else if(goal == "opt"){ - N = nrow(M) - d = ncol(M) - opt2 = (N/(N^(1/d)))^2 - G = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2))) - } - else if(goal == "sum_inv"){ - G = sumInvDistance(rbind(M,row)) - } - else if(goal == "avg_dist"){ - G = -sOptimality(rbind(M,row)) - } +# TODO: vectorize the inner loops, currently this method is very slow +cpSweepLHS = function(Integer N, Integer d, Integer reps, String goal, Double eps) +return (Matrix[Double] M) { + # optimal distance + opt2 = (N / (N^(1 / d)))^2 + + # generate initial random latin hypercube + M = matrix(0, rows = N, cols = d) + for (i in 1:d) { + M[, i] = sample(N, N) } + # initialize G according to the goal + G = evaluateGoalFullMatrix(M, goal) + pred = TRUE + iter = 0 + while (pred) { + iter = iter + 1 + if (iter >= reps) { + pred = FALSE + } - lhs = function( - Integer N = 10, - Integer d = 2, - String method = "random", - Integer reps = 1, - String goal = "maximin", - Double eps = 0.001, - Integer genetic_generations = 5, - Integer genetic_population_size = 10, - Double genetic_mutation_rate = 0.25, - String return_type = "coordinates" - ) - return (Matrix[Double] M) { - if (method == "random"){ - - G = 1/0 - M = matrix(0,rows = N, cols = d) - for (i in 1:reps) { - - newM = matrix(0,rows = N, cols = d) - for (j in 1:d) { - newM[,j] = sample(N,N); - } - - newG = evaluateGoalFullMatrix(newM,goal) - if(newG < G){ - M = newM - G = newG - } + # iterate over all columns and swap all pairs of rows in each column + for (j in 1:d) { + # initialize new goal vector, for easier calculation we flatten the 2D matrix to a 1D vector + # where the index of the improvement vector corresponds to the row pair (k,l) with k < l + # so the 1D vector index = l + (k-1)*N + new_goal_vector = matrix(0, rows = 1, cols = N * (N - 1)) + for (k in 1:(N - 1)) { + for (l in (k + 1):(N)) { + # add new goal to the goal vector according to the goal + newM = swap(M, j, k, l) + new_goal_vector[1, l + (k - 1) * N] = G - evaluateGoalFullMatrix(newM, goal) } } - - else if (method == "build"){ - # optimal distance - opt2 = (N / (N^(1/d)))^2 - - # available ints per column - available_coords = seq(1,N) %*% matrix(1,rows = 1, cols = d) - - # generate first point at random - first_sample = sample(N,d,TRUE) - M = t(first_sample) - - # replace used ints for each column, with last sample - # in subsequent runs we only look at the first N-1,N-2,... rows - replacement_mask = table(first_sample,seq(1,N),N,d) - available_coords = (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[N,])) * replacement_mask - - # for the next N -2 points: - # generate reps candidates and choose the candidate, - # for which the minimal distance to the already chosen points is maximal - # and update the available_coords matrix - - for (i in 2:(N-1)){ - #number of coordinates left to choose from - coords_left = N - i + 1 - # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv - G = 1/0 - - #initialize row vector and replacement mask for available coordinates - row = matrix(0,rows = 1, cols = d) - replacement_mask = matrix(0,rows = coords_left, cols = d) - - for(j in 1:reps){ - # random mask for coordinates generation - rand = table(sample(coords_left,d,TRUE),seq(1,N),N,d) - - # get coordinates and flatten to row vector - vals_matrix = rand * available_coords - row_candidate = matrix(1,rows = 1, cols = N) %*% vals_matrix - - # check if the candidate is better then the current best - # accoring to the goal and update the row vector and replacement mask accordingly - - newG = evaluateGoalNewRow(M,row_candidate,goal) - if(newG < G){ - row = row_candidate - G = newG - replacement_mask = rand - } - } - # add the best candidate to the matrix of chosen points - M = rbind(M,row) - # replace chosen coordinates with last, not yet 1picked coordinates - available_coords= (1 - replacement_mask) * available_coords + (matrix(1,rows = N, cols = 1) %*% (available_coords[coords_left,])) * replacement_mask - } - # for the last point there is only one choice left - M = rbind(M,available_coords[1,]) + # if there was an improvement, swap the rows with the best improvement + if (min(new_goal_vector) < G) { + # calculate k and l from the index of the goal vector + idx = as.scalar(rowIndexMax(new_goal_vector)) + row1 = ((idx - 1) %/% N) + 1 + row2 = ((idx - 1) %% N) + 1 + M = swap(M, j, row1, row2) } - # TODO: vectorize the inner loops, currently this method is very slow - else if (method == "cp_sweep"){ - # optimal distance - opt2 = (N / (N^(1/d)))^2 - - # generate initial random latin hypercube - M = matrix(0,rows = N, cols = d) - for (i in 1:d) { - M[,i] = sample(N,N); - } + } - # initialize G according to the goal - G = evaluateGoalFullMatrix(M,goal) - - pred = TRUE - iter = 0 - while (pred){ - iter = iter + 1; - if(iter >= reps){ - pred = FALSE; - } - #iterate over all columns and swap all pairs of rows in each column - for(j in 1:d){ - # initialize new goal vector, for easier calculation we flatten the 2D matrix to a 1D vector - # where the index of the improvement vector corresponds to the row pair (k,l) with k < l - # so the 1D vector index = l + (k-1)*N - new_goal_vector = matrix(0,rows = 1, cols = N*(N-1)) - for(k in 1:(N-1)){ - for(l in (k+1):(N)){ - # add new goal to the goal vector according to the goal - newM = swap(M,j,k,l) - new_goal_vector[1, l + (k-1)*N] = G - evaluateGoalFullMatrix(newM,goal) - } - } - #if there was an improvement, swap the rows with the best improvement - if(min(new_goal_vector) < G){ - #calculate k and l from the index of the goal vector - idx = as.scalar(rowIndexMax(new_goal_vector)) - row1 = ((idx - 1) %/% N) + 1 - row2 = ((idx - 1) %% N) + 1 - M = swap(M,j,row1,row2) - - } - } - newG = evaluateGoalFullMatrix(M,goal) - # stop if no significant improvement was made, otherwise update G - # we can use abs() here, because we only update M if there was an improvement - if(abs(G - newG) < eps){ - pred = FALSE - } - else{ - G = newG - } - } - } - # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html - else if (method == "genetic") { - population = matrix(0, rows=genetic_population_size*N, cols=d) - scores = matrix(0, rows=genetic_population_size, cols=1) - - # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix - for (p in 1:genetic_population_size) { - design = lhs(N = N, d = d, method = "random", reps = reps) - population_start = designStart(p, N) - population_end = designEnd(p, N) - population[population_start:population_end,] = design - - # 2. Calculate the S-optimality score for each design - scores[p,1] = sOptimality(design) - } + newG = evaluateGoalFullMatrix(M, goal) + # stop if no significant improvement was made, otherwise update G + # we can use abs() here, because we only update M if there was an improvement + if (abs(G - newG) < eps) { + pred = FALSE + } + else { + G = newG + } + } +} + +geneticLHS = function(Integer N, Integer d, Integer reps, String goal, Integer genetic_generations, Integer genetic_population_size, Double genetic_mutation_rate) +return (Matrix[Double] M) { + population = matrix(0, rows = genetic_population_size * N, cols = d) + scores = matrix(0, rows = genetic_population_size, cols = 1) + + # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix + for (p in 1:genetic_population_size) { + design = lhs(N = N, d = d, method = "random", reps = reps) + population_start = designStart(p, N) + population_end = designEnd(p, N) + population[population_start:population_end, ] = design + + # 2. Calculate the S-optimality score for each design + scores[p, 1] = sOptimality(design) + } - for (generation in 1:genetic_generations) { - # 3. Put the best design in the first position of survivors and throw away the worst half of the remaining population - survivors_cutoff = floor(genetic_population_size / 2) + 1 - survivors = matrix(0, rows=survivors_cutoff * N, cols=d) - for (survivor in 1:survivors_cutoff) { - current_best_index = as.scalar(rowIndexMax(t(scores))) - current_best_design = getDesignByIndex(population, current_best_index, N) - - # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again - survivors_start = designStart(survivor, N) - survivors_end = designEnd(survivor, N) - survivors[survivors_start:survivors_end,] = current_best_design - scores[current_best_index, 1] = -1/0 - } - - # 4. Create new designs from the survivors - new_population = matrix(0, rows=genetic_population_size * N, cols=d) - best_design = survivors[1:N,] - new_population[1:N,] = best_design - next_free_slot = 2 - - # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. - for (survivor_index in 2:survivors_cutoff) { - new_design = getDesignByIndex(survivors, survivor_index, N) - source_col = as.scalar(sample(d, 1)) - target_col = as.scalar(sample(d, 1)) - new_design[,target_col] = best_design[,source_col] - - # 5. Attempt to mutate the child design by switching two elements in a random column. - new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) - - # Store the generated child design in the new population - new_population_start = designStart(next_free_slot, N) - new_population_end = designEnd(next_free_slot, N) - new_population[new_population_start:new_population_end,] = new_design - next_free_slot = next_free_slot + 1 - } - - # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. - for (survivor_index in 2:survivors_cutoff) { - if (next_free_slot <= genetic_population_size) { - source_design = getDesignByIndex(survivors, survivor_index, N) - new_design = best_design - source_col = as.scalar(sample(d, 1)) - target_col = as.scalar(sample(d, 1)) - new_design[,target_col] = source_design[,source_col] - - # 5. Attempt to mutate the child design by switching two elements in a random column. - new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) - - # Store the generated child design in the new population - new_population_start = designStart(next_free_slot, N) - new_population_end = designEnd(next_free_slot, N) - new_population[new_population_start:new_population_end,] = new_design - next_free_slot = next_free_slot + 1 - } - } - - population = new_population - - # Score the new population - for (p in 1:genetic_population_size) { - design = getDesignByIndex(population, p, N) - scores[p,1] = sOptimality(design) - } - } + for (generation in 1:genetic_generations) { + # 3. Put the best design in the first position of survivors and throw away the worst half of the remaining population + survivors_cutoff = floor(genetic_population_size / 2) + 1 + survivors = matrix(0, rows = survivors_cutoff * N, cols = d) + for (survivor in 1:survivors_cutoff) { + current_best_index = as.scalar(rowIndexMax(t(scores))) + current_best_design = getDesignByIndex(population, current_best_index, N) + + # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again + survivors_start = designStart(survivor, N) + survivors_end = designEnd(survivor, N) + survivors[survivors_start:survivors_end, ] = current_best_design + scores[current_best_index, 1] = -1 / 0 + } - # Return the best design from the final population - best_index = as.scalar(rowIndexMax(t(scores))) - M = getDesignByIndex(population, best_index, N) - } + # 4. Create new designs from the survivors + new_population = matrix(0, rows = genetic_population_size * N, cols = d) + best_design = survivors[1:N, ] + new_population[1:N, ] = best_design + next_free_slot = 2 + + # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. + for (survivor_index in 2:survivors_cutoff) { + new_design = getDesignByIndex(survivors, survivor_index, N) + source_col = as.scalar(sample(d, 1)) + target_col = as.scalar(sample(d, 1)) + new_design[, target_col] = best_design[, source_col] + + # 5. Attempt to mutate the child design by switching two elements in a random column. + new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + + # Store the generated child design in the new population + new_population_start = designStart(next_free_slot, N) + new_population_end = designEnd(next_free_slot, N) + new_population[new_population_start:new_population_end, ] = new_design + next_free_slot = next_free_slot + 1 + } - if(return_type == "distribution"){ - r = rand(rows = N, cols = d, min = 0, max = 1/N) - M = M/N - r + # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. + for (survivor_index in 2:survivors_cutoff) { + if (next_free_slot <= genetic_population_size) { + source_design = getDesignByIndex(survivors, survivor_index, N) + new_design = best_design + source_col = as.scalar(sample(d, 1)) + target_col = as.scalar(sample(d, 1)) + new_design[, target_col] = source_design[, source_col] + + # 5. Attempt to mutate the child design by switching two elements in a random column. + new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + + # Store the generated child design in the new population + new_population_start = designStart(next_free_slot, N) + new_population_end = designEnd(next_free_slot, N) + new_population[new_population_start:new_population_end, ] = new_design + next_free_slot = next_free_slot + 1 } } + + population = new_population + + # Score the new population + for (p in 1:genetic_population_size) { + design = getDesignByIndex(population, p, N) + scores[p, 1] = sOptimality(design) + } + } + + # Return the best design from the final population + best_index = as.scalar(rowIndexMax(t(scores))) + M = getDesignByIndex(population, best_index, N) +} + +minLhsDistance = function(Matrix[Double] M) +return (Double min_dist) { + D = dist(M) + diag(matrix(1 / 0, rows = nrow(M), cols = 1)) + min_dist = min(D) +} + +swap = function(Matrix[Double] M, Integer col, Integer row1, Integer row2) +return (Matrix[Double] M) { + temp = M[row1, col] + M[row1, col] = M[row2, col] + M[row2, col] = temp +} + +sOptimality = function(Matrix[Double] M) +return (Double score) { + D = dist(M) + score = sum(D) / (nrow(M) * (nrow(M) - 1)) +} + +designStart = function(Integer index, Integer N) +return (Integer start) { + start = (index - 1) * N + 1 +} + +designEnd = function(Integer index, Integer N) +return (Integer end) { + end = index * N +} + +getDesignByIndex = function(Matrix[Double] population, Integer index, Integer N) +return (Matrix[Double] design) { + start = designStart(index, N) + end = designEnd(index, N) + design = population[start:end, ] +} + +attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) +return (Matrix[Double] design) { + if (as.scalar(rand(rows = 1, cols = 1, min = 0, max = 1)) < mutation_rate) { + column_to_mutate = as.scalar(sample(ncol(design), 1)) + rows_to_mutate = sample(nrow(design), 2) + row1 = as.scalar(rows_to_mutate[1, 1]) + row2 = as.scalar(rows_to_mutate[2, 1]) + design = swap(design, column_to_mutate, row1, row2) + } +} + +sumInvDistance = function(Matrix[Double] M) +return (Double sum_inv) { + inv_dist = 1 / lower.tri(target = dist(M), diag = FALSE, values = TRUE) + # replace infinite values with 0 and sum the inverse distances + inv_dist = replace(target = inv_dist, pattern = 1 / 0, replacement = 0) + sum_inv = sum(inv_dist) +} + +evaluateGoalFullMatrix = function(Matrix[Double] M, String goal) +return (Double G) { + if (goal == "maximin") { + G = -minLhsDistance(M) + } + else if (goal == "opt") { + N = nrow(M) + d = ncol(M) + opt = N / (N^(1 / d)) + G = abs(opt - minLhsDistance(M)) + } + else if (goal == "sum_inv") { + G = sumInvDistance(M) + } + else if (goal == "avg_dist") { + G = -sOptimality(M) + } +} + +evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal) + return(Double G){ + if(goal == "maximin"){ + G = -min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)) + } + else if(goal == "opt"){ + N = nrow(M) + d = ncol(M) + opt2 = (N/(N^(1/d)))^2 + G = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2))) + } + else if(goal == "sum_inv"){ + G = sumInvDistance(rbind(M,row)) + } + else if(goal == "avg_dist"){ + G = -sOptimality(rbind(M,row)) + } +} \ No newline at end of file From 5700701da96ad61fdf22556b4b9df7a35a227ad9 Mon Sep 17 00:00:00 2001 From: = <=> Date: Thu, 9 Jul 2026 18:57:16 +0200 Subject: [PATCH 09/20] fix buggy behaviour in cpSweep method and evaluateGoalNewRow function update docstring --- scripts/builtin/lhs.dml | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 5eed04670ca..a4b8d8b69b4 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -45,8 +45,8 @@ # genetic_population_size Size of the population ("genetic" method) # genetic_mutation_rate Mutation rate, between 0 and 1 ("genetic" method) # return_type How should the matrix be returned: -# - "coordinates": return the matrix of vectors with point integer point coordinates in the d-dimensional space -# - "distribution": return the matrix as vectors sampled from the d-dimensional uniform distribution in range 0-1 +# - "buckets": return the matrix with bucket ids +# - "numerical": return the matrix as vectors sampled from the d-dimensional uniform distribution in range 0-1 # ------------------------------------------------------------------------------------------------- # # OUTPUT: @@ -157,7 +157,7 @@ return (Matrix[Double] M) { # check if the candidate is better then the current best # according to the goal and update the row vector and replacement mask accordingly - newG = evaluateGoalNewRow(M, row_candidate, goal) + newG = evaluateGoalNewRow(M, row_candidate, goal, N) if (newG < G) { row = row_candidate G = newG @@ -189,7 +189,6 @@ return (Matrix[Double] M) { # initialize G according to the goal G = evaluateGoalFullMatrix(M, goal) - pred = TRUE iter = 0 while (pred) { @@ -197,40 +196,38 @@ return (Matrix[Double] M) { if (iter >= reps) { pred = FALSE } - + oldG = G # iterate over all columns and swap all pairs of rows in each column for (j in 1:d) { + # initialize new goal vector, for easier calculation we flatten the 2D matrix to a 1D vector # where the index of the improvement vector corresponds to the row pair (k,l) with k < l # so the 1D vector index = l + (k-1)*N - new_goal_vector = matrix(0, rows = 1, cols = N * (N - 1)) + new_goal_vector = matrix(1/0, rows = 1, cols = N * (N - 1)) for (k in 1:(N - 1)) { for (l in (k + 1):(N)) { # add new goal to the goal vector according to the goal newM = swap(M, j, k, l) - new_goal_vector[1, l + (k - 1) * N] = G - evaluateGoalFullMatrix(newM, goal) + new_goal_vector[1, l + (k - 1) * N] = evaluateGoalFullMatrix(newM, goal) } } # if there was an improvement, swap the rows with the best improvement if (min(new_goal_vector) < G) { # calculate k and l from the index of the goal vector - idx = as.scalar(rowIndexMax(new_goal_vector)) + idx = as.scalar(rowIndexMin(new_goal_vector)) row1 = ((idx - 1) %/% N) + 1 row2 = ((idx - 1) %% N) + 1 M = swap(M, j, row1, row2) + G = min(new_goal_vector) } } - newG = evaluateGoalFullMatrix(M, goal) # stop if no significant improvement was made, otherwise update G # we can use abs() here, because we only update M if there was an improvement - if (abs(G - newG) < eps) { + if (abs(oldG - G) < eps) { pred = FALSE } - else { - G = newG - } } } @@ -396,16 +393,15 @@ return (Double G) { } } -evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal) +evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal, Integer N) return(Double G){ if(goal == "maximin"){ G = -min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)) } else if(goal == "opt"){ - N = nrow(M) - d = ncol(M) - opt2 = (N/(N^(1/d)))^2 - G = abs(opt2 - min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2))) + d = ncol(row) + opt = (N/(N^(1/d))) + G = abs(opt - sqrt(min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)))) } else if(goal == "sum_inv"){ G = sumInvDistance(rbind(M,row)) From 267e636c27f953b88324d292d00e02c1e84a6e22 Mon Sep 17 00:00:00 2001 From: = <=> Date: Sat, 11 Jul 2026 14:10:46 +0200 Subject: [PATCH 10/20] change lhs to return the optimization objective value --- scripts/builtin/lhs.dml | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index a4b8d8b69b4..43068160173 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -52,7 +52,7 @@ # OUTPUT: # ------------------------------------------------------------------------------------------------------ # M A N x d matrix M -# Contains the coordinates of sampled points +# G Value of the given objective # ------------------------------------------------------------------------------------------------------ lhs = function( @@ -67,7 +67,7 @@ lhs = function( Double genetic_mutation_rate = 0.25, String return_type = "buckets" ) -return (Matrix[Double] M) { +return (Matrix[Double] M, Double G) { if (goal != "maximin" & goal != "opt" & goal != "sum_inv" & goal != "avg_dist") { stop("Goal should be one of: maximin, opt, sum_inv, avg_dist") } @@ -91,7 +91,8 @@ return (Matrix[Double] M) { else { stop("Method should be one of: random, build, cp_sweep, genetic") } - + + G = evaluateGoalFullMatrix(M,goal) if (return_type == "numerical") { r = rand(rows = N, cols = d, min = 0, max = 1 / N) M = M / N - r @@ -118,9 +119,6 @@ return (Matrix[Double] M) { buildLHS = function(Integer N, Integer d, Integer reps, String goal) return (Matrix[Double] M) { - # optimal distance - opt2 = (N / (N^(1 / d)))^2 - # available buckets per column free_buckets = seq(1, N) %*% matrix(1, rows = 1, cols = d) @@ -178,8 +176,6 @@ return (Matrix[Double] M) { # TODO: vectorize the inner loops, currently this method is very slow cpSweepLHS = function(Integer N, Integer d, Integer reps, String goal, Double eps) return (Matrix[Double] M) { - # optimal distance - opt2 = (N / (N^(1 / d)))^2 # generate initial random latin hypercube M = matrix(0, rows = N, cols = d) @@ -244,7 +240,7 @@ return (Matrix[Double] M) { population[population_start:population_end, ] = design # 2. Calculate the S-optimality score for each design - scores[p, 1] = sOptimality(design) + scores[p, 1] = avgDistance(design) } for (generation in 1:genetic_generations) { @@ -310,7 +306,7 @@ return (Matrix[Double] M) { # Score the new population for (p in 1:genetic_population_size) { design = getDesignByIndex(population, p, N) - scores[p, 1] = sOptimality(design) + scores[p, 1] = avgDistance(design) } } @@ -332,7 +328,7 @@ return (Matrix[Double] M) { M[row2, col] = temp } -sOptimality = function(Matrix[Double] M) +avgDistance = function(Matrix[Double] M) return (Double score) { D = dist(M) score = sum(D) / (nrow(M) * (nrow(M) - 1)) @@ -389,7 +385,7 @@ return (Double G) { G = sumInvDistance(M) } else if (goal == "avg_dist") { - G = -sOptimality(M) + G = -avgDistance(M) } } @@ -407,6 +403,6 @@ evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal G = sumInvDistance(rbind(M,row)) } else if(goal == "avg_dist"){ - G = -sOptimality(rbind(M,row)) + G = -avgDistance(rbind(M,row)) } } \ No newline at end of file From ff6454ca62ab07e9531000cd1d5a9e39e9723685 Mon Sep 17 00:00:00 2001 From: = <=> Date: Sat, 11 Jul 2026 14:11:59 +0200 Subject: [PATCH 11/20] add script to benchmark lhs againt random sampling by fitting a gmm to sampled data --- .../functions/builtin/lhs_gmm_benchmark.dml | 80 +++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml diff --git a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml new file mode 100644 index 00000000000..b72eb30e60f --- /dev/null +++ b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml @@ -0,0 +1,80 @@ +d = 5 +method = ifdef($method, "random") +goal = ifdef($goal, "maximin") +reps = ifdef($reps, 10) +seed = ifdef($seed, -1) +fname = ifdef($fname, "lhsExperimentOutput.csv") +Ntrain = list(10,15,20,30,50,75,100,200,500) +Ntest = 1000000 + +getParameterErrors = function(Matrix[Double] trueMu, Matrix[Double] trueSigma, Matrix[Double] mu, Matrix[Double] precCholesky) +return (Double muDistanceError, Double covMatrixError){ + # calculate euclidean distance of estimated and true mu + error = trueMu - mu + muDistanceError = sqrt(sum(error*error)) + + #calculate Frobenius norm of the difference between estimated and true sigma + estimatedSigma = inv(precCholesky %*% t(precCholesky)) + sigmaDiff = trueSigma - estimatedSigma + covMatrixError = sqrt(sum(sigmaDiff*sigmaDiff)) +} + +getMeanLogLikelihood = function(Matrix[Double] testData, Matrix[Double] mu, Matrix[Double] precCholesky) +return (Double meanLogLikelihood){ + # calculate the mean log likelihood of test data originating from the estimated distribution + y = testData %*% precCholesky - mu %*% precCholesky + log_prob = rowSums(y * y) + log_det_chol = sum(log(diag(t(precCholesky)))) + es_log_prob = -.5 * (ncol(testData) * log(2 * pi) + log_prob) + log_det_chol + meanLogLikelihood = mean(es_log_prob) +} + + +# setup true distribution parameters +trueMu = rand(rows = 1, cols = d, min = -10 , max = 10, seed = seed) +A = rand(rows = d, cols = d, min = 0, max = 5, seed = seed) +# ensures the matrix is symmetric and psd +trueSigma = A %*% t(A) + diag(matrix(0.01, rows = d, cols = d)) +L = cholesky(trueSigma) + +# generate test data +Z = rand(rows = Ntest, cols = 5, pdf = "normal", seed = seed) +testData = (Z %*% t(L)) + (matrix(1, rows = Ntest, cols = 1) %*% trueMu) + + +results = matrix(0,rows = length(Ntrain), cols = 9) +for(i in 1:length(Ntrain)){ + N = as.scalar(Ntrain[i]) + # random sampling + Z = rand(rows = N, cols = d, pdf = "normal", seed = seed) + # transform to the test multivariate normal distribution + trainDataRandom = (Z %*% t(L)) + (matrix(1, rows = N, cols = 1) %*% trueMu) + #fit the gmm with one component + [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataRandom, n_components = 1, model = "VVV") + [muDistanceErrorRandom, covMatrixErrorRandom] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) + trainLlRandom = getMeanLogLikelihood(trainDataRandom, mu, precCholesky) + testLlRandom = getMeanLogLikelihood(testData, mu, precCholesky) + # lhs sampling + [LHS,G] = lhs(N = N, d = d, reps = reps, method = method, goal = goal, return_type = "numerical") + + # transform the uniform data into N(0,1) oisdhdofisdfion + for(j in 1:d){ + for(k in 1:N){ + p = as.scalar(LHS[k,j]) + q = qnorm(target = p, mean = 0, std = 1) + LHS[k,j] = q + } + } + # transform data into the true multivariate normal distribution + trainDataLHS = (LHS %*% t(L)) + (matrix(1, rows = N, cols = 1) %*% trueMu) + [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataLHS, n_components = 1, model = "VVV") + [muDistanceErrorLHS, covMatrixErrorLHS] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) + trainLlLHS = getMeanLogLikelihood(trainDataLHS, mu, precCholesky) + testLlLHS = getMeanLogLikelihood(testData, mu, precCholesky) + results[i,] = cbind(as.matrix(N), as.matrix(muDistanceErrorRandom), as.matrix(covMatrixErrorRandom), + as.matrix(trainLlRandom), as.matrix(testLlRandom), + as.matrix(muDistanceErrorLHS), as.matrix(covMatrixErrorLHS), + as.matrix(trainLlLHS), as.matrix(testLlLHS)) +} + +write(results, fname, format="csv", sep=";"); \ No newline at end of file From 98830b60744c71854f538394259c6b6787fa39fd Mon Sep 17 00:00:00 2001 From: = <=> Date: Mon, 13 Jul 2026 21:59:28 +0200 Subject: [PATCH 12/20] add kl divergence to benchmarking metrics --- .../functions/builtin/lhs_gmm_benchmark.dml | 26 ++++++++++++++++--- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml index b72eb30e60f..cc016f07e4b 100644 --- a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml +++ b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml @@ -4,7 +4,7 @@ goal = ifdef($goal, "maximin") reps = ifdef($reps, 10) seed = ifdef($seed, -1) fname = ifdef($fname, "lhsExperimentOutput.csv") -Ntrain = list(10,15,20,30,50,75,100,200,500) +Ntrain = list(10,13,15,20,25,30,40,50,75,100) Ntest = 1000000 getParameterErrors = function(Matrix[Double] trueMu, Matrix[Double] trueSigma, Matrix[Double] mu, Matrix[Double] precCholesky) @@ -19,6 +19,22 @@ return (Double muDistanceError, Double covMatrixError){ covMatrixError = sqrt(sum(sigmaDiff*sigmaDiff)) } +getKLDivergence = function(Matrix[Double] trueMu, Matrix[Double] trueSigma, Matrix[Double] mu, Matrix[Double] precCholesky, Double d) +return (Double Dkl){ + # calculate kullback leibler divergence between true and estimate distribution + # Dkl(P||Q) + muDiff = trueMu - mu + invEstimatedSigma = precCholesky %*% t(precCholesky) + term1 = trace(invEstimatedSigma %*% trueSigma) + term2 = as.scalar(muDiff%*%invEstimatedSigma%*%t(muDiff)) + + logDetEst = -2*sum(log(diag(precCholesky))) + L = cholesky(trueSigma) + logDetTrue = 2*sum(log(diag(L))) + term3 = logDetEst - logDetTrue + Dkl = 0.5 * (term1 + term2 + term3 - d) +} + getMeanLogLikelihood = function(Matrix[Double] testData, Matrix[Double] mu, Matrix[Double] precCholesky) return (Double meanLogLikelihood){ # calculate the mean log likelihood of test data originating from the estimated distribution @@ -42,7 +58,7 @@ Z = rand(rows = Ntest, cols = 5, pdf = "normal", seed = seed) testData = (Z %*% t(L)) + (matrix(1, rows = Ntest, cols = 1) %*% trueMu) -results = matrix(0,rows = length(Ntrain), cols = 9) +results = matrix(0,rows = length(Ntrain), cols = 11) for(i in 1:length(Ntrain)){ N = as.scalar(Ntrain[i]) # random sampling @@ -52,6 +68,7 @@ for(i in 1:length(Ntrain)){ #fit the gmm with one component [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataRandom, n_components = 1, model = "VVV") [muDistanceErrorRandom, covMatrixErrorRandom] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) + DklRandom = getKLDivergence(trueMu, trueSigma, mu, precCholesky, d) trainLlRandom = getMeanLogLikelihood(trainDataRandom, mu, precCholesky) testLlRandom = getMeanLogLikelihood(testData, mu, precCholesky) # lhs sampling @@ -69,11 +86,12 @@ for(i in 1:length(Ntrain)){ trainDataLHS = (LHS %*% t(L)) + (matrix(1, rows = N, cols = 1) %*% trueMu) [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataLHS, n_components = 1, model = "VVV") [muDistanceErrorLHS, covMatrixErrorLHS] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) + DklLHS = getKLDivergence(trueMu, trueSigma, mu, precCholesky, d) trainLlLHS = getMeanLogLikelihood(trainDataLHS, mu, precCholesky) testLlLHS = getMeanLogLikelihood(testData, mu, precCholesky) - results[i,] = cbind(as.matrix(N), as.matrix(muDistanceErrorRandom), as.matrix(covMatrixErrorRandom), + results[i,] = cbind(as.matrix(N), as.matrix(muDistanceErrorRandom), as.matrix(covMatrixErrorRandom), as.matrix(DklRandom), as.matrix(trainLlRandom), as.matrix(testLlRandom), - as.matrix(muDistanceErrorLHS), as.matrix(covMatrixErrorLHS), + as.matrix(muDistanceErrorLHS), as.matrix(covMatrixErrorLHS), as.matrix(DklLHS), as.matrix(trainLlLHS), as.matrix(testLlLHS)) } From 91e03f0220ab137d3ee8c3a3f7de6906a82199e3 Mon Sep 17 00:00:00 2001 From: grzybkos <83309779+grzybkos@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:20:52 +0200 Subject: [PATCH 13/20] use the goal parameter for scoring in genetic LHS --- scripts/builtin/lhs.dml | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 43068160173..c1e99c3c61e 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -234,13 +234,13 @@ return (Matrix[Double] M) { # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix for (p in 1:genetic_population_size) { - design = lhs(N = N, d = d, method = "random", reps = reps) + design = lhs(N = N, d = d, method = "random", reps = reps, goal = goal) population_start = designStart(p, N) population_end = designEnd(p, N) population[population_start:population_end, ] = design - # 2. Calculate the S-optimality score for each design - scores[p, 1] = avgDistance(design) + # 2. Calculate the score for each design according to the goal + scores[p, 1] = evaluateGoalFullMatrix(design, goal) } for (generation in 1:genetic_generations) { @@ -248,14 +248,14 @@ return (Matrix[Double] M) { survivors_cutoff = floor(genetic_population_size / 2) + 1 survivors = matrix(0, rows = survivors_cutoff * N, cols = d) for (survivor in 1:survivors_cutoff) { - current_best_index = as.scalar(rowIndexMax(t(scores))) + current_best_index = as.scalar(rowIndexMin(t(scores))) current_best_design = getDesignByIndex(population, current_best_index, N) - # Put the current best design in the survivors matrix, then set its score to -inf so that it is not selected again + # Put the current best design in the survivors matrix, then set its score to +inf (worst) so that it is not selected again survivors_start = designStart(survivor, N) survivors_end = designEnd(survivor, N) survivors[survivors_start:survivors_end, ] = current_best_design - scores[current_best_index, 1] = -1 / 0 + scores[current_best_index, 1] = 1 / 0 } # 4. Create new designs from the survivors @@ -303,15 +303,15 @@ return (Matrix[Double] M) { population = new_population - # Score the new population + # Score the new population according to the goal for (p in 1:genetic_population_size) { design = getDesignByIndex(population, p, N) - scores[p, 1] = avgDistance(design) + scores[p, 1] = evaluateGoalFullMatrix(design, goal) } } # Return the best design from the final population - best_index = as.scalar(rowIndexMax(t(scores))) + best_index = as.scalar(rowIndexMin(t(scores))) M = getDesignByIndex(population, best_index, N) } From ab6398f771b1004ec0493204d12af514bda55a33 Mon Sep 17 00:00:00 2001 From: grzybkos <83309779+grzybkos@users.noreply.github.com> Date: Tue, 14 Jul 2026 14:32:29 +0200 Subject: [PATCH 14/20] adjust LHS script params and outputs --- src/test/scripts/functions/builtin/lhs.dml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml index dae502ab0c2..af08c661326 100644 --- a/src/test/scripts/functions/builtin/lhs.dml +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -24,5 +24,9 @@ d = $2 reps = $3 method = $4 goal = $5 -M = lhs(N = N,d = d,reps = reps, method = method, goal = goal) -write(M, $6) +generations = $6 +population_size = $7 +[M, G] = lhs(N = N,d = d,reps = reps, method = method, goal = goal, + genetic_generations = generations, genetic_population_size = population_size) +write(M, $8) +write(G, $9) From 8083156291a154545cacdce58548faf920435fae Mon Sep 17 00:00:00 2001 From: grzybkos <83309779+grzybkos@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:19:18 +0200 Subject: [PATCH 15/20] test all LHS method and parameter combinations, explain failures --- .../builtin/part1/BuiltinLHSTest.java | 67 ++++++++----------- 1 file changed, 29 insertions(+), 38 deletions(-) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index d03701a4c88..f9a5ce0d722 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -44,13 +44,14 @@ public void checkLatinHypercubeValidity(HashMap m,int N, int double val = m.get(new CellIndex(row, col)); int intVal = (int) val; - assertTrue(intVal >= 1 && intVal <= N); - - assertFalse(seen[intVal]); + assertTrue("value " + val + " at row " + row + " col " + col + " is outside the bucket range (1 to " + N + ")", + intVal >= 1 && intVal <= N); + + assertFalse("duplicate bucket " + intVal + " in column " + col, seen[intVal]); seen[intVal] = true; } for (int i = 1; i <= N; i++) { - assertTrue(seen[i]); + assertTrue("bucket " + i + " never used in column " + col, seen[i]); } } } @@ -60,55 +61,45 @@ public void setUp() { addTestConfiguration(TEST_NAME, new TestConfiguration(TEST_CLASS_DIR, TEST_NAME,new String[]{"B"})); } - - @Test - public void testRandomTwoDim() { - runLhsTest(5,2,5,"random","opt"); - } - - @Test - public void testRandomMultiDim() { - runLhsTest(10,4,5,"random","opt"); + @Test + public void testRandomTwoDim() { + runLhsTest(5, 2, 5, "random", "opt", 0, 0); } @Test - public void testImproved() { - runLhsTest(10,4,5,"build","opt"); - } - @Test - public void testSumInv() { - runLhsTest(10,4,5,"build","sum_inv"); - } - @Test - public void testMaximin() { - runLhsTest(10,4,5,"build","maximin"); + public void testAllMethodGoalCombinations() { + String[] methods = {"random", "build", "cp_sweep", "genetic"}; + String[] goals = {"maximin", "opt", "sum_inv", "avg_dist"}; + for (String method : methods) { + for (String goal : goals) { + try { + runLhsTest(10, 4, 5, method, goal, 3, 5); + } + catch (Throwable t) { + throw new AssertionError("failed for method=" + method + ", goal=" + goal + ": " + t.getMessage(), t); + } + } + } } - @Test - public void testCPsweep() { - runLhsTest(10,4,5,"cp_sweep","sum_inv"); - } - @Test - public void testGenetic() { - runLhsTest(10,4,5,"genetic","sum_inv"); - } - private void runLhsTest(int N,int d, int repetitions, String method, String goal) - { + private void runLhsTest(int N,int d, int repetitions, String method, String goal, int generations, int populationSize) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); - try - { + try { loadTestConfiguration(getTestConfiguration(TEST_NAME)); String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; - programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method, goal, output("C")}; + programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method, goal, + Integer.toString(generations), Integer.toString(populationSize), output("C"), output("G")}; - //execute test runTest(true, false, null, -1); - HashMap m = readDMLMatrixFromOutputDir("C"); checkLatinHypercubeValidity(m,N,d); + + HashMap g = readDMLScalarFromOutputDir("G"); + Double G = g.get(new CellIndex(1, 1)); + assertTrue("G should be a finite number but was: " + G, G != null && !G.isNaN() && !G.isInfinite()); } finally { rtplatform = platformOld; From 4e1420a80f5c8d6850db999df406ee9f2f594974 Mon Sep 17 00:00:00 2001 From: grzybkos <83309779+grzybkos@users.noreply.github.com> Date: Tue, 14 Jul 2026 17:47:00 +0200 Subject: [PATCH 16/20] refactor child generation in genetic LHS method --- scripts/builtin/lhs.dml | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index c1e99c3c61e..8c9f2dd2cc4 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -264,15 +264,10 @@ return (Matrix[Double] M) { new_population[1:N, ] = best_design next_free_slot = 2 - # 4.1: First: use each other survivor as base and replace one of its columns with a random column from the best design. + # 4.1: use each other survivor as base and transplant a column from the best design for (survivor_index in 2:survivors_cutoff) { - new_design = getDesignByIndex(survivors, survivor_index, N) - source_col = as.scalar(sample(d, 1)) - target_col = as.scalar(sample(d, 1)) - new_design[, target_col] = best_design[, source_col] - - # 5. Attempt to mutate the child design by switching two elements in a random column. - new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + survivor_design = getDesignByIndex(survivors, survivor_index, N) + new_design = makeGeneticChild(survivor_design, best_design, genetic_mutation_rate) # Store the generated child design in the new population new_population_start = designStart(next_free_slot, N) @@ -281,17 +276,11 @@ return (Matrix[Double] M) { next_free_slot = next_free_slot + 1 } - # 4.2: For the empty positions, use the best design as base and replace one column with one from a survivor. + # 4.2: for the empty positions, use the best design as base and transplant a column from a survivor for (survivor_index in 2:survivors_cutoff) { if (next_free_slot <= genetic_population_size) { - source_design = getDesignByIndex(survivors, survivor_index, N) - new_design = best_design - source_col = as.scalar(sample(d, 1)) - target_col = as.scalar(sample(d, 1)) - new_design[, target_col] = source_design[, source_col] - - # 5. Attempt to mutate the child design by switching two elements in a random column. - new_design = attemptGeneticMutation(new_design, genetic_mutation_rate) + survivor_design = getDesignByIndex(survivors, survivor_index, N) + new_design = makeGeneticChild(best_design, survivor_design, genetic_mutation_rate) # Store the generated child design in the new population new_population_start = designStart(next_free_slot, N) @@ -351,6 +340,18 @@ return (Matrix[Double] design) { design = population[start:end, ] } +# Creates a child design by transplanting a random column from a donor design +# into the base design and attempting a mutation. +makeGeneticChild = function(Matrix[Double] base, Matrix[Double] donor, Double mutation_rate) +return (Matrix[Double] child) { + source_col = as.scalar(sample(ncol(base), 1)) + target_col = as.scalar(sample(ncol(base), 1)) + child = base + child[, target_col] = donor[, source_col] + # Attempt to mutate the child design by switching two elements in a random column. + child = attemptGeneticMutation(child, mutation_rate) +} + attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) return (Matrix[Double] design) { if (as.scalar(rand(rows = 1, cols = 1, min = 0, max = 1)) < mutation_rate) { From da68d710bf4f879bc4afaef3d3b5dc9f1e2bbc26 Mon Sep 17 00:00:00 2001 From: grzybkos <83309779+grzybkos@users.noreply.github.com> Date: Tue, 14 Jul 2026 20:25:31 +0200 Subject: [PATCH 17/20] add seed parameter to LHS builtin and improve LHS tests --- scripts/builtin/lhs.dml | 62 +++++++++++-------- .../builtin/part1/BuiltinLHSTest.java | 41 +++++++++--- src/test/scripts/functions/builtin/lhs.dml | 7 ++- 3 files changed, 75 insertions(+), 35 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 8c9f2dd2cc4..5118122d532 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -47,6 +47,7 @@ # return_type How should the matrix be returned: # - "buckets": return the matrix with bucket ids # - "numerical": return the matrix as vectors sampled from the d-dimensional uniform distribution in range 0-1 +# seed The seed for making results reproducible. Defaults to -1, which creates random seeds. # ------------------------------------------------------------------------------------------------- # # OUTPUT: @@ -65,7 +66,8 @@ lhs = function( Integer genetic_generations = 5, Integer genetic_population_size = 10, Double genetic_mutation_rate = 0.25, - String return_type = "buckets" + String return_type = "buckets", + Integer seed = -1 ) return (Matrix[Double] M, Double G) { if (goal != "maximin" & goal != "opt" & goal != "sum_inv" & goal != "avg_dist") { @@ -76,17 +78,17 @@ return (Matrix[Double] M, Double G) { } if (method == "random") { - M = randomLHS(N, d, reps, goal) + M = randomLHS(N, d, reps, goal, seed) } else if (method == "build") { - M = buildLHS(N, d, reps, goal) + M = buildLHS(N, d, reps, goal, seed) } else if (method == "cp_sweep") { - M = cpSweepLHS(N, d, reps, goal, eps) + M = cpSweepLHS(N, d, reps, goal, eps, seed) } # Reference: https://bertcarnell.github.io/lhs/reference/geneticLHS.html else if (method == "genetic") { - M = geneticLHS(N, d, reps, goal, genetic_generations, genetic_population_size, genetic_mutation_rate) + M = geneticLHS(N, d, reps, goal, genetic_generations, genetic_population_size, genetic_mutation_rate, seed) } else { stop("Method should be one of: random, build, cp_sweep, genetic") @@ -94,19 +96,19 @@ return (Matrix[Double] M, Double G) { G = evaluateGoalFullMatrix(M,goal) if (return_type == "numerical") { - r = rand(rows = N, cols = d, min = 0, max = 1 / N) + r = rand(rows = N, cols = d, min = 0, max = 1 / N, seed = seed) M = M / N - r } } -randomLHS = function(Integer N, Integer d, Integer reps, String goal) +randomLHS = function(Integer N, Integer d, Integer reps, String goal, Integer seed) return (Matrix[Double] M) { G = 1 / 0 M = matrix(0, rows = N, cols = d) for (i in 1:reps) { newM = matrix(0, rows = N, cols = d) for (j in 1:d) { - newM[, j] = sample(N, N) + newM[, j] = sample(N, N, deriveSeed(seed, (i - 1) * d + j)) } newG = evaluateGoalFullMatrix(newM, goal) @@ -117,13 +119,13 @@ return (Matrix[Double] M) { } } -buildLHS = function(Integer N, Integer d, Integer reps, String goal) +buildLHS = function(Integer N, Integer d, Integer reps, String goal, Integer seed) return (Matrix[Double] M) { # available buckets per column free_buckets = seq(1, N) %*% matrix(1, rows = 1, cols = d) # generate first point at random - first_sample = sample(N, d, TRUE) + first_sample = sample(N, d, TRUE, deriveSeed(seed, 1)) M = t(first_sample) # replace used ints for each column, with last sample @@ -147,7 +149,7 @@ return (Matrix[Double] M) { for (j in 1:reps) { # random mask for bucket generation - rand = table(sample(num_free_buckets, d, TRUE), seq(1, N), N, d) + rand = table(sample(num_free_buckets, d, TRUE, deriveSeed(seed, i * reps + j)), seq(1, N), N, d) # get bucket and flatten to row vector vals_matrix = rand * free_buckets @@ -174,13 +176,13 @@ return (Matrix[Double] M) { } # TODO: vectorize the inner loops, currently this method is very slow -cpSweepLHS = function(Integer N, Integer d, Integer reps, String goal, Double eps) +cpSweepLHS = function(Integer N, Integer d, Integer reps, String goal, Double eps, Integer seed) return (Matrix[Double] M) { # generate initial random latin hypercube M = matrix(0, rows = N, cols = d) for (i in 1:d) { - M[, i] = sample(N, N) + M[, i] = sample(N, N, deriveSeed(seed, i)) } # initialize G according to the goal @@ -227,14 +229,16 @@ return (Matrix[Double] M) { } } -geneticLHS = function(Integer N, Integer d, Integer reps, String goal, Integer genetic_generations, Integer genetic_population_size, Double genetic_mutation_rate) +geneticLHS = function(Integer N, Integer d, Integer reps, String goal, Integer genetic_generations, Integer genetic_population_size, Double genetic_mutation_rate, Integer seed) return (Matrix[Double] M) { population = matrix(0, rows = genetic_population_size * N, cols = d) scores = matrix(0, rows = genetic_population_size, cols = 1) + population_seed_offsets = genetic_population_size * reps * d + # 1. Generate initial random Latin Hypercube designs and store them all in the population matrix for (p in 1:genetic_population_size) { - design = lhs(N = N, d = d, method = "random", reps = reps, goal = goal) + design = lhs(N = N, d = d, method = "random", reps = reps, goal = goal, seed = deriveSeed(seed, (p - 1) * reps * d)) population_start = designStart(p, N) population_end = designEnd(p, N) population[population_start:population_end, ] = design @@ -266,8 +270,9 @@ return (Matrix[Double] M) { # 4.1: use each other survivor as base and transplant a column from the best design for (survivor_index in 2:survivors_cutoff) { + child_base = population_seed_offsets + ((generation - 1) * genetic_population_size + next_free_slot) * 5 survivor_design = getDesignByIndex(survivors, survivor_index, N) - new_design = makeGeneticChild(survivor_design, best_design, genetic_mutation_rate) + new_design = makeGeneticChild(survivor_design, best_design, genetic_mutation_rate, deriveSeed(seed, child_base)) # Store the generated child design in the new population new_population_start = designStart(next_free_slot, N) @@ -279,8 +284,9 @@ return (Matrix[Double] M) { # 4.2: for the empty positions, use the best design as base and transplant a column from a survivor for (survivor_index in 2:survivors_cutoff) { if (next_free_slot <= genetic_population_size) { + child_base = population_seed_offsets + ((generation - 1) * genetic_population_size + next_free_slot) * 5 survivor_design = getDesignByIndex(survivors, survivor_index, N) - new_design = makeGeneticChild(best_design, survivor_design, genetic_mutation_rate) + new_design = makeGeneticChild(best_design, survivor_design, genetic_mutation_rate, deriveSeed(seed, child_base)) # Store the generated child design in the new population new_population_start = designStart(next_free_slot, N) @@ -323,6 +329,12 @@ return (Double score) { score = sum(D) / (nrow(M) * (nrow(M) - 1)) } +# Returns seed base + offset, but keeps the -1 (random) seed unchanged. +deriveSeed = function(Integer seed, Integer offset) +return (Integer s) { + s = ifelse(seed == -1, -1, seed + offset) +} + designStart = function(Integer index, Integer N) return (Integer start) { start = (index - 1) * N + 1 @@ -342,21 +354,21 @@ return (Matrix[Double] design) { # Creates a child design by transplanting a random column from a donor design # into the base design and attempting a mutation. -makeGeneticChild = function(Matrix[Double] base, Matrix[Double] donor, Double mutation_rate) +makeGeneticChild = function(Matrix[Double] base, Matrix[Double] donor, Double mutation_rate, Integer child_seed) return (Matrix[Double] child) { - source_col = as.scalar(sample(ncol(base), 1)) - target_col = as.scalar(sample(ncol(base), 1)) + source_col = as.scalar(sample(ncol(base), 1, deriveSeed(child_seed, 1))) + target_col = as.scalar(sample(ncol(base), 1, deriveSeed(child_seed, 2))) child = base child[, target_col] = donor[, source_col] # Attempt to mutate the child design by switching two elements in a random column. - child = attemptGeneticMutation(child, mutation_rate) + child = attemptGeneticMutation(child, mutation_rate, deriveSeed(child_seed, 3)) } -attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate) +attemptGeneticMutation = function(Matrix[Double] design, Double mutation_rate, Integer seed) return (Matrix[Double] design) { - if (as.scalar(rand(rows = 1, cols = 1, min = 0, max = 1)) < mutation_rate) { - column_to_mutate = as.scalar(sample(ncol(design), 1)) - rows_to_mutate = sample(nrow(design), 2) + if (as.scalar(rand(rows = 1, cols = 1, min = 0, max = 1, seed = seed)) < mutation_rate) { + column_to_mutate = as.scalar(sample(ncol(design), 1, deriveSeed(seed, 1))) + rows_to_mutate = sample(nrow(design), 2, deriveSeed(seed, 2)) row1 = as.scalar(rows_to_mutate[1, 1]) row2 = as.scalar(rows_to_mutate[2, 1]) design = swap(design, column_to_mutate, row1, row2) diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java index f9a5ce0d722..9f44d10775b 100644 --- a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinLHSTest.java @@ -29,6 +29,7 @@ import org.apache.sysds.common.Types.ExecMode; import org.apache.sysds.test.AutomatedTestBase; import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; public class BuiltinLHSTest extends AutomatedTestBase @@ -63,7 +64,7 @@ public void setUp() { @Test public void testRandomTwoDim() { - runLhsTest(5, 2, 5, "random", "opt", 0, 0); + runLhsTest(5, 2, 5, "random", "opt", 0, 0, -1); } @Test @@ -73,7 +74,7 @@ public void testAllMethodGoalCombinations() { for (String method : methods) { for (String goal : goals) { try { - runLhsTest(10, 4, 5, method, goal, 3, 5); + runLhsTest(10, 4, 5, method, goal, 3, 5, -1); } catch (Throwable t) { throw new AssertionError("failed for method=" + method + ", goal=" + goal + ": " + t.getMessage(), t); @@ -82,7 +83,33 @@ public void testAllMethodGoalCombinations() { } } - private void runLhsTest(int N,int d, int repetitions, String method, String goal, int generations, int populationSize) { + @Test + public void testReproducibility() { + String[] methods = {"random", "build", "cp_sweep", "genetic"}; + for (String method : methods) { + LhsResult r1 = runLhsTest(10, 4, 5, method, "maximin", 3, 5, 321); + LhsResult r2 = runLhsTest(10, 4, 5, method, "maximin", 3, 5, 321); + TestUtils.compareMatrices(r1.m(), r2.m(), 0, "first-run-" + method, "second-run-" + method); + } + } + + @Test + public void testGeneticOptimizesGoal() { + // "genetic" can't be worse than "random" with the same seed + // this primarily tests whether the scoring logic is correct + String[] goals = {"maximin", "opt", "sum_inv", "avg_dist"}; + for (String goal : goals) { + double gGenetic = runLhsTest(10, 4, 1, "genetic", goal, 10, 20, 1234).g(); + double gRandom = runLhsTest(10, 4, 1, "random", goal, 10, 20, 1234).g(); + assertTrue("genetic G=" + gGenetic + " should not be worse than random G=" + gRandom + " for goal " + goal, + gGenetic <= gRandom); + } + } + + private record LhsResult(HashMap m, double g) {} + + private LhsResult runLhsTest(int N,int d, int repetitions, String method, String goal, + int generations, int populationSize, int seed) { ExecMode platformOld = setExecMode(ExecMode.HYBRID); try { @@ -90,16 +117,16 @@ private void runLhsTest(int N,int d, int repetitions, String method, String goal String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; programArgs = new String[]{"-args", Integer.toString(N), Integer.toString(d), Integer.toString(repetitions), method, goal, - Integer.toString(generations), Integer.toString(populationSize), output("C"), output("G")}; + Integer.toString(generations), Integer.toString(populationSize), Integer.toString(seed), output("C"), output("G")}; runTest(true, false, null, -1); HashMap m = readDMLMatrixFromOutputDir("C"); checkLatinHypercubeValidity(m,N,d); - HashMap g = readDMLScalarFromOutputDir("G"); - Double G = g.get(new CellIndex(1, 1)); - assertTrue("G should be a finite number but was: " + G, G != null && !G.isNaN() && !G.isInfinite()); + double G = readDMLScalarFromOutputDir("G").get(new CellIndex(1, 1)); + assertTrue("G should be a finite number but was: " + G, !Double.isNaN(G) && !Double.isInfinite(G)); + return new LhsResult(m, G); } finally { rtplatform = platformOld; diff --git a/src/test/scripts/functions/builtin/lhs.dml b/src/test/scripts/functions/builtin/lhs.dml index af08c661326..95961a03cf9 100644 --- a/src/test/scripts/functions/builtin/lhs.dml +++ b/src/test/scripts/functions/builtin/lhs.dml @@ -26,7 +26,8 @@ method = $4 goal = $5 generations = $6 population_size = $7 +seed = $8 [M, G] = lhs(N = N,d = d,reps = reps, method = method, goal = goal, - genetic_generations = generations, genetic_population_size = population_size) -write(M, $8) -write(G, $9) + genetic_generations = generations, genetic_population_size = population_size, seed = seed) +write(M, $9) +write(G, $10) From 8db9be55b1a803ea85d679a56ad75a4723fab0b0 Mon Sep 17 00:00:00 2001 From: = <=> Date: Tue, 14 Jul 2026 23:15:23 +0200 Subject: [PATCH 18/20] fix sumInv to return inverse of squared distances fix minor bug in build method remove unnecessary metrics in benchmark script --- scripts/builtin/lhs.dml | 33 ++++++++++--------- .../functions/builtin/lhs_gmm_benchmark.dml | 28 ++++------------ 2 files changed, 24 insertions(+), 37 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index 5118122d532..e89c1eeb5af 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -34,11 +34,11 @@ # - for "random" method, the best sample is chosen from reps random samples # - for "build" method, reps candidates are generated for each point and the best candidate is chosen # - for "cp_sweep" method, reps sweeps are performed for each column -# - for "genetic" method, reps candidates are generate in the initial design +# - for "genetic" method, reps candidates are generated in the initial design # goal The goal to be minimized during the sampling procedure, possible values are: # - "maximin" - maximize the minimum distance between points (in practice the negative minimum distance is minimized) # - "opt" - minimize the difference between the minimum distance and the optimal distance N/(N^(1/d)) -# - "sum_inv" - minimize the sum of the inverse distances between points +# - "sum_inv" - minimize the sum of the inverse squared distances between points # - "avg_dist" - maximize avg euclidean distance between points (in practice the negative is minimized) # eps The threshold for stopping the cp_sweep method # genetic_generations Number of generations ("genetic" method) @@ -93,12 +93,13 @@ return (Matrix[Double] M, Double G) { else { stop("Method should be one of: random, build, cp_sweep, genetic") } - + G = evaluateGoalFullMatrix(M,goal) if (return_type == "numerical") { r = rand(rows = N, cols = d, min = 0, max = 1 / N, seed = seed) M = M / N - r } + } randomLHS = function(Integer N, Integer d, Integer reps, String goal, Integer seed) @@ -130,17 +131,17 @@ return (Matrix[Double] M) { # replace used ints for each column, with last sample # in subsequent runs we only look at the first N-1,N-2,... rows - replacement_mask = table(first_sample, seq(1, N), N, d) + replacement_mask = table(first_sample, seq(1, d), N, d) free_buckets = (1 - replacement_mask) * free_buckets + (matrix(1, rows = N, cols = 1) %*% (free_buckets[N, ])) * replacement_mask # for the next N - 2 points: - # generate reps candidates and choose the candidate, - # for which the minimal distance to the already chosen points is maximal + # generate reps candidates + # choose the candidate with the lowest G # and update the free_buckets matrix for (i in 2:(N - 1)) { # number of buckets left to choose from num_free_buckets = N - i + 1 - # initialize current best distance to 0 for maximin and to infinity for opt and sum_inv + # initialize G to infinity G = 1 / 0 # initialize row vector and replacement mask for free buckets @@ -149,9 +150,8 @@ return (Matrix[Double] M) { for (j in 1:reps) { # random mask for bucket generation - rand = table(sample(num_free_buckets, d, TRUE, deriveSeed(seed, i * reps + j)), seq(1, N), N, d) - - # get bucket and flatten to row vector + rand = table(sample(num_free_buckets, d, TRUE, deriveSeed(seed, i * reps + j)), seq(1, d), N, d) + # get buckets and flatten to row vector vals_matrix = rand * free_buckets row_candidate = matrix(1, rows = 1, cols = N) %*% vals_matrix @@ -175,7 +175,6 @@ return (Matrix[Double] M) { M = rbind(M, free_buckets[1, ]) } -# TODO: vectorize the inner loops, currently this method is very slow cpSweepLHS = function(Integer N, Integer d, Integer reps, String goal, Double eps, Integer seed) return (Matrix[Double] M) { @@ -185,7 +184,6 @@ return (Matrix[Double] M) { M[, i] = sample(N, N, deriveSeed(seed, i)) } - # initialize G according to the goal G = evaluateGoalFullMatrix(M, goal) pred = TRUE iter = 0 @@ -199,7 +197,7 @@ return (Matrix[Double] M) { for (j in 1:d) { # initialize new goal vector, for easier calculation we flatten the 2D matrix to a 1D vector - # where the index of the improvement vector corresponds to the row pair (k,l) with k < l + # where the index of the new goal vector corresponds to the row pair (k,l) with k < l # so the 1D vector index = l + (k-1)*N new_goal_vector = matrix(1/0, rows = 1, cols = N * (N - 1)) for (k in 1:(N - 1)) { @@ -222,7 +220,6 @@ return (Matrix[Double] M) { } # stop if no significant improvement was made, otherwise update G - # we can use abs() here, because we only update M if there was an improvement if (abs(oldG - G) < eps) { pred = FALSE } @@ -377,8 +374,10 @@ return (Matrix[Double] design) { sumInvDistance = function(Matrix[Double] M) return (Double sum_inv) { - inv_dist = 1 / lower.tri(target = dist(M), diag = FALSE, values = TRUE) - # replace infinite values with 0 and sum the inverse distances + row_sq = rowSums(M ^ 2) + ones = matrix(1, rows=nrow(M), cols=1) + squared_dist = (row_sq %*% t(ones)) + (ones %*% t(row_sq)) - (2 * M %*% t(M)) + inv_dist = 1 / lower.tri(target = squared_dist, diag = FALSE, values = TRUE) inv_dist = replace(target = inv_dist, pattern = 1 / 0, replacement = 0) sum_inv = sum(inv_dist) } @@ -405,11 +404,13 @@ return (Double G) { evaluateGoalNewRow = function (Matrix[Double] M, Matrix[Double] row, String goal, Integer N) return(Double G){ if(goal == "maximin"){ + # sqrt not needed G = -min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)) } else if(goal == "opt"){ d = ncol(row) opt = (N/(N^(1/d))) + # sqrt is necessary here G = abs(opt - sqrt(min(rowSums((M - (matrix(1, rows = nrow(M), cols = 1) %*% row))^2)))) } else if(goal == "sum_inv"){ diff --git a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml index cc016f07e4b..80bb72dfbea 100644 --- a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml +++ b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml @@ -7,18 +7,6 @@ fname = ifdef($fname, "lhsExperimentOutput.csv") Ntrain = list(10,13,15,20,25,30,40,50,75,100) Ntest = 1000000 -getParameterErrors = function(Matrix[Double] trueMu, Matrix[Double] trueSigma, Matrix[Double] mu, Matrix[Double] precCholesky) -return (Double muDistanceError, Double covMatrixError){ - # calculate euclidean distance of estimated and true mu - error = trueMu - mu - muDistanceError = sqrt(sum(error*error)) - - #calculate Frobenius norm of the difference between estimated and true sigma - estimatedSigma = inv(precCholesky %*% t(precCholesky)) - sigmaDiff = trueSigma - estimatedSigma - covMatrixError = sqrt(sum(sigmaDiff*sigmaDiff)) -} - getKLDivergence = function(Matrix[Double] trueMu, Matrix[Double] trueSigma, Matrix[Double] mu, Matrix[Double] precCholesky, Double d) return (Double Dkl){ # calculate kullback leibler divergence between true and estimate distribution @@ -58,7 +46,7 @@ Z = rand(rows = Ntest, cols = 5, pdf = "normal", seed = seed) testData = (Z %*% t(L)) + (matrix(1, rows = Ntest, cols = 1) %*% trueMu) -results = matrix(0,rows = length(Ntrain), cols = 11) +results = matrix(0,rows = length(Ntrain), cols = 7) for(i in 1:length(Ntrain)){ N = as.scalar(Ntrain[i]) # random sampling @@ -67,12 +55,11 @@ for(i in 1:length(Ntrain)){ trainDataRandom = (Z %*% t(L)) + (matrix(1, rows = N, cols = 1) %*% trueMu) #fit the gmm with one component [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataRandom, n_components = 1, model = "VVV") - [muDistanceErrorRandom, covMatrixErrorRandom] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) DklRandom = getKLDivergence(trueMu, trueSigma, mu, precCholesky, d) trainLlRandom = getMeanLogLikelihood(trainDataRandom, mu, precCholesky) testLlRandom = getMeanLogLikelihood(testData, mu, precCholesky) # lhs sampling - [LHS,G] = lhs(N = N, d = d, reps = reps, method = method, goal = goal, return_type = "numerical") + [LHS,G] = lhs(N = N, d = d, reps = reps, method = method, goal = goal, return_type = "numerical", seed = seed) # transform the uniform data into N(0,1) oisdhdofisdfion for(j in 1:d){ @@ -85,14 +72,13 @@ for(i in 1:length(Ntrain)){ # transform data into the true multivariate normal distribution trainDataLHS = (LHS %*% t(L)) + (matrix(1, rows = N, cols = 1) %*% trueMu) [labels, predict_prob, df, bic, mu, precCholesky, weight] = gmm(X = trainDataLHS, n_components = 1, model = "VVV") - [muDistanceErrorLHS, covMatrixErrorLHS] = getParameterErrors(trueMu, trueSigma, mu, precCholesky) DklLHS = getKLDivergence(trueMu, trueSigma, mu, precCholesky, d) - trainLlLHS = getMeanLogLikelihood(trainDataLHS, mu, precCholesky) + trainLlLHS = getMeanLogLikelihood(trainDataLHS, mu, precCholesky) testLlLHS = getMeanLogLikelihood(testData, mu, precCholesky) - results[i,] = cbind(as.matrix(N), as.matrix(muDistanceErrorRandom), as.matrix(covMatrixErrorRandom), as.matrix(DklRandom), - as.matrix(trainLlRandom), as.matrix(testLlRandom), - as.matrix(muDistanceErrorLHS), as.matrix(covMatrixErrorLHS), as.matrix(DklLHS), - as.matrix(trainLlLHS), as.matrix(testLlLHS)) + results[i,] = cbind(as.matrix(N), + as.matrix(DklRandom), as.matrix(trainLlRandom), as.matrix(testLlRandom), + as.matrix(DklLHS), as.matrix(trainLlLHS), as.matrix(testLlLHS) + ) } write(results, fname, format="csv", sep=";"); \ No newline at end of file From 20d877a2b5425a26248bbdc2d205e7edd4b2a722 Mon Sep 17 00:00:00 2001 From: = <=> Date: Wed, 15 Jul 2026 12:35:30 +0200 Subject: [PATCH 19/20] fix random sampling dimension and comment in benchmark script --- src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml index 80bb72dfbea..ad40ada930c 100644 --- a/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml +++ b/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml @@ -42,7 +42,7 @@ trueSigma = A %*% t(A) + diag(matrix(0.01, rows = d, cols = d)) L = cholesky(trueSigma) # generate test data -Z = rand(rows = Ntest, cols = 5, pdf = "normal", seed = seed) +Z = rand(rows = Ntest, cols = d, pdf = "normal", seed = seed) testData = (Z %*% t(L)) + (matrix(1, rows = Ntest, cols = 1) %*% trueMu) @@ -61,7 +61,7 @@ for(i in 1:length(Ntrain)){ # lhs sampling [LHS,G] = lhs(N = N, d = d, reps = reps, method = method, goal = goal, return_type = "numerical", seed = seed) - # transform the uniform data into N(0,1) oisdhdofisdfion + # transform the uniform data into N(0,1) for(j in 1:d){ for(k in 1:N){ p = as.scalar(LHS[k,j]) From 8bff6b41a593ff6c4dba0ff625e2134fce5635db Mon Sep 17 00:00:00 2001 From: = <=> Date: Wed, 15 Jul 2026 13:09:06 +0200 Subject: [PATCH 20/20] add stop condition for numerical arguments make build method work for N=2 --- scripts/builtin/lhs.dml | 76 +++++++++++++++++++++++++++-------------- 1 file changed, 50 insertions(+), 26 deletions(-) diff --git a/scripts/builtin/lhs.dml b/scripts/builtin/lhs.dml index e89c1eeb5af..806ece9f81d 100644 --- a/scripts/builtin/lhs.dml +++ b/scripts/builtin/lhs.dml @@ -76,6 +76,29 @@ return (Matrix[Double] M, Double G) { if (return_type != "buckets" & return_type != "numerical") { stop("Return type should be one of: buckets, numerical") } + if (N < 2) { + stop("N should be a positive integer larger than 1") + } + + if (d < 1) { + stop("d should be a positive integer") + } + + if (reps < 1) { + stop("reps should be a positive integer") + } + + if (method == "genetic" & genetic_generations < 1) { + stop("genetic_generations should be >= 1") + } + + if (method == "genetic" & genetic_population_size < 2) { + stop("genetic_population_size should be >= 2") + } + + if (method == "genetic" & (genetic_mutation_rate < 0 | genetic_mutation_rate > 1)) { + stop("genetic_mutation_rate should be between 0 and 1") + } if (method == "random") { M = randomLHS(N, d, reps, goal, seed) @@ -138,37 +161,38 @@ return (Matrix[Double] M) { # generate reps candidates # choose the candidate with the lowest G # and update the free_buckets matrix - for (i in 2:(N - 1)) { - # number of buckets left to choose from - num_free_buckets = N - i + 1 - # initialize G to infinity - G = 1 / 0 - - # initialize row vector and replacement mask for free buckets - row = matrix(0, rows = 1, cols = d) - replacement_mask = matrix(0, rows = num_free_buckets, cols = d) - - for (j in 1:reps) { - # random mask for bucket generation - rand = table(sample(num_free_buckets, d, TRUE, deriveSeed(seed, i * reps + j)), seq(1, d), N, d) - # get buckets and flatten to row vector - vals_matrix = rand * free_buckets - row_candidate = matrix(1, rows = 1, cols = N) %*% vals_matrix - - # check if the candidate is better then the current best - # according to the goal and update the row vector and replacement mask accordingly - newG = evaluateGoalNewRow(M, row_candidate, goal, N) - if (newG < G) { - row = row_candidate - G = newG - replacement_mask = rand + if (N>2){ + for (i in 2:(N - 1)) { + # number of buckets left to choose from + num_free_buckets = N - i + 1 + # initialize G to infinity + G = 1 / 0 + + # initialize row vector and replacement mask for free buckets + row = matrix(0, rows = 1, cols = d) + replacement_mask = matrix(0, rows = num_free_buckets, cols = d) + + for (j in 1:reps) { + # random mask for bucket generation + rand = table(sample(num_free_buckets, d, TRUE, deriveSeed(seed, i * reps + j)), seq(1, d), N, d) + # get buckets and flatten to row vector + vals_matrix = rand * free_buckets + row_candidate = matrix(1, rows = 1, cols = N) %*% vals_matrix + + # check if the candidate is better then the current best + # according to the goal and update the row vector and replacement mask accordingly + newG = evaluateGoalNewRow(M, row_candidate, goal, N) + if (newG < G) { + row = row_candidate + G = newG + replacement_mask = rand + } } - } - # add the best candidate to the matrix of chosen points M = rbind(M, row) # replace chosen buckets with last, not yet picked buckets free_buckets = (1 - replacement_mask) * free_buckets + (matrix(1, rows = N, cols = 1) %*% (free_buckets[num_free_buckets, ])) * replacement_mask + } } # for the last point there is only one bucket left