Mode choice estimation and skims updates

src/main/java/diary/RunRouter.java

@@ -45,7 +45,7 @@ public class RunRouter {
     public static void main(String[] args) throws IOException, FactoryException {
 
         if(args.length != 3) {
-            throw new RuntimeException("Program requires 2 arguments: \n" +
+            throw new RuntimeException("Program requires 3 arguments: \n" +
                     "(0) Properties file \n" +
                     "(1) Route data output (.csv) \n" +
                     "(2) Link data output (specify .gpkg OR write \"true\" to include links in csv file. Otherwise write \"false\") \n");

src/main/java/estimation/CoefficientsWriter.java

@@ -14,7 +14,7 @@ public class CoefficientsWriter {
     private final static Logger logger = Logger.getLogger(CoefficientsWriter.class);
     private final static String SEP = ",";
 
-    static void print(AbstractModelSpecification u, BFGS.Results results, double[] se, double[] t, double[] pVal, String[] sig, String filePath) {
+    static void print(AbstractModelSpecification u, BFGS.Results results,double ll0, double[] se, double[] t, double[] pVal, String[] sig, String filePath) {
 
         PrintWriter out = ioUtils.openFileForSequentialWriting(new File(filePath),false);
         assert out != null;
@@ -43,6 +43,7 @@ static void print(AbstractModelSpecification u, BFGS.Results results, double[] s
         out.println();
         out.println("Iterations: " + results.iterations);
         out.println();
+        out.println("LL0 = " + ll0);
         out.println("Start LL = " + results.llStart);
         out.println("Final LL = " + results.llOut);
 

src/main/java/estimation/MultinomialLogit.java

@@ -26,15 +26,18 @@ public static void run(AbstractModelSpecification u, int[] y, int k, double lamb
             throw new IllegalArgumentException("Invalid maximum number of iterations: " + maxIter);
         }
 
-        MultinomialLogitObjective objective = new MultinomialLogitObjective(u, y, k, lambda);
-
-        // Run the optimisation algorithm
+        // Get starting values
         double[] w = u.getStarting();
-        BFGS.Results results = BFGS.minimize(objective,u.getDynamicComponent(), w, tol, maxIter);
-        logger.info("finished estimation.");
 
-        // Get overlap with fastest path
+        // Create test objective just to check LL0
+        MultinomialLogitObjective test = new MultinomialLogitObjective(u, y, k, lambda);
+        double ll0 = test.f(new double[w.length]);
+        logger.info("LL0 = " + ll0);
 
+        // Run BFGS algrorithm
+        MultinomialLogitObjective objective = new MultinomialLogitObjective(u, y, k, lambda);
+        BFGS.Results results = BFGS.minimize(objective,u.getDynamicComponent(), w, tol, maxIter);
+        logger.info("finished estimation.");
 
         // Approximate variance-coviariance matrix (from BFGS method) – for debugging only
         // Matrix approxVarCov = Matrix.of(results.hessian);
@@ -49,7 +52,7 @@ public static void run(AbstractModelSpecification u, int[] y, int k, double lamb
         // Check if matrix is singular
         if(Arrays.stream(eigenvalues.wr).anyMatch(e -> e == 0)) {
             logger.error("Hessian matrix is singular! Cannot compute standard errors");
-            CoefficientsWriter.print(u,results,null,null,null,null,resultsFileName + ".txt");
+            CoefficientsWriter.print(u,results,ll0,null,null,null,null,resultsFileName + ".txt");
             CoefficientsWriter.write(u,results,null,null,null,null,resultsFileName);
         } else {
 
@@ -69,7 +72,7 @@ public static void run(AbstractModelSpecification u, int[] y, int k, double lamb
             String[] sig = sig(t);
 
             // Print results to screen and text file
-            CoefficientsWriter.print(u,results,se,t,pVal,sig,resultsFileName + ".txt");
+            CoefficientsWriter.print(u,results,ll0,se,t,pVal,sig,resultsFileName + ".txt");
 
             // Print iteration details to csv
             CoefficientsWriter.write(u,results,se,t,pVal,sig,resultsFileName + ".csv");

src/main/java/estimation/RunMnlManchester.java

@@ -1,10 +1,7 @@
 package estimation;
 
 import estimation.specifications.AbstractModelSpecification;
-import estimation.specifications.manchester.HBA;
-import estimation.specifications.manchester.HBD;
-import estimation.specifications.manchester.NHBO;
-import estimation.specifications.manchester.NHBW;
+import estimation.specifications.manchester.*;
 import gis.GisUtils;
 import gis.GpkgReader;
 import io.DiaryReader;
@@ -21,7 +18,6 @@
 import routing.Bicycle;
 import routing.travelTime.WalkTravelTime;
 import smile.classification.ClassLabels;
-import trip.Place;
 import trip.Trip;
 
 import java.io.IOException;

src/main/java/estimation/dynamic/RouteDataDynamic.java

@@ -285,11 +285,6 @@ private void computeRouteStats() {
                 .sorted().toArray();
         overlapStats = mode.toUpperCase() + " OVERLAP: " + getSummaryStats(overlap);
         logger.info(overlapStats);
-
-        // Runtime (for debugging only)
-        long endTime= System.currentTimeMillis();
-        double runtime = (endTime - startTime) / 1000.;
-        logger.info(mode.toUpperCase() + " Route stats computation runtime (s) = " + runtime);
     }
 
     public String getStats() {

src/main/java/estimation/specifications/manchester/HBD.java

@@ -42,9 +42,9 @@ List<String> sociodemographicVariables() {
     List<RouteAttribute> streetEnvironmentAttributesBike() {
         List<RouteAttribute> bike = new ArrayList<>();
         bike.add(new RouteAttribute("g_bike_grad", l -> Math.max(Math.min(Gradient.getGradient(l),0.5),0.)));
-        bike.add(new RouteAttribute("g_bike_vgvi", l -> Math.max(0.,0.81 - LinkAmbience.getVgviFactor(l))));
+//        bike.add(new RouteAttribute("g_bike_vgvi", l -> Math.max(0.,0.81 - LinkAmbience.getVgviFactor(l))));
         bike.add(new RouteAttribute("g_bike_stressLink", l -> LinkStress.getStress(l, TransportMode.bike)));
-        bike.add(new RouteAttribute("g_bike_stressLink_f","g_bike_stressLink", i -> value(i,"p.female") == 1));
+//        bike.add(new RouteAttribute("g_bike_stressLink_f","g_bike_stressLink", i -> value(i,"p.female") == 1));
         bike.add(new RouteAttribute("g_bike_stressLink_c","g_bike_stressLink", i -> value(i,"p.age_group_agg_5_14") == 1));
 //        bike.add(new RouteAttribute("g_bike_stressLink_o","g_bike_stressLink",i -> value(i,"p.age_65up") == 1));
         return bike;
@@ -54,13 +54,15 @@ List<RouteAttribute> streetEnvironmentAttributesWalk() {
         List<RouteAttribute> walk = new ArrayList<>();
 //        bike.add(new BuiltEnvironmentAttribute("grad",l -> Math.max(Math.min(Gradient.getGradient(l),0.5),0.),p -> true));
         walk.add(new RouteAttribute("g_walk_vgvi", l -> Math.max(0.,0.81 - LinkAmbience.getVgviFactor(l))));
+//        walk.add(new RouteAttribute("g_walk_vgvi_c", "g_walk_vgvi", i -> value(i,"p.age_group_agg_5_14") == 1));
+//        walk.add(new RouteAttribute("g_walk_vgvi_o", "g_walk_vgvi", i -> value(i,"p.age_65up") == 1));
         walk.add(new RouteAttribute("g_walk_stressJct", l -> JctStress.getStressProp(l,TransportMode.walk)));
 //        walk.add(new RouteAttribute("g_walk_stressJct_c", "g_walk_stressJct", i -> value(i,"p.age_group_agg_5_14") == 1));
 //        walk.add(new RouteAttribute("g_walk_stressJct_o", "g_walk_stressJct", i -> value(i,"p.age_65up") == 1));
         walk.add(new RouteAttribute("g_walk_speed", l -> Math.min(1.,l.getFreespeed() / 22.35)));
         walk.add(new RouteAttribute("g_walk_speed_c", "g_walk_speed", i -> value(i,"p.age_group_agg_5_14") == 1));
         walk.add(new RouteAttribute("g_walk_speed_o", "g_walk_speed", i -> value(i,"p.age_65up") == 1));
-        walk.add(new RouteAttribute("g_walk_aadt", l -> Math.min(1.,((int) l.getAttributes().getAttribute("aadt")) * 0.865/6000.)));
+//        walk.add(new RouteAttribute("g_walk_aadt", l -> Math.min(1.,((int) l.getAttributes().getAttribute("aadt")) * 0.865/6000.)));
 //        walk.add(new RouteAttribute("g_walk_aadt_c", "g_walk_aadt", i -> value(i,"p.age_group_agg_5_14") == 1));
 //        walk.add(new RouteAttribute("g_walk_aadt_o", "g_walk_aadt", i -> value(i,"p.age_65up") == 1));
 

src/main/java/estimation/specifications/manchester/HBW.java

@@ -19,7 +19,7 @@
 public class HBW extends MNL_Manchester {
 
     private final static List<String> SOCIODEMOGRAPHIC_VARIABLES = List.of(
-            "p.age_group_agg_15_24","p.age_group_agg_40_54","p.age_group_agg_55",
+            "p.age_group_agg_15_24","p.age_group_agg_40_54","p.age_55up",
             "p.female",
             "hh.cars_gr_0","hh.cars_gr_2","hh.cars_gr_3",
             "hh.income_agg_high");
@@ -41,7 +41,6 @@ List<RouteAttribute> streetEnvironmentAttributesBike() {
 //        bike.add(new RouteAttribute("g_bike_vgvi", l -> Math.max(0.,0.81 - LinkAmbience.getVgviFactor(l))));
         bike.add(new RouteAttribute("g_bike_stressLink", l -> LinkStress.getStress(l,TransportMode.bike)));
         bike.add(new RouteAttribute("g_bike_stressLink_f","g_bike_stressLink", i -> value(i,"p.female") == 1));
-//        bike.add(new BuiltEnvironmentAttribute("stressLink_c",l -> LinkStress.getStress(l,TransportMode.bike),p -> (boolean) p.getAttributes().getAttribute("p.under15")));
         return bike;
     }
 
@@ -60,9 +59,9 @@ List<RouteAttribute> streetEnvironmentAttributesWalk() {
     protected List<String> fixed() {
         List<String> fixed = coefficients().keySet().stream().filter(s -> (s.contains("carD"))).collect(Collectors.toList());
 //        fixed.add("b_carP_hh.income_agg_low");
-        fixed.add("b_carP_p.age_group_agg_55");
+        fixed.add("b_carP_p.age_55up");
         fixed.add("b_walk_p.age_group_agg_40_54");
-        fixed.add("b_walk_p.age_group_agg_55");
+        fixed.add("b_walk_p.age_55up");
 //        fixed.add("b_carP_t.full_purpose_HBO");
 //        fixed.add("b_bike_p.age_group_agg_40_69");
 //        fixed.add("b_bike_hh.income_agg_high");

src/main/java/estimation/specifications/manchester/MNL_Manchester.java

@@ -40,7 +40,6 @@ public MNL_Manchester(LogitData data, Trip[] trips, Set<SimpleFeature> OAs,
             bikeRouteData = new RouteDataDynamic(data.getIds(), trips,this, TransportMode.bike,OAs,netBike,vehBike,ttBike,attributesBike);
             walkRouteData = new RouteDataDynamic(data.getIds(), trips,this, TransportMode.walk,OAs,netWalk,vehWalk,ttWalk,attributesWalk);
         }
-//        System.exit(-1);
         initialiseDynamicUtilDeriv();
     }
 

src/main/java/network/CreateMatsimNetworkRoad.java

@@ -63,11 +63,11 @@ public static void main(String[] args) {
         // Create network links
         edges.forEach((id,edge) -> addLinkToNetwork(id,edge,net,fac));
 
-//        // Add volumes from events file
-//        NetworkUtils2.addSimulationVolumes(readSimulationVolumes(),net);
-//
-//        // Write crossing attributes
-//        NetworkUtils2.addCrossingAttributes(net);
+        // Add volumes from events file
+        NetworkUtils2.addSimulationVolumes(readSimulationVolumes(),net);
+
+        // Write crossing attributes
+        NetworkUtils2.addCrossingAttributes(net);
 
         // Identify disconnected links
         NetworkUtils2.identifyDisconnectedLinks(net,walk);

src/main/java/skim/RunSkimsPt.java

@@ -36,13 +36,14 @@ public static void main(String[] args) throws IOException, FactoryException {
 
         CalculateSkimMatrices skims = new CalculateSkimMatrices("/",numberOfThreads);
         skims.loadSamplingPointsFromFile(zonesFilename);
-        PTSkimMatrices.PtIndicators<String> ptIndicators = skims.calculatePTMatrices(transitNetworkFilePath,transitScheduleFilePath,28800,30600,config,(a,b) -> false);
+        PTSkimMatrices.PtIndicators<String> ptIndicators = skims.calculatePTMatrices(transitNetworkFilePath,transitScheduleFilePath,28800,30600,config,(a,b) -> b.getTransportMode().equals("bus"));
 
         // Choose which matrices to write
         Map<String, FloatMatrix<String>> matricesToWrite = new LinkedHashMap<>();
         matricesToWrite.put("travelTime",ptIndicators.travelTimeMatrix);
         matricesToWrite.put("accessTime",ptIndicators.accessTimeMatrix);
         matricesToWrite.put("egressTime",ptIndicators.egressTimeMatrix);
+        matricesToWrite.put("busTimeShare",ptIndicators.trainTravelTimeShareMatrix);
 
         // Write matrices
         OmxWriter.createOmxSkimMatrix(outputFile,skims.getCoordsPerZone().keySet(),matricesToWrite);