#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 plant i
continuous Cost;					// The objective function

int bv;							// Global utility variable for injury method
double best_so_far;				// Global utility variable for bound

// Definitions and declarations from A* search
#define NIL -99						
#define QUEUESIZE 400						 
int n_lnk[QUEUESIZE], predecessor[QUEUESIZE];
double f[QUEUESIZE];
int ste,end;

// The next three arrays hold state information
int next_decision[QUEUESIZE];						// The integer k
int decision_list[QUEUESIZE][M];						// bit0,...,bitk-1
int permutation[QUEUESIZE][M];						// i0,...,ik-1


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++)							// Demand must be met
            sigma(int i=0;i<M;i++) X[i][j] == 1.0;

	and(int i=0;i<M;i++)							// Capacity can’t be exceeded
	  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();				
		determine(bv);					// Branching variable
		if bv == M;
		then break;					// Exit the loop
		fix(Y[bv],X[bv]);
	}
}

determine(int & bv)
{
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
			then {						// number of users
				bv=i;
				cnt=xijcnt;
			}
		}
}

fix_common_variable_bindings(int ste)
{
int perm;
	and(int k=0;k<next_decision[ste];k++) {
		perm = permutation[ste][k];
		if wp(Y[perm]) == decision_list[ste][k];
		then Y[perm] == decision_list[ste][k];
	}
}

rebuild_node(int ste)
{
	and(int k=0;k<next_decision[ste];k++) 
		Y[permutation[ste][k]] == decision_list[ste][k];
	min: Cost; 					// Optimize for injury method
}

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

goal_node(int ste)
{
	// Set bv = branching variable or fail if eureka effect
	not determine_branching_variable(bv);			// As in Model 13.3
}

freeze_node(int ste)
{
	next_decision[ste] = M;						// All decisions made
	and(int j=0;j<M;j++) {
		decision_list[ste][j] = wp(Y[j]);
		permutation[ste][j] = j;					// Identity permutation
	}
	best_so_far = f[ste];						// Best known bound
}

land_and_doig_try_moves(int ste,int & end) 
{
	if land_and_doig_add_move(ste,end,1); 
	then end = end + 1;
	if land_and_doig_add_move(ste,end,0); 
	then end = end + 1;
}

land_and_doig_add_move(int ste,newstate,bit) 
{
	not not {			// Double negation to undo constraints
		Y[bv] == bit;
		min: Cost;						// Evaluate node
		f[newstate] = wp(Cost);				// Record evaluation
		f[newstate] < best_so_far;			// Check bound
		update_state(ste,newstate,bit);		// Record state
	}
}

update_state(int ste, newstate, bit)
{
int nd;
	nd = next_decision[ste];
	// Copy parent state
	and(int i=0;i<nd;i++){			
		permutation[newstate][i] = permutation[ste][i];
		decision_list[newstate][i] = decision_list[ste][i];
	}
	// Add new elements
	permutation[newstate][nd] = bv;								
	decision_list[newstate][nd] = bit;
	next_decision[newstate] = next_decision[ste]+1;
	if next_decision[newstate] == M;			// Goal node case
	then best_so_far = f[newstate];
}

init_queue_vars()
{
	ste = 0;				// First state to examine is root node
	end = 1;				// First empty slot on priority queue
	n_lnk[0] = NIL;			// First element has no ancestor
}

insert_states(int i,ste,end) 
{
	// end-i states to be added to priority queue
	and(int j=i;j<end;j++) {
		if n_lnk[ste] == NIL; 
		then splice(j,ste,NIL);
		else find_and_fill_slot(j);
	}
}

splice(int j,pre,post)
{
			n_lnk[pre] = j;
			n_lnk[j] = post;
}

find_and_fill_slot(int j)
{
int ptr,pre_ptr;

	pre_ptr = ste; 
	ptr = n_lnk[ste];				// Never insert before ste
	and(;n_lnk[pre_ptr] != NIL;){
		if f[j] < f[ptr]; 					
		then break;					// Insertion spot found
		pre_ptr = ptr;					// Move pre_ptr up to ptr
		ptr = n_lnk[ptr];				// Move ptr up to its successor
	}
	splice(j,pre_ptr,ptr);				// Insert the new state
}

initialize_bound()
{
	best_so_far = ub(Cost); 				// Bound from quick_and_dirty
}

output()			// Land and Doig version
{
	printf("ste is %d\n",ste);
	printf("f[ste] is %f\n",f[ste]);
	and(int k=0;k<M;k++) 
		if decision_list[ste][k] == 1;
		then printf("Build warehouse %d\n",permutation[ste][k]);
}

2lp_main()			// Land and Doig Search
{
	init_queue_vars();				// As for A* search
	data();						// As in Model 13.3
	setup();						// Ditto 
	quick_and_dirty();				// Ditto
	initialize_bound();				// Set initial bound

	and(int i=end;i<QUEUESIZE-2;i=end){
		if next_decision[ste] == M;						
		then {		// Optimal node at front of priority queue
				output();			// New version below
			return;			  
		}
		not not {
			fix_common_variable_bindings(ste);	// Good idea
			rebuild_node(ste);					// Node to be developed
			if goal_node(ste);					// No injuries ? 
			then freeze_node(ste);						// Eureka effect
			else {		// Generate children
				land_and_doig_try_moves(ste,end);
				insert_states(i,ste,end);					// As in A* search
				ste = n_lnk[ste];						// Resume search
			}
		ste != NIL; 			// Check queue not mined out
		}
	}
	printf("Out of space\n");
}