#define TOURS 65  
#define FLIGHTS 32
#define MAXNUM 7
#define N 10					// > (log base 2 (TOURS))  
#define LEFT 0
#define RIGHT 1

continuous Z,Tour[TOURS];

double cost[TOURS];			// Cost of tour

int card[FLIGHTS];				// Cardinality of set of tours for flight
int set[FLIGHTS][MAXNUM];
int trail[3*FLIGHTS]; 			// > (log base 2 (MAXNUM))* FLIGHTS
int top;						// Top of trail stack
int cnt; 						// Node counter
int itinerary[FLIGHTS][N];		// The vector itinerary[f] is 
							// stack for flight f
int depth[FLIGHTS];				// depth[f]-1 is top of stack for f


data()
{
#define NA -1
set ={34,4,50,40,22,33,61,
	8,62,21,19,32,34,1,
	51,3,63,13,25,35,2,
	31,5,14,29,36,6,NA,
	10,41,27,37,5,NA,NA,
	11,56,20,23,38,4,NA,
	12,14,52,42,39,32,NA,
	18,1,43,28,40,9,NA,
	21,13,44,28,41,6,NA,
	12,4,62,53,45,42,7,
	54,55,21,19,46,43,38,
	21,3,26,13,25,44,41,
	4,5,56,29,45,42,NA,
	47,55,27,46,44,NA,NA,
	11,31,20,23,47,41,NA,
	57,54,7,22,48,49,NA,
	48,32,16,28,49,6,NA,
	17,33,53,50,40,NA,NA,
	19,30,26,59,51,33,NA,
	15,44,16,49,52,46,NA,
	17,22,26,60,53,12,NA,
	58,14,64,8,2,54,36,
	51,11,21,9,2,55,50,
	27,23,26,3,5,56,37,
	14,15,4,9,57,34,NA,
	0,8,17,58,25,NA,NA,
	1,19,10,13,59,0,NA, 
	2,24,7,12,60,39,NA,
	8,0,16,18,61,22,NA,
	7,13,6,28,62,45,NA,
	2,24,17,8,12,63,50,
	28,31,26,9,49,64,60};
#undef NA

cost =	{7,3,7,2,5,4,6,4,6,2,2,3,4,
		5,3,3,3,2,3,3,3,2,2,2,3,8,
		8,5,6,9,10,11,7,2,2,6,3,3,6,
		3,6,7,2,7,2,3,3,6,3,2,4,3,
		3,2,2,3,2,4,3,3,3,6,2,2,4};

card = 	{7,7,7,6,5,6,6,6,6,7,7,
		7,6,5,6,6,6,5,6,6,6,7,
		7,7,6,5,6,6,6,6,7,7};
}

setup()
{
	and(int i=0;i<TOURS;i++) 
		Tour[i] <= 1; 				// Each tour can be used at most once
	
	and(int flt=0;flt<FLIGHTS;flt++)									// Cover each flight 
	  sigma(int i=0;i<card[flt];i++)Tour[set[flt][i]] == 1;

	// Objective function
	Z == sigma(int k=0;k<TOURS;k++) cost[k]*Tour[k]; 
}

injured_flight(int & be)
{
int nots;			// Number of tours
	nots = MAXNUM+1;
	be = FLIGHTS;
	and(int flt=0;flt<FLIGHTS;flt++)	// Loop through flights
		if {
			and(int i=0;i<card[flt];i++)		// Injured ?
				wp(Tour[set[flt][i]]) < 1.0;
			card[flt] < nots;				// Worse than previous ?
		}
		then {						// Worst injury found so far
			be = flt;						// This is element for branching
			nots = card[flt];				// Update number of tours
		}
	be < FLIGHTS;					// Injury found ?
}

compute_interval(int flt, int & m, & mid, & n)
{
	m = 0; n = card[flt]-1;			// Original interval
	and(int i=0;i<depth[flt];i++)		// Binary path
		if itinerary[flt][i] == LEFT; 		// of branches
		then n = m + (n - m)/2;			// left and right
		else m = m + (n - m)/2+1;

 	mid = m+(n-m)/2; 				// New midpoint
	depth[flt] = depth[flt]+1;			// Increment depth
}

tend_to_flight(int flt,k)
{
int m,mid,n;				// Storeback variables

	compute_interval(flt,m,mid,n);									// Find m,mid,n
	top = k;					// Increment top of stack
	trail[top] = flt;				// Push flt onto trail stack
	// Bias search towards heavier half
	if wp(sigma(int i=m;i<=mid;i++) Tour[set[flt][i]]) >
	    wp(sigma(int i=mid+1;i<=n;i++) Tour[set[flt][i]]);
	then left_then_right(flt,m,mid,n,k);
	else right_then_left(flt,m,mid,n,k);
}

left_then_right(int flt,m,mid,n,k)
{
	either left(flt,m,mid,n); 	// top is equal to k
	or {
		if top > k;			// Deep backtracking has occurred
		then {  				// Untrail required
			and(int t=top;t>k;t--)
				depth[trail[t]] = depth[trail[t]] - 1;
			top = k;			// Reset
		}
		right(flt,m,mid,n);
	}
}

left(int flt,m,mid,n)
{
	cnt = cnt+1;						// Count new node
	itinerary[flt][depth[flt]-1] = LEFT;
	shrink_sets(flt,m,mid,n,LEFT);		// Branch left
}

right(int flt,m,mid,n)
{
	cnt = cnt+1;						// Count new node
	itinerary[flt][depth[flt]-1] = RIGHT;
	shrink_sets(flt,m,mid,n,RIGHT);		// Branch right
}

shrink_sets(int flt,a,b,c,dir)
{
	if dir == LEFT;
	then and(int i=b+1;i<=c;i++)
		Tour[set[flt][i]] == 0; 			// Annihilate right
	else 	and(int i=a;i<=b;i++)
		Tour[set[flt][i]] == 0;				// Annihilate left
	min: Z;					// Re-optimize for injury method
}

right_then_left(int r,m,mid,n,int k)
{
	either right(r,m,mid,n);
	or {
		and(int t=top;t>k;t--)					// Untrail if needed
			depth[trail[t]] = depth[trail[t]] - 1;
		top = k;							// Re-enforce if needed
		left(r,m,mid,n);
	}
}

optimize()
{
int flight;

	min: Z;
	find_min: Z; 
	subject_to
		and(int k=0;;k++){				// k records depth of search
			absgap(1.0); 				// Each solution strictly better
			if injured_flight(flight); 
			then tend_to_flight(flight,k);
			else break;				// Eureka effect
		}
	output();
}

output()
{
	printf("There will be %.0f tours needed\n",
		sigma(int k=0;k<TOURS;k++) Tour[k]);
	printf("At a cost of %f\n",Z);
	printf("There were %d choice nodes generated\n",cnt);
}

2lp_main()
{
int flight;
 
	data();
	setup();  
	top = -1;				// Empty stack
 	optimize();
}

 