//  apply.c   "PSEUDO-CRYSTAL" NESTED SAMPLING APPLICATION
// (GNU General Public License software, (C) Sivia and Skilling 2006)
// Problem:     M switches s = 0 or 1, grouped in clusters of widths h.
//                 e.g. M=10, s = {0,0,0,1,1,1,1,0,0,1}
//                            h = {  3  ,   4   , 2 ,1}
// Inputs:
//  Prior(s)    is uniform, 1/2^N on each of 2^M states
//  Likelihood  is L(s) = exp( SUM h(h-1)/M )
// Outputs:
//  Evidence    is Z = SUM L(s) Prior(s)
//  Posterior   is P(s) = L(s) Prior(s)/ Z
//  Information is H = SUM P(s) log(P(s)/Prior(s))
/*__________________________________________________________________*/
#define n     1    // # Objects
#define MAX 800    // # iterates
#define M  1000    // # switches in this application
/*__________________________________________________________________*/
typedef struct
{
    char    s[M];  // state of switches
    double  logL;  // logLikelihood = ln prob(data | s)
    double  logWt; // log(Weight), adding to SUM(Wt) = Evidence Z
} Object;
/*__________________________________________________________________*/
double logLhood(   // logLikelihood function
char*  s)          // switches
{
    int    i, j;      // left and right counters
    double logL = 0;  // logLikelihood accumulator
    i = 0;                                // L.H. boundary
    for( j = 1; j < M; j++ )
        if( s[j] != s[j-1] )
        {                                 // R.H. boundary found
           logL += (j - i) * (j - i - 1); // cluster width h = j-i
           i = j;                         // reset L.H. boundary
        }
    logL += (j - i) * (j - i - 1);        // R.H. cluster
    return logL / M + sqrt(DBL_EPSILON) * UNIFORM;  // normalised
}                // jitter eliminates ties between likelihood values
/*__________________________________________________________________*/
void Prior(        // Set Object according to prior
Object* Obj)       // Object being set
{
    int  j;
    for( j = 0; j < M; j++ )
        Obj->s[j] = (int)(2 * UNIFORM) % 2;   // 0 or 1
    Obj->logL = logLhood(Obj->s);
}
/*__________________________________________________________________*/
void Explore(      // Evolve object within likelihood constraint
Object* Obj,       // Object being evolved
double  logLstar)  // Likelihood constraint L > Lstar
{
    int     m = 10 * M;   // MCMC counter (pre-judged # steps)
    int     try;          // Try flipping this switch
    double  logLtry;      // Trial loglikelihood
    for( ; m > 0; m-- )
    {
        try = (int)(M * UNIFORM) % M;      // random switch
        Obj->s[try] = 1 - Obj->s[try];     // try flipping
        logLtry = logLhood(Obj->s);        // trial loglikelihood
        if( logLtry > logLstar )
            Obj->logL = logLtry;           // accept
        else
            Obj->s[try] = 1 - Obj->s[try]; // reject
    }
}
/*__________________________________________________________________*/
void Results(      // Output nested sampling sequence
Object* Samples,   // Objects defining posterior
int     nest,      // # Samples
double  logZ)      // Evidence (= total weight = SUM[Samples] Weight)
{
    int  k;        // Sample counter
    for( k = 0; k < nest; k++ )
        printf("%7.2f %8.4f\n", -(k+1.) / n, Samples[k].logL);
}            // print log(enclosed prior mass), log(likelihood)