ccubes-cl/src/CCubes.c

#include <R.h>
#include <R_ext/RS.h>
#include <R_ext/Boolean.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"

#include "real.h"
#include "cl_setup.h"

#include "clccubes.h"

#include "config.h"
#include "logging.h"

// $ export BITS_PER_WORD=32 in the Terminal, before running R

SEXP CCubes(SEXP tt) {

    // simulate the R command:
    // system.file("ccubes.cl", package = "IEEE")

    const char *kernel_file = NULL;
    SEXP args = PROTECT(allocVector(VECSXP, 2));
    SET_VECTOR_ELT(args, 0, mkString("ccubes.cl"));
    SET_VECTOR_ELT(args, 1, mkString("IEEE"));

    SEXP names = PROTECT(allocVector(STRSXP, 2));
    SET_STRING_ELT(names, 1, mkChar("package"));
    Rf_setAttrib(args, R_NamesSymbol, names);

    SEXP r_result = PROTECT(R_tryEval(
        Rf_lang3(install("do.call"), mkString("system.file"), args),
        R_GlobalEnv,
        NULL
    ));

    Rprintf("Kernel file: '%s'\n", CHAR(STRING_ELT(r_result, 0)));
    kernel_file = CHAR(STRING_ELT(r_result, 0));
    UNPROTECT(3);

    // 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=true in the Terminal, before running R
    char *print_info = getenv("PRINT_INFO"); // Read from the PATH
    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;
    double elapsed_time;

    config_set_int("log", LOG_LEVEL_WARN);
    config_set_int("log:clccubes", LOG_LEVEL_DEBUG);
    config_set_int("log:ccubes", LOG_LEVEL_DEBUG);
    config_set_int("log:cl", LOG_LEVEL_DEBUG);

    char log[100] = {0};
    sprintf(log, "log-ccubes");
    config_set_string("out", log);

    struct ccubes_context *ctx = ccubes_create(kernel_file);
    if (ctx == NULL) {
        log_error("ccubes", "ccubes_do_tasks failed");
    }


    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 bits_needed = ceil(log2(max_value)); // Compute the necessary bits
    int value_bit_width = 1;
    while (value_bit_width < bits_needed) {
        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];

    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 = (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 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
    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 *) before R_Calloc() prefills all values with 0s

    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);
    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 previous (level k - 1), minimum number of PIs to solve the coverage matrix
    int solmin = 0;

    // the positions of the PIs solving the coverage matrix
    // 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;
        last_index[i] = 0;
    }

    Rboolean ON_set_covered = false;
    if (PRINT_INFO) {
        Rprintf("ON-set minterms: %d\n", posrows);
    }


    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);
        }
	log_debug("ccubes", "---k: %d\n", k);

	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 */
		int current_batch = n_tasks < n_tasks_batch ? n_tasks : n_tasks_batch;

		log_debug("ccubes", "Tasks %d - %d out of %d",
		    task, task + current_batch - 1, n_tasks);

		bool *redundant;
		bool *coverage;
		unsigned int *fixed_bits;
		unsigned int *value_bits;
		unsigned int *pichart_values;
		ccubes_do_tasks(ctx,
		    current_batch,
		    task,
		    k,
		    ninputs,
		    posrows,
		    negrows,
		    implicant_words,
		    value_bit_width,
		    value_bit_mask,
		    pichart_words,
		    estimPI,
		    nofvalues,
		    ON_set,
		    OFF_set,
		    p_implicants_pos,
		    p_implicants_val,
		    last_index,
		    p_covered,
		    p_pichart_pos,
		    redundant,
		    coverage,
		    fixed_bits,
		    value_bits,
		    pichart_values
			);

		for (int current_task = 0; current_task < current_batch; current_task++) {
			log_debug("ccubes", "Task %d", current_task);

			if (!ctx->h_redundant[current_task]) {
				/* LOG TASK */
				log_debug_raw("ccubes", "redundant[%d]: %d\n", current_task, ctx->h_redundant[current_task]);

				log_debug_raw("ccubes", "coverage[%d]:", current_task);
				for (int j = 0; j < posrows; j++) {
					log_debug_raw("ccubes", " %d",
					    ctx->h_coverage[current_task * posrows + j]);
				}
				log_debug_raw("ccubes", "\n");

				log_debug_raw("ccubes", "fixed_bits[%d]:", current_task);
				for (int j = 0; j < implicant_words; j++) {
					log_debug_raw("ccubes", " %d",
					    ctx->h_fixed_bits[current_task * implicant_words + j]);
				}
				log_debug_raw("ccubes", "\n");

				log_debug_raw("ccubes", "value_bits[%d]:", current_task);
				for (int j = 0; j < implicant_words; j++) {
					log_debug_raw("ccubes", " %d",
					    ctx->h_value_bits[current_task * implicant_words + j]);
				}
				log_debug_raw("ccubes", "\n");

				log_debug_raw("ccubes", "pichart_values[%d]:", current_task);
				for (int j = 0; j < pichart_words; j++) {
					log_debug_raw("ccubes", " %d",
					    ctx->h_pichart_values[current_task * pichart_words + j]);
				}
				log_debug_raw("ccubes", "\n");


				int covsum = 0;
				for (int i = 0; i < posrows; i++) {
					covsum += ctx->h_coverage[current_task * posrows + i];
				}
				// 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] = ctx->h_fixed_bits[current_task * implicant_words + w];
				    p_implicants_val[implicant_words * foundPI + w] = ctx->h_value_bits[current_task * implicant_words + w];
				}

				// 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] = ctx->h_pichart_values[current_task * pichart_words + w];
				    }

				    p_pichart[posrows * foundPI + r] = ctx->h_coverage[current_task * posrows + r];
				}

				++foundPI;

				// when needed, increase allocated memory
				if (foundPI / estimPI > 0.9) {
				    estimPI += 100000;
				    p_pichart =        R_Realloc(p_pichart,        posrows * estimPI,         int);
				    p_pichart_pos =    R_Realloc(p_pichart_pos,    pichart_words * estimPI,   unsigned int);
				    p_implicants_val = R_Realloc(p_implicants_val, implicant_words * 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 = foundPI; i < posrows * estimPI; i++) {
					p_pichart[i] = 0;
				    }
				    for (unsigned int i = foundPI; i < pichart_words * estimPI; i++) {
					p_pichart_pos[i] = 0U;
				    }
				    for (unsigned int i = foundPI; i < implicant_words * estimPI; i++) {
					p_implicants_val[i] = 0U;
					p_implicants_pos[i] = 0U;
				    }

				    if (PRINT_INFO) {
					multiplier++;
					Rprintf("%dx ", multiplier);
				    }
				}
			}
		}

		/* change to something less aggresive for reuse */
		ccubes_clean_up(ctx);
	}

        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 coverage matrix: %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;
        }
    }

    // Rprintf("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 / value_bit_width);                     // Word index within the implicant
            int bit_index = (r % (BITS_PER_WORD / value_bit_width)) * value_bit_width;  // Bit position within the word

                                                                           // (1U << bit_index) is just as good
            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
        }
    }

    ccubes_destroy(ctx);

    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 > n) return 0;
    if (k == 0 || k == n) return 1;

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

        // 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;
}