ZeroDivisionError: division by zero

[ali-akhavan89]

Hello,

I have a simple Random Walk model that doesn't have any division formulation inside Vensim. But it seems that pysd does so during truncation calculation. If the standard deviation of the Normal Random in Vensim is zero, pysd throws the following error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/cache.py:25, in Cache.__call__.<locals>.cached_func(*args)
     24 try:
---> 25     return self.data[func.__name__]
     26 except KeyError:

KeyError: 'shock'

During handling of the above exception, another exception occurred:

ZeroDivisionError                         Traceback (most recent call last)
Cell In[1], line 4
      2 model = pysd.read_vensim('RandomWalk.mdl')
      3 model = pysd.load('RandomWalk.py')
----> 4 stocks = model.run()
      5 stocks

File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/model.py:1551, in Model.run(self, params, return_columns, return_timestamps, initial_condition, final_time, time_step, saveper, reload, progress, flatten_output, cache_output, output_file)
   1548 self.output.set_capture_elements(self.capture_elements)
   1549 self.output.initialize(self)
-> 1551 self._integrate()
   1553 return self.output.collect(self, flatten_output)

File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/model.py:2339, in Model._integrate(self)
   2337 while self.time.in_bounds():
   2338     if self.time.in_return():
-> 2339         self.output.update(self)
   2341     self._integrate_step()
   2342     progressbar.update()

File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/output.py:483, in ModelOutput.update(self, model)
    478 def update(self, model):
    479     """
    480     Delegating the update of the results object and its elements
    481     to the appropriate handler.
    482     """
--> 483     self.handler.update(model)

File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/output.py:298, in DataFrameHandler.update(self, model)
    296 self.ds['time'].append(model.time.round())
    297 for key in self.capture_elements_step:
--> 298     self.ds[key].append(getattr(model.components, key)())
    300 self.__step += 1

File /opt/miniconda3/lib/python3.12/site-packages/pysd/py_backend/cache.py:27, in Cache.__call__.<locals>.cached_func(*args)
     25     return self.data[func.__name__]
     26 except KeyError:
---> 27     value = func(*args)
     28     self.data[func.__name__] = value
     29     return value

File ~/pysd/RandomWalk.py:134, in shock()
    126 @component.add(
    127     name="Shock",
    128     comp_type="Auxiliary",
   (...)
    131 )
    132 def shock():
    133     return stats.truncnorm.rvs(
--> 134         (-10000 - 0) / noisestd(),
    135         (10000 - 0) / noisestd(),
    136         loc=0,
    137         scale=noisestd(),
    138         size=(),
    139     )

ZeroDivisionError: division by zero

Here's the script I use:

import pysd
model = pysd.read_vensim('RandomWalk.mdl')
model = pysd.load('RandomWalk.py')
stocks = model.run()
stocks

I've also attached the .mdl and .py files.

RandomWalk.zip

Let me know if I can provide further information.

Best,
Ali

[ali-akhavan89]

I'm also a bit worried about the division by standard deviation in the code, and I'd appreciate if you could explain the reason behind that.

Thanks again.
Ali

[enekomartinmartinez]

Hi, thank you for opening this.

I don't understand the meaning of having a standard deviation equal to 0 with a random function.

What is Vensim doing in this case? Is it returning the h shift parameter value always? Is it ignoring the given std and giving a distribution with std equal to 1? This is not specified in their documentation, so if you can test that let me know please, I have no Vensim available anymore.

Vensim a lot of times allows operations that make no sense, e.g. allocation functions work with negative values and give an unexpected output without raising any error or warning. So the fact that works in Vensim doesn't mean that it is working well.

We are dividing by the standard deviation to have the same output as Vensim. Scipy's documentation gives more information about the distribution https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.truncnorm.html

In their note they have the equation

a, b = (a_trunc - loc) / scale, (b_trunc - loc) / scale

while Vensim takes a_trunc (m) and b_trunc (x) as parameters, scipy functions need a and b, so we need to remove loc (h) and divide by scale (r). But I would personally expect scale being > 0 always.

This transformation is tested in our test suite as we compared different distribution outputs by Vensim with the results of our translations.

[ali-akhavan89]

Thank you!

Vensim returns loc=mean if scale=std_dev=0, same as NumPy and SciPy.

I'm not an expert, but I would use the regular random generator function and check if the returned values are within bounds (min, max). If not, I'd replace them with another draw from the function.

Hope it helps,
Ali

[enekomartinmartinez]

Includying truncation after the generation of the values may be tricky when the variable has subscripts. I prefer using zidz functions when computing a and b for scipy, it will also keep the code clean. Could you confirm that the same happens with RANDOM EXPONENTIAL when r=0? This way I will make the fix for both distributions.
Thank you!

[ali-akhavan89]

Thank you. Yes, RANDOM EXPONENTIAL returns the value for shift when stretch = 0

[enekomartinmartinez]

Thank you, It will take me some days to do the fix, as I am currently traveling. I will try to make it the next week. Eneko 2025(e)ko urtarrilakren 27(a) 15:23:16 (CET)-(e)an, ali-akhavan89 ***@***.***>-(e)k hau idatzi zuen:

…

Thank you! Vensim returns `loc=mean` if `scale=std_dev=0`, same as NumPy and SciPy. I'm not an expert, but I would use the regular random generator function and check if the returned values are within bounds (min, max). If not, I'd replace them with another draw from the function. Hope it helps, Ali -- Reply to this email directly or view it on GitHub: #461 (comment) You are receiving this because you commented. Message ID: ***@***.***>