diff --git a/.gitignore b/.gitignore
index acdf0d2..babe4ca 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,4 +20,4 @@ todo-nosync.txt
 *.bin
 valgrind.txt
 scratch/
-net/
+net/*
diff --git a/README.md b/README.md
index e504c00..f8dbea9 100644
--- a/README.md
+++ b/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`.
diff --git a/branch-info.txt b/branch-info.txt
deleted file mode 100644
index b9d7ce6..0000000
--- a/branch-info.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-This branch adds the option to compute reduced-cost information on all network
-links (in-bush and out-of-bush) at each master iteration.  These values can be
-used to compute acceptance and rejection gaps for checking path flow
-consistency, or adjusting numerical precision on-the-fly.
-
-The first implementation and debugging will be on the serial build (runs are
-deterministic, so debugging is easier), before being ported to parallel.
diff --git a/branch.txt b/branch.txt
index 1e73779..d007f0a 100644
--- a/branch.txt
+++ b/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.
+
diff --git a/include/bush.h b/include/bush.h
index d9276dd..b954c02 100644
--- a/include/bush.h
+++ b/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);
diff --git a/include/convexcombination.h b/include/convexcombination.h
index 97f9e99..675e2ba 100644
--- a/include/convexcombination.h
+++ b/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);
diff --git a/include/networks.h b/include/networks.h
index 2490ae7..3802c4f 100644
--- a/include/networks.h
+++ b/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;
 
 
diff --git a/include/tap.h b/include/tap.h
index 1123c7c..d37559a 100644
--- a/include/tap.h
+++ b/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);
diff --git a/include/utils.h b/include/utils.h
index a7addb7..21bbca9 100644
--- a/include/utils.h
+++ b/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
diff --git a/src/bush.c b/src/bush.c
index 96849f4..c59d110 100644
--- a/src/bush.c
+++ b/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]);
 }
 
diff --git a/src/convexcombination.c b/src/convexcombination.c
index c26061a..04104a1 100644
--- a/src/convexcombination.c
+++ b/src/convexcombination.c
@@ -249,7 +249,7 @@ void updateLinks(network_type *network, bool rectifyFlows) {
 
     for (ij = 0; ij < network->numArcs; ij++) {
         network->arcs[ij].cost =
-                network->arcs[ij].calculateCost(&network->arcs[ij]);
+                network->arcs[ij].calculateCost(&network->arcs[ij], NO_FIXED_COST);
         network->arcs[ij].der =
                 network->arcs[ij].calculateDer(&network->arcs[ij]);
     }
@@ -282,12 +282,12 @@ double bisection(network_type *network, double **direction, int iteration,
         lambda = (lmax + lmin) / 2;
         der = 0;
         for (c = 0; c < network->numClasses; c++) {
-            changeFixedCosts(network, c);
+//            changeFixedCosts(network, c);
             for (ij = 0; ij < network->numArcs; ij++) {
                 /* Get total flow based on complete shift */
                 tempFlow = network->arcs[ij].flow 
                             + lambda * direction[ij][network->numClasses];
-                der += evaluateLinkCost(&network->arcs[ij], tempFlow)
+                der += evaluateLinkCost(&network->arcs[ij], tempFlow, c)
                         * direction[ij][c];
             }
         }
@@ -338,9 +338,9 @@ double NewtonIteration(network_type *network, double **direction, double lmin,
     int ij, c;
     double lambda = 0, numer = 0, denom = 0;
     for (c = 0; c < network->numClasses; c++) {
-        changeFixedCosts(network, c);
+//        changeFixedCosts(network, c);
         for (ij = 0; ij < network->numArcs; ij++) {
-            numer += network->arcs[ij].calculateCost(&network->arcs[ij]) 
+            numer += network->arcs[ij].calculateCost(&network->arcs[ij], c)
                     * direction[ij][c];
         }
     }
@@ -404,7 +404,7 @@ void AONdirection(network_type *network, double **direction,
 #endif
 
     for (c = 0; c < network->numClasses; c++) {
-        changeFixedCosts(network, c);
+//        changeFixedCosts(network, c);
 #if PARALLELISM
         for (r = 0; r < network->numZones; ++r) {
             args[r].clss = c;
@@ -465,7 +465,7 @@ void CFWdirection(network_type *network, double **direction,
     /* Now calculate parameters for conjugacy */
     numer = 0; denom = 0;
     for (c = 0; c < network->numClasses; c++) {
-        changeFixedCosts(network, c);
+//        changeFixedCosts(network, c);
         for (ij = 0; ij < network->numArcs; ij++) {
             numer += network->arcs[ij].der
                         * AON[ij][c] * oldDirection[ij][c];
@@ -519,7 +519,7 @@ void BFWdirection(network_type *network, double **direction,
     /* Now calculate parameters for conjugacy */
     numer = 0; denom = 0;
     for (c = 0; c < network->numClasses; c++) {
-        changeFixedCosts(network, c);
+//        changeFixedCosts(network, c);
         for (ij = 0; ij < network->numArcs; ij++) {
             numer += network->arcs[ij].der
                         * AON[ij][c] * oldOldDirection[ij][c];
@@ -534,7 +534,7 @@ void BFWdirection(network_type *network, double **direction,
 
     numer = 0; denom = 0;
     for (c = 0; c < network->numClasses; c++) {
-        changeFixedCosts(network, c);
+//        changeFixedCosts(network, c);
         for (ij = 0; ij < network->numArcs; ij++) {
             numer += network->arcs[ij].der
                         * AON[ij][c] * oldDirection[ij][c];
@@ -807,11 +807,11 @@ double CCaverageExcessCost(network_type *network, double tstt, double sptt) {
 }
 
 /* Evaluate a link performance function WITHOUT changing the underlying flow */
-double evaluateLinkCost(arc_type *arc, double flow) {
+double evaluateLinkCost(arc_type *arc, double flow, int c) {
     double oldFlow = arc->flow;
     double cost;
     arc->flow = flow;
-    cost = arc->calculateCost(arc);
+    cost = arc->calculateCost(arc, c);
     arc->flow = oldFlow;
     return cost;
 }
diff --git a/src/networks.c b/src/networks.c
index 75a753f..226ce55 100644
--- a/src/networks.c
+++ b/src/networks.c
@@ -500,28 +500,6 @@ void heapDijkstra(int origin, double *label, int *backnode, network_type
 }
 
 
-/*
-changeFixedCosts: Updates the network costs for the new specified toll and 
-distance factors.
-
-*/
-void changeFixedCosts(network_type *network, int class) {
-    int ij;
-    double oldCost;
-
-    for (ij = 0; ij < network->numArcs; ij++) {
-#ifdef PARALLELISM
-        pthread_mutex_lock(&network->arc_muts[ij]);
-#endif
-        oldCost = network->arcs[ij].fixedCost;
-        network->arcs[ij].fixedCost = network->arcs[ij].classCost[class];
-        network->arcs[ij].cost += (network->arcs[ij].fixedCost - oldCost);
-#ifdef PARALLELISM
-        pthread_mutex_unlock(&network->arc_muts[ij]);
-#endif
-    }
-}
-
 /*
 finalizeNetwork: After adding the links and nodes to the network struct, this
 function generates the forward and reverse star lists. 
diff --git a/src/parallel_bush.c b/src/parallel_bush.c
index 52f81e8..3a166b4 100644
--- a/src/parallel_bush.c
+++ b/src/parallel_bush.c
@@ -647,16 +647,16 @@ void checkFlows_par(network_type *network, bushes_type *bushes, int t_id) {
     deleteVector(flowCheck);
 }
 
-void exactCostUpdate_par(int ij, double shift, network_type *network, int class) {
-    network->arcs[ij].classFlow[class] += shift;
+void exactCostUpdate_par(int ij, double shift, network_type *network, int c) {
+    network->arcs[ij].classFlow[c] += shift;
     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]);
 }
 
-void classUpdate_par(int hi, int class, double shift,  network_type *network) {
+void classUpdate_par(int hi, int c, double shift,  network_type *network) {
     pthread_mutex_lock(&network->arc_muts[hi]);
-    network->arcs[hi].classFlow[class] -= shift;
+    network->arcs[hi].classFlow[c] -= shift;
     pthread_mutex_lock(&network->arc_muts[hi]);
 }
 
diff --git a/src/tap.c b/src/tap.c
index 0653237..4215ae6 100644
--- a/src/tap.c
+++ b/src/tap.c
@@ -30,7 +30,7 @@ double BeckmannFunction(network_type *network) {
     double Beckmann = 0;
     int ij, c;
     for (ij = 0; ij < network->numArcs; ij++) {
-        Beckmann += network->arcs[ij].calculateInt(&network->arcs[ij], FALSE);
+        Beckmann += network->arcs[ij].calculateInt(&network->arcs[ij], NO_FIXED_COST);
         for (c = 0; c < network->numClasses; c++) {
             Beckmann += network->arcs[ij].classFlow[c]
                        * network->arcs[ij].classCost[c];
@@ -120,13 +120,14 @@ double linkCost(arc_type *arc, cost_type costFunction) {
 
 /*
  * generalBPRcost -- Evaluates the BPR function for an arbitrary polynomial.
+ * If -1 is provided for second argument, fixed cost is excluded from calculation.
  */
-double generalBPRcost(struct arc_type *arc) {
+double generalBPRcost(struct arc_type *arc, int c) {
    if (arc->flow <= 0)
    // Protect against negative flow values and 0^0 errors
-       return arc->freeFlowTime + arc->fixedCost;
+       return arc->freeFlowTime + (c >= 0 ? arc->classCost[c] : 0);
 
-   return arc->fixedCost + arc->freeFlowTime *
+   return (c >= 0 ? arc->classCost[c] : 0) + arc->freeFlowTime *
        (1 + arc->alpha * pow(arc->flow / arc->capacity, arc->beta));
 }
 
@@ -148,19 +149,22 @@ double generalBPRder(struct arc_type *arc) {
 
 /*
  * generalBPRint -- Evaluate integral of an arbitrary polynomial BPR function.
+ * If -1 is provided for second argument, fixed cost is excluded from calculation.
  */
-double generalBPRint(struct arc_type *arc, bool includeFixedCost) {
+double generalBPRint(struct arc_type *arc, int c) {
    if (arc->flow <= 0) return 0; /* Protect against negative flow values and 
                                     0^0 errors */
-   return arc->flow * (includeFixedCost == TRUE ? arc->fixedCost : 0
+   return arc->flow * ((c >= 0 ? arc->classCost[c] : 0)
            + arc->freeFlowTime * 
            (1 + arc->alpha / (arc->beta + 1) * 
             pow(arc->flow / arc->capacity, arc->beta)));
 }
 
-/* linearBPRcost/der/int -- Faster implementation for linear BPR functions. */
-double linearBPRcost(struct arc_type *arc) {
-   return arc->fixedCost + arc->freeFlowTime *
+/* linearBPRcost/der/int -- Faster implementation for linear BPR functions.
+ * If -1 is provided for second argument, fixed cost is excluded from calculation of cost and integral.
+ * */
+double linearBPRcost(struct arc_type *arc, int c) {
+   return (c >= 0 ? arc->classCost[c] : 0) + arc->freeFlowTime *
        (1 + arc->alpha * arc->flow / arc->capacity);
 }
 
@@ -168,18 +172,19 @@ double linearBPRder(struct arc_type *arc) {
    return arc->freeFlowTime * arc->alpha / arc->capacity;
 }
 
-double linearBPRint(struct arc_type *arc, bool includeFixedCost) {
-   return arc->flow * (includeFixedCost == TRUE ? arc->fixedCost : 0
+double linearBPRint(struct arc_type *arc, int c) {
+   return arc->flow * ((c >= 0 ? arc->classCost[c] : 0)
            + arc->freeFlowTime * (1 + arc->flow*arc->alpha/arc->capacity/2));
 }
 
 /* quarticBPRcost/der/int -- Faster implementation for 4th-power BPR functions
+ * If -1 is provided for second argument, fixed cost is excluded from calculation of cost and integral.
  */
-double quarticBPRcost(struct arc_type *arc) {
+double quarticBPRcost(struct arc_type *arc, int c) {
    double y = arc->flow / arc->capacity;
    y *= y;
    y *= y;
-   return arc->fixedCost + arc->freeFlowTime * (1 + arc->alpha * y);
+   return (c >= 0 ? arc->classCost[c] : 0) + arc->freeFlowTime * (1 + arc->alpha * y);
 }
 
 double quarticBPRder(struct arc_type *arc) {
@@ -190,16 +195,19 @@ double quarticBPRder(struct arc_type *arc) {
 
 }
 
-double quarticBPRint(struct arc_type *arc, bool includeFixedCost) {
+double quarticBPRint(struct arc_type *arc, int c) {
    double y = arc->flow / arc->capacity;
    y *= y;
    y *= y;
-   return arc->flow * (includeFixedCost == TRUE ? arc->fixedCost : 0
+   return arc->flow * ((c >= 0 ? arc->classCost[c] : 0)
            + arc->freeFlowTime * (1 + arc->alpha * y / 5));
 }
 
-double conicCost(struct arc_type *arc) {
-    double time = arc->freeFlowTime + arc->fixedCost;
+/* conicCost/der/int -- Conic Delay calculations
+ * If -1 is provided for second argument, fixed cost is excluded from calculation of cost and integral.
+ */
+double conicCost(struct arc_type *arc, int c) {
+    double time = arc->freeFlowTime + (c >= 0 ? arc->classCost[c] : 0);
     double x = arc->flow/arc->capacity;
     double v_s = arc->flow/arc->saturationFlow;
 
@@ -253,11 +261,11 @@ double conicDer(struct arc_type *arc) {
     return der;
 }
 
-double conicInt(struct arc_type *arc, bool includeFixedCost) {
+double conicInt(struct arc_type *arc, int c) {
     fatalError("Conic integrals not yet implemented; set calculateBeckmann "
                "to FALSE in parameters.");
     /* Suppress compiler warning about unused arguments */
-    displayMessage(FULL_DEBUG, "%p%d", arc, (int) includeFixedCost);
+    displayMessage(FULL_DEBUG, "%p%d", arc, (int) c);
     return IS_MISSING;
 }
 
@@ -271,13 +279,11 @@ network struct -- this allows SPTT to be calculated without changing any values
 in the network, at the expense of a little more run time.
 */
 double SPTT(network_type *network) {
-    int r, j, c, originNode;
+    int r, j, originNode;
     double sptt = 0;
     declareVector(double, SPcosts, network->numNodes);
     for (r = 0; r < network->numOrigins; r++) {
         originNode = origin2node(network, r);
-        c = origin2class(network, r);
-        changeFixedCosts(network, c);
         BellmanFord_NoLabel(originNode, SPcosts, network, DEQUE, NULL, NULL);
         for (j = 0; j < network->numZones; j++) {
             sptt += network->demand[r][j] * SPcosts[j];
@@ -316,7 +322,7 @@ double classCost(network_type *network, int class, double timeFactor,
     for (ij = 0; ij < network->numArcs; ij++) {
         sum += network->arcs[ij].classFlow[class]
                * (timeFactor *
-                    (network->arcs[ij].cost - network->arcs[ij].fixedCost)
+                    (network->arcs[ij].cost - network->arcs[ij].classCost[class])
                   + distanceFactor * network->arcs[ij].length
                   + tollFactor * network->arcs[ij].classToll[class]);
     }
@@ -341,7 +347,7 @@ double classTravelTime(network_type *network, int class) {
 double classGeneralizedCost(network_type *network, int class) {
     int ij;
     double sum = 0;
-    changeFixedCosts(network, class);
+//    changeFixedCosts(network, class);
     for (ij = 0; ij < network->numArcs; ij++) {
         sum += network->arcs[ij].classFlow[class] * network->arcs[ij].cost;
     }
@@ -356,7 +362,7 @@ according to costFunction
 void updateAllCosts(network_type *network) {
     int i;
     for (i = 0; i < network->numArcs; i++)
-      network->arcs[i].cost=network->arcs[i].calculateCost(&network->arcs[i]);
+      network->arcs[i].cost=network->arcs[i].calculateCost(&network->arcs[i], NO_FIXED_COST);
 }
 
 /*