#define ITEMS 15					// Number of food items
#define LIMIT 3					// Number of servings allowed

// indices of nutitional categories
#define CALOR 0
#define PROT 1
#define FAT 2
#define SODIUM 3
#define VITA 4
#define VITC 5
#define VITB1 6
#define VITB2 7
#define NIACIN 8
#define CALCIUM 9
#define IRON 10

// Indices of food items
#define BURGER 0
#define LEANB 1
#define BIGB 2
#define FRIES 3
#define NUGGETS 4
#define HONEY 5
#define CSALAD 6
#define GSALAD 7
#define EGGMCM 8
#define CEREAL 9
#define YOGURT 10
#define MILK 11
#define OJ 12
#define GJ 13
#define AJ 14

continuous Serv[ITEMS];
continuous Cost;
continuous Fat, Calories;
continuous Variety, Y[ITEMS];
double prices[ITEMS];
int bv;
int food[ITEMS][IRON+1];
string iname[ITEMS];

data()
{
	// First Column of Table 11.1
	prices = { .59,1.79,1.65,.68,1.56,0.0,2.69,1.96,
				1.36,1.09,.63,.56,.88,.68,.68 };

	// Columns 2-12 of Table 11.1
	food = { 	255, 12, 9, 490, 4, 4, 20, 10, 20, 10, 15,
			320, 22, 10, 670, 10, 10, 25, 20, 35, 15, 20,
			500, 25, 26, 890, 6, 2, 30, 25, 35, 25, 20,
			220, 3, 12, 110, 0, 15,10, 0, 10, 0, 2,
			270, 20, 15, 580, 0, 0, 8, 8, 40, 0, 6,
			45, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
			170, 17, 9, 400, 100, 35, 20, 15, 20, 15, 8,
			50, 4, 2, 70, 90, 35, 6, 6, 2, 4, 8,
			280, 18, 11, 710, 10, 0, 30, 20, 20, 25, 15,
			90, 2, 1, 220, 20, 20, 20, 20, 20, 2, 20,
			105, 4, 1, 80, 2, 0, 2, 10, 2, 10, 0,
			110, 9, 2, 130, 10, 4, 8, 30, 0, 30, 0,
			80, 1, 0, 0, 0, 120, 10, 0, 0, 0, 0,
			80, 1, 0, 0, 0, 100, 4, 2, 2, 0, 0,
			90, 0, 0, 5, 0, 2, 2, 0, 0, 0, 4 
	};

	iname = {	"burger", "lean burger", "big burger", "fries", "nuggets", "honey",
				"chef salad", "garden salad", "egg sandwich", "cereal",
				"yogurt", "milk", "orange juice", "grapefruit juice", "apple juice"
			};
}

setup()
{
	and(int j=0;j<ITEMS;j++)  // Limit on number of servings 
		Serv[j] <= LIMIT;

	Serv[MILK] == Serv[CEREAL];			// Milk and cereal link

	Calories >= 2000;					// Minimum calories

	// At least 100% of these nutrients
	and(int i=VITA;i<=IRON;i++) 
		sigma(int j=0;j<ITEMS;j++)
			food[j][i]*Serv[j] >= 100; 	

	// At least 55 grams of protein
	sigma(int j=0;j<ITEMS;j++)food[j][PROT]*Serv[j] >= 55;

	// At most 3000 mg of sodium
	sigma(int j=0;j<ITEMS;j++)
		food[j][SODIUM]*Serv[j] <= 3000;	
 
	// Grams of fat
	Fat == 
	sigma(int j=0;j<ITEMS;j++) food[j][FAT]*Serv[j]; 
	// Total calories
	Calories == 
	sigma(int j=0;j<ITEMS;j++) food[j][CALOR]*Serv[j];
	// No more than 30% of calories from fat
	// given that each gram of fat has 9 calories
	9.0*Fat <= .30*Calories;

	// Objective function
	Cost == 
		sigma(int j=0;j<ITEMS;j++) prices[j]*Serv[j];	
}

find_least_expensive_meal()
{
	find_min: Cost; subject_to 
		and(int i=0;i<ITEMS;i++) 
			or(int k=0;k<=LIMIT;k++)
				Serv[i] == k;
}

extra_setup()
{
	and(int i=0;i<ITEMS;i++) {
		Y[i] <= 1.0;
		Y[i] <= Serv[i];
	}
	Variety == sigma(int i=0;i<ITEMS;i++) Y[i];
}

find_greatest_variety()
{
	find_max: Variety;
	subject_to {
		absgap(1.0);			// Fractional improvements impossible
		and(;;) {
			max: Variety;
			if determine_first_injured(bv);
			then split_in_two_halves(Serv[bv]);
			else break; 				// Eureka effect
		}
	}
}

split_in_two_halves(continuous X)
{
double cwp;

	cwp = wp(X);
	either X <= floor(cwp); 
	or X >= ceil(cwp);
}

#define fail 0==1

determine_first_injured(int & bv)
{
	and(int i=0;i<ITEMS;i++)
		if not integral(Serv[i]);
		then { 
			bv = i; 
			return; 
		}
	fail;
}

output()
{
	printf("Total cost is %.2f \n", Cost);
	printf("Total fat is %.0f \n", Fat);
	printf("Total calories are %.0f \n\n", Calories);

	and(int i=0;i<ITEMS;i++) {
		if Serv[i] == 0;
		then;
		else
		printf("Buy %.0f of item %s \n", Serv[i],iname[i]);
	}
}

2lp_main()
{
double cache;

	data();
	setup();
	if not {
		find_least_expensive_meal();				// As before
		cache = wp(Cost);
		printf("The best possible cost is %.2f\n",cache);
	}
	then { 
		printf("There is no solution\n");
		return;
	}
	extra_setup();
	Cost <= 1.5*cache;
	find_greatest_variety();
	output();
}
 