.Rbuildignore

@@ -1,6 +0,0 @@
-## Ignore travis config file
-^\.travis\.yml$
-\.vscode
-^.*\.Rproj$
-^\.Rproj\.user$
-^.*\.txt$

.gitignore

@@ -1,7 +0,0 @@
-*.js~
-*.marks
-*.history
-*.o
-*.so
-*.DS_Store
-instructiuni.txt

DESCRIPTION

@@ -1,21 +0,0 @@
-Package: IEEE
-Version: 1.0.0
-Title: For the IEEE article
-Authors@R: c(
-  person(given = "Adrian", family = "Dusa",
-  role = c("aut", "cre", "cph"),
-  email = "dusa.adrian@unibuc.ro",
-  comment = c(ORCID = "0000-0002-3525-9253"))
-  )
-URL: https://github.com/dusadrian/IEEE
-BugReports: https://github.com/dusadrian/IEEE/issues
-Depends: R (>= 3.6.0)
-Imports: LogicOpt, 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.
-  'QCA' is a methodology that bridges the qualitative and quantitative
-  divide in social science research. It uses a Boolean minimization
-  algorithm, resulting in a minimal causal configuration associated
-  with a given phenomenon.
-License: GPL (>= 3)

NAMESPACE

@@ -1,9 +0,0 @@
-importFrom("gurobi", "gurobi")
-importFrom("Matrix", "Matrix")
-importFrom("utils", "tail")
-importFrom("LogicOpt", "logicopt")
-export(gendat)
-export(make_infile)
-export(solvechart)
-useDynLib(IEEE, .registration = TRUE)
-

R/gendat.R

@@ -1,23 +0,0 @@
-
-gendat <- function(i, ncols, nrows) {
-    set.seed(i * nrows * ncols)
-
-    # all configurations BUT the outcome have to be unique
-    dat <- unique(
-        matrix(sample(0:1, (ncols - 1)*nrows, replace = TRUE), ncol = ncols - 1)
-    )
-
-    # now that configurations are unique, we can add the outcome
-    dat <- cbind(dat, sample(0:1, nrow(dat), replace = TRUE))
-
-    if (ncols > length(LETTERS) + 1) {
-        colnames(dat) <- c(
-            paste("I", seq(ncols - 1), sep = ""),
-            "OUT"
-        )
-    } else {
-        colnames(dat) <- c(LETTERS[seq(ncols - 1)], "OUT")
-    }
-
-    return(dat[order(dat[, ncols], decreasing = TRUE), ])
-}

R/make_infile.R

@@ -1,74 +0,0 @@
-
-make_infile <- function(dat, espname = "infile.esp", copy = FALSE, ...) {
-
-    on.exit(suppressWarnings(sink()))
-
-    dots <- list(...)
-    # number of outcomes
-    no <- ifelse(is.null(dots$no), 1, dots$no)
-
-    dat <- as.data.frame(dat)
-    inputs <- seq(ncol(dat) - no)
-    outputs <- setdiff(seq(ncol(dat)), inputs)
-    iname <- names(dat)[inputs]
-    oname <- names(dat)[outputs]
-    outcome <- dat[, outputs, drop = FALSE]
-    outcome[outcome == 0] <- "-"
-
-    # espname <- paste("infile", i, ncols, nrows, "esp", sep = ".")
-
-    sink(espname)
-    # on.exit(sink())
-    cat(paste(".i", length(inputs), "\n"))
-    cat(paste(".o", no, "\n"))
-    cat(paste(".ilb", paste(iname, collapse = " "), "\n"))
-    cat(paste(".ob", paste(oname, collapse = " "), "\n"))
-    cat(
-        paste(
-            paste(
-                apply(dat[, inputs], 1, paste, collapse = ""),
-                apply(outcome, 1, paste, collapse = ""),
-                sep = " "
-            ),
-            collapse = "\n"
-        )
-    )
-    cat("\n.e\n")
-    sink()
-
-
-    lo <- LogicOpt::logicopt(esp_file = espname, mode = "echo")[[1]]
-    outcome <- lo[, outputs, drop = FALSE]
-    outcome[] <- lapply(
-        outcome,
-        function(x) {
-            admisc::recode(x, "1=1; 0=-; else=0")
-        }
-    )
-
-    sink(espname)
-    cat(paste(".i", length(inputs), "\n"))
-    cat(paste(".o", no, "\n"))
-    cat(paste(".ilb", paste(iname, collapse = " "), "\n"))
-    cat(paste(".ob", paste(oname, collapse = " "), "\n"))
-    cat(
-        paste(
-            paste(
-                apply(lo[, inputs], 1, paste, collapse = ""),
-                apply(outcome, 1, paste, collapse = ""),
-                sep = " "
-            ),
-            collapse = "\n"
-        )
-    )
-    cat("\n.e\n")
-    sink()
-
-    if (copy) {
-        i <- dots$i
-        if (is.null(i)) {
-            i <- get("i", envir = parent.frame())
-        }
-        file.copy(espname, paste("infile", i, ncol(dat), nrow(dat), "esp", sep = "."))
-    }
-}

R/solvechart.R

@@ -1,133 +0,0 @@
-`solvechart` <- function(x, ...) {
-
-    dots <- list(...)
-
-    if (all(rowSums(x) > 0)) {
-
-        PIlayers <- attr(x, "PIlayers")
-
-        # number of:
-        PIs <- ncol(x)
-        poc <- nrow(x) # positive observed configurations (minterms)
-
-        # sink("timings.txt", append = TRUE)
-        # cat(sprintf("Minterms: %s, PIs: %s", poc, PIs))
-        # sink()
-
-        model <- list(
-            A = x,
-            modelsense = "min",
-            rhs = rep(1, poc),
-            sense = rep(">=", poc),
-            vtype = rep("B", PIs)
-        )
-
-        if (!is.null(PIlayers)) {
-            PIlayers <- PIlayers[PIlayers > 0]
-        }
-
-        if (length(PIlayers) < 2 || isFALSE(dots$multiobj)) {
-            model$obj <- rep(1, PIs)
-        } else {
-            # cat(paste(PIlayers, collapse = " "))
-
-            for (i in seq(length(PIlayers), 2)) {
-                PIlayers[i] <- PIlayers[i] - PIlayers[i - 1]
-            }
-
-            # cat(" -- ",  paste(PIlayers, collapse = " "), "\n")
-
-            value <- rep(
-                rev(2^seq(0, length(PIlayers) - 1)),
-                PIlayers
-            )
-
-            if (length(value) < PIs) {
-                value <- c(value, rep(tail(value, n = 1), PIs - length(value)))
-            } else if (PIs < length(value)) {
-                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 = ""),
-            #     col.names = FALSE
-            # )
-
-            # admisc::export(
-            #     as.data.frame(t(attr(x, "implicants"))),
-            #     file = paste("imp", sample(5000, 1), ".csv", sep = ""),
-            #     col.names = FALSE
-            # )
-
-            model$multiobj <- list(
-                list(
-                    objn = rep(1, PIs),
-                    priority = 1,
-                    weight = 1
-                ),
-                list(
-                    objn = -1 * value,
-                    priority = 0,
-                    weight = 1
-                )
-            )
-        }
-
-        tc <- admisc::tryCatchWEM(
-            solution <- gurobi::gurobi(
-                model,
-                params = list(
-                    OutputFlag = 0,
-                    LogToConsole = 0
-                )
-            )
-        )
-
-        if (is.null(tc$error)) {
-            # sink("timings.txt", append = TRUE)
-            # cat(sprintf(", Time: %s\n", round(solution$runtime, 3)))
-            # sink()
-
-            if (round(solution$objval[1]) < solution$objval[1]) {
-                # the weighted method did not yield precise results
-                temp <- sapply(solution$pool, function(x) {
-                    return(round(x$objval[1]) == x$objval[1])
-                })
-
-                if (any(temp)) {
-                    solution <- solution$pool[[which(temp)[1]]]$xn
-                } else {
-                    # There is still no ideintified solution with an integer number of PIs
-                    # Give up weigthing
-                    model$multiobj <- NULL
-
-                    # and return the uneighted solution
-                    model$obj <- rep(1, PIs)
-                    solution <- gurobi::gurobi(
-                        model,
-                        params = list(
-                            OutputFlag = 0,
-                            LogToConsole = 0
-                        )
-                    )$x
-                }
-            } else {
-                solution <- solution$x
-            }
-        } else {
-            solution <- lpSolve::lp(
-                "min",
-                rep(1, ncol(x)),
-                x,
-                ">=",
-                1,
-                int.vec = seq(nrow(x)),
-                all.bin = TRUE
-            )$solution
-        }
-
-        return(as.integer(which(solution > 0)))
-    }
-}

README.md

@@ -1,16 +0,0 @@
-BUILD
------
-cd ${HOME}
-R CMD build IEEE
-export LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu:/opt/gurobi1200/linux64/lib
-R CMD INSTALL IEEE_1.0.0.tar.gz
-
-TEST
-----
-R
---
-> library(IEEE)
-> compare(1,20,30,file="R-log.txt")
-
-$ cat R-log.txt
-1, 20, 30, 0.006, 0.357, "(!N&!P&!R) + (J&O&!S) + (B&!H&!Q) + (J&Q)", "(J&Q) + (!A&B&!O) + (F&O&!Q)", 3332, 233

ccubes.cl

@@ -109,9 +109,9 @@ nchoosek(int n, int k)
 #endif
 __kernel void
 ccubes_task(int k,
-		__global const int *nofvalues,			/* IN: RC */
-		__global const int *ON_set,			/* IN: RC */
-		__global const int *OFF_set,			/* IN: RC */
+		__global const real *nofvalues,			/* IN: RC */
+		__global const real *ON_set,			/* IN: RC */
+		__global const real *OFF_set,			/* IN: RC */
 		__global const unsigned int *p_implicants_pos,	/* IN: RC */
 		__global const unsigned int *p_implicants_val,	/* IN: RC */
 		__global const int *last_index,			/* IN: RC */

cl_setup.c

@@ -284,7 +284,6 @@ cl_build(struct cl_uctx uctx, cl_device_type dev,
 	FILE *kern_file = NULL;
 	char *kern_src = NULL;
 	size_t srcsz = 0;
-	size_t ret;
 
 	int type;
 
@@ -327,11 +326,7 @@ cl_build(struct cl_uctx uctx, cl_device_type dev,
 		result = CL_INVALID_VALUE;
 		goto err;
 	}
-	ret = fread(kern_src, 1, srcsz, kern_file);
-	if (ret != srcsz) {
-		log_warn("cl", "fread() failed!");
-		goto err;
-	}
+	fread(kern_src, 1, srcsz, kern_file);
 	kern_src[srcsz] = 0;
 
 	log_info("cl", "FILE DUMP BEGINS");

clccubes.c

@@ -133,9 +133,9 @@ err:
 
 int
 ccubes_alloc(struct ccubes_context *ctx,
-    int *nofvalues,			/* IN: RC */
-    int *ON_set,				/* IN: RC */
-    int *OFF_set,			/* IN: RC */
+    real *nofvalues,			/* IN: RC */
+    real *ON_set,				/* IN: RC */
+    real *OFF_set,			/* IN: RC */
     unsigned int *p_implicants_pos,	/* IN: RC */
     unsigned int *p_implicants_val,	/* IN: RC */
     int *last_index,			/* IN: RC */
@@ -161,24 +161,24 @@ ccubes_alloc(struct ccubes_context *ctx,
 	 * INPUTS
 	 */
 
-	/* __global const int *nofvalues, IN: RC */
+	/* __global const real *nofvalues, IN: RC */
 	ctx->nofvalues = clCreateBuffer(ctx->clctx->ctx,
 	    CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	    ctx->ninputs * sizeof(int), nofvalues, &rc);
+	    ctx->ninputs * sizeof(real), nofvalues, &rc);
 	if (rc != CL_SUCCESS) {
 		goto err;
 	}
-	/* __global const int *ON_set,	IN: RC */
+	/* __global const real *ON_set,	IN: RC */
 	ctx->ON_set = clCreateBuffer(ctx->clctx->ctx,
 	    CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	    ctx->posrows * ctx->ninputs * sizeof(int), ON_set, &rc);
+	    ctx->posrows * ctx->ninputs * sizeof(real), ON_set, &rc);
 	if (rc != CL_SUCCESS) {
 		goto err;
 	}
-	/*	__global const int *OFF_set, IN: RC */
+	/*	__global const real *OFF_set, IN: RC */
 	ctx->OFF_set = clCreateBuffer(ctx->clctx->ctx,
 	    CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	    ctx->ninputs * ctx->negrows * sizeof(int), OFF_set, &rc);
+	    ctx->ninputs * ctx->negrows * sizeof(real), OFF_set, &rc);
 	if (rc != CL_SUCCESS) {
 		goto err;
 	}
@@ -381,7 +381,7 @@ err:
 }
 
 int
-ccubes_do_tasks(int n_tasks,
+ccubes(int n_tasks,
     int n_tasks_off,
     int k,
     int ninputs,
@@ -391,9 +391,9 @@ ccubes_do_tasks(int n_tasks,
     int value_bit_width,
     int pichart_words,
     int estimPI,
-    int *nofvalues,			/* IN: RC */
-    int *ON_set,			/* IN: RC */
-    int *OFF_set,			/* IN: RC */
+    real *nofvalues,			/* IN: RC */
+    real *ON_set,			/* IN: RC */
+    real *OFF_set,			/* IN: RC */
     unsigned int *p_implicants_pos,	/* IN: RC */
     unsigned int *p_implicants_val,	/* IN: RC */
     int *last_index,			/* IN: RC */

clccubes.h

@@ -70,11 +70,10 @@ struct ccubes_context {
 
 	/* ND-Range */
 	size_t gws; 	/* global work size */
-	size_t goff; 	/* global offset */
 };
 
 int
-ccubes_do_tasks(int n_tasks,
+ccubes(int n_tasks,
     int n_tasks_off,
     int k,
     int ninputs,
@@ -84,9 +83,9 @@ ccubes_do_tasks(int n_tasks,
     int value_bit_width,
     int pichart_words,
     int estimPI,
-    int *nofvalues,			/* IN: RC */
-    int *ON_set,			/* IN: RC */
-    int *OFF_set,			/* IN: RC */
+    real *nofvalues,			/* IN: RC */
+    real *ON_set,			/* IN: RC */
+    real *OFF_set,			/* IN: RC */
     unsigned int *p_implicants_pos,	/* IN: RC */
     unsigned int *p_implicants_val,	/* IN: RC */
     int *last_index,			/* IN: RC */
@@ -98,4 +97,9 @@ ccubes_do_tasks(int n_tasks,
     unsigned int *pichart_values	/* OUT: RW */
 );
 
+int
+clccubes(struct ccubes_context *ccubesctx, cl_mem alpha0, cl_mem G, size_t signals,
+    size_t atoms, uint32_t sparsity, uint32_t pcoding, cl_mem gamma, cl_uint
+    num_events_in_wait_list, const cl_event *event_wait_list, cl_event *ev_ccubes);
+
 #endif

src/CCubes.c

@@ -1,643 +0,0 @@
-#include <R_ext/RS.h>
-#include <R_ext/Boolean.h>
-#include <R.h>
-#include <Rinternals.h>
-#include <Rmath.h>
-#include <R_ext/Rdynload.h>
-#include <sys/time.h>
-
-
-#include <stdbool.h>
-#include <math.h>
-#include "CCubes.h"
-
-#ifdef _OPENMP
-    #undef match
-    #include <omp.h>
-#endif
-
-#include "real.h"
-#include "cl_setup.h"
-
-#include "clccubes.h"
-
-#include "config.h"
-#include "logging.h"
-
-
-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;
-    if (BITS_PER_WORD != 8 && BITS_PER_WORD != 16 && BITS_PER_WORD != 32 && BITS_PER_WORD != 64) {
-        BITS_PER_WORD = 32;  // Default to 32
-    }
-
-    // $ export PRINT_INFO=1 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';
-
-    int multiplier = 0;
-    struct timeval start, end;
-    double elapsed_time;
-
-    config_set_int("log", LOG_LEVEL_WARN);
-    config_set_int("log:ccubes", LOG_LEVEL_WARN);
-    config_set_int("cl", LOG_LEVEL_DEBUG);
-
-    if (PRINT_INFO) {
-        Rprintf("--- START ---\n");
-        gettimeofday(&start, NULL); // Start time
-    }
-
-    int *p_tt = INTEGER(tt);
-    int ttrows = nrows(tt); // number of rows in the data matrix
-    int ninputs = ncols(tt) - 1; // number of inputs (columns - 1, the last one is the outcome)
-
-    // calculate the number of positive output rows (the ON set)
-    int posrows = 0;
-    for (int r = 0; r < ttrows; r++) {
-        posrows += p_tt[ninputs * ttrows + r];
-    }
-
-    // calculate the number of negative output rows (the OFF set)
-    int negrows = ttrows - posrows;
-
-    if (negrows == 0) {
-        // if there are no negative output rows, no PIs can be found
-        // all inputs will be completely minimized
-        return(R_NilValue);
-    }
-
-    // split the minterms in the ON and OFF set matrices
-    int ON_set[posrows * ninputs];
-    int OFF_set[ninputs * negrows];
-    int rowpos = 0, rowneg = 0;
-    int max_value = 0;
-
-    for (int r = 0; r < ttrows; r++) {
-        if (p_tt[ninputs * ttrows + r] == 1) { // positive
-            for (int c = 0; c < ninputs; c++) {
-                int value = p_tt[c * ttrows + r];
-                ON_set[c * posrows + rowpos] = value;
-                if (value > max_value) {
-                    max_value = value;
-                }
-            }
-            rowpos += 1;
-        }
-        else { // negative
-            for (int c = 0; c < ninputs; c++) {
-                int value = p_tt[c * ttrows + r];
-                OFF_set[c * negrows + rowneg] = value;
-                if (value > max_value) {
-                    max_value = value;
-                }
-            }
-            rowneg += 1;
-        }
-    }
-
-    int value_bit_width = 0;
-    while (max_value > 0) {
-        max_value >>= 1; // Shift right until no bits remain
-        value_bit_width++;
-    }
-
-    // calculate the number of values (biggest number) for each input
-    int nofvalues[ninputs];
-    int nofpi[ninputs];
-
-    for (int i = 0; i < ninputs; i++) {
-        nofvalues[i] = 0; // initialize
-        nofpi[i] = 0; // initialize
-
-        for (int r = 0; r < ttrows; r++) {
-            if (nofvalues[i] < p_tt[i * ttrows + r]) {
-                nofvalues[i] = p_tt[i * ttrows + r];
-            }
-        }
-
-        // add 1 because if the biggest number is 2 then it has three levels: 0, 1 and 2
-        nofvalues[i] += 1;
-
-        if (nofvalues[i] == 1) {
-            // no input ever has less than two values
-            nofvalues[i] = 2;
-        }
-    }
-
-    // preallocating for an estimated large number of 10000 found PIs
-    // this number will be iteratively increased when the found PIs reach the upper limit
-    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);
-
-    // 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
-    // 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
-    int last_index[posrows]; // descending order
-
-    // 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
-
-    int pichart_words = (posrows + BITS_PER_WORD - 1) / BITS_PER_WORD; // Words needed per PI chart 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);
-    unsigned int *p_implicants_val = (unsigned int *) R_Calloc(estimPI * implicant_words, unsigned int);
-
-    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 solmin = 0;
-
-    // the positions of the PIs solving the PI chart
-    // a vector which 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;
-        last_index[i] = 0;
-    }
-
-    Rboolean ON_set_covered = false;
-    if (PRINT_INFO) {
-        Rprintf("ON-set minterms: %d\n", posrows);
-        #ifdef _OPENMP
-            Rprintf("OpenMP enabled, %d workers\n", omp_get_max_threads());
-        #endif
-    }
-
-
-    int stop_counter = 0; // to stop if two consecutive levels of complexity yield no PIs
-    int k;
-    for (k = 1; k <= ninputs; k++) {
-        if (PRINT_INFO) {
-            Rprintf("---k: %d\n", k);
-        }
-
-	int n_tasks = 512;
-        for (int task = 0; task < nchoosek(ninputs, k); task+=n_tasks) {
-		bool *coverage;
-		unsigned int *fixed_bits;
-		unsigned int *value_bits;
-		unsigned int *pichart_values;
-		ccubes_do_tasks(n_tasks,
-		    task,
-		    k,
-		    ninputs,
-		    posrows,
-		    negrows,
-		    implicant_words,
-		    value_bit_width,
-		    pichart_words,
-		    estimPI,
-		    nofvalues,
-		    ON_set,
-		    OFF_set,
-		    p_implicants_pos,
-		    p_implicants_val,
-		    last_index,
-		    p_covered,
-		    p_pichart_pos,
-		    coverage,
-		    fixed_bits,
-		    value_bits,
-		    pichart_values
-			);
-	}
-
-        #ifdef _OPENMP
-            #pragma omp parallel for schedule(dynamic)
-        #endif
-
-        for (int task = 0; task < nchoosek(ninputs, k); task++) {
-            int tempk[k];
-            int x = 0;
-            int start_point = 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) {
-                    start_point -= nchoosek(ninputs - (x + 1), k - (i + 1));
-                    x++;
-                }
-                tempk[i] = x;
-                x++;
-            }
-
-            // allocate vectors of decimal row numbers for the positive and negative rows
-            int decpos[posrows];
-            int decneg[negrows];
-
-            // create the vector of multiple bases, useful when calculating the decimal representation
-            // of a particular combination of columns, for each row
-            int mbase[k];
-            mbase[0] = 1; // the first number is _always_ equal to 1, irrespective of the number of values in a certain input
-
-            // calculate the vector of multiple bases, for example if we have k = 3 (three inputs) with
-            // 2, 3 and 2 values then mbase will be [1, 2, 6] from: 1, 1 * 2 = 2, 2 * 3 = 6
-            for (int i = 1; i < k; i++) {
-                mbase[i] = mbase[i - 1] * nofvalues[tempk[i - 1]];
-            }
-
-            // calculate decimal numbers, using mbase, fills in decpos and decneg
-            for (int r = 0; r < posrows; r++) {
-                decpos[r] = 0;
-                for (int c = 0; c < k; c++) {
-                    decpos[r] += ON_set[tempk[c] * posrows + r] * mbase[c];
-                }
-            }
-
-            for (int r = 0; r < negrows; r++) {
-                decneg[r] = 0;
-                for (int c = 0; c < k; c++) {
-                    decneg[r] += OFF_set[tempk[c] * negrows + r] * mbase[c];
-                }
-            }
-
-
-            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
-
-            int found = 0;
-
-            // identifies all unique decimal rows, for the selected combination of k inputs
-            for (int r = 0; r < posrows; r++) {
-                int prev = 0;
-                Rboolean unique = true; // Rboolean flag, assume the row is unique
-                while (prev < found && unique) {
-                    unique = decpos[possible_rows[prev]] != decpos[r];
-                    prev++;
-                }
-
-                if (unique) {
-                    possible_rows[found] = r;
-                    possible_cover[found] = true;
-                    found++;
-                }
-            }
-
-            if (found > 0) {
-                // some of the ON set numbers are possible PIs (not found in the OFF set)
-                int frows[found];
-
-                // verify if this is a possible PI
-                // (if the same decimal number is not found in the OFF set)
-                for (int i = found - 1; i >= 0; i--) {
-                    int j = 0;
-                    while (j < negrows && possible_cover[i]) {
-                        if (decpos[possible_rows[i]] == decneg[j]) {
-                            possible_cover[i] = false;
-                            found--;
-                        }
-                        j++;
-                    }
-
-                    if (possible_cover[i]) {
-                        frows[found - i - 1] = possible_rows[i];
-                    }
-                }
-                // Rprintf("task: %d; rows: %d\n", task, found);
-
-                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];
-
-                    // using bit shifting, store the fixed bits and value bits
-                    unsigned int fixed_bits[implicant_words];
-                    unsigned int value_bits[implicant_words];
-
-                    for (int i = 0; i < implicant_words; i++) {
-                        fixed_bits[i] = 0U;
-                        value_bits[i] = 0U;
-                    }
-
-                    for (int c = 0; c < k; c++) {
-                        int value = ON_set[tempk[c] * posrows + frows[f]];
-                        tempc[c] = value + 1;
-
-                        int word_index = tempk[c] / BITS_PER_WORD;
-                        int bit_index = tempk[c] % BITS_PER_WORD;
-
-                        fixed_bits[word_index] |= 1U << bit_index;
-                        value_bits[word_index] |= (unsigned int)value << (bit_index * value_bit_width);
-                    }
-
-                    // check if the current PI is not redundant
-                    Rboolean redundant = false;
-
-                    int i = 0;
-                    while (i < prevfoundPI && !redundant) {
-                        // /*
-                        // - ck contains the complexity level for each of the previously found non-redundant PIs
-                        // - indx is a matrix containing the indexes of the columns where the values were stored
-                        // - a redundant PI is one for which all values from a previous PI are exactly the same:
-                        // 0 0 1 2 0, let's say previously found PI
-                        // which means a corresponding ck = 2 and a corresponding indx = [3, 4]
-                        // 0 0 1 2 1 is redundant because on both columns 3 and 4 the values are equal
-                        // therefore sumeq = 2 and it will be equal to v = 2 when reaching the complexity level ck = 2
-                        // */
-
-                        Rboolean is_subset = true; // Assume it's a subset unless proven otherwise
-
-                        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]) {
-                                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]) {
-                                is_subset = false;
-                                break;
-                            }
-                        }
-
-                        redundant = is_subset;
-
-                        i++;
-                    }
-
-                    if (redundant) continue;
-
-                    Rboolean coverage[posrows];
-                    int covsum = 0;
-                    unsigned int pichart_values[pichart_words];
-                    for (int w = 0; w < pichart_words; w++) {
-                        pichart_values[w] = 0U;
-                    }
-
-                    for (int r = 0; r < posrows; r++) {
-                        coverage[r] = decpos[r] == decpos[frows[f]];
-                        if (coverage[r]) {
-                            int word_index = r / BITS_PER_WORD;
-                            int bit_index = r % BITS_PER_WORD;
-                            pichart_values[word_index] |= (1U << bit_index);
-                        }
-                        covsum += coverage[r];
-                    }
-
-                    // verify row dominance
-                    int rd = 0;
-                    while (rd < last_index[covsum - 1] && !redundant) {
-
-                        bool dominated = true;
-                        for (int w = 0; w < pichart_words; w++) {
-                            if ((pichart_values[w] & p_pichart_pos[p_covered[rd] * pichart_words + w]) != pichart_values[w]) {
-                                dominated = false;
-                                break;
-                            }
-                        }
-
-                        redundant = dominated;
-                        rd++;
-                    }
-
-                    if (redundant) continue;
-
-
-                    // Rprintf("It is a prime implicant\n");
-                    // This operation first gets a new index to push in the global array in a concurrent way
-                    // Then adds the result there.
-                    // We could synchronize only the index and let the copy operation happen in parallel BUT this
-                    // creates a false sharing problem and the performance is down by several factors.
-
-                    #ifdef _OPENMP
-                        #pragma omp critical
-                    #endif
-                    {
-
-                        // push the PI information to the global arrays
-
-                        for (int i = foundPI; i > last_index[covsum - 1]; i--) {
-                            p_covered[i] = p_covered[i - 1];
-                        }
-
-                        p_covered[last_index[covsum - 1]] = foundPI;
-
-                        for (int l = 1; l < covsum; l++) {
-                            last_index[l - 1] += 1;
-                        }
-
-                        for (int w = 0; w < implicant_words; w++) {
-                            p_implicants_pos[implicant_words * foundPI + w] = fixed_bits[w];
-                            p_implicants_val[implicant_words * foundPI + w] = value_bits[w];
-                        }
-
-                        // populate the PI chart
-                        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];
-                            }
-
-                            p_pichart[posrows * foundPI + r] = coverage[r];
-                        }
-
-                        ++foundPI;
-
-                        // 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_implicants_val = R_Realloc(p_implicants_val, ninputs * estimPI, unsigned int);
-                            p_implicants_pos = R_Realloc(p_implicants_pos, ninputs * estimPI, unsigned int);
-                            p_covered =        R_Realloc(p_covered,        estimPI,           int);
-
-                            for (unsigned int i = old_size; i < posrows * estimPI; i++) {
-                                p_pichart[i] = 0;
-                            }
-                            for (unsigned int i = old_size; i < estimPI; i++) {
-                                p_pichart_pos[i] = 0U;
-                            }
-                            for (unsigned int i = old_size; i < ninputs * estimPI; i++) {
-                                p_implicants_val[i] = 0U;
-                                p_implicants_pos[i] = 0U;
-                            }
-
-                            if (PRINT_INFO) {
-                                multiplier++;
-                                Rprintf("%dx ", multiplier);
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        nofpi[k - 1] = foundPI;
-
-        if (foundPI > 0 && !ON_set_covered) {
-            Rboolean test_coverage = true;
-
-            int r = 0;
-            while (r < posrows && test_coverage) {
-
-                Rboolean minterm_covered = false;
-                int c = 0;
-
-                while (c < foundPI && !minterm_covered) {
-                    minterm_covered = p_pichart[c * posrows + r];
-                    c++;
-                }
-
-                test_coverage = minterm_covered;
-                r++;
-            }
-
-            ON_set_covered = test_coverage;
-        }
-
-        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;
-                Rprintf("Time taken finding %d PIs: %f sec.\n", foundPI, elapsed_time);
-                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];
-            }
-
-            SEXP PIlayers = PROTECT(allocVector(INTSXP, ninputs));
-            for (int i = 0; i < ninputs; i++) {
-                INTEGER(PIlayers)[i] = nofpi[i];
-            }
-            setAttrib(pic, install("PIlayers"), PIlayers);
-
-            // if this file is run directly using SHLIB, the following line is needed
-            // R_ParseEvalString("library(IEEE)", R_GlobalEnv);
-
-            SEXP pkg_env = PROTECT(R_FindNamespace(mkString("IEEE")));
-            SEXP solvechart = PROTECT(Rf_findVarInFrame(pkg_env, Rf_install("solvechart")));
-            SEXP evalinR = PROTECT(R_tryEval(Rf_lang2(solvechart, pic), pkg_env, NULL));
-
-            solmin = length(evalinR);
-            for (int i = 0; i < solmin; i++) {
-                indices[i] = INTEGER(evalinR)[i] - 1; // R is 1-based
-            }
-
-            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);
-                gettimeofday(&start, NULL);
-            }
-
-            if (solmin == prevsolmin) {
-                // the minimum number of PIs did not change in the current level of complexity
-                // we can safely retain the less complex PIs from the previous level
-                for (int i = 0; i < solmin; i++) {
-                    indices[i] = previndices[i];
-                }
-                stop_counter += 1;
-            }
-            else {
-                // this means solmin is in fact smaller than the previously found solmin
-                // or it is the very first time a solmin was found
-                // only here it makes sense to overwrite prevsolmin and previndices,
-                // otherwise they are just as good as the ones from the previous level
-
-                prevsolmin = solmin;
-                for (int i = 0; i < solmin; i++) {
-                    previndices[i] = indices[i];
-                }
-
-                stop_counter = 0;
-            }
-        }
-
-        prevfoundPI = foundPI;
-
-        // printf("stop_counter: %d\n", stop_counter);
-
-        // One level of complexity up, and the solution minima does not change
-        if (stop_counter > 0) {
-            // the search can stop
-            break;
-        }
-    }
-
-    // printf("solmin: %d\n", solmin);
-    if (PRINT_INFO) {
-        Rprintf("--- END ---\n");
-    }
-
-    SEXP sol = PROTECT(allocMatrix(INTSXP, solmin, ninputs));
-    int *p_sol = INTEGER(sol);
-
-    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 bit_index = r % BITS_PER_WORD;  // Bit position within the word
-
-            if (p_implicants_pos[indices[c] * implicant_words + word_index] & (1U << bit_index)) {
-                value = 1U + ((p_implicants_val[indices[c] * implicant_words + word_index] >> (bit_index * value_bit_width)) & ((1U << value_bit_width) - 1U));
-            }
-
-            p_sol[r * solmin + c] = value; // transposed
-        }
-    }
-
-    R_Free(p_pichart);
-    R_Free(p_implicants_val);
-    R_Free(p_implicants_pos);
-    R_Free(p_pichart_pos);
-    R_Free(p_covered);
-
-    UNPROTECT(1);
-    return (sol);
-}
-
-
-long long unsigned int nchoosek(
-    int n,
-    int k
-) {
-    if (k == 0 || k == n) return 1;
-    if (k == 1) return n;
-
-    long long unsigned int result = 1;
-
-    if (k > n - k) {
-        k = n - k;
-    }
-
-    for (int i = 0; i < k; i++) {
-        result = result * (n - i) / (i + 1);
-    }
-
-    return result;
-}

src/CCubes.h

@@ -1,16 +0,0 @@
-
-SEXP CCubes(
-    SEXP tt
-);
-
-void resize(
-    int **array,
-    const int rows,
-    const unsigned int newcols,
-    const unsigned int oldcols
-);
-
-long long unsigned int nchoosek(
-    int n,
-    int k
-);

src/Makevars

@@ -1,2 +0,0 @@
-PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS)
-PKG_LIBS = $(SHLIB_OPENMP_CFLAGS)

src/registerDynamicSymbol.c

@@ -1,8 +0,0 @@
-#include <R.h>
-#include <Rinternals.h>
-#include <R_ext/Rdynload.h>
-
-void R_init_IEEE(DllInfo* info) {
-  R_registerRoutines(info, NULL, NULL, NULL, NULL);
-  R_useDynamicSymbols(info, TRUE);
-}

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
-);