// apply.c  "ON/OFF SWITCHING" NESTED SAMPLING APPLICATION
// (GNU General Public License software, (C) Sivia and Skilling 2006)
// Problem:     5 switches s, with probabilistic data.
// Inputs:
//  Prior(s)    is {0.1, 0.2, 0.3, 0.2, 0.1} for the 5 switches.
//  Likelihood  is Pr(s[2] OR s[4]) = 0.2,  Pr(s[3] OR s[4]) = 0.9,
//                 Pr(s[1] AND s[4]) = 0.8.
// 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    100   // # Objects
#define MAX 1000   // # iterates
#define M      5   // # switches
#define Mplus  8   // power-of-2 >= M
#include "tree.c"  // binary tree procedures void PutRate, int GetRate
const double PrON[M] = {0.1, 0.2, 0.3, 0.2, 0.1};  // Prior(s is ON)
/*__________________________________________________________________*/
typedef struct
{
    char    s[M];          // switches
    double  Tree[2*Mplus]; // binary rate tree
    double  logL;          // logLikelihood = ln Prob(data | s)
    double  logWt;         // log(Weight), with SUM(Wt) = Evidence Z
} Object;
/*__________________________________________________________________*/
double logLhood(        // logLikelihood function
char*  s)               // switches
{
    double  L = 1.0;
    L *= (s[2] | s[4]) ? 0.2 : 0.8;
    L *= (s[3] | s[4]) ? 0.9 : 0.1;
    L *= (s[1] & s[4]) ? 0.8 : 0.2;
    return  log(L) + sqrt(DBL_EPSILON) * UNIFORM;
}                // jitter eliminates ties between likelihood values
/*__________________________________________________________________*/
void Prior(             // Set Object according to prior
Object* Obj)            // Object being set
{
    int  i;
// Initialise empty Tree of transition rates
    Obj->Tree[0] = Mplus;
    for( i = 1; i < 2*Mplus; i++ )
        Obj->Tree[i] = 0.0;
// Initialise object
    for( i = 0; i < M; i++ )
    {
        Obj->s[i] = (UNIFORM < PrON[i]);
        PutRate(Obj->Tree, i, Obj->s[i] ? 1.-PrON[i] : PrON[i]);
    }
    Obj->logL = logLhood(Obj->s);
}
/*__________________________________________________________________*/
void Explore(           // Evolve object within likelihood constraint
Object* Obj,            // Object being evolved
double  logLstar)       // Likelihood constraint L > Lstar
{
    double Interval = 30.0; // pre-judged time to equilibrate
    double t = 0.0;         // evolution time, initialized 0
    double logLtry;         // trial loglikelihood
    int    i;               // switch being flipped
    while( (t += -log(UNIFORM) / Obj->Tree[1]) < Interval )
    {
        i = GetRate(Obj->Tree);
        Obj->s[i] = 1 - Obj->s[i];          // trial state
        logLtry = logLhood(Obj->s);
        if( logLtry > logLstar )
        {                                   // accept
            Obj->logL = logLtry;
            PutRate(Obj->Tree, i,
                     Obj->s[i] ? 1.-PrON[i] : PrON[i]);
        }
        else
            Obj->s[i] = 1 - Obj->s[i];      // reject
    }
}
/*__________________________________________________________________*/
void Results(     // Posterior properties, here mean of s
Object* Samples,  // Objects defining posterior
int     nest,     // # Samples
double  logZ)     // Evidence (= total weight = SUM[Samples] Weight)
{
    double post[M] = {0,0,0,0,0}; // posterior prob(ON)
    double w;                     // Proportional weight
    int    i;                     // Sequence counter
    int    k;                     // Sample counter
    for( i = 0; i < nest; i++ )
    {
        w = exp(Samples[i].logWt - logZ);
        for( k = 0; k < M; k++ )
            if( Samples[i].s[k] )  post[k] += w;
    }
    for( k = 0; k < M; k++ )
        printf("%d: Posterior(ON)=%3.0f%%\n", k, 100.*post[k]);
}