[FlopCpp] Question on indexing constraints

[Michal Kaut]

Hello,

I forgot to attach the example file to my previous email, so I am 
sending it again - so please ignore the previous one. I apologize for 
any inconvenience...


I have been struggling to understand the indexing of FlopC++ constraints
and must admit I still do not understand it properly. What confuses me
most is that for each type of constraints I wanted to test, I was able
to find several "seemingly equivalent" formulations, all of which
compile, but only one of which gives the expected results.

I have tried to demonstrate it on the attached example, where the
constraints with a name differing only in number are supposed to do the
same thing, but only those that have "WORKS" in the comment actually
generate the expected result. (I should probably confess that I am an
AMPL user, so I have probably different understanding what an "index"
means, compared to the GAMS-influenced designers of FlopC++.)
- Could someone try to explain to me why the bad formulations do not
work, i.e. what they actually do? (They do something, as they compile.)
- Note that the example is only to illustrate indexing, the resulting
model is not supposed to make any sense.


In addition, I was rather disappointed by the compilation speed (i.e.
the lack of it): this very small example takes over a minute to compile
on mu 3GHz Pentium 4 machine with gcc 3.4.5 (mingw on WinXP).
Of this time, 39 seconds are spend on the last three sets of constraints
(getIndx1-3), i.e. the compile+link time is 26 seconds without those
constraints. Is this normal/expected behaviour?
(Another interesting thing is that adding each of the three sets of
constraints individually adds 5s, 10s and 10s respectively, i.e. 25
seconds in total, while adding all three adds 39 seconds to the
compilation time...)


And one more question: is it possible to get rid of the "eval subsetref"
messages? Even with this small example, I get hundreds of them and
setting model.silent() does not switch them off...


Thanks a lot in advance.

Regards,
Michal Kaut

-------------- next part --------------
// $Id$
#include "flopc.hpp"
using namespace flopc;
#include "OsiClpSolverInterface.hpp"

int main() {

	int i;

	MP_model &model	= MP_model::getDefaultModel();
	model.setSolver(new OsiClpSolverInterface);
	model.verbose();
//	model.silent();

	enum {a1, a2, a3, numA};
	enum {b1, b2, numB};
	enum {c1, c2, numC};

	MP_set A(numA);
	MP_set B(numB);
	MP_set C(numC);
	MP_subset<2> AB(A,B);
		AB.insert(a1, b1);
		AB.insert(a1, b2);
		AB.insert(a2, b1);
		AB.insert(a3, b2);
	MP_subset<1> As(A);
		As.insert(a1);
		As.insert(a3);

	MP_index a, b, c, ab;

	MP_data dA(A); // dA = {2,4,6}
		for (i=0; i<numA; i++)
			dA(i) = 2*(i+1);
	MP_data dB(B); // dB = {1,2}
		for (i=0; i<numB; i++)
			dB(i) = i+1;
	MP_data dAB(AB); // dAB = {{3,4},{5,.},{.,8}}
	MP_data dAB2(AB);
	MP_data dAB3(AB);
	MP_data dAB4(AB);
		// Using 2-dim set name and component set names -> WORKS
		dAB(AB(A,B)) = dA(A) + dB(B);
		// Using 2-dim set name and component indices -> WORKS
		dAB2(AB(a,b)) = dA(a) + dB(b);
		// Using component set names only
		dAB3(A,B) = dA(A) + dB(B);
		// Using only indices
		dAB4(a,b) = dA(a) + dB(b);
	MP_data dC(C); // dC = {5,6}
		for (i=0; i<numC; i++)
			dC(i) = i+5;
	MP_data dBC(B,C); // dBC = {{5,6},{10,12}}
	MP_data dBC2(B,C);
	MP_data dBC3(B,C);
		// Using set names only -> WORKS
		forall(B*C, dBC(B,C) = dB(B) * dC(C));
		// Using set names and indices
		forall(B(b)*C(c), dBC(b,c) = dB(b) * dC(c));
		// Using set names and indices, different syntax
		forall((B(b),C(c)), dBC(b,c) = dB(b) * dC(c));

	MP_variable xA(A);
	MP_variable xAB(AB);
	MP_variable xAC(A,C);
	MP_variable xABC(AB,C);
	// AMPL: for a given 'a' in 'A', want "sum{b in B : (a,b) in AB} x[a,b]"
	MP_constraint cA0(A); cA0.setName("cA0");
	MP_constraint cA1(A); cA1.setName("cA1");
	MP_constraint cA2(A); cA2.setName("cA2");
		// Using the two-dim set name -> WORKS
		cA0(A)  =  sum(AB(A,B), dAB(AB) * xAB(AB)) <= 1;
		// Using set names
		cA1(A)  =  sum(AB(A,B), dAB(A,B) * xAB(A,B)) <= 1;
		// Using indices
		cA2(a)  =  sum(AB(a,b), dAB(a,b) * xAB(a,b)) <= 1;
	//
	MP_constraint cAB1(AB); cAB1.setName("cAB1");
	MP_constraint cAB2(AB); cAB2.setName("cAB2");
	MP_constraint cAB3(AB); cAB3.setName("cAB3");
		// Using only set names -> WORKS
		cAB1(AB)  =  xAB(AB) <= sum(C, dC(C) * xABC(AB,C));
		// Using only indices
		cAB2(a,b)  =  xAB(a,b) <= sum(C(c), dC(c) * xABC(a,b,c));
		// Using mixture of set names and indices
		cAB3(a,b)  =  xAB(a,b) <= sum(C, dC(C) * xABC(a,b,C));
	//
	MP_constraint cAC1(A,C); cAC1.setName("cAC1");
	MP_constraint cAC2(A,C); cAC2.setName("cAC2");
		// Using set names -> WORKS
		cAC1(A,C)  =  (dA(A) + dC(C)) * xAC(A,C) <= xA(A);
		// Using indices -> Works in this case, not in another one
		cAC2(a,c)  =  (dA(a) + dC(c)) * xAC(a,c) <= xA(a);
	//
	MP_constraint cABC0(AB,C); cABC0.setName("cABC0");
	MP_constraint cABC1(AB,C); cABC1.setName("cABC1");
	MP_constraint cABC2(AB,C); cABC2.setName("cABC2");
	MP_constraint cABC3(AB,C); cABC3.setName("cABC3");
		// Using only set names, with AB(A,B) sub-indices -> WORKS
		cABC0(AB(A,B),C)  =  xAC(A,C) + xABC(AB,C) == dBC(B,C);
		// Using only set names
		cABC1(AB,C)  =  xAC(A,C) + xABC(AB,C) == dBC(B,C);
		// Using indices
		cABC2(a,b,c)  =  xAC(a,c) + xABC(a,b,c) == dBC(b,c);
		// Using a mixture of set names and indices
		cABC3(AB(a,b),C)  =  xAC(a,C) + xABC(a,b,C) == dBC(B,C);
	//
	MP_constraint cAs1(As); cAs1.setName("cAs1");
	MP_constraint cAs2(As); cAs2.setName("cAs2");
	MP_constraint cAs3(As); cAs3.setName("cAs3");
		// Using the superset name for indexing -> WORKS
		cAs1(As(A))  =  xA(A) <= dA(A);
		// Using the (sub)set name for indexing
		cAs2(As)  =  xA(As) <= dA(As);
		// Using the (sub)set + index for indexing -> WORKS
		cAs3(As(a))  =  xA(a) <= dA(a);
	//
	// These constraints use all variables -> get the variable indices...
	MP_constraint getIndx1; getIndx1.setName("getIndx1");
	MP_constraint getIndx2; getIndx2.setName("getIndx2");
	MP_constraint getIndx3; getIndx3.setName("getIndx3");
		// Using only sets (no indices) - WORKS
		getIndx1  =  sum(A, dA(A) * xA(A)) - sum(AB, dAB(AB) * xAB(AB))
								 + sum(A*C, xAC(A,C))
								 - sum(AB*C, dAB(AB) * dC(C) * xABC(AB,C)) <= 100;
		// Using indices; single index for set AB
		getIndx2  =  sum(A(a), dA(a) * xA(a)) - sum(AB(ab), dAB(ab) * xAB(ab))
								 + sum(A(a)*C(c), xAC(a,c))
								 - sum(AB(ab)*C(c), dAB(ab) * dC(c) * xABC(ab,c)) <= 100;
		// Using indices; double index for set AB
		getIndx3  =  sum(A(a), dA(a) * xA(a)) - sum(AB(a,b), dAB(a,b) * xAB(a,b))
								 + sum(A(a)*C(c), xAC(a,c))
								 - sum(AB(a,b)*C(c), dAB(a,b) * dC(c) * xABC(a,b,c)) <= 100;

	MP_expression ObjF = sum(AB, dAB(AB) * xAB(AB));
	model.setObjective(ObjF);
	model.attach();
	model->writeMps("sup_boat");


	/**********************************************************************/
	cout << "Set AB has " << AB.size() << " elements." << endl;
	dA.display("Data dA:");
	dB.display("Data dB:");
	dAB.display("Data dAB:");
	dAB2.display("Data dAB2:");
	dAB3.display("Data dAB3 [WRONG]:");
	dAB4.display("Data dAB4 [WRONG]:");
	dBC.display("Data dBC:");
	dBC2.display("Data dBC2 [WRONG]:");
	dBC3.display("Data dBC3 [WRONG]:");
	/**********************************************************************/

	model.minimize();

	return 0;
}