C++ BitSet::size方法代码示例

Delta delta( const BitSet& bit,
             const Scalar& default_value ) {
    Delta u( bit.size(), 1 );
    for( uint i = 0; i < bit.size(); ++i ) {
        if( bit[i]) {
            u.insert( i, 0 ) = default_value;
        }
    }
    return u;
}

FncKhunTucker::FncKhunTucker(const NormalizedSystem& sys, Function& _df, Function** _dg, const IntervalVector& current_box, const BitSet& active) :
								Fnc(1,1), nb_mult(0), // **tmp**
								n(sys.nb_var), sys(sys), df(_df), dg(active.size()), active(active),
								eq(BitSet::empty(sys.f_ctrs.image_dim())),
								ineq(BitSet::empty(sys.f_ctrs.image_dim())),
								bound_left(BitSet::empty(sys.nb_var)),
								bound_right(BitSet::empty(sys.nb_var)) {

	int l=1; // lambda0

	if (_dg==NULL && !active.empty()) {
		ibex_error("[FncKhunTucker] an unconstrained system cannot have active constraints");
	}

	unsigned int i=0; // index of a constraint in the active set
	for (BitSet::const_iterator c=active.begin(); c!=active.end(); ++c, ++i) {
		dg.set_ref(i,*_dg[c]);
		if (sys.ops[c]==EQ) eq.add(i);
		else ineq.add(i);
		//cout << " constraint n°" << c << " active\n";
	}
	l+=active.size();

	for (int j=0; j<sys.box.size(); j++) {
		if (current_box[j].lb() <= sys.box[j].lb()) {
			bound_left.add(j);
			//cout << " left bound n°" << j << " active\n";
			l++;
		}
		if (current_box[j].ub() >= sys.box[j].ub()) {
			bound_right.add(j);
			//cout << " right bound n°" << j << " active\n";
			l++;
		}

	}

	(int&) nb_mult = l;

	(int&) _nb_var = n + l;

	assert(_nb_var == sys.nb_var + eq.size() + ineq.size() + bound_left.size() + bound_right.size() + 1);

	(Dim&) _image_dim = Dim(_nb_var, 1);

}

void BronKerbosch::bkClassic(BitSet P, BitSet R, BitSet X)
{
  assert((P & X).none());
  assert((P & R).none());
  assert((R & X).none());

  // let's print P, R and X
  //std::cout << "P = ";
  //print(P, std::cout);
  //std::cout << ", R = ";
  //print(R, std::cout);
  //std::cout << ", X = ";
  //print(X, std::cout);
  //std::cout << std::endl;

  // Reports maximal cliques in P \cup R (but not in X)
  BitSet P_cup_X = P | X;
  if (P_cup_X.none())
  {
    report(R);
  }
  else
  {
    for (size_t i = 0; i < P.size(); ++i)
    {
      if (P[i])
      {
        Node v = _bitToNode[i];
        BitSet R_ = R;
        R_[_nodeToBit[v]] = 1;
        // report all maximal cliques in ( (P | N[v]) & R) \ (X & N[v]) )
        bkClassic(P & _bitNeighborhood[v], R_, X & _bitNeighborhood[v]);
        P[i] = 0;
        X[i] = 1;
      }
    }
  }

  /* Invariants:
   * - R is a clique
   * - Each node v in P is adjacent to all nodes in R
   *   (nodes in P are used to extend R, upon usage it's moved to X)
   * - Each node v in X is adjacent to all nodes in R
   *   (maximal cliques containing R \cup v have already been reported)
   */
}

void BronKerbosch::report(const BitSet& R)
{
  NodeVector clique;
  for (size_t i = 0; i < R.size(); ++i)
  {
    if (R[i])
    {
      if (g_verbosity >= VERBOSE_DEBUG)
      {
        std::cerr << " " << _g.id(_bitToNode[i]);
      }
      clique.push_back(_bitToNode[i]);
    }
  }
  _cliques.push_back(clique);
  if (g_verbosity >= VERBOSE_DEBUG)
    std::cerr << std::endl;
}

void BronKerbosch::printBitSet(const BitSet& S, std::ostream& out) const
{
  out << "{";
  bool first = true;
  for (size_t i = 0; i < S.size(); ++i)
  {
    if (S[i])
    {
      if (!first)
      {
        out << ", ";
      }
      else
      {
        first = false;
      }
      out << _g.id(_bitToNode[i]);
    }
  }
  out << "}";
}

int main()
{
  BitSet<8> a;
  a.set(0);
  a.set(1);
  a.set(2);

  BitSet<8> b;
  b.set(3);

  cout << a.intersects(b) << endl;

  b.set(1);

  cout << a.intersects(b) << endl;

  a.set(451);

  cout << a.size() << endl;

  return 0;
}

void
CLuceneIndexWriter::deleteEntry(const string& entry,
        lucene::index::IndexReader* reader) {
    int deleted = 0;
    wstring path(utf8toucs2(entry));
{
    Term t(systemlocation(), path.c_str());
    deleted += reader->deleteDocuments(&t);
    // if no file was deleted, no more can be deleted
    if (deleted == 0) return;
}
{
    Term t(parentlocation(), path.c_str());
    deleted += reader->deleteDocuments(&t);

    // if we have only deleted one file up to now, we cannot delete more
    if (deleted < 2) return;
}
{
    // delete all deeper nested files
    wstring v = utf8toucs2(entry+"/");
    Term* t = _CLNEW Term(parentlocation(), v.c_str());
    PrefixFilter* filter = _CLNEW PrefixFilter(t);
    BitSet* b = filter->bits(reader);
    _CLDELETE(filter);
    _CLDECDELETE(t);
    int32_t size = b->size();
    for (int id = 0; id < size; ++id) {
        if (b->get(id) && !reader->isDeleted(id)) {
            reader->deleteDocument(id);
            deleted++;
        }
    }
    _CLDELETE(b);
}
}

void FncKhunTucker::jacobian(const IntervalVector& x_lambda, IntervalMatrix& J, const BitSet& components, int v) const {

	if (components.size()!=n+nb_mult) {
		not_implemented("FncKhunTucker: 'jacobian' for selected components");
		//J.resize(n+nb_mult,n+nb_mult);
	}

	IntervalVector x=x_lambda.subvector(0,n-1);

	int lambda0=n;	// The index of lambda0 in the box x_lambda is nb_var.

	int l=lambda0; // mutipliers indices counter. The first multiplier is lambda0.

	// matrix corresponding to the "Hessian expression" lambda_0*d^2f+lambda_1*d^2g_1+...=0
	IntervalMatrix hessian=x_lambda[l] * df.jacobian(x,v<n? v : -1); // init
	if (v==-1 || v==l) J.put(0, l, df.eval_vector(x), false);

	l++;

	IntervalVector gx;
	if (!active.empty())
		gx = sys.f_ctrs.eval_vector(x,active);

	// normalization equation (init)
	J[lambda0].put(0,Vector::zeros(n));
	J[lambda0][lambda0]=1.0;

	IntervalVector dgi(n); // store dg_i([x]) (used in several places)

	for (BitSet::const_iterator i=ineq.begin(); i!=ineq.end(); ++i) {
		hessian += x_lambda[l] * dg[i].jacobian(x,v<n? v : -1);
		dgi=dg[i].eval_vector(x);
		J.put(0, l, dgi, false);

		J.put(l, 0, (x_lambda[l]*dgi), true);
		J.put(l, n, Vector::zeros(nb_mult), true);
		J[l][l]=gx[i];

		J[lambda0][l] = 1.0;

		l++;
	}

	for (BitSet::const_iterator i=ineq.begin(); i!=ineq.end(); ++i) {
		hessian += x_lambda[l] * dg[i].jacobian(x,v<n? v : -1);
		dgi=dg[i].eval_vector(x);
		J.put(0, l, dgi, false);

		J.put(l, 0, dgi, true);
		J.put(l, n, Vector::zeros(nb_mult), true);

		J[lambda0][l] = 2*x_lambda[l];

		l++;
	}

	for (BitSet::const_iterator v=bound_left.begin(); v!=bound_left.end(); ++v) {
		// this constraint does not contribute to the "Hessian expression"
		dgi=Vector::zeros(n);
		dgi[v]=-1.0;
		J.put(0, l, dgi, false);

		J.put(l, 0, (x_lambda[l]*dgi), true);
		J.put(l, n, Vector::zeros(nb_mult), true);
		J[l][l]=(-x[v]+sys.box[v].lb());

		J[lambda0][l] = 1.0;
		l++;
	}


	for (BitSet::const_iterator v=bound_right.begin(); v!=bound_right.end(); ++v) {
		// this constraint does not contribute to the "Hessian expression"
		dgi=Vector::zeros(n);
		dgi[v]=1.0;
		J.put(0, l, dgi, false);

		J.put(l, 0, (x_lambda[l]*dgi), true);
		J.put(l, n, Vector::zeros(nb_mult), true);
		J[l][l]=(x[v]-sys.box[v].ub());

		J[lambda0][l] = 1.0;
		l++;
	}

	assert(l==nb_mult+n);

	J.put(0,0,hessian);

}

IntervalVector FncKhunTucker::eval_vector(const IntervalVector& x_lambda, const BitSet& components) const {

	if (components.size()!=n+nb_mult) {
		not_implemented("FncKhunTucker: 'eval_vector' for selected components");
		//J.resize(n+nb_mult,n+nb_mult);
	}

	IntervalVector res(n+nb_mult);

	// Variables in x_lambda are organized as follows:
	// x - lambda0 - mu - lambda - beta

	IntervalVector x=x_lambda.subvector(0,n-1);

	int lambda0=n;	// The index of lambda0 in the box x_lambda is nb_var.

	int l=lambda0; // multipliers indices counter. The first multiplier is lambda0.

	// vector corresponding to the "gradient expression" lambda_0*dg + lambda_1*dg_1 + ... (init
	IntervalVector grad=x_lambda[l] * df.eval_vector(x); // init

	l++;

	IntervalVector gx;

	if (!active.empty()) {
		gx=sys.f_ctrs.eval_vector(x,active);
	}
	// normalization equation lambda_0 + ... = 1.0
	res[lambda0] = x_lambda[lambda0] - 1.0; // init

	for (BitSet::const_iterator i=ineq.begin(); i!=ineq.end(); ++i) {
		grad += x_lambda[l] * dg[i].eval_vector(x);
		res[l] = x_lambda[l] * gx[i];
		res[lambda0] += x_lambda[l];
		l++;
	}

	for (BitSet::const_iterator i=eq.begin(); i!=eq.end(); ++i) {
		grad += x_lambda[l] * dg[i].eval_vector(x);
		res[l] = gx[i];
		res[lambda0] += sqr(x_lambda[l]);
		l++;
	}

	for (BitSet::const_iterator v=bound_left.begin(); v!=bound_left.end(); ++v) {
		grad[v] -= x_lambda[l];
		res[l] = x_lambda[l] * (-x[v]+sys.box[v].lb());
		res[lambda0] += x_lambda[l];
		l++;
	}

	for (BitSet::const_iterator v=bound_right.begin(); v!=bound_right.end(); ++v) {
		grad[v] += x_lambda[l];
		res[l] = x_lambda[l] * (x[v]-sys.box[v].ub());
		res[lambda0] += x_lambda[l];
		l++;
	}

	assert(l==nb_mult+n);

	res.put(0, grad);

	return res;
}

CtcAcid::CtcAcid(const System& sys, const BitSet& cid_vars, Ctc& ctc, bool optim, int s3b, int scid,
		double var_min_width, double ct_ratio): Ctc3BCid (cid_vars,ctc,s3b,scid,cid_vars.size(),var_min_width),
		system(sys), nbcalls(0), nbctvar(0), ctratio(ct_ratio),  nbcidvar(0), nbtuning(0), optim(optim)  {
	// [gch] BNE check the argument "cid_vars.nb_set()" given to _3BCID
}

int LinearizerDuality::linearize(const IntervalVector& box, LPSolver& lp_solver, BoxProperties& prop)  {
	// ========= get active constraints ===========
	/* Using system cache seems not interesting. */
	//BxpSystemCache* cache=(BxpSystemCache*) prop[BxpSystemCache::get_id(sys)];
	BxpSystemCache* cache=NULL;
	//--------------------------------------------------------------------------

	BitSet* active;

	if (cache!=NULL) {
		active = &cache->active_ctrs();
	} else {
		active = new BitSet(sys.active_ctrs(box));
	}
	// ============================================

	size_t n = sys.nb_var;
	size_t m = sys.f_ctrs.image_dim();
	size_t n_total = n +  m*n;

	int nb_ctr=0; // number of inequalities added in the LP solver

//	BxpLinearRelaxArgMin* argmin=(BxpLinearRelaxArgMin*) prop[BxpLinearRelaxArgMin::get_id(sys)];
//
//	if (argmin && argmin->argmin()) {
//		pt=*argmin->argmin();
//	} else
		pt=box.mid();

	if (!active->empty()) {

		//IntervalMatrix J=cache? cache->active_ctrs_jacobian() : sys.f_ctrs.jacobian(box,active);
		//IntervalMatrix J=sys.f_ctrs.jacobian(box,*active);
		IntervalMatrix J(active->size(),n); // derivatives over the box
		sys.f_ctrs.hansen_matrix(box,pt,J,*active);

		if (J.is_empty()) {
			if (cache==NULL) delete active;
			return -1;
		}

		// the evaluation of the constraints in the point
		IntervalVector gx(sys.f_ctrs.eval_vector(pt,*active));
		if (gx.is_empty()) {
			if (cache==NULL) delete active;
			return 0;
		}


		int i=0; // counter of active constraints
		for (BitSet::iterator c=active->begin(); c!=active->end(); ++c, i++)  {

			if (!sys.f_ctrs.deriv_calculator().is_linear[c]) {
				for (size_t j=0; j<n; j++) {
					Vector row(n_total,0.0);
					row[j]=1;
					row[n + c*n +j]=1;

					double rhs = pt[j] - lp_solver.get_epsilon();

					lp_solver.add_constraint(row, LEQ, rhs);
					nb_ctr++;
				}
			}

			Vector row(n_total,0.0);
			row.put(0,J[i].lb());

			IntervalVector gl(J[i].lb());

			Vector diam_correctly_rounded = (IntervalVector(J[i].ub())-gl).lb();

			for (size_t j=0; j<n; j++) {
				if (diam_correctly_rounded[j]<0)
					ibex_error("negative diameter");
			}

			row.put(n + c*n,-diam_correctly_rounded);

			double rhs = (-gx[i] + (gl*pt)).lb()- lp_solver.get_epsilon();

			lp_solver.add_constraint(row, LEQ, rhs);
			nb_ctr++;
		}
	}

	if (cache==NULL) delete active;
	return nb_ctr;
}

	/*
	 * Create MyCtcExist. The number of variables
	 * is the number of bits set to "true" in the
	 * "vars" structure; this number is obtained
	 *  via vars.size().
	 */
	MyCtcExist(Ctc& c, const BitSet& vars, const IntervalVector& box_y) :
		Ctc(vars.size()), c(c), vars(vars), box_y(box_y) { }

void BronKerbosch::bkPivot(BitSet P, BitSet R, BitSet X)
{
  assert((P & X).none());
  assert((P & R).none());
  assert((R & X).none());

  // let's print P, R and X
  //std::cout << "P = ";
  //print(P, std::cout);
  //std::cout << ", R = ";
  //print(R, std::cout);
  //std::cout << ", X = ";
  //print(X, std::cout);
  //std::cout << std::endl;

  // Reports maximal cliques in P \cup R (but not in X)
  BitSet P_cup_X = P | X;
  if (P_cup_X.none())
  {
    report(R);
  }
  else
  {
    // choose a pivot u from (P | X) s.t |P & N(u)| is maximum, Tomita et al. (2006)
    size_t maxBitCount = 0;
    Node max_u = lemon::INVALID;
    for (size_t i = 0; i < P.size(); ++i)
    {
      if (P_cup_X[i])
      {
        Node u = _bitToNode[i];
        BitSet P_cap_Nu = P & _bitNeighborhood[u];
        size_t s = P_cap_Nu.count();
        if (s >= maxBitCount)
        {
          max_u = u;
          maxBitCount = s;
        }
      }
    }

    assert(max_u != lemon::INVALID);
    BitSet P_diff_Nu = P - _bitNeighborhood[max_u];
    for (size_t i = 0; i < P.size(); ++i)
    {
      if (P_diff_Nu[i])
      {
        Node v = _bitToNode[i];
        BitSet R_ = R;
        R_[_nodeToBit[v]] = 1;
        // report all maximal cliques in ( (P | N[v]) & R) \ (X & N[v]) )
        bkPivot(P & _bitNeighborhood[v], R_, X & _bitNeighborhood[v]);
        P[i] = 0;
        X[i] = 1;
      }
    }
  }

  /* Invariants:
   * - R is a clique
   * - Each node v in P is adjacent to all nodes in R
   *   (nodes in P are used to extend R, upon usage it's moved to X)
   * - Each node v in X is adjacent to all nodes in R
   *   (maximal cliques containing R \cup v have already been reported)
   */
}

Vector UnconstrainedLocalSearch::conj_grad(const Vector& gk, const Matrix& Bk, const Vector& xk, const Vector& x_gcp, const IntervalVector& region, const BitSet& I) {
	int hn = n-I.size(); // the restricted dimension

	//  cout << " [conj_grad] init x_gcp= " << x_gcp << endl;
	if (hn==0) return x_gcp;  // we are in a corner: nothing to do

	// ====================== STEP 3.0 : Initialization ======================
	Vector x=x_gcp; // next point, initialized to gcp

	// gradient of the quadratic model on zk1
	Vector r = -gk-Bk*(x-xk);

	double eta = get_eta(gk,xk,region,I);

	// Initialization of the conjuguate gradient restricted to the direction not(I[i])
	Vector hp(hn); // the restricted conjugate direction
	Vector hx(hn); // the restricted iterate
	Vector hr(hn); // the restricted gradient
	Vector hy(hn); // temporary vector
	Matrix hB(hn,hn); // the restricted hessian matrix
	IntervalVector hregion(hn); // the restricted region

	// initialization of \hat{B}:
	int p=0, q=0;
	for (int i=0; i<n; i++) {
		if (!I[i]) {
			for (int j=0; j<n; j++) {
				if (!I[j]) {
					hB[p][q] = Bk[i][j];
					q++;
				}
			}
			p++;
			q=0;
		}
	}

	// initialization of \hat{r} and \hat{region}
	p=0;
	for (int i=0; i<n; i++) {
		if (!I[i]) {
			hregion[p] = region[i];
			hr[p] = r[i];
			hx[p] = x[i];
			p++;
		}
	}

	double rho1 = 1.0; 					// norm of the restricted gradient at the previous iterate
	double rho2 = ::pow(hr.norm(),2); 	// norm of the restricted gradient

	// ====================== STEP 3.1 : Test for the required accuracy in c.g. iteration ======================
	bool cond = (rho2>::pow(eta,2));

	try {
		while (cond) {

			// ====================== STEP 3.2 : Conjugate gradient recurrences ======================

			// Update the restricted conjugate direction
			// \hat{p} = \hat{r} +beta* \hat{p}
			hp = hr + (rho2/rho1)*hp;

			// Update the temporary vector
			// \hat{y} = \hat{Bk}*\hat{p}
			hy = hB*hp;

			//  cout << " [conj_grad] current hr=" << hr << endl;
			//  cout << " [conj_grad] current hp=" << hp << endl;
			LineSearch ls(hregion,hx,hp,data,sigma);

			double alpha1=ls.alpha_max();
			//  cout << " [conj_grad] alpha1=" << alpha1 << endl;

			// first termination condition
			// we check if the hessian is "positive definite for \hat{p}"
			// otherwise the quadratic approximation is concave
			// and the solution if on the border of the region
			double aux = hp*hy;
			if ( aux <= 0) {
				cond = false;
				hx = ls.endpoint();
			}
			else {
				// second termination condition alpha2>alpha1
				double alpha2 = rho2/aux;

				if (alpha2>alpha1) {
					cond = false;
					hx = ls.endpoint();
				}
				else {
					// otherwise update x, r=\hat{r}, hy=y, rho1 and rho2= \hat{r}*\hat{r}=r*r
					hx += alpha2*hp;
					ls.proj(hx);
					hr -= alpha2*hy;
					rho1 = rho2;
					rho2 = hr*hr;
					cond = (rho2>(::pow(eta,2)));
				}
//.........这里部分代码省略.........

// returns number of elements in the bitset
static SB_INLINE uint32_t bitset_card(BitSet& bs) {
    return bs.size();
}