Merge HEAD on Github from AD.

DESCRIPTION

@@ -10,7 +10,7 @@ Authors@R: c(
 URL: https://github.com/dusadrian/IEEE
 BugReports: https://github.com/dusadrian/IEEE/issues
 Depends: R (>= 3.6.0)
-Imports: LogicOpt, Matrix, gurobi
+Imports: admisc, LogicOpt, lpSolve, Matrix, gurobi
 Description: An extensive set of functions to perform Qualitative Comparative Analysis:
   crisp sets ('csQCA'), temporal ('tQCA'), multi-value ('mvQCA')
   and fuzzy sets ('fsQCA'), using a GUI - graphical user interface.

NAMESPACE

@@ -1,7 +1,10 @@
 importFrom("gurobi", "gurobi")
+importFrom("lpSolve", "lp")
 importFrom("Matrix", "Matrix")
 importFrom("utils", "tail")
 importFrom("LogicOpt", "logicopt")
+importFrom("admisc", "recode", "export", "tryCatchWEM")
+export(compare)
 export(gendat)
 export(make_infile)
 export(solvechart)

R/compare.R

@@ -0,0 +1,90 @@
+compare <- function(
+    iter = 1:10,
+    ncols = 21:27,
+    nrows = 25:35,
+    file = "results.csv",
+    espresso = TRUE,
+    avg = 1
+) {
+    on.exit({
+        suppressWarnings(sink())
+    })
+
+    for (nc in ncols) {
+        for (nr in nrows) {
+            for (i in iter) {
+                dat <- gendat(i, nc, nr)
+                make_infile(dat)
+                ## $ sudo port install espresso
+                # lo <- system("espresso -Dexact -o eqntott infile.esp", intern = TRUE)[1]
+
+                if (espresso) {
+                    etc <- admisc::tryCatchWEM(
+                        es_time <- system.time(
+                            elo <- LogicOpt::logicopt(esp_file = "infile.esp")
+                        )
+                    )
+                }
+
+                ctc <- admisc::tryCatchWEM({
+                    cc_time <- NULL
+                    for (ii in seq(avg)) {
+                        cc_time <- c(
+                            cc_time,
+                            system.time(
+                                cc <- .Call("CCubes", dat, PACKAGE = "IEEE")
+                            )[3]
+                        )
+                    }
+                })
+
+                if (is.null(ctc)) {
+                    colnames(cc) <- colnames(dat)[-nc]
+                    if (espresso) {
+                        soles <- apply(
+                            elo[[1]][, -nc], 1,
+                            function(x) {
+                                x <- x[x != "-"]
+                                names(x)[x == 0] <- paste("!", names(x)[x == 0], sep = "")
+                                return(c(paste(names(x), collapse = "&"), length(x)))
+                            }
+                        )
+                    }
+
+                    solcc <- apply(cc, 1, function(x) {
+                        x <- x[x > 0] - 1
+                        names(x)[x == 0] <- paste("!", names(x)[x == 0], sep = "")
+                        return(c(paste(names(x), collapse = "&"), length(x)))
+                    })
+
+                    sink(file, append = TRUE)
+
+                    cat(
+                        sprintf(
+                        "%s, %s, %s, %s, %s, %s%s, %s%s\n",
+                        i, nc, nr,
+                        round(es_time[3], 3),
+                        round(min(cc_time), 3),
+                        ifelse(
+                            espresso,
+                            paste0("\"(", paste(soles[1, ], collapse = ") + ("), ")\", "),
+                            ""
+                        ),
+                        paste0("\"(", paste(solcc[1, ], collapse = ") + ("), ")\""),
+                        ifelse(
+                            espresso,
+                            paste(paste(soles[2, ], collapse = ""), ", ", sep = ""),
+                            ""
+                        ),
+                        paste(solcc[2, ], collapse = "")
+                        )
+                    )
+
+                    sink()
+                } else {
+                    cat(sprintf("i: %s, nc: %s, nr: %s\n", i, nc, nr))
+                }
+            }
+        }
+    }
+}

R/gendat.R

@@ -1,6 +1,6 @@
 
-gendat <- function(i, ncols, nrows) {
+gendat <- function(iter, ncols, nrows) {
-    set.seed(i * nrows * ncols)
+    set.seed(iter * nrows * ncols)
 
     # all configurations BUT the outcome have to be unique
     dat <- unique(

R/solvechart.R

@@ -48,7 +48,6 @@
                 value <- value[seq(PIs)]
             }
 
-            # cat("Saving the PI chart to a file\n")
             # admisc::export(
             #     as.data.frame(t(x)),
             #     file = paste("pic", sample(5000, 1), ".csv", sep = ""),

man/compare.Rd

@@ -0,0 +1,41 @@
+\name{compare}
+
+\alias{compare}
+
+\title{Compare CCubes results against Espresso}
+
+\description{
+Compare CCubes results against Espresso
+}
+
+\usage{
+compare(
+    iter = 1:10,
+    ncols = 21:27,
+    nrows = 25:35,
+    file = "results.csv",
+    espresso = TRUE,
+    avg = FALSE
+)
+}
+
+\arguments{
+    \item{iter}{Iterator}
+    \item{ncols}{Number of columns}
+    \item{nrows}{Number of rows}
+    \item{file}{Character, name / path to the file where results are written}
+    \item{espresso}{Boolean, also write Espresso results}
+    \item{avg}{Numeric, number of runs to average results over}
+}
+
+\author{
+Adrian Dusa
+}
+
+\examples{
+\dontrun{
+    compare(1, 7, 15)
+}
+}
+
+\keyword{functions}

man/gendat.Rd

@@ -0,0 +1,33 @@
+\name{gendat}
+
+\alias{gendat}
+
+\title{Generate a binary matrix}
+
+\description{
+Generate a random binary matrix for testing purposes.
+}
+
+\usage{
+gendat(iter, ncols, nrows)
+}
+
+\arguments{
+    \item{iter}{Iterator}
+    \item{ncols}{Number of columns}
+    \item{nrows}{Number of rows}
+}
+
+\value{
+A matrix containing binary values.
+}
+
+\author{
+Adrian Dusa
+}
+
+\examples{
+gendat(1, 7, 15)
+}
+
+\keyword{functions}

man/make_infile.Rd

@@ -0,0 +1,39 @@
+\name{make_infile}
+
+\alias{make_infile}
+
+\title{Create an infile for espresso}
+
+\description{
+Create a text file containing the data for the expresso minimizer.
+}
+
+\usage{
+make_infile(dat, espname = "infile.esp", copy = FALSE, ...)
+}
+
+\arguments{
+    \item{dat}{A binary matrix}
+    \item{espname}{Name of the text file}
+    \item{copy}{Logical, copy the file with a different name.}
+    \item{...}{Other arguments.}
+}
+
+\details{
+Different iterations generate different versions of the data.
+They are all temporarily saved into \dQuote{infile.esp} unless
+the argument \code{copy} is activated, which copies the infile
+into a file with a name containing the iterator, the number of
+columns and the number of rows.
+}
+
+\author{
+Adrian Dusa
+}
+
+\examples{
+dat <- gendat(1, 7, 15)
+make_infile(dat)
+}
+
+\keyword{functions}

man/solvechart.Rd

@@ -0,0 +1,43 @@
+\name{solvechart}
+
+\alias{solvechart}
+
+\title{Find the minimum for a prime implicants chart}
+
+\description{
+Solve a PI chart using a minimum number of prime implicants.
+}
+
+\usage{
+solvechart(x, ...)
+}
+
+\arguments{
+    \item{x}{A binary matrix.}
+    \item{...}{Other arguments.}
+}
+
+\details{
+A PI chart, in this package, is a binary matrix containing the prime implicants on the columns
+and the observed positive minterms on the rows.
+}
+
+\value{
+A vector of binary values
+}
+
+\author{
+Adrian Dusa
+}
+
+\examples{
+set.seed(12345)
+x <- matrix(
+    sample(0:1, prob = c(0.65, 0.35), size = 35, replace = TRUE),
+    nrow = 5
+)
+
+solvechart(x)
+}
+
+\keyword{functions}

src/CCubes.c

@@ -1,6 +1,6 @@
+#include <R.h>
 #include <R_ext/RS.h>
 #include <R_ext/Boolean.h>
-#include <R.h>
 #include <Rinternals.h>
 #include <Rmath.h>
 #include <R_ext/Rdynload.h>
@@ -24,11 +24,10 @@
 #include "config.h"
 #include "logging.h"
 
+// $ export BITS_PER_WORD=32 in the Terminal, before running R
 
 SEXP CCubes(SEXP tt) {
 
-    // $ export BITS_PER_WORD=32 in the Terminal, before running R
-
     // 32 bits per word, in bit shifting representation
     char *bits_per_word = getenv("BITS_PER_WORD"); // Read from the PATH
     int BITS_PER_WORD = bits_per_word ? atoi(bits_per_word) : 32;
@@ -36,9 +35,15 @@ SEXP CCubes(SEXP tt) {
         BITS_PER_WORD = 32; // Default to 32
     }
 
-    // $ export PRINT_INFO=1 in the Terminal, before running R
+    // $ export PRINT_INFO=true in the Terminal, before running R
     char *print_info = getenv("PRINT_INFO"); // Read from the PATH
-    Rboolean PRINT_INFO = print_info && print_info[0] == '1';
+    Rboolean PRINT_INFO = print_info && (
+        strcmp(print_info, "1") == 0 ||
+        strcmp(print_info, "TRUE") == 0 ||
+        strcmp(print_info, "true") == 0 ||
+        strcmp(print_info, "T") == 0 ||
+        strcmp(print_info, "t") == 0
+    );
 
     int multiplier = 0;
     struct timeval start, end;
@@ -112,12 +117,18 @@ SEXP CCubes(SEXP tt) {
         }
     }
 
-    int value_bit_width = 0;
+    int bits_needed = ceil(log2(max_value)); // Compute the necessary bits
-    while (max_value > 0) {
+    int value_bit_width = 1;
-        max_value >>= 1; // Shift right until no bits remain
+    while (value_bit_width < bits_needed) {
-        value_bit_width++;
+        value_bit_width *= 2; // Round up to the next power of 2
     }
 
+    if (value_bit_width > BITS_PER_WORD) {
+        BITS_PER_WORD = value_bit_width; // Adjust the bits per word
+    }
+
+    int value_bit_mask = (1U << value_bit_width) - 1U;
+
     // calculate the number of values (biggest number) for each input
     int nofvalues[ninputs];
     int nofpi[ninputs];
@@ -146,10 +157,10 @@ SEXP CCubes(SEXP tt) {
     int estimPI = 250000;
 
     // the index of the PIs, in descending order of their number of covered ON-set minterms
-    int *p_covered = R_Calloc(estimPI, int);
+    int *p_covered = (int *) R_Calloc(estimPI, int);
 
     // many PIs will have the same coverage, but they don't necessarily cover the same minterms
-    // to employ row dominance when solving the PI chart, we need to compare the coverage of the
+    // to employ row dominance when solving the coverage matrix, we need to compare the coverage of the
     // current PI with the coverage of the previous PIs. If this PI survives the comparison, its
     // coverage has to be added in the p_covered vector, and its order in the p_covered
     // vector, at the last index of the PI coverage with the same number of minterms
@@ -158,9 +169,9 @@ SEXP CCubes(SEXP tt) {
     // p_pichart = malloc(estimPI * posrows * sizeof(int));
     // memset(p_pichart, false, estimPI * posrows * sizeof(int));
     int *p_pichart = (int *) R_Calloc(estimPI * posrows, int);
-    // prefixing (int *) prefills in all values with 0s
+    // prefixing (int *) before R_Calloc() prefills all values with 0s
 
-    int pichart_words = (posrows + BITS_PER_WORD - 1) / BITS_PER_WORD; // Words needed per PI chart columns
+    int pichart_words = (posrows + BITS_PER_WORD - 1) / BITS_PER_WORD; // Words needed per coverage matrix columns
     unsigned int *p_pichart_pos = (unsigned int *) R_Calloc(estimPI * pichart_words, unsigned int);
     int implicant_words = (ninputs + BITS_PER_WORD - 1) / BITS_PER_WORD; // Words needed per PI representation
     unsigned int *p_implicants_pos = (unsigned int *) R_Calloc(estimPI * implicant_words, unsigned int);
@@ -168,15 +179,14 @@ SEXP CCubes(SEXP tt) {
 
     int prevfoundPI = 0;  // the number of previously found PIs
     int foundPI = 0;
-    int prevsolmin = 0;     // the minimum number of PIs that solve the PI chart
+    int prevsolmin = 0;   // the previous (level k - 1), minimum number of PIs to solve the coverage matrix
     int solmin = 0;
 
-    // the positions of the PIs solving the PI chart
+    // the positions of the PIs solving the coverage matrix
-    // a vector which can never be lengthier than the number of ON minterms (posrows)
+    // this vector can never be lengthier than the number of ON minterms (posrows)
     int previndices[posrows];
     int indices[posrows];
 
-
     for (int i = 0; i < posrows; i++) {
         previndices[i] = 0;
         indices[i] = 0;
@@ -200,6 +210,11 @@ SEXP CCubes(SEXP tt) {
         }
 
 	int n_tasks = nchoosek(ninputs, k);
+        if (n_tasks == 0) {
+            // overflow, too many tasks
+            return (R_NilValue);
+        }
+
 	int n_tasks_batch = 512;
         for (int task = 0; task < n_tasks; task+=n_tasks_batch) {
 		/* adjust if batch size is larger than total job size */
@@ -277,22 +292,22 @@ SEXP CCubes(SEXP tt) {
             #pragma omp parallel for schedule(dynamic)
         #endif
 
-        for (int task = 0; task < nchoosek(ninputs, k); task++) {
+        for (int task = 0; task < n_tasks; task++) {
             int tempk[k];
             int x = 0;
-            int start_point = task;
+            int combination = task;
 
-            // fill the combination for the current task
+            // fill the combination for the current task / combination number
             for (int i = 0; i < k; i++) {
-                while (nchoosek(ninputs - (x + 1), k - (i + 1)) <= start_point) {
+                while (nchoosek(ninputs - (x + 1), k - (i + 1)) <= combination) {
-                    start_point -= nchoosek(ninputs - (x + 1), k - (i + 1));
+                    combination -= nchoosek(ninputs - (x + 1), k - (i + 1));
                     x++;
                 }
                 tempk[i] = x;
                 x++;
             }
 
-            // allocate vectors of decimal row numbers for the positive and negative rows
+            // allocate vectors of decimal numbers for the ON-set and OFF-set rows
             int decpos[posrows];
             int decneg[negrows];
 
@@ -326,7 +341,7 @@ SEXP CCubes(SEXP tt) {
             int possible_rows[posrows];
 
             Rboolean possible_cover[posrows];
-            possible_cover[0] = true; // Rboolean flag, to be set with false if found among the OFF set
+            possible_cover[0] = true; // boolean flag, to be set with false if found among the OFF set
 
             int found = 0;
 
@@ -370,7 +385,6 @@ SEXP CCubes(SEXP tt) {
 
                 for (int f = 0; f < found; f++) {
 
-
                     // create a temporary vector of length k, containing the values from the initial ON set
                     // plus 1 (because 0 now signals a minimization, it becomes 1, and 1 becomes 2 etc.
                     int tempc[k];
@@ -388,11 +402,11 @@ SEXP CCubes(SEXP tt) {
                         int value = ON_set[tempk[c] * posrows + frows[f]];
                         tempc[c] = value + 1;
 
-                        int word_index = tempk[c] / BITS_PER_WORD;
+                        int word_index = tempk[c] / (BITS_PER_WORD / value_bit_width);
-                        int bit_index = tempk[c] % BITS_PER_WORD;
+                        int bit_index = (tempk[c] % (BITS_PER_WORD / value_bit_width)) * value_bit_width;
 
-                        fixed_bits[word_index] |= 1U << bit_index;
+                        fixed_bits[word_index] |= (value_bit_mask << bit_index);
-                        value_bits[word_index] |= (unsigned int)value << (bit_index * value_bit_width);
+                        value_bits[word_index] |= ((unsigned int)value << bit_index);
                     }
 
                     // check if the current PI is not redundant
@@ -414,13 +428,15 @@ SEXP CCubes(SEXP tt) {
 
                         for (int w = 0; w < implicant_words; w++) {
                             // If the new PI has values on positions outside the existing PI’s fixed positions, it’s not a subset
-                            if ((fixed_bits[w] & p_implicants_pos[i * implicant_words + w]) != p_implicants_pos[i * implicant_words + w]) {
+                            unsigned int index_mask = p_implicants_pos[i * implicant_words + w];
+
+                            if ((fixed_bits[w] & index_mask) != index_mask) {
                                 is_subset = false;
                                 break;
                             }
 
                             // then compare the value bits, if one or more values on those positions are different, it’s not a subset
-                            if ((value_bits[w] & p_implicants_val[i * implicant_words + w]) != p_implicants_val[i * implicant_words + w]) {
+                            if ((value_bits[w] & index_mask) != (p_implicants_val[i * implicant_words + w] & index_mask)) {
                                 is_subset = false;
                                 break;
                             }
@@ -479,7 +495,6 @@ SEXP CCubes(SEXP tt) {
                         #pragma omp critical
                     #endif
                     {
-
                         // push the PI information to the global arrays
 
                         for (int i = foundPI; i > last_index[covsum - 1]; i--) {
@@ -497,7 +512,7 @@ SEXP CCubes(SEXP tt) {
                             p_implicants_val[implicant_words * foundPI + w] = value_bits[w];
                         }
 
-                        // populate the PI chart
+                        // populate the coverage matrix
                         for (int r = 0; r < posrows; r++) {
                             for (int w = 0; w < pichart_words; w++) {
                                 p_pichart_pos[foundPI * pichart_words + w] = pichart_values[w];
@@ -510,21 +525,20 @@ SEXP CCubes(SEXP tt) {
 
                         // when needed, increase allocated memory
                         if (foundPI / estimPI > 0.9) {
-                            int old_size = estimPI;
                             estimPI += 100000;
                             p_pichart =        R_Realloc(p_pichart,        posrows * estimPI,         int);
-                            p_pichart_pos =    R_Realloc(p_pichart_pos,    estimPI,           unsigned int);
+                            p_pichart_pos =    R_Realloc(p_pichart_pos,    pichart_words * estimPI,   unsigned int);
-                            p_implicants_val = R_Realloc(p_implicants_val, ninputs * estimPI, unsigned int);
+                            p_implicants_val = R_Realloc(p_implicants_val, implicant_words * estimPI, unsigned int);
-                            p_implicants_pos = R_Realloc(p_implicants_pos, ninputs * estimPI, unsigned int);
+                            p_implicants_pos = R_Realloc(p_implicants_pos, implicant_words * estimPI, unsigned int);
                             p_covered =        R_Realloc(p_covered,        estimPI,                   int);
 
-                            for (unsigned int i = old_size; i < posrows * estimPI; i++) {
+                            for (unsigned int i = foundPI; i < posrows * estimPI; i++) {
                                 p_pichart[i] = 0;
                             }
-                            for (unsigned int i = old_size; i < estimPI; i++) {
+                            for (unsigned int i = foundPI; i < pichart_words * estimPI; i++) {
                                 p_pichart_pos[i] = 0U;
                             }
-                            for (unsigned int i = old_size; i < ninputs * estimPI; i++) {
+                            for (unsigned int i = foundPI; i < implicant_words * estimPI; i++) {
                                 p_implicants_val[i] = 0U;
                                 p_implicants_pos[i] = 0U;
                             }
@@ -565,7 +579,6 @@ SEXP CCubes(SEXP tt) {
         if (ON_set_covered) {
             // Rprintf("posrows: %d; foundPI: %d\n", posrows, foundPI);
 
-
             if (PRINT_INFO) {
                 gettimeofday(&end, NULL); // End time
                 elapsed_time = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1e6;
@@ -573,8 +586,6 @@ SEXP CCubes(SEXP tt) {
                 gettimeofday(&start, NULL);
             }
 
-
-
             SEXP pic = PROTECT(allocMatrix(INTSXP, posrows, foundPI));
             for (long long unsigned int i = 0; i < posrows * foundPI; i++) {
                 INTEGER(pic)[i] = p_pichart[i];
@@ -599,13 +610,12 @@ SEXP CCubes(SEXP tt) {
             }
 
             UNPROTECT(5);
+
             // Rprintf("solution minima: %d\n", solmin);
-
-
             if (PRINT_INFO) {
                 gettimeofday(&end, NULL);
                 elapsed_time = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1e6;
-                Rprintf("Time spent solving the PI chart: %f sec.\n", elapsed_time);
+                Rprintf("Time spent solving the coverage matrix: %f sec.\n", elapsed_time);
                 gettimeofday(&start, NULL);
             }
 
@@ -643,7 +653,7 @@ SEXP CCubes(SEXP tt) {
         }
     }
 
-    // printf("solmin: %d\n", solmin);
+    // Rprintf("solmin: %d\n", solmin);
     if (PRINT_INFO) {
         Rprintf("--- END ---\n");
     }
@@ -654,11 +664,13 @@ SEXP CCubes(SEXP tt) {
     for (int c = 0; c < solmin; c++) {
         for (int r = 0; r < ninputs; r++) {
             unsigned int value = 0;
-            int word_index = r / BITS_PER_WORD;  // Word index within the implicant
+            int word_index = r / (BITS_PER_WORD / value_bit_width);                     // Word index within the implicant
-            int bit_index = r % BITS_PER_WORD;  // Bit position within the word
+            int bit_index = (r % (BITS_PER_WORD / value_bit_width)) * value_bit_width;  // Bit position within the word
 
-            if (p_implicants_pos[indices[c] * implicant_words + word_index] & (1U << bit_index)) {
+                                                                           // (1U << bit_index) is just as good
-                value = 1U + ((p_implicants_val[indices[c] * implicant_words + word_index] >> (bit_index * value_bit_width)) & ((1U << value_bit_width) - 1U));
+            if (p_implicants_pos[indices[c] * implicant_words + word_index] & (value_bit_mask << bit_index)) {
+                value = (p_implicants_val[indices[c] * implicant_words + word_index] >> bit_index) & value_bit_mask;
+                value++; // 0 indicates a minimization, so we increment the value
             }
 
             p_sol[r * solmin + c] = value; // transposed
@@ -677,13 +689,12 @@ SEXP CCubes(SEXP tt) {
     return (sol);
 }
 
-
 long long unsigned int nchoosek(
     int n,
     int k
 ) {
+    if (k > n) return 0;
     if (k == 0 || k == n) return 1;
-    if (k == 1) return n;
 
     long long unsigned int result = 1;
 
@@ -692,7 +703,21 @@ long long unsigned int nchoosek(
     }
 
     for (int i = 0; i < k; i++) {
-        result = result * (n - i) / (i + 1);
+        // result = result * (n - i) / (i + 1);
+
+        // Check for potential overflow before multiplication
+        if (result > ULLONG_MAX / (n - i)) {
+            return 0; // Indicate overflow
+        }
+
+        result *= (n - i);
+
+        // Check for potential overflow before division
+        if (result % (i + 1) != 0) {
+            return 0; // Indicate overflow
+        }
+
+        result /= (i + 1);
     }
 
     return result;

src/ccubes.cl

@@ -147,12 +147,12 @@ ccubes_task(int k,
 
 	int tempk[NINPUTS]; /* max is tempk[ninputs] */
 	int x = 0;
-	int start_point = task;
+	int combination = task;
 
 	// fill the combination for the current task
 	for (int i = 0; i < k; i++) {
-		while (nchoosek(NINPUTS - (x + 1), k - (i + 1)) <= start_point) {
+		while (nchoosek(NINPUTS - (x + 1), k - (i + 1)) <= combination) {
-			start_point -= nchoosek(NINPUTS - (x + 1), k - (i + 1));
+			combination -= nchoosek(NINPUTS - (x + 1), k - (i + 1));
 			x++;
 		}
 		tempk[i] = x;
@@ -255,11 +255,11 @@ ccubes_task(int k,
 				int value = ON_set[tempk[c] * POSROWS + frows[f]];
 				tempc[c] = value + 1;
 
-				int word_index = tempk[c] / BITS_PER_WORD;
+				int word_index = tempk[c] / (BITS_PER_WORD / value_bit_width);
-				int bit_index = tempk[c] % BITS_PER_WORD;
+				int bit_index = (tempk[c] % (BITS_PER_WORD / value_bit_width)) * value_bit_width;
 
-				fixed_bits[word_index] |= 1U << bit_index;
+				fixed_bits[word_index] |= (value_bit_mask << bit_index);
-				value_bits[word_index] |= (unsigned int)value << (bit_index * VALUE_BIT_WIDTH);
+				value_bits[word_index] |= ((unsigned int)value << bit_index);
 			}
 
 			// check if the current PI is not redundant
@@ -281,13 +281,15 @@ ccubes_task(int k,
 
 				for (int w = 0; w < IMPLICANT_WORDS; w++) {
 					// If the new PI has values on positions outside the existing PI’s fixed positions, it’s not a subset
-					if ((fixed_bits[w] & p_implicants_pos[i * IMPLICANT_WORDS + w]) != p_implicants_pos[i * IMPLICANT_WORDS + w]) {
+					unsigned int index_mask = p_implicants_pos[i * implicant_words + w];
+
+					if ((fixed_bits[w] & index_mask) != index_mask) {
 						is_subset = false;
 						break;
 					}
 
 					// then compare the value bits, if one or more values on those positions are different, it’s not a subset
-					if ((value_bits[w] & p_implicants_val[i * IMPLICANT_WORDS + w]) != p_implicants_val[i * IMPLICANT_WORDS + w]) {
+					if ((value_bits[w] & index_mask) != (p_implicants_val[i * implicant_words + w] & index_mask)) {
 						is_subset = false;
 						break;
 					}

src/row_dominance.c

@@ -1,100 +0,0 @@
-
-#include <R_ext/Boolean.h>
-#include "row_dominance.h"
-
-
-void row_dominance(
-    int p_pichart[],
-    int p_implicants[],
-    int *p_ck,
-    int pirows,
-    unsigned int *foundPI,
-    int ninputs
-) {
-
-    unsigned int picols = *foundPI;
-
-
-    // int* survcols = (int *) R_Calloc (picols, int);
-
-    // for (int i = 0; i < picols; i++) {
-    //     survcols[i] = true;
-    // }
-
-    Rboolean survcols[picols];
-    int colsums[picols];
-    int sortcol[picols];
-    int temp;
-
-    for (unsigned int c = 0; c < picols; c++) {
-        colsums[c] = 0;
-
-        for (int r = 0; r < pirows; r++) {
-            colsums[c] += p_pichart[c * pirows + r];
-        }
-
-        sortcol[c] = c;
-        survcols[c] = true;
-    }
-
-    for (unsigned int c1 = 0; c1 < picols; c1++) {
-        for (unsigned int c2 = c1 + 1; c2 < picols; c2++) {
-            if (colsums[sortcol[c1]] < colsums[sortcol[c2]]) {
-                temp = sortcol[c1];
-                sortcol[c1] = sortcol[c2];
-                sortcol[c2] = temp;
-            }
-        }
-    }
-
-    for (unsigned int c1 = 0; c1 < picols; c1++) {
-        if (survcols[sortcol[c1]]) {
-            for (unsigned int c2 = c1 + 1; c2 < picols; c2++) {
-                if (survcols[sortcol[c2]]) {
-                    if (colsums[sortcol[c1]] > colsums[sortcol[c2]]) {
-
-                        Rboolean itcovers = true; // assume it's covered
-                        int r = 0;
-
-                        while (r < pirows && itcovers) {
-                            if (p_pichart[sortcol[c2] * pirows + r]) {
-                                itcovers = p_pichart[sortcol[c1] * pirows + r];
-                            }
-                            r++;
-                        }
-
-                        if (itcovers) {
-                            survcols[sortcol[c2]] = false;
-                            --(*foundPI);
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-
-    if (*foundPI < picols) {
-
-        // move (overwrite) all surviving columns towards the beginning
-        int s = 0;
-        for (unsigned int c = 0; c < picols; c++) {
-            if (survcols[c]) {
-                for (int r = 0; r < pirows; r++) {
-                    p_pichart[s * pirows + r] = p_pichart[c * pirows + r];
-                }
-
-                for (int r = 0; r < ninputs; r++) {
-                    p_implicants[s * ninputs + r] = p_implicants[c * ninputs + r];
-                }
-
-                // same with the vector positions
-                p_ck[s] = p_ck[c];
-
-                s++;
-            }
-        }
-    }
-
-    return;
-}

src/row_dominance.h

@@ -1,10 +0,0 @@
-#include <stdbool.h>
-
-void row_dominance(
-    int p_pichart[],
-    int p_implicants[],
-    int *p_ck,
-    int pirows,
-    unsigned int *foundPI,
-    int ninputs
-);