📄 bayesys3.c
字号:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: ShapeSolve
//
// Purpose: Latitude phi for which
// ShapeCumul(phi) / ShapeCumul(PI/2) = outward
// Equator is outward = 0, phi = 0
// North pole is outward = 1, phi = PI/2
//
// Method: FROWN(a) = ShapeCumul(a) is concave,
// FROWN''(a) < 0,
// so if a < aim, tangent meets aim at better lower bound
// SMILE(b) = - (ShapeCumul(PI/2) - ShapeCumul(b))^(1/N) is convex,
// SMILE''(b) > 0,
// so if b > aim, tangent meets aim at better upper bound
// Chop between lower and upper bounds
//
// History: JS 18 Sep 2003
//-----------------------------------------------------------------------------
static
void ShapeSolve(
int N, // I dimension of surface of sphere
double outward, // I poleward cumulant
double inwardN, // I (equator-ward cumulant)^(1/N)
double* cosp, // O cos(latitude)
double* sinp) // O sin(latitude)
{
static const double pibytwo = 1.57079632679489661923;
double frac = 1.0 / N;
double a; // left bound
double ya; // FROWN(a)
double ma; // FROWN'(a) > 0
double b; // right bound
double yb; // SMILE(b)
double mb; // SMILE'(b) > 0
double c; // central estimate
double yc; // ShapeCumul(c)
double mc; // ShapeCumul'(c)
double top; // ShapeCumul(PI/2)
double aim; // required value for FROWN
double bim; // required value for SMILE
int iter; // 5 chops suffice for full accuracy
if( N == 1 )
{
a = outward * pibytwo;
*cosp = cos(a);
*sinp = sin(a);
return;
}
if( N == 2 )
{
*cosp = sqrt(1.0 - outward * outward);
*sinp = outward;
return;
}
top = ShapeCumul(N, pibytwo, &mc);
a = outward * top;
b = a * a * (N - 2) / 6.0;
if( b * b <= DBL_EPSILON )
{
*sinp = a * (1.0 + b);
*cosp = sqrt(1.0 - *sinp * *sinp);
return;
}
a = inwardN * pow(N * top, frac);
b = a * a / (2.0 * N + 4.0);
if( b * b <= DBL_EPSILON )
{
*cosp = a * (1.0 - b);
*sinp = sqrt(1.0 - *cosp * *cosp);
return;
}
aim = outward * top;
bim = - pow(top - aim, frac);
a = 0.0; ya = 0.0; ma = 1.0;
b = pibytwo; yb = 0.0; mb = pow(frac, frac);
for( iter = 0; iter < 5; iter++ )
{
c = (a + b) / 2.0;
yc = ShapeCumul(N, c, &mc);
if( yc < aim )
{
a = c;
ya = yc; ma = mc;
a += (aim - ya) / ma;
}
else
{
b = c;
yb = - pow(top - yc, frac); mb = - mc * frac * yb / (top - yc);
b += (bim - yb) / mb;
}
}
// protected interpolation
ya = ShapeCumul(N, a, &ma) - aim;
yb = ShapeCumul(N, b, &mb) - aim;
if( ya >= yb )
c = (a + b) / 2.0; // unusual
else
{
c = (a * yb - b * ya) / (yb - ya); // linear
mc = (mb - ma) / (yb - ya);
c += 0.5 * mc * (b - c) * (c - a); // quadratic correction
}
if( c < a ) c = a; // final protection
if( c > b ) c = b; // final protection
*cosp = cos(c);
*sinp = sin(c);
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: ShapeCumul
//
// Purpose: phi N-1
// INTEGRAL cos (theta) d(theta)
// 0
//
// History: JS 18 Sep 2003
//-----------------------------------------------------------------------------
static
double ShapeCumul( // O integral
int N, // I dimension of surface of hemisphere
double phi, // I latitude
double* deriv) // O cos(phi)^(N-1)
{
double r, s, t, u, sinp, cosp, cos2;
int i;
sinp = sin(phi);
cosp = cos(phi);
cos2 = cosp * cosp;
s = t = u = 1.0;
for( i = (N & 1) + 2; i < N; i += 2 )
{
r = 1.0 / i;
t *= 1.0 - r;
u *= cos2;
s += t * u;
}
t *= N - 1.0;
s *= sinp;
u *= cosp;
if( N & 1 )
{
s = (phi + cosp * s) / 2.0;
u *= cosp;
t /= 2.0;
}
*deriv = u;
return s / t;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: ShapeIndex
//
// Purpose: Index vector y into increasing order
// y[p[0]] <= y[p[1]] <= .... <= y[p[n-1]]
//
// History: JS 27 Apr 1999, 27 Oct 2003
//-----------------------------------------------------------------------------
static
void ShapeIndex(
int N, // I dimension
double* y, // I vector being indexed by value
double* p) // O index (integer values)
{
int j, k, l, m, q;
double r;
if( N <= 1 )
{
p[0] = 0.0;
return;
}
for( j = 0; j < N; j++ )
p[j] = (double)j;
l = N / 2;
k = N - 1;
while( 1 )
{
if( l > 0 )
{
q = (int)p[--l];
r = y[q];
}
else
{
q = (int)p[k];
r = y[q];
p[k--] = p[0];
if( k == 0 )
{
p[0] = (double)q;
return;
}
}
m = l;
j = l + l + 1;
while( j <= k )
{
if( j < k && y[(int)p[j]] < y[(int)p[j+1]] )
j++;
if( r < y[(int)p[j]] )
{
p[m] = p[j];
m = j;
j = j + j + 1;
}
else
j = k + 1;
}
p[m] = (double)q;
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: BayesAlloc
//
// Purpose: Allocate memory for BayeSys (and MassInf)
//
// History: JS 16 Oct 2002, 1 Feb 2003, 13 Sep 2003
//-----------------------------------------------------------------------------
static
int BayesAlloc( // O 0, or -ve error
CommonStr* Common, // I(O) general information
ObjectStr* Objects, // (O) new ENSEMBLE of objects [ENSEMBLE]
Node* Links, // (O) new trees [ENSEMBLE]
Node* Ocean) // (O) new ocean [1]
{
int ENSEMBLE = Common->ENSEMBLE;
int Ndim = Common->Ndim;
int Valency = Common->Valency;
int Nsize = Ndim + Valency + 1;
ObjectStr* Object;
int j, k;
int CALLvalue = 0;
// Allow subsequent CALLOC memory allocation to fail gracefully
Common->offset = NULL;
Common->permute = NULL;
// Allocate
// Common
CALLOC(Common->offset, Ndim, unsigned)
CALLOC(Common->permute, Ndim, int)
// Objects
for( k = 0; k < ENSEMBLE; k++ )
{
Object = &Objects[k];
Object->xLabel = NULL;
Object->yLabel = NULL;
Object->xOrigin = NULL;
Object->yOrigin = NULL;
Object->xTry = NULL;
Object->yTry = NULL;
Object->work = NULL;
Object->Cubes = NULL;
CALLOC(Object->xLabel, Ndim, unsigned)
CALLOC(Object->yLabel, Ndim, unsigned)
CALLOC(Object->xOrigin, Ndim, unsigned)
CALLOC(Object->yOrigin, Ndim, unsigned)
CALLOC(Object->xTry, Ndim, unsigned)
CALLOC(Object->yTry, Ndim, unsigned)
CALLOC(Object->work, Ndim, unsigned)
Object->Nstore = Common->MinAtoms + 2;
CALLOC(Object->Cubes, Object->Nstore, double*)
for( j = 0; j < Object->Nstore; j++ )
CALLOC(Object->Cubes[j], Nsize, double)
CALL( SetLink(Ndim, &Links[k]) )
}
if( Valency )
CALL( FluxAlloc(Common, Objects) )
// Ocean
CALL( SetLink(Ndim, Ocean) )
Exit:
return CALLvalue;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: BayesInit
//
// Purpose: Initialise BayeSys (and MassInf)
//
// History: JS 16 Oct 2002, 1 Feb 2003
//-----------------------------------------------------------------------------
static
int BayesInit( // O 0, or -ve error
CommonStr* Common, // I O general information
ObjectStr* Objects) // O new ENSEMBLE of objects
{
unsigned* Rand = Common->Rand;
int ENSEMBLE = Common->ENSEMBLE;
int Ndim = Common->Ndim;
int Valency = Common->Valency;
int Nsize = Ndim + Valency + 1;
ObjectStr* Object;
int i, j, k;
int CALLvalue = 0;
// Initialise BayeSys
Common->Nbits = 0;
for( Rand[0] = 1; Rand[0]; Rand[0] <<= 1 )
Common->Nbits++;
CALL( RanInit(Rand, Common->Iseed) )
Common->Iseed = CALLvalue;
Common->cool = 0.0;
Common->Evidence = 0.0;
Common->Information = 0.0;
Common->Nsystem = 0;
Common->Success = 0.0;
Common->CPU = 0.0;
for( k = 0; k < ENSEMBLE; k++ )
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -