Skip to content

Commit

Permalink
Add age/gender/stress interactions to dynamic LCP estimation
Browse files Browse the repository at this point in the history
  • Loading branch information
@CorinStaves
CorinStaves committed Sep 24, 2024
1 parent 4446fed commit 033fcdf
Show file tree
Hide file tree
Showing 7 changed files with 186 additions and 103 deletions.
23 changes: 20 additions & 3 deletions src/main/java/estimation/CoefficientsWriter.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,10 @@ public class CoefficientsWriter {
private final static Logger logger = Logger.getLogger(CoefficientsWriter.class);
private final static String SEP = ",";

static void print(AbstractUtilitySpecification u, BFGS.Results results, double[] se, double[] t, double[] pVal, String[] sig) {
static void print(AbstractUtilitySpecification u, BFGS.Results results, double[] se, double[] t, double[] pVal, String[] sig, String filePath) {

PrintWriter out = ioUtils.openFileForSequentialWriting(new File(filePath),false);
assert out != null;

logger.info("ESTIMATION RESULTS AFTER " + results.iterations + " ITERATIONS:");
Iterator<String> coeffNames = u.getCoeffNames().iterator();
Expand All @@ -24,13 +27,27 @@ static void print(AbstractUtilitySpecification u, BFGS.Results results, double[]
String[] pAll = u.expand(pVal,"%.5f");
String[] sigAll = u.expand(sig);

System.out.printf("| %-30s | %-10s | %-7s | %-10s | %-10s |%n","COEFFICIENT NAME","VALUE","STD.ERR","T.TEST","P.VAL");
String header = String.format("| %-30s | %-10s | %-7s | %-10s | %-10s |%n","COEFFICIENT NAME","VALUE","STD.ERR","T.TEST","P.VAL");
out.print(header);
System.out.print(header);
int i = 0;
while(coeffNames.hasNext()) {
System.out.printf("| %-30s | % .7f | %-7s | %-10s | %-7s %-3s |%n",coeffNames.next(),finalResults[i],seAll[i],tAll[i],pAll[i],sigAll[i]);
String line = String.format("| %-30s | % .7f | %-7s | %-10s | %-7s %-3s |%n",coeffNames.next(),finalResults[i],seAll[i],tAll[i],pAll[i],sigAll[i]);
out.print(line);
System.out.print(line);
i++;
}

// Print key data
out.println();
out.println();
out.println("Iterations: " + results.iterations);
out.println();
out.println("Start LL = " + results.llStart);
out.println("Final LL = " + results.llOut);

out.close();
logger.info("Wrote these results to " + filePath);
}

// Write results to csv file
Expand Down
7 changes: 6 additions & 1 deletion src/main/java/estimation/LogitData.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,11 @@ public int[] getChoices() {
}

public double getValue(int row, String col) {
return this.predictors[row][colIndex.get(col)];
Integer colIdx = colIndex.get(col);
if(colIdx != null) {
return this.predictors[row][colIdx];
} else {
throw new RuntimeException("Column \"" + col + "\" in utility function but not found in database!");
}
}
}
8 changes: 4 additions & 4 deletions src/main/java/estimation/MultinomialLogit.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,11 @@ public static void run(AbstractUtilitySpecification u, int[] y, int k, double la
double[] pVal = pVal(t);
String[] sig = sig(t);

// Print results to screen
CoefficientsWriter.print(u,results,se,t,pVal,sig);
// Print results to screen and text file
CoefficientsWriter.print(u,results,se,t,pVal,sig,resultsFileName + ".txt");

// Print results to file
CoefficientsWriter.write(u,results,se,t,pVal,sig,resultsFileName);
// Print iteration details to csv
CoefficientsWriter.write(u,results,se,t,pVal,sig,resultsFileName + ".csv");
}
}

Expand Down
11 changes: 6 additions & 5 deletions src/main/java/estimation/RunMnlDynamic.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

import estimation.utilities.AbstractUtilitySpecification;
import estimation.utilities.MNL_Dynamic;
import estimation.utilities.MNL_Static;
import gis.GisUtils;
import gis.GpkgReader;
import io.DiaryReader;
Expand Down Expand Up @@ -42,7 +43,7 @@ public static void main(String[] args) throws IOException, FactoryException {
LogitData logitData = new LogitData(args[1],"choice","t.ID");
logitData.read();

// Read Boundary Shapefile
/* // Read Boundary Shapefile
logger.info("Reading boundary shapefile...");
Geometry boundary = GpkgReader.readNetworkBoundary();
Expand Down Expand Up @@ -75,7 +76,7 @@ public static void main(String[] args) throws IOException, FactoryException {
TravelTime ttBike = bicycle.getTravelTimeFast(networkBike);
// Deal with intrazonal trips – can remove after we get X/Y coordinates for TRADS)
Set<SimpleFeature> OAs = GisUtils.readGpkg("zones/2011/gm_oa.gpkg");
Set<SimpleFeature> OAs = GisUtils.readGpkg("zones/2011/gm_oa.gpkg");*/

// Organise classes
int[] y = logitData.getChoices();
Expand All @@ -85,12 +86,12 @@ public static void main(String[] args) throws IOException, FactoryException {
System.out.println("Identified " + k + " classes.");

// Utility function
AbstractUtilitySpecification u = new MNL_Dynamic(logitData,trip_data,OAs,networkBike,bike,ttBike,networkWalk,null,ttWalk);
// AbstractUtilityFunction u = new MNL_Static(logitData);
// AbstractUtilitySpecification u = new MNL_Dynamic(logitData,trip_data,OAs,networkBike,bike,ttBike,networkWalk,null,ttWalk);
AbstractUtilitySpecification u = new MNL_Static(logitData);


// Start model
MultinomialLogit.run(u,y,k,0,1e-10,500,"dynamic6.csv");
MultinomialLogit.run(u,y,k,0,1e-10,500,"estimation/results/static9");

logger.info("finished estimation.");
}
Expand Down
100 changes: 67 additions & 33 deletions src/main/java/estimation/dynamic/DynamicRouter.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
import org.matsim.core.network.NetworkUtils;
import org.matsim.core.router.FastDijkstraFactory;
import org.matsim.core.router.util.LeastCostPathCalculator;
import org.matsim.core.router.util.TravelDisutility;
import org.matsim.core.router.util.TravelTime;
import org.matsim.core.utils.misc.Counter;
import org.matsim.vehicles.Vehicle;
Expand All @@ -36,10 +37,11 @@

public class DynamicRouter implements DynamicUtilityComponent {

private final static boolean ENABLE_DYNAMIC_ROUTING = true;
private final static boolean ENABLE_DYNAMIC_ROUTING = false;
private final static Logger logger = Logger.getLogger(DynamicRouter.class);
private final int numberOfThreads;
final Trip[] trips;
final int[] tripDisutilityType;
final String mode;
final Network network;
final TravelTime tt;
Expand All @@ -48,13 +50,19 @@ public class DynamicRouter implements DynamicUtilityComponent {
final int gammaGradPos;
final int gammaVgviPos;
final int gammaStressLinkPos;
final int gammaStressLinkFemalePos;
final int gammaStressLinkUnder15Pos;
final int gammaStressJctPos;
final int gammaStressJctFemalePos;
final int gammaStressJctUnder15Pos;
final PathData pathData;


public DynamicRouter(Trip[] trips, AbstractUtilitySpecification u, String mode, Set<SimpleFeature> zones,
Network network, Vehicle vehicle, TravelTime tt,
String gammaGrad, String gammaVgvi, String gammaStressLink, String gammaStressJct) {
String gammaGrad, String gammaVgvi,
String gammaStressLink, String gammaStressLinkFemale, String gammaStressLinkUnder15,
String gammaStressJct, String gammaStressJctFemale, String gammaStressJctUnder15) {
this.trips = trips;
this.u = u;
this.mode = mode;
Expand All @@ -64,9 +72,25 @@ public DynamicRouter(Trip[] trips, AbstractUtilitySpecification u, String mode,
this.gammaGradPos = u.getCoeffPos(gammaGrad);
this.gammaVgviPos = u.getCoeffPos(gammaVgvi);
this.gammaStressLinkPos = u.getCoeffPos(gammaStressLink);
this.gammaStressLinkFemalePos = u.getCoeffPos(gammaStressLinkFemale);
this.gammaStressLinkUnder15Pos = u.getCoeffPos(gammaStressLinkUnder15);
this.gammaStressJctPos = u.getCoeffPos(gammaStressJct);
this.gammaStressJctFemalePos = u.getCoeffPos(gammaStressJctFemale);
this.gammaStressJctUnder15Pos = u.getCoeffPos(gammaStressJctUnder15);
this.numberOfThreads = Resources.instance.getInt(Properties.NUMBER_OF_THREADS);

// Set disutility types
tripDisutilityType = new int[trips.length];
for(int i = 0 ; i < trips.length ; i++) {
if(u.value(i,"p.age_group_agg_5_14") == 1) {
tripDisutilityType[i] = 2; // CHILD UNDER 15
} else if(u.value(i,"p.female") == 1) {
tripDisutilityType[i] = 1; // FEMALE
} else {
tripDisutilityType[i] = 0; // MALE
}
}

// Sort inter/intra-zonal trips, and compute fixed results for intrazonal
pathData = new PathData(trips.length);
computeIntrazonalTripData(pathData,trips,u,zones,network,mode);
Expand Down Expand Up @@ -130,38 +154,37 @@ public void update(double[] xVarOnly) {

private void updatePathData(double[] xVarOnly) {

// Get latest coefficients
// GET LATEST COEFFICIENTS
double[] x = u.expandCoeffs(xVarOnly);
double mGrad = x[gammaGradPos];
double mVgvi = x[gammaVgviPos];
double mStressLink = x[gammaStressLinkPos];
double mStressJct = x[gammaStressJctPos];

if(mGrad < 0) {
logger.warn("Negative gradient coefficient: " + mGrad + " Setting mGrad=0 for routing...");
mGrad = 0;
}
if(mVgvi < 0) {
logger.warn("Negative vgvi coefficient: " + mVgvi + ". Set mVgvi=0 for routing...");
mVgvi = 0;
}
if(mStressLink < 0) {
logger.warn("Negative link stress coefficient: " + mStressLink + ". Set mStressLink=0 for routing...");
mStressLink = 0;
}
if(mStressJct < 0) {
logger.warn("Negative junction stress coefficient: " + mStressJct + ". Set mStressJct=0 for routing...");
mStressJct = 0;
}

// Gradient & VGVI
double mGrad = zeroIfNegative(x[gammaGradPos],"gradient");
double mVgvi = zeroIfNegative(x[gammaVgviPos],"VGVI");

// Link stress + adjustments
double mStressLinkMale = zeroIfNegative(x[gammaStressLinkPos],"link stress (baseline)");
double mStressLinkFemale = zeroIfNegative(x[gammaStressLinkPos] + x[gammaStressLinkFemalePos],"link stress (female)");
double mStressLinkUnder15 = zeroIfNegative(x[gammaStressLinkPos] + x[gammaStressLinkUnder15Pos],"link stress (under 15)");

// Junction stress + adjustments
double mStressJctMale = zeroIfNegative(x[gammaStressJctPos],"junction stress (baseline)");
double mStressJctFemale = zeroIfNegative(x[gammaStressJctPos] + x[gammaStressJctFemalePos],"junction stress (female)");
double mStressJctUnder15 = zeroIfNegative(x[gammaStressJctPos] + x[gammaStressJctUnder15Pos],"junction stress (under 15)");

// Update all routable trips (multithreaded)
JibeDisutility2 disutility = new JibeDisutility2(network, vehicle, mode, tt, mGrad, mVgvi, mStressLink, mStressJct);
ConcurrentLinkedQueue<Integer> tripsQueue = IntStream.range(0, trips.length).boxed().collect(Collectors.toCollection(ConcurrentLinkedQueue::new));
TravelDisutility disutilityMale = new JibeDisutility2(network, vehicle, mode, tt, mGrad, mVgvi, mStressLinkMale, mStressJctMale);
TravelDisutility disutilityFemale = new JibeDisutility2(network, vehicle, mode, tt, mGrad, mVgvi, mStressLinkFemale, mStressJctFemale);
TravelDisutility disutilityUnder15 = new JibeDisutility2(network, vehicle, mode, tt, mGrad, mVgvi, mStressLinkUnder15, mStressJctUnder15);

Thread[] threads = new Thread[numberOfThreads];
Counter counter = new Counter("Routed ", " / " + trips.length + " trips.");
ConcurrentLinkedQueue<Integer> tripsQueue = IntStream.range(0, trips.length).boxed().collect(Collectors.toCollection(ConcurrentLinkedQueue::new));
for(int i = 0 ; i < numberOfThreads ; i++) {
LeastCostPathCalculator dijkstra = new FastDijkstraFactory(false).createPathCalculator(network, disutility, tt);
TripWorker worker = new TripWorker(tripsQueue,counter,mode,tt,vehicle,dijkstra,pathData);
LeastCostPathCalculator[] lcpCalculators = new LeastCostPathCalculator[3];
lcpCalculators[0] = new FastDijkstraFactory(false).createPathCalculator(network, disutilityMale, tt);
lcpCalculators[1] = new FastDijkstraFactory(false).createPathCalculator(network, disutilityFemale, tt);
lcpCalculators[2] = new FastDijkstraFactory(false).createPathCalculator(network, disutilityUnder15, tt);
TripWorker worker = new TripWorker(tripsQueue,counter,mode,tt,vehicle,lcpCalculators,tripDisutilityType,pathData);
threads[i] = new Thread(worker,"DynamicRouteUpdate-" + i);
threads[i].start();
}
Expand All @@ -176,6 +199,15 @@ private void updatePathData(double[] xVarOnly) {
}
}

private static double zeroIfNegative(double n, String name) {
if(n < 0) {
logger.warn("Negative " + name + " coefficient:" + n + ". Set to 0 for routing...");
return 0.;
} else {
return n;
}
}

public static class PathData {
final int length;
final boolean[] intrazonal;
Expand Down Expand Up @@ -205,20 +237,22 @@ private static class TripWorker implements Runnable {
private final ConcurrentLinkedQueue<Integer> tripIndices;
private final Counter counter;
private final Vehicle vehicle;
private final LeastCostPathCalculator pathCalculator;
private final LeastCostPathCalculator[] lcpCalculators;
private final int[] lcpCalculatorIdx;
private final TravelTime tt;
private final String mode;
private final PathData pathData;


public TripWorker(ConcurrentLinkedQueue<Integer> tripIndices, Counter counter, String mode,
TravelTime travelTime, Vehicle vehicle, LeastCostPathCalculator pathCalculator,PathData pathData) {
public TripWorker(ConcurrentLinkedQueue<Integer> tripIndices, Counter counter, String mode, TravelTime travelTime, Vehicle vehicle,
LeastCostPathCalculator[] lcpCalculators, int[] lcpCalculatorIdx, PathData pathData) {
this.tripIndices = tripIndices;
this.counter = counter;
this.mode = mode;
this.tt = travelTime;
this.vehicle = vehicle;
this.pathCalculator = pathCalculator;
this.lcpCalculators = lcpCalculators;
this.lcpCalculatorIdx = lcpCalculatorIdx;
this.pathData = pathData;
}

Expand All @@ -231,7 +265,7 @@ public void run() {
}
//this.counter.incCounter();
if(!pathData.intrazonal[i]) {
LeastCostPathCalculator.Path path = pathCalculator.calcLeastCostPath(pathData.originNodes[i], pathData.destinationNodes[i], 0., null, vehicle);
LeastCostPathCalculator.Path path = lcpCalculators[lcpCalculatorIdx[i]].calcLeastCostPath(pathData.originNodes[i], pathData.destinationNodes[i], 0., null, vehicle);
double tripGradient = 0.;
double tripVgvi = 0.;
double tripStressLink = 0.;
Expand Down
Loading

0 comments on commit 033fcdf

Please sign in to comment.