#define NT 12							// Number of non-terminal nodes
double signed_reduced_cost[NT];			// To record reduced costs

#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();
	tables(1,in,nnz_in);
	tables(0,out,nnz_out);
	new_setup();
	if not first_pass();
	then error_hndlr(-1);				// Report and exit
	and(int p=1;p<NT;p++)
		if signed_reduced_cost[p] < 0.0;
		then 
			if not {		// Negation will restore feasible region
				printf("Analyzing node %d\n",p);
				keep_the_very_same(p);
				parametric_analysis(Power[p],Power[0],.5,110,0.0,-1);
			}
			then error_hndlr(p);			// Report and exit
}

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

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

new_setup()
{
	and(int n=0;n<NODES;n++) {
		Power[n] <= node_cap[n];
		// Power[n] >= node_demand[n]; FREE UP LOWER BOUNDS
	}
	and(int arc=0;arc<ARCS;arc++)
		Flow[arc] <= arc_cap[arc];
	flows_in();
	flows_out();
}

first_pass()
{
	and(int n=0;n<NODES;n++)	
		Power[n] >= node_demand[n];
	min: Power[0];
	and(int i=1;i<NT;i++)
		signed_reduced_cost[i] = rc(Power[i]);
}

keep_the_very_same(int p)
{
	and(int n=0;n<NODES;n++)	
		if n != p; 
		then Power[n] >= node_demand[n];
}

optimize(continuous Z)
{    
	min: Z;
}

sample(double & r, & c, b,continuous R,Z,int dir)
{
	dir*R <= dir*b;
	optimize(Z);
	r = range(R);
	c = fabs(rc(R));

	printf("At b=%.2f,",b);
	printf("Obj: %.2f, Res: %.2f, Rd Cst: %.2f, Rg: %.2f\n",Z,R,c,r);

	local_output();			// Application specific
}

check_bound(continuous R, int dir, double limit)
{
	if dir == 1;
	then ub(R) >= limit;
	else lb(R) <= limit;
}

parametric_analysis(continuous X,Z, double delta, b, limit, int dir)
{
double c,r;

	if not check_bound(X,dir,limit);
	then {
		printf("Limit exceeds known bound\n");
		return;
	}

	and(;;) {		// Loop until break
		if not 
			sample(r,c,b,X,Z,dir);		//Do and undo constraints					
		then {						// Adjustment for infeasibility
			printf("Infeasibility detected at b=%f\n",b);
			b = b + dir*limit/10;			// Move b along
		}
		else {
			if c == 0.0;				// If reduced cost is 0.0
			then break;				// end analysis
			update(b,r,delta,dir);		// Revise b
		}
		if dir*b > dir*limit;				// If no room left
		then break;					// exit the loop and end analysis
	}
}

update(double & b,  r, delta, int dir)
{
	if r <= delta; 
	then b = b + dir*delta;
	else b = b + dir*r;		
}


error_hndlr(int error)
{
	printf("Program exited: error condition #%d\n",error);
	exit();
}

local_output()
{
	;		// Null statement
}


 