#ifndef INCLUDE_GUARD_ISING_MCMC_H
#define INCLUDE_GUARD_ISING_MCMC_H

#include "R_api.h"

/** Sample a single value from the Ising model via a Monte-Carlo Markov-Chain
 * method given the Ising model parameters in compact half-vectorized form.
 *
 *      f(x) = p0(params) exp(vech(x x')' params)
 *
 * @important requires to be wrapped between `GetRNGstate()` and `PutRNGState()`
 * calls. For example
 *
 * GetRNGstate()
 *  for (size_t sample = 0; sample < nr_samples; ++sample) {
 *      ising_mcmc_vech(warmup, dim, dim, params, &X[sample * dim]);
 *  }
 * PutRNGState()
 */
static inline void ising_mcmc_vech(
    const size_t warmup,
    const size_t dim, const double* params,
    int* X
) {
    // Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
    for (size_t i = 0, I = 0; i < dim; I += dim - i++) {
        double invProb = 1.0 + exp(-params[I]);
        X[i] = unif_rand() * invProb < 1.0;
    }

    // Skip the first samples of the Markov-Chain (warmup)
    for (size_t skip = 0; skip < warmup; ++skip) {
        // For every component
        for (size_t i = 0; i < dim; ++i) {
            // Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
            double log_odds = 0.0;
            // Tracks position in half-vectorized storage
            size_t J = i;
            // Sum all log-odds for `j < i` (two way interactions)
            for (size_t j = 0; j < i; J += dim - ++j) {
                log_odds += X[j] ? params[J] : 0.0;
            }
            // Allways add log-odds for `i = j` (self-interaction)
            log_odds += params[J];
            // Continue adding log-odds for `j > i` (two way interactions)
            J++;
            for (size_t j = i + 1; j < dim; ++j, ++J) {
                log_odds += X[j] ? params[J] : 0.0;
            }
            // Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
            X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
        }
    }
}

// Thread save version of `ising_mcmc` (using thread save PRNG)
static inline void ising_mcmc_vech_thrd(
    const size_t warmup,
    const size_t dim, const double* params,
    int* X, rng_seed_t seed
) {
    // Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
    for (size_t i = 0, I = 0; i < dim; I += dim - i++) {
        double invProb = 1.0 + exp(-params[I]);
        X[i] = unif_rand_thrd(seed) * invProb < 1.0;
    }

    // Skip the first samples of the Markov-Chain (warmup)
    for (size_t skip = 0; skip < warmup; ++skip) {
        // For every component
        for (size_t i = 0; i < dim; ++i) {
            // Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
            double log_odds = 0.0;
            // Tracks position in half-vectorized storage
            size_t J = i;
            // Sum all log-odds for `j < i` (two way interactions)
            for (size_t j = 0; j < i; J += dim - ++j) {
                log_odds += X[j] ? params[J] : 0.0;
            }
            // Allways add log-odds for `i = j` (self-interaction)
            log_odds += params[J];
            // Continue adding log-odds for `j > i` (two way interactions)
            J++;
            for (size_t j = i + 1; j < dim; ++j, ++J) {
                log_odds += X[j] ? params[J] : 0.0;
            }
            // Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
            X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
        }
    }
}

/** Sample a single value from the Ising model via a Monte-Carlo Markov-Chain
 * method given the Ising model parameters in symmetric matrix form `A`.
 *
 *      f(x) = prob_0(A) exp(x' A x)
 *
 * The relation to the natural parameters `params` as used in the half-vectorized
 * version `ising_mcmc_vech()`
 *
 *      f(x) = prob_0(params) exp(vech(x x')' params)
 *
 * is illustrated via an example with `dim = 5`;
 *
 *          p0                             [ p0   p1/2 p2/2  p3/2  p4/2  ]
 *          p1 p5                          [ p1/2 p5   p6/2  p7/2  p8/2  ]
 * params = p2 p6 p9             =>    A = [ p2/2 p6/2 p9    p10/2 p11/2 ]
 *          p3 p7 p10 p12                  [ p3/2 p7/2 p10/2 p12   p13/2 ]
 *          p4 p8 p11 p13 p14              [ p4/2 p8/2 p11/2 p13/2 p14   ]
 *
 * or the other way arroung given the symmetric matrix form `A` converted to
 * the natural parameters
 *
 *     [ a00 a10 a20 a30 a40 ]                   a00
 *     [ a10 a11 a21 a31 a41 ]                   2*a10 a11
 * A = [ a20 a21 a22 a32 a42 ]    =>    params = 2*a20 2*a21 a22
 *     [ a30 a31 a32 a33 a43 ]                   2*a30 2*a31 2*a32 a33
 *     [ a40 a41 a42 a43 a44 ]                   2*a40 2*a41 2*a42 2*a43 a44
 *
 */
static inline void ising_mcmc_mat(
    const size_t warmup,
    const size_t dim, const size_t ldA, const double* A,
    int* X
) {
    // Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
    for (size_t i = 0; i < dim; ++i) {
        X[i] = unif_rand() * (1.0 + exp(-A[i * (ldA + 1)])) < 1.0;
    }

    // Skip the first samples of the Markov-Chain (warmup)
    for (size_t skip = 0; skip < warmup; ++skip) {
        // For every component
        for (size_t i = 0; i < dim; ++i) {
            // Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
            double log_odds = 0.0;
            for (size_t j = 0; j < i; ++j) {
                log_odds += (2 * X[j]) * A[i * ldA + j];
            }
            log_odds += A[i * (ldA + 1)];
            for (size_t j = i + 1; j < dim; ++j) {
                log_odds += (2 * X[j]) * A[i * ldA + j];
            }
            // Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
            X[i] = unif_rand() * (1.0 + exp(-log_odds)) < 1.0;
        }
    }
}

// thread save version of `ising_mcmc_mat()` (using thread save PRNG)
static inline void ising_mcmc_mat_thrd(
    const size_t warmup,
    const size_t dim, const size_t ldA, const double* A,
    int* X, rng_seed_t seed
) {
    // Initialize elements `X_i ~ Bernoulli(P(X_i = 1 | X_-i = 0))`
    for (size_t i = 0; i < dim; ++i) {
        X[i] = unif_rand_thrd(seed) * (1.0 + exp(-A[i * (ldA + 1)])) < 1.0;
    }

    // Skip the first samples of the Markov-Chain (warmup)
    for (size_t skip = 0; skip < warmup; ++skip) {
        // For every component
        for (size_t i = 0; i < dim; ++i) {
            // Compute conditional probability `P(X_i = 1 | X_-i = x_-i)`
            double log_odds = 0.0;
            for (size_t j = 0; j < i; ++j) {
                log_odds += (2 * X[j]) * A[i * ldA + j];
            }
            log_odds += A[i * (ldA + 1)];
            for (size_t j = i + 1; j < dim; ++j) {
                log_odds += (2 * X[j]) * A[i * ldA + j];
            }
            // Update `i`th element `X_i ~ Bernoulli(P(X_i = 1 | X_-i = x_-i))`
            X[i] = unif_rand_thrd(seed) * (1.0 + exp(-log_odds)) < 1.0;
        }
    }
}


#endif /* INCLUDE_GUARD_ISING_MCMC_H */