C++ InfAlg类代码示例

/** Only maximization step included
 *
 *  This function performs maximization based on the evidence collected from tab file
 *  Since there's no belief propagation stage, it does not return proper likelihood.
 */
Real UnsafeEmalg::iterateWithoutEstep(MaximizationStep &mstep) {
    Real logZ = 0;
    Real likelihood = 0;

    // Expectation calculation
    int nProcessed = 0;

    for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
        InfAlg* clamped = _estep.clone();

        // Apply evidence
        for( Evidence::Observation::const_iterator i = e->begin(); i != e->end(); ++i )
            clamped->clamp( clamped->fg().findVar(i->first), i->second );

        mstep.addExpectations( *clamped );

        delete clamped;

        nProcessed++;

        if (nProcessed % 10 == 0) {
            cout << nProcessed << " number of samples processed" << endl;
        }
    }

    // Maximization of parameters
    mstep.maximize(_estep.fg());

    return likelihood;
}

string mapper(EMdata &dat) {
	FactorGraph fg = dat.fg;

	PropertySet infprops = getProps();

	InfAlg* inf = newInfAlg(INF_TYPE, fg, infprops);
	inf->init();

	// Read sample from file
	Evidence e;
	istringstream estream(dat.tabFile);
	e.addEvidenceTabFile(estream, fg);

	// Read EM specification
	istringstream emstream(dat.emFile);
	EMAlg em(e, *inf, emstream);

	// perform 1 iteration e-step
	Real likelihood = hadoop_iterate(em._msteps, em._evidence, em._estep);

	dat.likelihood = likelihood;
	dat.msteps = em._msteps;

	// Clean up
	delete inf;

	dat.emFile = ""; // reduce amount of serialization
	dat.tabFile = "";

	return emToString(dat);
}

int main() {
	// This is the INFERENCE/EM engine, derived from the 
    // libDAI example program (example_sprinkler_em)
    // (http://www.cs.ubc.ca/~murphyk/Bayes/bnintro.html)
    //
    // The factor graph file (input.fg) has to be generated first
	// and the data file (sprinkler.tab),
	// as well as the EM commands 

    // Read the factorgraph from the file
    FactorGraph Network;
    Network.ReadFromFile( "input.fg" );

    // Prepare junction-tree object for doing exact inference for E-step
    PropertySet infprops;
    infprops.set( "verbose", (size_t)1 );
    infprops.set( "updates", string("HUGIN") );
    infprops.set( "maxiter", string("1000") );
    infprops.set( "tol", string("0.00001") );
    infprops.set( "logdomain", true);
    infprops.set( "updates", string("SEQFIX") );
    InfAlg* inf = newInfAlg("BP", Network, infprops );
    inf->init();

    // Read sample from file
    Evidence e;
    ifstream estream( "input.tab" );
    e.addEvidenceTabFile( estream, Network );
    cerr << "Number of samples: " << e.nrSamples() << endl;

    // Read EM specification
    ifstream emstream( "input.em" );
    EMAlg em(e, *inf, emstream);

    // Iterate EM until convergence
    while( !em.hasSatisfiedTermConditions() ) {
        Real l = em.iterate();
        cerr << "Iteration " << em.Iterations() << " likelihood: " << l <<endl;
		Real c = inf->logZ();
        cerr << "Iteration infAlg " << em.Iterations() << " likelihood: " << c <<endl;
    }

    cout.precision(6);
    cout << inf->fg();

    delete inf;

    return 0;
}

EMAlg::EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &msteps_file )
  : _evidence(evidence), _estep(estep), _msteps(), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
{
    msteps_file.exceptions( std::istream::eofbit | std::istream::failbit | std::istream::badbit );
    size_t num_msteps = -1;
    msteps_file >> num_msteps;
    _msteps.reserve(num_msteps);
    for( size_t i = 0; i < num_msteps; ++i )
        _msteps.push_back( MaximizationStep( msteps_file, estep.fg() ) );
}

void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
    for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
        Permute &perm = i->second;
        VarSet &vs = _varsets[i->first];

        Factor b = alg.belief(vs);
        Prob p( b.states(), 0.0 );
        for( size_t entry = 0; entry < b.states(); ++entry )
            p[entry] = b[perm.convertLinearIndex(entry)]; // apply inverse permutation
        _estimation->addSufficientStatistics( p );
    }
}

void SharedParameters::collectExpectations( InfAlg &alg ) {
    for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
        Permute &perm = i->second;
        VarSet &vs = _varsets[i->first];

        Factor b = alg.belief(vs);
        Prob p( b.nrStates(), 0.0 );
        for( size_t entry = 0; entry < b.nrStates(); ++entry )
            p.set( entry, b[perm.convertLinearIndex(entry)] ); // apply inverse permutation
        (*_expectations) += p;
    }
}

Real EMAlg::iterate( MaximizationStep &mstep ) {
    Real logZ = 0;
    Real likelihood = 0;

    _estep.run();
    logZ = _estep.logZ();

    // Expectation calculation
    for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
        InfAlg* clamped = _estep.clone();
        // Apply evidence
        for( Evidence::Observation::const_iterator i = e->begin(); i != e->end(); ++i )
            clamped->clamp( clamped->fg().findVar(i->first), i->second );
        clamped->init();
        clamped->run();

        likelihood += clamped->logZ() - logZ;

        mstep.addExpectations( *clamped );

        delete clamped;
    }

    // Maximization of parameters
    mstep.maximize( _estep.fg() );

    return likelihood;
}

Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
    Factor Bi;
    Real maxdiff = 0;

    if( props.verbose >= 2 )
        cerr << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl;

    if( props.cavity == Properties::CavityType::UNIFORM )
        Bi = Factor(delta(i));
    else {
        InfAlg *cav = newInfAlg( name, *this, opts );
        cav->makeCavity( i );

        if( props.cavity == Properties::CavityType::FULL )
            Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
        else if( props.cavity == Properties::CavityType::PAIR ) {
            vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, false );
            for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
                Bi *= pairbeliefs[ij];
        } else if( props.cavity == Properties::CavityType::PAIR2 ) {
            vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, true );
            for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
                Bi *= pairbeliefs[ij];
        }
        maxdiff = cav->maxDiff();
        delete cav;
    }
    Bi.normalize();
    _cavitydists[i] = Bi;

    return maxdiff;
}

int main( int argc, char** argv ) {
    if( argc != 4 )
        usage("Incorrect number of arguments.");

    FactorGraph fg;
    fg.ReadFromFile( argv[1] );

    PropertySet infprops;
    infprops.set( "verbose", (size_t)0 );
    infprops.set( "updates", string("HUGIN") );
    InfAlg* inf = newInfAlg( "JTREE", fg, infprops );
    inf->init();

    Evidence e;
    ifstream estream( argv[2] );
    e.addEvidenceTabFile( estream, fg );

    cout << "Number of samples: " << e.nrSamples() << endl;
    for( Evidence::iterator ps = e.begin(); ps != e.end(); ps++ )
        cout << "Sample #" << (ps - e.begin()) << " has " << ps->size() << " observations." << endl;

    ifstream emstream( argv[3] );
    EMAlg em(e, *inf, emstream);

    while( !em.hasSatisfiedTermConditions() ) {
        Real l = em.iterate();
        cout << "Iteration " << em.Iterations() << " likelihood: " << l <<endl;
    }

    cout << endl << "Inferred Factor Graph:" << endl << "######################" << endl;
    cout.precision(12);
    cout << inf->fg();

    delete inf;

    return 0;
}

Real EM_estep(MaximizationStep &mstep, const Evidence &evidence, InfAlg &inf) {
	Real likelihood = 0;

	inf.run();
	Real logZ = inf.logZ();

	// Expectation calculation
	for (Evidence::const_iterator e = evidence.begin(); e != evidence.end(); ++e) {
		InfAlg* clamped = inf.clone();
		// Apply evidence
		for (Evidence::Observation::const_iterator i = e->begin(); i != e->end(); ++i)
			clamped->clamp(clamped->fg().findVar(i->first), i->second);
		clamped->init();
		clamped->run();

		likelihood += clamped->logZ() - logZ;

		mstep.addExpectations(*clamped);

		delete clamped;
	}

	return likelihood;
}

Real BBPCostFunction::evaluate( const InfAlg &ia, const vector<size_t> *stateP ) const {
    Real cf = 0.0;
    const FactorGraph &fg = ia.fg();

    switch( (size_t)(*this) ) {
        case CFN_BETHE_ENT: // ignores state
            cf = -ia.logZ();
            break;
        case CFN_VAR_ENT: // ignores state
            for( size_t i = 0; i < fg.nrVars(); i++ )
                cf += -ia.beliefV(i).entropy();
            break;
        case CFN_FACTOR_ENT: // ignores state
            for( size_t I = 0; I < fg.nrFactors(); I++ )
                cf += -ia.beliefF(I).entropy();
            break;
        case CFN_GIBBS_B:
        case CFN_GIBBS_B2:
        case CFN_GIBBS_EXP: {
            DAI_ASSERT( stateP != NULL );
            vector<size_t> state = *stateP;
            DAI_ASSERT( state.size() == fg.nrVars() );
            for( size_t i = 0; i < fg.nrVars(); i++ ) {
                Real b = ia.beliefV(i)[state[i]];
                switch( (size_t)(*this) ) {
                    case CFN_GIBBS_B:
                        cf += b;
                        break;
                    case CFN_GIBBS_B2:
                        cf += b * b / 2.0;
                        break;
                    case CFN_GIBBS_EXP:
                        cf += exp( b );
                        break;
                    default:
                        DAI_THROW(UNKNOWN_ENUM_VALUE);
                }
            }
            break;
        } case CFN_GIBBS_B_FACTOR:
          case CFN_GIBBS_B2_FACTOR:
          case CFN_GIBBS_EXP_FACTOR: {
            DAI_ASSERT( stateP != NULL );
            vector<size_t> state = *stateP;
            DAI_ASSERT( state.size() == fg.nrVars() );
            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
                size_t x_I = getFactorEntryForState( fg, I, state );
                Real b = ia.beliefF(I)[x_I];
                switch( (size_t)(*this) ) {
                    case CFN_GIBBS_B_FACTOR:
                        cf += b;
                        break;
                    case CFN_GIBBS_B2_FACTOR:
                        cf += b * b / 2.0;
                        break;
                    case CFN_GIBBS_EXP_FACTOR:
                        cf += exp( b );
                        break;
                    default:
                        DAI_THROW(UNKNOWN_ENUM_VALUE);
                }
            }
            break;
        } default:
            DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(*this));
    }
    return cf;
}

//.........这里部分代码省略.........
        for( size_t ev = 0; ev < evid.size(); ev++ )
            bestMPE[ev].state = stateVec( fgs[ev].nrVars(), 0 );
        vector<size_t> ev2go;
        ev2go.reserve( evid.size() );
        for( size_t ev = 0; ev < evid.size(); ev++ )
            ev2go.push_back( ev );

        // solve inference problems
        if( verbose )
            cout << "Solving inference problems..." << endl;
        bool first = true;
        size_t nrsolvers = 3;
        vector<size_t> nrsubsolvers( nrsolvers );
        nrsubsolvers[0] = 2;
        nrsubsolvers[1] = 1;
        nrsubsolvers[2] = 1;
        double MPEdamping = 0.49;
        // for each (sub)solver
        for( size_t solver = 0; solver < nrsolvers; solver++ ) {
            if( verbose )
                cout << "  Solver " << solver << endl;

            // for each evidence case
            size_t subsolver = 0;
            for( long _ev = 0; _ev < (long)ev2go.size(); ) {
                bool improved = false;
                size_t ev = ev2go[_ev];
                if( verbose )
                    cout << "    Evidence case " << ev << ", subsolver = " << subsolver << "..." << endl;

                // construct inference algorithm on clamped factor graph
                if( verbose )
                    cout << "      Constructing inference algorithm..." << endl;
                InfAlg *ia = NULL;
                double tic = toc();
                bool failed = false;
                try {
                    // construct
                    if( solver == 0 ) { // the quick one
                        double remtime = (UAI_time - (toc() - starttic)) * 0.9;
                        if( remtime < 1.0 )
                            remtime = 1.0 ;
                        double maxtime = remtime / (ev2go.size() - _ev);
                        if( verbose ) {
                            cout << "      Past time:     " << (toc() - starttic) << endl;
                            cout << "      Remaining time:" << remtime << endl;
                            cout << "      Allotted time: " << maxtime << endl;
                        }
                        string maxtimestr;
                        if( maxtime != INFINITY )
                            maxtimestr = ",maxtime=" + toString(maxtime);
                        // quick and dirty...
                        if( task == "MPE" )
                            ia = newInfAlgFromString( "BP[inference=MAXPROD,updates=SEQRND,logdomain=" + toString(subsolver) + ",tol=1e-9,maxiter=10000" + maxtimestr + ",damping=0.1,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
                        else {
                            if( couldBeGrid )
                                ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9,maxiter=10000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
                            else
                                ia = newInfAlgFromString( "BP[inference=SUMPROD,updates=SEQRND,logdomain=" + toString(subsolver) + ",tol=1e-9,maxiter=10000" + maxtimestr + ",damping=0.0,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
                        }
                    } else if( solver == 1 ) { // the exact one
                        string maxmemstr;
                        if( UAI_memory != 0 )
                            maxmemstr = ",maxmem=" + toString(UAI_memory);
                        if( task == "MPE" )
                            ia = newInfAlgFromString( "JTREE[inference=MAXPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );

 void doDAI() {
     double tic = toc();
     if( obj != NULL ) {
         obj->init();
         obj->run();
         time += toc() - tic;
         try {
             logZ = obj->logZ();
             has_logZ = true;
         } catch( Exception &e ) {
             has_logZ = false;
         }
         try {
             maxdiff = obj->maxDiff();
             has_maxdiff = true;
         } catch( Exception &e ) {
             has_maxdiff = false;
         }
         try {
             iters = obj->Iterations();
             has_iters = true;
         } catch( Exception &e ) {
             has_iters = false;
         }
         q = allBeliefs();
     };
 }

 vector<Factor> allBeliefs() {
     vector<Factor> result;
     for( size_t i = 0; i < obj->fg().nrVars(); i++ )
         result.push_back( obj->belief( obj->fg().var(i) ) );
     return result;
 }

 string identify() { 
     if( obj != NULL )
         return obj->identify(); 
     else
         return "NULL";
 }