+ * $ 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 + *
+ * $ 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/CranleyPattersonShiftDOE.java b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java
new file mode 100644
index 0000000..4b11dac
--- /dev/null
+++ b/src/main/java/com/jdoe/algorithms/CranleyPattersonShiftDOE.java
@@ -0,0 +1,66 @@
+package com.doe.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..34e3b09
--- /dev/null
+++ b/src/main/java/com/jdoe/algorithms/DoeFold.java
@@ -0,0 +1,95 @@
+package com.doe.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
+ * 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
+ * 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
+ * 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