Refactor cost calculations for specific classes by directly passing in class during cost

.gitignore

@@ -20,4 +20,4 @@ todo-nosync.txt
 *.bin
 valgrind.txt
 scratch/
-net/
+net/*

README.md

@@ -68,7 +68,7 @@ The key points are as follows:
 - The `master` branch will always point to a stable release, each of which will have an associated version number as tag, and a set of release notes.
 - If you are working on a new feature, **always** branch off of `develop` (**not** `master`).  The `develop` branch is used for features under development, experiments, or other changes.  Each branch should have a `branch.txt` file in the root directory briefly stating the purpose of the branch.
 - The only branches off of `master` will be to fix significant bugs as a `hotfix`.
-- When a new feature is complete, merge it back into `develop` (**not** `master`).  
+- When a new feature is complete, merge it back into `develop` (**not** `master`).
 
 This transition is currently in progress, identifying active branches which need to be merged, and other branches still with potential for development.
 When complete, the `master` branch will point to the v1.0 release, and other incomplete branches will be migrated to `develop`.

branch.txt

@@ -1 +1,4 @@
-This branch has code to write *path flows* corresponding to a bush-based solution.  No attempt is made to maximize entropy, just divide flow proportionally along bush paths.
+class-cost-refactor
+
+This branch aims to simplify calculations for the multiclass version of TAP.  By passing cost as an argument to the cost functions, this streamlines code (avoids loops through classes) and may possibly have performance and multithreading benefits.  Testing still underway.
+

include/bush.h

@@ -261,7 +261,7 @@ typedef struct algorithmBParameters_type{
                                  struct algorithmBParameters_type *);
    void     (*topologicalOrder)(int, network_type *, bushes_type *,
                                 struct algorithmBParameters_type *);
-   void     (*linkShiftB)(int, double, network_type *);
+   void     (*linkShiftB)(int, double, network_type *, int class);
 
 #if PARALLELISM
    int numThreads;
@@ -288,10 +288,8 @@ algorithmBParameters_type initializeAlgorithmBParameters();
 /* Main Algorithm B helper functions */
 void initializeAlgorithmB(network_type *network, bushes_type **bushes,
                           algorithmBParameters_type *parameters);
-void updateBatchBushes(network_type *network, bushes_type *bushes,
-                       int *lastClass, algorithmBParameters_type *parameters);
-void updateBatchFlows(network_type *network, bushes_type *bushes,
-                      int *lastClass, algorithmBParameters_type *parameters);
+void updateBatchBushes(network_type *network, bushes_type *bushes, algorithmBParameters_type *parameters);
+void updateBatchFlows(network_type *network, bushes_type *bushes, algorithmBParameters_type *parameters);
 void loadBatch(int batch, network_type *network, bushes_type **bushes,
                algorithmBParameters_type *parameters);
 void storeBatch(int batch, network_type *network, bushes_type *bushes,
@@ -360,7 +358,7 @@ bool isInBush(int origin, int ij, network_type *network, bushes_type *bushes);
 bool isMergeNode(int origin, int i, bushes_type *bushes);
 int pred2merge(int ij);
 int merge2pred(int m);
-void exactCostUpdate(int ij, double shift, network_type *network);
+void exactCostUpdate(int ij, double shift, network_type *network, int class);
 void linearCostUpdate(int ij, double shift, network_type *network);
 void noCostUpdate(int ij, double shift, network_type *network);
 void checkFlows(network_type *network, bushes_type *bushes);

include/convexcombination.h

@@ -86,7 +86,7 @@ void updateLinks(network_type *network, bool rectifyLinks);
 
 double CCrelativeGap(network_type *network, double tstt, double sptt);
 double CCaverageExcessCost(network_type *network, double tstt, double sptt);
-double evaluateLinkCost(arc_type *arc, double flow);
+double evaluateLinkCost(arc_type *arc, double flow, int class);
 
 double calculateConjugacy(network_type *network, double **direction1,
                           double **direction2);

include/networks.h

@@ -84,9 +84,9 @@ typedef struct arc_type {
     int     linkType;
 
     double  fixedCost; /* Reflects toll and distance */
-    double  (*calculateCost)(struct arc_type *arc);
+    double  (*calculateCost)(struct arc_type *arc, int);
     double  (*calculateDer)(struct arc_type *der);
-    double  (*calculateInt)(struct arc_type *der, bool);
+    double  (*calculateInt)(struct arc_type *der, int);
 } arc_type;
 
 

include/tap.h

@@ -43,18 +43,18 @@ double TSTT(network_type *network);
 void updateAllCosts(network_type *network);
 void updateAllCostDers(network_type *network);
 
-double generalBPRcost(struct arc_type *arc);
+double generalBPRcost(struct arc_type *arc, int c);
 double generalBPRder(struct arc_type *arc);
-double generalBPRint(struct arc_type *arc, bool includeFixedCost);
-double linearBPRcost(struct arc_type *arc);
+double generalBPRint(struct arc_type *arc, int c);
+double linearBPRcost(struct arc_type *arc, int c);
 double linearBPRder(struct arc_type *arc);
-double linearBPRint(struct arc_type *arc, bool includeFixedCost);
-double quarticBPRcost(struct arc_type *arc);
+double linearBPRint(struct arc_type *arc, int c);
+double quarticBPRcost(struct arc_type *arc, int c);
 double quarticBPRder(struct arc_type *arc);
-double quarticBPRint(struct arc_type *arc, bool includeFixedCost);
-double conicCost(struct arc_type *arc);
+double quarticBPRint(struct arc_type *arc, int c);
+double conicCost(struct arc_type *arc, int c);
 double conicDer(struct arc_type *arc);
-double conicInt(struct arc_type *arc, bool includeFixedCost);
+double conicInt(struct arc_type *arc, int c);
 
 
 int arcNumber(network_type *network, arc_type *arc);

include/utils.h

@@ -15,6 +15,7 @@
 #define EXIT_DEBUG -2
 
 #define IS_MISSING -1
+#define NO_FIXED_COST -1
 #define STRING_SIZE 20000
 
 #define PAUSE_ON_ERROR FALSE

src/bush.c

@@ -253,7 +253,7 @@ void initializeAlgorithmB(network_type *network, bushes_type **bushes,
                 network->arcs[ij].flow += (*bushes)->flow[ij];
                 network->arcs[ij].classFlow[c] += (*bushes)->flow[ij];
                 network->arcs[ij].cost =
-                    network->arcs[ij].calculateCost(&network->arcs[ij]);
+                    network->arcs[ij].calculateCost(&network->arcs[ij], NO_FIXED_COST);
             }
         }
         snprintf(batchFileName, 2*STRING_SIZE, "%s%d.bin",
@@ -297,9 +297,8 @@ void storeBatch(int batch, network_type *network, bushes_type *bushes,
 }
 
 void updateBatchBushes(network_type *network, bushes_type *bushes,
-                       int *lastClass, algorithmBParameters_type *parameters) {
+                       algorithmBParameters_type *parameters) {
 #if PARALLELISM
-    int c; 
     struct thread_args args[network->batchSize];
     for (int j = 0; j < network->batchSize; ++j) {
         args[j].id = j;
@@ -312,47 +311,32 @@ void updateBatchBushes(network_type *network, bushes_type *bushes,
     for (int j = 0; j < network->batchSize; ++j) {
         if (outOfOrigins(network, j) == TRUE) break;
         bushes->updateBush[j] = TRUE;
-        c = origin2class(network, j);
-        if (c != *lastClass) {
-            changeFixedCosts(network, c);
-        }
         thpool_add_work(thpool, (void (*)(void *)) updateBushPool,
                         (void*)&args[j]);
-        *lastClass = c;
     }
     thpool_wait(thpool);
 
     for (int j = 0; j < network->batchSize; ++j) {
         if (outOfOrigins(network, j) == TRUE) break;
         bushes->updateBush[j] = TRUE;
-        c = origin2class(network, j);
-        if (c != *lastClass) {
-            changeFixedCosts(network, c);
-        }
-        thpool_add_work(thpool, (void (*)(void *)) updateFlowsPool, 
+        thpool_add_work(thpool, (void (*)(void *)) updateFlowsPool,
                         (void*)&args[j]);
-        *lastClass = c;
     }
     thpool_wait(thpool);
 #else
-    int origin, c;
+    int origin;
     for (origin = 0; origin < network->batchSize; origin++) {
         if (outOfOrigins(network, origin) == TRUE) break;
         bushes->updateBush[origin] = TRUE;
-        c = origin2class(network, origin);
-        if (c != *lastClass) {
-            changeFixedCosts(network, c);
-        }
         updateBushB(origin, network, bushes, parameters);
         updateFlowsB(origin, network, bushes, parameters);
-        *lastClass = c;
     }
 #endif
 }
 
 void updateBatchFlows(network_type *network, bushes_type *bushes,
-                      int *lastClass, algorithmBParameters_type *parameters) {
-    int i, c;                       
+                      algorithmBParameters_type *parameters) {
+    int i;
     bool doneAny;
     for (i = 0; i < parameters->innerIterations; i++) {
         doneAny = FALSE;
@@ -369,10 +353,6 @@ void updateBatchFlows(network_type *network, bushes_type *bushes,
         for (int j = 0; j < network->batchSize; ++j) {
             if (outOfOrigins(network, j) == TRUE) break;
             if (bushes->updateBush[j] == FALSE) continue;
-            c = origin2class(network, j);
-            if (c != *lastClass) {
-                changeFixedCosts(network, c);
-            }
             thpool_add_work(thpool, (void (*)(void *)) updateFlowsPool,
                             (void*)&args[j]);
         }
@@ -386,12 +366,7 @@ void updateBatchFlows(network_type *network, bushes_type *bushes,
         for (origin = 0; origin < network->batchSize; origin++) {
             if (outOfOrigins(network, origin) == TRUE) break;
             if (bushes->updateBush[origin] == FALSE) continue;
-            c = origin2class(network, origin);
-            if (c != *lastClass) {
-                changeFixedCosts(network, c);
-            }
             doneAny |= updateFlowsB(origin,network,bushes,parameters);
-            *lastClass = c;
         }
 #endif
         if (doneAny == FALSE) break;
@@ -596,7 +571,7 @@ bool isInBush(int origin, int ij, network_type *network, bushes_type *bushes) {
  */
 double bushSPTT(network_type *network, bushes_type *bushes,
               algorithmBParameters_type *parameters) {
-    int b, r, ij, i, j, c, lastClass = IS_MISSING, originNode;
+    int b, r, ij, i, j, c, originNode;
     double frac, rc, acceptanceGap = 0, rejectionGap = INFINITY, consistency;
     double sptt = 0;
     if (parameters->calculateBins == TRUE) {
@@ -608,10 +583,6 @@ double bushSPTT(network_type *network, bushes_type *bushes,
     for (r = 0; r < network->batchSize; r++) {
         if (outOfOrigins(network, r) == TRUE) break;
         originNode = origin2node(network, r);
-        c = origin2class(network, r);
-        if (c != lastClass) {
-            changeFixedCosts(network, c);
-        }
         scanBushes(r, network, bushes, parameters, NO_LONGEST_PATH);
         BellmanFord_NoLabel(originNode, bushes->SPcost, network, DEQUE,
                             bushes->SPcost, bushes->bushOrder[r]);
@@ -692,19 +663,14 @@ double bushSPTT(network_type *network, bushes_type *bushes,
 }
 
 double bushTSTT(network_type *network, bushes_type *bushes) {
-    int r, ij, c, lastClass = IS_MISSING;
+    int r, ij;
     double tstt = 0;
     for (r = 0; r < network->batchSize; r++) {
         if (outOfOrigins(network, r) == TRUE) break;
-        c = origin2class(network, r);
-        if (c != lastClass) {
-            changeFixedCosts(network, c);
-        }
         calculateBushFlows(r, network, bushes);
         for (ij = 0; ij < network->numArcs; ij++) {
             tstt += bushes->flow[ij] * network->arcs[ij].cost;
         }
-        lastClass = c;
     }
     return tstt;
 }
@@ -763,14 +729,10 @@ double bushAEC(network_type *network, bushes_type *bushes,
 double bushMEC(network_type *network, bushes_type *bushes,
              algorithmBParameters_type *parameters) {
     double mec = 0;
-    int j, r, c, lastClass = IS_MISSING, originNode;
+    int j, r, originNode;
     for (r = 0; r < network->batchSize; r++) {
         if (outOfOrigins(network, r) == TRUE) break;
         originNode = origin2node(network, r);
-        c = origin2class(network, r);
-        if (c != lastClass) {
-            changeFixedCosts(network, c);
-        }
         scanBushes(r, network, bushes, parameters, LONGEST_USED_PATH);
         BellmanFord_NoLabel(originNode, bushes->SPcost, network, DEQUE,
                             bushes->SPcost, bushes->bushOrder[r]);
@@ -865,15 +827,10 @@ void deleteBushes(network_type *network, bushes_type *bushes) {
 */
 void initializeBushesB(network_type *network, bushes_type *bushes,
                        algorithmBParameters_type *parameters) {
-    int c, origin, lastClass = IS_MISSING;
+    int origin;
 
     for (origin = 0; origin < network->batchSize; origin++) {
         if (outOfOrigins(network, origin) == TRUE) break;
-        c = origin2class(network, origin);
-        if (c != lastClass) {
-            changeFixedCosts(network, c);
-        }
-        lastClass = c;
         /* createInitialBush also sets preds, bushOrder */
         parameters->createInitialBush(origin, network, bushes, parameters);
         calculateBushFlows(origin, network, bushes);
@@ -899,6 +856,8 @@ void scanBushes(int origin, network_type *network, bushes_type *bushes,
     }
 
     /* Ensure costs are up to date */
+    // TODO: updateAllCosts in the non-exactCostUpdate will not work in
+    // parallel with the current fixed cost mechanism
     if (parameters->linkShiftB != &exactCostUpdate) updateAllCosts(network);
 
     bushes->SPcost[origin2node(network, origin)] = 0;
@@ -1470,7 +1429,7 @@ void newtonFlowShift(int j, merge_type *merge, int origin,
         }
         bushes->flow[hi] -= shift;
         network->arcs[hi].classFlow[c] -= shift;
-        parameters->linkShiftB(hi, -shift, network);
+        parameters->linkShiftB(hi, -shift, network, c);
         i = network->arcs[hi].tail;      
     }
     i = j;
@@ -1485,7 +1444,7 @@ void newtonFlowShift(int j, merge_type *merge, int origin,
         }
         bushes->flow[hi] += shift;
         network->arcs[hi].classFlow[c] += shift;
-        parameters->linkShiftB(hi, shift, network);
+        parameters->linkShiftB(hi, shift, network, c);
         i = network->arcs[hi].tail;      
     }
 }
@@ -1666,9 +1625,9 @@ void printBush(int minVerbosity, int origin, network_type *network,
  * changing its flow, by explicitly recomputing the BPR function and its
  * derivative.
  */
-void exactCostUpdate(int ij, double shift, network_type *network) {
+void exactCostUpdate(int ij, double shift, network_type *network, int c) {
     network->arcs[ij].flow += shift;
-    network->arcs[ij].cost=network->arcs[ij].calculateCost(&network->arcs[ij]);
+    network->arcs[ij].cost=network->arcs[ij].calculateCost(&network->arcs[ij], c);
     network->arcs[ij].der = network->arcs[ij].calculateDer(&network->arcs[ij]);
 }