Sandbox

For playing around during development.


In [7]:
import pysd

In [1]:
#%pylab inline
import pysd
print(pysd.__version__)
# #print pysd.__file__
# import pandas as pd
# #import os
# # import xarray as xr
# import numpy as np
# import networkx as nx


0.9.0

In [3]:
l = [1, 2, 3]

while l:
    print(l[0])
    l.pop()


1
1
1

In [25]:
l = set([1, 2, 'burnside'])
j = set(['hi', 'bob'])

In [28]:
l.update(j)
l.difference_update(j)
l


Out[28]:
{1, 2, 'burnside'}

In [15]:
l.remove('burnside')

In [18]:
if l:
    print(l)


{1, 2, 'burnside'}

In [ ]:
x = pd.Series(np.arange(0,40,1))
# y = pd.Series(index=range(0,40))
# y[:10] = x[:10]**2 
# y[10:30] = -1*(x[10:30]-20)**2 + 2*x[10]**2
# y[30:40] = (x[30:40]-40)**2

y = -1*np.cos(x*2*np.pi/40)
plt.plot(x, y)

d1 = y.diff()
d2 = d1.diff()
inc_growth = (d1 > 0) & (d2 > 0)
dec_growth = (d1 > 0) & (d2 < 0)
inc_decline = (d1 < 0) & (d2 < 0)
dec_decline = (d1 < 0) & (d2 > 0)

(inc_growth[inc_growth].index.max() < 
 dec_growth[dec_growth].index.min() <
 dec_growth[dec_growth].index.max() <
 inc_decline[inc_decline].index.min() <
 inc_decline[inc_decline].index.max() <
 dec_decline[dec_decline].index.min())

In [3]:
model = pysd.read_vensim('tests/test-models/samples/teacup/teacup.mdl')
doc = model.doc()
doc


Out[3]:
Real Name Py Name Unit Lims Type Eqn Comment
0 Characteristic Time characteristic_time Minutes : (0.0, None) constant 10 How long will it take the teacup to cool 1/e o...
1 FINAL TIME final_time Minute : (None, None) constant 30 The final time for the simulation.
2 Heat Loss to Room heat_loss_to_room Degrees Fahrenheit/Minute : (None, None) component (Teacup Temperature - Room Temperature) / Char... This is the rate at which heat flows from the ...
3 INITIAL TIME initial_time Minute : (None, None) constant 0 The initial time for the simulation.
4 Room Temperature room_temperature Degrees Fahrenheit : (-459.67, None) constant 70 Put in a check to ensure the room temperature ...
5 SAVEPER saveper Minute : (0.0, None) component TIME STEP The frequency with which output is stored.
6 TIME STEP time_step Minute : (0.0, None) constant 0.125 The time step for the simulation.
7 Teacup Temperature teacup_temperature Degrees Fahrenheit : (32.0, 212.0) component INTEG ( -Heat Loss to Room, 180) The model is only valid for the liquid phase o...

In [3]:
pysd.testing.sample_pspace(model, samples=10)


Out[3]:
Room Temperature Characteristic Time
0 -366.160321 16.488710
1 94.430933 38.938432
2 -176.090179 2.204597
3 1409.317356 3.426968
4 -416.883537 6.072245
5 664.603176 0.747561
6 -237.021635 9.802032
7 247.807415 15.158609
8 -276.607301 3.798534
9 -46.465066 7.966876

In [3]:
doc


Out[3]:
Real Name Py Name Unit Type Comment
0 Characteristic Time characteristic_time Minutes [0,?] constant How long will it take the teacup to cool 1/e o...
4 FINAL TIME final_time Minute constant The final time for the simulation.
1 Heat Loss to Room heat_loss_to_room Degrees Fahrenheit/Minute component This is the rate at which heat flows from the ...
5 INITIAL TIME initial_time Minute constant The initial time for the simulation.
2 Room Temperature room_temperature Degrees Fahrenheit [-459.67,?] constant Put in a check to ensure the room temperature ...
6 SAVEPER saveper Minute [0,?] component The frequency with which output is stored.
7 TIME STEP time_step Minute [0,?] constant The time step for the simulation.
3 Teacup Temperature teacup_temperature Degrees Fahrenheit [32,212] component The model is only valid for the liquid phase o...

In [5]:
set(doc[doc['Type'] == 'constant']['Real Name']) - {'FINAL TIME', 'INITIAL TIME', 'TIME STEP'}


Out[5]:
{'Characteristic Time', 'Room Temperature'}

In [160]:
import pyDOE
from scipy.stats.distributions import uniform, expon, truncnorm, norm
unit_lhd = pyDOE.lhs(2, samples=10, criterion='c')
lhd[:,0] = expon(loc=0, scale=4).ppf(unit_lhd[:,0])
lhd


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-160-e972c553e028> in <module>()
      2 from scipy.stats.distributions import uniform, expon, truncnorm, norm
      3 unit_lhd = pyDOE.lhs(2, samples=10, criterion='c')
----> 4 lhd[:,0] = expon(loc=0, scale=4).ppf(unit_lhd[:,0])
      5 lhd

ValueError: could not broadcast input array from shape (10) into shape (50000)

In [162]:
expon(loc=-1, scale=-1).ppf(unit_lhd[:,0])


Out[162]:
array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan])

In [149]:
lhd = pyDOE.lhs(2, samples=50000, criterion='c')
lhd[:,0] = expon(loc=0, scale=4).ppf(lhd[:,0])
lhd


Out[149]:
array([[ 1.13396347,  0.75099   ],
       [ 2.10951415,  0.51567   ],
       [ 4.38609788,  0.58637   ],
       ..., 
       [ 0.0576535 ,  0.70835   ],
       [ 2.52211701,  0.32877   ],
       [ 1.63385841,  0.44381   ]])

In [151]:
np.isinf(-np.inf)


Out[151]:
True

In [140]:
a = pd.DataFrame(index=range(10))
a['hi'] = 5
a


Out[140]:
hi
0 5
1 5
2 5
3 5
4 5
5 5
6 5
7 5
8 5
9 5

In [132]:
model = pysd.read_vensim('tests/test-models/tests/variable_ranges/test_variable_ranges.mdl')
model.doc()['Real Name'].values


Out[132]:
array(['Unbounded Variable Value Above 0',
       'Unbounded Variable Value Below 0', 'Both Bounds Above 0',
       'Both Bounds Below 0', 'Both Bounds Identical',
       'Unbounded Variable Value 0', 'Lower Bound 0',
       'Lower Bound Above 0', 'Lower Bound Below 0', 'Upper Bound 0',
       'Upper Bound Above 0', 'Upper Bound Below 0', 'Outflow',
       'Lower Bounded Stock', 'Broken Upper Bounded Variable',
       'FINAL TIME', 'INITIAL TIME', 'SAVEPER', 'TIME STEP'], dtype=object)

In [133]:
df = pd.DataFrame(columns=['Unbounded Variable Value Above 0',
       'Unbounded Variable Value Below 0', 'Both Bounds Above 0',
       'Both Bounds Below 0', 'Both Bounds Identical',
       'Unbounded Variable Value 0', 'Lower Bound 0',
       'Lower Bound Above 0', 'Lower Bound Below 0', 'Upper Bound 0',
       'Upper Bound Above 0', 'Upper Bound Below 0'], index=range(10), data=0)
for i, row in df.iterrows():
    a = row
    break
    print(row)

In [110]:
for row in df.itertuples():
    a = row
    break
    print(row)

In [116]:
import scipy.stats

In [118]:
data = scipy.random.normal(loc=10, size=100)
scipy.stats.kstest(rvs=data, cdf='norm')


Out[118]:
KstestResult(statistic=0.99999999999999978, pvalue=0.0)

In [131]:
data = scipy.random.normal(loc=10, size=50000)
a = scipy.stats.kstest(rvs=data, cdf='uniform', args=(10,100))
a.pvalue


Out[131]:
0.0

If there are two ends defined, use a uniform distribution between them If only one end is defined, use exponential distribution with current location set as mean, other end as location If neither end is defined, use normal distribution with current value as mean, magnitude of current value as scale.


In [76]:
def sample_pspace(model=None, bounds=None, param_list=None, samples=100):
    """
    In iterator over various locations in the parameter space,
    useful for

    Parameters
    ----------
    model
    bounds
    param_list
    samples: int
        How many samples to include in the iterator?

    Returns
    -------
    params: dictionary
        dictionary of 'variable':'value' pairs
    """
    if isinstance(bounds, pd.DataFrame):
        bounds = bounds.set_index('Real Name')
    elif bounds is None:
        
    elif isinstance(bounds, str):
        if bounds.split('.')[-1] in ['xls', 'xlsx']:
            bounds = pd.read_excel(bounds, sheetname='Bounds', index_col='Real Name')
        elif bounds.split('.')[-1] == 'csv':
            bounds = pd.read_csv(bounds, index_col='Real Name')
        elif bounds.split('.')[-1] == 'tab':
            bounds = pd.read_csv(bounds, sep='\t', index_col='Real Name')
        else:
            raise ValueError('Unknown file type: bounds')
    else:
        raise ValueError('Unknown type: bounds')

    if param_list is None:
        param_list = list(bounds.index)

    pyDOE.lhs()


4
4
/Users/houghton/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:9: DeprecationWarning: inspect.getargspec() is deprecated, use inspect.signature() or inspect.getfullargspec()

In [29]:
model = pysd.read_vensim('tests/test-models/samples/pendulum/Single_Pendulum.mdl')

In [32]:
res = model.run()
res.plot()


Out[32]:
<matplotlib.axes._subplots.AxesSubplot at 0x11d860d30>

In [34]:
t_series = res['Angular Position']
#def is_periodic(t_series):
fft = np.fft.fft(t_series)
plt.plot(fft)


/Users/houghton/anaconda/lib/python3.6/site-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Out[34]:
[<matplotlib.lines.Line2D at 0x11da67cc0>]

In [38]:
from scipy import signal
xs = np.arange(0, 4*np.pi, 0.05)
data = np.sin(xs)
plt.plot(xs, data)
peakind = signal.find_peaks_cwt(data, np.arange(1,10))
peakind, xs[peakind], data[peakind]


Out[38]:
(array([ 32, 158]), array([ 1.6,  7.9]), array([ 0.9995736 ,  0.99894134]))

In [42]:
peakind = signal.find_peaks_cwt(t_series, np.arange(1,10))
multiple_peaks = len(peakind) > 1
multiple_peaks


Out[42]:
True

In [64]:
segment = int(np.mean(np.diff(peakind))/8)
segment


Out[64]:
8

In [66]:
measure = np.mean([t_series.autocorr(i) for i in range(1,segment)])
measure > .75


Out[66]:
True

In [17]:
a = pd.DataFrame(np.ones((3,3)))
a['hi'] = 3
a
(a.index==2).argmax()


Out[17]:
2

In [22]:
model = pysd.read_vensim('tests/test-models/samples/teacup/teacup.mdl')
model.doc()


Out[22]:
Real Name Py Name Unit Comment
0 Characteristic Time characteristic_time Minutes
1 Heat Loss to Room heat_loss_to_room Degrees/Minute This is the rate at which heat flows from the ...
2 Room Temperature room_temperature
3 Teacup Temperature teacup_temperature Degrees
4 FINAL TIME final_time Minute The final time for the simulation.
5 INITIAL TIME initial_time Minute The initial time for the simulation.
6 SAVEPER saveper Minute [0,?] The frequency with which output is stored.
7 TIME STEP time_step Minute [0,?] The time step for the simulation.

In [26]:
model = pysd.read_vensim('tests/test-models/samples/teacup/teacup.mdl')
result0 = model.run()
result1 = model.run(params={'Room Temperature': 1000})
result2 = model.run()
model.reload()
result3 = model.run()

self.assertTrue((result0 == result3).all().all())
self.assertFalse((result0 == result1).all().all())
self.assertTrue((result1 == result2).all().all())


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-26-18127db6e222> in <module>()
      2 result1 = model.run(params={'Room Temperature': 1000})
      3 result2 = model.run()
----> 4 model.reload()
      5 result3 = model.run()
      6 

AttributeError: 'Model' object has no attribute 'reload'

In [ ]:
model.

In [21]:
model.run()


Out[21]:
Characteristic Time FINAL TIME Heat Loss to Room INITIAL TIME Room Temperature SAVEPER TIME TIME STEP Teacup Temperature Time
0.000 10 30 11.000000 0 70 0.125 0.000 0.125 180.000000 0.000
0.125 10 30 10.862500 0 70 0.125 0.125 0.125 178.625000 0.125
0.250 10 30 10.726719 0 70 0.125 0.250 0.125 177.267188 0.250
0.375 10 30 10.592635 0 70 0.125 0.375 0.125 175.926348 0.375
0.500 10 30 10.460227 0 70 0.125 0.500 0.125 174.602268 0.500
0.625 10 30 10.329474 0 70 0.125 0.625 0.125 173.294740 0.625
0.750 10 30 10.200356 0 70 0.125 0.750 0.125 172.003556 0.750
0.875 10 30 10.072851 0 70 0.125 0.875 0.125 170.728511 0.875
1.000 10 30 9.946940 0 70 0.125 1.000 0.125 169.469405 1.000
1.125 10 30 9.822604 0 70 0.125 1.125 0.125 168.226037 1.125
1.250 10 30 9.699821 0 70 0.125 1.250 0.125 166.998212 1.250
1.375 10 30 9.578573 0 70 0.125 1.375 0.125 165.785734 1.375
1.500 10 30 9.458841 0 70 0.125 1.500 0.125 164.588413 1.500
1.625 10 30 9.340606 0 70 0.125 1.625 0.125 163.406057 1.625
1.750 10 30 9.223848 0 70 0.125 1.750 0.125 162.238482 1.750
1.875 10 30 9.108550 0 70 0.125 1.875 0.125 161.085501 1.875
2.000 10 30 8.994693 0 70 0.125 2.000 0.125 159.946932 2.000
2.125 10 30 8.882260 0 70 0.125 2.125 0.125 158.822595 2.125
2.250 10 30 8.771231 0 70 0.125 2.250 0.125 157.712313 2.250
2.375 10 30 8.661591 0 70 0.125 2.375 0.125 156.615909 2.375
2.500 10 30 8.553321 0 70 0.125 2.500 0.125 155.533210 2.500
2.625 10 30 8.446404 0 70 0.125 2.625 0.125 154.464045 2.625
2.750 10 30 8.340824 0 70 0.125 2.750 0.125 153.408244 2.750
2.875 10 30 8.236564 0 70 0.125 2.875 0.125 152.365641 2.875
3.000 10 30 8.133607 0 70 0.125 3.000 0.125 151.336071 3.000
3.125 10 30 8.031937 0 70 0.125 3.125 0.125 150.319370 3.125
3.250 10 30 7.931538 0 70 0.125 3.250 0.125 149.315378 3.250
3.375 10 30 7.832394 0 70 0.125 3.375 0.125 148.323936 3.375
3.500 10 30 7.734489 0 70 0.125 3.500 0.125 147.344886 3.500
3.625 10 30 7.637808 0 70 0.125 3.625 0.125 146.378075 3.625
... ... ... ... ... ... ... ... ... ... ...
26.375 10 30 0.773966 0 70 0.125 26.375 0.125 77.739657 26.375
26.500 10 30 0.764291 0 70 0.125 26.500 0.125 77.642911 26.500
26.625 10 30 0.754737 0 70 0.125 26.625 0.125 77.547375 26.625
26.750 10 30 0.745303 0 70 0.125 26.750 0.125 77.453032 26.750
26.875 10 30 0.735987 0 70 0.125 26.875 0.125 77.359869 26.875
27.000 10 30 0.726787 0 70 0.125 27.000 0.125 77.267871 27.000
27.125 10 30 0.717702 0 70 0.125 27.125 0.125 77.177023 27.125
27.250 10 30 0.708731 0 70 0.125 27.250 0.125 77.087310 27.250
27.375 10 30 0.699872 0 70 0.125 27.375 0.125 76.998719 27.375
27.500 10 30 0.691123 0 70 0.125 27.500 0.125 76.911235 27.500
27.625 10 30 0.682484 0 70 0.125 27.625 0.125 76.824844 27.625
27.750 10 30 0.673953 0 70 0.125 27.750 0.125 76.739534 27.750
27.875 10 30 0.665529 0 70 0.125 27.875 0.125 76.655289 27.875
28.000 10 30 0.657210 0 70 0.125 28.000 0.125 76.572098 28.000
28.125 10 30 0.648995 0 70 0.125 28.125 0.125 76.489947 28.125
28.250 10 30 0.640882 0 70 0.125 28.250 0.125 76.408823 28.250
28.375 10 30 0.632871 0 70 0.125 28.375 0.125 76.328712 28.375
28.500 10 30 0.624960 0 70 0.125 28.500 0.125 76.249604 28.500
28.625 10 30 0.617148 0 70 0.125 28.625 0.125 76.171483 28.625
28.750 10 30 0.609434 0 70 0.125 28.750 0.125 76.094340 28.750
28.875 10 30 0.601816 0 70 0.125 28.875 0.125 76.018161 28.875
29.000 10 30 0.594293 0 70 0.125 29.000 0.125 75.942934 29.000
29.125 10 30 0.586865 0 70 0.125 29.125 0.125 75.868647 29.125
29.250 10 30 0.579529 0 70 0.125 29.250 0.125 75.795289 29.250
29.375 10 30 0.572285 0 70 0.125 29.375 0.125 75.722848 29.375
29.500 10 30 0.565131 0 70 0.125 29.500 0.125 75.651312 29.500
29.625 10 30 0.558067 0 70 0.125 29.625 0.125 75.580671 29.625
29.750 10 30 0.551091 0 70 0.125 29.750 0.125 75.510912 29.750
29.875 10 30 0.544203 0 70 0.125 29.875 0.125 75.442026 29.875
30.000 10 30 0.537400 0 70 0.125 30.000 0.125 75.374001 30.000

241 rows × 10 columns


In [ ]:


In [6]:
d = {'a':1, 'b':2}

In [10]:
for key, val in d.items():
    print(key, val)


a 1
b 2

In [11]:
import time
time.time()


Out[11]:
1499715812.290142

In [15]:
class MyObj(object):
    def __init__(self):
        self.inittime = time.time()
        
    def reload_self(self):
        self.__init__()

In [16]:
a = MyObj()
a.inittime


Out[16]:
1499715987.2452679

In [17]:
a.reload_self()
a.inittime


Out[17]:
1499715988.093183

In [ ]: