File non_monotonic_forecasting.py

Michel Bierlaire, EPFL Fri Jul 25 2025, 17:27:35

Forecasting with a MDCEV model and the “non-monotonic utility” specification.

Example: non monotonic utility
Forecasting observation 0 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '8.52', 2: '163', 3: '322', 4: '5.73'} objective 210, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '8.52', 2: '163', 3: '322', 4: '5.73'} objective 210, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '6.16', 2: '223', 3: '265', 4: '5.34'} objective 249, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '6.16', 2: '223', 3: '265', 4: '5.35'} objective 249, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '17.9', 2: '299', 3: '173', 4: '10.2'} objective 188, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '17.9', 2: '299', 3: '173', 4: '10.2'} objective 188, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '46.1', 2: '320', 3: '128', 4: '6.18'} objective 171, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '46.1', 2: '320', 3: '128', 4: '6.17'} objective 171, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '16.7', 2: '212', 3: '253', 4: '18.4'} objective 168, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '16.7', 2: '212', 3: '253', 4: '18.4'} objective 168, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '14.3', 2: '372', 3: '109', 4: '5.37'} objective 187, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '14.3', 2: '372', 3: '109', 4: '5.37'} objective 187, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '52.5', 2: '124', 3: '30.1', 4: '294'} objective 195, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '52.5', 2: '123', 3: '30.1', 4: '294'} objective 195, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '18', 2: '324', 3: '152', 4: '5.9'} objective 221, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '18', 2: '324', 3: '152', 4: '5.9'} objective 221, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '35.7', 2: '291', 3: '169', 4: '5.1'} objective 227, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '35.7', 2: '291', 3: '169', 4: '5.1'} objective 227, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '1.38', 2: '12', 3: '484', 4: '3.1'} objective 420, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '1.38', 2: '12', 3: '484', 4: '3.1'} objective 420, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 1 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '1.82', 2: '121', 3: '17', 4: '360'} objective 255, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '1.83', 2: '121', 3: '17', 4: '360'} objective 255, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '59.8', 2: '348', 3: '83.7', 4: '8.72'} objective 166, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '59.4', 2: '348', 3: '84.3', 4: '8.71'} objective 166, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '14.5', 2: '57.9', 3: '423', 4: '4.11'} objective 257, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '14.5', 2: '57.9', 3: '423', 4: '4.11'} objective 257, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '2.63', 2: '470', 3: '22.5', 4: '4.52'} objective 270, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '2.63', 2: '470', 3: '22.5', 4: '4.51'} objective 270, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '18.8', 2: '43.5', 3: '429', 4: '8.41'} objective 227, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '18.8', 2: '43.5', 3: '429', 4: '8.41'} objective 227, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '16.4', 2: '437', 3: '44.1', 4: '2.93'} objective 254, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '16.4', 2: '437', 3: '44.1', 4: '2.93'} objective 254, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '11.6', 2: '324', 3: '156', 4: '8.26'} objective 220, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '11.6', 2: '324', 3: '156', 4: '8.26'} objective 220, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '14.3', 2: '389', 3: '93.4', 4: '3.42'} objective 258, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '14.3', 2: '389', 3: '93.4', 4: '3.41'} objective 258, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '9.09', 2: '401', 3: '84.1', 4: '5.79'} objective 218, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '9.09', 2: '401', 3: '84.1', 4: '5.79'} objective 218, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '19.9', 2: '352', 3: '119', 4: '9.42'} objective 163, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '19.9', 2: '352', 3: '119', 4: '9.41'} objective 163, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with brute force algorithm: 20.2 seconds
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with analytical algorithm: 1.03 seconds
                 1             2            3            4
count  2000.000000  2.000000e+03  2000.000000  2000.000000
mean     39.598668  3.267214e+02   122.819661    10.860305
std      79.603989  1.411482e+02   125.865887    31.240736
min       0.000000  3.640836e-15     0.000000     0.000000
25%       6.887531  2.312174e+02    31.750637     4.066402
50%      13.874335  3.747594e+02    70.409042     5.882995
75%      30.098289  4.426312e+02   172.528683     8.805091
max     497.780408  5.000000e+02   500.000000   482.881548
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     39.602338   326.719112   122.818632    10.859918
std      79.609622   141.150621   125.866995    31.243921
min       0.000000     0.000000     0.000000     0.000000
25%       6.886657   231.826845    31.749006     4.065986
50%      13.873640   374.762290    70.406599     5.882881
75%      30.103048   442.632262   172.489746     8.807932
max     497.783463   500.000000   500.000000   482.850818
                 1             2            3            4
count  2000.000000  2.000000e+03  2000.000000  2000.000000
mean     38.561604  3.233693e+02   127.299435    10.769694
std      78.690955  1.438931e+02   128.870804    31.042549
min       0.000000  1.370070e-15     0.000000     0.000000
25%       6.448176  2.145525e+02    31.190340     3.906875
50%      12.450833  3.725790e+02    73.637691     5.842609
75%      27.639689  4.437999e+02   187.630390     8.913430
max     487.783138  5.000000e+02   499.616097   486.839276
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     38.561587   323.361076   127.305445    10.771892
std      78.693250   143.900485   128.877816    31.050409
min       0.000000     0.000000     0.000000     0.000000
25%       6.449612   214.741469    31.189704     3.908586
50%      12.449893   372.572041    73.635679     5.842149
75%      27.649573   443.800831   187.630687     8.919349
max     487.784147   500.000000   499.616172   486.818113

import sys
import time

import numpy as np
import pandas as pd
from IPython.display import display

import biogeme.biogeme_logging as blog
from biogeme.database import Database
from biogeme.results_processing import EstimationResults
from non_monotonic_specification import the_non_monotonic
from process_data import database

logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example: non monotonic utility')

result_file = 'saved_results/non_monotonic.yaml'
try:
    results = EstimationResults.from_yaml_file(filename=result_file)
except FileNotFoundError as e:
    print(e)
    print(f'File {result_file} is missing.')
    sys.exit()

the_non_monotonic.estimation_results = results

# %
# We apply the model only on the first two rows of the database.
two_rows_of_database: Database = database.extract_rows([10, 11])
# %
budget_in_hours = 500

# %
# # Validation

# %
# As the implementation is still experimental, we compare the result obtained by the bruteforce algorithm and
# the analytical algorithm for a few draws.

# Note that minor discrepancies between the outcome of the two algorithms are likely to occur, due to numerical
# imprecision, inevitable in finite arithmetic.

# However, if there are major differences, it should be reported.

# %
number_of_draws = 10

# %
# We generate the draws
epsilons = [
    np.random.gumbel(
        loc=0, scale=1, size=(number_of_draws, the_non_monotonic.number_of_alternatives)
    )
    for _ in range(two_rows_of_database.num_rows())
]

# %
# We first compare the results obtained from the brute force and the analytical algorithms, for each draw.
the_non_monotonic.validate_forecast(
    database=two_rows_of_database, total_budget=budget_in_hours, epsilons=epsilons
)

# %
# # Forecasting
# We use a larger number of draws to obtain the forecast.

# %
number_of_draws = 2000

# %
# We generate the draws
epsilons = the_non_monotonic.generate_epsilons(
    number_of_observations=two_rows_of_database.num_rows(),
    number_of_draws=number_of_draws,
)
# %
# First, the brute force algorithm.
start_time = time.time()
optimal_consumptions_brute_force: list[pd.DataFrame] = the_non_monotonic.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=True,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with brute force algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Then, the analytical algorithm.
start_time = time.time()
optimal_consumptions_analytical: list[pd.DataFrame] = the_non_monotonic.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=False,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with analytical algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Results for the first observation, brute force method
display(optimal_consumptions_brute_force[0].describe())

# %
# Results for the first observation, analytical method
display(optimal_consumptions_analytical[0].describe())

# %
# Results for the second observation, brute force method
display(optimal_consumptions_brute_force[1].describe())

# %
# Results for the second observation, analytical method
display(optimal_consumptions_analytical[1].describe())