"""File netherlands_model.py

Model used for aggregation and forecasting

Michel Bierlaire
Fri Jul  1 09:54:05 2022

"""

import biogeme.database as db
import pandas as pd
from biogeme.expressions import Beta, Variable, log, exp

df = pd.read_csv('netherlands.dat', sep='\t')
database = db.Database('netherlands', df)

sp = Variable('sp')
rail_ivtt = Variable('rail_ivtt')
rail_acc_time = Variable('rail_acc_time')
rail_egr_time = Variable('rail_egr_time')
car_ivtt = Variable('car_ivtt')
car_walk_time = Variable('car_walk_time')
car_cost = Variable('car_cost')
rail_cost = Variable('rail_cost')
choice = Variable('choice')

purpose = Variable('purpose')
business = purpose == 1
other = purpose == 0

gender = Variable('gender')
female = gender == 1
male = gender == 0

exclude = sp != 0
database.remove(exclude)

rail_time = rail_ivtt + rail_acc_time + rail_egr_time
car_time = car_ivtt + car_walk_time

DUTCH_GUILDERS_TO_EUROS = 0.44378022
car_cost_euro = car_cost * DUTCH_GUILDERS_TO_EUROS
rail_cost_euro = rail_cost * DUTCH_GUILDERS_TO_EUROS

beta_time_car = Beta('beta_time_car', 0, None, None, 0)

asc_car_business = Beta('asc_car_business', 0, None, None, 0)
beta_cost_business = Beta('beta_cost_business', 0, None, None, 0)
beta_time_rail_business = Beta('beta_time_rail_business', 0, None, None, 0)

asc_car_other = Beta('asc_car_other', 0, None, None, 0)
beta_cost_other = Beta('beta_cost_other', 0, None, None, 0)
beta_time_rail_other = Beta('beta_time_rail_other', 0, None, None, 0)

asc_car_male = Beta('asc_car_male', 0, None, None, 0)
asc_car_female = Beta('asc_car_female', 0, None, None, 0)

asc_car = (
    business * asc_car_business
    + other * asc_car_other
    + male * asc_car_male
    + female * asc_car_female
)

beta_cost = business * beta_cost_business + other * beta_cost_other
beta_time_rail = business * beta_time_rail_business + other * beta_time_rail_other

v_car = asc_car + beta_cost * car_cost_euro + beta_time_car * car_time
v_rail = beta_cost * rail_cost_euro + beta_time_rail * rail_time
prob_car = 1 / (1 + exp(v_rail - v_car))
prob_rail = 1 - prob_car
prob_observation = prob_car * (choice == 0) + prob_rail * (choice == 1)
logprob = log(prob_observation)