📄 bayesys3.c
字号:
signed char kType; // type of 3rd object, -1=WRITE, 0=Unused, 1=read
int iEarly; // min start time # for 1st object i
int jEarly; // min start time # for 2nd object j
int kEarly; // min start time # for 3rd object k
int iLate; // max start time # for 1st object i
int jLate; // max start time # for 2nd object j
int kLate; // max start time # for 3rd object k
int slave; // assigned processor
int CPU; // CPU time in this operation
int Success; // # successes in this operation
} OperStr; // Operation
/*********************************/
/* Internal constants and macros */
/*********************************/
//
// CALLOC(p,n,t) allocates vector p[0:n-1] of type t
// REALLOC(p,n,t) (re-)allocates vector p[0:n-1] of type t
// FREE(p) frees CALLOC or REALLOC or NULL vector p[0:*], sets p=NULL
// CALL catches negative error codes
// PLUS(x,y) log(exp(x)+exp(y))
// SWAP(x,y) exchange
//
#undef CALLOC
#define CALLOC(p,n,t) {p=NULL;\
if((n)>0&&!(p=(t*)calloc((size_t)(n),sizeof(t))))\
{CALLvalue=E_MALLOC;goto Exit;}/*printf("%p %d\n",p,(size_t)(n)*sizeof(t));*/}
#undef REALLOC
#define REALLOC(p,n,t) {/*printf("%p -1\n",p);*/\
if((n)>0&&!(p=(t*)realloc(p,(size_t)(n)*sizeof(t))))\
{CALLvalue=E_MALLOC;goto Exit;}/*printf("%p %d\n",p,(size_t)(n)*sizeof(t));*/}
#undef FREE
#define FREE(p) {if(p){/*printf("%p -1\n",p);*/(void)free((void*)p);} p=NULL;}
#undef CALL
#define CALL(x) {if( (CALLvalue = (x)) < 0 ) goto Exit;}
#undef PLUS
#define PLUS(x,y) ((x)>(y)?(x)+log(1.+exp((y)-(x))):(y)+log(1.+exp((x)-(y))))
#undef SWAP
#define SWAP(a,b) {swap = a; a = b; b = swap;}
/***********************/
/* Internal prototypes */
/***********************/
// BayeSys3 main procedures
static void ShapeSolve (int, double, double, double*, double*);
static double ShapeCumul (int, double, double*);
static void ShapeIndex (int, double*, double*);
static int BayesAlloc (CommonStr*, ObjectStr*, Node*, Node*);
static int BayesInit (CommonStr*, ObjectStr*);
static int PriorInit (CommonStr*, ObjectStr*, Node*, int, double*);
static int Control (CommonStr*, ObjectStr*, int, double*, double*);
static double Copies (double*, int);
static int Anneal (CommonStr*, ObjectStr*, int, double);
static void Sort2 (int*, double*, int);
static int FillOcean (CommonStr*, ObjectStr*, Node*, int);
static int MCMCengines (CommonStr*, ObjectStr*, Node*);
static void BayesFree (CommonStr*, ObjectStr*, Node*, Node*);
// Engines
static int DoOperations(OperStr*, CommonStr*, ObjectStr*, Node*, int);
static int Do1operation(OperStr*, CommonStr*, ObjectStr*, Node*, int);
static int SetPrior (OperStr*, CommonStr*, ObjectStr*, Node*);
static int CopyObject (OperStr*, CommonStr*, ObjectStr*);
static int MCMCsetup (OperStr*, CommonStr*, ObjectStr*, Node*);
static int LifeStory1 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int LifeStory2 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int Chameleon1 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int Chameleon2 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int Leapfrog1 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int Leapfrog2 (OperStr*, CommonStr*, ObjectStr*, Node*);
static int GuidedWalk (OperStr*, CommonStr*, ObjectStr*, Node*);
static int MCMCexit (OperStr*, CommonStr*, ObjectStr*, Node*);
// Event library
static double Birth (CommonStr*, int);
static double Death (CommonStr*, int);
static int SafeCube (CommonStr*, ObjectStr*);
// Atom control
static int Empty (CommonStr*, ObjectStr*);
static int Try1 (double*, CommonStr*, ObjectStr*);
static int Try2 (double*, double*, CommonStr*, ObjectStr*);
static int Insert1 (CommonStr*, ObjectStr*, Node*);
static int Insert2 (CommonStr*, ObjectStr*, Node*);
static int Delete1 (int,double*,double*, CommonStr*, ObjectStr*, Node*);
// Label library
static void Topology (CommonStr*);
static void CubetoLabel (unsigned*, ObjectStr*, CommonStr*);
static void LabeltoCube (ObjectStr*, unsigned*, CommonStr*);
static void RanLabel (unsigned*, int, unsigned*);
static void CopyLabel (unsigned*, unsigned*, int);
static void NewLabel (unsigned*, unsigned*, unsigned*, int,
int, int, unsigned*);
static int Outside (unsigned*, unsigned*, unsigned*, int);
// MassInf ancillary library
static int FluxAlloc (CommonStr*, ObjectStr*);
static int FluxInit (CommonStr*, ObjectStr*);
static int FluxCalib0 (CommonStr*, ObjectStr*, Node*);
static void FluxCalib (CommonStr*, ObjectStr*);
static void FluxFree (CommonStr*, ObjectStr*);
// MassInf outer procedures
static void FluxEmpty (double*, CommonStr*, ObjectStr*);
static int FluxTry1 (double*, CommonStr*, ObjectStr*);
static int FluxTry2 (double*, double*, CommonStr*, ObjectStr*);
static int FluxInsert1 (double*, CommonStr*, ObjectStr*);
static int FluxInsert2 (double*, CommonStr*, ObjectStr*);
static int FluxDelete1 (double*, CommonStr*, ObjectStr*);
// MassInf application procedures for user Footprints
static int AtomBits (double*, CommonStr*, int*, double*, int*, int*);
static int Overlap1 (int*, int*, int*, int);
static int Overlap2 (int*, int*, int*, int*, int*, int*, int);
static void SetIndex (int, int*, int*, int*, int*, int*, int*);
static void InsBits (double*, double*, int*, int*, double*, int);
static void DelBits (double*, double*, int*, int*, double*, int);
// MassInf application procedures for quadratic chisquared
static void SetGrad1 (double*, double*, double*, int, int*, int*, double*,
double*, double*);
static void SetGrad2 (double*, double*, double*, double*, int, int*, int*,
double*, int*, int*, double*, double*, double*,
double*, double*, double*, int*);
// MassInf application procedures for quadratic probabilities
static double GaussLhood (int, double*, double*, double*);
static double GaussTry1 (int, double, double, double, int, double*, double*);
static void GaussTry2 (int, double, double, double, int, double*, double*,
double*, double*, double*, int*, double*, double*);
static void GaussInsert1(int, Rand_t, double, double, double, int, double*,
double*, double*);
static void GaussInsert2(int, Rand_t, double, double, double, int, double*,
double*, double*, double*, double*, int*, double*,
double*);
static double GaussLhood1 (int, double*, double*, double*, int);
static double GaussLhood2 (int, double*, double*, double*, double*, double*,
int*, double*, double*);
// without Valency ....
static double Gauss1Marginal(int, double, double, double, double, double);
static double Gauss2Marginal(int, double, double, double, double, double,
double, double, double);
static double Gauss1Sample (int, Rand_t, double, double, double, double,
double);
static void Gauss2Sample (int, Rand_t, double, double, double, double,
double, double, double, double, double*, double*);
// without ProbON ....
static double gauss1marginal(int, double, double, double, double);
static double gauss2marginal(int, double, double, double, double, double,
double, double);
static double gauss1sample (int, Rand_t, double, double, double, double);
static void gauss2sample (int, Rand_t, double, double, double, double,
double, double, double, double*, double*);
// MassInf application procedures for Poisson probabilities
static double PoissLhood (int, double*, double*, double*);
static int PoissTry1 (int, double, double, double, double*, double*,
double*, int*, int, int*, int*, double*, double*);
static int PoissTry2 (int, double, double, double, double*, double*,
double*, int*, double*, int, int*, int*, double*,
int*, int*, double*, int*, double*, double*);
static int PoissInsert1 (int, Rand_t, double, double, double, double*,
double*, double*, int*, int, int*, int*, double*,
double*);
static int PoissInsert2 (int, Rand_t, double, double, double, double*,
double*, double*, int*, double*, int, int*, int*,
double*, int*, int*, double*, int*, double*,
double*);
static double PoissLhood1 (double*, double*, double*, int, int*, int*,
double*, double*, int);
static double PoissLhood2 (double*, double*, double*, double*, int, int*,
int*, double*, int*, int*, double*, int*,
double*, double*);
// without Valency ....
static int Poiss1Marginal(int, double, double, double, double*, double*,
double*, int*, int, int*, double*, double*);
static int Poiss2Marginal(int, double, double, double, double*, double*,
double*, int*, double*, int, int*, double*, int,
int*, double*, double*, double*);
static int Poiss1Sample (int, Rand_t, double, double, double, double*,
double*, double*, int*, int, int*, double*,
double*);
static int Poiss2Sample (int, Rand_t, double, double, double, double*,
double*, double*, int*, double*, int, int*,
double*, int, int*, double*, double*, double*);
// without ProbON ....
static int Poisson1 (int, Rand_t, double, double, double*, double*,
double*, int*, int, int*, double*, double*);
static int Poisson2 (int, Rand_t, double, double, double*, double*,
double*, int*, double*, int, int*, double*, int,
int*, double*, double*, double*);
static double poiss1lhood (double, double*, double*, double*, int, int*,
double*);
static void poiss2lhood (double, double, double*, double*, double*,
double*, int, int*, double*, int, int*, double*,
double*, double*);
//=============================================================================
//
// BayeSys3 main procedures
// main
// |
//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// BayeSys3
// | \_____________________________________________________________
// | | | | | | |
// | PriorInit Control MCMCengines BayesAlloc BayesFree FillOcean
// | \ Anneal / \ BayesInit
// | \ | / \ /
// | DoOperations Topology
// | |
// | Do1operation
// | _____________|______________________________________________
// | | | | | | |
// | CopyObject SetPrior MCMCsetup LifeStory1 LifeStory2 MCMCexit
// | | | Chameleon1 /|
// | | | Chameleon2 / |
// | | | Leapfrog1 / |
// | | | Leapfrog2 / |
// | | | GuidedWalk / |
// | | | \ / |
// | Empty Empty Insert1 Insert2
// | Insert1 Insert1 Delete1 Try2
// | Try1 Try1
// |
// | Empty Insert1 Try1 Delete1 Insert2 Try2
// | | | \_________|_\______|_\_________|_\_________|_\___
// | | | | | | | |
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
// | | | | | | | |
// | UserEmpty UserInsert1 UserTry1 UserDelete1 UserInsert2 UserTry2 |
// | |
// | |
// UserMonitor UserFoot
//=============================================================================
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: BayeSys3
//
// Purpose: Perform Markov chain Monte Carlo algorithm for
// Bayesian analysis or maximisation over unit hypercube,
// using an ensemble of one or more objects.
//
// History: JS BayeSys1 3 Mar 1999 - 3 Nov 2000
// BayeSys2 25 Apr 2001 - 10 Sep 2001
// BayeSys3 12 Jan 2002 - 3 Feb 2004
//-----------------------------------------------------------------------------
//
int BayeSys3( // O +ve finish code from UserMonitor (or error)
CommonStr* Common, // I O General information
ObjectStr* Objects) // O ENSEMBLE of objects
{
static const double BIG = DBL_MAX / 3.0000; // for annealing protection
static const int NOCEAN = 50000; // max # atoms in Ocean, for atom-widths
int NTABLE = 1000; // current size of evidence table
#undef NEXTRA
#define NEXTRA 12 /* # likelihoods in prior settings */
int ENSEMBLE = Common->ENSEMBLE;
int Valency = Common->Valency;
double Lextra[NEXTRA]; // Old likelihoods for stability
Node* Links = NULL; // Linked lists of atom labels
Node Ocean[1]; // Linked list of accumulated atoms
double cold; // Old reciprocal temperature
double dcool; // Cooling increment
double Lbar; // <logL>
double* Ltable = NULL; // table of strictly increasing logL values
double* ctable = NULL; // table of strictly increasing cool values
int* ntable = NULL; // table of occupation numbers for averaging
double Evid; // SUM[1..] Ltable[i] * (ctable[i]-ctable[i-1])
int j; // table counter
int k; // object counter
int finish; // exit code from UserMonitor
double t; // temporary
int CALLvalue = 0;
#if DEBUG
printf("WARNING: DEBUG switched on\n");
#endif
// Check input parameters
if( Common->ENSEMBLE < 1 ) return E_BAYESYS_PARMS;
if( Common->MinAtoms < 1 ) return E_BAYESYS_PARMS;
if( Common->MaxAtoms && Common->MaxAtoms < Common->MinAtoms )
return E_BAYESYS_PARMS;
if( Common->MaxAtoms < Common->MinAtoms && Common->Alpha == 0.0 )
return E_BAYESYS_PARMS;
if( Common->Rate <= 0.0 ) return E_BAYESYS_PARMS;
if( Common->Ndim < 1 ) return E_BAYESYS_PARMS;
if( Valency )
{
if( Common->MassInf != 0 && Common->MassInf != 1
&& Common->MassInf != 2 && Common->MassInf != 3
&& Common->MassInf != 100 && Common->MassInf != 101 )
return E_MASSINF_PARMS;
if( Common->Ndata < 0 ) return E_MASSINF_PARMS;
if( ! Common->Data ) return E_MASSINF_PARMS;
if( ! Common->Acc ) return E_MASSINF_PARMS;
if( Common->ProbON <= 0.0 ) return E_MASSINF_PARMS;
if( Common->MassInf >= 100 )
for( k = 0; k < Common->Ndata; k++ )
{
if( Common->Acc[k] < 0. || Common->Data[k]+Common->Acc[k] < 0. )
return E_MASSINF_DATA;
if( Common->Acc[k] == 0.0 && Common->Data[k] != 0.0 )
return E_MASSINF_DATA;
}
}
// Allocate memory
CALLOC(Ltable, NTABLE, double)
CALLOC(ctable, NTABLE, double)
CALLOC(ntable, NTABLE, int)
CALLOC(Links, ENSEMBLE, Node)
CALL( BayesAlloc(Common, Objects, Links, Ocean) )
// Initialise system
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -