"""Flips a coin with the given probability.

"""Convert Cartesian coordinates to polar.

x = np.asarray(x)

y = np.asarray(y)

# TODO: use hypot?

rho = np.sqrt(x**2 + y**2)

theta = np.arctan2(y, x)

if z is None:

return theta, rho

else:

"""Convert polar coordinates to Cartesian.

if hasattr(theta, 'units'):

if theta.units == UNITS.degree:

theta = theta.to(UNITS.radian)

if theta.units != UNITS.radian:

msg = """In pol2cart, theta must be either a number or

raise ValueError(msg)

x = rho * np.cos(theta)

y = rho * np.sin(theta)

if z is None:

return x, y

else:

"""Returns an array of evenly-spaced values in the interval [start, stop].

underride(options, dtype=np.float64)

# see if either of the arguments has units

units = getattr(start, 'units', None)

units = getattr(stop, 'units', units)

array = np.linspace(start, stop, num, **options)

if units:

array = array * units

"""Returns an array of evenly-spaced values in an interval.

if stop is None:

stop = start

start = 0

underride(options, endpoint=False)

# see if any of the arguments has units

units = getattr(start, 'units', None)

units = getattr(stop, 'units', units)

units = getattr(step, 'units', units)

n = np.round((stop - start) / step)

if options['endpoint']:

n += 1

array = np.full(int(n), magnitude(step))

array[0] = magnitude(start)

array = np.cumsum(array)

if units:

array = array * units

"""Returns the magnitude of a Quantity or number.

"""Returns the magnitude of a Quantity or number, or sequence.

try:

return [magnitude(elt) for elt in x]

except TypeError: # not iterable

"""Returns the units of a Quantity or number.

"""Removes units from the values in a Series.

res = copy(series)

print(type(res))

for label, value in res.iteritems():

res[label] = magnitude(value)

"""Apply units to `x`, if necessary.

if isinstance(x, Quantity):

return x.to(units)

else:

"""Find the parameters that yield the best fit for the data.

# if any of the params are quantities, strip the units

x0 = [magnitude(x) for x in params]

# override `full_output` so we get a message if something goes wrong

options['full_output'] = True

# run leastsq

with units_off():

best_params, cov_x, infodict, mesg, ier = leastsq(error_func,

x0=x0, args=args, **options)

details = ModSimSeries(infodict)

details.set(cov_x=cov_x, mesg=mesg, ier=ier)

# if we got a Params object, we should return a Params object

if isinstance(params, Params):

best_params = Params(Series(best_params, params.index))

# return the best parameters and details

"""Finds the input value that minimizes `min_func`.

# try:

# print(bounds[0])

# min_func(bounds[0], *args)

# except Exception as e:

# msg = """Before running scipy.integrate.min_bounded, I tried

# running the slope function you provided with the

# initial conditions in system and t=0, and I got

# the following error:"""

# logger.error(msg)

# raise(e)

underride(options, xatol=1e-3)

# TODO: Do we need to remove units from bounds?

with units_off():

res = minimize_scalar(min_func,

bracket=bounds,

bounds=bounds,

args=args,

method='bounded',

options=options)

if not res.success:

msg = """scipy.optimize.minimize_scalar did not succeed.

The message it returned is %s""" % res.message

raise Exception(msg)

"""Finds the input value that maximizes `max_func`.

def min_func(*args):

return -max_func(*args)

res = min_bounded(min_func, bounds, *args, **options)

# we have to negate the function value before returning res

res.fun = -res.fun

"""Finds the input value that minimizes `min_func`.

underride(options, tol=1e-3)

with units_off():

res = scipy.optimize.minimize(min_func, x0, *args, **options)

"""Integrates an ordinary differential equation.

# make sure `system` contains `ts`

if not hasattr(system, 'ts'):

msg = """It looks like `system` does not contain `ts`

raise ValueError(msg)

# make sure `system` contains `init`

if not hasattr(system, 'init'):

msg = """It looks like `system` does not contain `init`

raise ValueError(msg)

# make the system parameters available as globals

unpack(system)

# try running the slope function with the initial conditions

try:

slope_func(init, ts[0], system)

except Exception as e:

msg = """Before running scipy.integrate.odeint, I tried

logger.error(msg)

raise(e)

# when odeint calls slope_func, it should pass `system` as

# the third argument. To make that work, we have to make a

# tuple with a single element and pass the tuple to odeint as `args`

args = (system,)

# now we're ready to run `odeint` with `init` and `ts` from `system`

with units_off():

array = odeint(slope_func, list(init), ts, args, **options)

# the return value from odeint is an array, so let's pack it into

# a TimeFrame with appropriate columns and index

frame = TimeFrame(array, columns=init.index, index=ts, dtype=np.float64)

"""Computes a numerical solution to a differential equation.

# make sure `system` contains `init`

if not hasattr(system, 'init'):

msg = """It looks like `system` does not contain `init`

raise ValueError(msg)

# make sure `system` contains `t_end`

if not hasattr(system, 't_end'):

msg = """It looks like `system` does not contain `t_end`

raise ValueError(msg)

# make the system parameters available as globals

unpack(system)

# the default value for t_0 is 0

t_0 = getattr(system, 't_0', 0)

# try running the slope function with the initial conditions

# try:

# slope_func(init, t_0, system)

# except Exception as e:

# msg = """Before running scipy.integrate.solve_ivp, I tried

# running the slope function you provided with the

# initial conditions in `system` and `t=t_0` and I got

# the following error:"""

# logger.error(msg)

# raise(e)

# wrap the slope function to reverse the arguments and add `system`

f = lambda t, y: slope_func(y, t, system)

def wrap_event(event):

"""Wrap the event functions.

wrapped = lambda t, y: event(y, t, system)

wrapped.terminal = getattr(event, 'terminal', True)

wrapped.direction = getattr(event, 'direction', 0)

return wrapped

# wrap the event functions so they take the right arguments

events = options.pop('events', [])

try:

events = [wrap_event(event) for event in events]

except TypeError:

events = wrap_event(events)

# remove dimensions from the initial conditions.

# we need this because otherwise `init` gets copied into the

# results array along with its units

# try:

# y_0 = [magnitude(x) for x in init]

# except TypeError:

# y_0 = [magnitude(init)]

y_0 = [magnitude(x) for x in init]

# run the solver

with units_off():

bunch = solve_ivp(f, [t_0, t_end], y_0, events=events, **options)

# separate the results from the details

y = bunch.pop('y')

t = bunch.pop('t')

details = ModSimSeries(bunch)

# pack the results into a TimeFrame

results = TimeFrame(np.transpose(y), index=t, columns=init.index)

"""Return the roots of the (non-linear) equations

# make sure we can run the given function with x0

try:

func(x0, *args)

except Exception as e:

msg = """Before running scipy.optimize.fsolve, I tried

logger.error(msg)

raise(e)

# make the tolerance more forgiving than the default

underride(options, xtol=1e-6)

x0 = magnitude(x0)

# run fsolve

with units_off():

result = scipy.optimize.fsolve(func, x0, args=args, **options)

"""Find the labels where the series passes through value.

interp = InterpolatedUnivariateSpline(series.index, series-value)

"""Creates an interpolation function.

# TODO: add error checking for nonmonotonicity

if sum(series.index.isnull()):

msg = """The Series you passed to interpolate contains

raise ValueError(msg)

# make the interpolate function extrapolate past the ends of

# the range, unless `options` already specifies a value for `fill_value`

underride(options, fill_value='extrapolate')

# call interp1d, which returns a new function object

interp_func = interp1d(series.index, series.values, **options)

units = getattr(series, 'units', None)

if units:

return lambda x: Quantity(interp_func(x), units)

else:

"""Interpolate the inverse function of a Series.

inverse = Series(series.index, index=series.values)

T = interpolate(inverse, **options)

"""Make the names in `series` available as globals.

# TODO: Make this a context manager, so the syntax is

# with series:

# and maybe even add an __exit__ that copies changes back

frame = inspect.currentframe()

caller = frame.f_back

"""Prints the source code for a given object.

"""Add key-value pairs to d only if key is not in d.

if d is None:

d = {}

for key, val in options.items():

d.setdefault(key, val)

"""Makes line plots.

# TODO: add lines to REPLOT_CACHE

x, y, style = parse_plot_args(*args, **options)

if isinstance(x, pd.DataFrame) or isinstance(y, pd.DataFrame):

raise ValueError("modsimpy.plot can't handle DataFrames.")

if x is None:

if isinstance(y, np.ndarray):

x = np.arange(len(y))

if isinstance(y, pd.Series):

x = y.index

y = y.values

x = magnitudes(x)

y = magnitudes(y)

underride(options, linewidth=3, alpha=0.6)

if style is not None:

lines = plt.plot(x, y, style, **options)

else:

lines = plt.plot(x, y, **options)

"""

try:

label = options['label']

except KeyError:

raise ValueError('To use replot, you must provide a label argument.')

axes = plt.gca()

key = (axes, label)

if key not in REPLOT_CACHE:

lines = plot(*args, **options)

if len(lines) != 1:

raise ValueError('Replot only works with a single plotted element.')

REPLOT_CACHE[key] = lines[0]

return lines

line = REPLOT_CACHE[key]

x, y, style = parse_plot_args(*args, **options)

line.set_xdata(x)

"""Parse the args the same way plt.plot does."""

x = None

y = None

style = None

if len(args) == 1:

y = args[0]

elif len(args) == 2:

if isinstance(args[1], str):

y, style = args

else:

x, y = args

elif len(args) == 3:

x, y, style = args

"""Makes a contour plot from a DataFrame.

x = df.columns

y = df.index

X, Y = np.meshgrid(x, y)

cs = plt.contour(X, Y, df, **options)

"""Save the current figure.

print('Saving figure to file', filename)

"""Decorate the current axes.

loc = options.pop('loc', 'best')

if options.pop('legend', True):

legend(loc=loc)

plt.gca().set(**options)

"""Draws a legend only if there is at least one labeled item.

underride(options, loc='best')

ax = plt.gca()

handles, labels = ax.get_legend_handles_labels()

"""Removes some labels from the legend.

ax = plt.gca()

handles, labels = ax.get_legend_handles_labels()

handle_list, label_list = [], []

for handle, label in zip(handles, labels):

if label not in bad_labels:

handle_list.append(handle)

label_list.append(label)

def __init__(self, fig, axes_seq):

self.fig = fig

self.axes_seq = axes_seq

self.current_axes_index = 0

def current_axes(self):

return self.axes_seq(self.current_axes_index)

# TODO: consider making SubPlots iterable

def next_axes(self):

self.current_axes_index += 1

figsize = {(2, 1): (8, 8),

(3, 1): (8, 10)}

key = nrows, ncols

default = (8, 5.5)

width, height = figsize.get(key, default)

plt.subplot(nrows, ncols, plot_number, **options)

fig = plt.gcf()

fig.set_figwidth(width)

"""Modified version of a Pandas Series,

def __init__(self, *args, **kwargs):

"""Initialize a Series.

if args or kwargs:

underride(kwargs, copy=True)

super().__init__(*args, **kwargs)

else:

super().__init__([], dtype=np.float64)

def _repr_html_(self):

"""Returns an HTML representation of the series.

df = pd.DataFrame(self.values, index=self.index, columns=['values'])

return df._repr_html_()

def __copy__(self, deep=True):

series = super().copy(deep=deep)

return self.__class__(series)

copy = __copy__

def set(self, **kwargs):

"""Uses keyword arguments to update the Series in place.

for name, value in kwargs.items():

self[name] = value

@property

def dt(self):

"""Intercept the Series accessor object so we can use `dt`

return self.loc['dt']

@property

def T(self):

"""Intercept the Series accessor object so we can use `T`

"""Contains system variables and their values.

def __init__(self, *args, **kwargs):

"""Initialize the series.

if len(args) == 0:

super().__init__(list(kwargs.values()), index=kwargs)

elif len(args) == 1:

super().__init__(*args, copy=True)

self.set(**kwargs)

else:

msg = '__init__() takes at most one positional argument'

"""Contains state variables and their values.

"""Represents the condition of a system.

xs = np.asarray(seq)

diff = np.ediff1d(xs, np.nan)

if isinstance(seq, Series):

return Series(diff, seq.index)