#!/usr/bin/perl -w
# Generate a passenger-centric LP to optimize an airline flight scedule
# Rowin Andruscavage
# ENCE 667 SP02 Dr. S. Gabriel
#
# 5/15/2003  Add input/output metrics for ENPM642 trade study
# 8/12/2005  Modify for arcology simulation transit optimization

use schedopt.pm

## if a input filename is provided, read initialization from it
if $argv[0] {
    print "unable to read $_";
}

# otherwise use the parametric initialization routines
else { &parametricInit };


## Multicriteria optimization loops
#@opt_bases = ("numacdev", "numunsrv", "profit");
@opt_bases = ("profit");  # perform a single run
foreach $opt_base (@opt_bases) {

    ## re-sort @opt_bases array so the current entry goes first
    @opt_bases2 = @opt_bases;
    until ($opt_base eq $opt_bases2[0]) {
	push(@opt_bases2, shift(@opt_bases2));
	## actually, we can just get rid of half or our cases
#	shift(@opt_bases2);
    }
#    print join(':', @opt_bases) . "\n";
#    print join(':', @opt_bases2). "\n";
    
    foreach $opt_subcrit (@opt_bases2) {
	if ($opt_subcrit eq $opt_base) {
	    $dir = "$opt_base";
	} else {
	    $dir = "$opt_base/$opt_subcrit";
	}
	mkdir($dir);


##
## Tweak base optimal solutions as constraints
##


	## number of points to try to relax values at
	## correspondes to resolution of the tradeoff curve
	$numpts = 1;
	$relax = 0;
	$first = 0;
	if (!($opt_subcrit eq $opt_base)) {
	    $relax = 1;  ## Turn on the relaxation loop
	    $first = 1;
## read the optimal value of the base optimization
	    %baseresults = &readresults("$opt_base/output.csv");
	    $basevalue = $baseresults{$opt_base};
	}	    

	$value = "";
	foreach (my $relx=0.0; $relx<1.0; $relx+=1.0/$numpts) {

	if ($relax == 1) {
# profit
	    if ($opt_base eq "profit") {
		my $newvalue = $basevalue - 2 * $basevalue * $relx;
		if ($value == $newvalue) {
		    unless ($first == 1) {
			next;  # no change, skip to next $relx
		    }
		}
		$value = $newvalue;
		$relaxbaseopt = "profit >= $value;\n";
	    }
	    
# number of unserved people
	    elsif ($opt_base eq "numunsrv") {    
		my $newvalue = int($basevalue + ($population - $basevalue) 
			     * $relx + 0.5);
		if ($value == $newvalue) {
		    unless ($first == 1) {
			next;  # no change, skip to next $relx
		    }
		}
		$value = $newvalue;
		$relaxbaseopt = "numunsrv <= $value; \n";
	    }
	    
# deviation from standard fleet size
	    elsif ($opt_base eq "numacdev") {
		my $maxdev = ($desirednumac > ($maxaircraft - $desirednumac) ? 
		    $desirednumac : $maxaircraft - $desirednumac);
		my $newvalue = int($basevalue + $maxdev * $relx + 0.5);
		if ($value == $newvalue) {
		    unless ($first == 1) {
			next;  # no change, skip to next $relx
		    }
		}
		$value = $newvalue;
		$relaxbaseopt = "numacdev <= $value; \n";
	    }
	}


##
##	LP formulation starts here
##

	
%input = (
	  population => "$population",
	  r_hub_node => "$r_hub_node",
	  node_maxthroughput => "$node_maxthroughput",
	  hub_maxthroughput => "$hub_maxthroughput",
	  numnodes => "$numnodes",
	  numhubs => "$numhubs",
	  r_pop_numnodes => "$r_pop_numnodes",
	  maxaircraft => "$maxaircraft",
	  desirednumac => "$desirednumac",
	  r_pop_aircraft => "$r_pop_aircraft",
	  time_incr => "$nt",
	  defsc => "$defsc",
	  fc => "$fc[0]",
	  aac => "$aac[0]",
	  fcap => "$fcap[0]",
	  );

print "\n\n--------------------------------------------------------------
Creating lp in $dir ...\n";
open (LPM, ">$dir/sched.lpm");


##
## output metrics
##

# profit
if ($opt_subcrit eq "profit" || ($opt_base eq "profit" && $opt_subcrit eq "profit")) {
    print LPM "max: profit;\n";
}

# number of unserved people
elsif ($opt_subcrit eq "numunsrv" || ($opt_base eq "numunsrv" && $opt_subcrit eq "numunsrv")) {    
    print LPM "min: numunsrv ;  /*  number of unserved passengers  */\n";
}

# deviation from standard fleet size
elsif ($opt_subcrit eq "numacdev" || ($opt_base eq "numacdev" && $opt_subcrit eq "numacdev")) {
    print LPM "min: numacdev;  /*  number of aircraft in service */\n";
}

	## relax base optimal value constraints
	if ($relax == 1) {
	    print LPM $relaxbaseopt;
	}
	

## FIXME:  currently has a limitation of 99 nodes and 9 time steps,
## since glpk can't process variables longer than 8 characters.

    print LPM "
/*
 Three airports, three time increments
 p[tiijjkk] : people from node i going to node j at time t with final destination k
 f[tsiijj] : flights of type s from node i going to node j at time t
 pp[tiikk] : pool of people at node i at time t whose destination is node k
 ap[tsii] :  pool of aircraft of type s at node i at time t
*/\n\n";


print LPM "/* Profit = passenger revenue - total flight cost 
     - (aircraft acquisition cost * number of aircraft)
*/
tpr - tfc - $aac[0] * numac = profit ;\n";
$t = $nt;
for ($k=1; $k<=$nn; $k++) {
    &mkdd($k);
    print LPM "+ $sc[$k-1] pp$t$k$k ";
}
print LPM "- tpr = 0 ; \n\n";


print LPM "/* Total flight cost */\n";
for ($t=1; $t<=$nt; $t++) {
    for ($i=1; $i<=$nn; $i++) {
	for ($j=1; $j<=$nn; $j++) {
	    if ($i != $j) {
		for ($s=1; $s<=$ns; $s++) {
		    &mkdd($i,$j);
		    print LPM "+ $fc[$s-1] f$t$s$i$j ";
		}
	    }
	}
    }
    print LPM "\n";
} print LPM "- tfc = 0 ; \n\n";


print LPM "/* Initial distribution of demand */\n";
$t=0;
for ($i=1; $i<=$nn; $i++) {
    for ($k=1; $k<=$nn; $k++) {
	if ($i != $k) {
	    &mkdd($i,$k);
	    # initial demand for travel from/to
	    print LPM "pp$t$i$k = $initdemand[$i-1][$k-1] ; \n";
	} else {
	    &mkdd($k);
	    # initial pool of passengers at their final destination
	    print LPM "pp$t$k$k = 0 ; \n";
	}
    }
} print LPM "\n";

print LPM "/* unserved persons remaining at time = $nt */\n";
$t=$nt;
for ($i=1; $i<=$nn; $i++) {
    for ($k=1; $k<=$nn; $k++) {
	if ($i != $k) {
	    &mkdd($i,$k);
	    print LPM "+ pp$t$i$k ";
	} 
    ## else unconstrained final pool of passengers at their final destination
    }
} print LPM " = numunsrv ; \n\n";


print LPM "/* initial aircraft starting positions = final positions */\n";
$t = 0;
for ($i=1; $i<=$nn; $i++) {
    if ($i != $j) {
	for ($s=1; $s<=$ns; $s++) {
	    &mkdd($i);
	    print LPM "ap$t$s$i = ap$nt$s$i ; \n";
	}
    }
} print LPM "\n";


print LPM "/* max aircraft */\n";
$t=0;
for ($i=1; $i<=$nn; $i++) {
    if ($i != $j) {
	for ($s=1; $s<=$ns; $s++) {
	    &mkdd($i);
	    print LPM "+ ap$t$s$i ";
	}
    }
} print LPM " - numac = 0 ; 
numac <= $maxaircraft ; 
spareacp + spareacm - numacdev = 0;
spareacp >= 0;
spareacm >= 0;
numac + spareacp - spareacm = $desirednumac;
\n\n";


print LPM "/* Apply flight capacities */\n";
for ($t=1; $t<=$nt; $t++) {
    for ($i=1; $i<=$nn; $i++) {
	for ($j=1; $j<=$nn; $j++) {
	    if ($i != $j) {	# no flights to nowhere
		for ($k=1; $k<=$nn; $k++) {
		    if ($i != $k) {	# no one leaves their destination
			&mkdd($i,$j,$k);
			print LPM  "+ p$t$i$j$k ";
		    }
		}
		for ($s=1; $s<=$ns; $s++) {
		    &mkdd($i,$j);
		    print LPM "- $fcap[$s-1] f$t$s$i$j ";
		}
		print LPM "< 0 ; \n";
	    }
	}
    }
    print LPM "\n";
} print LPM "\n";


print LPM "/* Apply airport capacities */\n";
for ($t=1; $t<=$nt; $t++) {
    for ($i=1; $i<=$nn; $i++) {
	for ($j=1; $j<=$nn; $j++) {
	    if ($i != $j) {	# no flights to nowhere
		$tmo = $t - 1;
		for ($s=1; $s<=$ns; $s++) {
		    &mkdd($i,$j);
		    print LPM "+ f$t$s$i$j + f$tmo$s$j$i ";
		}
		print LPM " < $ncap[$i-1] ; \n";
	    }
	}
    }
    print LPM "\n";
} print LPM "\n";


# currently assumes flight duration $dur = 1
print LPM "/* Apply aircraft movements */\n";

for ($s=1; $s<=$ns; $s++) {
    for ($t=1; $t<=$nt; $t++) {
	for ($i=1; $i<=$nn; $i++) {
	    for ($j=1; $j<=$nn; $j++) {
		if ($i != $j) {	# no flights to nowhere
		    $tmo = $t - 1;  # previous timestep
		    &mkdd($i,$j);
		    print LPM "- f$t$s$i$j + f$tmo$s$j$i ";
		}
	    }
	    &mkdd($i);
	    print LPM "+ ap$tmo$s$i - ap$t$s$i = 0 ; \n";
	} 
    }
} print LPM "\n\n";


print LPM "/* Apply time 0 constraints - no flights */\n";
$t=0;
for ($i=1; $i<=$nn; $i++) {
    for ($j=1; $j<=$nn; $j++) {
	if ($i != $j) {
	    for ($s=1; $s<=$ns; $s++) {
		&mkdd($i,$j);
		print LPM "f$t$s$i$j = 0 ; \n";
	    }
	}
    }
} print LPM "\n\n";


# Assumes $dur = 1;
print LPM "/* Apply passenger movements */\n";
for ($t=1; $t<=$nt; $t++) {
    $tmo = $t - 1;  # previous timestep
    for ($k=1; $k<=$nn; $k++) {
	for ($i=1; $i<=$nn; $i++) {
	    if ($i != $k) {	# no passengers leave their destination
		for ($j=1; $j<=$nn; $j++) {
		    if ($j != $i) {	# no passengers going nowhere
			# passengers leaving
			&mkdd($i,$j,$k);
			print LPM "- p$t$i$j$k ";
			if ($j != $k) { # no passengers leave their destination
			    # passengers arriving
			    &mkdd($i,$j,$k);
			    print LPM "+ p$tmo$j$i$k ";
			}
		    }
		}
		# update passenger pools
		&mkdd($i,$k);
		print LPM "+ pp$tmo$i$k - pp$t$i$k = 0 ; \n";
	    } else {  # pool of people who arrived at their destination
		for ($j=1; $j<=$nn; $j++) {
		    if ($j != $k) {	# no passengers going nowhere
			# number of passengers arriving destination
			&mkdd($j,$k);
			print LPM "+ p$tmo$j$k$k ";
		    }
		}
		# passenger pool who arrived at their destination
		&mkdd($k);
		print LPM "+ pp$tmo$k$k - pp$t$k$k = 0 ; \n";
	    }
	}
    } 
} print LPM "\n\n";


print LPM "/* Apply time 0 constraints - no passenger movements */\n";
$t=0;
for ($i=1; $i<=$nn; $i++) {
    for ($j=1; $j<=$nn; $j++) {
	if ($i != $j) {		# no passengers going nowhere
	    for ($k=1; $k<=$nn; $k++) {
		if ($i != $k) {	# no passengers leaving their destination
		    &mkdd($i,$j,$k);
		    print LPM "p$t$i$j$k = 0 ; \n";
		}
	    }
	}
    } 
} print LPM "\n\n";


print LPM "/* Apply integer constraint to flights */\n";
print LPM "int ";
for ($t=1; $t<=$nt; $t++) {
    for ($i=1; $i<=$nn; $i++) {
	for ($j=1; $j<=$nn; $j++) {
	    if ($i != $j) {
		for ($s=1; $s<=$ns; $s++) {
		    &mkdd($i,$j);
		    print LPM "f$t$s$i$j";
		    unless (($t==$nt) && ($i==$nn) && ($j==$nn-1) && ($s==$ns)) {
			print LPM " , ";
		    }
		}
	    }
	}
    } 
} print LPM " ;\n\n";

close (LPM);

##
##	End of LP formulation
##


## DEBUG  
#	if( 0 ) {

print "Solving LP ...\n";
system("lp2mps < $dir/sched.lpm > $dir/sched.mps; glpsol -o $dir/sched.ogl $dir/sched.mps");

## Make a human-readable results summary for debugging
print "Summarizing results from $dir/sched.ogl to $dir/results.txt for\n";
system("
echo \"Results\" > $dir/results.txt
grep profit $dir/sched.ogl >> $dir/results.txt
grep tfc $dir/sched.ogl >> $dir/results.txt
grep tpr $dir/sched.ogl >> $dir/results.txt
grep numac $dir/sched.ogl >> $dir/results.txt
grep numunsrv $dir/sched.ogl >> $dir/results.txt

echo \"Flight schedule:\" >> $dir/results.txt
grep f[0-9] $dir/sched.ogl | grep \" [1-9] \" | sort -k 2 >> $dir/results.txt

echo \"Passenger pools\" >> $dir/results.txt
grep pp[0-9] $dir/sched.ogl | sort -k 2 >> $dir/results.txt

echo \"Aircraft pools\" >> $dir/results.txt
grep ap[0-9] $dir/sched.ogl | grep \" [1-9] \" | sort -k 2 >> $dir/results.txt
");

print "Collecting results from $dir/sched.ogl to $dir/output.csv ";
open(RESULTS, "$dir/sched.ogl");
@rawresults = <RESULTS>;
close(RESULTS);

@results = (grep(/profit/, @rawresults), 
	    grep(/tfc/, @rawresults),
	    grep(/tpr/, @rawresults),
	    grep(/numac/, @rawresults),
	    grep(/numacdev/, @rawresults),
	    grep(/numunsrv/, @rawresults));
# print join("\n", @results);

%output = ();
foreach $line (@results) {
    @field = split(/\s+/,$line);
#    print join(':', @field) . "\n";
    $output{$field[2]} = $field[3];
}


	if ($relax == 0 || ($relax == 1 && $first == 1)) {
	    open (OUTPUT, ">$dir/output.csv");
	    ## print line of column numbers
	    print OUTPUT "#";
	    my $xiter = 0;
	    foreach $key (keys %input) {
		$xiter++;
		print OUTPUT "$xiter,";
	    }
	    foreach $key (keys %output) {
		$xiter++;
		print OUTPUT "$xiter,";
	    }
	    print OUTPUT "\n";

	    ## print line of keys
	    print OUTPUT "#";
	    foreach $key (keys %input) {
		print OUTPUT "$key,";
	    }
	    foreach $key (keys %output) {
		print OUTPUT "$key,";
	    }
	    if ($relax == 1) {
		chomp $relaxbaseopt;
		print OUTPUT $relaxbaseopt; 
	    }
	    print OUTPUT "\n";
	    close (OUTPUT);
	} 
	
	open (OUTPUT, ">>$dir/output.csv");
	foreach $key (keys %input) {
	    print OUTPUT "$input{$key},";
	}
	foreach $key (keys %output) {
	    print OUTPUT "$output{$key},";
	}
	if ($relax == 1) {
	    $first = 0;  ## not first relax loop anymore
 ## record how this line was constrained
	    print OUTPUT "$value"; 
	}
	## no trailing comma
	print OUTPUT "\n";
	close (OUTPUT);

## only loop through once if not playing with relaxed constraints
	unless ($relax == 1) {
	    last;  
	}
    }   # foreach $relx 
}	# foreach $opt_subcrit
}	# foreach $opt_base


#}  ## DEBUG

## change numeric vars to double digit format
sub mkdd {
    foreach (@_) {
	if ($_ > 99) { 
	    print STDERR "Overflow: $_ can't be expressed in two digits!\n";
	    exit; 
	}
	$_ = sprintf("%02d", $_);
    }
    return @_;
}


## Reads first two lines of a csv file (name given in $_[0]) and
## returns a hash of its keys => values
sub readresults {
    open (BASERESULTS, $_[0]);
    my @baseresults = <BASERESULTS>;
    close (BASERESULTS);
    
    ## skip first line are just column numbers
    my @bkeys = split(/,/, $baseresults[1]);
    my @bvals = split(/,/, $baseresults[2]);
    
    my %baseresults = ();
    for (my $x=0; $x<$#bkeys; $x++) {
	$baseresults{$bkeys[$x]} = $bvals[$x];
    }
    return %baseresults;
}
