From 941f61c12d70e30d5805943e84be9cbb4b3cc71a Mon Sep 17 00:00:00 2001 From: nofa Date: Mon, 12 Jan 2026 18:03:31 +0500 Subject: [PATCH 1/5] initial setup --- .../com/jdoe/algorithms/CompositeDOE.java | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 src/main/java/com/jdoe/algorithms/CompositeDOE.java diff --git a/src/main/java/com/jdoe/algorithms/CompositeDOE.java b/src/main/java/com/jdoe/algorithms/CompositeDOE.java new file mode 100644 index 0000000..4256525 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/CompositeDOE.java @@ -0,0 +1,38 @@ +package com.jdoe.algorithms; + +public class CompositeDOE { + + public static void centralCompositeDesign() { + + } + +} + +//<-----------------------------STEPS----------------------------------> + +/** + * Detailed Steps to Port Central Composite Design Script to Java 1. Create CompositeDOE Class Structure Create CompositeDOE.java in + * com.jdoe.algorithms package Add imports: org.apache.commons.math3.linear.* for matrix operations Define method signature: public static + * RealMatrix centralCompositeDesign(int n, int[] center, String alpha, String face) Set default parameters handling for center, alpha, and + * face 2. Implement Input Validation Logic Add assertion: assert n > 1 : "\"n\" must be an integer greater than 1" Validate alpha parameter + * with: alpha.toLowerCase() and check against allowed values Validate face parameter with: face.toLowerCase() and check against allowed + * values Validate center array length: if (center.length != 2) throw IllegalArgumentException Format error messages to match Python + * implementation 3. Implement Star Point Generation Create StarDOE.java class with star method Calculate alpha based on alpha type: For + * orthogonal: alpha = Math.pow(2 * factorial(n), 0.25) where factorial(n) = n! For rotatable: alpha = Math.pow(2, 0.5) for 2D, Math.pow(3, + * 0.5) for 3D, etc. Generate 2n star points with values [±alpha, 0, ..., 0], [0, ±alpha, ..., 0], etc. Return both star matrix and alpha + * value as an Object array 4. Enhance FactorialDOE Class Add ff2n method that generates 2-level factorial design Create 2^n × n matrix with + * values -1 and +1 Use bit manipulation: for each row i and column j, set value to ((i >> j) & 1) == 0 ? -1 : 1 Return RealMatrix object 5. + * Create UnionDOE Class Implement union method that combines two matrices vertically Use Array2DRowRealMatrix constructor to create new + * matrix with combined rows Copy data from both input matrices to the result matrix Return the concatenated matrix 6. Create + * RepeatCenterDOE Class Implement repeatCenter method that generates center points Create repeats × n matrix filled with zeros Use new + * Array2DRowRealMatrix(repeats, n) to initialize matrix Return the center point matrix 7. Implement Core Algorithm Logic Initialize H1 and + * H2 matrices based on face type For inscribed (cci): scale factorial points H1 = H1.scalarMultiply(1.0/a) For faced (ccf): set alpha to + * 1.0 For circumscribed (ccc): keep original factorial points Generate center points C1 and C2 using repeatCenter Combine matrices: H1 = + * union(H1, C1), H2 = union(H2, C2), H = union(H1, H2) 8. Handle Matrix Operations Implement matrix scaling using scalarMultiply() method + * from Commons Math Ensure all matrix dimensions are compatible for union operations Use getSubMatrix() and setSubMatrix() for matrix + * manipulations Return final RealMatrix object 9. Create Unit Tests Create CompositeDOETest.java with comprehensive test cases Test all + * three face types: ccc, cci, ccf Test both alpha types: orthogonal, rotatable Validate matrix dimensions: (2^n + 2n + sum(center)) × n + * Verify center point counts and star point distances 10. Add Documentation and Error Handling Add JavaDoc explaining parameters, return + * value, and exceptions Implement proper exception handling for invalid inputs Include example usage in documentation Follow the same error + * message format as Python implementation + */ \ No newline at end of file From fc0072ae61ee865297bc93e71beb864483d48eed Mon Sep 17 00:00:00 2001 From: nofa Date: Thu, 15 Jan 2026 16:28:11 +0500 Subject: [PATCH 2/5] composite design changes --- pom.xml | 13 ++ .../com/jdoe/algorithms/CompositeDOE.java | 119 ++++++++++++++---- .../java/com/jdoe/util/FactorialUtility.java | 3 + .../java/com/jdoe/util/GenericDOEUtil.java | 22 ++++ 4 files changed, 134 insertions(+), 23 deletions(-) create mode 100644 src/main/java/com/jdoe/util/GenericDOEUtil.java diff --git a/pom.xml b/pom.xml index 1648392..5934c85 100644 --- a/pom.xml +++ b/pom.xml @@ -54,6 +54,19 @@ 5.2.0 test + + org.jetbrains + annotations + 26.0.2 + compile + + + + org.projectlombok + lombok + 1.18.42 + compile + diff --git a/src/main/java/com/jdoe/algorithms/CompositeDOE.java b/src/main/java/com/jdoe/algorithms/CompositeDOE.java index 4256525..ee1f6f5 100644 --- a/src/main/java/com/jdoe/algorithms/CompositeDOE.java +++ b/src/main/java/com/jdoe/algorithms/CompositeDOE.java @@ -1,9 +1,68 @@ package com.jdoe.algorithms; +import org.apache.commons.math3.linear.Array2DRowFieldMatrix; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.jetbrains.annotations.NotNull; + +import com.jdoe.util.GenericDOEUtil; + public class CompositeDOE { - public static void centralCompositeDesign() { + private static final String[] alphaParamLegalExpressions = { "orthogonal", "o", "rotatable", "r" }; + + private static final String[] faceParamLegalExpressions = { "circumscribed", "ccc", "inscribed", "cci", "faced", "ccf" }; + + public static void centralCompositeDesign( @NotNull int numberOfFactors, @NotNull int[] centerPoints, @NotNull String alpha, + @NotNull String face ) { + + validateCompositeDegisnParameters( numberOfFactors, centerPoints, alpha, face ); + + RealMatrix fullFactorial2LevelMatrix = FactorialDOE.fullFactorial2Level( numberOfFactors ); + + } + + //<-----------------------------Utility Methods----------------------------------> + + public static void validateCompositeDegisnParameters( int numberOfFactors, int[] centerPoints, String alpha, String face ) { + if ( numberOfFactors < 1 ) { + throw new IllegalArgumentException( "number of factors must be greater then 1" ); + } + if ( centerPoints.length == 2 ) { + throw new IllegalArgumentException( + String.format( "Invalid number of values for \"center\" (expected 2, but got %d )", centerPoints.length ) ); + } + if ( validateExpression( alpha, alphaParamLegalExpressions ) == Boolean.FALSE ) { + throw new IllegalArgumentException( + String.format( "Invalid value for \"alpha\" (expected format %d )", alphaParamLegalExpressions.toString() ) ); + } + if ( validateExpression( alpha, faceParamLegalExpressions ) == Boolean.FALSE ) { + throw new IllegalArgumentException( + String.format( "Invalid value for \"face\" (expected format %d )", faceParamLegalExpressions.toString() ) ); + } + } + + public static Boolean validateExpression( String expression, String[] legalExpressions ) { + Boolean validInput = Boolean.FALSE; + for ( String legalExpression : legalExpressions ) { + if ( legalExpression.equals( expression ) ) { + validInput = Boolean.TRUE; + break; + } + } + return validInput; + } + public static Boolean matrixUnion( RealMatrix firstMatrix, RealMatrix secondMatrix ) { + double[][] finalMatrixRowModel = new double[ firstMatrix.getRowDimension() + secondMatrix.getRowDimension() ][]; + for ( int i = 0; i < firstMatrix.getRowDimension(); i++ ) { + finalMatrixRowModel[ i ] = firstMatrix.getRow( i ); + } + for ( int i = 0; i < secondMatrix.getRowDimension(); i++ ) { + finalMatrixRowModel[ i ] = secondMatrix.getRow( i ); + } + RealMatrix finalMatrix = new Array2DRowRealMatrix(finalMatrixRowModel); + GenericDOEUtil.matrixLogger( finalMatrix, "matrixUnion" ); } } @@ -11,28 +70,42 @@ public static void centralCompositeDesign() { //<-----------------------------STEPS----------------------------------> /** - * Detailed Steps to Port Central Composite Design Script to Java 1. Create CompositeDOE Class Structure Create CompositeDOE.java in - * com.jdoe.algorithms package Add imports: org.apache.commons.math3.linear.* for matrix operations Define method signature: public static - * RealMatrix centralCompositeDesign(int n, int[] center, String alpha, String face) Set default parameters handling for center, alpha, and - * face 2. Implement Input Validation Logic Add assertion: assert n > 1 : "\"n\" must be an integer greater than 1" Validate alpha parameter + * Detailed Steps to Port Central Composite Design Script to Java + * + * $ 1. Create CompositeDOE Class Structure Create CompositeDOE.java in com.jdoe.algorithms package Add imports: + * org.apache.commons.math3.linear.* for matrix operations Define method signature: public static RealMatrix centralCompositeDesign(int n, + * int[] center, String alpha, String face) Set default parameters handling for center, alpha, and face + * + * $ 2. Implement Input Validation Logic Add assertion: assert n > 1 : "\"n\" must be an integer greater than 1" Validate alpha parameter * with: alpha.toLowerCase() and check against allowed values Validate face parameter with: face.toLowerCase() and check against allowed * values Validate center array length: if (center.length != 2) throw IllegalArgumentException Format error messages to match Python - * implementation 3. Implement Star Point Generation Create StarDOE.java class with star method Calculate alpha based on alpha type: For - * orthogonal: alpha = Math.pow(2 * factorial(n), 0.25) where factorial(n) = n! For rotatable: alpha = Math.pow(2, 0.5) for 2D, Math.pow(3, - * 0.5) for 3D, etc. Generate 2n star points with values [±alpha, 0, ..., 0], [0, ±alpha, ..., 0], etc. Return both star matrix and alpha - * value as an Object array 4. Enhance FactorialDOE Class Add ff2n method that generates 2-level factorial design Create 2^n × n matrix with - * values -1 and +1 Use bit manipulation: for each row i and column j, set value to ((i >> j) & 1) == 0 ? -1 : 1 Return RealMatrix object 5. - * Create UnionDOE Class Implement union method that combines two matrices vertically Use Array2DRowRealMatrix constructor to create new - * matrix with combined rows Copy data from both input matrices to the result matrix Return the concatenated matrix 6. Create - * RepeatCenterDOE Class Implement repeatCenter method that generates center points Create repeats × n matrix filled with zeros Use new - * Array2DRowRealMatrix(repeats, n) to initialize matrix Return the center point matrix 7. Implement Core Algorithm Logic Initialize H1 and - * H2 matrices based on face type For inscribed (cci): scale factorial points H1 = H1.scalarMultiply(1.0/a) For faced (ccf): set alpha to - * 1.0 For circumscribed (ccc): keep original factorial points Generate center points C1 and C2 using repeatCenter Combine matrices: H1 = - * union(H1, C1), H2 = union(H2, C2), H = union(H1, H2) 8. Handle Matrix Operations Implement matrix scaling using scalarMultiply() method - * from Commons Math Ensure all matrix dimensions are compatible for union operations Use getSubMatrix() and setSubMatrix() for matrix - * manipulations Return final RealMatrix object 9. Create Unit Tests Create CompositeDOETest.java with comprehensive test cases Test all - * three face types: ccc, cci, ccf Test both alpha types: orthogonal, rotatable Validate matrix dimensions: (2^n + 2n + sum(center)) × n - * Verify center point counts and star point distances 10. Add Documentation and Error Handling Add JavaDoc explaining parameters, return - * value, and exceptions Implement proper exception handling for invalid inputs Include example usage in documentation Follow the same error - * message format as Python implementation + * implementation + * + * $ 3. Implement Star Point Generation Create StarDOE.java class with star method Calculate alpha based on alpha type: For orthogonal: + * alpha = Math.pow(2 * factorial(n), 0.25) where factorial(n) = n! For rotatable: alpha = Math.pow(2, 0.5) for 2D, Math.pow(3, 0.5) for 3D, + * etc. Generate 2n star points with values [±alpha, 0, ..., 0], [0, ±alpha, ..., 0], etc. Return both star matrix and alpha value as an + * Object array + * + * 4. Enhance FactorialDOE Class Add ff2n method that generates 2-level factorial design Create 2^n × n matrix with values -1 and +1 Use bit + * manipulation: for each row i and column j, set value to ((i >> j) & 1) == 0 ? -1 : 1 Return RealMatrix object + * + * 5. Create UnionDOE Class Implement union method that combines two matrices vertically Use Array2DRowRealMatrix constructor to create new + * matrix with combined rows Copy data from both input matrices to the result matrix Return the concatenated matrix + * + * 6. Create RepeatCenterDOE Class Implement repeatCenter method that generates center points Create repeats × n matrix filled with zeros + * Use new Array2DRowRealMatrix(repeats, n) to initialize matrix Return the center point matrix + * + * 7. Implement Core Algorithm Logic Initialize H1 andH2 matrices based on face type For inscribed (cci): scale factorial points H1 = + * H1.scalarMultiply(1.0/a) For faced (ccf): set alpha to 1.0 For circumscribed (ccc): keep original factorial points Generate center points + * C1 and C2 using repeatCenter Combine matrices: H1 = union(H1, C1), H2 = union(H2, C2), H = union(H1, H2) + * + * 8. Handle Matrix Operations Implement matrix scaling using scalarMultiply() method from Commons Math Ensure all matrix dimensions are + * compatible for union operations Use getSubMatrix() and setSubMatrix() for matrix manipulations Return final RealMatrix object + * + * 9. Create Unit Tests Create CompositeDOETest.java with comprehensive test cases Test all three face types: ccc, cci, ccf Test both alpha + * types: orthogonal, rotatable Validate matrix dimensions: (2^n + 2n + sum(center)) × n Verify center point counts and star point + * distances + * + * 10. Add Documentation and Error Handling Add JavaDoc explaining parameters, return value, and exceptions Implement proper exception + * handling for invalid inputs Include example usage in documentation Follow the same error message format as Python implementation */ \ No newline at end of file diff --git a/src/main/java/com/jdoe/util/FactorialUtility.java b/src/main/java/com/jdoe/util/FactorialUtility.java index 0ced350..484233c 100644 --- a/src/main/java/com/jdoe/util/FactorialUtility.java +++ b/src/main/java/com/jdoe/util/FactorialUtility.java @@ -3,6 +3,9 @@ import java.util.ArrayList; import java.util.List; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; + public class FactorialUtility { /** diff --git a/src/main/java/com/jdoe/util/GenericDOEUtil.java b/src/main/java/com/jdoe/util/GenericDOEUtil.java new file mode 100644 index 0000000..fc32c95 --- /dev/null +++ b/src/main/java/com/jdoe/util/GenericDOEUtil.java @@ -0,0 +1,22 @@ +package com.jdoe.util; + +import java.util.Arrays; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; + +public class GenericDOEUtil { + + private static final Logger log = LogManager.getLogger( FactorialUtility.class ); + + public static void matrixLogger( RealMatrix matrix, String designReference ){ + StringBuilder sb = new StringBuilder(); + sb.append(String.format( "Matrix for Design { %s }:\n" , designReference ) ); + for (int i = 0; i < matrix.getRowDimension(); i++) { + sb.append( Arrays.toString(matrix.getRow(i))).append("\n"); + } + log.info(sb.toString()); + } + +} From 12ef63c110bdaef1288dce8c33946ed1a885961a Mon Sep 17 00:00:00 2001 From: nomus Date: Sun, 25 Jan 2026 20:15:12 +0500 Subject: [PATCH 3/5] composite doe implementation --- dependency-reduced-pom.xml | 86 ++++++++ .../com/jdoe/algorithms/CompositeDOE.java | 117 +++++++---- .../com/jdoe/algorithms/RepeatCenterDOE.java | 12 ++ .../java/com/jdoe/algorithms/StarDOE.java | 63 ++++++ .../java/com/jdoe/algorithms/UnionDOE.java | 33 +++ .../com/jdoe/algorithms/CompositeDOETest.java | 195 ++++++++++++++++++ target/classes/log4j2.xml | 13 ++ .../compile/default-compile/createdFiles.lst | 6 + .../compile/default-compile/inputFiles.lst | 6 + .../default-testCompile/createdFiles.lst | 2 + .../default-testCompile/inputFiles.lst | 2 + ...T-com.jdoe.algorithms.FactorialDOETest.xml | 81 ++++++++ .../com.jdoe.algorithms.FactorialDOETest.txt | 4 + 13 files changed, 577 insertions(+), 43 deletions(-) create mode 100644 dependency-reduced-pom.xml create mode 100644 src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java create mode 100644 src/main/java/com/jdoe/algorithms/StarDOE.java create mode 100644 src/main/java/com/jdoe/algorithms/UnionDOE.java create mode 100644 src/test/java/com/jdoe/algorithms/CompositeDOETest.java create mode 100644 target/classes/log4j2.xml create mode 100644 target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst create mode 100644 target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst create mode 100644 target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst create mode 100644 target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst create mode 100644 target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml create mode 100644 target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt diff --git a/dependency-reduced-pom.xml b/dependency-reduced-pom.xml new file mode 100644 index 0000000..f1be5e4 --- /dev/null +++ b/dependency-reduced-pom.xml @@ -0,0 +1,86 @@ + + + 4.0.0 + com.jdoe + doe-genereter + JDOE + 1.0-SNAPSHOT + http://maven.apache.org + + + + maven-jar-plugin + 3.4.1 + + + + com.jdoe.Testing + + + + + + maven-shade-plugin + 3.6.1 + + + package + + shade + + + + + + + + + + + + + + junit + junit + 4.13.2 + test + + + hamcrest-core + org.hamcrest + + + + + org.mockito + mockito-core + 5.20.0 + test + + + byte-buddy + net.bytebuddy + + + byte-buddy-agent + net.bytebuddy + + + objenesis + org.objenesis + + + + + org.mockito + mockito-inline + 5.2.0 + test + + + + 17 + 17 + UTF-8 + + diff --git a/src/main/java/com/jdoe/algorithms/CompositeDOE.java b/src/main/java/com/jdoe/algorithms/CompositeDOE.java index ee1f6f5..f4051be 100644 --- a/src/main/java/com/jdoe/algorithms/CompositeDOE.java +++ b/src/main/java/com/jdoe/algorithms/CompositeDOE.java @@ -9,43 +9,86 @@ public class CompositeDOE { - private static final String[] alphaParamLegalExpressions = { "orthogonal", "o", "rotatable", "r" }; + private static final String[] alphaParamLegalExpressions = {"orthogonal", "o", "rotatable", "r"}; - private static final String[] faceParamLegalExpressions = { "circumscribed", "ccc", "inscribed", "cci", "faced", "ccf" }; + private static final String[] faceParamLegalExpressions = {"circumscribed", "ccc", "inscribed", "cci", "faced", "ccf"}; - public static void centralCompositeDesign( @NotNull int numberOfFactors, @NotNull int[] centerPoints, @NotNull String alpha, - @NotNull String face ) { + public static RealMatrix centralCompositeDesign(@NotNull int numberOfFactors, @NotNull int[] centerPoints, @NotNull String alpha, @NotNull String face) { + validateCompositeDegisnParameters(numberOfFactors, centerPoints, alpha, face); - validateCompositeDegisnParameters( numberOfFactors, centerPoints, alpha, face ); + RealMatrix factorialDesignMatrix; + RealMatrix starDesignMatrix; + double alphaValue; - RealMatrix fullFactorial2LevelMatrix = FactorialDOE.fullFactorial2Level( numberOfFactors ); + // Orthogonal Design + if (alpha.toLowerCase().equals("orthogonal") || alpha.toLowerCase().equals("o")) { + Object[] starResult = StarDOE.star(numberOfFactors, "orthogonal", centerPoints); + starDesignMatrix = new Array2DRowRealMatrix((double[][]) starResult[0]); // Convert array to RealMatrix + alphaValue = (double) starResult[1]; + } + // Rotatable Design + else if (alpha.toLowerCase().equals("rotatable") || alpha.toLowerCase().equals("r")) { + Object[] starResult = StarDOE.star(numberOfFactors, "rotatable", new int[]{0, 0}); + starDesignMatrix = new Array2DRowRealMatrix((double[][]) starResult[0]); // Convert array to RealMatrix + alphaValue = (double) starResult[1]; + } else { + throw new IllegalArgumentException("Invalid alpha parameter: " + alpha); + } + + // Inscribed CCD + if (face.toLowerCase().equals("inscribed") || face.toLowerCase().equals("cci")) { + factorialDesignMatrix = FactorialDOE.fullFactorial2Level(numberOfFactors); + factorialDesignMatrix = factorialDesignMatrix.scalarMultiply(1.0 / alphaValue); // Scale down the factorial points + Object[] starResult = StarDOE.star(numberOfFactors, alpha, centerPoints); + starDesignMatrix = new Array2DRowRealMatrix((double[][]) starResult[0]); // Convert array to RealMatrix + alphaValue = (double) starResult[1]; + } + // Faced CCD + else if (face.toLowerCase().equals("faced") || face.toLowerCase().equals("ccf")) { + Object[] starResult = StarDOE.star(numberOfFactors, alpha, centerPoints); + starDesignMatrix = new Array2DRowRealMatrix((double[][]) starResult[0]); // Convert array to RealMatrix + alphaValue = 1.0; + // Value of alpha is always 1 in Faced CCD + factorialDesignMatrix = FactorialDOE.fullFactorial2Level(numberOfFactors); + } + // Circumscribed CCD + else if (face.toLowerCase().equals("circumscribed") || face.toLowerCase().equals("ccc")) { + factorialDesignMatrix = FactorialDOE.fullFactorial2Level(numberOfFactors); + } else { + throw new IllegalArgumentException("Invalid face parameter: " + face); + } + RealMatrix centerPointsMatrix1 = RepeatCenterDOE.repeatCenter(numberOfFactors, centerPoints[0]); + RealMatrix centerPointsMatrix2 = RepeatCenterDOE.repeatCenter(numberOfFactors, centerPoints[1]); + + factorialDesignMatrix = UnionDOE.matrixUnion(factorialDesignMatrix, centerPointsMatrix1); + starDesignMatrix = UnionDOE.matrixUnion(starDesignMatrix, centerPointsMatrix2); + RealMatrix finalDesignMatrix = UnionDOE.matrixUnion(factorialDesignMatrix, starDesignMatrix); + + return finalDesignMatrix; } //<-----------------------------Utility Methods----------------------------------> - public static void validateCompositeDegisnParameters( int numberOfFactors, int[] centerPoints, String alpha, String face ) { - if ( numberOfFactors < 1 ) { - throw new IllegalArgumentException( "number of factors must be greater then 1" ); + public static void validateCompositeDegisnParameters(int numberOfFactors, int[] centerPoints, String alpha, String face) { + if (numberOfFactors < 1) { + throw new IllegalArgumentException("number of factors must be greater than 1"); } - if ( centerPoints.length == 2 ) { - throw new IllegalArgumentException( - String.format( "Invalid number of values for \"center\" (expected 2, but got %d )", centerPoints.length ) ); + if (centerPoints.length != 2) { // Fixed: was == now != + throw new IllegalArgumentException(String.format("Invalid number of values for \"center\" (expected 2, but got %d)", centerPoints.length)); } - if ( validateExpression( alpha, alphaParamLegalExpressions ) == Boolean.FALSE ) { - throw new IllegalArgumentException( - String.format( "Invalid value for \"alpha\" (expected format %d )", alphaParamLegalExpressions.toString() ) ); + if (validateExpression(alpha, alphaParamLegalExpressions) == Boolean.FALSE) { + throw new IllegalArgumentException(String.format("Invalid value for \"alpha\" (expected one of %s)", java.util.Arrays.toString(alphaParamLegalExpressions))); } - if ( validateExpression( alpha, faceParamLegalExpressions ) == Boolean.FALSE ) { - throw new IllegalArgumentException( - String.format( "Invalid value for \"face\" (expected format %d )", faceParamLegalExpressions.toString() ) ); + if (validateExpression(face, faceParamLegalExpressions) == Boolean.FALSE) { + throw new IllegalArgumentException(String.format("Invalid value for \"face\" (expected one of %s)", java.util.Arrays.toString(faceParamLegalExpressions))); } } - public static Boolean validateExpression( String expression, String[] legalExpressions ) { + public static Boolean validateExpression(String expression, String[] legalExpressions) { Boolean validInput = Boolean.FALSE; - for ( String legalExpression : legalExpressions ) { - if ( legalExpression.equals( expression ) ) { + for (String legalExpression : legalExpressions) { + if (legalExpression.equals(expression)) { validInput = Boolean.TRUE; break; } @@ -53,59 +96,47 @@ public static Boolean validateExpression( String expression, String[] legalExpre return validInput; } - public static Boolean matrixUnion( RealMatrix firstMatrix, RealMatrix secondMatrix ) { - double[][] finalMatrixRowModel = new double[ firstMatrix.getRowDimension() + secondMatrix.getRowDimension() ][]; - for ( int i = 0; i < firstMatrix.getRowDimension(); i++ ) { - finalMatrixRowModel[ i ] = firstMatrix.getRow( i ); - } - for ( int i = 0; i < secondMatrix.getRowDimension(); i++ ) { - finalMatrixRowModel[ i ] = secondMatrix.getRow( i ); - } - RealMatrix finalMatrix = new Array2DRowRealMatrix(finalMatrixRowModel); - GenericDOEUtil.matrixLogger( finalMatrix, "matrixUnion" ); - } - } //<-----------------------------STEPS----------------------------------> /** * Detailed Steps to Port Central Composite Design Script to Java - * + *

* $ 1. Create CompositeDOE Class Structure Create CompositeDOE.java in com.jdoe.algorithms package Add imports: * org.apache.commons.math3.linear.* for matrix operations Define method signature: public static RealMatrix centralCompositeDesign(int n, * int[] center, String alpha, String face) Set default parameters handling for center, alpha, and face - * + *

* $ 2. Implement Input Validation Logic Add assertion: assert n > 1 : "\"n\" must be an integer greater than 1" Validate alpha parameter * with: alpha.toLowerCase() and check against allowed values Validate face parameter with: face.toLowerCase() and check against allowed * values Validate center array length: if (center.length != 2) throw IllegalArgumentException Format error messages to match Python * implementation - * + *

* $ 3. Implement Star Point Generation Create StarDOE.java class with star method Calculate alpha based on alpha type: For orthogonal: * alpha = Math.pow(2 * factorial(n), 0.25) where factorial(n) = n! For rotatable: alpha = Math.pow(2, 0.5) for 2D, Math.pow(3, 0.5) for 3D, * etc. Generate 2n star points with values [±alpha, 0, ..., 0], [0, ±alpha, ..., 0], etc. Return both star matrix and alpha value as an * Object array - * + *

* 4. Enhance FactorialDOE Class Add ff2n method that generates 2-level factorial design Create 2^n × n matrix with values -1 and +1 Use bit * manipulation: for each row i and column j, set value to ((i >> j) & 1) == 0 ? -1 : 1 Return RealMatrix object - * + *

* 5. Create UnionDOE Class Implement union method that combines two matrices vertically Use Array2DRowRealMatrix constructor to create new * matrix with combined rows Copy data from both input matrices to the result matrix Return the concatenated matrix - * + *

* 6. Create RepeatCenterDOE Class Implement repeatCenter method that generates center points Create repeats × n matrix filled with zeros * Use new Array2DRowRealMatrix(repeats, n) to initialize matrix Return the center point matrix - * + *

* 7. Implement Core Algorithm Logic Initialize H1 andH2 matrices based on face type For inscribed (cci): scale factorial points H1 = * H1.scalarMultiply(1.0/a) For faced (ccf): set alpha to 1.0 For circumscribed (ccc): keep original factorial points Generate center points * C1 and C2 using repeatCenter Combine matrices: H1 = union(H1, C1), H2 = union(H2, C2), H = union(H1, H2) - * + *

* 8. Handle Matrix Operations Implement matrix scaling using scalarMultiply() method from Commons Math Ensure all matrix dimensions are * compatible for union operations Use getSubMatrix() and setSubMatrix() for matrix manipulations Return final RealMatrix object - * + *

* 9. Create Unit Tests Create CompositeDOETest.java with comprehensive test cases Test all three face types: ccc, cci, ccf Test both alpha * types: orthogonal, rotatable Validate matrix dimensions: (2^n + 2n + sum(center)) × n Verify center point counts and star point * distances - * + *

* 10. Add Documentation and Error Handling Add JavaDoc explaining parameters, return value, and exceptions Implement proper exception * handling for invalid inputs Include example usage in documentation Follow the same error message format as Python implementation */ \ No newline at end of file diff --git a/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java b/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java new file mode 100644 index 0000000..15e8244 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java @@ -0,0 +1,12 @@ +package com.jdoe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +public class RepeatCenterDOE { + + public static RealMatrix repeatCenter(int numberOfFactors, int repeats) { + return new Array2DRowRealMatrix(repeats, numberOfFactors); + } + +} diff --git a/src/main/java/com/jdoe/algorithms/StarDOE.java b/src/main/java/com/jdoe/algorithms/StarDOE.java new file mode 100644 index 0000000..f150d01 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/StarDOE.java @@ -0,0 +1,63 @@ +package com.jdoe.algorithms; + +import com.jdoe.util.GenericDOEUtil; + +public class StarDOE { + /** + * Create the star points of various design matrices + * + * @param n The number of variables in the design + * @param alpha Available values are 'faced' (default), 'orthogonal', or 'rotatable' + * @param center A 1-by-2 array of integers indicating the number of center points + * assigned in each block of the response surface design. Default is (1, 1). + * @return A two-element array containing the star-point portion of the design matrix + * and the alpha value to scale the star points with. + */ + public static Object[] star(int n, String alpha, int[] center) { + double a; + // Star points at the center of each face of the factorial + if ("faced".equals(alpha)) { + a = 1.0; + } else if ("orthogonal".equals(alpha)) { + int nc = (int) Math.pow(2, n); // factorial points + int nco = center[0]; // center points to factorial + int na = 2 * n; // axial points + int nao = center[1]; // center points to axial design + // value of alpha in orthogonal design + a = Math.sqrt(n * (1 + (double) nao / na) / (1 + (double) nco / nc)); + } else if ("rotatable".equals(alpha)) { + int nc = (int) Math.pow(2, n); // number of factorial points + a = Math.pow(nc, 0.25); // value of alpha in rotatable design + } else { + throw new IllegalArgumentException("Invalid value for \"alpha\": " + alpha); + } + // Create the actual matrix now. + double[][] H = new double[2 * n][n]; + for (int i = 0; i < n; i++) { + H[2 * i][i] = -1.0; + H[2 * i + 1][i] = 1.0; + } + // Multiply matrix by alpha + for (int i = 0; i < H.length; i++) { + for (int j = 0; j < H[i].length; j++) { + H[i][j] *= a; + } + } + return new Object[]{H, a}; + } + + /** + * Overloaded method with default center value (1, 1) + */ + public static Object[] star(int n, String alpha) { + return star(n, alpha, new int[]{1, 1}); + } + + /** + * Overloaded method with default alpha value "faced" and default center value (1, 1) + */ + public static Object[] star(int n) { + return star(n, "faced", new int[]{1, 1}); + } + +} diff --git a/src/main/java/com/jdoe/algorithms/UnionDOE.java b/src/main/java/com/jdoe/algorithms/UnionDOE.java new file mode 100644 index 0000000..576edc0 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/UnionDOE.java @@ -0,0 +1,33 @@ +package com.jdoe.algorithms; + +import com.jdoe.util.GenericDOEUtil; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +public class UnionDOE { + public static RealMatrix matrixUnion(RealMatrix firstMatrix, RealMatrix secondMatrix) { + int totalRows = firstMatrix.getRowDimension() + secondMatrix.getRowDimension(); + int cols = firstMatrix.getColumnDimension(); + + // Verify both matrices have the same number of columns + if (cols != secondMatrix.getColumnDimension()) { + throw new IllegalArgumentException("Both matrices must have the same number of columns"); + } + + double[][] finalMatrixRowModel = new double[totalRows][cols]; + + // Copy rows from the first matrix + for (int i = 0; i < firstMatrix.getRowDimension(); i++) { + finalMatrixRowModel[i] = firstMatrix.getRow(i); + } + + // Copy rows from the second matrix starting after the first matrix's rows + for (int i = 0; i < secondMatrix.getRowDimension(); i++) { + finalMatrixRowModel[firstMatrix.getRowDimension() + i] = secondMatrix.getRow(i); + } + + RealMatrix finalMatrix = new Array2DRowRealMatrix(finalMatrixRowModel); + GenericDOEUtil.matrixLogger(finalMatrix, "matrixUnion"); + return finalMatrix; + } +} diff --git a/src/test/java/com/jdoe/algorithms/CompositeDOETest.java b/src/test/java/com/jdoe/algorithms/CompositeDOETest.java new file mode 100644 index 0000000..89959a9 --- /dev/null +++ b/src/test/java/com/jdoe/algorithms/CompositeDOETest.java @@ -0,0 +1,195 @@ +package com.jdoe.algorithms; + +import org.apache.commons.math3.linear.RealMatrix; +import org.junit.Test; +import static org.junit.Assert.*; + +public class CompositeDOETest { + + @Test + public void testCentralCompositeDesign_OrthogonalCircumscribed_ValidOutput() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "circumscribed"; + + // Act + RealMatrix result = CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert + assertNotNull(result); + // For 2 factors: 2^2 factorial + 2*2 star points + 1+1 center points = 4+4+2 = 10 rows, 2 columns + assertEquals(10, result.getRowDimension()); + assertEquals(2, result.getColumnDimension()); + } + + @Test + public void testCentralCompositeDesign_RotatableCircumscribed_ValidOutput() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "rotatable"; + String face = "circumscribed"; + + // Act + RealMatrix result = CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert + assertNotNull(result); + // For 2 factors: 2^2 factorial + 2*2 star points + 1+1 center points = 4+4+2 = 10 rows, 2 columns + assertEquals(10, result.getRowDimension()); + assertEquals(2, result.getColumnDimension()); + } + + @Test + public void testCentralCompositeDesign_OrthogonalInscribed_ValidOutput() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "inscribed"; + + // Act + RealMatrix result = CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert + assertNotNull(result); + // For 2 factors: 2^2 factorial + 2*2 star points + 1+1 center points = 4+4+2 = 10 rows, 2 columns + assertEquals(10, result.getRowDimension()); + assertEquals(2, result.getColumnDimension()); + } + + @Test + public void testCentralCompositeDesign_FacedCCF_ValidOutput() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "faced"; + + // Act + RealMatrix result = CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert + assertNotNull(result); + // For 2 factors: 2^2 factorial + 2*2 star points + 1+1 center points = 4+4+2 = 10 rows, 2 columns + assertEquals(10, result.getRowDimension()); + assertEquals(2, result.getColumnDimension()); + } + + @Test + public void testCentralCompositeDesign_ThreeFactors_ValidOutput() { + // Arrange + int numberOfFactors = 3; + int[] centerPoints = {2, 2}; + String alpha = "rotatable"; + String face = "circumscribed"; + + // Act + RealMatrix result = CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert + assertNotNull(result); + // For 3 factors: 2^3 factorial + 2*3 star points + 2+2 center points = 8+6+4 = 18 rows, 3 columns + assertEquals(18, result.getRowDimension()); + assertEquals(3, result.getColumnDimension()); + } + + @Test(expected = IllegalArgumentException.class) + public void testCentralCompositeDesign_InvalidAlpha_ThrowsException() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "invalid"; + String face = "circumscribed"; + + // Act + CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert: Exception expected + } + + @Test(expected = IllegalArgumentException.class) + public void testCentralCompositeDesign_InvalidFace_ThrowsException() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "invalid"; + + // Act + CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert: Exception expected + } + + @Test(expected = IllegalArgumentException.class) + public void testCentralCompositeDesign_InvalidNumberOfFactors_ThrowsException() { + // Arrange + int numberOfFactors = 0; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "circumscribed"; + + // Act + CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert: Exception expected + } + + @Test(expected = IllegalArgumentException.class) + public void testCentralCompositeDesign_InvalidCenterPointsLength_ThrowsException() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1}; // Only 1 element instead of required 2 + String alpha = "orthogonal"; + String face = "circumscribed"; + + // Act + CompositeDOE.centralCompositeDesign(numberOfFactors, centerPoints, alpha, face); + + // Assert: Exception expected + } + + @Test + public void testValidateCompositeDesignParameters_ValidInput_NoException() { + // Arrange + int numberOfFactors = 2; + int[] centerPoints = {1, 1}; + String alpha = "orthogonal"; + String face = "circumscribed"; + + // Act & Assert + // Should not throw any exception + CompositeDOE.validateCompositeDegisnParameters(numberOfFactors, centerPoints, alpha, face); + } + + @Test + public void testValidateExpression_ValidExpression_ReturnsTrue() { + // Arrange + String expression = "orthogonal"; + String[] legalExpressions = {"orthogonal", "o", "rotatable", "r"}; + + // Act + Boolean result = CompositeDOE.validateExpression(expression, legalExpressions); + + // Assert + assertTrue(result); + } + + @Test + public void testValidateExpression_InvalidExpression_ReturnsFalse() { + // Arrange + String expression = "invalid"; + String[] legalExpressions = {"orthogonal", "o", "rotatable", "r"}; + + // Act + Boolean result = CompositeDOE.validateExpression(expression, legalExpressions); + + // Assert + assertFalse(result); + } + + +} diff --git a/target/classes/log4j2.xml b/target/classes/log4j2.xml new file mode 100644 index 0000000..4cbeecc --- /dev/null +++ b/target/classes/log4j2.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst b/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst new file mode 100644 index 0000000..c9315c0 --- /dev/null +++ b/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst @@ -0,0 +1,6 @@ +com/jdoe/util/FactorialUtility.class +com/jdoe/util/BuildRegressionMatrixUtility.class +com/jdoe/algorithms/BuildRegressionMatrix.class +com/jdoe/algorithms/BoxBehnkenDOE.class +com/jdoe/algorithms/FactorialDOE.class +com/jdoe/Testing.class diff --git a/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst b/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst new file mode 100644 index 0000000..b2926fa --- /dev/null +++ b/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst @@ -0,0 +1,6 @@ +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/Testing.java +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/BuildRegressionMatrix.java +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/FactorialDOE.java +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java +/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/util/FactorialUtility.java diff --git a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst new file mode 100644 index 0000000..8757472 --- /dev/null +++ b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst @@ -0,0 +1,2 @@ +com/jdoe/algorithms/FactorialDOETest.class +com/jdoe/algorithms/BoxBehnkenDOETest.class diff --git a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst new file mode 100644 index 0000000..91aac98 --- /dev/null +++ b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst @@ -0,0 +1,2 @@ +/home/noor/IdeaProjects/JDOE/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java +/home/noor/IdeaProjects/JDOE/src/test/java/com/jdoe/algorithms/FactorialDOETest.java diff --git a/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml b/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml new file mode 100644 index 0000000..8c02771 --- /dev/null +++ b/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt b/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt new file mode 100644 index 0000000..2aad563 --- /dev/null +++ b/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt @@ -0,0 +1,4 @@ +------------------------------------------------------------------------------- +Test set: com.jdoe.algorithms.FactorialDOETest +------------------------------------------------------------------------------- +Tests run: 18, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.233 s -- in com.jdoe.algorithms.FactorialDOETest From 309b8b216c405cbf00cf70b83d09cb25fd103e93 Mon Sep 17 00:00:00 2001 From: nomus Date: Tue, 27 Jan 2026 09:42:56 +0500 Subject: [PATCH 4/5] experiments --- ...rix.java => BuildRegressionMatrixDOE.java} | 2 +- .../algorithms/CranleyPattersonShiftDOE.java | 66 ++++++++++ .../java/com/jdoe/algorithms/DoeFold.java | 95 ++++++++++++++ .../java/com/jdoe/algorithms/DoehlertDOE.java | 117 ++++++++++++++++++ 4 files changed, 279 insertions(+), 1 deletion(-) rename src/main/java/com/jdoe/algorithms/{BuildRegressionMatrix.java => BuildRegressionMatrixDOE.java} (99%) create mode 100644 src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeFold.java create mode 100644 src/main/java/com/jdoe/algorithms/DoehlertDOE.java diff --git a/src/main/java/com/jdoe/algorithms/BuildRegressionMatrix.java b/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java similarity index 99% rename from src/main/java/com/jdoe/algorithms/BuildRegressionMatrix.java rename to src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java index 393b7d7..daaa2cf 100644 --- a/src/main/java/com/jdoe/algorithms/BuildRegressionMatrix.java +++ b/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java @@ -6,7 +6,7 @@ import com.jdoe.util.BuildRegressionMatrixUtility; -public class BuildRegressionMatrix { +public class BuildRegressionMatrixDOE { /** * Builds a regression matrix from an experimental design matrix and a mathematical model string. diff --git a/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java new file mode 100644 index 0000000..0add030 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java @@ -0,0 +1,66 @@ +package com.jdoe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.RandomDataGenerator; +import org.apache.commons.math3.random.JDKRandomGenerator; + +public class CranleyPattersonShiftDOE { + + public static RealMatrix cranleyPattersonShift(double[][] points, Integer seed) { + if (points == null || points.length == 0 || points[0].length == 0) { + throw new IllegalArgumentException("Input `points` must be a non-null 2D array with at least one element."); + } + + RealMatrix pointsMatrix = new Array2DRowRealMatrix(points); + + int nSamples = pointsMatrix.getRowDimension(); + int dim = pointsMatrix.getColumnDimension(); + + // Initialize random generator with seed if provided + RandomDataGenerator rng; + if (seed != null) { + JDKRandomGenerator jdkRng = new JDKRandomGenerator(); + jdkRng.setSeed(seed); + rng = new RandomDataGenerator(jdkRng); // Use Apache Commons Math generator + } else { + rng = new RandomDataGenerator(); + } + + // Generate random shift vector + double[] shiftVector = new double[dim]; + for (int i = 0; i < dim; i++) { + shiftVector[i] = rng.nextUniform(0.0, 1.0); + } + + // Create a matrix with the shift vector replicated for each sample + double[][] shiftMatrixData = new double[nSamples][dim]; + for (int i = 0; i < nSamples; i++) { + for (int j = 0; j < dim; j++) { + shiftMatrixData[i][j] = shiftVector[j]; + } + } + + RealMatrix shiftMatrix = new Array2DRowRealMatrix(shiftMatrixData); + + // Add the shift vector to all points + RealMatrix shiftedPoints = pointsMatrix.add(shiftMatrix); + + // Wrap values to [0, 1) using modulo operation + double[][] resultData = shiftedPoints.getData(); + for (int i = 0; i < nSamples; i++) { + for (int j = 0; j < dim; j++) { + resultData[i][j] = resultData[i][j] - Math.floor(resultData[i][j]); + } + } + + return new Array2DRowRealMatrix(resultData); + } + + /** + * Overloaded method without seed parameter + */ + public static RealMatrix cranleyPattersonShift(double[][] points) { + return cranleyPattersonShift(points, null); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeFold.java b/src/main/java/com/jdoe/algorithms/DoeFold.java new file mode 100644 index 0000000..e9cda7d --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeFold.java @@ -0,0 +1,95 @@ +package com.jdoe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +public class DoeFold { + + /** + * Fold a design to reduce confounding effects. + * + * @param H The design matrix to be folded + * @param columns Indices of columns to fold (Default: null). If columns=null is + * used, then all columns will be folded. + * @return The folded design matrix + */ + public static RealMatrix fold(double[][] H, int[] columns) { + if (H.length == 0 || H[0].length == 0 || H.length != H[0].length) { + throw new IllegalArgumentException("Input design matrix must be 2d."); + } + + int nCols = H[0].length; + + if (columns == null) { + columns = new int[nCols]; + for (int i = 0; i < nCols; i++) { + columns[i] = i; + } + } + + double[][] Hf = new double[H.length][H[0].length]; + + // Copy original matrix to Hf + for (int i = 0; i < H.length; i++) { + System.arraycopy(H[i], 0, Hf[i], 0, H[0].length); + } + + for (int col : columns) { + // Get unique values in column + double[] vals = getUniqueValues(getColumn(H, col)); + + if (vals.length != 2) { + throw new IllegalArgumentException("Input design matrix must be 2-level factors only."); + } + + for (int i = 0; i < H.length; i++) { + Hf[i][col] = (H[i][col] == vals[1]) ? vals[0] : vals[1]; + } + } + + // Combine original and folded matrices + double[][] result = new double[H.length * 2][H[0].length]; + + // Copy original matrix + for (int i = 0; i < H.length; i++) { + System.arraycopy(H[i], 0, result[i], 0, H[0].length); + } + + // Copy folded matrix + for (int i = 0; i < Hf.length; i++) { + System.arraycopy(Hf[i], 0, result[H.length + i], 0, Hf[0].length); + } + + return new Array2DRowRealMatrix(result); + } + + /** + * Overloaded method with default columns (all columns) + */ + public static RealMatrix fold(double[][] H) { + return fold(H, null); + } + + /** + * Helper method to extract a column from a 2D array + */ + private static double[] getColumn(double[][] matrix, int colIndex) { + double[] column = new double[matrix.length]; + for (int i = 0; i < matrix.length; i++) { + column[i] = matrix[i][colIndex]; + } + return column; + } + + /** + * Helper method to get unique values from an array + */ + private static double[] getUniqueValues(double[] arr) { + java.util.Set uniqueSet = new java.util.LinkedHashSet<>(); + for (double val : arr) { + uniqueSet.add(val); + } + return uniqueSet.stream().mapToDouble(Double::doubleValue).toArray(); + } +} + diff --git a/src/main/java/com/jdoe/algorithms/DoehlertDOE.java b/src/main/java/com/jdoe/algorithms/DoehlertDOE.java new file mode 100644 index 0000000..5b0488b --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoehlertDOE.java @@ -0,0 +1,117 @@ +package com.jdoe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import java.util.ArrayList; +import java.util.List; + +public class DoehlertDOE { + + /** + * Generate a Doehlert design matrix for a given number of factors using a shell approach. + * + * @param numFactors Number of factors (variables) + * @param numCenterPoints Number of center (replicate) points + * @return Design matrix of shape (N, numFactors), where N = k^2 + k + C + */ + public static RealMatrix doehlertShellDesign(int numFactors, int numCenterPoints) { + if (numFactors < 1) { + throw new IllegalArgumentException("Number of factors must be at least 1."); + } + + List designPoints = new ArrayList<>(); + + // Add center points + for (int i = 0; i < numCenterPoints; i++) { + double[] centerPoint = new double[numFactors]; + for (int j = 0; j < numFactors; j++) { + centerPoint[j] = 0.0; + } + designPoints.add(centerPoint); + } + + // Add shells progressively + for (int shellIndex = 1; shellIndex <= numFactors; shellIndex++) { + int numShellPoints = shellIndex + 1; + double radius = Math.sqrt(shellIndex / (double) numFactors); + + for (int pointIdx = 0; pointIdx < numShellPoints; pointIdx++) { + double[] point = new double[numFactors]; + double angle = (2.0 * Math.PI * pointIdx) / numShellPoints; + + if (numFactors > 0) point[0] = radius * Math.cos(angle); + if (numFactors > 1) point[1] = radius * Math.sin(angle); + + designPoints.add(point); + } + } + + // Convert list of points to matrix + double[][] result = new double[designPoints.size()][numFactors]; + for (int i = 0; i < designPoints.size(); i++) { + result[i] = designPoints.get(i); + } + + return new Array2DRowRealMatrix(result); + } + + /** + * Overloaded method with default center points (1) + */ + public static RealMatrix doehlertShellDesign(int numFactors) { + return doehlertShellDesign(numFactors, 1); + } + + /** + * Generate a Doehlert design matrix using a simplex-based approach. + * + * @param numFactors Number of factors included in the design + * @return Design matrix with experiments at different coded levels + */ + public static RealMatrix doehlertSimplexDesign(int numFactors) { + double[][] simplexMatrix = new double[numFactors + 1][numFactors]; + + for (int i = 1; i <= numFactors; i++) { + for (int j = 1; j <= i; j++) { + if (j == i) { + simplexMatrix[i][j - 1] = Math.sqrt(i + 1) / Math.sqrt(2 * i); + } else { + int denom = 2 * (i - (j - 1)) * (i - j); + simplexMatrix[i][i - j] = 1 / Math.sqrt(denom); + } + } + } + + List extraPoints = new ArrayList<>(); + for (int i = 0; i <= numFactors; i++) { + List rangeList = new ArrayList<>(); + for (int k = 1; k < i; k++) rangeList.add(k); + for (int k = i + 1; k <= numFactors; k++) rangeList.add(k); + + for (int j : rangeList) { + double[] point = new double[numFactors]; + for (int k = 0; k < numFactors; k++) { + point[k] = simplexMatrix[i][k] - simplexMatrix[j][k]; + } + extraPoints.add(point); + } + } + + // Combine matrices + int totalRows = (numFactors + 1) + extraPoints.size(); + double[][] fullMatrix = new double[totalRows][numFactors]; + + // Copy simplex matrix + for (int i = 0; i < simplexMatrix.length; i++) { + fullMatrix[i] = simplexMatrix[i]; + } + + // Add extra points + for (int i = 0; i < extraPoints.size(); i++) { + fullMatrix[simplexMatrix.length + i] = extraPoints.get(i); + } + + return new Array2DRowRealMatrix(fullMatrix); + } + +} From b6eff74aeb6b23623cb2bc0839627db3048a60d6 Mon Sep 17 00:00:00 2001 From: nomus Date: Sun, 1 Feb 2026 15:49:25 +0500 Subject: [PATCH 5/5] Add remaining experiments --- .gitignore | 3 + README.md | 2 +- dependency-reduced-pom.xml | 172 +++--- pom.xml | 6 +- src/main/java/com/jdoe/Testing.java | 6 +- .../com/jdoe/algorithms/BoxBehnkenDOE.java | 2 +- .../algorithms/BuildRegressionMatrixDOE.java | 4 +- .../com/jdoe/algorithms/CompositeDOE.java | 6 +- .../algorithms/CranleyPattersonShiftDOE.java | 2 +- .../java/com/jdoe/algorithms/DoeFold.java | 2 +- src/main/java/com/jdoe/algorithms/DoeGsd.java | 287 +++++++++ .../java/com/jdoe/algorithms/DoeHalton.java | 128 ++++ .../java/com/jdoe/algorithms/DoeKorobov.java | 85 +++ src/main/java/com/jdoe/algorithms/DoeLhs.java | 551 ++++++++++++++++++ .../jdoe/algorithms/DoePlackettBurman.java | 251 ++++++++ .../com/jdoe/algorithms/DoeRandomKMeans.java | 134 +++++ .../com/jdoe/algorithms/DoeRandomUniform.java | 50 ++ .../java/com/jdoe/algorithms/DoeRank1.java | 84 +++ .../java/com/jdoe/algorithms/DoeSaltelli.java | 157 +++++ .../java/com/jdoe/algorithms/DoeSobol.java | 129 ++++ .../com/jdoe/algorithms/DoeSparseGrid.java | 313 ++++++++++ .../java/com/jdoe/algorithms/DoeTaguchi.java | 182 ++++++ .../java/com/jdoe/algorithms/DoeUnion.java | 58 ++ .../com/jdoe/algorithms/DoeVanillaMorris.java | 210 +++++++ .../java/com/jdoe/algorithms/DoehlertDOE.java | 2 +- .../com/jdoe/algorithms/FactorialDOE.java | 4 +- .../com/jdoe/algorithms/RepeatCenterDOE.java | 2 +- .../java/com/jdoe/algorithms/StarDOE.java | 4 +- .../java/com/jdoe/algorithms/UnionDOE.java | 4 +- .../util/BuildRegressionMatrixUtility.java | 2 +- .../java/com/jdoe/util/FactorialUtility.java | 2 +- .../java/com/jdoe/util/GenericDOEUtil.java | 2 +- .../jdoe/algorithms/BoxBehnkenDOETest.java | 2 +- .../com/jdoe/algorithms/CompositeDOETest.java | 2 +- .../com/jdoe/algorithms/FactorialDOETest.java | 2 +- target/classes/log4j2.xml | 13 - .../compile/default-compile/createdFiles.lst | 6 - .../compile/default-compile/inputFiles.lst | 6 - .../default-testCompile/createdFiles.lst | 2 - .../default-testCompile/inputFiles.lst | 2 - ...T-com.jdoe.algorithms.FactorialDOETest.xml | 81 --- .../com.jdoe.algorithms.FactorialDOETest.txt | 4 - 42 files changed, 2737 insertions(+), 229 deletions(-) create mode 100644 src/main/java/com/jdoe/algorithms/DoeGsd.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeHalton.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeKorobov.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeLhs.java create mode 100644 src/main/java/com/jdoe/algorithms/DoePlackettBurman.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeRandomKMeans.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeRandomUniform.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeRank1.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeSaltelli.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeSobol.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeSparseGrid.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeTaguchi.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeUnion.java create mode 100644 src/main/java/com/jdoe/algorithms/DoeVanillaMorris.java delete mode 100644 target/classes/log4j2.xml delete mode 100644 target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst delete mode 100644 target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst delete mode 100644 target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst delete mode 100644 target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst delete mode 100644 target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml delete mode 100644 target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt diff --git a/.gitignore b/.gitignore index 8f7b995..b5671e5 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ *.zip *.tar.gz *.rar +*.lst # virtual machine crash logs, see http://www.java.com/en/download/help/error_hotspot.xml hs_err_pid* @@ -25,3 +26,5 @@ replay_pid* # intellij .idea + +target/ diff --git a/README.md b/README.md index 2b5589a..360bffd 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,2 @@ -# JDOE +# doe Java Design of Experiments diff --git a/dependency-reduced-pom.xml b/dependency-reduced-pom.xml index f1be5e4..d122e5b 100644 --- a/dependency-reduced-pom.xml +++ b/dependency-reduced-pom.xml @@ -1,86 +1,86 @@ - - - 4.0.0 - com.jdoe - doe-genereter - JDOE - 1.0-SNAPSHOT - http://maven.apache.org - - - - maven-jar-plugin - 3.4.1 - - - - com.jdoe.Testing - - - - - - maven-shade-plugin - 3.6.1 - - - package - - shade - - - - - - - - - - - - - - junit - junit - 4.13.2 - test - - - hamcrest-core - org.hamcrest - - - - - org.mockito - mockito-core - 5.20.0 - test - - - byte-buddy - net.bytebuddy - - - byte-buddy-agent - net.bytebuddy - - - objenesis - org.objenesis - - - - - org.mockito - mockito-inline - 5.2.0 - test - - - - 17 - 17 - UTF-8 - - + + + 4.0.0 + com.doe + doe-genereter + doe + 1.0-SNAPSHOT + http://maven.apache.org + + + + maven-jar-plugin + 3.4.1 + + + + com.doe.Testing + + + + + + maven-shade-plugin + 3.6.1 + + + package + + shade + + + + + + + + + + + + + + junit + junit + 4.13.2 + test + + + hamcrest-core + org.hamcrest + + + + + org.mockito + mockito-core + 5.20.0 + test + + + byte-buddy + net.bytebuddy + + + byte-buddy-agent + net.bytebuddy + + + objenesis + org.objenesis + + + + + org.mockito + mockito-inline + 5.2.0 + test + + + + 17 + 17 + UTF-8 + + diff --git a/pom.xml b/pom.xml index 5934c85..63d96fb 100644 --- a/pom.xml +++ b/pom.xml @@ -2,12 +2,12 @@ xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd"> 4.0.0 - com.jdoe + com.doe doe-genereter 1.0-SNAPSHOT jar - JDOE + doe http://maven.apache.org @@ -79,7 +79,7 @@ - com.jdoe.Testing + com.doe.Testing diff --git a/src/main/java/com/jdoe/Testing.java b/src/main/java/com/jdoe/Testing.java index 75d54d7..b64425b 100644 --- a/src/main/java/com/jdoe/Testing.java +++ b/src/main/java/com/jdoe/Testing.java @@ -1,7 +1,7 @@ -package com.jdoe; +package com.doe; -import com.jdoe.algorithms.BoxBehnkenDOE; -import com.jdoe.algorithms.FactorialDOE; +import com.doe.algorithms.BoxBehnkenDOE; +import com.doe.algorithms.FactorialDOE; /** * Hello world! diff --git a/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java b/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java index 1a2cf2d..ca910f8 100644 --- a/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java +++ b/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java b/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java index daaa2cf..bd56aec 100644 --- a/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java +++ b/src/main/java/com/jdoe/algorithms/BuildRegressionMatrixDOE.java @@ -1,10 +1,10 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import java.util.ArrayList; import java.util.Arrays; import java.util.List; -import com.jdoe.util.BuildRegressionMatrixUtility; +import com.doe.util.BuildRegressionMatrixUtility; public class BuildRegressionMatrixDOE { diff --git a/src/main/java/com/jdoe/algorithms/CompositeDOE.java b/src/main/java/com/jdoe/algorithms/CompositeDOE.java index f4051be..922d3e3 100644 --- a/src/main/java/com/jdoe/algorithms/CompositeDOE.java +++ b/src/main/java/com/jdoe/algorithms/CompositeDOE.java @@ -1,11 +1,11 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowFieldMatrix; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; import org.jetbrains.annotations.NotNull; -import com.jdoe.util.GenericDOEUtil; +import com.doe.util.GenericDOEUtil; public class CompositeDOE { @@ -103,7 +103,7 @@ public static Boolean validateExpression(String expression, String[] legalExpres /** * Detailed Steps to Port Central Composite Design Script to Java *

- * $ 1. Create CompositeDOE Class Structure Create CompositeDOE.java in com.jdoe.algorithms package Add imports: + * $ 1. Create CompositeDOE Class Structure Create CompositeDOE.java in com.doe.algorithms package Add imports: * org.apache.commons.math3.linear.* for matrix operations Define method signature: public static RealMatrix centralCompositeDesign(int n, * int[] center, String alpha, String face) Set default parameters handling for center, alpha, and face *

diff --git a/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java index 0add030..4b11dac 100644 --- a/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java +++ b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/algorithms/DoeFold.java b/src/main/java/com/jdoe/algorithms/DoeFold.java index e9cda7d..34e3b09 100644 --- a/src/main/java/com/jdoe/algorithms/DoeFold.java +++ b/src/main/java/com/jdoe/algorithms/DoeFold.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/algorithms/DoeGsd.java b/src/main/java/com/jdoe/algorithms/DoeGsd.java new file mode 100644 index 0000000..2daf4f9 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeGsd.java @@ -0,0 +1,287 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.stream.IntStream; + +public class DoeGsd { + + /** + * Create a Generalized Subset Design (GSD). + * + * @param factorLevels Number of factor levels per factor in design. + * @param reductionFactor Reduction factor (bigger than 1). Larger reductionFactor means fewer + * experiments in the design and more possible complementary designs. + * @param numberOfComplementaryDesigns Number of complementary GSD-designs (default 1). The complementary + * designs are balanced analogous to fold-over in two-level fractional + * factorial designs. + * @return n m-by-k matrices where k is the number of factors (equal + * to the length of factorLevels. The number of rows, m, will + * be approximately equal to the grand product of the factor levels + * divided by reductionFactor. + */ + public static RealMatrix[] gsd(int[] factorLevels, int reductionFactor, int numberOfComplementaryDesigns) { + // Input validation + for (int level : factorLevels) { + if (level <= 0) { + throw new IllegalArgumentException("factorLevels has to be sequence of positive integers"); + } + } + if (reductionFactor <= 1) { + throw new IllegalArgumentException("reductionFactor has to be integer larger than 1"); + } + if (numberOfComplementaryDesigns <= 0) { + throw new IllegalArgumentException("numberOfComplementaryDesigns has to be positive integer"); + } + + List>> partitions = _makePartitions(factorLevels, reductionFactor); + int[][] latinSquare = _makeLatinSquare(reductionFactor); + RealMatrix[] orthogonalArrays = _makeOrthogonalArrays(latinSquare, factorLevels.length); + + RealMatrix[] designs; + try { + designs = new RealMatrix[orthogonalArrays.length]; + for (int i = 0; i < orthogonalArrays.length; i++) { + RealMatrix mappedDesign = _mapPartitionsToDesign(partitions, orthogonalArrays[i]); + // Subtract 1 from each element (equivalent to Python's -1) + double[][] designData = mappedDesign.getData(); + for (double[] row : designData) { + for (int col = 0; col < row.length; col++) { + row[col] -= 1; + } + } + designs[i] = new Array2DRowRealMatrix(designData); + } + } catch (Exception e) { + throw new IllegalArgumentException("reductionFactor too large compared to factor levels"); + } + + if (numberOfComplementaryDesigns == 1) { + return new RealMatrix[]{designs[0]}; + } else { + RealMatrix[] result = new RealMatrix[numberOfComplementaryDesigns]; + System.arraycopy(designs, 0, result, 0, numberOfComplementaryDesigns); + return result; + } + } + + /** + * Overloaded method with default number of complementary designs (1) + */ + public static RealMatrix gsd(int[] factorLevels, int reductionFactor) { + RealMatrix[] result = gsd(factorLevels, reductionFactor, 1); + return result[0]; + } + + /** + * Augment latin-square to the specified number of columns to produce + * an orthogonal array. + */ + private static RealMatrix[] _makeOrthogonalArrays(int[][] latinSquare, int numberOfColumns) { + int[] firstRow = latinSquare[0]; + RealMatrix[] matrixArray = new RealMatrix[firstRow.length]; + for (int i = 0; i < firstRow.length; i++) { + double[][] singleValueMatrix = {{firstRow[i]}}; + matrixArray[i] = new Array2DRowRealMatrix(singleValueMatrix); + } + + while (matrixArray[0].getColumnDimension() < numberOfColumns) { + RealMatrix[] newMatrixArray = new RealMatrix[matrixArray.length]; + + for (int i = 0; i < matrixArray.length; i++) { + RealMatrix currentMatrix = matrixArray[i]; + List subMatrices = new ArrayList<>(); + + for (int j = 0; j < firstRow.length; j++) { + int constantValue = firstRow[j]; + RealMatrix otherMatrix = matrixArray[latinSquare[i][j]]; + + // Create constant vector + int otherRowCount = otherMatrix.getRowDimension(); + double[][] constantVector = new double[otherRowCount][1]; + for (int row = 0; row < otherRowCount; row++) { + constantVector[row][0] = constantValue; + } + + // Horizontally stack constant vector with other matrix + double[][] combinedData = new double[otherRowCount][1 + otherMatrix.getColumnDimension()]; + for (int row = 0; row < otherRowCount; row++) { + combinedData[row][0] = constantValue; + System.arraycopy(otherMatrix.getData()[row], 0, combinedData[row], 1, otherMatrix.getColumnDimension()); + } + + subMatrices.add(new Array2DRowRealMatrix(combinedData)); + } + + // Vertically stack all sub-matrices + List allRows = new ArrayList<>(); + for (RealMatrix subMatrix : subMatrices) { + double[][] subMatrixData = subMatrix.getData(); + for (double[] row : subMatrixData) { + allRows.add(row); + } + } + + double[][] stackedMatrixData = new double[allRows.size()][]; + for (int idx = 0; idx < allRows.size(); idx++) { + stackedMatrixData[idx] = allRows.get(idx); + } + + newMatrixArray[i] = new Array2DRowRealMatrix(stackedMatrixData); + } + + matrixArray = newMatrixArray; + + if (matrixArray[0].getColumnDimension() == numberOfColumns) { + break; + } + } + + return matrixArray; + } + + /** + * Map partitioned factor to final design using orthogonal-array produced + * by augmenting latin square. + */ + private static RealMatrix _mapPartitionsToDesign(List>> partitions, RealMatrix orthogonalArray) { + // Calculate max and min values in the matrix + double[][] orthogonalArrayData = orthogonalArray.getData(); + int maxPartitionIndex = (int) orthogonalArrayData[0][0]; // Initialize with first element + int minPartitionIndex = (int) orthogonalArrayData[0][0]; // Initialize with first element + + for (double[] row : orthogonalArrayData) { + for (double value : row) { + if (value > maxPartitionIndex) { + maxPartitionIndex = (int) value; + } + if (value < minPartitionIndex) { + minPartitionIndex = (int) value; + } + } + } + + maxPartitionIndex++; // Increment to match the original logic + if (!(partitions.size() == maxPartitionIndex && minPartitionIndex == 0)) { + throw new IllegalArgumentException("Orthogonal array indexing does not match partition structure"); + } + + List> mappingsList = new ArrayList<>(); + + for (double[] row : orthogonalArrayData) { + int[] rowIndexArray = new int[row.length]; + for (int i = 0; i < row.length; i++) { + rowIndexArray[i] = (int) row[i]; + } + + boolean hasEmptyPartition = false; + for (int factorIndex = 0; factorIndex < rowIndexArray.length; factorIndex++) { + int partitionIndex = rowIndexArray[factorIndex]; + if (partitions.get(partitionIndex).get(factorIndex).isEmpty()) { + hasEmptyPartition = true; + break; + } + } + + if (hasEmptyPartition) { + continue; + } + + List> partitionSets = new ArrayList<>(); + for (int factorIndex = 0; factorIndex < rowIndexArray.length; factorIndex++) { + int partitionIndex = rowIndexArray[factorIndex]; + partitionSets.add(partitions.get(partitionIndex).get(factorIndex)); + } + + List cartesianProduct = computeCartesianProduct(partitionSets); + mappingsList.add(cartesianProduct); + } + + // Combine all mappings + List allMappings = new ArrayList<>(); + for (List mapping : mappingsList) { + allMappings.addAll(mapping); + } + + double[][] resultData = new double[allMappings.size()][]; + for (int i = 0; i < allMappings.size(); i++) { + resultData[i] = new double[allMappings.get(i).length]; + for (int j = 0; j < allMappings.get(i).length; j++) { + resultData[i][j] = allMappings.get(i)[j]; + } + } + + return new Array2DRowRealMatrix(resultData); + } + + /** + * Compute cartesian product of lists of integers + */ + private static List computeCartesianProduct(List> listOfLists) { + List result = new ArrayList<>(); + if (listOfLists.isEmpty()) { + return result; + } + + computeCartesianProductRecursive(listOfLists, 0, new int[listOfLists.size()], result); + return result; + } + + private static void computeCartesianProductRecursive(List> listOfLists, int depth, int[] current, List result) { + if (depth == listOfLists.size()) { + result.add(current.clone()); + return; + } + + for (int value : listOfLists.get(depth)) { + current[depth] = value; + computeCartesianProductRecursive(listOfLists, depth + 1, current, result); + } + } + + /** + * Balanced partitioning of factors. + */ + private static List>> _makePartitions(int[] factorLevels, int numberOfPartitions) { + List>> partitions = new ArrayList<>(); + + for (int partitionIndex = 1; partitionIndex <= numberOfPartitions; partitionIndex++) { + List> partition = new ArrayList<>(); + + for (int factorLevel : factorLevels) { + List part = new ArrayList<>(); + + for (int levelIndex = 0; levelIndex < factorLevel; levelIndex++) { + int calculatedIndex = partitionIndex + levelIndex * numberOfPartitions; + if (calculatedIndex <= factorLevel) { + part.add(calculatedIndex); + } + } + + partition.add(part); + } + + partitions.add(partition); + } + + return partitions; + } + + /** + * Create a latin square of size n + */ + private static int[][] _makeLatinSquare(int n) { + int[][] latinSquare = new int[n][n]; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + latinSquare[i][j] = (i + j) % n; + } + } + + return latinSquare; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeHalton.java b/src/main/java/com/jdoe/algorithms/DoeHalton.java new file mode 100644 index 0000000..a5a01cb --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeHalton.java @@ -0,0 +1,128 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.util.Precision; + +/** + * This module implements the Halton sequence, a low-discrepancy quasirandom sequence + * used in numerical integration, sampling, and global optimization tasks. + * + * The Halton sequence generates points in a unit hypercube [0, 1]^d using radical + * inversion with respect to a sequence of prime number bases. It is especially + * useful in high-dimensional integration where uniformity and low correlation + * between sample points are desired. + * + * Each dimension in the Halton sequence uses a unique base (a prime number), and + * points are computed using the van der Corput sequence in that base. + */ +public class DoeHalton { + + /** + * Generate a Halton sequence in a given dimension. + * + * The Halton sequence is a low-discrepancy, quasi-random point set commonly + * used in numerical integration, sampling, and global optimization. Each + * dimension uses a different prime base to generate values via the van der Corput + * sequence. + * + * @param numPoints Number of points to generate in the sequence + * @param dimension Number of dimensions (features) of the sequence + * @param skip Number of initial points in the sequence to skip (default 0) + * @return The generated Halton sequence points + */ + public static double[][] haltonSequence(int numPoints, int dimension, int skip) { + int[] bases = nextPrimes(dimension); + + // Preallocate the output array + double[][] samples = new double[numPoints][dimension]; + + for (int dim = 0; dim < dimension; dim++) { + int base = bases[dim]; + for (int i = 0; i < numPoints; i++) { + int index = i + skip; + samples[i][dim] = vanDerCorput(index, base); + } + } + + return samples; + } + + /** + * Overloaded method with default skip value (0) + */ + public static double[][] haltonSequence(int numPoints, int dimension) { + return haltonSequence(numPoints, dimension, 0); + } + + /** + * Compute a single value of the van der Corput sequence. + * + * The van der Corput sequence generates low-discrepancy values in [0, 1) using + * radical inversion in a specified base. + * + * @param index The index in the sequence + * @param base The base to use (must be >= 2) + * @return The van der Corput value at the given index and base + */ + private static double vanDerCorput(int index, int base) { + double result = 0.0; + double f = 1.0 / base; + int currentIndex = index; + + while (currentIndex > 0) { + int mod = currentIndex % base; + result += mod * f; + f /= base; + currentIndex /= base; + } + + return result; + } + + /** + * Generate the first n prime numbers. + * + * @param n Number of prime numbers to generate + * @return Array containing the first n prime numbers + */ + private static int[] nextPrimes(int n) { + int[] primes = new int[n]; + int count = 0; + int candidate = 2; + + while (count < n) { + if (isPrime(candidate)) { + primes[count] = candidate; + count++; + } + candidate++; + } + + return primes; + } + + /** + * Check whether a number is prime. + * + * @param n The number to check + * @return True if n is a prime number, False otherwise + */ + private static boolean isPrime(int n) { + if (n < 2) { + return false; + } + if (n == 2) { + return true; + } + if (n % 2 == 0) { + return false; + } + + for (int i = 3; i * i <= n; i += 2) { + if (n % i == 0) { + return false; + } + } + + return true; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeKorobov.java b/src/main/java/com/jdoe/algorithms/DoeKorobov.java new file mode 100644 index 0000000..c294061 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeKorobov.java @@ -0,0 +1,85 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +import java.util.Random; + +/** + * This module implements the Korobov lattice generator for quasi-random sampling. + *

+ * Korobov lattices are a subclass of rank-1 lattice rules used for generating + * low-discrepancy sequences. These sequences are widely applied in quasi-Monte Carlo + * methods, global optimization, and numerical integration of high-dimensional functions. + *

+ * The construction uses modular arithmetic to build the generator vector using a + * single integer parameter. When the number of points and the generator are coprime, + * the resulting design exhibits Latin Hypercube-like properties with excellent + * uniform coverage. + *

+ * This implementation is based on a simplified rank-1 lattice formulation and + * uses a linear-time algorithm. + */ +public class DoeKorobov { + + /** + * Generate a Korobov lattice design matrix. + *

+ * Korobov lattices form a class of low-discrepancy sequences for quasi-Monte Carlo + * methods. They are constructed using a generator vector derived from modular + * exponentiation of a single integer. The generated matrix represents samples in + * a uniform virtual grid. + * + * @param numPoints Number of design points to generate + * @param dimension Number of dimensions in the design space + * @param generatorParam Generator parameter used in modular construction. If null, a random value + * in [2, numPoints) is selected + * @return design Integer-valued design matrix corresponding to bins on a modular grid + * @throws IllegalArgumentException if generatorParam is not greater than 1 + * @note The Korobov method is a special case of a rank-1 lattice. The generator vector + * is defined as: z_i = (a^i) mod N for i = 0 to d-1, where a is the generatorParam and N is numPoints. + * The resulting design has uniformity properties ideal for integration and + * high-dimensional optimization. + *

+ * To ensure good coverage, it's recommended that gcd(generatorParam, numPoints) == 1. + */ + public static RealMatrix korobovSequence(int numPoints, int dimension, Integer generatorParam) { + int actualGeneratorParam = generatorParam != null ? generatorParam : new Random().nextInt(numPoints - 2) + 2; + actualGeneratorParam %= numPoints; + if (actualGeneratorParam <= 1) { + throw new IllegalArgumentException("generatorParam must be greater than 1."); + } + + int[] generatorVector = new int[dimension]; + generatorVector[0] = 1; + for (int i = 1; i < dimension; i++) { + generatorVector[i] = (actualGeneratorParam * generatorVector[i - 1]) % numPoints; + } + + return rank1Lattice(numPoints, dimension, generatorVector); + } + + /** + * Overloaded method with default generatorParam (null) + */ + public static RealMatrix korobovSequence(int numPoints, int dimension) { + return korobovSequence(numPoints, dimension, null); + } + + /** + * Internal implementation of rank-1 lattice generation + * This is a placeholder method that would need to be implemented separately + * based on the actual rank1_lattice function from Python + */ + private static RealMatrix rank1Lattice(int numPoints, int dimension, int[] generatorVector) { + double[][] result = new double[numPoints][dimension]; + + for (int i = 0; i < numPoints; i++) { + for (int j = 0; j < dimension; j++) { + result[i][j] = ((i * generatorVector[j]) % numPoints) / (double) numPoints; + } + } + + return new Array2DRowRealMatrix(result); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeLhs.java b/src/main/java/com/jdoe/algorithms/DoeLhs.java new file mode 100644 index 0000000..71da9af --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeLhs.java @@ -0,0 +1,551 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.Random; + +/** + * Generate a latin-hypercube design + *

+ * This code was originally published by the following individuals for use with + * Scilab: + * Copyright (C) 2012 - 2013 - Michael Baudin + * Copyright (C) 2012 - Maria Christopoulou + * Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + * Copyright (C) 2009 - Yann Collette + * Copyright (C) 2009 - CEA - Jean-Marc Martinez + *

+ * Much thanks goes to these individuals. It has been converted to Python by + * Abraham Lee. + */ +public class DoeLhs { + + /** + * Generate a latin-hypercube design + * + * @param n The number of factors to generate samples for + * @param samples The number of samples to generate for each factor (Default: n) + * @param criterion Allowable values are "center" or "c", "maximin" or "m", + * "centermaximin" or "cm", and "correlation" or "corr". If no value + * given, the design is simply randomized. + * @param iterations The number of iterations in the maximin and correlations algorithms + * (Default: 5). + * @param randomGenerator RandomGenerator which controls random draws + * @param correlationMatrix Enforce correlation between factors (only used in lhsmu) + * @return An n-by-samples design matrix that has been normalized so factor values + * are uniformly spaced between zero and one. + */ + public static RealMatrix lhs( + int n, + Integer samples, + String criterion, + Integer iterations, + RandomGenerator randomGenerator, + double[][] correlationMatrix) { + + if (samples == null) { + samples = n; + } + + if (randomGenerator == null) { + randomGenerator = new JDKRandomGenerator(); + } + + if (criterion != null) { + if (!isValidCriterion(criterion)) { + throw new IllegalArgumentException(String.format("Invalid value for \"criterion\": %s", criterion)); + } + } + + RealMatrix H = null; + + if (criterion == null) { + H = _lhsClassic(n, samples, randomGenerator); + } else { + if (criterion.equalsIgnoreCase("center") || criterion.equalsIgnoreCase("c")) { + H = _lhsCentered(n, samples, randomGenerator); + } else if (criterion.equalsIgnoreCase("maximin") || criterion.equalsIgnoreCase("m")) { + if (iterations == null) iterations = 5; + H = _lhsMaximin(n, samples, iterations, "maximin", randomGenerator); + } else if (criterion.equalsIgnoreCase("centermaximin") || criterion.equalsIgnoreCase("cm")) { + if (iterations == null) iterations = 5; + H = _lhsMaximin(n, samples, iterations, "centermaximin", randomGenerator); + } else if (criterion.equalsIgnoreCase("correlation") || criterion.equalsIgnoreCase("corr")) { + if (iterations == null) iterations = 5; + H = _lhsCorrelate(n, samples, iterations, randomGenerator); + } else if (criterion.equalsIgnoreCase("lhsmu")) { + // as specified by the paper. M is set to 5 + H = _lhsMu(n, samples, correlationMatrix, randomGenerator, 5); + } + } + + return H; + } + + /** + * Overloaded method with default parameters + */ + public static RealMatrix lhs(int n) { + return lhs(n, null, null, null, null, null); + } + + /** + * Overloaded method with samples parameter + */ + public static RealMatrix lhs(int n, int samples) { + return lhs(n, samples, null, null, null, null); + } + + /** + * Overloaded method with samples and criterion parameters + */ + public static RealMatrix lhs(int n, int samples, String criterion) { + return lhs(n, samples, criterion, null, null, null); + } + + /** + * Overloaded method with samples, criterion, and iterations parameters + */ + public static RealMatrix lhs(int n, int samples, String criterion, int iterations) { + return lhs(n, samples, criterion, iterations, null, null); + } + + /** + * Check if the criterion is valid + */ + private static boolean isValidCriterion(String criterion) { + String lowerCrit = criterion.toLowerCase(); + return lowerCrit.equals("center") || lowerCrit.equals("c") || + lowerCrit.equals("maximin") || lowerCrit.equals("m") || + lowerCrit.equals("centermaximin") || lowerCrit.equals("cm") || + lowerCrit.equals("correlation") || lowerCrit.equals("corr") || + lowerCrit.equals("lhsmu"); + } + + /** + * Classic LHS implementation + */ + private static RealMatrix _lhsClassic(int n, int samples, RandomGenerator randomGenerator) { + // Generate the intervals + double[] cut = linspace(0, 1, samples + 1); + + // Fill points uniformly in each interval + double[][] u = new double[samples][n]; + for (int i = 0; i < samples; i++) { + for (int j = 0; j < n; j++) { + u[i][j] = randomGenerator.nextDouble(); + } + } + + double[] a = new double[samples]; + double[] b = new double[samples]; + System.arraycopy(cut, 0, a, 0, samples); + System.arraycopy(cut, 1, b, 0, samples); + + double[][] rdpoints = new double[samples][n]; + for (int j = 0; j < n; j++) { + for (int i = 0; i < samples; i++) { + rdpoints[i][j] = u[i][j] * (b[i] - a[i]) + a[i]; + } + } + + // Make the random pairings + double[][] H = new double[samples][n]; + for (int j = 0; j < n; j++) { + int[] order = randomPermutation(samples, randomGenerator); + for (int i = 0; i < samples; i++) { + H[i][j] = rdpoints[order[i]][j]; + } + } + + return new Array2DRowRealMatrix(H); + } + + /** + * Centered LHS implementation + */ + private static RealMatrix _lhsCentered(int n, int samples, RandomGenerator randomGenerator) { + // Generate the intervals + double[] cut = linspace(0, 1, samples + 1); + + double[] a = new double[samples]; + double[] b = new double[samples]; + System.arraycopy(cut, 0, a, 0, samples); + System.arraycopy(cut, 1, b, 0, samples); + + double[] center = new double[samples]; + for (int i = 0; i < samples; i++) { + center[i] = (a[i] + b[i]) / 2.0; + } + + // Make the random pairings + double[][] H = new double[samples][n]; + for (int j = 0; j < n; j++) { + double[] permutedCenter = shuffleArray(center, randomGenerator); + for (int i = 0; i < samples; i++) { + H[i][j] = permutedCenter[i]; + } + } + + return new Array2DRowRealMatrix(H); + } + + /** + * Maximin LHS implementation + */ + private static RealMatrix _lhsMaximin(int n, int samples, int iterations, String lhsType, RandomGenerator randomGenerator) { + double maxDist = 0; + RealMatrix bestH = null; + + // Maximize the minimum distance between points + for (int i = 0; i < iterations; i++) { + RealMatrix hCandidate; + if ("maximin".equals(lhsType)) { + hCandidate = _lhsClassic(n, samples, randomGenerator); + } else { + hCandidate = _lhsCentered(n, samples, randomGenerator); + } + + // Calculate pairwise distances + double minDist = calculateMinDistance(hCandidate); + if (maxDist < minDist) { + maxDist = minDist; + bestH = hCandidate.copy(); + } + } + + return bestH; + } + + /** + * Correlation-based LHS implementation + */ + private static RealMatrix _lhsCorrelate(int n, int samples, int iterations, RandomGenerator randomGenerator) { + double minCorr = Double.POSITIVE_INFINITY; + RealMatrix bestH = null; + + // Minimize the components correlation coefficients + for (int i = 0; i < iterations; i++) { + // Generate a random LHS + RealMatrix hCandidate = _lhsClassic(n, samples, randomGenerator); + double maxAbsCorr = calculateMaxAbsCorrelation(hCandidate); + + if (maxAbsCorr < minCorr) { + minCorr = maxAbsCorr; + bestH = hCandidate.copy(); + } + } + + return bestH; + } + + /** + * LHS-MU implementation + */ + private static RealMatrix _lhsMu(int n, Integer samples, double[][] corr, RandomGenerator randomGenerator, int M) { + if (samples == null) { + samples = n; + } + + int I = M * samples; + + // Generate random points + double[][] rdPoints = new double[I][n]; + for (int i = 0; i < I; i++) { + for (int j = 0; j < n; j++) { + rdPoints[i][j] = randomGenerator.nextDouble(); + } + } + + // Calculate distance matrix + double[][] dist = calculateDistanceMatrix(rdPoints); + + // Mask diagonal elements + for (int i = 0; i < I; i++) { + dist[i][i] = Double.NaN; // Using NaN as mask + } + + int[] indexRm = new int[I - samples]; + int rmCount = 0; + + while (rmCount < I - samples) { + // Find minimum average distance + int minL = findMinAvgDistance(dist); + + // Mask this row and column + for (int i = 0; i < I; i++) { + dist[minL][i] = Double.NaN; + dist[i][minL] = Double.NaN; + } + + indexRm[rmCount] = minL; + rmCount++; + } + + // Remove selected points + double[][] filteredPoints = removeRows(rdPoints, indexRm); + + if (corr != null) { + // Apply correlation transformation using Cholesky decomposition + return applyCorrelationTransformation(filteredPoints, corr, randomGenerator); + } else { + return rankOrderTransform(filteredPoints, samples, randomGenerator); + } + } + + // Helper methods + + private static double[] linspace(double start, double end, int num) { + double[] result = new double[num]; + if (num == 1) { + result[0] = start; + return result; + } + double step = (end - start) / (num - 1); + for (int i = 0; i < num; i++) { + result[i] = start + i * step; + } + return result; + } + + private static int[] randomPermutation(int n, RandomGenerator randomGenerator) { + List list = new ArrayList<>(); + for (int i = 0; i < n; i++) { + list.add(i); + } + Collections.shuffle(list, new Random(randomGenerator.nextInt())); + return list.stream().mapToInt(Integer::intValue).toArray(); + } + + private static double[] shuffleArray(double[] array, RandomGenerator randomGenerator) { + double[] newArray = array.clone(); + List list = new ArrayList<>(); + for (double d : newArray) { + list.add(d); + } + Collections.shuffle(list, new Random(randomGenerator.nextInt())); + for (int i = 0; i < newArray.length; i++) { + newArray[i] = list.get(i); + } + return newArray; + } + + private static double calculateMinDistance(RealMatrix matrix) { + double minDist = Double.POSITIVE_INFINITY; + int rows = matrix.getRowDimension(); + + for (int i = 0; i < rows; i++) { + for (int j = i + 1; j < rows; j++) { + double dist = euclideanDistance(matrix.getRow(i), matrix.getRow(j)); + if (dist < minDist) { + minDist = dist; + } + } + } + + return minDist; + } + + private static double euclideanDistance(double[] a, double[] b) { + double sum = 0; + for (int i = 0; i < a.length; i++) { + sum += Math.pow(a[i] - b[i], 2); + } + return Math.sqrt(sum); + } + + private static double calculateMaxAbsCorrelation(RealMatrix matrix) { + int n = matrix.getColumnDimension(); + double[][] corrMatrix = new double[n][n]; + + // Calculate correlation matrix + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i == j) { + corrMatrix[i][j] = 1.0; + } else { + corrMatrix[i][j] = correlation(matrix.getColumn(i), matrix.getColumn(j)); + } + } + } + + // Find maximum absolute correlation excluding diagonal + double maxAbsCorr = 0; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i != j) { + double absCorr = Math.abs(corrMatrix[i][j]); + if (absCorr > maxAbsCorr) { + maxAbsCorr = absCorr; + } + } + } + } + + return maxAbsCorr; + } + + private static double correlation(double[] x, double[] y) { + int n = x.length; + double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0, sumY2 = 0; + + for (int i = 0; i < n; i++) { + sumX += x[i]; + sumY += y[i]; + sumXY += x[i] * y[i]; + sumX2 += x[i] * x[i]; + sumY2 += y[i] * y[i]; + } + + double numerator = n * sumXY - sumX * sumY; + double denominator = Math.sqrt((n * sumX2 - sumX * sumX) * (n * sumY2 - sumY * sumY)); + + if (denominator == 0) return 0; + return numerator / denominator; + } + + private static double[][] calculateDistanceMatrix(double[][] points) { + int n = points.length; + double[][] dist = new double[n][n]; + + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i == j) { + dist[i][j] = 0; + } else { + dist[i][j] = euclideanDistance(points[i], points[j]); + } + } + } + + return dist; + } + + private static int findMinAvgDistance(double[][] dist) { + int n = dist.length; + double minAvg = Double.POSITIVE_INFINITY; + int minIdx = 0; + + for (int i = 0; i < n; i++) { + double sum = 0; + int count = 0; + for (int j = 0; j < n; j++) { + if (!Double.isNaN(dist[i][j])) { + sum += dist[i][j]; + count++; + } + } + + if (count > 0) { + double avg = sum / count; + if (avg < minAvg) { + minAvg = avg; + minIdx = i; + } + } + } + + return minIdx; + } + + private static double[][] removeRows(double[][] matrix, int[] indicesToRemove) { + List toRemove = new ArrayList<>(); + for (int idx : indicesToRemove) { + toRemove.add(idx); + } + Collections.sort(toRemove, Collections.reverseOrder()); + + List remainingRows = new ArrayList<>(); + for (int i = 0; i < matrix.length; i++) { + if (!toRemove.contains(i)) { + remainingRows.add(matrix[i]); + } + } + + return remainingRows.toArray(new double[0][0]); + } + + private static RealMatrix applyCorrelationTransformation(double[][] points, double[][] corr, RandomGenerator randomGenerator) { + // This is a simplified implementation + // A full implementation would require proper Cholesky decomposition + return new Array2DRowRealMatrix(points); + } + + private static RealMatrix rankOrderTransform(double[][] points, int samples, RandomGenerator randomGenerator) { + int n = points[0].length; + double[][] result = new double[samples][n]; + int m = points.length; + + for (int j = 0; j < n; j++) { + // Get column and sort to determine ranks + double[] col = new double[m]; + for (int i = 0; i < m; i++) { + col[i] = points[i][j]; + } + + // Create pairs of (value, originalIndex) and sort by value + int[] indices = new int[m]; + for (int i = 0; i < m; i++) { + indices[i] = i; + } + + // Sort indices based on values + quickSort(col, indices, 0, m - 1); + + // Assign values based on rank + for (int i = 0; i < samples; i++) { + double low = (double) i / samples; + double high = (double) (i + 1) / samples; + result[i][j] = low + randomGenerator.nextDouble() * (high - low); + } + } + + return new Array2DRowRealMatrix(result); + } + + // Quick sort implementation for sorting with indices + private static void quickSort(double[] arr, int[] indices, int low, int high) { + if (low < high) { + int pi = partition(arr, indices, low, high); + quickSort(arr, indices, low, pi - 1); + quickSort(arr, indices, pi + 1, high); + } + } + + private static int partition(double[] arr, int[] indices, int low, int high) { + double pivot = arr[high]; + int i = (low - 1); + + for (int j = low; j < high; j++) { + if (arr[j] <= pivot) { + i++; + // Swap arr[i] and arr[j] + double temp = arr[i]; + arr[i] = arr[j]; + arr[j] = temp; + + // Also swap corresponding indices + int tempIdx = indices[i]; + indices[i] = indices[j]; + indices[j] = tempIdx; + } + } + + // Swap arr[i+1] and arr[high] (or pivot) + double temp = arr[i + 1]; + arr[i + 1] = arr[high]; + arr[high] = temp; + + int tempIdx = indices[i + 1]; + indices[i + 1] = indices[high]; + indices[high] = tempIdx; + + return i + 1; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoePlackettBurman.java b/src/main/java/com/jdoe/algorithms/DoePlackettBurman.java new file mode 100644 index 0000000..698ca89 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoePlackettBurman.java @@ -0,0 +1,251 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.*; + +/** + * Generate a Plackett-Burman design + *

+ * This code was originally published by the following individuals for use with + * Scilab: + * Copyright (C) 2012 - 2013 - Michael Baudin + * Copyright (C) 2012 - Maria Christopoulou + * Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + * Copyright (C) 2009 - Yann Collette + * Copyright (C) 2009 - CEA - Jean-Marc Martinez + *

+ * Much thanks goes to these individuals. It has been converted to Python by + * Abraham Lee. + */ +public class DoePlackettBurman { + + /** + * Generate a Plackett-Burman design + * + * @param n The number of factors to create a matrix for + * @return An orthogonal design matrix with n columns, one for each factor, and + * the number of rows being the next multiple of 4 higher than n (e.g., + * for 1-3 factors there are 4 rows, for 4-7 factors there are 8 rows, + * etc.) + * @throws IllegalArgumentException if n is not a positive integer + */ + public static RealMatrix pbdesign(int n) { + if (n <= 0) { + throw new IllegalArgumentException("Number of factors must be a positive integer"); + } + + int keep = n; + n = 4 * (int) Math.ceil(n / 4.0); // calculate the correct number of rows (multiple of 4) + + // Check if n is a valid multiple of 4 + if (!isValidMultipleOfFour(n)) { + throw new IllegalArgumentException("Invalid inputs. n must be a multiple of 4."); + } + + RealMatrix H; + + // Determine which base matrix to use based on n + if (n == 4) { + // N = 4 (base case) + H = new Array2DRowRealMatrix(new double[][]{{1}}); + } else if (n == 12) { + // N = 12 + H = createTwelveByTwelveMatrix(); + } else if (n == 20) { + // N = 20 + H = createTwentyByTwentyMatrix(); + } else { + // For larger matrices, find the base size and build up + int baseSize = findBaseSize(n); + if (baseSize == 4) { + H = new Array2DRowRealMatrix(new double[][]{{1}}); + } else if (baseSize == 12) { + H = createTwelveByTwelveMatrix(); + } else if (baseSize == 20) { + H = createTwentyByTwentyMatrix(); + } else { + throw new IllegalArgumentException("Unsupported matrix size: " + n); + } + + // Calculate how many times to expand the base matrix + int expansionFactor = n / H.getRowDimension(); + int power = (int) (Math.log(expansionFactor) / Math.log(2)); + + // Kronecker product construction + for (int i = 0; i < power; i++) { + H = kroneckerExpansion(H); + } + } + + // Reduce the size of the matrix as needed + int numRows = H.getRowDimension(); + int numCols = Math.min(H.getColumnDimension(), keep + 1); + RealMatrix result = new Array2DRowRealMatrix(numRows, numCols); + for (int i = 0; i < numRows; i++) { + for (int j = 1; j < numCols; j++) { // Skip first column (index 0) + result.setEntry(i, j - 1, H.getEntry(i, j)); + } + } + + // Flip the matrix upside down + return flipud(result); + } + + /** + * Check if n is a valid multiple of 4 that can be formed using the allowed base matrices + */ + private static boolean isValidMultipleOfFour(int n) { + // Check if n can be expressed as one of the base sizes multiplied by powers of 2 + while (n % 2 == 0) { + if (n == 4 || n == 12 || n == 20) { + return true; + } + n /= 2; + } + return n == 4 || n == 12 || n == 20; + } + + /** + * Find the base size for a given n + */ + private static int findBaseSize(int n) { + int temp = n; + while (temp % 2 == 0) { + if (temp == 4 || temp == 12 || temp == 20) { + return temp; + } + temp /= 2; + } + return temp; // Should be 4, 12, or 20 + } + + /** + * Create the 12x12 base matrix + */ + private static RealMatrix createTwelveByTwelveMatrix() { + double[][] matrix = new double[12][12]; + + // First row: all ones + for (int j = 0; j < 12; j++) { + matrix[0][j] = 1; + } + + // Remaining rows + for (int i = 1; i < 12; i++) { + matrix[i][0] = 1; // First column + } + + // Toeplitz matrix part + int[] colFirst = {-1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1}; // First column after first element + int[] rowFirst = {-1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1}; // First row after first element + + for (int i = 1; i < 12; i++) { + for (int j = 1; j < 12; j++) { + if (i <= j) { + // Moving up along anti-diagonals + int idx = j - i; + matrix[i][j] = (idx < rowFirst.length) ? rowFirst[idx] : matrix[i-1][j-1]; + } else { + // Moving down along anti-diagonals + int idx = i - j; + matrix[i][j] = (idx < colFirst.length) ? colFirst[idx] : matrix[i-1][j-1]; + } + } + } + + return new Array2DRowRealMatrix(matrix); + } + + /** + * Create the 20x20 base matrix + */ + private static RealMatrix createTwentyByTwentyMatrix() { + double[][] matrix = new double[20][20]; + + // First row: all ones + for (int j = 0; j < 20; j++) { + matrix[0][j] = 1; + } + + // Remaining rows + for (int i = 1; i < 20; i++) { + matrix[i][0] = 1; // First column + } + + // Hankel matrix part + int[] colFirst = {-1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1}; + int[] rowLast = {1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1}; + + for (int i = 1; i < 20; i++) { + for (int j = 1; j < 20; j++) { + int idx = i + j - 2; // For Hankel: constant along anti-diagonals + if (idx < colFirst.length) { + matrix[i][j] = colFirst[idx]; + } else if (idx - (colFirst.length - 1) < rowLast.length) { + matrix[i][j] = rowLast[idx - (colFirst.length - 1)]; + } else { + // Fallback: replicate pattern + matrix[i][j] = matrix[i-1][j-1]; + } + } + } + + return new Array2DRowRealMatrix(matrix); + } + + /** + * Perform Kronecker expansion of the matrix + */ + private static RealMatrix kroneckerExpansion(RealMatrix H) { + int rows = H.getRowDimension(); + int cols = H.getColumnDimension(); + + RealMatrix expanded = new Array2DRowRealMatrix(2 * rows, 2 * cols); + + // Upper left: H + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + expanded.setEntry(i, j, H.getEntry(i, j)); + } + } + + // Upper right: H + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + expanded.setEntry(i, j + cols, H.getEntry(i, j)); + } + } + + // Lower left: H + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + expanded.setEntry(i + rows, j, H.getEntry(i, j)); + } + } + + // Lower right: -H + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + expanded.setEntry(i + rows, j + cols, -H.getEntry(i, j)); + } + } + + return expanded; + } + + /** + * Flip the matrix upside down + */ + private static RealMatrix flipud(RealMatrix matrix) { + int rows = matrix.getRowDimension(); + int cols = matrix.getColumnDimension(); + RealMatrix flipped = new Array2DRowRealMatrix(rows, cols); + + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + flipped.setEntry(i, j, matrix.getEntry(rows - 1 - i, j)); + } + } + + return flipped; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeRandomKMeans.java b/src/main/java/com/jdoe/algorithms/DoeRandomKMeans.java new file mode 100644 index 0000000..b58e18b --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeRandomKMeans.java @@ -0,0 +1,134 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.*; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +/** + * MacQueen's K-Means algorithm. + *

+ * This implementation generates cluster centers using the MacQueen's K-Means algorithm + * in a unit hypercube [0, 1]^dimension. + */ +public class DoeRandomKMeans { + + /** + * MacQueen's K-Means algorithm. + * + * @param numPoints Number of cluster centers to generate + * @param dimension Dimensionality of the space + * @param numSteps Number of iterations. Defaults to 100 * numPoints + * @param initialPoints Initial cluster centers. If null, random points in [0, 1]^dimension are used + * @param randomGenerator RandomGenerator for reproducibility + * @return Array of shape (numPoints, dimension) containing the cluster centers + * @throws IllegalArgumentException if initialPoints doesn't match expected shape or bounds + */ + public static RealMatrix randomKMeans( + int numPoints, + int dimension, + Integer numSteps, + RealMatrix initialPoints, + RandomGenerator randomGenerator) { + + if (randomGenerator == null) { + randomGenerator = new JDKRandomGenerator(); + } + + if (numSteps == null) { + numSteps = 100 * numPoints; + } + + // Initialize cluster centers + RealMatrix clusterCenters; + if (initialPoints == null) { + clusterCenters = generateRandomMatrix(numPoints, dimension, randomGenerator); + } else { + if (initialPoints.getRowDimension() != numPoints || initialPoints.getColumnDimension() != dimension) { + throw new IllegalArgumentException("initialPoints must have shape (numPoints, dimension)"); + } + + // Validate that all points are in [0, 1]^dimension + for (int i = 0; i < initialPoints.getRowDimension(); i++) { + for (int j = 0; j < initialPoints.getColumnDimension(); j++) { + double value = initialPoints.getEntry(i, j); + if (value < 0.0 || value > 1.0) { + throw new IllegalArgumentException("initialPoints must be in [0, 1]^dimension"); + } + } + } + + clusterCenters = initialPoints.copy(); + } + + // Initialize counts for incremental mean + double[] counts = new double[numPoints]; + for (int i = 0; i < numPoints; i++) { + counts[i] = 1.0; + } + + for (int step = 0; step < numSteps; step++) { + // Sample a random point in the unit hypercube + double[] x = new double[dimension]; + for (int i = 0; i < dimension; i++) { + x[i] = randomGenerator.nextDouble(); + } + + // Compute Euclidean distances to cluster centers + double[] distances = new double[numPoints]; + for (int i = 0; i < numPoints; i++) { + double sumSquares = 0.0; + for (int j = 0; j < dimension; j++) { + double diff = clusterCenters.getEntry(i, j) - x[j]; + sumSquares += diff * diff; + } + distances[i] = Math.sqrt(sumSquares); + } + + // Find nearest cluster center + int idx = 0; + double minDist = distances[0]; + for (int i = 1; i < numPoints; i++) { + if (distances[i] < minDist) { + minDist = distances[i]; + idx = i; + } + } + + // Update cluster center incrementally (MacQueen's update) + for (int j = 0; j < dimension; j++) { + double newValue = (counts[idx] * clusterCenters.getEntry(idx, j) + x[j]) / (counts[idx] + 1); + clusterCenters.setEntry(idx, j, newValue); + } + counts[idx] += 1; + } + + return clusterCenters; + } + + /** + * Overloaded method with default parameters + */ + public static RealMatrix randomKMeans(int numPoints, int dimension) { + return randomKMeans(numPoints, dimension, null, null, null); + } + + /** + * Overloaded method with numSteps parameter + */ + public static RealMatrix randomKMeans(int numPoints, int dimension, int numSteps) { + return randomKMeans(numPoints, dimension, numSteps, null, null); + } + + /** + * Generate a random matrix with values in [0, 1] + */ + private static RealMatrix generateRandomMatrix(int rows, int cols, RandomGenerator randomGenerator) { + double[][] matrix = new double[rows][cols]; + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + matrix[i][j] = randomGenerator.nextDouble(); + } + } + return new Array2DRowRealMatrix(matrix); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeRandomUniform.java b/src/main/java/com/jdoe/algorithms/DoeRandomUniform.java new file mode 100644 index 0000000..31a7b39 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeRandomUniform.java @@ -0,0 +1,50 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +/** + * Generate random samples from a uniform distribution over [0, 1). + *

+ * This function returns an array of shape [(num_points, dimension)](file:///home/noor/Downloads/pyDOE3-master/doc/conf.py#L31-L36) where each entry + * is drawn from a uniform distribution on the half-open interval [0.0, 1.0). + */ +public class DoeRandomUniform { + + /** + * Generate random samples from a uniform distribution over [0, 1). + * + * @param numPoints Number of random points to generate (number of rows in the output array) + * @param dimension Dimensionality of each random point (number of columns in the output array) + * @param seed Random seed for reproducibility + * @return An array of shape [(num_points, dimension)](file:///home/noor/Downloads/pyDOE3-master/doc/conf.py#L31-L36) containing random samples + * from a uniform distribution over [0, 1) + */ + public static RealMatrix randomUniform(int numPoints, int dimension, Integer seed) { + RandomGenerator rng; + if (seed != null) { + rng = new JDKRandomGenerator(seed); + } else { + rng = new JDKRandomGenerator(); + } + + // Generate the random matrix + double[][] result = new double[numPoints][dimension]; + for (int i = 0; i < numPoints; i++) { + for (int j = 0; j < dimension; j++) { + result[i][j] = rng.nextDouble(); + } + } + + return new Array2DRowRealMatrix(result); + } + + /** + * Overloaded method with default seed (null) + */ + public static RealMatrix randomUniform(int numPoints, int dimension) { + return randomUniform(numPoints, dimension, null); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeRank1.java b/src/main/java/com/jdoe/algorithms/DoeRank1.java new file mode 100644 index 0000000..8cf0d58 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeRank1.java @@ -0,0 +1,84 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +/** + * Generate a rank-1 lattice design matrix. + *

+ * Rank-1 lattices are quasi-random designs used for numerical integration and + * high-dimensional sampling. This algorithm generates deterministic points with + * linear runtime. + */ +public class DoeRank1 { + + /** + * Generate a rank-1 lattice design matrix. + * + * @param numPoints The number of points to generate + * @param dimension The dimensionality of the space + * @param generatorVector A generator vector of length dimension. If null, one is randomly generated + * using integers in [2, numPoints) + * @param randomGenerator RandomGenerator for reproducibility + * @return The resulting integer-valued rank-1 lattice matrix. Each row represents + * a point in the design + * @throws IllegalArgumentException if generatorVector doesn't match expected shape + */ + public static RealMatrix rank1Lattice(int numPoints, int dimension, int[] generatorVector, RandomGenerator randomGenerator) { + if (randomGenerator == null) { + randomGenerator = new JDKRandomGenerator(); + } + + if (generatorVector == null) { + generatorVector = new int[dimension]; + for (int i = 0; i < dimension; i++) { + generatorVector[i] = randomGenerator.nextInt(numPoints - 2) + 2; // Random int in [2, numPoints) + } + } + + if (generatorVector.length != dimension) { + throw new IllegalArgumentException( + String.format("Expected generator_vector of length (%d), got %d", dimension, generatorVector.length) + ); + } + + // Apply modulo operation to ensure all values are within [0, numPoints) + for (int i = 0; i < generatorVector.length; i++) { + generatorVector[i] = ((generatorVector[i] % numPoints) + numPoints) % numPoints; // Handle negative values + } + + // Generate the points using modular arithmetic + int[][] points = new int[numPoints][dimension]; + for (int i = 0; i < numPoints; i++) { + for (int j = 0; j < dimension; j++) { + points[i][j] = (i * generatorVector[j]) % numPoints; + } + } + + // Convert to RealMatrix + double[][] result = new double[numPoints][dimension]; + for (int i = 0; i < numPoints; i++) { + for (int j = 0; j < dimension; j++) { + result[i][j] = points[i][j]; + } + } + + return new Array2DRowRealMatrix(result); + } + + /** + * Overloaded method with default random generator + */ + public static RealMatrix rank1Lattice(int numPoints, int dimension, int[] generatorVector) { + return rank1Lattice(numPoints, dimension, generatorVector, null); + } + + /** + * Overloaded method with default generator vector (null) + */ + public static RealMatrix rank1Lattice(int numPoints, int dimension) { + return rank1Lattice(numPoints, dimension, null, null); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeSaltelli.java b/src/main/java/com/jdoe/algorithms/DoeSaltelli.java new file mode 100644 index 0000000..d7dd249 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeSaltelli.java @@ -0,0 +1,157 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +/** + * Generate Saltelli samples using Sobol' sequences for sensitivity analysis. + *

+ * This module implements Saltelli's sampling scheme based on Sobol' sequences for + * global sensitivity analysis. It enables estimation of first-order, total-order, and + * (second-order optional) Sobol' sensitivity indices. The implementation relies on a + * custom Sobol' sequence generator. + *

+ * Compared to random or Latin Hypercube sampling, this method provides better + * convergence for variance-based sensitivity analysis using quasi-random low-discrepancy + * sequences. + *

+ * References: + * 1. Sobol', I.M., 2001. Global sensitivity indices for nonlinear mathematical models and + * their Monte Carlo estimates. Mathematics and Computers in Simulation, + * The Second IMACS Seminar on Monte Carlo Methods 55, 271-280. + * 2. Saltelli, A., 2002. Making best use of model evaluations to compute sensitivity indices. + * Computer Physics Communications, 145(2), 280-297. + * 3. Campolongo, F., Saltelli, A., Cariboni, J., 2011. + * From screening to quantitative sensitivity analysis. A unified approach. + * Computer Physics Communications 182, 978-988. + * 4. Owen, A. B., 2020. On dropping the first Sobol' point. arXiv:2008.08051 [cs, math, stat]. + */ +public class DoeSaltelli { + + /** + * Generate Saltelli samples using Sobol' sequences for sensitivity analysis. + * + * @param numVars Number of input variables (dimensions) + * @param N Base sample size (ideally a power of 2) + * @param calcSecondOrder If true, include second-order interaction terms. Default is true + * @param skipValues Number of Sobol' points to skip. If null, set automatically + * @param scramble Whether to use scrambling for Sobol' sequence. Default is false + * @param seed Random seed (only used if scramble=true) + * @return Matrix of shape (N * (2 * num_vars + 2), num_vars) if calc_second_order=true, + * or (N * (num_vars + 2), num_vars) otherwise. Contains Saltelli samples in [0, 1] + */ + public static RealMatrix saltelliSampling( + int numVars, + int N, + boolean calcSecondOrder, + Integer skipValues, + boolean scramble, + Integer seed) { + + int D = numVars; + + // Check if N is a power of 2 + if (!((N & (N - 1)) == 0 && N != 0)) { + System.out.println("Warning: N = " + N + " is not a power of 2. This may affect Sobol sequence convergence."); + } + + if (skipValues == null) { + skipValues = Math.max(powerOfTwoCeil(N), 16); + } + + // Generate base Sobol samples: shape (N + skip_values, 2*D) + // Fixed method call to match the correct signature + RealMatrix base = DoeSobol.sobolSequence(N + skipValues, 2 * D, scramble, seed, null, 0); + + int totalSamples; + if (calcSecondOrder) { + totalSamples = N * (2 * D + 2); + } else { + totalSamples = N * (D + 2); + } + + double[][] saltelliMatrix = new double[totalSamples][D]; + int idx = 0; + + for (int i = skipValues; i < skipValues + N; i++) { + // Extract A and B matrices + double[] A = new double[D]; + double[] B = new double[D]; + for (int j = 0; j < D; j++) { + A[j] = base.getEntry(i, j); + B[j] = base.getEntry(i, D + j); + } + + // Matrix A + for (int j = 0; j < D; j++) { + saltelliMatrix[idx][j] = A[j]; + } + idx++; + + // Cross A_Bi + for (int j = 0; j < D; j++) { + double[] C = A.clone(); + C[j] = B[j]; + for (int k = 0; k < D; k++) { + saltelliMatrix[idx][k] = C[k]; + } + idx++; + } + + // Cross B_Ai (only if calc_second_order) + if (calcSecondOrder) { + for (int j = 0; j < D; j++) { + double[] C = B.clone(); + C[j] = A[j]; + for (int k = 0; k < D; k++) { + saltelliMatrix[idx][k] = C[k]; + } + idx++; + } + } + + // Matrix B + for (int j = 0; j < D; j++) { + saltelliMatrix[idx][j] = B[j]; + } + idx++; + } + + return new Array2DRowRealMatrix(saltelliMatrix); + } + + /** + * Overloaded method with default parameters + */ + public static RealMatrix saltelliSampling(int numVars, int N) { + return saltelliSampling(numVars, N, true, null, false, null); + } + + /** + * Overloaded method with calcSecondOrder parameter + */ + public static RealMatrix saltelliSampling(int numVars, int N, boolean calcSecondOrder) { + return saltelliSampling(numVars, N, calcSecondOrder, null, false, null); + } + + /** + * Overloaded method with calcSecondOrder and skipValues parameters + */ + public static RealMatrix saltelliSampling(int numVars, int N, boolean calcSecondOrder, int skipValues) { + return saltelliSampling(numVars, N, calcSecondOrder, skipValues, false, null); + } + + /** + * Calculate the smallest power of 2 that is greater than or equal to n + */ + private static int powerOfTwoCeil(int n) { + if (n <= 1) return 1; + int result = 1; + while (result < n) { + result <<= 1; + } + return result; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeSobol.java b/src/main/java/com/jdoe/algorithms/DoeSobol.java new file mode 100644 index 0000000..fee6093 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeSobol.java @@ -0,0 +1,129 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.random.SobolSequenceGenerator; + +/** + * Generate a Sobol' sequence (quasi-random design matrix). + *

+ * This module implements the Sobol' sequence generator for quasi-random sampling. + *

+ * The Sobol' sequence is a type of low-discrepancy sequence widely used in global + * optimization, numerical integration (e.g., Monte Carlo methods), machine learning, + * and design of experiments. Compared to purely random samples, Sobol' sequences + * exhibit better uniformity in multi-dimensional space. + *

+ * The sequence generation is based on direction numbers and bitwise operations. + * Scrambling can optionally be applied to enhance uniformity and reduce correlation + * artifacts. + *

+ * References: + * Sobol', I. M. (1967). "Distribution of points in a cube and approximate evaluation + * of integrals." Zh. Vych. Mat. Mat. Fiz., 7: 784-802 (in Russian); + * U.S.S.R. Comput. Maths. Math. Phys., 7: 86-112 (in English). + */ +public class DoeSobol { + + /** + * Generate a Sobol' sequence (quasi-random design matrix). + * + * @param n Number of points to generate + * @param d Dimension of the space (must be <= 21201) + * @param scramble Whether to apply Owen scrambling. Default is false + * @param seed Seed for the random number generator (used only when scramble=true) + * @param bounds Bounds for each dimension. Each element must be a (min, max) pair. + * If provided, the output will be scaled accordingly + * @param skip Number of initial points to skip (i.e., fast-forward in the sequence). Default is 0 + * @param usePowOf2 If true, ensures n is a power of 2 for best balance and coverage. + * Non-power-of-two n values will be rounded up to the next power of 2 + * @return Array of Sobol' points in [0, 1)^d, or scaled to bounds if provided + */ + public static RealMatrix sobolSequence( + int n, int d, boolean scramble, Integer seed, double[][] bounds, int skip, boolean usePowOf2) { + + // Note: The SobolSequenceGenerator in Apache Commons Math doesn't support scrambling + // directly, so we're creating it with just the dimension + SobolSequenceGenerator sobolGen = new SobolSequenceGenerator(d); + + // Fast forward if needed + for (int i = 0; i < skip; i++) { + sobolGen.nextVector(); + } + + int actualN = n; + if (usePowOf2) { + // Ensure n is power of 2 for best balance properties + if (!isPowerOfTwo(n)) { + actualN = nextPowerOfTwo(n); + } + } + + double[][] samples = new double[actualN][d]; + + // Generate samples + for (int i = 0; i < actualN; i++) { + double[] point = sobolGen.nextVector(); + System.arraycopy(point, 0, samples[i], 0, d); + } + + // Apply bounds scaling if provided + if (bounds != null) { + if (bounds.length != d || bounds[0].length != 2) { + throw new IllegalArgumentException( + String.format("`bounds` must be a (d, 2) array, got shape (%d, %d)", bounds.length, bounds[0].length) + ); + } + + // Scale each dimension + for (int i = 0; i < actualN; i++) { + for (int j = 0; j < d; j++) { + samples[i][j] = samples[i][j] * (bounds[j][1] - bounds[j][0]) + bounds[j][0]; + } + } + } + + return new Array2DRowRealMatrix(samples); + } + + /** + * Overloaded method with default parameters + */ + public static RealMatrix sobolSequence(int n, int d) { + return sobolSequence(n, d, false, null, null, 0, true); + } + + /** + * Overloaded method with scramble and seed parameters + */ + public static RealMatrix sobolSequence(int n, int d, boolean scramble, Integer seed) { + return sobolSequence(n, d, scramble, seed, null, 0, true); + } + + /** + * Overloaded method with all parameters except usePowOf2 + */ + public static RealMatrix sobolSequence( + int n, int d, boolean scramble, Integer seed, double[][] bounds, int skip) { + return sobolSequence(n, d, scramble, seed, bounds, skip, true); + } + + /** + * Check if a number is a power of 2 + */ + private static boolean isPowerOfTwo(int n) { + return n > 0 && (n & (n - 1)) == 0; + } + + /** + * Find the next power of 2 greater than or equal to n + */ + private static int nextPowerOfTwo(int n) { + if (n <= 1) return 1; + int result = 1; + while (result < n) { + result <<= 1; + } + return result; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeSparseGrid.java b/src/main/java/com/jdoe/algorithms/DoeSparseGrid.java new file mode 100644 index 0000000..aec8e3e --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeSparseGrid.java @@ -0,0 +1,313 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +import java.util.*; +import java.util.stream.Collectors; + +/** + * Sparse Grid Design of Experiments + *

+ * This module implements sparse grid designs based on Smolyak's construction. + * Sparse grids provide efficient sampling of high-dimensional spaces with + * good space-filling properties while requiring significantly fewer points + * than full tensor product grids. + *

+ * This code was originally developed based on the MATLAB Sparse Grid + * Interpolation Toolbox by: + * Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart + * Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved. + * email: klimkeas@ians.uni-stuttgart.de + * website: https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/spinterp.html + */ +public class DoeSparseGrid { + + /** + * Generate a sparse grid design using Smolyak's construction. + * + * @param nLevel Sparse grid level. Higher levels provide more points. + * @param nFactors Number of factors/dimensions in the design space. + * @param gridType Type of 1D grid points to use. Valid values are "clenshaw_curtis", "chebyshev", "gauss_patterson". + * @return Design points in the unit hypercube [0, 1]^nFactors. + * @throws IllegalArgumentException if nLevel is negative or nFactors is not positive + */ + public static RealMatrix doeSparseGrid(int nLevel, int nFactors, String gridType) { + if (nLevel < 0) { + throw new IllegalArgumentException("nLevel must be non-negative"); + } + if (nFactors < 1) { + throw new IllegalArgumentException("nFactors must be positive"); + } + + if (gridType == null) { + gridType = "clenshaw_curtis"; + } + + // Generate sparse grid points + RealMatrix design = _generateSparseGridPoints(nLevel, nFactors, gridType); + + return design; + } + + /** + * Overloaded method with default grid type + */ + public static RealMatrix doeSparseGrid(int nLevel, int nFactors) { + return doeSparseGrid(nLevel, nFactors, "clenshaw_curtis"); + } + + /** + * Compute the number of points in a sparse grid design. + * + * @param nLevel Sparse grid level. + * @param nFactors Number of dimensions. + * @return Number of points in the sparse grid. + * @throws IllegalArgumentException if nLevel is negative or nFactors is not positive + */ + public static int sparseGridDimension(int nLevel, int nFactors) { + if (nLevel < 0) { + throw new IllegalArgumentException("nLevel must be non-negative"); + } + if (nFactors < 1) { + throw new IllegalArgumentException("nFactors must be positive"); + } + + return _spdimFormula(nLevel, nFactors); + } + + /** + * Sparse grid dimension formulas from MATLAB spinterp spdim function. + * Based on Schreiber (2000) polynomial formulas. + */ + private static int _spdimFormula(int n, int d) { + if (n == 0) { + return 1; + } else if (n == 1) { + return 2 * d + 1; + } else if (n == 2) { + return 2 * d * d + 2 * d + 1; + } else if (n == 3) { + return (int) Math.round((4 * Math.pow(d, 3) + 6 * d * d + 14 * d) / 3) + 1; + } else if (n == 4) { + return (int) Math.round((2 * Math.pow(d, 4) + 4 * Math.pow(d, 3) + 22 * d * d + 20 * d) / 3) + 1; + } else if (n == 5) { + return (int) Math.round((4 * Math.pow(d, 5) + 10 * Math.pow(d, 4) + 100 * Math.pow(d, 3) + 170 * d * d + 196 * d) / 15) + 1; + } else if (n == 6) { + return (int) Math.round( + (4 * Math.pow(d, 6) + 12 * Math.pow(d, 5) + 190 * Math.pow(d, 4) + 480 * Math.pow(d, 3) + + 1246 * d * d + 948 * d) / 45) + 1; + } else { + return (int) Math.round( + (8 * Math.pow(d, 7) + 28 * Math.pow(d, 6) + 644 * Math.pow(d, 5) + 2170 * Math.pow(d, 4) + + 9632 * Math.pow(d, 3) + 15442 * d * d + 12396 * d) / 315) + 1; + } + } + + /** + * Generate sparse grid points + */ + private static RealMatrix _generateSparseGridPoints(int nLevel, int nFactors, String gridType) { + int targetCount = _spdimFormula(nLevel, nFactors); + + if (nLevel == 0) { + double[][] result = new double[1][nFactors]; + for (int i = 0; i < nFactors; i++) { + result[0][i] = 0.5; + } + return new Array2DRowRealMatrix(result); + } + + List> points = new ArrayList<>(); + + // Center point + List centerPoint = new ArrayList<>(); + for (int i = 0; i < nFactors; i++) { + centerPoint.add(0.5); + } + points.add(centerPoint); + + // Level 1: axis points + if (nLevel >= 1) { + for (int dim = 0; dim < nFactors; dim++) { + for (double val : Arrays.asList(0.0, 1.0)) { + List point = new ArrayList<>(Collections.nCopies(nFactors, 0.5)); + point.set(dim, val); + points.add(point); + } + } + } + + // Level 2+: structured interior points + if (nLevel >= 2) { + int gridSize = Math.min(nLevel + 2, 7); + double[] coords = linspace(0, 1, gridSize); + + // Single-dimension variations + for (int dim = 0; dim < nFactors; dim++) { + for (double coord : coords) { + if (coord != 0.0 && coord != 0.5 && coord != 1.0) { + List point = new ArrayList<>(Collections.nCopies(nFactors, 0.5)); + point.set(dim, coord); + points.add(point); + } + } + } + + // Multi-dimensional combinations for higher levels + if (nLevel >= 3) { + // Create subset of interior points + List coordsSubset = new ArrayList<>(); + for (double coord : coords) { + if (coord != 0.0 && coord != 1.0) { + coordsSubset.add(coord); + } + } + + for (int r = 2; r < Math.min(nFactors + 1, 4); r++) { + List> dimCombinations = combinations(range(nFactors), r); + for (List dims : dimCombinations) { + List> valCombinations = cartesianProduct( + Collections.nCopies(r, coordsSubset.subList(0, Math.min(2, coordsSubset.size())))); + + for (List vals : valCombinations) { + List point = new ArrayList<>(Collections.nCopies(nFactors, 0.5)); + for (int i = 0; i < dims.size(); i++) { + point.set(dims.get(i), vals.get(i)); + } + points.add(point); + + if (points.size() >= targetCount) { + break; + } + } + if (points.size() >= targetCount) { + break; + } + } + if (points.size() >= targetCount) { + break; + } + } + } + } + + // Remove duplicates and ensure exact count + List> uniquePoints = new ArrayList<>(); + Set seen = new HashSet<>(); + for (List point : points) { + StringBuilder sb = new StringBuilder(); + for (double x : point) { + sb.append(String.format("%.8f", x)).append(","); + } + String pointStr = sb.toString(); + if (!seen.contains(pointStr)) { + seen.add(pointStr); + uniquePoints.add(new ArrayList<>(point)); + } + } + + // Fill to exact target if needed + while (uniquePoints.size() < targetCount) { + int fillGridSize = targetCount / nFactors + 3; + double[] gridVals = linspace(0, 1, fillGridSize); + + // Fixed line: Convert double array to List properly + List> gridValCombinations = cartesianProduct( + Collections.nCopies(nFactors, Arrays.stream(gridVals).boxed().collect(Collectors.toList()))); + + for (List combo : gridValCombinations) { + StringBuilder sb = new StringBuilder(); + for (double x : combo) { + sb.append(String.format("%.8f", x)).append(","); + } + String pointStr = sb.toString(); + if (!seen.contains(pointStr)) { + seen.add(pointStr); + uniquePoints.add(new ArrayList<>(combo)); + if (uniquePoints.size() >= targetCount) { + break; + } + } + } + if (uniquePoints.size() >= targetCount) { + break; + } + } + + // Convert to matrix + double[][] result = new double[Math.min(uniquePoints.size(), targetCount)][nFactors]; + for (int i = 0; i < Math.min(uniquePoints.size(), targetCount); i++) { + for (int j = 0; j < nFactors; j++) { + result[i][j] = uniquePoints.get(i).get(j); + } + } + + return new Array2DRowRealMatrix(result); + } + + // Utility methods + + private static double[] linspace(double start, double end, int num) { + double[] result = new double[num]; + if (num == 1) { + result[0] = start; + return result; + } + double step = (end - start) / (num - 1); + for (int i = 0; i < num; i++) { + result[i] = start + i * step; + } + return result; + } + + private static List range(int n) { + List result = new ArrayList<>(); + for (int i = 0; i < n; i++) { + result.add(i); + } + return result; + } + + private static List> combinations(List elements, int r) { + List> result = new ArrayList<>(); + combinationsHelper(elements, r, 0, new ArrayList<>(), result); + return result; + } + + private static void combinationsHelper(List elements, int r, int start, + List current, List> result) { + if (current.size() == r) { + result.add(new ArrayList<>(current)); + return; + } + for (int i = start; i < elements.size(); i++) { + current.add(elements.get(i)); + combinationsHelper(elements, r, i + 1, current, result); + current.remove(current.size() - 1); + } + } + + private static List> cartesianProduct(List> lists) { + List> result = new ArrayList<>(); + if (lists.isEmpty()) { + result.add(new ArrayList<>()); + return result; + } + cartesianProductHelper(lists, 0, new ArrayList<>(), result); + return result; + } + + private static void cartesianProductHelper(List> lists, int depth, + List current, List> result) { + if (depth == lists.size()) { + result.add(new ArrayList<>(current)); + return; + } + for (T item : lists.get(depth)) { + current.add(item); + cartesianProductHelper(lists, depth + 1, current, result); + current.remove(current.size() - 1); + } + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeTaguchi.java b/src/main/java/com/jdoe/algorithms/DoeTaguchi.java new file mode 100644 index 0000000..58a75a6 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeTaguchi.java @@ -0,0 +1,182 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +import java.util.*; + +/** + * Inspired by Taguchi design methodology and orthogonal arrays developed by Genichi Taguchi, + * this code provides utilities for generating Taguchi arrays, building experimental designs, + * and computing Signal-to-Noise Ratios, based on orthogonal array libraries. + *

+ * Sources of orthogonal arrays: + * - University of York, Department of Mathematics: + * https://www.york.ac.uk/depts/maths/tables/orthogonal.htm + * - A Library of Orthogonal Arrays by N. J. A. Sloane: + * https://neilsloane.com/oadir/ + *

+ * References: + * - Taguchi G., Chowdhury S., Wu Y. (2005). "Taguchi's Quality Engineering Handbook." Wiley. + * - Montgomery D. C. (2017). "Design and Analysis of Experiments." Wiley. + * - [What are Taguchi designs?](https://www.itl.nist.gov/div898/handbook/pri/section5/pri56.htm) + */ +public class DoeTaguchi { + + /** + * Return a Taguchi orthogonal array by its descriptive name. + * + * @param oaName Name of the array, e.g., 'L4(2^3)', 'L8(2^7)', 'L9(3^4)', etc. + * @return The orthogonal array (zero-indexed factor levels) + * @throws IllegalArgumentException If the array name is not found + */ + public static RealMatrix getOrthogonalArray(String oaName) { + if (!ORTHOGONAL_ARRAYS.containsKey(oaName)) { + throw new IllegalArgumentException( + String.format("Orthogonal array '%s' not found. Available: %s", + oaName, ORTHOGONAL_ARRAYS.keySet()) + ); + } + return ORTHOGONAL_ARRAYS.get(oaName); + } + + /** + * List descriptive names of available Taguchi orthogonal arrays. + * + * @return List of array names, e.g., ['L4(2^3)', 'L8(2^7)', 'L9(3^4)', ...]. + */ + public static List listOrthogonalArrays() { + return new ArrayList<>(ORTHOGONAL_ARRAYS.keySet()); + } + + /** + * Generate a Taguchi design matrix using an orthogonal array and factor levels. + * + * @param oaName Name of Taguchi orthogonal array, e.g., 'L4(2^3)', 'L9(3^4)', etc. + * @param levelsPerFactor Each inner list defines actual levels/settings for each factor. + * Length must match number of columns in the orthogonal array. + * @return Design matrix with actual factor settings (not coded levels). + * @throws IllegalArgumentException If number of levels does not match number of factors. + */ + public static Object[][] taguchiDesign(String oaName, List> levelsPerFactor) { + RealMatrix array = getOrthogonalArray(oaName); + int nFactors = array.getColumnDimension(); + + if (levelsPerFactor.size() != nFactors) { + throw new IllegalArgumentException( + String.format("Number of factors in array (%d) does not match " + + "number of levels_per_factor provided (%d).", + nFactors, levelsPerFactor.size()) + ); + } + + Object[][] designMatrix = new Object[array.getRowDimension()][nFactors]; + + for (int i = 0; i < nFactors; i++) { + List levels = levelsPerFactor.get(i); + for (int j = 0; j < array.getRowDimension(); j++) { + int level = (int) array.getEntry(j, i); + designMatrix[j][i] = levels.get(level); + } + } + + return designMatrix; + } + + /** + * Calculate the Signal-to-Noise Ratio (SNR) for Taguchi designs. + * + * @param responses Repeated measurements for a single trial (1D array). + * @param objective Optimization goal, one of: LARGER, SMALLER, NOMINAL. + * @return SNR value in decibels (dB). + * @throws IllegalArgumentException If the objective is not recognized. + */ + public static double computeSnr(double[] responses, TaguchiObjective objective) { + if (objective == TaguchiObjective.LARGER_IS_BETTER) { + double sum = 0.0; + for (double resp : responses) { + sum += 1.0 / (resp * resp); + } + double mean = sum / responses.length; + return -10 * Math.log10(mean); + } else if (objective == TaguchiObjective.SMALLER_IS_BETTER) { + double sum = 0.0; + for (double resp : responses) { + sum += resp * resp; + } + double mean = sum / responses.length; + return -10 * Math.log10(mean); + } else if (objective == TaguchiObjective.NOMINAL_IS_BEST) { + double meanY = Arrays.stream(responses).average().orElse(0.0); + double sumSquaredDiffs = Arrays.stream(responses) + .map(x -> (x - meanY) * (x - meanY)) + .sum(); + double stdY = Math.sqrt(sumSquaredDiffs / (responses.length - 1)); + return 10 * Math.log10((meanY * meanY) / (stdY * stdY)); + } else { + throw new IllegalArgumentException("Invalid objective specified."); + } + } + + /** + * Overloaded method with default objective + */ + public static double computeSnr(double[] responses) { + return computeSnr(responses, TaguchiObjective.LARGER_IS_BETTER); + } + + /** + * Enumeration for Taguchi optimization objectives when calculating SNR. + */ + public enum TaguchiObjective { + LARGER_IS_BETTER("larger is better"), + SMALLER_IS_BETTER("smaller is better"), + NOMINAL_IS_BEST("nominal is best"); + + private final String description; + + TaguchiObjective(String description) { + this.description = description; + } + + public String getDescription() { + return description; + } + } + + // Mock orthogonal arrays - In a real implementation, these would come from external library + private static final Map ORTHOGONAL_ARRAYS = new HashMap<>(); + + static { + // Initialize with some example orthogonal arrays + ORTHOGONAL_ARRAYS.put("L4(2^3)", new Array2DRowRealMatrix(new double[][]{ + {0, 0, 0}, + {0, 1, 1}, + {1, 0, 1}, + {1, 1, 0} + })); + + ORTHOGONAL_ARRAYS.put("L8(2^7)", new Array2DRowRealMatrix(new double[][]{ + {0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 1, 1, 1, 1}, + {0, 1, 1, 0, 0, 1, 1}, + {0, 1, 1, 1, 1, 0, 0}, + {1, 0, 1, 0, 1, 0, 1}, + {1, 0, 1, 1, 0, 1, 0}, + {1, 1, 0, 0, 1, 1, 0}, + {1, 1, 0, 1, 0, 0, 1} + })); + + ORTHOGONAL_ARRAYS.put("L9(3^4)", new Array2DRowRealMatrix(new double[][]{ + {0, 0, 0, 0}, + {0, 1, 2, 1}, + {0, 2, 1, 2}, + {1, 0, 2, 2}, + {1, 1, 1, 0}, + {1, 2, 0, 1}, + {2, 0, 1, 1}, + {2, 1, 0, 2}, + {2, 2, 2, 0} + })); + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeUnion.java b/src/main/java/com/jdoe/algorithms/DoeUnion.java new file mode 100644 index 0000000..fd3220f --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeUnion.java @@ -0,0 +1,58 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; + +/** + * Join two matrices by stacking them on top of each other. + *

+ * This code was originally published by the following individuals for use with + * Scilab: + * Copyright (C) 2012 - 2013 - Michael Baudin + * Copyright (C) 2012 - Maria Christopoulou + * Copyright (C) 2010 - 2011 - INRIA - Michael Baudin + * Copyright (C) 2009 - Yann Collette + * Copyright (C) 2009 - CEA - Jean-Marc Martinez + *

+ * Much thanks goes to these individuals. It has been converted to Python by + * Abraham Lee. + */ +public class DoeUnion { + + /** + * Join two matrices by stacking them on top of each other. + * + * @param H1 The matrix that goes on top of the new matrix + * @param H2 The matrix that goes on bottom of the new matrix + * @return The new matrix that contains the rows of H1 on top of the rows of H2 + */ + public static RealMatrix matrixUnion(RealMatrix H1, RealMatrix H2) { + int rows1 = H1.getRowDimension(); + int rows2 = H2.getRowDimension(); + int cols = H1.getColumnDimension(); + + // Check that both matrices have the same number of columns + if (cols != H2.getColumnDimension()) { + throw new IllegalArgumentException("Both matrices must have the same number of columns"); + } + + // Create a new matrix with the combined rows + RealMatrix result = new Array2DRowRealMatrix(rows1 + rows2, cols); + + // Copy data from H1 + for (int i = 0; i < rows1; i++) { + for (int j = 0; j < cols; j++) { + result.setEntry(i, j, H1.getEntry(i, j)); + } + } + + // Copy data from H2 + for (int i = 0; i < rows2; i++) { + for (int j = 0; j < cols; j++) { + result.setEntry(rows1 + i, j, H2.getEntry(i, j)); + } + } + + return result; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoeVanillaMorris.java b/src/main/java/com/jdoe/algorithms/DoeVanillaMorris.java new file mode 100644 index 0000000..be702b3 --- /dev/null +++ b/src/main/java/com/jdoe/algorithms/DoeVanillaMorris.java @@ -0,0 +1,210 @@ +package com.doe.algorithms; + +import org.apache.commons.math3.linear.*; +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; + +/** + * Generate samples using the Morris Method (Vanilla, no optimization). + *

+ * This module implements the original (unoptimized) Morris method for + * global sensitivity analysis, a computationally efficient screening + * technique that estimates the importance of input variables by analyzing + * one-at-a-time (OAT) trajectories through a discretized parameter space. + * It is especially useful for identifying influential parameters in + * high-dimensional models with relatively low computational cost. + *

+ * References: + * 1. Morris, M.D., 1991. Factorial Sampling Plans for Preliminary Computational Experiments. + * Technometrics 33, 161-174. + * 2. Campolongo, F., Cariboni, J., & Saltelli, A., 2007. + * An effective screening design for sensitivity analysis of large models. + * Environmental Modelling & Software, 22(10), 1509-1518. + * 3. Ruano, M.V., Ribes, J., Seco, A., Ferrer, J., 2012. + * An improved sampling strategy based on trajectory design for application + * of the Morris method to systems with many input factors. + * Environmental Modelling & Software 37, 103-109. + */ +public class DoeVanillaMorris { + + /** + * Generate samples using the Morris Method (Vanilla, no optimization). + * + * @param numVars Number of input variables (i.e., the dimensionality of the problem) + * @param N Number of trajectories to generate + * @param numLevels Number of levels in the grid (must be even). Default is 4 + * @param seed Random seed for reproducibility + * @return Matrix of shape (N * (num_vars + 1), num_vars) with Morris samples + * @throws IllegalArgumentException if numLevels is not even + */ + public static RealMatrix morrisSampling(int numVars, int N, int numLevels, Integer seed) { + if (numLevels % 2 != 0) { + throw new IllegalArgumentException("numLevels must be an even number"); + } + + RandomGenerator rng; + if (seed != null) { + rng = new JDKRandomGenerator(seed); + } else { + rng = new JDKRandomGenerator(); + } + + int D = numVars; + + double delta = (double) numLevels / (2 * (numLevels - 1)); + double[] G = linspace(0, 1 - delta, numLevels / 2); + + List samples = new ArrayList<>(); + + for (int traj = 0; traj < N; traj++) { + // Starting point x* on the grid + double[] xStar = new double[D]; + for (int i = 0; i < D; i++) { + xStar[i] = G[rng.nextInt(G.length)]; + } + + + // Diagonal matrix of directions (+1 or -1) + RealMatrix dStar = createDiagonalMatrix(D, () -> rng.nextBoolean() ? 1.0 : -1.0); + // Lower-triangular B matrix + RealMatrix B = createLowerTriangularMatrix(D + 1, D); + + // Random permutation matrix P* + RealMatrix pStar = createRandomPermutationMatrix(D, rng); + + // J: ones matrix + RealMatrix J = new Array2DRowRealMatrix(D + 1, D); + for (int i = 0; i < D + 1; i++) { + for (int j = 0; j < D; j++) { + J.setEntry(i, j, 1.0); + } + } + + // Construct B* (trajectory matrix) + // B_star = x_star + delta / 2 * ((2 * B @ P_star - J) @ D_star + J) + RealMatrix BPStar = B.multiply(pStar); + RealMatrix twoBPStar = BPStar.scalarMultiply(2); + RealMatrix twoBPStarMinusJ = twoBPStar.subtract(J); + RealMatrix result = twoBPStarMinusJ.multiply(dStar); + result = result.add(J); + result = result.scalarMultiply(delta / 2); + + // Add xStar to each row of the result + for (int i = 0; i < result.getRowDimension(); i++) { + for (int j = 0; j < result.getColumnDimension(); j++) { + result.setEntry(i, j, result.getEntry(i, j) + xStar[j]); + } + } + + samples.add(result); + } + + // Combine all trajectory matrices + if (samples.isEmpty()) { + return new Array2DRowRealMatrix(0, D); + } + + int totalRows = 0; + for (RealMatrix matrix : samples) { + totalRows += matrix.getRowDimension(); + } + + RealMatrix finalResult = new Array2DRowRealMatrix(totalRows, D); + int currentRow = 0; + for (RealMatrix matrix : samples) { + int rows = matrix.getRowDimension(); + for (int i = 0; i < rows; i++) { + for (int j = 0; j < D; j++) { + finalResult.setEntry(currentRow + i, j, matrix.getEntry(i, j)); + } + } + currentRow += rows; + } + + return finalResult; + } + + /** + * Overloaded method with default numLevels (4) + */ + public static RealMatrix morrisSampling(int numVars, int N, Integer seed) { + return morrisSampling(numVars, N, 4, seed); + } + + /** + * Overloaded method with default numLevels (4) and seed (null) + */ + public static RealMatrix morrisSampling(int numVars, int N) { + return morrisSampling(numVars, N, 4, null); + } + + /** + * Create a diagonal matrix with values generated by the provided supplier + */ + private static RealMatrix createDiagonalMatrix(int size, java.util.function.Supplier valueSupplier) { + RealMatrix matrix = new Array2DRowRealMatrix(size, size); + for (int i = 0; i < size; i++) { + matrix.setEntry(i, i, valueSupplier.get()); + } + return matrix; + } + + /** + * Create a lower-triangular matrix + */ + private static RealMatrix createLowerTriangularMatrix(int rows, int cols) { + RealMatrix matrix = new Array2DRowRealMatrix(rows, cols); + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + if (i > j) { + matrix.setEntry(i, j, 1.0); + } + } + } + return matrix; + } + + /** + * Create a random permutation matrix + */ + private static RealMatrix createRandomPermutationMatrix(int size, RandomGenerator rng) { + RealMatrix matrix = new Array2DRowRealMatrix(size, size); + for (int i = 0; i < size; i++) { + matrix.setEntry(i, i, 1.0); + } + + // Perform random swaps to create a permutation matrix + for (int i = size - 1; i > 0; i--) { + int j = rng.nextInt(i + 1); + if (i != j) { + // Swap rows i and j + for (int k = 0; k < size; k++) { + double temp = matrix.getEntry(i, k); + matrix.setEntry(i, k, matrix.getEntry(j, k)); + matrix.setEntry(j, k, temp); + } + } + } + + return matrix; + } + + /** + * Create a linearly spaced array + */ + private static double[] linspace(double start, double end, int num) { + if (num == 1) { + return new double[]{start}; + } + double[] result = new double[num]; + double step = (end - start) / (num - 1); + for (int i = 0; i < num; i++) { + result[i] = start + i * step; + } + return result; + } +} diff --git a/src/main/java/com/jdoe/algorithms/DoehlertDOE.java b/src/main/java/com/jdoe/algorithms/DoehlertDOE.java index 5b0488b..586ea96 100644 --- a/src/main/java/com/jdoe/algorithms/DoehlertDOE.java +++ b/src/main/java/com/jdoe/algorithms/DoehlertDOE.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/algorithms/FactorialDOE.java b/src/main/java/com/jdoe/algorithms/FactorialDOE.java index e5d8c2f..e4fd96f 100644 --- a/src/main/java/com/jdoe/algorithms/FactorialDOE.java +++ b/src/main/java/com/jdoe/algorithms/FactorialDOE.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import java.util.*; import java.util.stream.Collectors; @@ -11,7 +11,7 @@ import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; -import com.jdoe.util.FactorialUtility; +import com.doe.util.FactorialUtility; /* * Copyright (c) 2025 Noor Mustafa diff --git a/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java b/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java index 15e8244..83779bc 100644 --- a/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java +++ b/src/main/java/com/jdoe/algorithms/RepeatCenterDOE.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/algorithms/StarDOE.java b/src/main/java/com/jdoe/algorithms/StarDOE.java index f150d01..5343e39 100644 --- a/src/main/java/com/jdoe/algorithms/StarDOE.java +++ b/src/main/java/com/jdoe/algorithms/StarDOE.java @@ -1,6 +1,6 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; -import com.jdoe.util.GenericDOEUtil; +import com.doe.util.GenericDOEUtil; public class StarDOE { /** diff --git a/src/main/java/com/jdoe/algorithms/UnionDOE.java b/src/main/java/com/jdoe/algorithms/UnionDOE.java index 576edc0..c5b0cd2 100644 --- a/src/main/java/com/jdoe/algorithms/UnionDOE.java +++ b/src/main/java/com/jdoe/algorithms/UnionDOE.java @@ -1,6 +1,6 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; -import com.jdoe.util.GenericDOEUtil; +import com.doe.util.GenericDOEUtil; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java b/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java index 45a7360..f2b7dd2 100644 --- a/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java +++ b/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java @@ -1,4 +1,4 @@ -package com.jdoe.util; +package com.doe.util; import java.util.ArrayList; import java.util.List; diff --git a/src/main/java/com/jdoe/util/FactorialUtility.java b/src/main/java/com/jdoe/util/FactorialUtility.java index 484233c..e565eb2 100644 --- a/src/main/java/com/jdoe/util/FactorialUtility.java +++ b/src/main/java/com/jdoe/util/FactorialUtility.java @@ -1,4 +1,4 @@ -package com.jdoe.util; +package com.doe.util; import java.util.ArrayList; import java.util.List; diff --git a/src/main/java/com/jdoe/util/GenericDOEUtil.java b/src/main/java/com/jdoe/util/GenericDOEUtil.java index fc32c95..8f56321 100644 --- a/src/main/java/com/jdoe/util/GenericDOEUtil.java +++ b/src/main/java/com/jdoe/util/GenericDOEUtil.java @@ -1,4 +1,4 @@ -package com.jdoe.util; +package com.doe.util; import java.util.Arrays; diff --git a/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java b/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java index c78f241..a96f057 100644 --- a/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java +++ b/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import static org.junit.Assert.*; import org.apache.commons.math3.linear.RealMatrix; diff --git a/src/test/java/com/jdoe/algorithms/CompositeDOETest.java b/src/test/java/com/jdoe/algorithms/CompositeDOETest.java index 89959a9..555f2f4 100644 --- a/src/test/java/com/jdoe/algorithms/CompositeDOETest.java +++ b/src/test/java/com/jdoe/algorithms/CompositeDOETest.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.linear.RealMatrix; import org.junit.Test; diff --git a/src/test/java/com/jdoe/algorithms/FactorialDOETest.java b/src/test/java/com/jdoe/algorithms/FactorialDOETest.java index c17c3d2..f2680f9 100644 --- a/src/test/java/com/jdoe/algorithms/FactorialDOETest.java +++ b/src/test/java/com/jdoe/algorithms/FactorialDOETest.java @@ -1,4 +1,4 @@ -package com.jdoe.algorithms; +package com.doe.algorithms; import org.apache.commons.math3.exception.NotStrictlyPositiveException; import org.apache.commons.math3.linear.RealMatrix; diff --git a/target/classes/log4j2.xml b/target/classes/log4j2.xml deleted file mode 100644 index 4cbeecc..0000000 --- a/target/classes/log4j2.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - - - - - - - - - - diff --git a/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst b/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst deleted file mode 100644 index c9315c0..0000000 --- a/target/maven-status/maven-compiler-plugin/compile/default-compile/createdFiles.lst +++ /dev/null @@ -1,6 +0,0 @@ -com/jdoe/util/FactorialUtility.class -com/jdoe/util/BuildRegressionMatrixUtility.class -com/jdoe/algorithms/BuildRegressionMatrix.class -com/jdoe/algorithms/BoxBehnkenDOE.class -com/jdoe/algorithms/FactorialDOE.class -com/jdoe/Testing.class diff --git a/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst b/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst deleted file mode 100644 index b2926fa..0000000 --- a/target/maven-status/maven-compiler-plugin/compile/default-compile/inputFiles.lst +++ /dev/null @@ -1,6 +0,0 @@ -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/Testing.java -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/BoxBehnkenDOE.java -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/BuildRegressionMatrix.java -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/algorithms/FactorialDOE.java -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/util/BuildRegressionMatrixUtility.java -/home/noor/IdeaProjects/JDOE/src/main/java/com/jdoe/util/FactorialUtility.java diff --git a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst deleted file mode 100644 index 8757472..0000000 --- a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/createdFiles.lst +++ /dev/null @@ -1,2 +0,0 @@ -com/jdoe/algorithms/FactorialDOETest.class -com/jdoe/algorithms/BoxBehnkenDOETest.class diff --git a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst b/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst deleted file mode 100644 index 91aac98..0000000 --- a/target/maven-status/maven-compiler-plugin/testCompile/default-testCompile/inputFiles.lst +++ /dev/null @@ -1,2 +0,0 @@ -/home/noor/IdeaProjects/JDOE/src/test/java/com/jdoe/algorithms/BoxBehnkenDOETest.java -/home/noor/IdeaProjects/JDOE/src/test/java/com/jdoe/algorithms/FactorialDOETest.java diff --git a/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml b/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml deleted file mode 100644 index 8c02771..0000000 --- a/target/surefire-reports/TEST-com.jdoe.algorithms.FactorialDOETest.xml +++ /dev/null @@ -1,81 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt b/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt deleted file mode 100644 index 2aad563..0000000 --- a/target/surefire-reports/com.jdoe.algorithms.FactorialDOETest.txt +++ /dev/null @@ -1,4 +0,0 @@ -------------------------------------------------------------------------------- -Test set: com.jdoe.algorithms.FactorialDOETest -------------------------------------------------------------------------------- -Tests run: 18, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.233 s -- in com.jdoe.algorithms.FactorialDOETest