#define MONTHS 6		// Number of periods
#define KINDS 5			// Kinds of vegetable oils
#define VATREV 150		// Selling price per vat
#define STORCOST 5		// Storage cost per vat
#define INITSTOR 500		// Initial storage in vats
#define ENDSTOR 500		// Final storage in vats
#define STORCAP 1000		// Monthly storage capacity

#define MODES 3			// Purchase, Process or Store
#define PURC 0
#define PROC 1
#define STOR 2

#define LIMITA 250		// Limit on group A production
#define LIMITB 200		// Limit on group B production
#define A 2
#define B 3
#define CPS 3
#define LBHARD 3.0		// Lower bound on average hardness
#define UBHARD 6.0		// Upper bound on average hardness

#define CORN 0
#define PEANUT 1
#define SAFFLOWER 2
#define SESAME 3
#define SUNFLOWER 4
#define EPSILON .000001

continuous Plan[MONTHS][MODES][KINDS];
continuous Revenue;

double price[MONTHS][KINDS];			// Table 3.2
double hardness[KINDS];				// Table 3.3
double cache;
int cnt;
2lp_main()
{
	data(price,hardness);		// As in Model 3.2

	first_stor(Plan[0]);
	last_stor(Plan[MONTHS-1][STOR]);

	and(int i=0;i<MONTHS;i++){
		storage_caps(Plan[i][STOR]);
		blend_constraints(Plan[i][PROC]);
	}

	and(int i=1;i<MONTHS;i++)
		link(Plan[i-1][STOR],Plan[i][PURC],
			Plan[i][PROC],Plan[i][STOR]);

	find_revenue(Plan,price,Revenue);
	hill_climb();

	find_max: Revenue; 
	subject_to 
		smart_search();

	output(Revenue,Plan);
	printf("The node count is %d\n",cnt);
}

data (double price[][KINDS],double hardness[])
{
	price = {	115, 110, 130, 120, 110,
			115,  90, 110, 130, 130,
			 95, 100, 130, 140, 110,
			125, 120, 120, 110, 120,
			105, 110, 150, 120, 100,
			135,  80, 140, 100,  90 
	};

	hardness = { 5.0, 4.2, 2.0, 6.1, 8.8 };
}

first_stor (continuous FirstMonth[MODES][KINDS])
{
	and(int j=0;j<KINDS;j++)
		FirstMonth[STOR][j] ==
				INITSTOR - FirstMonth[PROC][j] + FirstMonth[PURC][j];
}

last_stor (continuous LastMonth[KINDS])
{
	and(int j=0;j<KINDS;j++)
		LastMonth[j] == ENDSTOR;
}

storage_caps (continuous Stored[KINDS])
{
	and(int j=0;j<KINDS;j++)
		Stored[j] <= STORCAP;
}

blend_constraints (continuous Oils[KINDS])
{
	sigma(int j=0;j<=A;j++) Oils[j] <= LIMITA;
	sigma(int j=B;j<KINDS;j++) Oils[j] <= LIMITB;
	sigma(int i=0;i<KINDS;i++) hardness[i]*Oils[i] 
			>= LBHARD*sigma(int j=0;j<KINDS;j++) Oils[j];
	sigma(int i=0;i<KINDS;i++) hardness[i]*Oils[i]
			<= UBHARD*sigma(int j=0;j<KINDS;j++) Oils[j];
}

link (continuous OldSt[KINDS],Purc[KINDS],Proc[KINDS],NewSt[KINDS])
{
	and(int k=0;k<KINDS;k++)
		innerlink(OldSt[k],Purc[k],Proc[k],NewSt[k]);
}

innerlink (continuous OldSt,Purc,Proc,NewSt)
{
	OldSt + Purc == Proc + NewSt;   
}

find_revenue (continuous Plan[MONTHS][MODES][KINDS],
		double price[MONTHS][KINDS], continuous Revenue)
{
	Revenue ==
			sigma(int i=0;i<MONTHS;i++)
				sigma(int j=0;j<KINDS;j++)
					( VATREV*Plan[i][PROC][j] - 
			  			price[i][j]*Plan[i][PURC][j] - 
			  					STORCOST*Plan[i][STOR][j] );
}

output (continuous  Revenue, Plan[MONTHS][MODES][KINDS])
{
string months[MONTHS],modes[MODES],kinds[KINDS];													

	months = {"Jan","Feb","Mar","Apr","May","June"};
	kinds = {"Corn","Peanut","Safflower","Sesame","Sunflower"};
	modes = {"Buy","Process","Store"};

	printf("Maximum revenue is %.2f\n",Revenue);

	and(int m=0;m<MONTHS;m++) {
		printf("\nIn the month of %s\n",months[m]);
		and(int k=0;k<KINDS;k++) {
			printf(" %s oil: ",kinds[k]);
			and(int t=0;t<MODES;t++)
				printf(" %s %.0f tons\t",modes[t],Plan[m][t][k]);
			printf("\n");
		}
		if  m == 0; then
		printf("Make a sale of %.0f tons\n",
			sigma(int k=0;k<KINDS;k++) 
			 (INITSTOR+Plan[m][PURC][k]-Plan[m][STOR][k]));
		else
		printf("Make a sale of %.0f tons\n",
			sigma(int k=0;k<KINDS;k++) 
			 (Plan[m-1][STOR][k]
				+ Plan[m][PURC][k] - Plan[m][STOR][k]));
	}
}

dont_use_oil(int k, continuous Oils[])
{
	Oils[k] == 0;
	if k==CORN;
	then {		// Excluding CORN means excluding the last two
		Oils[SESAME] == 0; 
		Oils[SUNFLOWER] == 0; 
	}
}

zero(int i,j,m, continuous Oils[])
{
	if i == 0; 						// Enough zeroes accounted for
	then return;
	j >= i;						// Check room left to find enough zeroes
	if wp(Oils[m]) < EPSILON; 		// Another zero found ?
	then zero(i-1,j-1,m+1,Oils);		// A recursion
	else zero(i,j-1,m+1,Oils);			// Ditto
}

#define THRESHOLD 20
use_oil(int k,continuous Oils[])
{
	Oils[k] >= THRESHOLD;		// Use Oil[k]
	if {
		k==2;					// SAFFLOWER		
		not zero(1,2,0,Oils);			// All previous oils nonzero ?
	}
	then {					// Last two oils cannot be used
		Oils[k+1]==0;
		Oils[k+2]==0;
	}
	if {
		k==3; 					// SESAME
		not zero(2,3,0,Oils);			// Two previous oils nonzero?
	}
	then Oils[k+1]==0;			// Last oil cannot be used
}

hill_climb()
{
double cache;
	if not {
		and(int i=0;i<MONTHS;i++)
			and(int k=0;k<KINDS;k++)
				steeper_gradient(i,k,Plan[i][PROC]);
		max: Revenue;
		cache = wp(Revenue);
	}
	then {
		printf("The hillclimb failed\n"); 
	}
	else Revenue >= cache; 			// Bluff for branch-and-bound
}

steeper_gradient(int i,k,continuous Oils[])
{
double u, uu;

	if not { 
		use_oil(k,Oils); 
		max: Revenue; 
		u = wp(Revenue); 
	}
	then u = 0;

	if not { 
		dont_use_oil(k,Oils); 
		max: Revenue; 
		uu = wp(Revenue); 
	}
	then uu = 0;

	if u >= uu;
	then use_oil(k,Oils);
	else dont_use_oil(k,Oils);
}

#define EPSILON .000001
a_ok(int j, continuous Oils[]) 
{
	if j == KINDS;
	then return;

	and(int i=j;i<KINDS;i++)
		if wp(Oils[i]) > EPSILON; 						// Condition 2
		then wp(Oils[i]) >= THRESHOLD - EPSILON;

	if {
		j==0; 			// Condition 3
		wp(Oils[j]) < EPSILON;
	}
	then and(int k=CPS;k<KINDS;k++) wp(Oils[k]) < EPSILON;
	else zero(2,KINDS,0,Oils); // Condition 1
}

find_injury(int & eye, & current_month, int depth[])
{
	and(int i=current_month; ;i=(i+1)%MONTHS) {
			if not a_ok(depth[i],Plan[i][PROC]);
			then {
				eye = i; 
				current_month=(i+1)%MONTHS;		// New month
				return;
			}
		 current_month  != (i+1)%MONTHS;				// Circled round ?
	}
}

#define GAP 1.0
#define ZERO 0.0

smart_search() 
{
int top;							// Top of trail stack
int trail[MONTHS*KINDS];			// Trail stack
int depth[MONTHS];					// Depth of injury repair
int current_month;					// For circular loop
int eye;							// Storeback variable

	
	top = -1;
	and(int i=0;i<MONTHS;i++) 
		depth[i] = 0;
	current_month = 0;					
	absgap(GAP); 				// Avoid virtually duplicate solutions
	and(int m=0;;m++) {
		penalty_computation(depth,GAP);		// To tighten model
		if find_injury(eye,current_month,depth);
		then fix(m,eye,Plan[eye][PROC],depth,trail,top);	
		else break;
	}
	printf("Solution at %f\n",Revenue);
	cache = wp(Revenue);
}

fix(int m,i,continuous Oils[],int depth[],trail[],&top)
{
int j;

	j = depth[i];
	depth[i] = depth[i]+1;						// Push stacks
	top = m;									// Same as top = top+1
	trail[top] = i;
	if wp(Oils[j]) >= THRESHOLD;
	then {
		either { cnt = cnt+1; use_oil(j,Oils); }
		or {   cnt = cnt+1;
			and(int t=top;t>m;t--) 				// Clean up stacks
				depth[trail[t]] = depth[trail[t]]-1;
			top = m;
			dont_use_oil(j,Oils);
		}
	}
	else {
		either { cnt = cnt+1; dont_use_oil(j,Oils);}
		or {  cnt=cnt+1;
			and(int t=top;t>m;t--) 				// Clean up stacks
				depth[trail[t]] = depth[trail[t]]-1;
			top = m;
			use_oil(j,Oils);
		}
	 }
}

ineligible(continuous Oil)
{
	c_either lb(Oil) != 0.0; 			// Oil already used ?
	or ub(Oil) == 0.0; 				// Oil already eliminated ?
}

penalty_computation(int depth[], double gap)
{
double cwp,src,crev;
int eye[MONTHS*KINDS];
int jay[MONTHS*KINDS];
int dir[MONTHS*KINDS];
int apex;
	apex = 0;
	max: Revenue; 				// Linear optimization for reduced costs
return;
	crev = wp(Revenue);			// Record linear optimum
	and(int i=0;i<MONTHS;i++) 
		and(int j=depth[i];j<KINDS;j++) {
			if ineligible(Plan[i][PROC][j]); 
			then continue; 					// Move on if ineligible
			cwp = wp(Plan[i][PROC][j]); 			// Current value
			src =  rc(Plan[i][PROC][j]); 			// Reduced cost
			if crev + THRESHOLD*src < cache + gap; 	// Bingo									
			then {		// Record variable
				eye[apex] = i;
				jay[apex] = j;
				dir[apex] = 0;						// To fix at lower bound
				apex = apex + 1;
			}
			else {  
				if crev - src*cwp < cache + GAP;		// Bingo
				then {		// Record variable
					eye[apex] = i;
					jay[apex] = j;
					dir[apex] = 1;		// To fix threshold bound
					apex = apex + 1;
				}
			}
		}
	and(int k=0;k<apex;k++)			// Use results
		if dir[k] == 0; 
		then Plan[eye[k]][PROC][jay[k]] == ZERO;
		else Plan[eye[k]][PROC][jay[k]] >= THRESHOLD;

	max: Revenue; 				// For injury method
}