Claude/translate iboss matlab to r o3 p4 o

Results/IBOSS_10M_NoK.csv

@@ -0,0 +1,51 @@
+Problem,n,m,Time,N,Obj
+Generator_10M_1_1.mat,10,10000000,4.364975,20,118.035786923113
+Generator_10M_1_10.mat,10,10000000,4.415632,20,121.831054336905
+Generator_10M_1_2.mat,10,10000000,4.505327,20,110.802208841728
+Generator_10M_1_3.mat,10,10000000,4.228981,20,97.4913487020384
+Generator_10M_1_4.mat,10,10000000,4.705958,20,82.1043845416311
+Generator_10M_1_5.mat,10,10000000,4.574209,20,106.555204708658
+Generator_10M_1_6.mat,10,10000000,4.610081,20,108.86236072323
+Generator_10M_1_7.mat,10,10000000,3.92479,20,111.10427408281
+Generator_10M_1_8.mat,10,10000000,4.576194,20,93.6717505708583
+Generator_10M_1_9.mat,10,10000000,4.181577,20,110.717917078181
+Generator_10M_2_1.mat,20,10000000,8.750429,40,233.649143646439
+Generator_10M_2_10.mat,20,10000000,8.158203,40,222.825589452468
+Generator_10M_2_2.mat,20,10000000,8.807856,40,247.316723304714
+Generator_10M_2_3.mat,20,10000000,8.476802,40,254.457593902342
+Generator_10M_2_4.mat,20,10000000,10.161736,40,164.960401257846
+Generator_10M_2_5.mat,20,10000000,8.482211,40,243.461899439601
+Generator_10M_2_6.mat,20,10000000,8.212102,40,217.628554872989
+Generator_10M_2_7.mat,20,10000000,9.864593,40,165.400233367178
+Generator_10M_2_8.mat,20,10000000,8.57601,40,194.125927173195
+Generator_10M_2_9.mat,20,10000000,9.808802,40,107.201125654082
+Generator_10M_3_1.mat,30,10000000,15.302226,60,306.816219066008
+Generator_10M_3_10.mat,30,10000000,15.242893,60,317.166229884767
+Generator_10M_3_2.mat,30,10000000,12.867234,60,355.865285326097
+Generator_10M_3_3.mat,30,10000000,11.732713,60,319.497459116525
+Generator_10M_3_4.mat,30,10000000,13.125791,60,281.424307835681
+Generator_10M_3_5.mat,30,10000000,12.649949,60,351.491558145388
+Generator_10M_3_6.mat,30,10000000,14.279315,60,404.277708765693
+Generator_10M_3_7.mat,30,10000000,14.9663,60,246.22990820306
+Generator_10M_3_8.mat,30,10000000,12.522791,60,375.02442093901
+Generator_10M_3_9.mat,30,10000000,12.813444,60,297.784534134187
+Generator_10M_4_1.mat,40,10000000,19.456385,80,444.850494589041
+Generator_10M_4_10.mat,40,10000000,15.861178,80,394.727145110651
+Generator_10M_4_2.mat,40,10000000,20.12846,80,389.352127798575
+Generator_10M_4_3.mat,40,10000000,16.918422,80,485.463994752531
+Generator_10M_4_4.mat,40,10000000,20.478763,80,401.684716464614
+Generator_10M_4_5.mat,40,10000000,20.605866,80,361.279450955749
+Generator_10M_4_6.mat,40,10000000,18.132969,80,506.501663997589
+Generator_10M_4_7.mat,40,10000000,20.856812,80,306.851493561119
+Generator_10M_4_8.mat,40,10000000,17.976415,80,412.531378745288
+Generator_10M_4_9.mat,40,10000000,17.790051,80,478.146575574834
+Generator_10M_5_1.mat,50,10000000,23.612641,100,636.010051571282
+Generator_10M_5_10.mat,50,10000000,25.422309,100,400.458913087078
+Generator_10M_5_2.mat,50,10000000,21.757456,100,578.079879913931
+Generator_10M_5_3.mat,50,10000000,23.239555,100,625.266206915978
+Generator_10M_5_4.mat,50,10000000,21.219201,100,542.017150596086
+Generator_10M_5_5.mat,50,10000000,22.66709,100,553.879731314336
+Generator_10M_5_6.mat,50,10000000,22.476265,100,514.870391508125
+Generator_10M_5_7.mat,50,10000000,23.882625,100,642.125141185357
+Generator_10M_5_8.mat,50,10000000,22.921065,100,539.3704097554
+Generator_10M_5_9.mat,50,10000000,22.515005,100,578.920599348739

iboss_Exact_D_Optimal_NoK.R

@@ -0,0 +1,97 @@
+# A Column Generation Approach to Exact Experimental Design
+# Authors: Selin Ahipasaoglu, Stefano Cipolla, Jacek Gondzio
+# arXiv ID: 2507.03210 (July 2025)
+#
+# R translation of matlab/iboss_Exact_D_Optimal_NoK.m
+# Uses the original R IBOSS package (iboss.od with compiled C++ backend)
+#
+# Usage:
+#   From the project root directory:
+#     Rscript iboss_Exact_D_Optimal_NoK.R
+#   Or within an R session:
+#     source("iboss_Exact_D_Optimal_NoK.R")
+#
+# Prerequisites:
+#   install.packages("R.matlab")  # for reading .mat files
+#   devtools::install(".")        # install the IBOSS package from source
+
+# --- Load required packages ---
+library(R.matlab)
+
+# Load the IBOSS package.
+# Option 1: If installed as a package
+# library(IBOSS)
+# Option 2: If running from the IBOSS project root (development mode)
+devtools::load_all(".")
+
+# --- Configuration ---
+# The path on which all the data are:
+problems_path <- "/scratch/sc9c23/Dataset/Gen1"
+
+# Find all the problems and store their names
+files <- sort(list.files(problems_path, pattern = "^Generator_10M_.*\\.mat$",
+                         full.names = FALSE))
+
+# N offsets
+NN <- c(0)
+
+# Total number of experiments
+expN <- length(files) * length(NN)
+
+# --- Initialize results data frame ---
+resultsIBOSS <- data.frame(
+  Problem = character(expN),
+  n       = numeric(expN),
+  m       = numeric(expN),
+  Time    = numeric(expN),
+  N       = numeric(expN),
+  Obj     = numeric(expN),
+  stringsAsFactors = FALSE
+)
+
+# --- Main experiment loop ---
+exp_num <- 0
+
+for (k in seq_along(files)) {
+  # Load .mat file and extract the data matrix
+  mat_data <- readMat(file.path(problems_path, files[k]))
+  X <- mat_data$genmatrix
+
+  cat(files[k], "\n")
+
+  # X is (n x m): n = number of features (rows), m = number of observations (columns)
+  n <- nrow(X)
+  m <- ncol(X)
+
+  for (ll in seq_along(NN)) {
+    cat(ll, "\n")
+
+    N <- 2 * n + NN[ll]
+    exp_num <- exp_num + 1
+
+    # Run IBOSS: transpose X so that iboss.od receives (observations x features)
+    # Y = ones(m,1) is a dummy response (only index selection matters here)
+    time_start <- proc.time()
+    fit <- iboss.od(t(X), rep(1, m), N)
+    time_IBOSS <- (proc.time() - time_start)[["elapsed"]]
+
+    # Select the columns of X corresponding to the chosen observation indices
+    X_S <- X[, fit$index, drop = FALSE]
+
+    # D-optimality objective: log(det(X_S %*% t(X_S)))
+    obj_IBOSS <- log(det(X_S %*% t(X_S)))
+
+    # Record results
+    resultsIBOSS$Problem[exp_num] <- files[k]
+    resultsIBOSS$n[exp_num]       <- n
+    resultsIBOSS$m[exp_num]       <- m
+    resultsIBOSS$Time[exp_num]    <- time_IBOSS
+    resultsIBOSS$N[exp_num]       <- N
+    resultsIBOSS$Obj[exp_num]     <- obj_IBOSS
+  }
+}
+
+# --- Write results to CSV ---
+write.csv(resultsIBOSS, file.path("Results", "IBOSS_10M_NoK.csv"), row.names = FALSE)
+
+cat("Results written to Results/IBOSS_10M_NoK.csv\n")

matlab/compile_mex.m

@@ -0,0 +1,14 @@
+% COMPILE_MEX  Compile the MEX files for iboss_od_mex.
+%
+%   Run this script once before using iboss_od_mex.
+%   The compiled binaries will be placed in the current directory.
+
+fprintf('Compiling getIdx_mex.cpp ...\n');
+mex(fullfile(fileparts(mfilename('fullpath')), 'getIdx_mex.cpp'), ...
+    '-outdir', fileparts(mfilename('fullpath')));
+
+fprintf('Compiling getIdxR_mex.cpp ...\n');
+mex(fullfile(fileparts(mfilename('fullpath')), 'getIdxR_mex.cpp'), ...
+    '-outdir', fileparts(mfilename('fullpath')));
+
+fprintf('Done. MEX files ready.\n');

matlab/getIdxR_mex.cpp

@@ -0,0 +1,114 @@
+/*
+ * getIdxR_mex.cpp  –  MATLAB MEX port of getIdxR() from src/getIdx.cpp
+ *
+ * Find the indices of the r smallest and r largest values in a vector,
+ * excluding a set of already-selected indices.
+ * Uses std::nth_element for O(n) selection (same algorithm as the R/C++ version).
+ *
+ * Usage (from MATLAB, after mex-compiling):
+ *   idx = getIdxR_mex(r, z, del)
+ *
+ *   r   – scalar integer, number of extremes per tail
+ *   z   – double column vector (n x 1)
+ *   del – double vector of 1-based indices to exclude
+ *
+ * Returns:
+ *   idx – (2*r x 1) double vector of 1-based indices
+ *         [r smallest indices; r largest indices]
+ */
+
+#include "mex.h"
+#include <algorithm>
+
+void mexFunction(int nlhs, mxArray *plhs[],
+                 int nrhs, const mxArray *prhs[])
+{
+    /* ---- input validation ---- */
+    if (nrhs != 3)
+        mexErrMsgIdAndTxt("getIdxR:nrhs",
+                          "Three inputs required: r, z, del.");
+    if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]))
+        mexErrMsgIdAndTxt("getIdxR:notDouble", "z must be a real double vector.");
+
+    int r = (int)mxGetScalar(prhs[0]);
+    double *z = mxGetPr(prhs[1]);
+    int n = (int)mxGetNumberOfElements(prhs[1]);
+    double *del_dbl = mxGetPr(prhs[2]);
+    int m = (int)mxGetNumberOfElements(prhs[2]);
+
+    /* ---- convert del to int and sort ---- */
+    int *del = new int[m];
+    for (int i = 0; i < m; i++)
+        del[i] = (int)del_dbl[i];
+    std::sort(del, del + m);
+
+    /* ---- build filtered values (exclude del indices) ---- */
+    double *y = new double[n - m];
+    int j = 0, k = 0;
+    for (int i = 0; i < n; i++) {
+        if (j >= m) {
+            y[k++] = z[i];
+        } else if (del[j] != i + 1) {
+            y[k++] = z[i];
+        } else {
+            j++;
+        }
+    }
+
+    /* ---- find r-th smallest value (lower threshold) ---- */
+    std::nth_element(y, y + r - 1, y + n - m);
+    double yrl = y[r - 1];
+
+    /* ---- find r-th largest value (upper threshold) ---- */
+    j = 0; k = 0;
+    for (int i = 0; i < n; i++) {
+        if (j >= m) {
+            y[k++] = -z[i];
+        } else if (del[j] != i + 1) {
+            y[k++] = -z[i];
+        } else {
+            j++;
+        }
+    }
+    std::nth_element(y, y + r - 1, y + n - m);
+    double yru = -y[r - 1];
+
+    delete[] y;
+
+    /* ---- collect indices (excluding del) ---- */
+    int *locl = new int[r];
+    int *locu = new int[r];
+    int jl = 0, ju = 0;
+    j = 0;
+
+    for (int i = 0; i < n; i++) {
+        if (j >= m) {
+            if (z[i] <= yrl && jl < r)
+                locl[jl++] = i + 1;
+            if (z[i] >= yru && ju < r)
+                locu[ju++] = i + 1;
+        } else if (del[j] != i + 1) {
+            if (z[i] <= yrl && jl < r)
+                locl[jl++] = i + 1;
+            if (z[i] >= yru && ju < r)
+                locu[ju++] = i + 1;
+        } else {
+            j++;
+        }
+        if (jl >= r && ju >= r)
+            break;
+    }
+
+    /* ---- build output ---- */
+    plhs[0] = mxCreateDoubleMatrix(2 * r, 1, mxREAL);
+    double *idx = mxGetPr(plhs[0]);
+
+    for (int i = 0; i < r; i++) {
+        idx[i]     = (double)locl[i];
+        idx[r + i] = (double)locu[i];
+    }
+
+    delete[] locl;
+    delete[] locu;
+    delete[] del;
+}

matlab/getIdx_mex.cpp

@@ -0,0 +1,81 @@
+/*
+ * getIdx_mex.cpp  –  MATLAB MEX port of getIdx() from src/getIdx.cpp
+ *
+ * Find the indices of the r smallest and r largest values in a vector.
+ * Uses std::nth_element for O(n) selection (same algorithm as the R/C++ version).
+ *
+ * Usage (from MATLAB, after mex-compiling):
+ *   idx = getIdx_mex(r, z)
+ *
+ *   r  – scalar integer, number of extremes per tail
+ *   z  – double column vector (n x 1)
+ *
+ * Returns:
+ *   idx – (2*r x 1) double vector of 1-based indices
+ *         [r smallest indices; r largest indices]
+ */
+
+#include "mex.h"
+#include <algorithm>
+
+void mexFunction(int nlhs, mxArray *plhs[],
+                 int nrhs, const mxArray *prhs[])
+{
+    /* ---- input validation ---- */
+    if (nrhs != 2)
+        mexErrMsgIdAndTxt("getIdx:nrhs",
+                          "Two inputs required: r (integer scalar), z (double vector).");
+    if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]))
+        mexErrMsgIdAndTxt("getIdx:notDouble", "z must be a real double vector.");
+
+    int r = (int)mxGetScalar(prhs[0]);
+    double *z = mxGetPr(prhs[1]);
+    int n = (int)mxGetNumberOfElements(prhs[1]);
+
+    if (r < 1 || 2 * r > n)
+        mexErrMsgIdAndTxt("getIdx:badR",
+                          "r must satisfy 1 <= r <= n/2.");
+
+    /* ---- find r-th smallest value (lower threshold) ---- */
+    double *y = new double[n];
+    for (int i = 0; i < n; i++)
+        y[i] = z[i];
+
+    std::nth_element(y, y + r - 1, y + n);
+    double yrl = y[r - 1];
+
+    /* ---- find r-th largest value (upper threshold) ---- */
+    for (int i = 0; i < n; i++)
+        y[i] = -z[i];
+
+    std::nth_element(y, y + r - 1, y + n);
+    double yru = -y[r - 1];
+
+    delete[] y;
+
+    /* ---- collect indices ---- */
+    int *locl = new int[r];
+    int *locu = new int[r];
+    int jl = 0, ju = 0;
+
+    for (int i = 0; i < n; i++) {
+        if (z[i] <= yrl && jl < r)
+            locl[jl++] = i + 1;          /* 1-based for MATLAB */
+        if (z[i] >= yru && ju < r)
+            locu[ju++] = i + 1;
+        if (jl >= r && ju >= r)
+            break;
+    }
+
+    /* ---- build output ---- */
+    plhs[0] = mxCreateDoubleMatrix(2 * r, 1, mxREAL);
+    double *idx = mxGetPr(plhs[0]);
+
+    for (int i = 0; i < r; i++) {
+        idx[i]     = (double)locl[i];
+        idx[r + i] = (double)locu[i];
+    }
+
+    delete[] locl;
+    delete[] locu;
+}