bayesys3.c

贝叶斯算法：盲分离技术

        Object = &Objects[k];
        j = Ranint(Rand);
        if( j < 0 )
            j = ~j;
        CALL( RanInit(Object->Rand, j) )
        Objects[k].Natoms = 0;
        for( j = 0; j < Object->Nstore; j++ )
            for( i = 0; i < Nsize; i++ )
                Object->Cubes[j][i] = 0.0;
        Object->reset = 1;
    }
    Topology(Common);
    if( Valency )
        CALL( FluxInit(Common, Objects) )
Exit:
    return CALLvalue;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  PriorInit
// 
// Purpose:   Initialise prior objects and extra loglikelihood values
// 
// History:   JS        16 Oct 2002, 1 Feb 2003
//-----------------------------------------------------------------------------
static
int PriorInit(       //   O  0, or -ve error
CommonStr* Common,   // I O  general information
ObjectStr* Objects,  //   O  ENSEMBLE of objects
Node*      Links,    //  (O) empty on exit
int        NEXTRA,   // I    # extra loglikelihood values
double*    Lextra)   //   O  extra loglikelihood values   [NEXTRA]
{
    int       ENSEMBLE   = Common->ENSEMBLE;
    OperStr*  Operations = NULL; // list of operations   [ENSEMBLE]
    OperStr*  Oper;              // & individual operation
    int       i, j, k;
    int       CALLvalue  = 0;

// Pre-calibrate objects
    if( Common->Valency )
        CALL( FluxCalib0(Common, Objects, Links) )
    CALLOC(Operations, ENSEMBLE, OperStr)
// Find O(10) preliminary prior objects, preferably all with different L
    Operations->engine = 0;       // SetPrior
    Operations->iType  = -1;    Operations->i = 0;
    Operations->jType  = 0;
    Operations->kType  = 0;
    for( k = 0; k < NEXTRA; k++ )
    {
        for( i = 0; i < NEXTRA; i++ )
        {                     // try to ensure all L different but don't insist
            CALL( DoOperations(Operations, Common, Objects, Links, 1) )
            for( j = 0; j < k; j++ )
                if( Objects->Lhood == Lextra[j] )
                    break;
            if( j == k )
                break;
        }
        Lextra[k] = Objects->Lhood;
    }
// Set ENSEMBLE of prior objects
    for( i = 0; i < ENSEMBLE; i++ )
    {
        Oper = &Operations[i];
        Oper->engine =  0;               // SetPrior
        Oper->iType  = -1;    Oper->i = i;
        Oper->jType  =  0;
        Oper->kType  =  0;
    }
    CALL( DoOperations(Operations, Common, Objects, Links, ENSEMBLE) )
Exit:
    FREE(Operations)
    return CALLvalue;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Control
//
// Purpose:   Rate-limited allowable cooling.
//            Aim for about Rate copy operations per object
//            (see Anneal for documentation on this).
//            
//            Guard against accidental near-coalescence of likelihood values
//            (which would allow arbitrarily large cooling), by including
//            artificial extra loglikelihood values, perhaps from old ensembles
//
// History:   John Skilling   28 Jan 2002, 25 Mar 2002, 18 Jul 2002
//                            17 Sep 2002  Rate pertains to most likely object
//                            17 Dec 2002  Lextra kept different
//                            20 Aug 2003  Separate Sort2 call
//                             3 Feb 2004  Use <# copies> not #copies(Lmax)
//-----------------------------------------------------------------------------
static 
int Control(        //   O  0, or -ve error
CommonStr* Common,  // I    General information
ObjectStr* Objects, // I    Sample objects, with Lhood               [ENSEMBLE]
int        Nextra,  // I    # extra objects
double*    Lextra,  // I O  Extra objects for stability                [Nextra]
double*    dcool)   //   O  cooling increment allowed by Rate
{
    double    Rate     = Common->Rate;
    int       ENSEMBLE = Common->ENSEMBLE;
    unsigned* Rand     = Common->Rand;
    double    Lmid;           // central likelihood value
    double    R;              // Rate, adjusted for finite allowable # copies
    int       N;              // total # likelihoods
    int       i, j;           // counters
    int       chop;           // binary chop counter
    double    a, copya;       // binary chop x and y (left)
    double    b, copyb;       // binary chop x and y (mid)
    double    c, copyc;       // binary chop x and y (right)
    double*   wa;             // binary chop data vector (left)
    double*   wb;             // binary chop data vector (mid)
    double*   wc;             // binary chop data vector (right)
    double*   work1  = NULL;  // workspace
    double*   work2  = NULL;  // workspace
    double*   work3  = NULL;  // workspace
    double*   Lvalue = NULL;  // workspace for sorted likelihoods
    double*   swap;           // exchange pointer
    double    Lmax;           // largest logL
    int       CALLvalue = 0;

// Default if all L equal or other anomaly
    *dcool = 0.0;

// Read all likelihoods
    N = ENSEMBLE + Nextra;
    CALLOC(Lvalue, N, double)
    CALLOC(work1, N, double)
    CALLOC(work2, N, double)
    CALLOC(work3, N, double)
    for( j = 0; j < ENSEMBLE; j++ )
        Lvalue[j] = Objects[j].Lhood;
    for( ; j < N; j++ )
        Lvalue[j] = Lextra[j-ENSEMBLE];
    if( N <= 1 )
        return 0;

// Offset loglikelihoods to MAX=0 for safety
    Lmax = Lvalue[0];
    for( i = 1; i < N; i++ )
        if( Lmax < Lvalue[i] )
            Lmax = Lvalue[i];
    for( i = 0; i < N; i++ )
        Lvalue[i] -= Lmax;
// Max # possible copies = max # recipients (per sample)
    R = N;
    for( i = 0; i < N; i++ )
        if( Lvalue[i] == 0.0 )
            R -= 1.0;
    R /= N;
    if( R > 0.0 )
    {
// Ensure # copies < max possible (the 3. is for rough backward compatibility)
        R = 1.0 / (3. / Rate + 1.0 / R);
// Guess initial value for dcool from normally distributed Lvalues
        a = b = 0.0;
        for( i = 0; i < N; i++ )
            a += Lvalue[i];
        a /= N;
        for( i = 0; i < N; i++ )
        {
            c = Lvalue[i] - a;
            b += c * c;
        }
        b = sqrt(b / N);
        a = b = 2.5066 * R / b;
// Bracket dcool by interval [a,b] a factor of 2 wide
        wa = work1;   wb = work2;
        for( i = 0; i < N; i++ )
            wa[i] = wb[i] = exp(a * Lvalue[i]);
        copya = copyb = Copies(wa, N);
        if( copya < R )
        {
            do
            {
                a = b;    copya = copyb;
                swap = wa;    wa = wb;    wb = swap;
                b = 2.0 * a;    
                for( i = 0; i < N; i++ )
                    wb[i] = wa[i] * wa[i];
                copyb = Copies(wb, N);
            }
            while( copyb < R );
        }
        else
        {
            do
            {
                b = a;    copyb = copya;
                swap = wb;    wb = wa;    wa = swap;
                a = b / 2.0;
                for( i = 0; i < N; i++ )
                    wa[i] = sqrt(wb[i]);
                copya = Copies(wa, N);
            }
            while( copya > R );
        }
// Binary chop to accuracy 1 in 1000
        wc = work3;
        for( chop = 0; chop < 10; chop++ )
        {
            c = (a + b) / 2;
            for( i = 0; i < N; i++ )
                wc[i] = sqrt(wa[i] * wb[i]);
            copyc = Copies(wc, N);
            if( copyc < R )
            {
                a = c;    copya = copyc;
                swap = wa;    wa = wc;    wc = swap;
            }
            else
            {
                b = c;    copyb = copyc;
                swap = wb;    wb = wc;    wc = swap;
            }
        }
// Final linear interpolation
        *dcool = (b * (R - copya) + a * (copyb - R)) / (copyb - copya);
    }

    if( Nextra )
    {
// Keep Lextra up-to-date enough, trying to keep all different
        for( j = 0; j < Nextra; j++ )
        {
            Lmid = Objects[Rangrid(Rand, ENSEMBLE)].Lhood;
            for( i = 0; i < Nextra; i++ )
                if( Lextra[i] == Lmid )
                    break;
            if( i == Nextra )
            {
                Lextra[Rangrid(Rand, Nextra)] = Lmid;
                break;
            }
        }
    }

// Phenomenological correction for finite ENSEMBLE, reducing effective dcool
// (Undo correction at top of Anneal)
    *dcool *= (N - 1.0) / N;
Exit:
    FREE(work3);
    FREE(work2);
    FREE(work1);
    FREE(Lvalue);
    return  CALLvalue;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Copies
// 
// Purpose:   <# copy operations per object> for weights w
//
//              = SUM | w - <w> |  /  2 SUM w
//
// History:   JS     3 Feb 2004
//-----------------------------------------------------------------------------
static
double Copies(
double*  w,      // I   weights
int      N)      // I   dimension
{
    double sumw, wbar, copy;
    int    i;

    sumw = 0.0;
    for( i = 0; i < N; i++ )
        sumw += w[i];
    wbar = sumw / N;
    copy = 0.0;
    for( i = 0; i < N; i++ )
        copy += fabs(w[i] - wbar);
    return  copy / (2.0 * sumw);
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Anneal
//
// Purpose:   Anneal ensemble by dcool, by copying from source objects "src"
//            to destination objects "dest", as uniformly as possible
//            consistently with excess/deficit real multiplicities N:-
//
//               N(src)  = exp(dcool * (Lhood[src] - Lbar)) - 1.0
//                                                   for Lhood[src] > Lbar only
//               N(dest) = 1.0 - exp(dcool * (Lhood[dest] - Lbar))
//                                                  for Lhood[dest] < Lbar only
//            where Lbar is adjusted for equality in
//               # copy operations = SUM[src] N(src) = SUM[dest] N(dest).
//
//            A source will be copied either n or n+1 times, where
//               n = (int)N(src) = 0,1,2,...
//            A destination will be overwritten, with probability N(dest).
//
// History:   John Skilling   28 Jan 2002, 25 Mar 2002, 18 Jul 2002
//                            17 Sep 2002  Return # copies
//                             3 Feb 2004  Internal vector allocations
//-----------------------------------------------------------------------------
static
int Anneal(         //   O  #copies diagnostic, or -ve error
CommonStr* Common,  // I    CopyObject info
ObjectStr* Objects, // I    Sample objects                           [ENSEMBLE]
int        Nextra,  // I    # extra objects, for phenomenological de-correction
double     dcool)   // I    Cooling increment
{
    int       ENSEMBLE   = Common->ENSEMBLE;
    unsigned* Rand       = Common->Rand;
    OperStr*  Operations = NULL;             // list of operations   [ENSEMBLE]
    OperStr*  Oper;           // & individual operation
    double    Lmid;           //