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

void CtcLMI::mkSolver(ibex::MatrixArray& matrices, const IntervalVector &box)
{
    // On modifie F0 en fonction des nouvelles bornes de boite
    for (int i = 0; i < box.size(); i++)
    {
        matrices[0][1 + box.size() + 2 * i][1 + box.size() + 2 * i] = box[i].lb();
        matrices[0][1 + box.size() + 2 * i + 1][1 + box.size() + 2 * i + 1] = -box[i].ub();
    }

    solver.setParameterType(SDPA::PARAMETER_DEFAULT);
    solver.setParameterMaxIteration(100);
    //    solver.setNumThreads(4);
    //    solver.printParameters(stdout);

    solver.inputConstraintNumber(matrices.size());
    solver.inputBlockNumber(1);

    solver.inputBlockSize(1, matrices[0].nb_cols());
    solver.inputBlockType(1, SDPA::SDP);

    solver.initializeUpperTriangleSpace();
    //    solver.inputCVec(1, 0);
    //    solver.inputCVec(2, -1);

    for (int i = 0; i < matrices.size(); i++)
        for (int j = 1; j < matrices[i].nb_cols() + 1; j++)
            for (int k = 1; k < matrices[i].nb_rows() + 1; k++)
                solver.inputElement(i, 1, k, j, matrices[i][k - 1][j - 1]);
}

void gauss_seidel(const IntervalMatrix& A, const IntervalVector& b, IntervalVector& x, double ratio) {
	int n=(A.nb_rows());
	assert(n == (A.nb_cols())); // throw NotSquareMatrixException();
	assert(n == (x.size()) && n == (b.size()));

	double red;
	Interval old, proj, tmp;

	do {
		red = 0;
		for (int i=0; i<n; i++) {
			old = x[i];
			proj = b[i];

			for (int j=0; j<n; j++)	if (j!=i) proj -= A[i][j]*x[j];
			tmp=A[i][i];

			bwd_mul(proj,tmp,x[i]);

			if (x[i].is_empty()) { x.set_empty(); return; }

			double gain=old.rel_distance(x[i]);
			if (gain>red) red=gain;
		}
	} while (red >= ratio);
}

void Function::gradient(const IntervalVector& x, IntervalVector& g) const {
	assert(g.size()==nb_var());
	assert(x.size()==nb_var());
	Gradient().gradient(*this,x,g);
//	if (!df) ((Function*) this)->df=new Function(*this,DIFF);
//	g=df->eval_vector(x);
}

bool inflating_gauss_seidel(const IntervalMatrix& A, const IntervalVector& b, IntervalVector& x, double min_dist, double mu_max) {
	int n=(A.nb_rows());
	assert(n == (A.nb_cols()));
	assert(n == (x.size()) && n == (b.size()));
	assert(min_dist>0);
	//cout << " ====== inflating Gauss-Seidel ========= " << endl;
	double red;
	IntervalVector xold(n);
	Interval proj;
	double d=DBL_MAX; // Hausdorff distances between 2 iterations
	double dold;
	double mu; // ratio of dist(x_k,x_{k-1)) / dist(x_{k-1},x_{k-2}).
	do {
		dold = d;
		xold = x;
		for (int i=0; i<n; i++) {
			proj = b[i];
			for (int j=0; j<n; j++)	if (j!=i) proj -= A[i][j]*x[j];
			x[i] = proj/A[i][i];
		}
		d=distance(xold,x);
		mu=d/dold;
		//cout << "  x=" << x << " d=" << d << " mu=" << mu << endl;
	} while (mu<mu_max && d>min_dist);
	//cout << " ======================================= " << endl;
	return (mu<mu_max);

}

void SepInverse::separate(IntervalVector& xin, IntervalVector& xout){

	assert(xin.size()==f.nb_var() && xout.size() == f.nb_var());


	xin &= xout;
	Domain tmp=f.eval_domain(xin);
	yin.init(Interval::ALL_REALS);
	yout.init(Interval::ALL_REALS);
	id->backward(tmp, yin);
	id->backward(tmp, yout);

	s.separate(yin, yout);

	if( yin.is_empty())
		xin.set_empty();
	else
		tmp = id->eval_domain(yin);
		f.backward(tmp, xin);

	if( yout.is_empty())
		xout.set_empty();
	else
		tmp = id->eval_domain(yout);
		f.backward(tmp, xout);

}

void midpoint(const IntervalVector& x, Eigen::VectorXd* x_mid)
{
	x_mid->resize(x.size());
	for (int i = 0; i < x.size(); ++i) {
		(*x_mid)[i] = (x[i].get_upper() + x[i].get_lower()) / 2.0;
	}
}

std::string print_mma(const IntervalVector& iv) {
	std::string res = "{";
	for(int i = 0; i < iv.size()-1; ++i) {
		res += "{" + std::to_string(iv[i].lb()) + ", " + std::to_string(iv[i].ub()) + "}, ";
	}
	res += "{" + std::to_string(iv[iv.size()-1].lb()) + ", " + std::to_string(iv[iv.size()-1].ub()) + "}}";
	return res;
}

bool TestIbex::almost_eq(const IntervalVector& y_actual, const IntervalVector& y_expected, double err) {
	if (y_actual.size()!=y_actual.size()) return false;
	if (y_actual.is_empty() && y_expected.is_empty()) return true;

	for (int i=0; i<y_actual.size(); i++) {
		if (!almost_eq(y_actual[i], y_expected[i],err)) return false;
	}

	return true;
}

void SetNodeReg::cut(const IntervalVector& box) {
	assert(is_leaf());
	unsigned int cutvar;
	if(father == NULL)
		cutvar = 0;
	else
		var = (father->var+1)%box.size();
	pt = box[var].mid();
	left = new SetNodeReg(status,this,(var+1)%box.size());
	right = new SetNodeReg(status,this,(var+1)%box.size());

}

Affine2Vector::Affine2Vector(const IntervalVector& x, bool b) :
		_n(x.size()),
		_vec(new Affine2[x.size()]) {
	if (!b) {
		for (int i = 0; i < x.size(); i++) {
			_vec[i] = Affine2(x[i]);
		}
	} else {
		for (int i = 0; i < x.size(); i++) {
			_vec[i] = Affine2(x.size(), i + 1, x[i]);
		}
	}
}

void hansen_bliek(const IntervalMatrix& A, const IntervalVector& B, IntervalVector& x) {
	int n=A.nb_rows();
	assert(n == A.nb_cols()); // throw NotSquareMatrixException();
	assert(n == x.size() && n == B.size());

	Matrix Id(n,n);
	for (int i=0; i<n; i++)
		for (int j=0; j<n; j++) {
			Id[i][j]   = i==j;
		}

	IntervalMatrix radA(A-Id);
	Matrix InfA(Id-abs(radA.lb()));
	Matrix M(n,n);

	real_inverse(InfA, M);  // may throw SingularMatrixException...

	for (int i=0; i<n; i++)
		for (int j=0; j<n; j++)
			if (! (M[i][j]>=0.0)) throw NotInversePositiveMatrixException();

	Vector b(B.mid());
	Vector delta = (B.ub())-b;

	Vector xstar = M * (abs(b)+delta);

	double xtildek, xutildek, nuk, max, min;

	for (int k=0; k<n; k++) {
		xtildek = (b[k]>=0) ? xstar[k] : xstar[k] + 2*M[k][k]*b[k];
		xutildek = (b[k]<=0) ? -xstar[k] : -xstar[k] + 2*M[k][k]*b[k];

		nuk = 1/(2*M[k][k]-1);
		max = nuk*xutildek;
		if (max < 0) max = 0;

		min = nuk*xtildek;
		if (min > 0) min = 0;

		/* compute bounds of x(k) */
		if (xtildek >= max) {
			if (xutildek <= min) x[k] = Interval(xutildek,xtildek);
			else x[k] = Interval(max,xtildek);
		} else {
			if (xutildek <= min) x[k] = Interval(xutildek,min);
			else { x.set_empty(); return; }
		}
	}
}

int Interset(IntervalVector x, IntervalVector y) {
assert(x.size() == y.size());
int res = 2;
int i;
for(i = 0;i<x.size();i++) {
	if (!x[i].is_subset(y[i])) {
		res = 1;
		break;}
	}
while(i<x.size()) {
    if(!x.intersects(y)) {
		res = 0;
		return res;
	}
	i++;}
return res;}

/*
 * Restricted propagation procedure. Used when a direction is solved, but no new q-intersection is found.
 * Same as propagate, except that upper bounds are not recorded (because we did not find a curr_qinter to work with).
 */
void propagate_no_ub(const Array<IntervalVector>& boxes, IntStack ***dirboxes, int dimension, bool left, IntervalVector& hull_qinter, vector<BitSet *>& nogoods) {
	
	unsigned int n = hull_qinter.size();
	unsigned int p = boxes.size();
	
	IntervalVector b(n);
	
	/* Create the new nogood */
	BitSet *newNogood = new BitSet(0,p-1,BitSet::empt);
	
	/* We iterate through the boxes to propagate bounds */
	/* This does not seem to be optimal ; maybe we should study directly the directions' lists ? */
	for (int i=0; i<p; i++) {
		b = boxes[i];
		
		/* Check if the box can be added to the new nogood */
		if ((left && (b[dimension].lb() < hull_qinter[dimension].lb())) || (!left && (b[dimension].ub() > hull_qinter[dimension].ub()))) {
			newNogood->add(i);	
		}
		
		/* Check if the box is strictly outside our current q-inter hull */
		if ((left && (b[dimension].ub() < hull_qinter[dimension].lb())) || (!left && (b[dimension].lb() > hull_qinter[dimension].ub()))) {
			for (int k=dimension+1; k<n; k++) {
				if (dirboxes[k][0]->contain(i)) dirboxes[k][0]->remove(i);
				if (dirboxes[k][1]->contain(i)) dirboxes[k][1]->remove(i);
			}
			
			continue;
		}
	}
	
	nogoods.push_back(newNogood);
}

void Function::hansen_matrix(const IntervalVector& box, IntervalMatrix& H) const {
	int n=nb_var();
	int m=expr().dim.vec_size();

	assert(H.nb_cols()==n);
	assert(box.size()==n);
	assert(expr().dim.is_vector());
	assert(H.nb_rows()==m);

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

	// test!
//	int tab[box.size()];
//	box.sort_indices(false,tab);
//	int var;

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

}

IntervalVector Set::inflate_one_float(const IntervalVector& box) {
	IntervalVector inflated=box;
	for (int i=0;i<box.size();i++) {
		double lb=box[i].lb();
		double ub=box[i].ub();
		inflated[i]=Interval(lb==NEG_INFINITY? lb : previous_float(lb), ub==POS_INFINITY? ub : next_float(ub));
	}
	return inflated;
}