#define ARCS 23
#define NODES 18

int from_to[ARCS][2];
double arc_cap[ARCS];
double cons[ARCS];
double node_cap[NODES];
double node_demand[NODES];
continuous Flow[ARCS], Power[NODES];
int in[ARCS],nnz_in[NODES];
int out[ARCS],nnz_out[NODES];

2lp_main()
{
	data();				// From tables 4.3 and 4.4
	tables(1,in,nnz_in);		// Build auxiliary tables
	tables(0,out,nnz_out);
	setup();				// Node constraints
	min: Power[0];			// Minimize power produced at plant
	output();				// Left to reader
}

data()
{
	from_to = { 0, 1, 0, 2, 0, 3, 1, 4, 1, 5, 2, 5, 2, 10, 3, 6,
			3, 7, 3, 8, 4, 12, 4, 9, 5, 9, 6, 11, 7, 11,
			7, 15, 8, 16, 8, 17, 9, 12, 9, 13, 10, 14, 11, 14, 11, 15 };

	arc_cap = { 400.0, 400.0, 550.0, 330.0, 230.0, 340.0, 610.0, 330.0, 
			250.0, 320.0, 140.0, 350.0, 250.0, 300.0, 300.0, 250.0, 
			500.0, 220.0, 311.0, 233.0, 190.0, 278.0, 175.0 };

	cons = { .88, .9, .9, .92, .88, .89, .96, .99, .96, .92, .88, .89, .98, .99,
			.96, .98, .96, .92, .88, .89, .97, .99, .92 };

	node_cap = { 1200.0, 700.0, 700.0, 650.0, 600.0, 400.0, 520.0,
				625.0, 775.0, 621.0, 632.0, 798.0, 223.0, 310.0,
				421.0, 300.0, 254.0, 285.0 };

	node_demand = { 0.0, 110.0, 110.0, 110.0, 110.0, 110.0, 110.0, 110.0, 110.0,  
				110.0, 110.0,110.0, 222.0, 100.0, 120.0, 120.0, 111.0, 107.0 };

	in = { 0, 1, 2, 3, 4, 5, 7, 8, 9, 11, 12, 6, 13, 14, 10, 18, 19,
		20, 21, 15, 22, 16, 17 };

	out = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
			17, 18, 19, 20, 21, 22 };
}
 
tables(int dir,indices[],nnz[]) 		// Dir will be 1 for in and 0 for out
{
int end;
	end = 0;
	and(int n=0;n<NODES;n++) 
		and(int arc=0;arc<ARCS;arc++)
			if from_to[arc][dir] == n;
			then {
				indices[end] = arc;
				nnz[n] = nnz[n]+1;
				end = end + 1;
			}
}

setup()
{
	and(int n=0;n<NODES;n++) {
		Power[n] <= node_cap[n];
		Power[n] >= node_demand[n];
	}

	and(int arc=0;arc<ARCS;arc++)
		Flow[arc] <= arc_cap[arc];

	flows_in();
	flows_out();
}

flows_in()
{
int begin, end;
	end = 0; 
	and(int n=0;n<NODES;n++) { 
		begin = end; 
		end = end + nnz_in[n]; 
		if begin < end; 
		then		// Flow in is sum of flow along incoming arcs
			 Power[n] == 
			 sigma(int k=begin;k<end;k++) 
				cons[in[k]]*Flow[in[k]];
	}
}

flows_out()
{
int begin,end;

	end = 0; 
	and(int n=0;n<NODES;n++) { 
		begin = end;
		end = end + nnz_out[n]; 
		if begin < end; 
		then		// Flow out is equal to flow reaching node
			Power[n] == 
			sigma(int k=begin;k<end;k++) Flow[out[k]];
	}
}

output()
{
	and(int n=0;n<NODES;n++) 
		printf("N[%d] is %.2f\n",n,Power[n]);
	printf("\n");
	and(int n=0;n<ARCS;n++)
		printf("Flow[%d] is %.2f\n",n,Flow[n]);
}
 