# 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 vehicle 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 3D array
     for ($t=0; $t<$nt; $t++) {
 	for ($i=0; $i<$nn; $i++) {
 	    for ($k=0; $k<$nn; $k++) {
 		# demand from node $i to node $k
 		$initdemand[$t][$i][$k] = 0;
 	    }
 	}
     }


    ## $tg transitgraph 3D array
    for ($s=0; $s<$ns; $s++) {   # different transit connectivity matrix for each vehicle type
	for ($i=0; $i<$nn; $i++) {	# source node
	    for ($j=0; $j<$nn; $j++) {	# destination node
		$transitgraph[$s][$i][$j] = 0;  # 0 = no connection, 1 = number of timesteps
	    }
	}
    }

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

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


    ## Vehicle setup variables:

    ## flight cost 3D array (type, source, destination)
     for ($s=0; $s<$ns; $s++) {
 	for ($i=0; $i<$nn; $i++) {
 	    for ($j=0; $j<$nn; $j++) {
 		# cost for aircraft type $s from node $i to node $j
 		$fc[$s][$i][$j] = 1;
 	    }
 	}
	$vac[$s] =	0;	# additional vehicle acquisition / lease cost by vehicle type 
	@fcap[$s] =     0;	# Max seats per flight by vehicle type
     }


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


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

    ## Solution variables
    $optcriteria = "max: numsrv;";

    $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 =~ /^maxvehicle/) { $maxvehicle = parseVar($line) }
	if ($line =~ /^desirednumveh/) { $desirednumveh = 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 =~ /^transitgraph/) { parseSparseArray($line, \@transitgraph); }
	if ($line =~ /^ncap/) { parseSparseArray($line, \@ncap); }
	if ($line =~ /^fcost/) { parseSparseArray($line, \@fc); }
	if ($line =~ /^vac/) { parseSparseArray($line, \@vac); }
	if ($line =~ /^fcap/) { parseSparseArray($line, \@fcap); }
	if ($line =~ /^sc/) { parseSparseArray($line, \@sc); }
	if ($line =~ /^startpos/) { parseSparseArray($line, \@startpos); }
	else { next; }		# skip unrecognized lines

    }


}


## 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 vehicle 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 station nodes 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 vehicle types	[0 < x < 10]
    $r_pop_numnodes = $population / $nn;

# population / number of vehicles ratio
    $maxvehicle = 7;
    $desirednumveh = 2;
    $r_pop_vehicle = $population / $maxvehicle;
    
    
## 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 ($t=0; $t<1; $t++) {	# Only set for $t=0 for now
	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[$t][$i][$k] = int($defdemand * $r_hub_node**2 + 0.5);
		    } elsif ($i < $nh || $k < $nh) {  # node to/from hub
			$initdemand[$t][$i][$k] = int($defdemand * $r_hub_node + 0.5);
		    } else {  # node to node
			$initdemand[$t][$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 station node, arrivals + departures per timestep
    for ($i=0; $i<$nn; $i++) {
	for ($s=0; $s<$ns; $s++) {
	    if ($i<$nh) {
		$ncap[$i][$s] = $hub_maxthroughput;
	    } else {
		$ncap[$i][$s] = $node_maxthroughput;
	    }
	}
    }
    
## FIXME:  currently hardcoded to use only 1 type of vehicle!
## parameter for each vehicle type s 
# $dur = 1;  # flight duration between nodes in timesteps
    @fc =	( 17000, 42000 );	# Flight cost by vehicle type 
    @vac =	( 10000, 40000 );	# Vehicle acquisition cost by vehicle type 
    @fcap = ( 200,	550 );	# Max seats per flight by vehicle type

    $cyclicschedule = 1;	## Vehicle must end up at the same starting points
    
## Number of vehicles 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 = "max: numsrv ;  /*  number of served passengers  */\n";
#    $optcriteria = "min: numvehdev;  /*  deviation from number of vehicle in service */\n";



    
}  # sub parametricInit 



## TODO: Add concept of time - nodes can be different distances apart,
## and have timesteps correspond to a length of time to varying
## resolutions.  This can be accomplished by introducing "non-station
## nodes" in between links or to bypass stations. 
## Bypass stations also means that we want transfers to take a
## configurable number of timesteps, rather than assuming it can all
## take place in a single timestep.

## 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";    
    print LPM "
/*
 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
 ppi[_t_i_k] : pool of people initialized to depart node i at time t with destination node k
 ap[_t_s_i] :  pool of vehicle of type s at node i at time t
*/\n\n";

    print LPM "/* Profit = passenger revenue - total flight cost */\n";
#    print LPM "tpr - tfc = profit ;\n";   # doesn't do good for some reason
    $t = $nt;
    for ($k=0; $k<$nn; $k++) {
	print LPM "+ $sc[$k] pp\_$t\_$k\_$k ";
    }
    print LPM "- tpr = 0 ; \n\n";


    ## TODO:  add "sitting in place" cost
    print LPM "/* Total flight cost */\n";
    for ($t=1; $t<=$nt; $t++) {   # skip t=0, since we don't have
				  # control of initialization parameters
	for ($i=0; $i<$nn; $i++) {
	    for ($j=0; $j<$nn; $j++) {
		if ($i != $j) {
		    for ($s=0; $s<$ns; $s++) {
			if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			    print LPM "+ $fc[$s][$i][$j] f\_$t\_$s\_$i\_$j ";
			}
		    }
		}
	    }
	}
	print LPM "\n";
    } 
    print LPM "- tfc = 0 ; \n\n";

    
    print LPM "/* Initial distribution of demand */\n";
    for ($i=0; $i<$nn; $i++) {
	for ($k=0; $k<$nn; $k++) {
	    if ($i != $k) {
		for ($t=0; $t<$nt; $t++) {
		    # initial demand for travel from/to
		    print LPM "ppi\_$t\_$i\_$k = $initdemand[$t][$i][$k] ; \n";
		}
	    }
	    ## Removed these variables entirely
	    #else {
	    #    # initial pool of passengers at their final destination at $t=0
	    #    print LPM "pp\_0\_$k\_$k = 0 ; \n";
	    #}
	}
    } print LPM "\n";
    
    print LPM "/* served persons */\n";
    $t=$nt;  # at the final timestep under consideration
    for ($i=0; $i<$nn; $i++) {
	print LPM "+ pp\_$t\_$i\_$i ";
    } 
    print LPM " = numsrv ; \n\n";
    
    print LPM "/* sum of waiting persons who are not at their destinations over every timestep */\n";
    for ($t=1; $t<=$nt; $t++) {  # we have no control over $t=0
	for ($i=0; $i<$nn; $i++) {
	    for ($k=0; $k<$nn; $k++) {
		if ($i != $k) {
		    print LPM "+ pp\_$t\_$i\_$k ";
		}
	    } 
	}
	print LPM "\n";
    }
    print LPM " = numwaiting ; \n\n";

    
    
    print LPM "/* initial vehicle starting positions */\n";
    for ($t=0; $t<=$nt; $t++) {
	for ($i=0; $i<$nn; $i++) {
	    for ($s=0; $s<$ns; $s++) {
		## Only set if provided
		if ($startpos[$t][$s][$i]) {
		    print LPM "ap\_$t\_$s\_$i = $startpos[$t][$s][$i] ; \n";
		}
	    }
	}
    } print LPM "\n";


    ## for repeatable cyclic demand & schedules
    if ($cyclicschedule) {
	print LPM "/* final positions = initial vehicle starting positions */\n";
	$t = 0;
	for ($i=0; $i<$nn; $i++) {
	    for ($s=0; $s<$ns; $s++) {
		print LPM "ap\_$t\_$s\_$i = ap\_$nt\_$s\_$i ; \n";
	    }
	} print LPM "\n";
    }

    print LPM "/* max vehicles */\n";
    $t=0;
    for ($i=0; $i<$nn; $i++) {
	for ($s=0; $s<$ns; $s++) {
	    print LPM "+ ap\_$t\_$s\_$i ";
	}
    } print LPM " - numveh = 0 ; 
numveh <= $maxvehicle ; 
sparevehp + sparevehm - numvehdev = 0;
sparevehp >= 0;
sparevehm >= 0;
numveh + sparevehp - sparevehm = $desirednumveh;
\n\n";
    
    
    print LPM "/* Apply flight capacities */\n";
    for ($t=1; $t<=$nt; $t++) {  # no control over $t=0
	for ($i=0; $i<$nn; $i++) {
	    for ($j=0; $j<$nn; $j++) {
		if ($i != $j) {	# no flights to nowhere
		    for ($k=0; $k<$nn; $k++) {
			if ($i != $k) {	# no one leaves their destination
			    print LPM  "+ p\_$t\_$i\_$j\_$k ";
			}
		    }
		    for ($s=0; $s<$ns; $s++) {
			if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			    print LPM "- $fcap[$s] f\_$t\_$s\_$i\_$j ";
			}
		    }
		    print LPM "<= 0 ; \n";
		}
	    }
	}
	print LPM "\n";
    } print LPM "\n";
    
    
    print LPM "/* Apply station berthing capacities */\n";
    for ($t=1; $t<=$nt; $t++) { 
	for ($i=0; $i<$nn; $i++) {
	    for ($s=0; $s<$ns; $s++) {
		for ($j=0; $j<$nn; $j++) {
		    if ($i != $j) {	# no flights to nowhere
			$tmo = $t - 1;
			if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			    print LPM "+ f\_$t\_$s\_$i\_$j + f\_$tmo\_$s\_$j\_$i ";
			}
		    }
		}
		print LPM " < $ncap[$i][$s] ; \n";
	    }
	}
	print LPM "\n";
    } print LPM "\n";
    
    
    # currently assumes flight duration $dur = 1
    print LPM "/* Apply vehicle movements */\n";
    
    for ($s=0; $s<$ns; $s++) {
	for ($t=1; $t<=$nt; $t++) {
	    for ($i=0; $i<$nn; $i++) {
		$tmo = $t - 1;  # previous timestep
		for ($j=0; $j<$nn; $j++) {
		    if ($i != $j) {	# no flights to nowhere
			if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			    print LPM "- f\_$t\_$s\_$i\_$j + f\_$tmo\_$s\_$j\_$i ";
			}
		    }
		}
		print LPM "+ ap\_$tmo\_$s\_$i - ap\_$t\_$s\_$i = 0 ; \n";
	    } 
	}
    } print LPM "\n\n";
    

    ## TODO: make this an initialization variable
    print LPM "/* Apply time 0 constraints - no flights */\n";
    $t=0;
    for ($i=0; $i<$nn; $i++) {
	for ($j=0; $j<$nn; $j++) {
	    if ($i != $j) {
		for ($s=0; $s<$ns; $s++) {
		    if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			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=0; $k<$nn; $k++) {
	    for ($i=0; $i<$nn; $i++) {
		if ($i != $k) {	# no passengers leave their destination
		    for ($j=0; $j<$nn; $j++) {
			if ($j != $i) {	# no passengers going nowhere
			    # passengers leaving
			    print LPM "- p\_$t\_$i\_$j\_$k ";
			    if ($j != $k) { # no passengers leave their destination
				# passengers arriving
				unless ($tmo == 0) {
				    print LPM "+ p\_$tmo\_$j\_$i\_$k ";
				}
			    }
			} 
		    }
		    # update passenger pools
		    unless ($tmo == 0) {
			print LPM "+ pp\_$tmo\_$i\_$k ";
		    }
		    print LPM "+ ppi\_$tmo\_$i\_$k - pp\_$t\_$i\_$k = 0 ; \n";
		} else {  # pool of people who arrived at their destination $i == $k
		    for ($j=0; $j<$nn; $j++) {
			if ($j != $k) {
			    # number of passengers arriving destination
			    unless ($tmo == 0) {
				print LPM "+ p\_$tmo\_$j\_$k\_$k ";
			    }
			}
		    }
		    # preserve passenger pool who already had arrived at their destination
		    unless ($tmo == 0) { print LPM "+ pp\_$tmo\_$i\_$k "; }
		    print LPM "- pp\_$t\_$i\_$k = 0 ; \n";
		}
	    }
	} 
    } print LPM "\n\n";
    
    
    ## TODO: make this an initialization variable
    print LPM "/* Apply time 0 constraints - no passenger movements */\n";
    $t=0;
    for ($i=0; $i<$nn; $i++) {
	for ($j=0; $j<$nn; $j++) {
	    if ($i != $j) {		# no passengers going nowhere
		for ($k=0; $k<$nn; $k++) {
		    if ($i != $k) {	# no passengers leaving their destination
			print LPM "p\_$t\_$i\_$j\_$k = 0 ; \n";
		    }
		}
	    }
	} 
    } print LPM "\n\n";
    
    
    print LPM "/* Apply integer constraint to flights */\n";
    print LPM "int ";
    my $first = 1;
    for ($t=1; $t<=$nt; $t++) {   # we have no control over $t=0
	for ($i=0; $i<$nn; $i++) {
	    for ($j=0; $j<$nn; $j++) {
		if ($i != $j) {
		    for ($s=0; $s<$ns; $s++) {
			if ($transitgraph[$s][$i][$j] > 0) {  # apply connectivity constraints
			    unless ($first) {	# don't print a comma before first entry
				print LPM ", ";
			    }
			    $first = 0; 
			    print LPM "f\_$t\_$s\_$i\_$j";
			}
		    }
		}
	    }
	} 
    } 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 numveh $dir/sched.ogl >> $dir/results.txt
grep numsrv $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 \"Vehicle pools\" >> $dir/results.txt
grep ap[0-9] $dir/sched.ogl | grep \" [1-9] \" | sort -k 2 >> $dir/results.txt
");
}

return 1;
