In [ ]:
#%%
"""File 19individualLevelParameters

:author: Michel Bierlaire, EPFL
:date: Wed Aug 26 14:56:49 2020

Calculation of the individual level parameters for model 05normalMixture
"""

import os
import pickle
import pandas as pd
import biogeme.database as db
import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import Beta, DefineVariable, bioDraws, MonteCarlo

# Read the data
df = pd.read_csv('swissmetro.dat', sep='\t')
database = db.Database('swissmetro', df)

# The following statement allows you to use the names of the variable
# as Python variable.
globals().update(database.variables)

# Here we use the "biogeme" way for backward compatibility
exclude = ((PURPOSE != 1) * (PURPOSE != 3) + (CHOICE == 0)) > 0
database.remove(exclude)

# Values of the parameters estimated by the model 05normalMixture
ASC_CAR = Beta('ASC_CAR', 0.137, None, None, 0)
ASC_TRAIN = Beta('ASC_TRAIN', -0.402, None, None, 0)
ASC_SM = Beta('ASC_SM', 0, None, None, 1)
B_COST = Beta('B_COST', -1.28, None, None, 0)

# Define a random parameter, normally distributed, designed to be used
# for Monte-Carlo simulation
B_TIME = Beta('B_TIME', -2.26, None, None, 0)

# It is advised not to use 0 as starting value for the following parameter.
B_TIME_S = Beta('B_TIME_S', 1.65, None, None, 0)
B_TIME_RND = B_TIME + B_TIME_S * bioDraws('B_TIME_RND', 'NORMAL')

# Definition of new variables
SM_COST = SM_CO * (GA == 0)
TRAIN_COST = TRAIN_CO * (GA == 0)

# Definition of new variables: adding columns to the database
CAR_AV_SP = DefineVariable('CAR_AV_SP', CAR_AV * (SP != 0), database)
TRAIN_AV_SP = DefineVariable('TRAIN_AV_SP', TRAIN_AV * (SP != 0), database)
TRAIN_TT_SCALED = DefineVariable('TRAIN_TT_SCALED', TRAIN_TT / 100.0, database)
TRAIN_COST_SCALED = DefineVariable(
    'TRAIN_COST_SCALED', TRAIN_COST / 100, database
)
SM_TT_SCALED = DefineVariable('SM_TT_SCALED', SM_TT / 100.0, database)
SM_COST_SCALED = DefineVariable('SM_COST_SCALED', SM_COST / 100, database)
CAR_TT_SCALED = DefineVariable('CAR_TT_SCALED', CAR_TT / 100, database)
CAR_CO_SCALED = DefineVariable('CAR_CO_SCALED', CAR_CO / 100, database)

# Definition of the utility functions
V1 = ASC_TRAIN + B_TIME_RND * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED
V2 = ASC_SM + B_TIME_RND * SM_TT_SCALED + B_COST * SM_COST_SCALED
V3 = ASC_CAR + B_TIME_RND * CAR_TT_SCALED + B_COST * CAR_CO_SCALED

# Associate utility functions with the numbering of alternatives
V = {1: V1, 2: V2, 3: V3}

# Associate the availability conditions with the alternatives
av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}

# Conditional to B_TIME_RND, we have a logit model (called the kernel)
prob_chosen = models.logit(V, av, CHOICE)

numerator = MonteCarlo(B_TIME_RND * prob_chosen)
denominator = MonteCarlo(prob_chosen)

simulate = {
    'Numerator': numerator,
    'Denominator': denominator,
    'Choice': CHOICE,
}

pickleFile = '19individualLevelParameters.pickle'
if os.path.isfile(pickleFile):
    with open(pickleFile, 'rb') as f:
        sim = pickle.load(f)
else:
    biosim = bio.BIOGEME(database, simulate, numberOfDraws=100000)
    sim = biosim.simulate()
    sim['Individual-level parameters'] = sim['Numerator'] / sim['Denominator']
    with open(pickleFile, 'wb') as f:
        pickle.dump(sim, f)

htmlFile = '19individualLevelParameters.html'
with open(htmlFile, 'w') as h:
    print(sim.to_html(), file=h)
print(f'Simulated values available in {htmlFile}')