diff --git a/EnKF.py b/EnKF.py
new file mode 100644
index 0000000..3a3dd35
--- /dev/null
+++ b/EnKF.py
@@ -0,0 +1,250 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+this class is for EnKF operations
+
+@author: cesny
+"""
+import numpy as np
+from utils import setCTMVehicles, forwardCTMPropagation, CTMcreateInitialEnsemble, VmaxCreateInitialEnsemble, VmaxtoCritDen, cellToLength, lengthToCell
+
+
+class EnKF:
+  '''
+  this class is used to implement EnKF operations
+  for implementation details, this code follows Evensen2003 
+  '''
+  def __init__(self, obsError, modelError, sampleSize, stateDim, obsDim, m=None, assimilatedDensities=None, H=None, EnKFtype='CTM', nonLinearObs=False, droneLoc = None, trafficNet=None, droneDenObsError=None):
+    self.obsError = obsError  # specifies standard dev. of observ. white noise
+    self.modelError = modelError  # specifies standard dev. of model white noise
+    self.sampleSize = sampleSize  # number of ensemble members
+    self.stateDim = stateDim  # dimension of an ensemble member
+    self.obsDim = obsDim  # dimension of observation vector
+    self.EnKFtype = EnKFtype  # could be 'CTM' or 'Vmax'
+    self.nonLinearObs = nonLinearObs
+    self.assimDen = assimilatedDensities  # a list with two elements, assim. den upstream and assim den downstream
+    self.m = m  # this is a function that takes as input the state vector as a list [vmax1, vmax2] and returns the model prediction of the parameters [v1, v2]
+    self.droneLoc = droneLoc  # stores the drone location (linkID, cell) tuple
+    self.H = H  # create the matrix (vector) H if it is available
+    self.locToCell = dict()  # maps the tuple (link, cell) to cell between 0 and 40
+    self.cellToLoc = dict()
+    self.trafficNet = trafficNet
+    self.droneDenObsError = droneDenObsError
+    # store data!
+    self.storePropEnsembles = list()
+    self.storeAhat = list()
+    self.storeA = list()
+    self.storeKalman = list()
+    self.storeD = list()
+    self.storeDmA = list()
+    self.storeInvPart = list()
+    self.storeCovPart = list()
+    self.storeAhatPrime = list()
+
+  def createLocToCell(self):
+    '''
+    creates the mapping between (linkid, cell) to
+    global cellindex between (0,40) across all links
+    '''
+    cellindex = 0
+    for linkID in self.trafficNet.linkDict:
+      if linkID is not 9:
+        for cellkey, cell in enumerate(self.trafficNet.linkDict[linkID].cells):
+          self.locToCell[(linkID, cellkey)] = cellindex
+          self.cellToLoc[cellindex] = (linkID, cellkey)
+          cellindex += 1
+    return None
+  
+  def genModErrorMatrix(self):
+    '''
+    generates a matrix of perturbations
+    that will be used to perturb CTM forecasts
+    '''
+    self.modelErrorMatrix = np.random.normal(loc=0.0, scale=self.modelError, size=(self.stateDim, self.sampleSize))  # for general multivariate normal, this could have been generated using  numpy.random.multivariate_normal(mean, cov)
+    return None
+  
+  def genObsErrorMatrix(self):
+    '''
+    generate a matrix of perturbations
+    that will be used to perturb sensor
+    observations
+    Adjusts for randomness based on 
+    congestion state
+    note for future only considering drone obs in embedded EnKFs
+    '''
+    if (self.EnKFtype is 'CTM') and (self.droneLoc is not None):
+      droneCell = self.locToCell[self.droneLoc]
+      self.obsErrorMatrix = np.random.normal(loc=0.0, scale=self.obsError, size=(self.obsDim, self.sampleSize))  # for general multivariate normal, this could have been generated using  numpy.random.multivariate_normal(mean, cov)
+      self.obsErrorMatrix[droneCell] = np.random.normal(loc=0.0, scale=self.droneDenObsError, size=self.sampleSize)  # lower error at location of  drone!
+    else:
+      self.obsErrorMatrix = np.random.normal(loc=0.0, scale=self.obsError, size=(self.obsDim, self.sampleSize))
+    return None
+      
+  def getUpdatedEnsembles(self):
+    '''
+    get the updated ensembles in list of lists format
+    needed for CTM processing, each sublist is a list
+    of densities
+    '''
+    Atranspose = np.transpose(self.A)
+    listofLists = Atranspose.tolist()
+    return listofLists
+  
+  def getMean(self):
+    '''
+    adds some layer of protection
+    '''
+    return self.mean
+  
+  def getP(self):
+    '''
+    return P, scales by 1/(N-1) as needed
+    in rest of code here this scale is 
+    ignored since it cancels out
+    keeps self.P intact
+    '''
+    P = (1.0/(self.sampleSize-1)) * self.P
+    return P
+    
+  def getPostDist(self):
+    '''
+    updates mean and covariance matrix based on
+    observations
+    '''
+    if self.nonLinearObs is False:
+      self.A = self.A + np.dot(self.K, self.D - np.dot(self.H, self.A))
+      scaleMatrix = np.full((self.sampleSize, self.sampleSize), 1.0/self.sampleSize)
+      self.Abar = np.dot(self.A, scaleMatrix)
+      self.mean = self.Abar[:,0]
+      self.P = self.P - np.dot(self.K, np.dot(self.H, self.P))
+      
+    elif self.nonLinearObs is True:
+      self.A = self.A + np.dot(self.K, self.D - self.Ahat)
+      self.storeDmA.append(self.D - self.Ahat)
+      scaleMatrix = np.full((self.sampleSize, self.sampleSize), 1.0/self.sampleSize)
+      self.Abar = np.dot(self.A, scaleMatrix)
+      self.mean = self.Abar[:,0]
+      self.Aprime = self.A - self.Abar
+      self.P = np.dot(self.Aprime, np.transpose(self.Aprime))
+      self.storeA.append(self.A)
+    return self.mean, self.P
+    
+  def getKalmanGain(self):
+    '''
+    computes the Kalman gain
+    '''
+    if self.nonLinearObs is False:
+      temp1 = np.dot(self.P, np.transpose(self.H))
+      temp2 = np.dot(self.H, self.P)
+      temp2 = np.dot(temp2, np.transpose(self.H))
+      temp2 = temp2 + self.R
+      temp2 = np.linalg.inv(temp2)
+      self.K = np.dot(temp1, temp2)
+    elif self.nonLinearObs is True:
+      Ahat0 = self.m(self.A[0],self.assimDen[0])  # Warning, works specifically with problem at hand
+      Ahat1 = self.m(self.A[1], self.assimDen[1])
+      self.Ahat = np.stack((Ahat0,Ahat1)) # compute Ahat through the m function
+      self.storeAhat.append(self.Ahat)
+      scaleMatrix = np.full((self.sampleSize, self.sampleSize), 1.0/self.sampleSize)
+      self.Ahatbar = np.dot(self.Ahat, scaleMatrix)
+      self.AhatPrime = self.Ahat - self.Ahatbar
+      self.storeAhatPrime.append(self.AhatPrime)
+      temp1 = np.dot(self.Aprime, np.transpose(self.AhatPrime))
+      self.storeCovPart.append(temp1)
+      temp2 = np.dot(self.AhatPrime, np.transpose(self.AhatPrime))
+      temp2 = temp2 + self.R
+      temp2 = np.linalg.inv(temp2)
+      self.storeInvPart.append(temp2)
+      self.K = np.dot(temp1, temp2)
+      self.storeKalman.append(self.K)
+    return None
+  
+  def getPriorDist(self):
+    '''
+    generate prior t+1|t
+    updates mean and covariance based on
+    predictions
+    '''
+    scaleMatrix = np.full((self.sampleSize, self.sampleSize), 1.0/self.sampleSize)
+    self.Abar = np.dot(self.A, scaleMatrix)
+    self.mean = self.Abar[:,0]
+    self.Aprime = self.A - self.Abar
+    self.P = np.dot(self.Aprime, np.transpose(self.Aprime))
+    # self.P = (1.0/(self.sampleSize-1)) * self.P
+    return self.mean, self.P
+  
+  def addModelNoise(self, forecasts):
+    '''
+    forecasts is a list of lists, where sublists
+    are model updates for different ensemble
+    members
+    --------
+    NOTE:
+    forecasts do not include noise, this method
+    adds the noise!
+    --------
+    returns noisy forecasts
+    '''
+    self.A = np.array(forecasts)
+    self.A = np.transpose(self.A)
+    self.genModErrorMatrix()
+    self.A = self.A + self.modelErrorMatrix
+    self.storePropEnsembles.append(self.A)
+    return None
+
+  def addObsNoise(self, observations):
+    '''
+    observations is a single list of
+    sensor measurements
+    adds obs noise to observations
+    returns noisy observations
+    if model has different congestion states
+    adjusts observation noise accordingly
+    in genObsErrorMatrix
+    '''
+    tempList = list()
+    for sample in range(self.sampleSize):
+      tempList.append(observations)
+    self.D = np.array(tempList)
+    self.D = np.transpose(self.D)
+    self.genObsErrorMatrix()
+    self.D = self.D + self.obsErrorMatrix
+    self.storeD.append(self.D)
+    return None
+  
+  def getObsCov(self):
+    '''
+    get the observation covariance matrix
+    '''
+    self.R = np.dot(self.obsErrorMatrix, np.transpose(self.obsErrorMatrix))
+    # self.R = (1.0/(self.sampleSize-1)) * self.R
+    return None
+  
+  def EnKFStep(self, forecasts, observations):
+    '''
+    implements one EnKF forecast and
+    assimilation step and updates
+    mean and covariance (get posterior)
+    then samples and returns next step
+    assimilated states
+    '''
+    self.createLocToCell()
+    self.addModelNoise(forecasts)
+    self.addObsNoise(observations)
+    self.getObsCov()
+    self.getPriorDist()
+    self.getKalmanGain()
+    self.getPostDist()
+    return self.getUpdatedEnsembles()
+
+
+
+
+
+
+
+
+
+    
+    
+  
diff --git a/findPath.py b/findPath.py
new file mode 100644
index 0000000..b769934
--- /dev/null
+++ b/findPath.py
@@ -0,0 +1,187 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+this module is for implementing one step lookahead control 
+
+determines drone location at next step based on information maximizing
+control algorithm
+
+@author: cesny
+"""
+import numpy as np
+from utils import setCTMVehicles, forwardCTMPropagation, CTMcreateInitialEnsemble, VmaxCreateInitialEnsemble, m, VmaxtoCritDen, cellToLength, lengthToCell
+
+
+
+
+class findPath:
+  '''
+  this class is for determining next drone 
+  location based on A-optimal control
+  '''
+  def __init__(self, location, time, trafficNet, EnKFCTM, EnKFV, timeHorizon=None, weight=0.5):
+    self.location = location  # current drone location, defined as a tuple (linkID, cell), cell count from zero
+    self.time = time  # current time
+    self.timeHorizon = timeHorizon  # time horizon to do MPC (number of timeSteps), set dynamically till drone visits all cells in each path, can assign otherwise
+    self.trafficNet = trafficNet  # the traffic network class
+    self.EnKFCTM = EnKFCTM  # the ensemble kalman filtering class for densities
+    self.EnKFV = EnKFV  # the ensemble kalman filtering class, use to predict covariance matrix
+    self.cellToLoc = dict()
+    self.locToCell = dict()  # maps the tuple (link, cell) to cell between 0 and 40
+    self.weight = weight  # this is the weight of vmax vs densities trace, weight corresponds to vmax, (1-w) corresponds to weight of densities trace
+  
+  def createLocToCell(self):
+    '''
+    creates the mapping between (linkid, cell) to
+    global cellindex between (0,40) across all links
+    '''
+    cellindex = 0
+    for linkID in self.trafficNet.linkDict:
+      if linkID is not 9:
+        for cellkey, cell in enumerate(self.trafficNet.linkDict[linkID].cells):
+          self.locToCell[(linkID, cellkey)] = cellindex
+          self.cellToLoc[cellindex] = (linkID, cellkey)
+          cellindex += 1
+    return None
+  
+  def generateDronePaths(self):
+    '''
+    gets the possible paths the drone can take based
+    on its current location. A path is a collection of
+    cells. Returns self.paths a dict of dicts determining the
+    cellid at every time step for possible paths
+    '''
+    self.dronePaths = dict()
+    self.dronePaths['left'] = dict()
+    self.dronePaths['right'] = dict()
+    droneCell = self.locToCell[self.location]
+    leftPath = list(range(0,droneCell+1))
+    leftPath.reverse()
+    rightPath = list(range(droneCell, 40))
+    # assuming at every time step the drone can move one cell, determine its path across time
+    for key, lcell in enumerate(leftPath):
+      if self.time+key <= 449:  # limited by time horizon
+        self.dronePaths['left'][self.time+key] = lcell
+    for key, rcell in enumerate(rightPath):
+      if self.time+key <= 449:
+        self.dronePaths['right'][self.time+key] = rcell
+    return None
+  
+  def getObservations(self):
+    '''
+    use traffic network to simulate paths based on current
+    condition and return expected observations for the given
+    paths (left path and right path). Observations is a list over
+    the 40 cells.
+    '''
+    self.pathObservations = dict()  # a dictionary with mean traffic state at every time step from now as predicted by propagated CTM ensembles
+    self.pathObservations['left'] = dict()
+    self.pathObservations['right'] = dict()
+    loadRange = max(len(self.dronePaths['left']), len(self.dronePaths['right']))
+    storeResults = dict()
+    CTMensembles = self.EnKFCTM.getUpdatedEnsembles()  # current ensembles
+    for lr in range(loadRange):
+      CTMensembles = forwardCTMPropagation(self.time + lr, self.trafficNet, CTMensembles)
+      storeResults[self.time + lr] =  [float(sum(col))/len(col) for col in zip(*CTMensembles)]  # store average of propagated ensembles as expected observed true state
+    
+    for time in self.dronePaths['left']:
+      self.pathObservations['left'][time] = storeResults[time]
+    for time in self.dronePaths['right']:
+      self.pathObservations['right'][time] = storeResults[time]
+    return None
+  
+  def getCovarianceMatrices(self):
+    '''
+    use the "observations" (from propagated ensembles) to determine
+    the final covariance matrix on the densities and vmax
+    covariance on vmax set manually based on actual conditions, covariance
+    of densities computed from EnKF with CTM forward simulation results as 
+    expected observations, drone observations accounted for by lower observation
+    error
+    '''
+    self.finalCovariancesCTM = dict()
+    self.finalCovariancesVmax = dict()
+    # densities covariance matrix
+    CTMensemblesLeft = self.EnKFCTM.getUpdatedEnsembles()
+    CTMensemblesRight = self.EnKFCTM.getUpdatedEnsembles()
+    VmaxEnsemblesLeft = self.EnKFV.getUpdatedEnsembles()
+    VmaxEnsemblesRight = self.EnKFV.getUpdatedEnsembles()
+    for time in self.dronePaths['left']:
+      CTMensemblesLeft = forwardCTMPropagation(time, self.trafficNet, CTMensemblesLeft)
+      self.EnKFCTM.droneLoc = self.cellToLoc[self.dronePaths['left'][time]]  # update drone location according to base policy, used for precise observations
+      CTMensemblesLeft = self.EnKFCTM.EnKFStep(CTMensemblesLeft, self.pathObservations['left'][time])  # update the ensembles
+    self.finalCovariancesCTM['left'] = self.EnKFCTM.getP()  
+    
+    for time in self.dronePaths['right']:
+      CTMensemblesRight = forwardCTMPropagation(time, self.trafficNet, CTMensemblesRight)
+      self.EnKFCTM.droneLoc = self.cellToLoc[self.dronePaths['right'][time]]  # update drone location according to base policy, used for precise observations
+      CTMensemblesRight = self.EnKFCTM.EnKFStep(CTMensemblesRight, self.pathObservations['right'][time])
+    self.finalCovariancesCTM['right'] = self.EnKFCTM.getP()
+    
+    # update the uf estimates
+    self.EnKFV.obsError = 10   
+    self.EnKFV.obsDim = 1
+    EnKFVmean = self.EnKFV.getMean()
+    # VfObserved = [20.0]
+    self.EnKFV.nonLinearObs = False
+    # do left
+    self.EnKFV.H = np.array([[1.0, 0.0]])
+    print(' ')
+    print('pre-update left: ', np.transpose(np.array(VmaxEnsemblesLeft))[:,0:3])
+    VmaxEnsemblesLeft = self.EnKFV.EnKFStep(VmaxEnsemblesLeft, [EnKFVmean[0]])
+    self.finalCovariancesVmax['left'] = self.EnKFV.getP()
+    print('left path: ', self.finalCovariancesVmax['left'], np.transpose(np.array(VmaxEnsemblesLeft))[:,0:3])
+    print(' ')
+    # do right
+    self.EnKFV.H = np.array([[0.0, 1.0]])
+    print('pre-update right: ', np.transpose(np.array(VmaxEnsemblesRight))[:,0:3])
+    VmaxEnsemblesRight = self.EnKFV.EnKFStep(VmaxEnsemblesRight, [EnKFVmean[1]])
+    self.finalCovariancesVmax['right'] = self.EnKFV.getP()
+    print('right path: ', self.finalCovariancesVmax['right'], np.transpose(np.array(VmaxEnsemblesRight))[:,0:3])
+    print(' ')
+    return None
+    
+  def getObjective(self):
+    '''
+    determines the objective function as a weighted form
+    of the trace of the covariance matrices for densities and vmax
+    corrects for dimension and scale differences
+    returns objective of left path and objective of right path in a dict
+    '''
+    self.ObjectiveVal = dict()
+    self.ObjectiveVal['left'] = (self.weight * np.trace(self.finalCovariancesVmax['left']) / 2.0) + ((1-self.weight) * np.trace(self.finalCovariancesCTM['left']) / 40.0)
+    self.ObjectiveVal['right'] = (self.weight * np.trace(self.finalCovariancesVmax['right']) / 2.0) + ((1-self.weight) * np.trace(self.finalCovariancesCTM['right']) / 40.0)
+    print(' ')
+    print('objective left: ', np.trace(self.finalCovariancesVmax['left']) / 2.0, np.trace(self.finalCovariancesCTM['left']) / 40.0)
+    print('objective right: ', np.trace(self.finalCovariancesVmax['right']) / 2.0, np.trace(self.finalCovariancesCTM['right']) / 40.0)
+    return None
+  
+  def moveDrone(self, direction):
+    '''
+    changes the location of the drone
+    '''
+    if direction == 'left':
+      currentCell = self.locToCell[self.location]
+      if currentCell != 0:
+        newLoc = currentCell - 1
+        self.location =  self.cellToLoc[newLoc]
+        
+    elif direction == 'right':
+      currentCell = self.locToCell[self.location]
+      if currentCell != 39:
+        newLoc = currentCell + 1
+        self.location =  self.cellToLoc[newLoc]
+    return self.location
+  
+  def updateLocation(self):
+    '''
+    determines the next step for the drone
+    go left or go right, updates self.location
+    '''
+    self.createLocToCell()
+    self.generateDronePaths()
+    self.getObservations()
+    self.getCovarianceMatrices()
+    self.getObjective()
+    minKey = min(self.ObjectiveVal, key=self.ObjectiveVal.get)
+    return self.moveDrone(minKey)
\ No newline at end of file
diff --git a/link.py b/link.py
new file mode 100644
index 0000000..1b157b0
--- /dev/null
+++ b/link.py
@@ -0,0 +1,162 @@
+# -*- coding: utf-8 -*-
+"""
+abstract link class
+includes methods for updating parameters using parameter EnKF estimates
+
+@author: cesny
+"""
+
+
+# link parameters, here common for all links, can be different, read in readLinks
+critDen = 100.0 # (veh/km)
+linkType = 'CTM'  # specifies the link type
+ffs = 100.0  # free flow speed in km/hr
+jamDen = 300.0  # veh/km
+qcap = ffs * critDen  # veh/hr  (assuming triangular FD)
+bws = (ffs * critDen) / (jamDen - critDen)  # backward wave speed  (km/hr)
+length = 1.5  # km
+timeStep = 10  # seconds  (needed to define num cells in case of CTM, delx=ufdelt)
+params = {'critDen':critDen, 'ffs':ffs, 'jamDen': jamDen, 'qcap': qcap, 
+          'bws':bws, 'length':length, 'timeStep': timeStep, 'linkType': 'CTM'}
+
+
+class Link:
+  """
+  this class is for link objects
+  upstreamPathCount is a dict of dicts, time: subdict
+  subdict is a dict of paths and their flow values
+  """
+  def __init__(self, linkID, unode, dnode, params):
+    self.ID = linkID
+    self.unode = unode  # upstream node object
+    self.dnode = dnode
+    self.params = params
+    self.params['qcap'] = params['ffs'] * params['critDen']  # cap in veh/hr
+    self.params['bws'] = (self.params['ffs'] * self.params['critDen']) / (self.params['jamDen'] - self.params['critDen'])  # backward wave speed in km/hr
+    self._upstreamCounts = dict()  # records cumulative counts across time steps
+    self._downstreamCounts = dict()
+    self.inFlow = 0
+    self.outFlow = 0  # current inFlow and outFlow into the link
+    
+  def calculateSendingFlow(self, time, timeStep):
+    """
+    overwrite in subclass
+    """
+    pass
+  
+  def calculateReceivingFlow(self, time, timeStep):
+    """
+    overwrite in subclass
+    """
+    pass
+  
+  def linkUpdate(self, time):
+    """
+    performs any internal calculations, puts flows on links, and 
+    removes flows from downstream end of link
+    """
+    self.flowIn(time)
+    self.flowOut(time)
+    return None
+  
+  def upstreamCount(self, time):
+    """
+    return the cumulative entries to a link up to time t
+    """
+    if time < 0:
+      return 0
+    return self._upstreamCounts(time)
+  
+  def downstreamCount(self, time):
+    """
+    returns the cumulative exists from a link
+    """
+    if time < 0:
+      return 0
+    return self._downstreamCounts(time)
+  
+  def vehiclesOnLink(self, time):
+    """
+    returns vehicles on link
+    """
+    return self.upstreamCount(time) - self.downstreamCount(time)
+  
+  def flowIn(self, time):
+    """
+    adds flow to link based on inFlows from 
+    transitionFlows
+    """
+    self._upstreamCounts[time] = self.inFlow  # I believe that this should be +=self.inflow but for our purposes it does not matter
+    return None
+  
+  def flowOut(self, time):
+    """
+    removes flow from the downstream end of the link
+    based on outFlow from transitionFlows
+    """
+    self._downstreamCounts[time] = self.outFlow
+    return None
+  
+  def linkDensity(self, time=None):
+    """
+    overwriten in cell transmission model to return
+    density at every cell
+    """
+    return self.vehiclesOnLink(time) / self.params['length']  # density in veh/km
+  
+  def updateVmaxCritDen(self, newVmax, newCritDen):
+    """
+    update the maximum speed
+    note assuming that the backward wave speed is fixed
+    """
+    if newVmax > 110:
+      print('.. WARNING! CFL condition violated ...')
+    self.params['ffs'] = newVmax
+    self.params['critDen'] = newCritDen
+    self.params['qcap'] = newVmax * newCritDen
+    # self.params['bws'] = (newVmax * self.params['critDen']) / (self.params['jamDen'] - self.params['critDen'])  # remains fixed across iterations
+    return None
+  
+      
+      
+      
+      
+      
+      
+      
+      
+      
+      
+      
+      
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+    
+    
+    
+    
+    
+    
+
+
+
+
+
diff --git a/linkModel.py b/linkModel.py
new file mode 100644
index 0000000..72bfb8e
--- /dev/null
+++ b/linkModel.py
@@ -0,0 +1,113 @@
+# -*- coding: utf-8 -*-
+"""
+cell transmission model
+
+@author: cesny
+"""
+import math
+from link import Link
+
+class Cell:
+  def __init__(self, capacity, maxVehicles, delta, timeStep, ffs):
+    self.capacity = capacity / 3600.0  # you have to divide by 3600 for veh/sec
+    self.maxVehicles = maxVehicles  # maxVehicles per cell
+    self.length = ffs * timeStep / 3600.0  # length in km
+    self.delta = delta
+    self.vehicles = 0
+    self.timeStep = timeStep
+    
+  def calculateSendingFlow(self):
+    return min(self.vehicles, self.capacity * self.timeStep)
+  
+  def calculateReceivingFlow(self):
+    return min(self.delta * (self.maxVehicles - self.vehicles), self.capacity * self.timeStep)
+  
+  def removeVehicles(self, numVehicles):
+    self.vehicles -= numVehicles
+    return None
+  
+  def addVehicles(self, numVehicles):
+    self.vehicles += numVehicles
+    return None
+  
+  def cellDensity(self):
+    return float(self.vehicles) /  self.length  # density in veh/km
+    
+
+class CTM(Link):
+  def __init__(self, linkID, unode, dnode, params):
+    Link.__init__(self, linkID, unode, dnode, params)
+    # create a list of cells
+    self.cells = list()
+    #  freeflowtime is equiv to length/delx*delt
+    #self.freeFlowTime = int((self.params['length']/(self.params['ffs'] / 3600.0) + self.params['timeStep'] - 1) / self.params['timeStep'])
+    # timeStep is used as a parameter since it affects cell size
+    self.numCells = math.ceil(self.params['length']/(self.params['ffs']*(1.0/3600)*self.params['timeStep']))
+    for c in range(self.numCells):
+      newCell = Cell(self.params['qcap'], self.params['jamDen'] * self.params['length'] / self.numCells,
+                     self.params['bws'] / self.params['ffs'], self.params['timeStep'], self.params['ffs'])
+      self.cells.append(newCell)
+    
+  def calculateSendingFlow(self, time, timeStep):
+    """
+    overwrites Link method
+    """
+    return self.cells[-1].calculateSendingFlow()
+  
+  def calculateReceivingFlow(self, time, timeStep):
+    """
+    overwrites Link method
+    """
+    return self.cells[0].calculateReceivingFlow()
+  
+  def linkUpdate(self, time):
+    """
+    performs any internal calculations, puts flows on links, and 
+    removes flows from downstream end of link
+    inFlow and outFlow were calculated based on the t-1
+    sending and receiving flows
+    """    
+    cellTransitionFlow = list()
+    # deal with intermediate cells
+    for c in range(0, self.numCells-1):
+      cellSendingFlow = self.cells[c].calculateSendingFlow()
+      cellReceivingFlow = self.cells[c+1].calculateReceivingFlow()
+      cellTransitionFlow.append(min(cellSendingFlow, cellReceivingFlow))
+    for c in range(0, self.numCells-1):
+      self.cells[c].removeVehicles(cellTransitionFlow[c])
+      self.cells[c+1].addVehicles(cellTransitionFlow[c])
+    self.flowIn(time)
+    self.flowOut(time)
+    return None
+      
+  def flowIn(self, time):
+    Link.flowIn(self, time)
+    self.cells[0].addVehicles(self.inFlow)
+    return None
+    
+  def flowOut(self, time):
+    Link.flowOut(self, time)
+    self.cells[len(self.cells)-1].removeVehicles(self.outFlow)
+    return None
+  
+  def linkDensity(self, time=None):
+    """
+    overwrites link method to return density at every cell
+    """
+    listofDensities = list()
+    for cell in self.cells:
+      listofDensities.append(cell.cellDensity())
+    return listofDensities
+  
+  def updateVmaxCritDen(self, newVmax, newCritDen):
+    """
+    overwrites link method to update cells
+    """
+    Link.updateVmaxCritDen(self, newVmax, newCritDen)
+    # now update cells, keep length and number of cells the same
+    for cell in self.cells:
+      cell.capacity = self.params['qcap'] / 3600.0 # already updated for link
+      cell.delta = self.params['bws'] / self.params['ffs']  # already upated in link 
+    return None
+
+      
\ No newline at end of file
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..f598d48
--- /dev/null
+++ b/main.py
@@ -0,0 +1,178 @@
+
+"""
+primary UAV navigation and estimation framework
+
+@author: cesny
+"""
+
+from network import Network
+import numpy as np
+from EnKF import EnKF
+import copy as cp
+import matplotlib.pyplot as plt
+from findPath import findPath
+from utils import readData, setCTMVehicles, forwardCTMPropagation, CTMcreateInitialEnsemble, VmaxCreateInitialEnsemble, m, VmaxtoCritDen, createLocToCell, cellToLength, lengthToCell
+
+ 
+
+if __name__ == '__main__':
+  # create traffic network
+  simTime = 4490
+  simTimeStep = 10
+  totalTimeSteps = range(int(np.ceil(float(simTime)/simTimeStep)) + 1)
+  linkfile = 'VISSIMnetwork/links.txt'
+  nodefile = 'VISSIMnetwork/nodes.txt'
+  demandfile = 'VISSIMnetwork/demand6600.txt'
+  trafficNet = Network(simTime, simTimeStep, nodefile, linkfile, demandfile)
+  LocToCell = createLocToCell(trafficNet)
+  # laod data observations from VISSIM
+  denData, spData = readData('data/model_001_Link Segment Results-6600.att')
+  
+  # true incident Vf
+  trueVf=20.0
+  
+  # UAV path planning weight lambda
+  pathWeight=1.0
+  
+  # set initial UAV location, link 5 cell 0
+  droneLocation = (5,0)
+  
+  # specify parameters for CTM-EnKF
+  totalcells = 0
+  for linkID in trafficNet.linkDict:
+    if linkID is not 9:
+      totalcells += trafficNet.linkDict[linkID].numCells
+  HCTM = np.array(np.identity(totalcells).tolist())
+  CTMobsSTDV = 10  # standard deviation veh/km
+  CTMdrObsSTDV = 2 # standard deviaiton veh/km of drone observations
+  CTMmodSTDV = 5  # standard deviation veh/km
+  CTMstateDim = 40  # 40 cells with monitored densities
+  CTMobsDim = 40  # 40 cells with monitored densities
+  CTMensembles = 100  # number of ensembles in EnKF
+  EnKFCTM = EnKF(obsError=CTMobsSTDV, modelError=CTMmodSTDV, sampleSize=CTMensembles, stateDim=CTMstateDim, obsDim=CTMobsDim, H=HCTM, droneLoc=droneLocation, trafficNet=trafficNet, droneDenObsError=CTMdrObsSTDV)
+  
+  # specify parameters for velocity EnKF when velocities are observed
+  VobsSTDV = 5 # obs standard deviation!
+  VmodSTDV = 5  # model random walk standard dev. in km/hr
+  VstateDim = 2  # two incident prone regions
+  VobsDimV = 2  # observing velocities on those two regions
+  Vensembles = 100  # number of ensembles in EnKF
+  EnKFV = EnKF(obsError=VobsSTDV, modelError=VmodSTDV, sampleSize=Vensembles, stateDim=VstateDim, obsDim=VobsDimV, m=m, assimilatedDensities=[0,0], EnKFtype='Vmax', nonLinearObs=True, droneLoc=droneLocation, trafficNet=trafficNet)
+  
+  # specify parameters for velocity EnKF when we obtain direct uf observations
+  VdrObsSTDV = 10  # error of observing the true uf value
+  VobsDimVf = 1  # when the drone observes it only does so at one location
+  
+  # creates initial ensembles
+  CTMensembles = CTMcreateInitialEnsemble(CTMstateDim, CTMensembles, CTMmodSTDV)
+  VmaxEnsembles = VmaxCreateInitialEnsemble(VstateDim,Vensembles,VmodSTDV)
+  
+  # store data
+  firstIncidentDen = list()
+  secIncidentDen = list()
+  Vmax1Estimated = list()
+  Vmax2Estimated = list()
+  critDenEstimated = list()
+  listofTime = list()
+  VmaxlistofTime = list()
+  storeDenTotal=list()
+  objective = list()
+  velObj = list()
+  storeDenInc = list()
+  storeDroneLocation = list()
+  droneLocCell=list()
+  droneLocKm=list()
+  
+  # simulate
+  for time in totalTimeSteps:  # cell indices 6 and 32 for inc1 and inc2, respectively (i.e., those are the incident prone locations)
+    CTMensembles = forwardCTMPropagation(time, trafficNet, CTMensembles)  # propagate ensembles using CTM
+    CTMensembles = EnKFCTM.EnKFStep(CTMensembles, denData[time])  # data assimilation, get updated density ensembles from EnKF
+    firstIncidentDen.append(EnKFCTM.mean[6])  # add best estimate of den in incident location to list
+    secIncidentDen.append(EnKFCTM.mean[32])
+    storeDenTotal.append(EnKFCTM.mean)
+    listofTime.append(time*10/60.0)
+    if (time % 30 == 0):  # if you observe a velocity measurement
+      print('velocity observation')
+      VmaxlistofTime.append(time)
+      EnKFV.nonLinearObs = True  # nonlinear relationship when velocity is observed
+      EnKFV.H = None  # there is no H matrix
+      # get assimilated densities
+      assimDensities = [EnKFCTM.mean[6], EnKFCTM.mean[32]]  # current density, use in nonlinear observations to determine vmax
+      EnKFV.assimDen = assimDensities
+      storeDenInc.append(assimDensities)  # stores the denisities at incident-prone locations when velocity measurements are collected
+      # adjust observations
+      EnKFV.obsError = VobsSTDV
+      EnKFV.obsDim = VobsDimV
+      VmaxEnsembles = EnKFV.EnKFStep(VmaxEnsembles, spData[time])  # propagate Vmax ensembles using random walk and assimilate observed data
+      # store results
+      Vmax1Estimated.append(EnKFV.mean[0])
+      Vmax2Estimated.append(EnKFV.mean[1])
+      critDenEstimated.append(VmaxtoCritDen(EnKFV.mean))
+      # update traffic parameters
+      trafficNet.updateVmaxCritDen(EnKFV.mean, VmaxtoCritDen(EnKFV.mean))  # Now the parameters are updates for the links of interest at the locations of interest
+      # store objective
+      Vobj = np.trace(EnKFV.getP())
+      print('post velocity trace: ', EnKFV.getP())
+      velObj.append(Vobj)
+      
+    if (droneLocation[0] == 2) or (droneLocation[0] == 7):  # if UAV at incident location
+      if time in VmaxlistofTime:
+        print("UAV at incident location when velocity is measured")
+      VmaxlistofTime.append(time)
+      EnKFV.nonLinearObs = False  # linear relationship when drone at incident location (direct uf observation)
+      if droneLocation[0] == 2:
+        EnKFV.H = np.array([[1.0, 0.0]])
+      if droneLocation[0] == 7:
+        EnKFV.H = np.array([[0.0, 1.0]])
+      # adjust observations
+      EnKFV.obsError = VdrObsSTDV
+      EnKFV.obsDim = VobsDimVf
+      VfObserved = [trueVf]
+      VmaxEnsembles = EnKFV.EnKFStep(VmaxEnsembles, VfObserved)  # the observed uf is only the one at the incident location, propagate Vmax ensembles using random walk and assimilate observed data
+      # store results
+      Vmax1Estimated.append(EnKFV.mean[0])
+      Vmax2Estimated.append(EnKFV.mean[1])
+      critDenEstimated.append(VmaxtoCritDen(EnKFV.mean))
+      # update traffic parameters
+      trafficNet.updateVmaxCritDen(EnKFV.mean, VmaxtoCritDen(EnKFV.mean))  # Now the parameters are updates for the links of interest at the locations of interest
+      # store objective
+      Vobj = np.trace(EnKFV.getP())
+      print('drone uf obs. at: ', droneLocation, EnKFV.getP())
+      print('updated ensembles: ', droneLocation, np.transpose(np.array(VmaxEnsembles))[:,0:3])
+      print(' ')
+      velObj.append(Vobj)  # store the variation in the trace of the EnKF-V (uf)
+    
+    # UAV path planning, note the objective below is NOT along candidate paths! It is the current instantaneous objective!
+    obj= (pathWeight * np.trace(EnKFV.getP()) / 2.0) + ((1-pathWeight) * np.trace(EnKFCTM.getP()) / 40.0)  # path planning objective
+    objective.append(obj)
+    # update the UAV location and update filters
+    print('pre-find path ensembles: ', np.transpose(np.array(VmaxEnsembles))[:,0:3])  # sanity check
+    explorePath = findPath(location=droneLocation, time=time, trafficNet=cp.deepcopy(trafficNet), EnKFCTM=cp.deepcopy(EnKFCTM), EnKFV=cp.deepcopy(EnKFV), weight=pathWeight)
+    droneLocation = explorePath.updateLocation()
+    print('post-find path ensembles: ', np.transpose(np.array(VmaxEnsembles))[:,0:3])  # sanity check
+    print('post-find path ensembles from EnKF: ', np.transpose(np.array(EnKFV.getUpdatedEnsembles()))[:,0:3])  # sanity check
+    EnKFCTM.droneLoc = droneLocation  # update UAV loc. in CTM EnKF
+    EnKFV.droneLoc = droneLocation  # update UAV loc. in uf EnKF
+    print(' ')
+    print(' ... ')
+    print('drone currently at: ', droneLocation)
+    storeDroneLocation.append(droneLocation)
+    droneLocCell.append(LocToCell[droneLocation])
+
+  # determine position of UAV in km from start of road
+  droneLocKm.append(cellToLength(storeDroneLocation, trafficNet))
+  
+  
+  ############# data
+  # create list of observations at critical time steps for sanity checks
+  spdict = dict()
+  for time in totalTimeSteps:
+    if time % 30 == 0:
+        spdict[time]=spData[time]
+  
+  ############# plotting
+  VmaxTrue = list()
+  for time in VmaxlistofTime:
+    VmaxTrue.append(20)
+  
+  
\ No newline at end of file
diff --git a/network.py b/network.py
new file mode 100644
index 0000000..c3afa66
--- /dev/null
+++ b/network.py
@@ -0,0 +1,240 @@
+# -*- coding: utf-8 -*-
+"""
+implements network methods
+including network loading
+
+@author: cesny
+"""
+
+import nodeModel
+import linkModel
+import numpy as np
+
+'''
+simulation parameters used in VISSIM model
+simTime = 3600
+simTimeStep = 10
+'''
+
+def linkFactory(aClass, linkID, unode, dnode, params):
+    """
+    This function is a factory
+    Classes are first class objects so they
+    can be passed into this function
+    :param aClass: A class object
+    :param pargs: parameters
+    :param kargs: key-value parameters
+    :return: initiated type specific class
+    """
+    return aClass(linkID, unode, dnode, params)
+  
+def nodeFactory(aClass, nodeID, nodeModel, fstar, rstar):
+    """
+    This function is a factory
+    Classes are first class objects so they
+    can be passed into this function
+    :param aClass: A class object
+    :param pargs: parameters
+    :param kargs: key-value parameters
+    :return: initiated type specific class
+    """
+    return aClass(nodeID, nodeModel, fstar, rstar)
+
+
+class Network:
+  """
+  This is a general network class for connecting links and nodes
+  """
+  def __init__(self, simTime, simStep, nodefile, linkfile, demandfile):
+    self.simTime = simTime  # time horizon
+    self.timeStep = simStep  # simulation time step
+    self.totalTimesteps = range(int(np.ceil(float(self.simTime)/self.timeStep)) + 1)
+    self.ODs = dict()  # create an OD dictionary, each OD is a class that stores demand
+    self.nodeDict = dict()  # dictionary of nodes
+    self.linkDict = dict()  # dictionary of links
+    self._setupNetwork(nodefile, linkfile, demandfile)
+  
+  def _setupNetwork(self, nodefile, linkfile, demandfile):
+    """
+    sets up the network
+    """
+    self.readNodes(nodefile)
+    self.readLinks(linkfile)
+    self.readDemand(demandfile)
+    return None
+  
+  def updateVmaxCritDen(self, newVmaxlist, newCritDenlist):
+    """
+    updates the critical densities in the network
+    call before calling update Vmax
+    """
+    self.linkDict[2].updateVmaxCritDen(newVmaxlist[0], newCritDenlist[0])
+    self.linkDict[7].updateVmaxCritDen(newVmaxlist[1], newCritDenlist[1])
+    return None
+  
+  
+  def readNodes(self, nfile):
+    """
+    reads node file, returns None
+    """
+    try:
+      with open(nfile) as nf:
+        print('... reading node file ...')
+        nf.readline()
+        for line in nf:
+          data = line.strip().split('\t')
+          if data[2] != '[]':
+            fstar = data[2].strip('[]').split(',')
+          else:
+            fstar = []
+          if data[3] != '[]':
+            rstar = data[3].strip('[]').split(',')
+          else:
+            rstar = []
+          self.nodeDict[int(data[0])] = nodeFactory(getattr(nodeModel, data[1]),\
+                          int(data[0]), data[1], fstar, rstar)
+          
+    except IOError as ioerr:
+      print('... error reading node file: ' + str(ioerr) + ' ...')
+    return None
+  
+  def readLinks(self, lfile):
+    """
+    lfile: link text file
+    network: network class
+    params: dictionary with model parameters, may be read from text file
+    reads links, returns None
+    """
+    if self.nodeDict == {}:
+      raise Exception('... read nodes before links ...')
+    
+    try:
+      with open(lfile) as lf:
+        print('... reading link file ...')
+        lf.readline()  # skip the first line
+        for line in lf:
+          data = line.strip().split('\t')
+          params = dict()
+          params['linkType'] = data[1]
+          params['length'] = float(data[4])
+          params['ffs'] = float(data[5])
+          params['critDen'] = float(data[6])
+          params['jamDen'] = float(data[7])
+          params['timeStep'] = self.timeStep
+          source = int(data[2])  # upstream node
+          sink = int(data[3])  # downstream node
+          sourceObject = self.nodeDict[source]
+          sinkObject = self.nodeDict[sink]
+          # create the link using params, and specify upstream downstream nodes
+          self.linkDict[int(data[0])] = linkFactory(getattr(linkModel, data[1]),\
+                          int(data[0]), sourceObject, sinkObject, params )
+    except IOError as ioerr:
+      print('... error reading link file: ' + str(ioerr) + ' ...')
+
+    self.setNodeAdjacency()
+    return None  
+  
+  def setNodeAdjacency(self):
+    """
+    puts link objects in node classes upstreamLinks and downstreamLinks
+    """
+    for nodeID in self.nodeDict:
+      node = self.nodeDict[nodeID]
+      for linkID in node.fstar:
+        node.downstreamLinks.append(self.linkDict[linkID])
+      for linkID in node.rstar:
+        node.upstreamLinks.append(self.linkDict[linkID])
+    return None
+  
+  def readDemand(self, dFile):
+    """
+    read OD file, tab delimited demand file
+    make sure time starts from zero
+    """
+    if self.nodeDict == {}:
+      raise Exception('... read nodes before links ...')
+    
+    listofOrigins = list()
+    try:
+      with open(dFile) as df:
+        print('... reading demand file ...')
+        df.readline()  # skips the first line
+        for line in df:
+          data = line.strip().split('\t')
+          time = int(data[0])  # based format of text file
+          if data[1] != '[]':
+            origins = data[1].strip('[]').split(',')
+            demand = data[2].strip('[]').split(',')
+            for key, origin in enumerate(origins):
+              self.nodeDict[int(origin)].demandRates[time] = float(demand[key])
+              if int(origin) not in listofOrigins:
+                listofOrigins.append(int(origin))
+          
+    except IOError as ioerr:
+      print('... error reading OD file: ' + str(ioerr) + ' ...')
+      
+    self._completeZoneSetup(listofOrigins)
+    return None
+  
+  def _completeZoneSetup(self, listofOrigins):
+    """
+    assigns demandRates of zero to zones when there's no flow
+    """
+    for nodeID in listofOrigins:
+      node = self.nodeDict[nodeID]
+      for time in self.totalTimesteps:
+        if time not in node.demandRates:
+          node.demandRates[time] = 0    
+    return None
+  
+  def loadNetworkStep(self, time):
+    """
+    implements the network loading algorithm for one time step
+    and returns the densities on the cells
+    """
+    for nodeID in self.nodeDict:
+      self.nodeDict[nodeID].nodeUpdate(time, self.timeStep)
+    for linkID in self.linkDict:
+      self.linkDict[linkID].linkUpdate(time + 1)
+    
+    linkDensities = list()
+    listofLinkDensities = list()
+    for linkID in self.linkDict:
+      if linkID is not 9:
+        listofDensitiesPerLink = self.linkDict[linkID].linkDensity()
+        listofLinkDensities.append(listofDensitiesPerLink)
+        for linkdenval in listofDensitiesPerLink:
+          linkDensities.append(linkdenval)
+    return linkDensities, listofLinkDensities
+  
+  def networkLoading(self):
+    """
+    the full network loading algorithm, returns dictionary of
+    dictionaries where outer keys are time, inner key is link, and value
+    is the densities on the cells
+    """
+    self.netLoadingResults = dict()
+    for time in self.totalTimesteps:
+      self.netLoadingResults[time] = self.loadNetworkStep(time)
+      
+    return self.netLoadingResults
+  
+  def resetCounts(self):
+    for linkID in self.linkDict:
+      link = self.linkDict[linkID]
+      link._upstreamCounts = dict()
+      link._downstreamCounts = dict()
+      link.inFlow = 0
+      link.outFlow = 0
+    return None
+
+if __name__ == '__main__':
+    simTime = 4490
+    simTimeStep = 10
+    totalTimeSteps = range(int(np.ceil(float(simTime)/simTimeStep)) + 1)
+    linkfile = 'VISSIMnetwork/links.txt'
+    nodefile = 'VISSIMnetwork/nodes.txt'
+    demandfile = 'VISSIMnetwork/demand6600.txt'
+    trafficNet = Network(simTime, simTimeStep, nodefile, linkfile, demandfile)
+    results=trafficNet.networkLoading()
+    print(results)
\ No newline at end of file
diff --git a/node.py b/node.py
new file mode 100644
index 0000000..5d5c170
--- /dev/null
+++ b/node.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+abstract node class
+
+@author: cesny
+"""
+
+class Node:
+  """
+  this class is for Node objects
+  """
+  def __init__(self, nodeID, nodeModel, fstar, rstar):
+    self.ID = nodeID
+    self.model = nodeModel
+    self.downstreamLinks = list()
+    self.upstreamLinks = list()
+    self.fstar = list()
+    self.rstar = list()
+    self._processStars(fstar, rstar)
+    
+  def _processStars(self, fstar, rstar):
+    """
+    creates fstar rstar from input files
+    """
+    for star in fstar:
+      intStar = int(star)
+      self.fstar.append(intStar)
+      
+    for star in rstar:
+      intStar = int(star)
+      self.rstar.append(intStar)
+    return None
+    
+  def nodeUpdate(self, time, timeStep):
+    """
+    moves flow by calculating sending/receiving flows and
+    calling the calculateTransitionFlows to get transitionFlow dict
+    this method also updates link counts
+    """
+    sendingFlow = dict()
+    receivingFlow = dict()
+    for inLink in self.upstreamLinks:
+      sendingFlow[inLink.ID] = inLink.calculateSendingFlow(time, timeStep)
+    for outLink in self.downstreamLinks:
+      receivingFlow[outLink.ID] = outLink.calculateReceivingFlow(time, timeStep)
+    # add transitiionFlows as a node attribute
+    self.transitionFlows = self.calculateTransitionFlows(sendingFlow, receivingFlow)
+    linkInflows = dict()
+    linkOutflows = dict()
+    for inLink in self.transitionFlows:
+      linkOutflows[inLink] = sum(self.transitionFlows[inLink].values())
+      for outLink in self.transitionFlows[inLink]:
+        linkInflows[outLink] = linkInflows.setdefault(outLink, 0.0) + self.transitionFlows[inLink][outLink] 
+    for inLink in self.upstreamLinks:
+      inLink.outFlow = linkOutflows[inLink.ID]
+    for outLink in self.downstreamLinks:
+      outLink.inFlow = linkInflows[outLink.ID]
+    # now make sure you use those inflows and outflows to updateLinks for time+1
+    return None
+    
+  def calculateTransitionFlows(self, sendingFlow, receivingFlow, proportions=None):
+    """
+    returns transition flow, which is a dict of dicts, keys are pair of inc
+    out links and vals are num of vehicles that can move during time interval
+    {a:{b:4, c:6}} a, b, and c are links, 4 & 6 are transitionFlows
+    """
+    pass
+  
diff --git a/nodeModel.py b/nodeModel.py
new file mode 100644
index 0000000..781c158
--- /dev/null
+++ b/nodeModel.py
@@ -0,0 +1,113 @@
+# -*- coding: utf-8 -*-
+"""
+node models
+
+@author: cesny
+"""
+
+from node import Node
+
+class Zone(Node):
+  def __init__(self, nodeID, nodeModel, fstar, rstar):
+    Node.__init__(self, nodeID, nodeModel, fstar, rstar)
+    self.__initialize()
+    
+  def __initialize(self):
+    """
+    sets subtype of node and creates
+    demandRates if origin
+    """
+    if len(self.fstar) == 0:
+      # Node is destination
+      self.subType = 'Destination'
+    if len(self.rstar) == 0:
+      self.subType = 'Origin'
+      self.demandRates = dict()  # create a demand rates dict if it's an origin!
+    if (len(self.rstar) != 0) and (len(self.fstar) != 0):
+      raise Exception('... wrong zone initialization ...')
+    if (len(self.rstar) >= 1) and (len(self.fstar) >= 1):
+      raise Exception('... code note ready for multiple ins or outs for des or orig ...') 
+    return None
+  
+  def calculateTransitionFlows(self, sendingFlow, receivingFlow, proportions=None):
+    """
+    raises exception
+    """
+    raise Exception('... load vehicles for Zones error in calculate transition flows ...')
+    return None
+  
+  def nodeUpdate(self, time, timeStep):
+    """
+    override the nodeUpdate method to update inFlows and outFlows
+    based on loading vehicles or removing them
+    """
+    if self.subType == 'Origin':
+      for outLink in self.downstreamLinks:
+        outLink.inFlow = self.demandRates[time] * (1.0/3600) * timeStep  # demand is in veh/hr, have to convert to timestep
+    if self.subType == 'Destination':
+      for inLink in self.upstreamLinks:  # if it's the destination, all the link from upstream gets out
+        inLink.outFlow = inLink.calculateSendingFlow(time, timeStep)
+    return None
+  
+    
+
+
+class SeriesNode(Node):
+  def __init__(self, nodeID, nodeModel, fstar, rstar):
+    Node.__init__(self, nodeID, nodeModel, fstar, rstar)
+  
+  def calculateTransitionFlows(self, sendingFlow, receivingFlow, proportions=None):
+    """
+    sendingFlow and receivingFlow are dictionaries
+    sendingFlow = {uplink1: val, uplink2: val..}
+    """
+    transitionFlows = dict()
+    for inLinkID in self.rstar:
+      transitionFlows[inLinkID] = dict()
+    
+    upLink = self.rstar[0]
+    downLink = self.fstar[0]
+    transitionFlows[upLink][downLink] = min(sendingFlow[upLink], receivingFlow[downLink])
+    return transitionFlows
+  
+
+class DivergeNode(Node):
+  def __init__(self, nodeID, nodeModel, fstar, rstar):
+    self.proportions = dict()  # a dictionary with props of flow to links, considering that proportions are fixed
+    Node.__init__(self, nodeID, nodeModel, fstar, rstar)
+  
+  def _processStars(self, fstar, rstar):
+    """
+    overrides node processStars to create proportions if they are static
+    """
+    for star in rstar:
+      intStar = int(star)
+      self.rstar.append(intStar)
+    self.proportions[self.rstar[0]] = dict()
+    for star in fstar:
+      link, prop = star.split(':')
+      intLink = int(link)
+      floatProp = float(prop)
+      self.fstar.append(intLink)
+      self.proportions[self.rstar[0]][intLink] = floatProp    
+    return None
+  
+  def calculateTransitionFlows(self, sendingFlow, receivingFlow, proportions=None):
+    """
+    sendingFlow and receivingFlow are dictionaries
+    sendingFlow = {uplink1: val, uplink2: val..}
+    """
+    transitionFlows = dict()
+    for inLinkID in self.rstar:
+      transitionFlows[inLinkID] = dict()
+    inLink = self.rstar[0]
+    thetas = list()
+    for outLink in self.fstar:
+      if (sendingFlow[inLink] != 0) and (self.proportions[inLink][outLink] != 0):
+        thetas.append(receivingFlow[outLink] / ((self.proportions[inLink][outLink])*(sendingFlow[inLink])))
+    thetas.append(1.0)
+    theta = min(thetas)
+    for outLink in self.fstar:
+      transitionFlows[inLink][outLink] = theta * self.proportions[inLink][outLink] * sendingFlow[inLink]
+    return transitionFlows
+
diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000..4dfdb14
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,171 @@
+# -*- coding: utf-8 -*-
+"""
+this contains utilities
+
+@author: cesny
+"""
+from network import Network
+import numpy as np
+import copy as cp
+
+
+def readData(textfile):
+  """
+  reads the text file
+  returns two lists with density and velocity
+  measured at a specific point
+  """
+  try:
+    with open(textfile) as tf:
+      density = dict()
+      speed = dict()
+      for line in tf:
+        data = line.strip().split(';')
+        first_number = 0
+        try:
+          first_number = float(data[0])
+        except:
+          pass
+        if first_number == 1:
+          roadid_str_list = data[2].split('-')
+          roadid = (float(roadid_str_list[0]), float(roadid_str_list[1]), \
+                      float(roadid_str_list[2]))
+          time = data[1].split('-')
+          timeStep = float(time[0]) / 10
+          density.setdefault(timeStep,[])
+          speed.setdefault(timeStep,[])
+          thresholds = [120.0, 370.0, 620.0, 870.0, 1120.0, 1370.0]  # thresholds are used so that you take a measurement in each cell
+          if (roadid[0] != 9) and (roadid[1] in thresholds):
+            density[timeStep].append(float(data[3]))
+            if (roadid[0] == 2 and roadid[1] == 120.0) or (roadid[0] == 7 and roadid[1] == 120.0):  # storing first cell only
+              if (float(data[3]) != 0.0):  #density is not zero
+                speed[timeStep].append(float(data[5]))
+              else:
+                speed[timeStep].append(100.0)
+  except IOError as ioerr:
+    print('failed to read' + str(ioerr))
+  return density, speed
+
+
+def setCTMVehicles(trafficNet, EnKFensemble):
+  '''
+  sets the vehicles in CTM traffic model
+  based on one particular ensemble output
+  an ensemble is a list of 40 CTM densities
+  across cells, set ramp cells to zero to
+  allow vehicles to exit
+  '''
+  cellindex = 0
+  for linkID in trafficNet.linkDict:
+    if linkID is not 9:
+      for cell in trafficNet.linkDict[linkID].cells:
+        EnKFveh = EnKFensemble[cellindex] * cell.length  # convert density to vehicles
+        cell.vehicles = EnKFveh  # assign EnKF determined vehicles to cell
+        cellindex += 1
+    if linkID is 9:
+      for cell in trafficNet.linkDict[linkID].cells:
+        cell.vehicles = 0
+  return trafficNet
+
+
+def forwardCTMPropagation(time, trafficNet, EnKFensembles):
+  '''
+  propagates a set of EnKF ensembles forward
+  each EnKF ensemble is a list of 40 density values across
+  cells
+  flow is propagated using cell transmission model
+  return a list of lists with propagated ensembles
+  '''
+  propagatedEnsembles = list()
+  for ensemble in EnKFensembles:
+    trafficNet = setCTMVehicles(trafficNet, ensemble)  # set current ensembles
+    propagatedEnsembles.append(trafficNet.loadNetworkStep(time)[0])  # propagate the ensembles forward using model
+  return propagatedEnsembles
+
+
+def CTMcreateInitialEnsemble(stateDim, ensembleSize, modSTDV, bestguess=20):
+  '''
+  creates an initial ensemble around a best guess
+  estimate of the state
+  best guess is a single number for best guess average
+  densities on the states
+  '''
+  return np.random.normal(loc=bestguess, scale=modSTDV, size=(ensembleSize,stateDim)).tolist()
+
+  
+def VmaxCreateInitialEnsemble(VstateDim,Vensembles,VmodSTDV, bestguess=80):
+  '''
+  creates an initial ensemble for vmax
+  '''
+  return np.random.normal(loc=bestguess, scale=VmodSTDV, size=(Vensembles,VstateDim)).tolist()
+
+  
+def m(vmax, rho):
+  '''
+  this function is used to get the model predicted
+  velocity measurement from vmax (i.e. model predicted diagnostic variable
+  for nonlinear EnKF observation matrix)
+  original back wave based on vmax=100, rhocr=80, jamDen = 300 - can be
+  fitted by measuring traffic considtion under normal conditions
+  '''
+  v=np.full(vmax.shape,0.0)
+  for key, vmaxVal in enumerate(vmax):  
+    rhoCritical = (80.0*100*300)/(vmaxVal*(300 - 80) + 80*100)  # update rho critical given current vmax
+    if rho < rhoCritical:
+      v[key] = vmaxVal
+    elif rho > rhoCritical:
+      v[key] = float(vmaxVal*(rhoCritical)*(300.0 - rho))/(rho*(300.0 - rhoCritical))
+  return v  # generate the predicted observation, this will be used to update vmax
+
+
+def VmaxtoCritDen(vmax):
+  '''
+  defines the relationship based on maintaining uncongested
+  backwave
+  '''
+  rhoCrit=np.full(vmax.shape,0.0)
+  for key, vmaxVal in enumerate(vmax):  
+    rhoCritical = (80.0*100*300)/(vmaxVal*(300 - 80) + 80*100)
+    rhoCrit[key] = rhoCritical
+  return rhoCrit
+
+
+def createLocToCell(trafficNet):
+  '''
+  creates the mapping between (linkid, cell) to
+  global cellindex between (0,40) across all links
+  '''
+  cellindex = 0
+  locToCell = dict()
+  cellToLoc = dict()
+  for linkID in trafficNet.linkDict:
+    if linkID is not 9:
+      for cellkey, cell in enumerate(trafficNet.linkDict[linkID].cells):
+        locToCell[(linkID, cellkey)] = cellindex
+        cellToLoc[cellindex] = (linkID, cellkey)
+        cellindex += 1
+  return locToCell
+
+
+def cellToLength(listofLocations, trafficNet):
+  '''
+  returns UAV location on the road network in km
+  across time, input list of UAV location (linkID, linkCell)
+  '''
+  DronePositions = list()
+  LocToCell = createLocToCell(trafficNet)
+  for location in listofLocations:
+    cellLocation = LocToCell[location]
+    physicalLocation = (5.0/18)*cellLocation + (5.0/36)
+    DronePositions.append(physicalLocation)
+  return DronePositions
+
+
+def lengthToCell(locKm):
+  '''
+  from position in km to cell number
+  '''
+  cell_loc = list()
+  for loc in locKm:
+     cell_loc.append(int(loc*(18.0/5)-0.5))
+  return cell_loc
\ No newline at end of file