# Perl module for schedule optimization generation


sub initFile {   ## 1 argument : string indicating filename of initialization data
    my $infilename = $_[0];   ## filename to load passed as arguement
    open INFILE, $infilename;
    my @data = <INFILE>;
    close INFILE;


    ## first pass through initialization file to grab important
    ## variables needed to initialize array variables properly
    foreach my $line (@data) {
	if ($line =~ /^\#/) { next; }	# skip comments
	if ($line =~ /^nn/) { $nn = parseVar($line) }  # number of nodes
	if ($line =~ /^ns/) { $ns = parseVar($line) }  # number of aircraft types / sizes
	if ($line =~ /^nt/) { $nt = parseVar($line) }  # number of transfer timesteps you want this solved in
	else { next; }		# skip unrecognized lines
    }


    ##
    ## Variables that need to be loaded by file:
    ##


    ## $initdemand 2D array
    ## TODO: allow additional demand to be added at later timesteps
    for ($i=0; $i<$nn; $i++) {
	for ($k=0; $k<$nn; $k++) {
	    # demand from node $i to node $k
	    $initdemand[$i][$k] = 0;
	}
    }

#    &print2Dlol("initdemand", @initdemand);	#DEBUG

    for ($i=0; $i<$nn; $i++) {
	$ncap[$i] = 1;	# vehicle capacity for each node per timestep
    }


    ## Vehicle setup variables:

    @fc =	( 10, 40 );	# Flight cost by aircraft type 
    @aac =	( 10, 40 );	# Aircraft acquisition cost by aircraft type 
    @fcap =     ( 10, 50 );	# Max seats per flight by aircraft type

    ## revenue per passenger delivered to destination
    for ($i=0; $i<$nn; $i++) {
	$sc[$i] = 0;
    }


    $maxaircraft = "Inf";	# total number of vehicles possibly in operation simultaneously
    $desirednumac = 0;	# preferred number of vehicles in operation simultaneously

    ## Solution variables
    $optcriteria = "min: numunsrv;";

    $dir = "output";  # default directory for output files, may be overridden
    mkdir($dir);


    ## parse initialization file
    foreach my $line (@data) {
	if ($line =~ /^\#/) { next; }	# skip comments
	if ($line =~ /^maxaircraft/) { $maxaircraft = parseVar($line) }
	if ($line =~ /^desirednumac/) { $desirednumac = parseVar($line) }
	if ($line =~ /^cyclicschedule/) { $cyclicschedule = parseVar($line) }
	if ($line =~ /^optcriteria/) { $optcriteria = parseVar($line) }
	if ($line =~ /^dir/) { $dir = parseVar($line) }
	if ($line =~ /^initdemand/) { parseSparseArray($line, \@initdemand); }
	if ($line =~ /^ncap/) { parseSparseArray($line, \@ncap); }
	if ($line =~ /^fc/) { parseSparseArray($line, \@fc); }
	if ($line =~ /^aac/) { parseSparseArray($line, \@aac); }
	if ($line =~ /^fcap/) { parseSparseArray($line, \@fcap); }
	if ($line =~ /^sc/) { parseSparseArray($line, \@sc); }
	if ($line =~ /^startpos/) { parseSparseArray($line, \@startpos); }
	else { next; }		# skip unrecognized lines

    }


#    &print2Dlol("initdemand", @initdemand);	#DEBUG

}


## take a line containing "name = value\n" and return just the value
sub parseVar {
    my ($name, $var) = split(/\s*=\s*/, $_[0]);
    chomp $var;
    print "Found $name = $var\n"; ## DEBUG
    return $var;
}

## take a line containing "name_i_j" = value\n" and an array reference
## and set the elements of the array reference [i][j]
sub parseSparseArray {
    my ($name, $var) = split(/\s*=\s*/, $_[0]);
    chomp $var;
    my $aref = $_[1];
    my @loc = split(/_/, $name);
    my $depth = $#loc;
    $name = shift(@loc);

#    print "Found $depth - dimensional $name = $var\n"; ## DEBUG
#    foreach my $iter (@loc) {
#	print " $iter ";
#    }

#    &print2Dlol("ar", @$aref);
    &recSetArr($aref, $var, @loc);
    return 1;
}

## recursively set multidimensional array elements
## array reference, a variable, and an array of subindexes
sub recSetArr {
    my $aref = shift;
    my $var = shift;
#    print "  Depth = $#_ = $_[0] \n";	#DEBUG
    if ($#_ <= 0) {
	$$aref[$_[0]] = $var;    my @ar = @$aref;

	return 1;
    }
    else {
	my $ind = shift @_;
	&recSetArr($$aref[$ind], $var, @_);
	return 0;
    }
}


## For debugging purposes
sub print2Dlol {
    my $name = shift;	# name text for display
    my @aref = @_;	# array reference
    print "\ns--------------------\n";
    print "$name = \n $aref[0][0]" ;
    for my $i ( 0 .. $#aref) {
	my $saref = $aref[$i];
	my $n = @$saref -1;
	for $j (0 .. $n) {
	    print " \[$i $j = $aref[$i][$j] \]";
	}
	print "\n";
    }
    print "\ne--------------------\n";
}



## read text containing a multidimensional array "[[x x][x x]]" and
## return and array
## performed recursively, so it works regardless of the dimension of the array
## FIXME: this is horribly broken -- needs to be rewritten using
## ungreedy regexps instead of greedy index($string)
sub parseArray {
#    my @el = split(/[0-9]*\.[0-9]*\s+[0-9]*\.[0-9]*/, $_[0]);  ## probe structure of braces
#   my @num = split(/\]*\s*\[*/, $_[0];		## store numbers only
#    my $dim = $#el[0];
#    foreach $seg (@el) {
#	
#    }

    ## inside loop, no more dimensions to process
    if (!(substr($_[0],0,1) eq '[')) {
	return split(/\s+/, $_[0]);
    }
    ## recursively string together subarrays
    else {
	my @ar = ();
	my $cb = ']';
	# working copy with front brace stripped
#	my $artext = $_[0];
	my $artext = substr($_[0], 1);
	chop $artext;  # remove trailing brace
	while (length($artext) > 0 ) {	# until we reach the end of the string 
#	    print "\$artext " . length($artext) . " = \'$artext\'\n";	#DEBUG
	    $artext =~ m/^(\[+)/ ;  # match one or more opening braces
	    $cb = $1;	# a number of opening braces
	    # turned into the same number of closing braces
	    $cb =~ s/\[/\]/g;	
	    $closepos = index($artext, $cb);
#	    print "\$cb = $cb; artextlength: " . length($artext) . "\n";	#DEBUG
	    my $subar = substr($artext, 0, $closepos);
	    print "\$subar = \'$subar\'\n";
	    push (@ar, &parseArray($subar));
	    $artext = substr($artext, $closepos);
	}
    }
    return @ar;
}



sub parametricInit {

##
##	Parametric Initialization Constants
##

## FIXME:  should read numbers from an input file rather than hardcode them in this script

# The traveling population
    $population = 2500;
    
    
##
##	Inputs
##

## hub / node ratio
    $r_hub_node = 1;	# range from 1.0 <= x
    
# number of aircraft that can arrive/leave in one timestep
# corresponds to the number of gates / loading docks
    $node_maxthroughput = 3;	
    $hub_maxthroughput = $node_maxthroughput * $r_hub_node;
    
# population / number of nodes ratio
    $nn = 4;  # number of nodes (airports) in the system [0 < x < 100] 
    $nh = 1;  # number of hub nodes		[0 < x < $nn]
    $nt = 4;  # number of time increments	[1 < x < 10]
    $ns = 1;  # number of aircraft types	[0 < x < 10]
    $r_pop_numnodes = $population / $nn;

# population / number of vehicles ratio
    $maxaircraft = 7;
    $desirednumac = 2;
    $r_pop_aircraft = $population / $maxaircraft;
    
    
## initial demand for travel to each node
# scale population per node relative to node size - hubs get more
# passengers proportional to hub/node ratio
    
    $hoph2h = $nh * ($nh - 1);  # hub to hub segments
    $hopn2n = ($nn - $nh) * ($nn - $nh - 1);	# node to node segments
    $hoph2n = $nn * ($nn - 1) - $hoph2h - $hopn2n;	# hub to/from node segments
## nominal ordinary none to node demand
    $defdemand = $population / ( $hoph2h * $r_hub_node**2 
				 + $hoph2n * $r_hub_node 
				 + $hopn2n);
    
    for ($i=0; $i<$nn; $i++) {
	for ($k=0; $k<$nn; $k++) {
	    if ($i != $k) {
		# scale population relative to hub size
		## int($x + 0.5) is a cheesy substitute for perl's lack of
		## a round() function
		if ($i < $nh && $k < $nh) {  # hub to hub
		    $initdemand[$i][$k] = int($defdemand * $r_hub_node**2 + 0.5);
		} elsif ($i < $nh || $k < $nh) {  # node to/from hub
		    $initdemand[$i][$k] = int($defdemand * $r_hub_node + 0.5);
		} else {  # node to node
		    $initdemand[$i][$k] = int($defdemand + 0.5);
		}
	    }
	}
    }
    
## Revenue per passenger delivered src to dest
    $defsc = 200;
    for ($i=0; $i<$nn; $i++) {
	$sc[$i] = $defsc;
    }
    
## Runway capacity at each airport node, landings + takeoffs per timestep
    for ($i=0; $i<$nn; $i++) {
	if ($i<$nh) {
	    $ncap[$i] = $hub_maxthroughput;
	} else {
	    $ncap[$i] = $node_maxthroughput;
	}
    }
    
## FIXME:  currently hardcoded to use only 1 type of aircraft!
## parameter for each aircraft type s 
# $dur = 1;  # flight duration between nodes in timesteps
    @fc =	( 17000, 42000 );	# Flight cost by aircraft type 
    @aac =	( 10000, 40000 );	# Aircraft acquisition cost by aircraft type 
    @fcap = ( 200,	550 );	# Max seats per flight by aircraft type

    $cyclicschedule = 1;	## Aircraft must end up at the same starting points
    
## Number of aircraft of each type (row) starting at each node (col)
#$defsupply = 3;
#for ($s=0; $s<$ns; $s++) {
#    for ($i=0; $i<$nn; $i++) {
#	$initsupply[$s][$i] = $defsupply;
#    }
#}
    $dir = "output";  # default directory for output files, may be overridden
    mkdir($dir);
    
    ## Set an optimization criteria
#    $optcriteria = "max: profit;    /*  maximize efficiency  */\n";
    $optcriteria = "min: numunsrv ;  /*  number of unserved passengers  */\n";
#    $optcriteria = "min: numacdev;  /*  deviation from number of aircraft in service */\n";



    
}  # sub parametricInit 





## TODO: Add topologies, so network is not always fully connected -
## remove vars for pathways that don't exist

## TODO: Add concept of time - nodes can be different distances apart,
## and have timesteps correspond to a length of time to varying
## resolutions

## TODO: Add ability to insert more demand at times other than 0

## TODO: Add ability to run multimodal graphs by allowing certain
## nodes to only service certain vehicle types

## TODO: Initialization with vehicles in transit


##
##	LP formulation starts here
##
sub writeModel {

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

    print LPM $optcriteria, "\n";
    
    
## FIXME:  currently hardcoded with a limitation of 99 nodes and 9 time steps,
## since glpk can't process variables longer than 8 characters.  Need
## to use separaters instead 
    print LPM "
/*
 Three airports, three time increments
 p[_t_i_j_k] : people from node i going to node j at time t with final destination k
 f[_t_s_i_j] : flights of type s from node i going to node j at time t
 pp[_t_i_k] : pool of people at node i at time t whose destination is node k
 ap[_t_s_i] :  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 */\n";
    for ($t=0; $t<=$nt; $t++) {
	for ($i=1; $i<=$nn; $i++) {
	    for ($s=1; $s<=$ns; $s++) {
		## Only set if provided
		if ($startpos[$t][$s][$i]) {
		    &mkdd($i);
		    print LPM "ap$t$s$i = $startpos[$t][$s][$i] ; \n";
		}
	    }
	}
    } print LPM "\n";


    if ($cyclicschedule) {
	print LPM "/* final positions = initial aircraft starting positions */\n";
	$t = 0;
	for ($i=1; $i<=$nn; $i++) {
	    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++) {
    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++) {
	    $tmo = $t - 1;  # previous timestep
	    for ($j=1; $j<=$nn; $j++) {
		if ($i != $j) {	# no flights to nowhere
		    &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
##

} # sub writeModel




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

}


sub summarizeResults {
## 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
");
}



## 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 @_;
}

 
return 1;
