//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//             Bayesian Inference
// 
// Filename:   bayesys3.c
//
// Purpose:    Obtain sample objects from posterior atomic distribution.
//
// Dedication: To my intellectual ancestors  the late Edwin T Jaynes,  and
//             Steve Gull, to my descendant Sibusiso Sibisi,  and the many
//             colleagues and friends  over the  past quarter-century  who
//             have inspired and encouraged the development of these ideas.
//
//             John Skilling, Kenmare, Ireland, November 2003
//             email: skilling@eircom.net
//=============================================================================
/*
    Copyright (c) 1999-2003 Maximum Entropy Data Consultants Ltd,
                            114c Milton Road, Cambridge CB4 1XE, England

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include "license.txt"
*/
//=============================================================================
//
//  An OBJECT is a combination of a-priori-equivalent "atoms":
//  statisticians call this construction a "mixture model".
// 
//  ATOMS
//  The number of atoms N in an object can range from MIN to MAX.
//  MIN >= 1 is required to avoid the algorithmically special null object N=0.
//  Generally, MAX >= MIN is required (with MAX=MIN allowed), but
//  MAX=infinity (absence of limit) can be specified with MAX=0.
//  [ Technically, this coding allows at most 1 atom per identifiable
//    location, so "infinity" means #(locations), perhaps 2^32 . ]
//
//  There are three standard priors Pr(N), selected by sign(Alpha).
//  Alpha = 0.
//      The prior is UNIFORM in MIN <= N <= MAX (and MAX must be finite).
//  Alpha > 0.
//      The prior is POISSON or BINOMIAL.
//      N  ~  Alpha +- sqrt(Alpha), subject to MIN and MAX
//  Alpha < 0.
//      The prior is GEOMETRIC.
//      N  ~  |Alpha| +- |Alpha|, subject to MIN and MAX
//
//  COORDINATES
//  The BayeSys prior for locating an atom in space is flat over the hypercube
//  (0,1)^Ndim  of dimension  Ndim .  In fact, the hypercube is treated as
//  wraparound-continuous, so is more properly described as a "hyper-torus".
//  BayeSys gives you coordinate vectors "double Cube[Ndim]" with values in
//  (0,1), but you can transform these into other "Coord" if you wish.
//  
//  The BayeShape procedure lets you modify some or all of your coordinates
//  into nine alternative configurations:
//
//   Shape  Description           Prob(Coord[i])                    Range of i
//
//     0   Permutation         uniform on integers Perm(0...N-1)     0...N-1
//     1   +ve orthant         exp(-x) in x > 0                      0...N-1
//     2   Simplex volume      uniform in SUM(x) < 1                 0...N-1
//     3   Simplex surface     uniform on SUM(x) = 1                  0...N
//     4   Ordered             uniform in 0 < x[0] < x[1] < ... < 1  0...N-1
//     5   Bell                Normal(0,1) in -inf < x < inf         0...N-1
//     6   Sphere volume       uniform in SUM(x^2) < 1               0...N-1
//     7   Hemisphere surface  uniform on SUM(x^2) = 1, x[N]>0        0...N
//     8   Sphere surface      uniform on SUM(x^2) = 1                0...N
//
//  For the interior of an awkward shape, allow a circumscribing volume,
//  but return the value 0 from UserTry1 and UserTry2 (instead of the
//  usual positive acceptance code) for any location outside the shape.
//  The ensemble will keep within the domain, because BayeSys only accepts
//  points for which you give strictly positive return codes.
//
//  Thus the unit disc  x^2 + y^2 < 1 could be programmed in several ways:
// (a) Use BayeShape with Shape=6 which gives the N=2 disc interior directly.
// (b) Use N=2 hypercube (= unit square) but expand it to (-1,1) by setting
//        x = 2*Cube[0] - 1,  y = 2*Cube[1] - 1
//     and reject any point outside the disc by returning 0 from UserTry1
//     and UserTry2.
// (c) Use hypercube to yield a unit square, but spread the prior measure
//     uniformly through the unit disc with your own transformation, e.g.
//        radius^2 = Cube[0],  angle = 2 PI Cube[1] .
//
//  FLUXES
//  The Ndim coordinates can be supplemented by  Valency  intensities or
//  "fluxes" z, being attributes for which the joint distribution
//                  Prior(z).Likelihood(D|z,...)
//  is integrable and can be sampled from.
//  The MassInf library incorporated in BayeSys provides Flux procedures for
//  its priors and linear data.
//
//  DISPLAY 
//  The coordinates and intensities are supplemented by a guidance width, being
//      log(fraction of hypercube volume that atom might plausibly range over).
//  This may help you to produce smooth displays from the atomic objects.
//
//  The number of atoms in an object is supplied to you as Natoms and
//  each atom's attributes are supplied to you consecutively as
//    double Cube[ 0,1,...,Ndim-1,  Ndim,...,Ndim+Valency-1,  Ndim+Valency ]
//              coordinates in (0,1)     intensities            log(width)
//=============================================================================
//
//  ENSEMBLE
//  The program uses an ensemble of sample objects.  Using several objects is
//  usually recommended, partly because objects that seem to be getting stuck
//  are overwritten by more successful ones which reduces the risk of failure,
//  and partly because the geometrical exploration engines only work with
//  several objects.
//
//=============================================================================
//
//  METHOD
//  Internally in the program, hypercube coordinates are mapped to an
//  extended-integer label whose range fills the hypercube but which preserves
//  a degree of locality: small changes in this integer will necessarily
//  correspond to small changes in coordinates.  
//  The user can use the Method parameter to control the style of this mapping,
//  and the operation of the various internal engines that control the 
//  evolution of the ensemble.
//
//     Method = 0    is the simplest algorithm, mapping hypercube coordinates
//                   onto extended-integer position labels by simple raster,
//                   and using the "LifeStory1" diffusion engine to create,
//                   destroy, and move atoms.
//                   This is very likely to be enhanced by switching on
//                   various bits of the Method integer.
//
// if( Method & 1 ), hypercube coordinates are mapped to position labels along
//                   a space-filling Hilbert curve.
//                   This is generally recommended, but unnecessary if Ndim=1.
//
// if( Method & 2 ), the "LifeStory" diffusion engine will include interactions
//                   with atoms' neighbours, otherwise not.
//                   This is generally recommended for its extra power.
//
// if( Method & 4 ), the algorithm includes the Chameleon1 engine, which lets
//                   atoms jump from one ensemble object to another, without
//                   changing position.
//                   This is ineffective if there is only one ensemble object.
//
// if( Method & 8 ), the algorithm includes the Chameleon2 engine, which lets
//                   atoms from different ensemble objects exchange position. 
//                   This is ineffective if there is only one ensemble object.
//
// if( Method & 16), the algorithm will include the Leapfrog1 engine, which
//                   inverts atom positions with respect to another atom from
//                   the same or a different ensemble object.
//                   This is designed to assist when the posterior distribution
//                   in more than one dimension is highly non-spherical.
//
// if( Method & 32), the algorithm will include the Leapfrog2 engine, which
//                   reflects atom positions with respect to two other atoms
//                   from the same or different ensemble objects.
//                   This is designed to assist when the posterior distribution
//                   in more than one dimension is highly non-spherical.
//
// if( Method & 64), the algorithm will include the GuidedWalk engine, which
//                   moves atoms along directions parallel to displacements
//                   between neighbours, thus following the local shape of
//                   the posterior distribution.
//                   This may supersede the Leapfrog engines.
//
//  ________________________________________________________________________
// |                                                                        |
// | GuidedWalk Leapfrog2 Leapfrog1 Chameleon2 Chameleon1 LifeStory2 Hilbert| 1
// |    off        off       off       off       off      LifeStory1 raster | 0
// |________________________________________________________________________|
//       64        32        16         8         4          2         1 Method
//
//  In the interest of overall efficiency, the algorithm tries to equalise the
//  computation time between its engines so that even if all the optional
//  engines did nothing, the computation cost of including them would be
//  limited to the number of engines (currently a factor of 6).
//
//  Author generally recommends switching everything in by
//                   Method = 127 (equivalently -1) .
//  
//=============================================================================
//
// The "Massive Inference" (MassInf) option is provided for applications where
// atoms have intensities or "fluxes" about which the data are linear.
//
// The PRIOR for the number of atoms and for their location is as for BayeSys3.
//
// Prior on flux z of atom is Pr(z) = ProbON * P(z) + (1 - ProbON) * delta(z) ,
// i.e. each flux is expected to be distributed as P(z) with probability ProbON
//      otherwise it is switched off with zero value.
// Common->ProbON supplies ProbON, and the "units" decimal digit of the switch
// Common->MassInf defines the shape P(z) of the prior according to
//
// (0) "monkeys"              P(z) = delta(z-q),   i.e. z = q = constant
//                     
// (1) "positive"             P(z) = exp(-z/q) / q    in z > 0
//    
// (2) "positive/negative"    P(z) = exp(-|z|/q) / 2q
//                                         2     2               2
// (3) "Gaussian"             P(z) = exp(-z / 2 q ) / sqrt(2 pi q )
//
// When flux can be positive or negative, the "positive/negative" prior is
// more tolerant of dynamic range than is the Gaussian prior.
// In each case the flux unit q is fixed at the positive value supplied
// in Common->FluxUnit0, or (if <= 0) kept close to its most probable value.
//
// An atom can have Valency (=1,2,...) independent fluxes, provided
// (a) all atoms have the same valency,
// (b) the flux unit is the same for each,
// (c) the mock-data footprints of the different valencies of an atom
//     do not overlap (no common cells), and
// (d) (if using LifeStory2 engine) for any pair of atoms, each footprint
//     of one overlaps no more than one valency footprint of the other.
//
//
// To specify the LIKELIHOOD (used as its log), define
//
//                    f     = object = sum of atoms
//
//                    x[j]  = positional coordinate of atom j
//
//                    z[j]  = flux(es) of atom j
//
//                    Footprint(x) = mock data from unit flux at x
//
//                    Mock  = SUM[j] z[j] Footprint(x[j]) = mock data of object
//
// Likeliood Pr(Data | f) can be one of
//
// (0xx) "chisquared", derived from
//
//                    Data  = signal,
//
//                    Acc   = 1/sigma, can be 0,   (sigma = standard deviation)
//
//                    M     = # measured data, for which Acc[k] > 0.
//                               M                         2
//                    Z     = PRODUCT[k] sqrt(2 PI sigma[k] )
//                                                      2          2
//                    Chisq = SUM[k] (Mock[k] - Data[k]) / sigma[k]
//                                    -1   -Chisq / 2
//                    Pr(Data | f) = Z    e
//
// (1xx) "Poisson", for positive problems with "monkey" or "positive" prior;
//
//                    Data  = counts, can be floating-point but must be >= 0.0
//
//                    Acc   = extraneous "background", which may be small but
//                            must be strictly +ve wherever Data > 0.0
//                                               -F[k]     D[k]
//                    Pr(Data | f) = PRODUCT[k] e      F[k]    / D[k]!
//
//                    where  F = Mock + Acc, D = Data + Acc,  x! = GAMMA(x+1) :
//                    this form of likelihood avoids singularity when active
//                    cells happen to be empty of mock data, and reduces to
//                    standard Poisson when the background -> 0
//
// according to the "hundreds" digit of the switch Common->MassInf.
//
// For example, Common->MassInf =  1  is Gaussian data with "positive" prior;
//              Common->MassInf = 101  is Poisson data with "positive" prior.
//=============================================================================
// History:
// MassInf1 v1.12-2.36  1998-2000
// BayeSys1 code        1999-2000
// BayeSys2 v1.02-1.05  1 May - 10 Sep 2001
// BayeSys3 v1.01-1.32  2 Jan 2002 - 4 Feb 2003
//          v2.00      10 Feb 2003
//          v2.01      11 Jul 2003  Fixed-q0 option to fix FluxUnit0 in MassInf
//          v2.02      19 Jul 2003  MassInf options extended to allow flux=0
//          v2.10      12 Aug 2003  MassInf needs FluxUnit0 from user
//          v2.20      14 Aug 2003  Remove any overlay in MassInf user's "bits"
//          v3.00      20 Aug 2003  Poisson data as MassInf option
//          v3.01      12 Sep 2003  MassInf control of flux=0 by ProbON
//          v3.02      18 Sep 2003  BayeShape coordinate options
//          v3.03      11 Oct 2003  "Hilbert" name used instead of Peano
//          v3.04      20 Oct 2003  Display width of atom roughy halved
//          v3.05      27 Oct 2003  BayeShape permutation option
//          v3.10      29 Nov 2003  Free software under GNU LGPL
//          v3.11       3 Feb 2004  Control uses <#copies>.  Better Evid & Info
//-----------------------------------------------------------------------------
// 
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include "bayesys3.h"
#include "random.h"
#include "hilbert.h"

#undef  PARALLEL
#define PARALLEL 0   // # slaved parallel processors [0 = none]
#undef  DEBUG
#define DEBUG    0   // check Cube/Link consistency  [0 = off]

/***********************/
/* Internal structures */
/***********************/
typedef struct
{
    int            i;        // 1st object
    int            j;        // [2nd object]
    int            k;        // [3rd object]
    signed char    engine;   // operation type
    signed char    iType;    // type of 1st object, -1=WRITE, 0=Unused, 1=read
    signed char    jType;    // type of 2nd object, -1=WRITE, 0=Unused, 1=read