<repo_name>jackdreilly/commroute<file_sep>/src/ctm.py
from cr_utils import Dumpable
from cr_network import Link
from demand import FlowNetwork

class FundamentalDiagram(Dumpable):
  """docstring for FundamentalDiagram"""

  def __init__(self, v, w, rho_max, q_max):
    super(FundamentalDiagram, self).__init__()
    self.v = float(v)
    self.w = float(w)
    self.rho_max = float(rho_max)
    self.q_max = float(q_max)

  def jsonify(self):
    """docstring for jsonify"""
    return {
      'v': self.v,
      'w': self.w,
      'rho_max': self.rho_max,
      'q_max': self.q_max
    }

  def rho_cong(self, flow):
    return self.rho_max - float(flow)/self.w

  def flow_cong(self, rho):
    return self.w * (self.rho_max - rho)

  def rho_ff(self, flow):
    return float(flow) / self.v

  def flow_ff(self, rho):
    return self.v * rho

  def flow(self, rho):
    return min(self.q_max, self.flow_cong(rho), self.flow_ff(rho))

  @classmethod
  def load_with_json_data(cls, data):
    return cls(
      v=data['v'],
      w=data['w'],
      q_max=data['q_max'],
      rho_max=data['rho_max']
    )


class CTMLink(Link):
  """docstring for CTMLink"""

  def __init__(self, l, fd, flow=0.0, *args, **kwargs):
    super(CTMLink, self).__init__(*args, **kwargs)
    self.l = l
    self.fd = fd
    self.flow = flow

  def jsonify(self):
    json = super(CTMLink, self).jsonify()
    json.update({
      'l': self.l,
      'fd': self.fd.jsonify(),
      'flow': self.flow
    })
    return json


  def d3_value(self):
    return self.l


class DensityCTMLink(CTMLink):
  """docstring for DensityCTMLink"""

  def __init__(self, rho=0.0, *args, **kwargs):
    super(DensityCTMLink, self).__init__(*args, **kwargs)
    self.rho = rho



  def jsonify(self):
    """docstring for jsonify"""
    json = super(DensityCTMLink, self).jsonify()
    json.update({
      'rho': self.rho
    })
    return json

  @classmethod
  def additional_kwargs(cls, data):
    return dict(
      l=data['l'],
      rho=data['rho'],
      flow=data['flow'],
      fd=FundamentalDiagram.load_with_json_data(data['fd'])
    )


class DensityCTMNetwork(FlowNetwork):
  link_class = DensityCTMLink

  def __init__(self):
    super(DensityCTMNetwork, self).__init__()


  def tt_free_flow(self, route):
    return sum(link.l / link.fd.v for link in route.get_links)

if __name__ == '__main__':
  print 'hei'
  <file_sep>/src/static_op.py
from cr_optimize import SimpleOptimizeMixIn
from ctm import *
from cvxpy import variable, eq, geq
from demand import ODDemand, RouteDemand

class LagrangianStaticProblem(FlowNetwork, SimpleOptimizeMixIn):
  """docstring for StaticProblem"""

  def cvxify(self):
    """docstring for cvxify"""
    self.cache_props()
    for link in self.get_links():
      link.v_flow = variable(name='flow: {0}'.format(link.name))
    for (source, sink), routes in self.od_routes.iteritems():
      for i, route in enumerate(routes):
        route.v_flow = variable(name='rf: o: {0}, d: {1} [{2}]'.format(source.name, sink.name, i))

class LagrangianConstrained(LagrangianStaticProblem):

  def con_junc(self):
    return [eq(sum(link.v_flow for link in junction.in_links),
               sum(link.v_flow for link in junction.out_links))
            for junction in self.junctions]

  def con_od_flows(self):
    od_demands = filter(lambda dem: isinstance(dem, ODDemand), self.demands)
    route_demands = filter(lambda dem: isinstance(dem, RouteDemand), self.demands)

    def route_flows(link):
      flow = 0.0
      link_routes = link.routes()
      for route in link_routes:
        flow += route.v_flow
      for dem in route_demands:
        if dem.route in link_routes:
          flow += dem.flow
      return flow

    def od_flows(source, sink):
      flow = 0
      for route in self.od_routes[source, sink]:
        flow += route.v_flow
      return flow

    return [
           eq(route_flows(link), link.v_flow)
           for link in self.get_links()
           ] + [
    eq(od_flows(dem.source, dem.sink), dem.flow)
    for dem in od_demands
    ] + [
      geq(route.v_flow, 0.0) for route in self.all_routes()
    ]

  def constraints(self):
    return super(LagrangianConstrained, self).constraints() + self.con_od_flows()<file_sep>/src/static_ctm.py
from cvxpy import geq, leq, quad_over_lin, variable, hstack, minimize, eq
from cvxpy import max as cvx_max
from cr_optimize import SimpleOptimizeMixIn
from cr_utils import flatten
from ctm import DensityCTMNetwork
from static_op import LagrangianConstrained
from numpy import inf

__author__ = 'jdr'


class CTMStaticProblem(DensityCTMNetwork, LagrangianConstrained):
  def cvxify(self):
    super(CTMStaticProblem, self).cvxify()
    self.cache_props()
    for link in self.get_links():
      link.v_dens = variable(name='dens: {0}'.format(link.name))
      # this next line isn't useful right now
      # link.d3_value = lambda: link.v_dens.value


class CTMConstrained(CTMStaticProblem):
  def con_ctm(self):
    return list(flatten([geq(link.v_flow, 0),
                         leq(link.v_flow, link.fd.q_max),
                         leq(link.v_flow, link.fd.v * link.v_dens),
                         leq(link.v_flow, link.fd.w * (link.fd.rho_max - link.v_dens)),
                         geq(link.v_dens, 0),
                         leq(link.v_dens, link.fd.rho_max),
                         ] for link in self.get_links()))

  def constraints(self):
    return super(CTMConstrained, self).constraints() + self.con_ctm()

  def check_feasible_constraints(self):
    return list(flatten(
      [geq(link.v_flow, link.flow), geq(link.v_dens, link.rho),leq(link.v_flow, link.flow), leq(link.v_dens, link.rho)]
        for link in self.get_links()
    ))

  def check_feasible(self):
    class FeasibleProgram(SimpleOptimizeMixIn):
      outer = self
      def constraints(self):
        return self.outer.constraints() + self.outer.check_feasible_constraints()

      def cvxify(self):
        self.outer.cvxify()

    op = FeasibleProgram()
    prog  = op.get_program()
    prog.show()
    return not (prog.solve(quiet=True) == inf)

class ComplianceConstrained(CTMConstrained):
  def route_tt_heuristic(self, route):
    def link_tt_heuristic(link):
      ff = link.l / link.fd.v
      q_max = (link.l / link.fd.q_max) * link.v_dens
      rho_hat = link.fd.rho_max - link.v_dens
      cong = link.l / link.fd.w * (quad_over_lin(link.fd.rho_max ** .5, rho_hat) - 1)
      return cvx_max(hstack([ff, q_max, cong]))

    return sum(map(link_tt_heuristic, route.links))

  def route_tt_real(self, route):
    def link_tt(link):
      if link.v_dens.value < link.fd.rho_max * 10e-8:
        return link.l / link.fd.v
      return link.l * link.v_dens.value / link.v_flow.value

    return sum([
    link_tt(link)
    for link in route.get_links
    ])

  def con_route_tt(self):
    return [
    leq(self.route_tt_heuristic(route), 10000.0)
    for route in self.all_routes()
    ]

  def constraints(self):
    return super(ComplianceConstrained, self).constraints() + self.con_route_tt()


class MinTTT(CTMStaticProblem):
  def objective(self):
    return minimize(sum(link.l * link.v_dens for link in self.get_links()))


class MinTTTComplianceProblem(MinTTT, ComplianceConstrained):
  def __init__(self):
    super(MinTTTComplianceProblem, self).__init__()

  def objective(self):
    return MinTTT.objective(self)

  def constraints(self):
    return ComplianceConstrained.constraints(self)


class MinTTTLagrangianCTMProblem(MinTTT, CTMConstrained):
  def __init__(self):
    super(MinTTTLagrangianCTMProblem, self).__init__()

  def objective(self):
    return MinTTT.objective(self)

  def constraints(self):
    return CTMStaticProblem.constraints(self)<file_sep>/src/cr_utils.py
__author__ = 'jdr'
from simplejson import dump, load
from itertools import chain

class Dumpable(object):
  def dump(self, fn):
    """docstring for dump"""
    with open(fn, 'w') as fn:
      dump(self.jsonify(), fn, indent=2)

  def jsonify(self):
    raise NotImplementedError('need to jsonify!')

  @classmethod
  def load(cls, fn, **kwargs):
    with open(fn, 'r') as fn:
      data = load(fn)
    return cls.load_with_json_data(data, **kwargs)

  @classmethod
  def load_with_json_data(cls, data, **kwargs):
    raise NotImplementedError('need to load object!')


def flatten(listOfLists):
  "Flatten one level of nesting"
  return chain.from_iterable(listOfLists)
<file_sep>/src/experiments.py
from demand import RouteDemand, ODDemand

__author__ = 'jdr'
from static_ctm import *
from ctm import FundamentalDiagram
from ctm import DensityCTMLink
from cr_network import Junction

def exp_1_create():
  """
  create the base network
  """
  v = 1.0
  w = .5
  rho_max = 3
  q_max = 1.0
  l = 1.0

  fd = FundamentalDiagram(
    v=v,
    w=w,
    rho_max=rho_max,
    q_max=q_max
  )

  net = DensityCTMNetwork()
  source = DensityCTMLink(
    net=net,
    name='source',
    l=l,
    fd=fd,
    rho=1.0
  )
  left = DensityCTMLink(
    net=net,
    name='left',
    l=l,
    fd=fd,
    rho=1.0
  )
  right = DensityCTMLink(
    net=net,
    name='right',
    l=2 * l,
    fd=fd,
    rho=1.0
  )
  sink = DensityCTMLink(
    net=net,
    name='sink',
    l=l,
    fd=fd,
    rho=1.0
  )
  junctions = [
    Junction([source], [left, right]),
    Junction([left, right], [sink]),
    ]
  for junction in junctions:
    net.add_junction(junction)

  net.dump("networks/exps/exp1/net.json")


def exp_1_demands():
  net = DensityCTMNetwork.load('networks/exps/exp1/net.json')
  left = net.route_by_names([
    'source', 'left', 'sink'
  ])
  route_demand = RouteDemand(left, .2)
  od_demand = ODDemand(net.link_by_name('source'), net.link_by_name('sink'), .8)
  net.demands.extend([route_demand, od_demand])
  net.dump('networks/exps/exp1/net_w_demand.json')


def exp_1_opt():
  '''
  @return: nothing, demonstrate

  show that you can add constraints on the fly, and that this increases the objective when you force large densities
  '''
  net = MinTTTComplianceProblem.load('networks/exps/exp1/net_w_demand.json')
  prog = net.get_program()
  prog.show()
  prog.solve()
  print 'obnjective', prog.objective.value
  left = net.link_by_name('left')
  prog.constraints.extend([
    geq(left.v_dens, left.fd.v * left.fd.q_max * 1.2)
  ])
  prog.show()
  prog.solve()
  print 'obnjective', prog.objective.value
  for link in net.get_links():
    link.flow = link.v_flow.value
    link.rho = link.v_dens.value
  net.d3ize()
  print net.check_feasible()


def exp_1_nash():
  net = MinTTTComplianceProblem.load('networks/exps/exp1/net_w_demand.json')
  left = net.link_by_name('left')
  right = net.link_by_name('right')
  left.flow = .75
  left.rho = left.fd.rho_cong(left.flow)
  right.flow = 1 - left.flow
  right.rho = right.fd.rho_ff(right.flow)
  source = net.link_by_name('source')
  sink = net.link_by_name('sink')
  source.flow = 1.0
  source.rho = source.fd.rho_ff(source.flow)
  sink.flow = 1.0
  sink.rho = source.fd.rho_ff(sink.flow)
  net.dump('networks/exps/exp1/net_nash.json')

def exp_1_nash_feasible():
  net = MinTTTComplianceProblem.load('networks/exps/exp1/net_nash.json')
  print net.check_feasible()

def main():
#  exp_1_create()
#  exp_1_demands()
#  exp_1_opt()
#  exp_1_nash()
  exp_1_nash_feasible()

if __name__ == '__main__':
  main()

<file_sep>/src/cr_optimize.py
from cr_utils import Dumpable
from cr_network import CRNetwork
from demand import Demand
from cvxpy import program, minimize

__author__ = 'jdr'

class OptimizeMixIn(object):

  def constraints(self):
    raise NotImplementedError("implement me")

  def objective(self):
    raise NotImplementedError("implement me")

  def cvxify(self):
    pass

  def get_program(self):
    self.cvxify()
    return program(self.objective(), self.constraints())

  def solve_problem(self):
    program = self.get_program()
    program.solve()
    return program


class Feasible:

  def objective(self):
    return minimize(0.0)


class Unconstrained:

  def constraints(self):
    return []

class SimpleOptimizeMixIn(Feasible, Unconstrained, OptimizeMixIn):
  pass