BAMresearch / probeye

A general framework for setting up parameter estimation problems.
MIT License
5 stars 5 forks source link

Output sensor order in forward model #103

Closed ppp0l closed 1 year ago

ppp0l commented 1 year ago

Hello everyone, I recently noticed an unexpected behaviour when using the scipy/MC solvers. Summarising: when defining the forward model class, one needs to have a consistent order for the output sensor in the interface() method and in the response() output dictionary; otherwise the solver runs, but gives completely wrong results. This is not specified anywhere in the documentation and I find it counterintuitive: having a dictionary, I would not expect the order of the keys to play any role.

I briefly discussed this with @danielandresarcones and we found that instead of using the dictionary structure, the scipy solver transforms the response() output dictionary directly into a numpy array without any checks. Being both interface() and response() user-defined, there is no guarantee of order consistency, but there is no warning/error for this. We also noted that the ordering of the experimental data dictionary does not play any role, since it is reordered upon being given to the forward model. So the experimental data ordering is always consistent with the sensor definition order.

Given this, we thought of three possible solutions:

What do you think? Do you have any other ideas?

Below is a python script as an example: you have two forward models, identical but for the output sensor order in the response() output dictionary, who behave differently. The problem for the second forward model lies in the inconsistency in the ordering of the sensors between interface() and response().

import numpy as np
import math 

from probeye.definition.inverse_problem import InverseProblem
from probeye.definition.distribution import Uniform
from probeye.definition.forward_model import ForwardModelBase
from probeye.definition.sensor import Sensor
from probeye.definition.likelihood_model import GaussianLikelihoodModel

from probeye.inference.scipy.solver import MaxLikelihoodSolver

# dummy linear model
def linearModel(par1, par2, pos) :
    return pos * par1 + par2

# true parameter values
par1True = 42
par2True = 16

# correct forward model
class CorrectDummyModel(ForwardModelBase) :
    def interface(self) :
        self.parameters = ["par1", "par2"]
        self.input_sensors= []
        self.output_sensors= [
            Sensor("S1", std_model="sigma"),
            Sensor("S2", std_model="sigma"),
        ]

    def response(self, inp: dict) :

        res1 = linearModel(inp['par1'], inp['par2'], 10)
        res2 = linearModel(inp['par1'], inp['par2'], -10)

        return {'S1' : res1, 'S2': res2, }

# incorrect forward model
class WrongDummyModel(ForwardModelBase) :
    def interface(self) :
        self.parameters = ["par1Wrong", "par2Wrong"]
        self.input_sensors= []
        self.output_sensors= [
            Sensor("S1", std_model="sigma"),
            Sensor("S2", std_model="sigma"),
        ]

    def response(self, inp: dict) :

        res1 = linearModel(inp['par1Wrong'], inp['par2Wrong'], 10)
        res2 = linearModel(inp['par2Wrong'], inp['par2Wrong'], -10)

        return {'S2': res2, 'S1' : res1, }

# # problem definition

problem = InverseProblem('Linear problem', print_header=False)

# # parameters

# parameters for the correct model
problem.add_parameter(
    "par1",
    prior=Uniform( low = -250.0, high = 250.0),
    tex="$p_1$",
)

problem.add_parameter(
    "par2",
    prior=Uniform( low = -100.0, high = 100.0),
    tex="$p_2$",
)

# parameters for the incorrect model
problem.add_parameter(
    "par1Wrong",
    prior=Uniform( low = -250.0, high = 250.0),
    tex="$p_{1W}$",
)

problem.add_parameter(
    "par2Wrong",
    prior=Uniform( low = -100.0, high = 100.0),
    tex="$p_{2W}$",
)

problem.add_parameter(
    "sigma",
    tex=r"$\sigma$",
    info="Std of zero-mean Gaussian noise model",
    prior=Uniform(low=0.0, high=20.0),
    domain="(0, +oo)"
)

# experimental data 

seed = 1
np.random.seed(seed)

res_test = dict()

res_test['S1'] = linearModel(par1True, par2True, 10)  + np.random.normal(loc=0.0, scale=1)
res_test['S2'] = linearModel(par1True, par2True, -10) + np.random.normal(loc=0.0, scale=1)

problem.add_experiment(
    name='Correct',
    sensor_data= res_test,
)

problem.add_experiment(
    name='Wrong',
    sensor_data= res_test,
)

# add forward models

problem.add_forward_model( CorrectDummyModel("RM"), experiments = ["Correct"] )

problem.add_forward_model( WrongDummyModel("WM"), experiments = ["Wrong"] )

# likelihoods
problem.add_likelihood_model(
    GaussianLikelihoodModel(experiment_name="Correct", model_error="additive")
)
problem.add_likelihood_model(
    GaussianLikelihoodModel(experiment_name="Wrong", model_error="additive")
)

scipy_solver = MaxLikelihoodSolver(problem, show_progress=False)
max_like_data = scipy_solver.run( solver_options = {
                                                    "maxiter" : 5000,
                                                    "maxfev" : 10000,
                                                    })
joergfunger commented 1 year ago

I would opt for sorting (don't think that sorting the dict with the references (the actual vectors are IMO not copied), or even use an ordered dict in python