CFactorGraph CompressInterface::createCFactorGraph() {

	createVarClustering();
	createFacClustering();

	// create lifted fg here
	vector<CFactor> superFacs;

	// clusterIdx => facIdx
	for (map<size_t,size_t>::iterator facIter = _facRepr.begin(); facIter != _facRepr.end(); facIter++) {
		VarSet superVarSet;

		foreach (const dai::BipartiteGraph::Neighbor &tmpVar, _cfg.nbF(facIter->second)) {
			Var varCluster = _varRepr[_varColorVec[tmpVar]];
			if (!superVarSet.contains(varCluster)) {
				superVarSet |= Var(varCluster);
			}
		}

		CFactor superFac = CFactor(superVarSet, _cfg.factor(facIter->second).p());
		superFac.sigma() = _cfg.factor(facIter->second).sigma();
		superFac.position() = _cfg.factor(facIter->second).position();

		superFac.counts() = createCounts(facIter->second, superVarSet);
		superFacs.push_back(superFac);
	}

	return CFactorGraph(superFacs);
}

Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
    size_t I = nbV(i)[_I];
    Factor piet = _pancakes[i];

    // recalculate _pancake[i]
    VarSet Ivars = factor(I).vars();
    Factor A_I;
    for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
        if( var(i) != *k )
            A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).marginal( Ivars / var(i), false );
    if( Ivars.size() > 1 )
        A_I ^= (1.0 / (Ivars.size() - 1));
    Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).marginal( Ivars / var(i), false );
    Factor quot = A_I / A_Ii;
    if( props.damping != 0.0 )
        quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);

    piet *= quot / _phis[i][_I].normalized();
    _phis[i][_I] = quot.normalized();

    piet.normalize();

    if( piet.hasNaNs() ) {
        cerr << name() << "::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
        hasNaNs = true;
    }

    return piet;
}

/**
* @brief Returns all local variables, including parameters.
*/
VarSet VarsHandler::getLocalVars() const {
	VarSet result;
	for (const auto &p : localVars) {
		result.insert(p.second);
	}
	return result;
}

void MakeHOIFG( size_t N, size_t M, size_t k, Real sigma, FactorGraph &fg ) {
    vector<Var> vars;
    vector<Factor> factors;

    vars.reserve(N);
    for( size_t i = 0; i < N; i++ )
        vars.push_back(Var(i,2));

    for( size_t I = 0; I < M; I++ ) {
        VarSet vars;
        while( vars.size() < k ) {
            do {
                size_t newind = (size_t)(N * rnd_uniform());
                Var newvar = Var(newind, 2);
                if( !vars.contains( newvar ) ) {
                    vars |= newvar;
                    break;
                }
            } while( 1 );
        }
        factors.push_back( RandomFactor( vars, sigma ) );
    }

    fg = FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
}

/** \param N number of variables
 *  \param M number of factors
 *  \param k number of variables that each factor depends on
 *  \param beta standard-deviation of Gaussian log-factor entries
 */
FactorGraph createHOIFG( size_t N, size_t M, size_t k, Real beta ) {
    vector<Var> vars;
    vector<Factor> factors;

    vars.reserve(N);
    for( size_t i = 0; i < N; i++ )
        vars.push_back(Var(i,2));

    for( size_t I = 0; I < M; I++ ) {
        VarSet vars;
        while( vars.size() < k ) {
            do {
                size_t newind = (size_t)(N * rnd_uniform());
                Var newvar = Var(newind, 2);
                if( !vars.contains( newvar ) ) {
                    vars |= newvar;
                    break;
                }
            } while( 1 );
        }
        factors.push_back( createFactorExpGauss( vars, beta ) );
    }

    return FactorGraph( factors.begin(), factors.end(), vars.begin(), vars.end(), factors.size(), vars.size() );
}

VarSet WhySemiring::getVars() const {
  VarSet res;
  for (auto ws : val) {
    res.insert(ws.begin(), ws.end());
  }
  return res;
}

Index *IncLit::index(Grounder *, Formula *, VarSet &bound)
{
	VarSet vars;
	VarVec bind;
	var_->vars(vars);
	std::set_difference(vars.begin(), vars.end(), bound.begin(), bound.end(), std::back_insert_iterator<VarVec>(bind));
	bound.insert(bind.begin(), bind.end());
	return new IncIndex(this);
}

void findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
    for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
        VarSet ind = fg.delta( *in );
        if( (newcl.size()) >= 2 && ind.contains( root ) ) {
            allcl.insert( newcl | *in );
        }
        else if( length > 1 )
            findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
    }
}

Factor LC::belief (const VarSet &ns) const {
    if( ns.size() == 0 )
        return Factor();
    else if( ns.size() == 1 )
        return beliefV( findVar( *(ns.begin()) ) );
    else {
        DAI_THROW(BELIEF_NOT_AVAILABLE);
        return Factor();
    }
}

void SharedParameters::setPermsAndVarSetsFromVarOrders() {
    if( _varorders.size() == 0 )
        return;
    DAI_ASSERT( _estimation != NULL );

    // Construct the permutation objects and the varsets
    for( FactorOrientations::const_iterator foi = _varorders.begin(); foi != _varorders.end(); ++foi ) {
        VarSet vs;
        _perms[foi->first] = calculatePermutation( foi->second, vs );
        _varsets[foi->first] = vs;
        DAI_ASSERT( _estimation->probSize() == vs.nrStates() );
    }
}

Factor ExactInf::belief( const VarSet &ns ) const {
    if( ns.size() == 0 )
        return Factor();
    else if( ns.size() == 1 ) {
        return beliefV( findVar( *(ns.begin()) ) );
    } else {
        size_t I;
        for( I = 0; I < nrFactors(); I++ )
            if( factor(I).vars() >> ns )
                break;
        if( I == nrFactors() )
            DAI_THROW(BELIEF_NOT_AVAILABLE);
        return beliefF(I).marginal(ns);
    }
}

WhySemiring::WhySemiring(std::string str_val)
{
  if(0 == str_val.compare("1")) {
    this->val = WhySet();
    val.insert(VarSet());
  }
  else if (0 == str_val.compare("0")) {
    this->val = WhySet();
  }
  else {
    this->val = WhySet();
    VarSet tmp = VarSet();
    tmp.insert(Var::GetVarId(str_val));
    val.emplace(tmp);
  }
}

/**
* @brief Precomputes @c funcInfo->isAlwaysModifiedBeforeRead for @c traversedFunc.
*/
void OptimFuncInfoCFGTraversal::precomputeAlwaysModifiedVarsBeforeRead() {
	// Initialization.
	funcInfo->varsAlwaysModifiedBeforeRead.clear();
	// Global variables which are read during the computation.
	VarSet readVars;

	// Currently, we only traverse the function's body up to the first compound
	// statement. Moreover, we only consider global variables as the computed
	// piece of information is useless for local variables.
	// TODO Use a CFG traversal for this to improve the analysis.
	ShPtr<Statement> stmt(traversedFunc->getBody());
	while (stmt) {
		if (stmt->isCompound()) {
			// As we are not using a CFG traversal, this is the end of the
			// computation.
			break;
		}

		ShPtr<ValueData> stmtData(va->getValueData(stmt));

		// Handle directly read variables.
		addToSet(stmtData->getDirReadVars(), readVars);

		// Handle function calls (indirectly accessed variables).
		for (auto i = stmtData->call_begin(), e = stmtData->call_end(); i != e; ++i) {
			ShPtr<CallInfo> callInfo(cio->computeCallInfo(*i, traversedFunc));
			for (const auto &var : globalVars) {
				if (callInfo->mayBeRead(var)) {
					readVars.insert(var);
				}
			}
		}

		// Handle directly written variables.
		for (auto i = stmtData->dir_written_begin(), e = stmtData->dir_written_end();
				i != e; ++i) {
			if (hasItem(globalVars, *i) && !hasItem(readVars, *i)) {
				// This global variable is modified before read.
				funcInfo->varsAlwaysModifiedBeforeRead.insert(*i);
			}
		}

		// TODO What about indirectly accessed variables?

		stmt = stmt->getSuccessor();
	}
}

std::vector<std::map<size_t, int> > AnyPositionCnfCompress::createCounts(size_t &gndFactor, VarSet &superVarSet) {
	// create zero entries for each position
	map<long, map<size_t, int> > countMap;
	foreach (const dai::BipartiteGraph::Neighbor &tmpVar, _cfg.nbF(gndFactor)) {
		Var liftedVar = _varRepr[_varColorVec[tmpVar]];
		size_t pos = find(_cfg.factor(gndFactor).sigma().begin(), _cfg.factor(gndFactor).sigma().end(), tmpVar.iter) - _cfg.factor(gndFactor).sigma().begin();
		countMap[liftedVar.label()][pos] = 0;
	}

	vector<map<size_t, int> > counts;
	size_t posCount;
	size_t negCount;
	for (vector<Var>::const_iterator iter = superVarSet.begin(); iter < superVarSet.end(); iter++) {
		posCount = 0;
		negCount = 0;

		foreach(const dai::BipartiteGraph::Neighbor tmpFac, _cfg.nbV(_cfg.findVar(*iter))) {
			if (_facRepr[_facColorVec[tmpFac]] == gndFactor) {
				size_t pos = find(_cfg.factor(tmpFac).sigma().begin(), _cfg.factor(tmpFac).sigma().end(), tmpFac.dual) - _cfg.factor(tmpFac).sigma().begin();
				double res = log(_cfg.factor(tmpFac).states() - _zeroStates[tmpFac]) /  log(2);
				size_t nrPosLiterals = size_t(res);
				bool sign = (pos < nrPosLiterals);
				if (sign) {
					posCount++;
				} else {
					negCount++;
				}
			}
		}
		map<size_t, int>::iterator posIter=countMap[iter->label()].begin();

		if (posCount > 0) {
			posIter->second = posCount;
		}

		if (negCount > 0) {
			for (size_t j=0;j<posCount; j++, posIter++) {}
			posIter->second = negCount;
		}

		counts.push_back(countMap[iter->label()]);
	}

	return counts;
}

void RIFActVarBind::getVars(VarSet &vars, const bool include_act_var) const
    throw() {
  if (include_act_var) {
    vars.insert(this->var);
  }
  if (this->frame != NULL) {
    this->frame->getVars(vars);
  }
}

std::vector<std::map<size_t, int> > CompressInterface::createCounts(size_t &gndFactor, VarSet &superVarSet) {
	// create zero entries for each position
	map<long, map<size_t, int> > countMap;
	foreach (const dai::BipartiteGraph::Neighbor &tmpVar, _cfg.nbF(gndFactor)) {
		Var liftedVar = _varRepr[_varColorVec[tmpVar]];
		size_t pos = find(_cfg.factor(gndFactor).sigma().begin(), _cfg.factor(gndFactor).sigma().end(), tmpVar.iter) - _cfg.factor(gndFactor).sigma().begin();
		countMap[liftedVar.label()][pos] = 0;
	}

	vector<map<size_t, int> > counts;
	for (vector<Var>::const_iterator iter = superVarSet.begin(); iter < superVarSet.end(); iter++) {
		foreach(const dai::BipartiteGraph::Neighbor gndFac, _cfg.nbV(_cfg.findVar(*iter))) {
			if (_facRepr[_facColorVec[gndFac]] == gndFactor) {
				size_t pos = find(_cfg.factor(gndFac).sigma().begin(), _cfg.factor(gndFac).sigma().end(), gndFac.dual) - _cfg.factor(gndFac).sigma().begin();
				countMap[iter->label()][pos]++;
			}
		}
		counts.push_back(countMap[iter->label()]);
	}
	return counts;
}

void Function::hansen_matrix(const IntervalVector& box, IntervalMatrix& H, const VarSet& set) const {
	int n=set.nb_var;
	int m=image_dim();

	assert(H.nb_cols()==n);
	assert(box.size()==nb_var());
	assert(H.nb_rows()==m);

	IntervalVector var_box=set.var_box(box);
	IntervalVector param_box=set.param_box(box);

	IntervalVector x=var_box.mid();
	IntervalMatrix J(m,n);

	for (int var=0; var<n; var++) {
		//var=tab[i];
		x[var]=var_box[var];
		jacobian(set.full_box(x,param_box),J,set);
		H.set_col(var,J.col(var));
	}
}

void Function::jacobian(const IntervalVector& box, IntervalMatrix& J, const VarSet& set) const {

	assert(J.nb_cols()==set.nb_var);
	assert(box.size()==nb_var());
	assert(J.nb_rows()==image_dim());

	IntervalVector g(nb_var());

	// calculate the gradient of each component of f
	for (int i=0; i<image_dim(); i++) {
		(*this)[i].gradient(box,g);
		J.set_row(i,set.var_box(g));
	}
}

void RIFTerm::getVars(VarSet &vars) const throw() {
  if (this->type == VARIABLE) {
    RIFVar *v = (RIFVar*)this->state;
    RIFVar v2 = *v;
    vars.insert(v2);
  } else if (this->type == FUNCTION) {
    func_state *f = (func_state*) this->state;
    if (f->args != NULL) {
      RIFTerm *mark = f->args->dptr();
      RIFTerm *end = mark + f->args->size();
      for (; mark != end; ++mark) {
        mark->getVars(vars);
      }
    }
  }
}

void ParityAggrLit::index(Grounder *g, Groundable *gr, VarSet &bound)
{
	(void)g;
	if(assign_)
	{
		VarSet vars;
		VarVec bind;
		lower_->vars(vars);
		std::set_difference(vars.begin(), vars.end(), bound.begin(), bound.end(), std::back_insert_iterator<VarVec>(bind));
		if(bind.size() > 0)
		{
			bound.insert(bind.begin(), bind.end());
			return;
		}
	}
	gr->instantiator()->append(new MatchIndex(this));
}