/** Methods computing or estimating the second moment of the Ising model given * the natural parameters of the p.m.f. in exponential family form. That is * f(x) = p0(params) exp(vech(x x')' params) * where `params` are the natural parameters. The scaling constant `p0(params)` * is also the zero event `x = (0, ..., 0)` probability. * * Three different version are provided. * - exact: Exact method computing the true second moment `E[X X']` given the * parameters `params`. This method has a runtime of `O(2^dim)` where * `dim` is the dimension of the random variable `X`. This means it is * infeasible for bigger (like 25 and above) values if `dim` * - MC: Monte-Carlo method used in case the exact method is not feasible. * - MC Thrd: Multi-Threaded Monte-Carlo method. * * The natural parameters `params` of the Ising model is a `dim (dim + 1) / 2` * dimensional numeric vector containing the log-odds. The indexing schema into * the parameters corresponding to the log-odds of single or two way interactions * with indicex `i, j` are according to the half-vectorization schema * * i = 4'th row => Symmetry => Half-Vectorized Storage * * [ * * * * * * ] [ * * * * * * ] [ * * * 15 * * ] * [ * * * * * * ] [ * * * * * * ] [ * * 12 16 * ] * [ * * * * * * ] [ * * * * * * ] [ * 8 * 17 ] * [ 3 8 12 15 16 17 ] [ 3 8 12 15 * * ] [ 3 * * ] * [ * * * * * * ] [ * * * 16 * * ] [ * * ] * [ * * * * * * ] [ * * * 17 * * ] [ * ] * */ #include // DBL_MAX #include "R_api.h" #include "bit_utils.h" #include "int_utils.h" #include "ising_MCMC.h" #ifndef __STDC_NO_THREADS__ #include #endif //////////////////////////////////////////////////////////////////////////////// // TODO: read // // https://developer.r-project.org/Blog/public/2019/04/18/common-protect-errors/ // and // https://developer.r-project.org/Blog/public/2018/12/10/unprotecting-by-value/index.html // as well as do the following PROPERLLY // TODO: make specialized version using the parameters in given sym mat form // // similar to `ising_sample()` // //////////////////////////////////////////////////////////////////////////////// // void ising_m2_exact_mat( // const size_t dim, const size_t ldA, const double* A, // double* M2 // ) { // double sum_0 = 1.0; // const uint32_t max_event = (uint32_t)1 << dim; // (void)memset(M2, 0, dim * dim * sizeof(double)); // for (uint32_t X = 1; X < max_event; ++X) { // // Evaluate quadratic form `X' A X` using symmetry of `A` // double XtAX = 0.0; // for (uint32_t Y = X; Y; Y &= Y - 1) { // const int i = bitScanLS32(Y); // XtAX += A[i * (ldA + 1)]; // for (uint32_t Z = Y & (Y - 1); Z; Z &= Z - 1) { // XtAX += 2 * A[i * ldA + bitScanLS32(Z)]; // } // } // const double prob_X = exp(XtAX); // sum_0 += prob_X; // for (uint32_t Y = X; Y; Y &= Y - 1) { // const int i = bitScanLS32(Y); // for (uint32_t Z = Y; Z; Z &= Z - 1) { // M2[i * dim + bitScanLS32(Z)] += prob_X; // } // } // } // const double prob_0 = 1.0 / sum_0; // for (size_t j = 0; j < dim; ++j) { // for (size_t i = 0; i < j; ++i) { // M2[j * dim + i] = M2[i * dim + j]; // } // for (size_t i = j; i < dim; ++i) { // M2[i * dim + j] *= prob_0; // } // } // } double ising_m2_exact(const size_t dim, const double* params, double* M2) { // Accumulator of sum of all (unscaled) probabilities double sum_0 = 1.0; // max event (actually upper bound) const uint32_t max_event = (uint32_t)1 << dim; // Length of parameters const size_t len = dim * (dim + 1) / 2; // Initialize M2 to zero (void)memset(M2, 0, len * sizeof(double)); // Iterate all `2^dim` binary vectors `X` for (uint32_t X = 1; X < max_event; ++X) { // Dot product `` double dot_X = 0.0; // all bits in `X` for (uint32_t Y = X; Y; Y &= Y - 1) { const int i = bitScanLS32(Y); const int I = (i * (2 * dim - 1 - i)) / 2; // add single and two way interaction log odds for (uint32_t Z = Y; Z; Z &= Z - 1) { dot_X += params[I + bitScanLS32(Z)]; } } // (unscaled) probability of current event `X` const double prob_X = exp(dot_X); sum_0 += prob_X; // Accumulate set bits probability for the first and second moment `E[X X']` for (uint32_t Y = X; Y; Y &= Y - 1) { const int i = bitScanLS32(Y); const int I = (i * (2 * dim - 1 - i)) / 2; for (uint32_t Z = Y; Z; Z &= Z - 1) { M2[I + bitScanLS32(Z)] += prob_X; } } } // Finish by scaling with zero event probability (Ising p.m.f. scaling constant) const double prob_0 = 1.0 / sum_0; for (size_t i = 0; i < len; ++i) { M2[i] *= prob_0; } return log(prob_0); } // Monte-Carlo method using Gibbs sampling for the second moment double ising_m2_MC( /* options */ const size_t nr_samples, const size_t warmup, /* dimension */ const size_t dim, /* vectors */ const double* params, double* M2 ) { // Length of `params` and `M2` const size_t len = (dim * (dim + 1)) / 2; // Dirty trick (reuse output memory for precise counting) _Static_assert(sizeof(uint64_t) == sizeof(double), "Dirty trick fails!"); uint64_t* counts = (uint64_t*)M2; // Initialize counts to zero (void)memset(counts, 0, len * sizeof(uint64_t)); // Allocate memory to store Markov-Chain value int* X = (int*)R_alloc(dim, sizeof(int)); // Accumulator for Monte-Carlo estimate for zero probability `P(X = 0)` double accum = 0.0, max_mdot_X = -DBL_MAX; // Create/Update R's internal PRGN state GetRNGstate(); // Spawn chains for (size_t sample = 0; sample < nr_samples; ++sample) { // Draw random sample `X ~ Ising(params)` ising_mcmc_vech(warmup, dim, params, X); // Accumulate component counts and compute (unscaled) probability of `X` uint64_t* count = counts; double dot_X = 0.0; for (size_t i = 0; i < dim; ++i) { const size_t I = (i * (2 * dim - 1 - i)) / 2; for (size_t j = i; j < dim; ++j) { *(count++) += X[i] & X[j]; dot_X += (X[i] & X[j]) ? params[I + j] : 0.0; } } if (-dot_X > max_mdot_X) { accum *= exp(max_mdot_X + dot_X); max_mdot_X = -dot_X; } accum += exp(-dot_X - max_mdot_X); } // Write R's internal PRNG state back (if needed) PutRNGstate(); // Compute means from counts for (size_t i = 0; i < len; ++i) { M2[i] = (double)counts[i] / (double)nr_samples; } // Prob. of zero event (Ising p.m.f. scaling constant) return log(accum) + max_mdot_X - log(2) * (double)dim - log((double)nr_samples); } // in case the compile supports standard C threads #ifndef __STDC_NO_THREADS__ // Worker thread to calling thread communication struct typedef struct thrd_data { rng_seed_t seed; // Pseudo-Random-Number-Generator seed/state value size_t nr_samples; // Nr. of samples this worker handles size_t warmup; // Monte-Carlo Chain burne-in length size_t dim; // Random variable dimension const double* params; // Ising model parameters int* X; // Working memory to store current binary sample uint32_t* counts; // (output) count of single and two way interactions double log_sum_exp; // (output) Monte-Carlo inbetween value for `log(P(X = 0))` } thrd_data_t; // Worker thread function int thrd_worker(thrd_data_t* data) { // Extract data as thread local variables (for convenience) const size_t dim = data->dim; int* X = data->X; // Initialize counts to zero (void)memset(data->counts, 0, (dim * (dim + 1) / 2) * sizeof(uint32_t)); // Accumulator for Monte-Carlo estimate for zero probability `P(X = 0)` double accum = 0.0, max_mdot_X = -DBL_MAX; // Spawn Monte-Carlo Chains (one foe every sample) for (size_t sample = 0; sample < data->nr_samples; ++sample) { // Draw random sample `X ~ Ising(params)` ising_mcmc_vech_thrd(data->warmup, dim, data->params, X, data->seed); // Accumulate component counts and compute (unscaled) probability of `X` uint32_t* count = data->counts; double dot_X = 0.0; for (size_t i = 0; i < dim; ++i) { const size_t I = (i * (2 * dim - 1 - i)) / 2; for (size_t j = i; j < dim; ++j) { *(count++) += X[i] & X[j]; dot_X += (X[i] & X[j]) ? data->params[I + j] : 0.0; } } if (-dot_X > max_mdot_X) { accum *= exp(max_mdot_X + dot_X); max_mdot_X = -dot_X; } accum += exp(-dot_X - max_mdot_X); } data->log_sum_exp = log(accum) + max_mdot_X; return 0; } double ising_m2_MC_thrd( /* options */ const size_t nr_samples, const size_t warmup, const size_t nr_threads, /* dimension */ const size_t dim, /* vectors */ const double* params, double* M2 ) { // Length of `params` and `M2` const size_t len = (dim * (dim + 1)) / 2; // Allocate working and self memory for worker threads int* Xs = (int*)R_alloc(dim * nr_threads, sizeof(int)); uint32_t* counts = (uint32_t*)R_alloc(len * nr_threads, sizeof(uint32_t)); thrd_t* threads = (thrd_t*)R_alloc(nr_threads, sizeof(thrd_data_t)); thrd_data_t* threads_data = (thrd_data_t*)R_alloc(nr_threads, sizeof(thrd_data_t)); // Provide instruction wor worker and dispatch them for (size_t tid = 0; tid < nr_threads; ++tid) { // Every thread needs its own PRNG seed! init_seed(threads_data[tid].seed); // divide work among workers (more or less) equaly with the first worker // (tid = 0) having the additional remainder of divided work (`nr_samples`) threads_data[tid].nr_samples = nr_samples / nr_threads + ( tid ? 0 : nr_samples % nr_threads ); threads_data[tid].warmup = warmup; threads_data[tid].dim = dim; threads_data[tid].params = params; // Every worker gets disjoint working memory threads_data[tid].X = &Xs[dim * tid]; threads_data[tid].counts = &counts[len * tid]; // dispatch the worker thrd_create(&threads[tid], (thrd_start_t)thrd_worker, (void*)&threads_data[tid]); } // Join (Wait for) all worker for (size_t tid = 0; tid < nr_threads; ++tid) { thrd_join(threads[tid], NULL); } // Accumulate worker results into first (tid = 0) worker counts result // while determining the log sum exp maximum double max = threads_data[0].log_sum_exp; for (size_t tid = 1; tid < nr_threads; ++tid) { for (size_t i = 0; i < len; ++i) { counts[i] += counts[tid * len + i]; } max = (max < threads_data[tid].log_sum_exp) ? threads_data[tid].log_sum_exp : max; } // Accum all `log(P(X = 0))` via "LogSumExp" trick into final result double accum = 0; for (size_t tid = 0; tid < nr_threads; ++tid) { accum += exp(threads_data[tid].log_sum_exp - max); } // convert discreat counts into means for (size_t i = 0; i < len; ++i) { M2[i] = (double)counts[i] / (double)nr_samples; } // Log of Prob. of zero event (Ising p.m.f. scaling constant) return log(accum) + max - log(2) * (double)dim - log((double)nr_samples); } #endif /* !__STDC_NO_THREADS__ */ /* Implementation of `.Call` entry point */ //////////////////////////////////////////////////////////////////////////////// // TODO: make specialized version using the parameters in given sym mat form // // similar to `ising_sample()` // //////////////////////////////////////////////////////////////////////////////// extern SEXP R_ising_m2( SEXP _params, SEXP _use_MC, SEXP _nr_samples, SEXP _warmup, SEXP _nr_threads ) { // Extract and validate arguments size_t protect_count = 0; if (!Rf_isReal(_params)) { _params = PROTECT(Rf_coerceVector(_params, REALSXP)); protect_count++; } // Depending on parameters given in symmetric matrix form or natural // or natrual parameters in the half-vectorized memory layout size_t dim; double* params; if (Rf_isMatrix(_params)) { if (Rf_nrows(_params) != Rf_ncols(_params)) { Rf_error("Invalid 'params' value, exected square matrix"); } // Convert to natural parameters dim = Rf_nrows(_params); params = (double*)R_alloc(dim * (dim + 1) / 2, sizeof(double)); double* A = REAL(_params); for (size_t j = 0, I = 0; j < dim; ++j) { for (size_t i = j; i < dim; ++i, ++I) { params[I] = (1 + (i != j)) * A[i + j * dim]; } } } else { // Get Ising random variable dimension `dim` from "half-vectorized" parameters // vector. That is, compute `dim` from `length(params) = dim (dim + 1) / 2`. dim = invTriag(Rf_length(_params)); if (!dim) { Rf_error("Expected parameter vector of length `p (p + 1) / 2` where" " `p` is the dimension of the random variable"); } params = REAL(_params); } // Determin the method to use and validate if possible const int use_MC = Rf_asLogical(_use_MC); if (use_MC == NA_LOGICAL) { Rf_error("Invalid 'use_MC' value, expected ether TRUE or FALSE"); } // Allocate result vector SEXP _M2 = PROTECT(Rf_allocVector(REALSXP, dim * (dim + 1) / 2)); ++protect_count; // asside computed log of zero event probability (log of recibrocal partition // function), the log scaling factor for the Ising model p.m.f. double log_prob_0 = -1.0; if (use_MC) { // Convert and validate arguments for the Monte-Carlo methods const size_t nr_samples = asUnsigned(_nr_samples); if (nr_samples == 0 || nr_samples == NA_UNSIGNED) { Rf_error("Invalid 'nr_samples' value, expected pos. integer"); } const size_t warmup = asUnsigned(_warmup); if (warmup == NA_UNSIGNED) { Rf_error("Invalid 'warmup' value, expected non-negative integer"); } const size_t nr_threads = asUnsigned(_nr_threads); if (nr_threads == 0 || nr_threads > 256) { Rf_error("Invalid 'nr_thread' value"); } if (nr_threads == 1) { // Single threaded Monte-Carlo method log_prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); } else { // Multi-Threaded Monte-Carlo method if provided, otherwise use // the single threaded version with a warning #ifdef __STDC_NO_THREADS__ Rf_warning("Multi-Threading NOT supported, using fallback."); log_prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); #else log_prob_0 = ising_m2_MC_thrd( nr_samples, warmup, nr_threads, dim, params, REAL(_M2) ); #endif } } else { // Exact method (ignore other arguments), only validate dimension if (25 < dim) { Rf_error("Dimension '%d' too big for exact method (max 24)", dim); } // and call the exact method log_prob_0 = ising_m2_exact(dim, params, REAL(_M2)); } // Set log-lokelihood as an attribute to the computed second moment SEXP _log_prob_0 = PROTECT(Rf_ScalarReal(log_prob_0)); ++protect_count; Rf_setAttrib(_M2, Rf_install("log_prob_0"), _log_prob_0); // release SEPXs to the garbage collector UNPROTECT(protect_count); return _M2; }