#define M 6						// Number of warehouses
#define N 30						// Number of supermarkets

double c[M][N];					// Objective function coefficients
double fc[M];						// Fixed cost of running warehouse i
double kap[M];						// Capacity of warehouse i
double d[N];						// Demand of supermarket j

continuous X[M][N];					// The amount of j’s demand met by i
continuous Y[M];					// Fuzzy boolean for warehouse i
continuous Cost;					// The objective function

data()
{

fc = { 10500, 18500, 30500, 17500, 24500 };

kap = { 19000, 27000, 29000, 39000, 23000, 44000 };

d = {
2146,
187,
672,
133 ,
2131,
559,
2370,
1089,
2133,
2132,
2495,
904,
1466,
2143,
3615,
564,
226,
3016,
253,
2195,
2138,
807,
551,
304,
814,
3337,
4368,
577,
482,
495
};

c = {
6739,
3204,
4914,
32372,
1715,
6421,
81972,
33391,
2020,
1459,
141015,
17684,
38207,
1953,
17181,
25640,
7031,
78453,
9452,
8597,
1581,
23170,
12087,
4883,
24063,
4124,
281463,
11056,
8585,
12480,


3727,
4673,
13451,
372672,
9745,
12055,
97602,
60774,
54470,
7146,
38011,
39723,
16111,
16981,
168663,
57109,
15576,
2542,
34056,
7095,
10355,
5457,
26409,
29982,
2152,
23701,
28499,
26544,
2480,
1995,

205925,
32069,
42477,
5044,
36054,
35602,
10492,
92515,
12441,
14113,
2030,
48702,
19877,
12851,
39682,
12148,
406770,
22113,
22449,
25455,
11116,
13346,
35106,
229188,
18070,
18181,
73603,
63568,
65177,
8618,

70728,
52917,
20714,
32575,
210766,
66703,
18481,
3928,
34221,
11999,
7650,
3845,
19622,
21024,
1577,
16197,
43134,
6370,
1869,
1402,
104130,
15322,
15319,
4089,
25399,
25154,
6305,
36644,
7754,
10500,

1326,
36072,
9670,
10822,
24603,
8180,
325852,
11424,
14122,
22151,
8229,
7880,
25927,
203364,
12049,
11400,
59561,
27330,
52117,
6428,
39587,
32225,
15620,
23312,
169251,
53124,
14368,
3020,
24448,
7886,

5219,
2396,
13876,
29681,
1061,
10383,
65767,
16770,
1324,
869,
12638,
8429,
15832,
3428,
16297,
15763,
2542,
27445,
3542,
7254,
693,
26166,
3801,
8930,
11050,
5611,
253234,
5582,
7458,
19069
};

}

setup()
{
      and(int j=0;j<N;j++)							// Supermarket constraints
            sigma(int i=0;i<M;i++) X[i][j] == 1.0;

	and(int i=0;i<M;i++)							// Tight warehouse constraints
	  sigma(int j=0;j<N;j++) d[j]*X[i][j] <= kap[i]*Y[i];

	Cost ==
		sigma(int i=0;i<M;i++) fc[i]*Y[i]				// Fixed costs
		+
		sigma(int i=0;i<M;i++)						// Proportional costs
			sigma(int j=0;j<N;j++) c[i][j]*X[i][j];
}

quick_and_dirty()		// Compute and enforce bluff
{
double scratch;

	min: Cost;
	// Add up proportional costs
	scratch = sigma(int i=0;i<M;i++)
				sigma(int j=0;j<N;j++) c[i][j]*wp(X[i][j]);

	// Add full fixed costs for any partially open warehouse
	and(int i=0;i<M;i++) 
		scratch = scratch + fc[i]*ceil(Y[i]);

	Cost <= scratch;			// A bluff

	printf("Quick and dirty solution is %f\n",scratch);
}

fix(continuous Y, X[])		// Called as fix(Y[i],X[i])
{
	either 		
		Y == 1.0;				// Open the warehouse
	or {
		Y == 0.0;				// Do not open the warehouse
		and(int j=0;j<N;j++)
			X[j] == 0.0;			// Good programming practice
	}
}

add_cuts()
{
int flag;

	flag = 1;					// Turn flag on
	and(;flag != 0;){			// Loop while flag is on
		flag = 0;					// Turn flag off
		and(int i=0;i<M;i++)
			and(int j=0;j<N;j++)
				if wp(X[i][j]) > wp(Y[i]);		// Violated cut?
				then {
					X[i][j] <= Y[i];				// Add cut
					flag = 1;					// Turn flag on
					min: Cost;				// Reoptimize
				}
	}
}

search()
{
int bv;

	and(;;) {
		min: Cost;					// Optimize linear relaxation
		add_cuts();						// Add the cuts
		if determine_branching_variable(bv);
		then fix(Y[bv],X[bv]);				// Repair injured variable
		else break;						// Exit the loop
	}
}

determine_branching_variable(int & bv)	
// Selects injured warehouse with greatest number of users
{
int cnt, xijcnt;

	cnt = -1;
	bv = M;					// Unchanged if all Y[i] integral
	and(int i=0;i<M;i++)
		if not integral(Y[i]);
		then {
			xijcnt = 0;
			and(int k=0;k<N;k++)				// Count users
				if wp(X[i][k]) > 0;
				then xijcnt = xijcnt+1;
			if xijcnt > cnt;						// Select i with largest number of users
			then {
				bv=i;
				cnt=xijcnt;
			}
		}
	bv < M;
}

2lp_main()
{
	data();
	setup();
	quick_and_dirty();			// Compute and enforce bluff

	find_min: Cost; 
	subject_to 
		search();

	and(int i=0;i<M;i++)
		if wp(Y[i]) == 1;
		then printf("\nUse warehouse %d",i);
}
 