Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistent sampled values when mixing continuous and discrete distributions #401

Closed
mlchodkowski opened this issue Aug 21, 2023 · 16 comments · Fixed by #402
Closed

Inconsistent sampled values when mixing continuous and discrete distributions #401

mlchodkowski opened this issue Aug 21, 2023 · 16 comments · Fixed by #402

Comments

@mlchodkowski
Copy link

Bug Report

Description:

I encountered a scenario within my use case where I need to sample both floating-point and integer numbers for my model. In the process of defining a params dictionary, I've noticed an inconsistency when combining distributions of type cp.Normal and cp.DiscreteUniform within the vary dictionary. Specifically, when utilizing the SC sampler, integer values are being sampled correctly, but they are incorrectly identified as being of type class <float>. This misidentification triggers constraints related to integer types within the params dictionary.

Furthermore, when employing PCE and MC samplers, I've observed a behavior where the samplers inconsistently sample both float and integer values. This inconsistency persists despite using the cp.DiscreteUniform distribution class (returned values like 18.996 or 24.52345). The reported values, despite being constrained to integers, are still indicated as class <float>. This behavior is inconsistent with the expected behavior of the distribution.

The expectation is that when using cp.DiscreteUniform distribution, integer values should be sampled consistently and correctly identified as integer types, not misclassified as floating-point types.

Steps to Reproduce:

First case

  1. Prepare model provided in Additional context section.
  2. Setup params dictionary:
params = {
    "steamT": {"type": "integer", "default": 600},
    "steamP": {"type": "integer", "default": 30},
    "liquidT": {"type": "float", "default": 10},
    "liquidP": {"type": "float", "default": 20},
    "liquidV": {"type": "integer", "default": 2},
    "tankV": {"type": "integer", "default": 15},
}
  1. Create a sample vary dictionary with example data like this:
vary = {
    "steamT": cp.DiscreteUniform(100, 600),
    "steamP": cp.DiscreteUniform(5, 30),
    "liquidT": cp.Normal(15, 3),
    "liquidP": cp.Normal(10, 2),
    "liquidV": cp.DiscreteUniform(1, 2),
    # "tankV": cp.Normal(15, 2) # Tank V is const
}
  1. Set actions, sampler and execute:
encoder = uq.encoders.GenericEncoder(
    template_fname="config.template", delimiter="$", target_filename="input.json"
)
decoder = uq.decoders.JSONDecoder(
    target_filename="output.json", output_columns=["finalT", "finalP"]
)
execute = ExecuteLocal(f"python3 {os.getcwd()}/steam-water-equilibrium.py input.json")
actions = Actions(CreateRunDirectory("/tmp"), Encode(encoder), execute, Decode(decoder))
campaign = uq.Campaign(name="steam-water", params=params, actions=actions)

### Sampler set
campaign.set_sampler(uq.sampling.SCSampler(vary=vary, polynomial_order=1))
with QCGPJPool(
    template=EasyVVUQParallelTemplate(),
    template_params={
        "venv": r"/home/<username>/.pyenv/versions/3.10.7/envs/steam-water"
    },
) as qcgpj:
    campaign.execute(pool=qcgpj).collate()
  1. Observe error:
>>> RuntimeError: Error when verifying the following new run:
{'steamT': 225.0, 'steamP': 11.0, 'liquidT': 12.0, 'liquidP': 8.0, 'liquidV': 1.0, 'tankV': 15}
Identified errors were:
{'liquidV': ['must be of integer type'], 'steamP': ['must be of integer type'], 'steamT': ['must be of integer type']}

As you can see, the values are correctly sampled as integers, but returned as class<float>. What's interesting is that when setting polynomial_order=1, some parameters, despite being declared as cp.Normal, were also sampled as integers (still returned as class<float>) When setting higher polynomial_order, values with distribution set as cp.Normal or other continuous distribution type are more likely to be floating point numbers like 18.2334.

Second case

  1. Repeat first case up to set_sampler() step.
  2. Set actions, sampler and execute:
### Sampler set
campaign.set_sampler(uq.sampling.PCESampler(vary=vary, polynomial_order=1))
with QCGPJPool(
    template=EasyVVUQParallelTemplate(),
    template_params={
        "venv": r"/home/<username>/.pyenv/versions/3.10.7/envs/steam-water"
    },
) as qcgpj:
    campaign.execute(pool=qcgpj).collate()
  1. Observe error:
>>> RuntimeError: Error when verifying the following new run:
{'steamT': 205.3738767467413, 'steamP': 10.000279551072868, 'liquidT': 12.0, 'liquidP': 8.0, 'liquidV': 0.9999999999999998, 'tankV': 15}
Identified errors were:
{'liquidV': ['must be of integer type'], 'steamP': ['must be of integer type'], 'steamT': ['must be of integer type']}

As you can see here the values were once again returned as class<float> but this time they were not integers? steamT is an integer class, but was sampled as 205.3738767467413, the same goes for steamP variable and liquidV. The only varaibles sampled as integers were ones declared as float.

Environment:

Operating System: [e.g., Ubuntu 22.04 - WSL1]
Python: 3.10.7
EasyVVUQ: currently newest on Pypi, 1.2.
Chaospy version: 4.3.2
Venv manager - pyenv

Additional Context:

Model that I'm using:

#!/usr/bin/env python3

# These code contains material and example problems in thermodynamics.
# This code was used from https://kyleniemeyer.github.io/computational-thermo/
# under Creative Commons Attribution 4.0 International License.

# Steam equilibrating with liquid water
######################################
# Accident scenario: steam leaks into a rigid, insulated tank that is partially filled with water.
# The steam and liquid water are not initially at thermal equilibrium, though they are at the same pressure.
# The steam is at temperature 600 C and pressure 20 MPa.
# The liquid water is initially at 40 C and pressure 20 MPa.
# The total volume of the tank is 10 m^3
# And the volume of the liquid water initially in the tank is 1 m^3.

# Eventually, the contents of the tank reach a uniform temperature and pressure.
# The tank is well-insulated and rigid.

# Problem: Determine the final temperature and pressure of the water in the tank.

import dataclasses
import json
import sys

import numpy as np
import cantera as ct
from pint import UnitRegistry


@dataclasses.dataclass
class Config:
    steamT: float  # Initial steam temperature (degC)
    steamP: float  # Initial steam pressure (MPa)
    liquidT: float  # Initial liquid temperature (degC)
    liquidP: float  # Initial liquid pressure (MPa)
    liquidV: float  # Liquid volume initially in the tank (m^3)
    tankV: float  # Volume of the tank (m^3)

    @staticmethod
    def decode(cfg):
        return Config(**cfg)


if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("ERROR: Config file not provided")
        sys.exit(1)
    with open(sys.argv[1]) as fh:
        config = json.load(fh, object_hook=Config.decode)

    # Setup UnitRegistry to handle unit conversion and manipulation
    ureg = UnitRegistry()
    Q_ = ureg.Quantity

    # Setup Cantera environment providing steam and liquid
    steam = ct.Water()
    temp_steam1 = Q_(config.steamT, "degC")
    pres1 = Q_(config.steamP, "MPa")
    steam.TP = temp_steam1.to("K").magnitude, pres1.to("Pa").magnitude
    steam()

    liquid = ct.Water()
    temp_liquid1 = Q_(config.liquidT, "degC")
    pres_liquid1 = Q_(config.liquidP, "MPa")
    liquid.TP = temp_liquid1.to("K").magnitude, pres_liquid1.to("Pa").magnitude
    liquid()

    # "Calculate the mass of liquid water in the tank, and then determine the volume and mass of steam"
    vol_tank = Q_(config.tankV, "meter^3")
    vol_liquid1 = Q_(config.liquidV, "meter^3")
    mass_liquid1 = vol_liquid1 / Q_(liquid.v, "m^3/kg")
    print(f"Mass of liquid at state 1: {mass_liquid1: .2f}")
    vol_steam1 = vol_tank - vol_liquid1
    mass_steam1 = vol_steam1 / Q_(steam.v, "m^3/kg")
    print(f"Mass of steam at state 1: {mass_steam1: .2f}")
    mass_1 = mass_liquid1 + mass_steam1
    print(f"Total mass of system: {mass_1: .2f}")

    mass_2 = mass_1
    spec_vol2 = vol_tank / mass_2
    print(f"Specific volume of state 2: {spec_vol2: .2e}")
    int_energy2 = (Q_(liquid.u, "J/kg") * mass_liquid1 + Q_(steam.u, "J/kg") * mass_steam1) / mass_2
    int_energy2.ito("kilojoule/kg")
    print(f"Internal energy of state 2: {int_energy2: .2f}")
    water_equilibrium = ct.Water()
    water_equilibrium.UV = int_energy2.to("J/kg").magnitude, spec_vol2.to("m^3/kg").magnitude
    water_equilibrium()

    # "At equilibrium, the tank contains a mixture of saturated liquid and vapor, with temperature and pressure"
    final_temperature = Q_(water_equilibrium.T, "K")
    final_pressure = Q_(water_equilibrium.P, "Pa")

    print(f"Final temperature: {final_temperature: .2f}")
    print(f"Final pressure: {final_pressure.to(ureg.MPa): .2f}")

    with open("output.json", "w+") as fh:
        content = json.dumps({"finalT": final_temperature.magnitude, "finalP": final_pressure.to(ureg.MPa).magnitude})
        fh.write(content)

My template file config.template:

{
  "steamT": $steamT,
  "steamP": $steamP,
  "liquidT": $liquidT,
  "liquidP": $liquidP,
  "liquidV": $liquidV,
  "tankV": $tankV
}

Reproducibility:

  • The bug happens consistently.

Related Issues:

If you found any similar or related issues on the repository, provide links here.

@wedeling
Copy link
Collaborator

Hi @mlchodkowski, thanks for this bug report. I was aware that this can create issues, and some time ago added a quick hack in the SC sampler such that it would actually return integer values. I'll have a look soon

@wedeling
Copy link
Collaborator

@mlchodkowski I'm working in this now, sorry it took so long, was on parental leave due to the birth of my son amongst other things. As a first comment, I think the issue is that you can only have one type in a numpy array. So even if chaospy will give a integer distribution, numpy turns it into a float once combined with the other parameters. A quick workaround that should work for the SC sampler is

campaign = uq.Campaign(name="steam-water", params=params, actions=actions,
                       work_dir=WORK_DIR, verify_all_runs=False)

The verify_all_runs flag just turns off the type check, and since the SC sampler gives integer values this should work. I will try to make a less hacky solution that converts the floats once they are passed from the sample array (sampler.xi_d in the SC case) to the solver. The PCE error you flag is probably due to an incorrectly selected quad rule. Will get back to you later today.

@wedeling
Copy link
Collaborator

@mlchodkowski I think it should be fixed now, without the workaround mentioned above. If you pull the issue401 branch (see #402) you can check. If so, I'll merge it into the dev branch.

@wedeling
Copy link
Collaborator

@mlchodkowski I went ahead and merged the pull request. You can now check your implementation by just pulling the latest dev. If there are still issues let me know.

@mlchodkowski
Copy link
Author

Hi @wedeling,

first of all, thanks for handling this problem. I've successfully pulled the changes into my local repo. AFAIK the SC sampler now works correctly, but I'm still testing the changes on my models.

Right now I've encountered a problem with PCE sampler. AFAIK the invalid quadrature rule that I probably use should now be automatically replaced when using DiscreteUniform distribution. However I've encountered a bug like this:

File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/easyvvuq/sampling/pce.py", line 257, in __init__
    self._nodes, self._weights = cp.generate_quadrature(order=polynomial_order,
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/frontend.py", line 135, in generate_quadrature
    return sparse_grid(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/sparse_grid.py", line 81, in sparse_grid
    x_lookup, w_lookup = _construct_lookup(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/sparse_grid.py", line 146, in _construct_lookup
    (abscissas,), weights = chaospy.generate_quadrature(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/frontend.py", line 152, in generate_quadrature
    abscissas, weights = _generate_quadrature(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/frontend.py", line 248, in _generate_quadrature
    abscissas, weights = quad_function(order, dist, **parameters)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/discrete.py", line 48, in discrete
    order = numpy.where(growth, numpy.where(order > 0, 2**order, 0), order)
ValueError: invalid literal for int() with base 10: 'false'

Additional context (this is my internal JSON format, which is further casted to a proper vary dict):

"sampled_params": [
    {
      "name": "initial_dose",
      "sampling_space": {
        "param_type": "integer",
        "default": 1000
      },
      "distribution": {
        "dist_type": "discrete_uniform",
        "lower": 100,
        "upper": 1000
      }
    },
    {
      "name": "k12",
      "sampling_space": {
        "param_type": "float",
        "default": 0.2
      },
      "distribution": {
        "dist_type": "uniform",
        "lower": 0.3,
        "upper": 0.7
      }
    },
    {
      "name": "k21",
      "sampling_space": {
        "param_type": "float",
        "default": 0.1
      },
      "distribution": {
        "dist_type": "uniform",
        "lower": 0.3,
        "upper": 0.7
      }
    },
    {
      "name": "k10",
      "sampling_space": {
        "param_type": "float",
        "default": 0.15
      },
      "distribution": {
        "dist_type": "uniform",
        "lower": 0.3,
        "upper": 0.7
      }
    },
    {
      "name": "time_steps",
      "sampling_space": {
        "param_type": "integer",
        "default": 100
      }
    },
    {
      "name": "time_interval",
      "sampling_space": {
        "param_type": "integer",
        "default": 1
      }
    }
  ],
  "method": {
    "method_type": "pce",
    "polynomial_order": "3"
  }

I've verified that the problem is not in this format, nor the casting-to-vary step.

@wedeling wedeling reopened this Oct 2, 2023
@wedeling
Copy link
Collaborator

wedeling commented Oct 2, 2023

I've reopened the issue. Are you using a sparse grid or a standard tensor grid?

@mlchodkowski
Copy link
Author

I don't use sparse grid (default arg sparse=False for PCE)

@wedeling
Copy link
Collaborator

wedeling commented Oct 2, 2023

It's a bit hard to see from here, could you post your vary dict and the polynomial order you use?

@mlchodkowski
Copy link
Author

mlchodkowski commented Oct 4, 2023

Sure, here it is:

params = {
    "initial_dose": {"type": "integer", "default": 100},
    "k12": {"type": "float", "default": 0.2},
    "k21": {"type": "float", "default": 0.1},
    "k10": {"type": "float", "default": 0.15},
    "time_steps": {"type": "integer", "default": 100},
    "time_interval": {"type": "integer", "default": 1},
}

vary = {
    "initial_dose": cp.DiscreteUniform(100, 1000),
    "k12": cp.Uniform(0.3, 0.7),
    "k21": cp.Uniform(0.3, 0.7),
    "k10": cp.Uniform(0.3, 0.7)
}

campaign.set_sampler(uq.sampling.PCESampler(vary=vary, polynomial_order=3))

EDIT: And here is the model that I'm using:

#!/usr/bin/env python3

import dataclasses
import json
import random
import sys

import numpy as np
import pandas as pd


@dataclasses.dataclass
class Config:
    # Parameters
    initial_dose: int  # Initial drug dose
    k12: float  # Rate constant from central to peripheral compartment
    k21: float  # Rate constant from peripheral to central compartment
    k10: float  # Elimination rate constant
    time_steps: int  # Number of time steps
    time_interval: int  # Time interval between steps

    @staticmethod
    def decode(cfg):
        return Config(**cfg)


def stochastic_term(amp):
    return amp * random.uniform(0.1, 1)


if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("ERROR: Config file not provided")
        sys.exit(1)
    with open(sys.argv[1]) as fh:
        config = json.load(fh, object_hook=Config.decode)

    amp = 0.15
    # Initialize concentration arrays for central and peripheral compartments
    central_concentration = np.zeros(config.time_steps)
    peripheral_concentration = np.zeros(config.time_steps)

    central_concentration[0] = config.initial_dose

    # Simulation loop
    for t in range(1, config.time_steps):
        # Calculate rate of change for each compartment
        dC1_dt = config.k21 * peripheral_concentration[t - 1] - config.k12 * central_concentration[
            t - 1
        ] * stochastic_term(amp)
        dC2_dt = (
            config.k12 * central_concentration[t - 1]
            - config.k21 * peripheral_concentration[t - 1]
            - config.k10 * peripheral_concentration[t - 1]
        ) * stochastic_term(amp)

        # Update concentrations using Euler's method
        central_concentration[t] = central_concentration[t - 1] + dC1_dt * config.time_interval
        peripheral_concentration[t] = peripheral_concentration[t - 1] + dC2_dt * config.time_interval

    # Plot the results
    time = np.arange(0, config.time_steps * config.time_interval, config.time_interval)

    df = pd.DataFrame({"time": time, "central": central_concentration, "peripheral": peripheral_concentration})
    df.to_csv("output.csv", index=False)

@mlchodkowski
Copy link
Author

@wedeling Hi, I wanted to inform you that I've encountered the same bug also with SC sampler. The vary is the same as above, sampler is SCSampler with polynomial_order=3. I've tested this on you latest commit (0d27a63) to 'dev' branch.

Error:

  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/easyvvuq/sampling/stochastic_collocation.py", line 106, in __init__
    self.check_max_quad_level()
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/easyvvuq/sampling/stochastic_collocation.py", line 247, in check_max_quad_level
    xi_i, wi_i = cp.generate_quadrature(order + j,
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/frontend.py", line 152, in generate_quadrature
    abscissas, weights = _generate_quadrature(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/frontend.py", line 248, in _generate_quadrature
    abscissas, weights = quad_function(order, dist, **parameters)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/quadrature/clenshaw_curtis.py", line 66, in clenshaw_curtis
    order = numpy.where(growth, numpy.where(order > 0, 2**order, 0), order)
ValueError: invalid literal for int() with base 10: 'false'

@wedeling
Copy link
Collaborator

OK, will have a look today.

@wedeling
Copy link
Collaborator

It seems to work on my end, I get results like this (1st-order Sobol indices):

Figure_2

I used the following script, when storing your model as model.py in the same directory

import os
import easyvvuq as uq
import chaospy as cp
import matplotlib.pyplot as plt

from easyvvuq.actions import CreateRunDirectory, Encode, ExecuteLocal, Decode, Actions

WORK_DIR = '/tmp'
OUT_COLS = ["peripheral","central"]

params = {
    "initial_dose": {"type": "integer", "default": 100},
    "k12": {"type": "float", "default": 0.2},
    "k21": {"type": "float", "default": 0.1},
    "k10": {"type": "float", "default": 0.15},
    "time_steps": {"type": "integer", "default": 100},
    "time_interval": {"type": "integer", "default": 1},
}

vary = {
    "initial_dose": cp.DiscreteUniform(100, 1000),
    "k12": cp.Uniform(0.3, 0.7),
    "k21": cp.Uniform(0.3, 0.7),
    "k10": cp.Uniform(0.3, 0.7)
}

encoder = uq.encoders.GenericEncoder(
    template_fname="config.template", delimiter="$", target_filename="input.json")
decoder = uq.decoders.SimpleCSV(target_filename="output.csv", output_columns=OUT_COLS)
execute = ExecuteLocal(f"python3 {os.getcwd()}/model.py input.json")
actions = Actions(CreateRunDirectory(WORK_DIR, flatten=True), Encode(encoder), execute, Decode(decoder))
campaign = uq.Campaign(name="test", params=params, actions=actions, work_dir=WORK_DIR)

# sampler = uq.sampling.SCSampler(vary=vary, polynomial_order=3)
sampler = uq.sampling.PCESampler(vary=vary, polynomial_order=3)

campaign.set_sampler(sampler)

campaign.execute().collate()
data_frame = campaign.get_collation_result()

# analysis = uq.analysis.SCAnalysis(sampler=sampler, qoi_cols=OUT_COLS)
analysis = uq.analysis.PCEAnalysis(sampler=sampler, qoi_cols=OUT_COLS)
campaign.apply_analysis(analysis)

results = campaign.get_last_analysis()

sobols = results.sobols_first()

fig = plt.figure(figsize=[8, 4])
ax1 = fig.add_subplot(121, title=OUT_COLS[0], xlabel="time", ylabel="1st order Sobol index")
ax2 = fig.add_subplot(122, title=OUT_COLS[1], xlabel="time", sharey=ax1)

for param in vary.keys():
    ax1.plot(sobols[OUT_COLS[0]][param], label=param)
    ax2.plot(sobols[OUT_COLS[1]][param], label=param)

plt.legend(loc=0)
plt.tight_layout()
plt.show()

The config.template file is as follows, also stored in the same directory:

{
  "initial_dose": $initial_dose,
  "k12": $k12,
  "k21": $k21,
  "k10": $k10,
  "time_steps": $time_steps,
  "time_interval": $time_interval
}

Could you perhaps run this, and pinpoint at what line it goes wrong?

@mlchodkowski
Copy link
Author

Hi, I've recreated my virtualenv in which i installed the newest easyvvuq. I've run your code from above and I got this error in campaign.apply_analysis(analysis):

/mnt/storage_2/project_data/grant_571/muqsa-app/.venv/lib/python3.10/site-packages/chaospy/descriptives/correlation/pearson.py:46: RuntimeWarning: invalid value encountered in divide
  return numpy.where(vvar > 0, cov/vvar, 0)
Error LinAlgError for peripheral when computing cp.QoI_Dist()

@mlchodkowski
Copy link
Author

mlchodkowski commented Oct 11, 2023

@wedeling I have succesfully localized the bug. The original error was caused by wrong type set in my my pydantic's models for my internal json file. I'm really sorry for the trouble. The error is now fixed. The implementation of SC and PCE samplers allowed for passing the boolean parameters as str's ('false' and 'true' instead of False and True). This was causing the invalid literal for int() with base 10: 'false' on my side.

However the above error (Error LinAlgError) still persist in some cases + there is a new one.

/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/chaospy/descriptives/correlation/pearson.py:46: RuntimeWarning: invalid value encountered in divide
  return numpy.where(vvar > 0, cov/vvar, 0)
/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/muqsa/runner.py:557: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
  result.fillna("", inplace=True)
Traceback (most recent call last):
  File "/opt/exp_soft/local/skylake/python/3.10.7/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/opt/exp_soft/local/skylake/python/3.10.7/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/muqsa/runner.py", line 580, in <module>
    r = Runner(sys.argv[1]).run()
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/muqsa/runner.py", line 64, in run
    self.analyse(campaign, decoder.output_columns, vary, um.method)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/muqsa/runner.py", line 365, in analyse
    Runner.analyse_in_campaign(campaign, qoi_columns, method.__root__.method_type)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/muqsa/runner.py", line 572, in analyse_in_campaign
    results.plot_sobols_treemap(q, filename=f"results/{q}_sobols_treemap_{postfix}.svg")
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/easyvvuq/analysis/results.py", line 464, in plot_sobols_treemap
    squarify.plot(sizes=values, label=keys, color=colors, ax=ax, pad=True)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 244, in plot
    rects = padded_squarify(normed, 0, 0, norm_x, norm_y)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 147, in padded_squarify
    rects = squarify(sizes, x, y, dx, dy)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 136, in squarify
    return layout(current, x, y, dx, dy) + squarify(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 128, in squarify
    while i < len(sizes) and worst_ratio(sizes[:i], x, y, dx, dy) >= worst_ratio(
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 86, in worst_ratio
    for rect in layout(sizes, x, y, dx, dy)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 48, in layout
    layoutrow(sizes, x, y, dx, dy) if dx >= dy else layoutcol(sizes, x, y, dx, dy)
  File "/home/users/mchodkowski/grant_571/project_data/muqsa-app/.venv/lib64/python3.10/site-packages/squarify/__init__.py", line 38, in layoutcol
    height = covered_area / dx
ZeroDivisionError: float division by zero

I suspect the issue is in the results that my model generates some bad or incompatible results. Maybe it's worth to handle this exception?

@wedeling
Copy link
Collaborator

@mlchodkowski no worries, glad you got it working. The LinAlgError seems to be related to the Sobol indices. Perhaps the Sobol indices are computed somewhere where there is zero variance (perhaps a boundary condition or something?). In this case you would get a divide by zero. I could look into catching this exception.

@wedeling wedeling closed this as completed Nov 9, 2023
@wedeling
Copy link
Collaborator

wedeling commented Nov 9, 2023

Closing this

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants