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

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 InHC4Revise::ibwd(const Function& f, const Domain& y, IntervalVector& x, const IntervalVector& xin) {

	Eval e;

	if (!xin.is_empty()) {
		e.eval(f,xin);

		assert(!f.expr().deco.d->is_empty());

		for (int i=0; i<f.nb_nodes(); i++)
			*f.node(i).deco.p = *f.node(i).deco.d;
	}
	else {
		for (int i=0; i<f.nb_nodes(); i++)
			f.node(i).deco.p->set_empty();
	}

	e.eval(f,x);

	assert(!f.expr().deco.d->is_empty());

	*f.expr().deco.d = y;

	try {

		f.backward<InHC4Revise>(*this);

		f.read_arg_domains(x);

	} catch(EmptyBoxException&) {
		x.set_empty();
	}
}

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 InHC4Revise::iproj(const Function& f, const Domain& y, IntervalVector& x) {

	for (int i=0; i<f.nb_nodes(); i++)
		f.node(i).deco.p->set_empty();

	Eval().eval(f,x);

	*f.expr().deco.d = y;

	try {
		f.backward<InHC4Revise>(*this);

		if (f.all_args_scalar()) {
			int j;
				for (int i=0; i<f.nb_used_vars; i++) {
					j=f.used_var[i];
					x[j]=f.arg_domains[j].i();
				}

			}
		else
			load(x,f.arg_domains,f.nb_used_vars,f.used_var);

	} catch(EmptyBoxException&) {
		x.set_empty();
	}
}

void CtcPolytopeHull::contract(IntervalVector& box) {

	if (!(limit_diam_box.contains(box.max_diam()))) return;
	// is it necessary?  YES (BNE) Soplex can give false infeasible results with large numbers
	//       	cout << " box before LR " << box << endl;


	try {
		// Update the bounds the variables
		mylinearsolver->initBoundVar(box);

		//returns the number of constraints in the linearized system
		int cont = lr.linearization(box,mylinearsolver);

		if(cont<1)  return;
		optimizer(box);

		//	mylinearsolver->writeFile("LP.lp");
		//		system ("cat LP.lp");
		//		cout << " box after  LR " << box << endl;
		mylinearsolver->cleanConst();


	}
	catch(EmptyBoxException&) {
		box.set_empty(); // empty the box before exiting in case of EmptyBoxException
		mylinearsolver->cleanConst();
		throw EmptyBoxException();
	}

}

//-------------------------------------------------------------------------------------------------------------
void CtcPixelMap::grid_to_world(IntervalVector& box) {
    for(unsigned int i = 0; i < I.ndim; i++) {
        box[i] &= Interval(pixel_coords[2*i], pixel_coords[2*i+1]+1) * I.leaf_size_[i] + I.origin_[i];
        if(box[i].is_empty()){
            box.set_empty();
            return;
        }   
    }
}

void CtcInteger::contract(IntervalVector& box) {
	for (int i=0; i<nb_var; i++) {
		if (!is_int[i]) continue;
		proj_integer(box[i]);
		if (box[i].is_empty()) {
			box.set_empty();
			throw EmptyBoxException();
		}
	}
}

void CtcIn::contract(IntervalVector& box) {

	// it's simpler here to use direct computation, but
	// we could also have used CtcFwdBwd

	try {
		HC4Revise().proj(_f,_d,box);
	} catch (EmptyBoxException& e) {
		box.set_empty();
		throw e;
	}
}

void CtcSegment::contract(IntervalVector &box) {

	if (nb_var==6) {
		ctc_f->contract(box);
		if (box.is_empty()) return;
		ctc_g->contract(box);
	}
	else {
		X_with_params[0] = box[0];
		X_with_params[1] = box[1];
		//        X_with_params=cart_prod(box,)

		ctc_f->contract(X_with_params);
		if (X_with_params.is_empty()) { box.set_empty(); return; }

		ctc_g->contract(X_with_params);
		if (X_with_params.is_empty()) { box.set_empty(); return; }

		//        box = X_with_params.subvector(0,1);
		box[0] = X_with_params[0];
		box[1] = X_with_params[1];
	}
}

void CtcFirstOrderTest::contract(IntervalVector& box, ContractContext& context) {
	if(box.size() == 2) {
		return;
	}
	BxpNodeData* node_data = (BxpNodeData*) context.prop[BxpNodeData::id];
	if(node_data == nullptr) {
		ibex_error("CtcFirstOrderTest: BxpNodeData must be set");
	}
	vector<IntervalVector> gradients;
	for (int i = 0; i < system_.normal_constraints_.size() - 1; ++i) {
		if (!system_.normal_constraints_[i].isSatisfied(box)) {
			gradients.push_back(system_.normal_constraints_[i].gradient(box));
		}
	}
	for (int i = 0; i < system_.sic_constraints_.size(); ++i) {
		if (!system_.sic_constraints_[i].isSatisfied(box, node_data->sic_constraints_caches[i])) {
			gradients.push_back(system_.sic_constraints_[i].gradient(box, node_data->sic_constraints_caches[i]));
		}
	}

	// Without the goal variable
	IntervalMatrix matrix(nb_var - 1, gradients.size() + 1);
	matrix.set_col(0, system_.goal_function_->gradient(box.subvector(0, nb_var - 2)));
	for (int i = 0; i < gradients.size(); ++i) {
		matrix.set_col(i + 1, gradients[i].subvector(0, nb_var - 2));
	}
	bool testfailed = true;
	if (matrix.nb_cols() == 1) {
		if (matrix.col(0).contains(Vector::zeros(nb_var - 2))) {
			testfailed = false;
		}
	} else {
		int* pr = new int[matrix.nb_rows()];
		int* pc = new int[matrix.nb_cols()];
		IntervalMatrix LU(matrix.nb_rows(), matrix.nb_cols());
		testfailed = true;
		try {
			interval_LU(matrix, LU, pr, pc);
		} catch(SingularMatrixException&) {
			testfailed = false;
		}
		delete[] pr;
		delete[] pc;
	}
	if(testfailed) {
		box.set_empty();
	}

}

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; }
		}
	}
}

void CtcEmpty::contract(IntervalVector& box) {
	BoolInterval t= pdc.test(box);
	switch (t) {
	case YES:
		box.set_empty();
		set_flag(FIXPOINT);
		break;
	case NO:
		// The constraint is inactive only if
		// the test is monotonous w.r.t. inclusion
		//	set_flag(INACTIVE);
		break;
	default: break;
	}
}

void CtcNotIn::contract(IntervalVector& box) {

	// it's simpler here to use direct computation, but
	// we could also have used CtCunion of two CtcFwdBwd

	IntervalVector savebox(box);
	try {
		HC4Revise().proj(f,d1,box);
	} catch (EmptyBoxException& ) {box.set_empty(); }
	try {
		HC4Revise().proj(f,d2,savebox);
	} catch (EmptyBoxException& ) {savebox.set_empty(); }

	box |= savebox;
	if (box.is_empty()) throw EmptyBoxException();

}

bool HC4Revise::proj(const Domain& y, IntervalVector& x) {
	eval.eval(x);
	//std::cout << "forward:" << std::endl; f.cf.print(d);

	bool is_inner=false;

	try {
		is_inner = backward(y);

		d.read_arg_domains(x);

		return is_inner;

	} catch(EmptyBoxException&) {
		x.set_empty();
		return false;
	}
}

void Gradient::gradient(const Array<Domain>& d2, IntervalVector& gbox) {
	assert(f.expr().dim.is_scalar());

	_eval.eval(d2);

	// outside definition domain -> empty gradient
	if (d.top->is_empty()) { gbox.set_empty(); return; }

	gbox.clear();

	g.write_arg_domains(gbox);

	f.forward<Gradient>(*this);

	g.top->i()=1.0;

	f.backward<Gradient>(*this);

	g.read_arg_domains(gbox);
}