From 143c29f407d580249f809e2e8dc113ba01c75775 Mon Sep 17 00:00:00 2001 From: Paul Irofti Date: Wed, 26 Mar 2025 15:15:57 +0200 Subject: [PATCH] =?UTF-8?q?Integrate=20R=20package=20from=20Adrian=20Du?= =?UTF-8?q?=C8=99a?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .Rbuildignore | 6 + .gitignore | 7 + DESCRIPTION | 21 ++ NAMESPACE | 9 + R/gendat.R | 23 ++ R/make_infile.R | 74 ++++ R/solvechart.R | 133 ++++++++ src/CCubes.c | 643 +++++++++++++++++++++++++++++++++++ src/CCubes.h | 16 + src/Makevars | 2 + ccubes.cl => src/ccubes.cl | 6 +- cl_setup.c => src/cl_setup.c | 7 +- cl_setup.h => src/cl_setup.h | 0 clccubes.c => src/clccubes.c | 26 +- clccubes.h => src/clccubes.h | 14 +- config.c => src/config.c | 0 config.h => src/config.h | 0 logging.c => src/logging.c | 0 logging.h => src/logging.h | 0 queue.h => src/queue.h | 0 real.h => src/real.h | 0 src/registerDynamicSymbol.c | 8 + src/row_dominance.c | 100 ++++++ src/row_dominance.h | 10 + tree.h => src/tree.h | 0 25 files changed, 1079 insertions(+), 26 deletions(-) create mode 100644 .Rbuildignore create mode 100644 .gitignore create mode 100644 DESCRIPTION create mode 100755 NAMESPACE create mode 100644 R/gendat.R create mode 100644 R/make_infile.R create mode 100644 R/solvechart.R create mode 100644 src/CCubes.c create mode 100755 src/CCubes.h create mode 100644 src/Makevars rename ccubes.cl => src/ccubes.cl (98%) rename cl_setup.c => src/cl_setup.c (98%) rename cl_setup.h => src/cl_setup.h (100%) rename clccubes.c => src/clccubes.c (95%) rename clccubes.h => src/clccubes.h (82%) rename config.c => src/config.c (100%) rename config.h => src/config.h (100%) rename logging.c => src/logging.c (100%) rename logging.h => src/logging.h (100%) rename queue.h => src/queue.h (100%) rename real.h => src/real.h (100%) create mode 100644 src/registerDynamicSymbol.c create mode 100644 src/row_dominance.c create mode 100755 src/row_dominance.h rename tree.h => src/tree.h (100%) diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..91db4ad --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,6 @@ +## Ignore travis config file +^\.travis\.yml$ +\.vscode +^.*\.Rproj$ +^\.Rproj\.user$ +^.*\.txt$ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f16e9bf --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +*.js~ +*.marks +*.history +*.o +*.so +*.DS_Store +instructiuni.txt \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..915afe7 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,21 @@ +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) diff --git a/NAMESPACE b/NAMESPACE new file mode 100755 index 0000000..3faae86 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,9 @@ +importFrom("gurobi", "gurobi") +importFrom("Matrix", "Matrix") +importFrom("utils", "tail") +importFrom("LogicOpt", "logicopt") +export(gendat) +export(make_infile) +export(solvechart) +useDynLib(IEEE, .registration = TRUE) + diff --git a/R/gendat.R b/R/gendat.R new file mode 100644 index 0000000..033d8b0 --- /dev/null +++ b/R/gendat.R @@ -0,0 +1,23 @@ + +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), ]) +} diff --git a/R/make_infile.R b/R/make_infile.R new file mode 100644 index 0000000..ef1ac69 --- /dev/null +++ b/R/make_infile.R @@ -0,0 +1,74 @@ + +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 = ".")) + } +} diff --git a/R/solvechart.R b/R/solvechart.R new file mode 100644 index 0000000..725df39 --- /dev/null +++ b/R/solvechart.R @@ -0,0 +1,133 @@ +`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))) + } +} diff --git a/src/CCubes.c b/src/CCubes.c new file mode 100644 index 0000000..05e8e64 --- /dev/null +++ b/src/CCubes.c @@ -0,0 +1,643 @@ +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include "CCubes.h" + +#ifdef _OPENMP + #undef match + #include +#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; +} diff --git a/src/CCubes.h b/src/CCubes.h new file mode 100755 index 0000000..8e1bde5 --- /dev/null +++ b/src/CCubes.h @@ -0,0 +1,16 @@ + +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 +); diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..0ee3615 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CFLAGS) diff --git a/ccubes.cl b/src/ccubes.cl similarity index 98% rename from ccubes.cl rename to src/ccubes.cl index cf0f1de..d62125c 100644 --- a/ccubes.cl +++ b/src/ccubes.cl @@ -109,9 +109,9 @@ nchoosek(int n, int k) #endif __kernel void ccubes_task(int k, - __global const real *nofvalues, /* IN: RC */ - __global const real *ON_set, /* IN: RC */ - __global const real *OFF_set, /* IN: RC */ + __global const int *nofvalues, /* IN: RC */ + __global const int *ON_set, /* IN: RC */ + __global const int *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 */ diff --git a/cl_setup.c b/src/cl_setup.c similarity index 98% rename from cl_setup.c rename to src/cl_setup.c index b5e52c9..cbea775 100644 --- a/cl_setup.c +++ b/src/cl_setup.c @@ -284,6 +284,7 @@ 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; @@ -326,7 +327,11 @@ cl_build(struct cl_uctx uctx, cl_device_type dev, result = CL_INVALID_VALUE; goto err; } - fread(kern_src, 1, srcsz, kern_file); + ret = fread(kern_src, 1, srcsz, kern_file); + if (ret != srcsz) { + log_warn("cl", "fread() failed!"); + goto err; + } kern_src[srcsz] = 0; log_info("cl", "FILE DUMP BEGINS"); diff --git a/cl_setup.h b/src/cl_setup.h similarity index 100% rename from cl_setup.h rename to src/cl_setup.h diff --git a/clccubes.c b/src/clccubes.c similarity index 95% rename from clccubes.c rename to src/clccubes.c index 92c631c..f13453a 100644 --- a/clccubes.c +++ b/src/clccubes.c @@ -133,9 +133,9 @@ err: int ccubes_alloc(struct ccubes_context *ctx, - real *nofvalues, /* IN: RC */ - real *ON_set, /* IN: RC */ - real *OFF_set, /* IN: RC */ + int *nofvalues, /* IN: RC */ + int *ON_set, /* IN: RC */ + int *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 real *nofvalues, IN: RC */ + /* __global const int *nofvalues, IN: RC */ ctx->nofvalues = clCreateBuffer(ctx->clctx->ctx, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - ctx->ninputs * sizeof(real), nofvalues, &rc); + ctx->ninputs * sizeof(int), nofvalues, &rc); if (rc != CL_SUCCESS) { goto err; } - /* __global const real *ON_set, IN: RC */ + /* __global const int *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(real), ON_set, &rc); + ctx->posrows * ctx->ninputs * sizeof(int), ON_set, &rc); if (rc != CL_SUCCESS) { goto err; } - /* __global const real *OFF_set, IN: RC */ + /* __global const int *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(real), OFF_set, &rc); + ctx->ninputs * ctx->negrows * sizeof(int), OFF_set, &rc); if (rc != CL_SUCCESS) { goto err; } @@ -381,7 +381,7 @@ err: } int -ccubes(int n_tasks, +ccubes_do_tasks(int n_tasks, int n_tasks_off, int k, int ninputs, @@ -391,9 +391,9 @@ ccubes(int n_tasks, int value_bit_width, int pichart_words, int estimPI, - real *nofvalues, /* IN: RC */ - real *ON_set, /* IN: RC */ - real *OFF_set, /* IN: RC */ + int *nofvalues, /* IN: RC */ + int *ON_set, /* IN: RC */ + int *OFF_set, /* IN: RC */ unsigned int *p_implicants_pos, /* IN: RC */ unsigned int *p_implicants_val, /* IN: RC */ int *last_index, /* IN: RC */ diff --git a/clccubes.h b/src/clccubes.h similarity index 82% rename from clccubes.h rename to src/clccubes.h index 1be063f..6c72361 100644 --- a/clccubes.h +++ b/src/clccubes.h @@ -70,10 +70,11 @@ struct ccubes_context { /* ND-Range */ size_t gws; /* global work size */ + size_t goff; /* global offset */ }; int -ccubes(int n_tasks, +ccubes_do_tasks(int n_tasks, int n_tasks_off, int k, int ninputs, @@ -83,9 +84,9 @@ ccubes(int n_tasks, int value_bit_width, int pichart_words, int estimPI, - real *nofvalues, /* IN: RC */ - real *ON_set, /* IN: RC */ - real *OFF_set, /* IN: RC */ + int *nofvalues, /* IN: RC */ + int *ON_set, /* IN: RC */ + int *OFF_set, /* IN: RC */ unsigned int *p_implicants_pos, /* IN: RC */ unsigned int *p_implicants_val, /* IN: RC */ int *last_index, /* IN: RC */ @@ -97,9 +98,4 @@ ccubes(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 diff --git a/config.c b/src/config.c similarity index 100% rename from config.c rename to src/config.c diff --git a/config.h b/src/config.h similarity index 100% rename from config.h rename to src/config.h diff --git a/logging.c b/src/logging.c similarity index 100% rename from logging.c rename to src/logging.c diff --git a/logging.h b/src/logging.h similarity index 100% rename from logging.h rename to src/logging.h diff --git a/queue.h b/src/queue.h similarity index 100% rename from queue.h rename to src/queue.h diff --git a/real.h b/src/real.h similarity index 100% rename from real.h rename to src/real.h diff --git a/src/registerDynamicSymbol.c b/src/registerDynamicSymbol.c new file mode 100644 index 0000000..bd87209 --- /dev/null +++ b/src/registerDynamicSymbol.c @@ -0,0 +1,8 @@ +#include +#include +#include + +void R_init_IEEE(DllInfo* info) { + R_registerRoutines(info, NULL, NULL, NULL, NULL); + R_useDynamicSymbols(info, TRUE); +} diff --git a/src/row_dominance.c b/src/row_dominance.c new file mode 100644 index 0000000..0b2c93a --- /dev/null +++ b/src/row_dominance.c @@ -0,0 +1,100 @@ + +#include +#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; +} diff --git a/src/row_dominance.h b/src/row_dominance.h new file mode 100755 index 0000000..0ac39fb --- /dev/null +++ b/src/row_dominance.h @@ -0,0 +1,10 @@ +#include + +void row_dominance( + int p_pichart[], + int p_implicants[], + int *p_ck, + int pirows, + unsigned int *foundPI, + int ninputs +); diff --git a/tree.h b/src/tree.h similarity index 100% rename from tree.h rename to src/tree.h