text:
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

import Utilities.KPDF as KPDF
import statsmodels.nonparametric.kernel_density as smkde
import statsmodels.nonparametric.bandwidths as smbw

import seaborn as sns
sns.set_context("paper")

xrvs = [stats.norm(loc=1.0, scale=0.2), stats.norm(loc=3.0, scale=0.1)]
yrvs = [stats.norm(loc=0.0, scale=0.2), stats.norm(loc=1.0, scale=0.3)]

N = 3000
M = N / (len(xrvs) + len(yrvs))

rvs = []
for xrv in xrvs:
    for yrv in yrvs:
        rvs.append(np.hstack([
            xrv.rvs(size=(M, 1)),
            yrv.rvs(size=(M, 1))]))
rvs = np.vstack(rvs)

def pdf(x, y):
    return (sum([rv.pdf(x) for rv in xrvs]) *
            sum([rv.pdf(y) for rv in yrvs]))

def plot(ax, Z, title):
    ax.scatter(rvs[:, 0], rvs[:, 1], alpha=0.3, s=1, color='black')
    ax.imshow(Z, aspect=1, origin='lower',
              cmap=sns.light_palette((210, 90, 60),
                                     input='husl',
                                     as_cmap=True),
              extent=(rvs[:, 0].min(),
                      rvs[:, 0].max(),
                      rvs[:, 1].min(),
                      rvs[:, 1].max()))
    ax.contour(x, y, Z, linewidth=1)
    ax.set_title(title)

x_flat = np.r_[rvs[:, 0].min():rvs[:, 0].max():128j] #pylint: disable=E1127
y_flat = np.r_[rvs[:, 1].min():rvs[:, 1].max():128j] #pylint: disable=E1127
x, y = np.meshgrid(x_flat, y_flat)
grid = np.hstack([x.reshape(-1, 1), y.reshape(-1, 1)])

n, d = rvs.shape
bw_kpdf = KPDF.MPDFOptimumBandwidth(rvs)
bw_scott = n ** (-1.0 / (d + 4))
bw_silverman = (n * (d + 2) / 4.) ** (-1. / (d + 4))

scale = np.array([[x_flat.ptp(), y_flat.ptp()]])
print('scale: %s' % scale)

print('bandwidth kpdf=%f scott=%f silverman=%f' %
      (bw_kpdf, bw_scott, bw_silverman))

fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(8, 11))
axes = axes.flat

p = pdf(x, y)
p = p / p.max()
plot(axes[0], p, 'True density')

kde = stats.kde.gaussian_kde(rvs.T)
print(kde.covariance)
z = kde(grid.T)
z = z.reshape(x.shape) / z.max()
plot(axes[1], z, 'Scipy ($\ell_2$ norm: %.3f)' % np.linalg.norm((p - z).flat))

bw = bw_kpdf
# bw = bw_scott

w = KPDF.MPDFGaussian(rvs, grid, bw_kpdf)
w = w.reshape(x.shape) / w.max()
plot(axes[2], w, 'KPDF bw:kpdf ($\ell_2$ norm: %.3f)' % np.linalg.norm((p -
     w).flat))

w = KPDF.MPDFGaussian(rvs, grid, bw_kpdf / 2)
w = w.reshape(x.shape) / w.max()
plot(axes[3], w, 'KPDF bw:kpdf/2 ($\ell_2$ norm: %.3f)' % np.linalg.norm((
    p - w).flat))

w = KPDF.MPDFGaussian(rvs, grid, bw_scott)
w = w.reshape(x.shape) / w.max()
plot(axes[4], w, 'KPDF bw:scott ($\ell_2$ norm: %.3f)' % np.linalg.norm((p -
     w).flat))

w = KPDF.MPDFGaussian(rvs, grid, bw_scott / 2)
w = w.reshape(x.shape) / w.max()
plot(axes[5], w, 'KPDF bw:scott/2 ($\ell_2$ norm: %.3f)' % np.linalg.norm((
    p - w).flat))

dens = smkde.KDEMultivariate(rvs, 'cc', bw='cv_ml')
print "SM bandwidth (cv_ml): " + repr(dens.bw)
w = dens.pdf(grid)
w = w.reshape(x.shape) / w.max()
plot(axes[6], w, 'SM bw:CVML ($\ell_2$ norm: %.3f)' % np.linalg.norm((p -
     w).flat))

dens = smkde.KDEMultivariate(rvs, 'cc', bw='cv_ls')
print "SM bandwidth (cv_ls): " + repr(dens.bw)
w = dens.pdf(grid)
w = w.reshape(x.shape) / w.max()
plot(axes[7], w, 'SM bw:CVLS ($\ell_2$ norm: %.3f)' % np.linalg.norm((p -
     w).flat))

plt.savefig('fig3.png')
plt.show()

max_stars_repo_name: GlobalParametrics/TCRM

max_issues_repo_name: GlobalParametrics/TCRM

max_forks_repo_name: GlobalParametrics/TCRM

reason: import numpy,import scipy,import statsmodels


text:
from astropy.convolution import Box1DKernel
import emcee
import numpy as np
from scipy import integrate, optimize, signal

from .utils import acf, get_noise, smooth

def Vrng_Basri2011(y):
    """Basri et al. 2011, AJ, 141, 20"""
    vrng = np.percentile(y, 95) - np.percentile(y, 5)
    return vrng

def HFrms_Basri2011(y):
    """Basri et al. 2011, AJ, 141, 20"""
    yf = smooth(y, Box1DKernel(width=4))
    resid = y - yf
    hfrms = np.std(resid)
    return hfrms

def Sph_Mathur2014(y, dy=None):
    """Mathur et al. 2014, A&A, 562, 124"""
    if dy is None:
        noise = get_noise(y)
    else:
        noise = np.min(dy)
    sph = np.sqrt(np.std(y) ** 2 - noise ** 2)
    return sph

def Sph_timeseries_Mathur2014(y, t, prot, dy=None, k=5):
    """Mathur et al. 2014, A&A, 562, 124"""
    window = k * prot
    mint = np.min(t)
    maxt = np.max(t)
    t_sph, y_sph = [], []
    N = int(max(1, (maxt - mint) // window))
    for i in range(N):
        indices = (t >= mint + i * window) & (t < mint + (i + 1) * window)
        yi = y[indices]
        if dy is not None:
            dyi = dy[indices]
        else:
            dyi = None
        try:
            y_sph = np.append(y_sph, Sph_Mathur2014(yi, dyi))
            t_sph = np.append(t_sph, mint + (i + 0.5) * window)
        except ValueError:
            continue
    t_sph = t_sph[np.isfinite(y_sph)]
    y_sph = y_sph[np.isfinite(y_sph)]
    return t_sph, y_sph

def Contrast_Mathur2014(y, t, prot, dy=None, k=5):
    """Mathur et al. 2014, A&A, 562, 124"""
    to, yo = Sph_timeseries_Mathur2014(y, t, prot, dy, k)
    sph = Sph_Mathur2014(y, dy)
    nlo = yo < sph
    nhi = ~nlo
    C = yo[nhi].mean() / yo[nlo].mean()
    return C

def m_k(y, t, prot, dy=None, k=5):
    to, yo = Sph_timeseries_Mathur2014(y, t, prot, dy, k)
    return yo.std() / yo.mean()

def Reff_He2015(y):
    """He et al. 2015, ApJS, 221, 18"""
    reff = np.sqrt(8) * np.std(y)
    return reff

def iAC_He2015(y):
    """He et al. 2015, ApJS, 221, 18"""
    N = y.size // 2
    lags, ryy = acf(y, maxlag=N)
    iac = integrate.simps(np.abs(ryy), lags) / N
    return iac

def SDR_Basri2018(t, y, prot):
    """Basri & Nguyen 2018, ApJ, 863, 190"""
    if prot > 8:
        yf = smooth(y, Box1DKernel(width=prot / 8))
    else:
        yf = y.copy()
    dips, _ = signal.find_peaks(-yf, width=4)
    seps = np.diff(t[dips])
    singledips = seps > 0.75 * prot
    doubledips = ~singledips
    singlemode = seps[singledips].sum()
    doublemode = seps[doubledips].sum()
    sdr = np.log10(singlemode / doublemode)
    return sdr

def tau_ACF_Santos2021(lags, ryy, prot, random_seed=None, use_prior=False):
    """Santos et al. 2021"""
    rng = np.random.default_rng(random_seed)
    maxima, _ = signal.find_peaks(ryy)
    idx = np.where(np.diff(np.diff(np.sign(np.diff(ryy[maxima])))) == -2)[0][0] + 2
    max_lag = maxima[idx]
    nflux_j = lags.size - np.arange(max_lag)
    ryy = ryy[:max_lag]
    lags = lags[:max_lag]
    sigma_p = 0.2 * prot
    sigma_j = 1 / np.sqrt(nflux_j)

    def predict(theta):
        tau, a, b, y0, P = theta
        phi = 2 * np.pi * lags / P
        return (1 - lags / tau) * (a * np.cos(phi) + b * np.cos(2 * phi) + y0)

    def log_prior(theta):
        tau, a, b, y0, P = theta
        log_tau = np.log10(tau)
        if not (
            (-1 < log_tau < 3)
            and (0.01 < a < 1.5)
            and (0.001 < b < 1)
            and (-1 < y0 < 1)
        ):
            return -np.inf
        chi = ((P - prot) / sigma_p) ** 2
        norm_const = -0.5 * np.log(2 * np.pi) - 0.5 * np.log(sigma_p)
        lp = norm_const - 0.5 * chi
        return lp

    def nll(theta):
        y = predict(theta)
        chi = np.sum(np.square((ryy - y) / sigma_j))
        norm_const = -0.5 * max_lag * np.log(2 * np.pi) - 0.5 * np.log(sigma_j).sum()
        return 0.5 * chi - norm_const

    def logprob(theta):
        return log_prior(theta) - nll(theta)

    n_walkers = 25
    n_dim = 5
    n_steps = 10_000
    burn = 8_000
    if use_prior:
        p0 = np.vstack(
            [
                10 ** rng.uniform(-1, 3, n_walkers),
                rng.uniform(0.01, 1.5, n_walkers),
                rng.uniform(0.001, 1, n_walkers),
                rng.uniform(-1, 1, n_walkers),
                rng.normal(prot, sigma_p, n_walkers),
            ]
        ).T
    else:
        x0 = [5 * prot, 0.75, 0.5, 0.0, prot]
        bounds = [
            (0.1, 1000),
            (0.01, 1.5),
            (0.001, 1),
            (-1, 1),
            (0.2 * prot, 1.8 * prot),
        ]
        res = optimize.minimize(fun=nll, x0=x0, method="L-BFGS-B", bounds=bounds)
        p0 = res.x + 1e-2 * rng.standard_normal((n_walkers, n_dim))
    sampler = emcee.EnsembleSampler(n_walkers, n_dim, logprob)
    sampler.run_mcmc(p0, n_steps, progress=True)
    samples = sampler.get_chain(discard=burn, flat=True)
    tau = np.percentile(samples[:, 0], [16, 50, 84], axis=0)
    theta_MAP = np.median(samples, axis=0)
    reduced_chi = np.sum(np.square((ryy - predict(theta_MAP)) / sigma_j)) / (
        max_lag - n_dim
    )
    return tau, theta_MAP, reduced_chi

ext: py

max_stars_repo_name: dioph/edastra

max_issues_repo_name: dioph/edastra

max_forks_repo_name: dioph/edastra

reason: import numpy,from scipy,from astropy


text:
__author__ = 'Izzat'
import wx
import wx.lib.scrolledpanel as scrolled
from core.symoro import SIMPLE, TREE, CLOSED_LOOP
from sympy import Expr, Symbol
#TODO: PROG_NAME

class DialogConversion(wx.Dialog):
    def __init__(self, prefix, robo, par_dict, parent=None):
        super(DialogConversion, self).__init__(parent)
        self.robo = robo
        self.par_dict = par_dict
        self.SetTitle(prefix + ": Enter numerical values")
        self.construct_sym()
        self.init_ui()

    def has_syms(self):
        return len(self.syms) > 0

    def construct_sym(self):
        params = self.robo.get_geom_head()[4:]
        q_vec = self.robo.q_vec
        self.syms = set()
        for par in params:
            for i in range(1, self.robo.NF):
                val = self.robo.get_val(i, par)
                if val in q_vec:
                    continue
                if isinstance(val, Expr):
                    for at in val.atoms(Symbol):
                        self.syms.add(at)

    def init_ui(self):
        p = scrolled.ScrolledPanel(self, -1)
        vbox = wx.BoxSizer(wx.VERTICAL)
        self.widgets = {}
        for par in self.syms:
            if par in self.par_dict:
                val = str(self.par_dict[par])
            else:
                val = 1.
            hor_sizer = wx.BoxSizer(wx.HORIZONTAL)
            label = wx.StaticText(p, label=str(par), size=(60, -1),
                                  style=wx.ALIGN_RIGHT)
            hor_sizer.Add(label, 0, wx.ALIGN_CENTER_VERTICAL)
            hor_sizer.AddSpacer(5)
            txt_box = wx.TextCtrl(p, size=(120, -1), value=str(val))
            hor_sizer.Add(txt_box)
            vbox.Add(hor_sizer, 0, wx.ALL | wx.ALIGN_CENTER_HORIZONTAL, 5)
            self.widgets[str(par)] = txt_box

        p.SetSizer(vbox)
        p.SetAutoLayout(1)
        p.SetupScrolling()

        main_sizer = wx.BoxSizer(wx.VERTICAL)
        main_sizer.Add(p, 1, wx.ALL | wx.EXPAND, 2)
        hor_sizer = wx.BoxSizer(wx.HORIZONTAL)
        ok_btn = wx.Button(self, wx.ID_OK, "OK")
        ok_btn.Bind(wx.EVT_BUTTON, self.OnOK)
        cancel_btn = wx.Button(self, wx.ID_CANCEL, "Cancel")
        cancel_btn.Bind(wx.EVT_BUTTON, self.OnCancel)
        hor_sizer.Add(ok_btn, 0, wx.ALL, 15)
        hor_sizer.Add(cancel_btn, 0, wx.ALL, 15)
        main_sizer.Add(hor_sizer, 0, wx.ALIGN_CENTER_HORIZONTAL, 0)
        self.SetSizer(main_sizer)
        #if n_syms == 0:
        #   self.EndModal(wx.ID_OK)

    def OnOK(self, _):
        self.EndModal(wx.ID_OK)

    def OnCancel(self, _):
        self.EndModal(wx.ID_CANCEL)

    def get_values(self):
        result = {}
        for par in self.syms:
            result[par] = self.widgets[str(par)].Value
        return result

class DialogDefinition(wx.Dialog):
    def __init__(self, prefix, name, nl, nj, structure, is_mobile, parent=None):
        super(DialogDefinition, self).__init__(parent, style=wx.SYSTEM_MENU |
                                wx.CAPTION | wx.CLOSE_BOX | wx.CLIP_CHILDREN)
        self.init_ui(name, nl, nj, is_mobile, structure)
        self.SetTitle(prefix + ": New robot definition")

    def init_ui(self, name, nl, nj, is_mobile, structure):
        main_sizer = wx.BoxSizer(wx.VERTICAL)

        #title
        label_main = wx.StaticText(self, label="Robot definition")
        main_sizer.Add(label_main, 0, wx.ALL | wx.ALIGN_CENTER_HORIZONTAL, 25)

        #grid
        grid = wx.GridBagSizer(15, 15)
        grid.Add(wx.StaticText(self, label='Name of the robot:'), pos=(0, 0),
                 flag=wx.BOTTOM | wx.ALIGN_LEFT, border=2)
        grid.Add(wx.TextCtrl(self, size=(92, -1), name='name', value=name),
                 pos=(0, 1))
        label = wx.StaticText(self, label='Number of moving links:')
        grid.Add(label, pos=(1, 0),
                 flag=wx.BOTTOM | wx.TOP | wx.ALIGN_LEFT, border=2)
        self.spin_links = wx.SpinCtrl(self, size=(92, -1),
                                      value=str(nl), min=1)
        self.spin_links.Bind(wx.EVT_SPINCTRL, self.OnSpinNL)
        grid.Add(self.spin_links, pos=(1, 1))

        label = wx.StaticText(self, label='Number of joints:')
        grid.Add(label, pos=(2, 0),
                 flag=wx.BOTTOM | wx.TOP | wx.ALIGN_LEFT, border=2)
        self.spin_joints = wx.SpinCtrl(self, size=(92, -1),
                                       value=str(nj), min=0)
        grid.Add(self.spin_joints, pos=(2, 1))

        grid.Add(wx.StaticText(self, label='Type of structure'), pos=(3, 0),
                 flag=wx.BOTTOM | wx.TOP | wx.ALIGN_LEFT, border=2)
        self.cmb_structure = wx.ComboBox(self, size=(92, -1),
                                         name='structure', style=wx.CB_READONLY,
                                         choices=[SIMPLE, TREE, CLOSED_LOOP],
                                         value=structure)
        grid.Add(self.cmb_structure, pos=(3, 1))
        self.cmb_structure.Bind(wx.EVT_COMBOBOX, self.OnTypeChanged)
        self.OnTypeChanged(None)

        main_sizer.Add(grid, 0, wx.ALL | wx.ALIGN_CENTER_HORIZONTAL, 15)
        self.ch_is_mobile = wx.CheckBox(self, label=' Is mobile')
        self.ch_is_mobile.Value = is_mobile
        self.ch_keep_geo = wx.CheckBox(self, label=' Keep geometric parameters')
        self.ch_keep_geo.Value = True
        self.ch_keep_dyn = wx.CheckBox(self, label=' Keep dynamic parameters')
        self.ch_keep_dyn.Value = True
        self.ch_keep_base = wx.CheckBox(self, label=' Keep base parameters')
        self.ch_keep_base.Value = True
        main_sizer.Add(self.ch_is_mobile, 0,
                       wx.LEFT | wx.RIGHT | wx.ALIGN_LEFT, 15)
        main_sizer.Add(self.ch_keep_geo, 0,
                       wx.TOP | wx.LEFT | wx.ALIGN_LEFT, 15)
        main_sizer.Add(self.ch_keep_dyn, 0,
                       wx.TOP | wx.LEFT | wx.ALIGN_LEFT, 15)
        main_sizer.Add(self.ch_keep_base, 0,
                       wx.TOP | wx.LEFT | wx.ALIGN_LEFT, 15)
        hor_sizer = wx.BoxSizer(wx.HORIZONTAL)
        ok_btn = wx.Button(self, wx.ID_OK, "OK")
        ok_btn.Bind(wx.EVT_BUTTON, self.OnOK)
        cancel_btn = wx.Button(self, wx.ID_CANCEL, "Cancel")
        cancel_btn.Bind(wx.EVT_BUTTON, self.OnCancel)
        hor_sizer.Add(ok_btn, 0, wx.ALL, 25)
        hor_sizer.Add(cancel_btn, 0, wx.ALL, 25)
        main_sizer.Add(hor_sizer)

        self.SetSizerAndFit(main_sizer)

    def OnOK(self, _):
        self.EndModal(wx.ID_OK)

    def OnCancel(self, _):
        self.EndModal(wx.ID_CANCEL)

    def OnTypeChanged(self, _):
        if self.cmb_structure.GetSelection() == 2:
            self.spin_joints.Enable(True)
        else:
            self.spin_joints.Enable(False)
            self.OnSpinNL(None)

    def OnSpinNL(self, _):
        self.spin_joints.SetRange(int(self.spin_links.Value), 100)
        if self.cmb_structure.GetSelection() != 2:
            self.spin_joints.Value = self.spin_links.Value

    def get_values(self):
        name = self.FindWindowByName('name').Value
        nl = int(self.spin_links.Value)
        nj = int(self.spin_joints.Value)
        return {'init_pars': (name, nl, nj, 2*nj - nl,
                               self.ch_is_mobile.Value,
                               self.cmb_structure.Value),
                'keep_geo': self.ch_keep_geo.Value,
                'keep_dyn': self.ch_keep_dyn.Value,
                'keep_base': self.ch_keep_base.Value}

reason: from sympy


text:
import os
import sys
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../src/')))
from permhull import *
from evenpermhull import *

def test_perms():
    from math import factorial
    # tests that generated groups are correct size
    for n in range(2,7):
        assert len(list(symmetric_group(n))) == factorial(n)
        assert len(list(alternating_group(n))) == factorial(n)//2

def test_cycle_types():
    assert len(list(even_cycle_types(6))) == 6
    assert len(list(even_cycle_types(7))) == 8
    num_partitions = [0,1,2,3,5,7,11,15,22,30,42,
                      56,77,101,135,176]
    num_even_cycle_types = [0,1,1,2,3,4,6,8]
    for n in range(1,15):
        cycles = cycle_types(n)
        assert len(list(cycles)) == num_partitions[n]
        assert all(len(lst)==n for lst in cycles)
        for i in range(len(list(cycles))):
            for j in range(i+1, len(list(cycles))):
                assert cycles[i] != cycles[j]
        even_cycles = even_cycle_types(n)
        for i in range(len(list(even_cycles))):
            for j in range(i+1, len(list(even_cycles))):
                assert not np.all(even_cycles[i] != even_cycles[j])
        assert all(len(lst)==n for lst in even_cycle_types(n))

def test_all_comb_coeff():
    from scipy.special import binom
    # checks that all tuples of convex coefficients sum to 1
    eps = 1e-14
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(2,1,.049))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(2,1,.1))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(2,1,.12))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(3,1,.1))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(3,1,.12))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(4,1,.1))
    assert all(abs(sum(c)-1)< eps for c in all_comb_coeff(4,1,.12))

def test_in_region():

    # PM_3
    x_ranges, lines = pm_boundary(3)
    assert in_region((.49 + .0j), x_ranges, lines)
    assert in_region((.49 + .2j), x_ranges, lines)
    assert not in_region((.49 - .5j), x_ranges, lines)
    assert not in_region((-.51 + .49j), x_ranges, lines)

    # PM_4
    x_ranges, lines = pm_boundary(4)
    assert in_region((.49 + .0j), x_ranges, lines)
    assert in_region((.49 + .2j), x_ranges, lines)
    assert not in_region((.49 - .52j), x_ranges, lines)
    assert in_region((-.51 + .49j), x_ranges, lines)

if __name__ == "__main__":
    test_perms()
    test_cycle_types()
    test_all_comb_coeff()
    test_in_region()
    print("All tests pass")

max_stars_repo_name: cptq/ds-spectra

max_issues_repo_name: cptq/ds-spectra

max_forks_repo_name: cptq/ds-spectra

reason: from scipy


text:
"""
Population dynamics module tests
"""

import os
import numpy as np
from scipy.ndimage import gaussian_filter
import popdyn as pd
from popdyn import summary
import shutil
import time

def some_random_k(shape, factor):
    k = np.zeros(shape=shape, dtype='float32')
    np.random.seed(2);
    i = np.random.randint(0, shape[0], 5)
    np.random.seed(2);
    j = np.random.randint(0, shape[1], 5)
    k[(i, j)] = 100000
    k = gaussian_filter(k, 10)
    np.random.seed(0);
    k += np.random.normal(0, 1000, shape)
    k = gaussian_filter(k, 10) * factor
    k[k < 0] = 0
    return k

# Domain
shape = (90, 160)

# Species
#================================================================================

# Stark
starks = pd.Species('Stark')

stark_male_infant = pd.AgeGroup('Stark', 'Infant', 'male', 0, 3)
stark_female_infant = pd.AgeGroup('Stark', 'Infant', 'female', 0, 3)

stark_male_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'male', 4, 12)
stark_female_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'female', 4, 12)

stark_male_adult = pd.AgeGroup('Stark', 'Adult', 'male', 13, 50)
stark_female_adult = pd.AgeGroup('Stark', 'Adult', 'female', 13, 50)

# Lannister
lannister = pd.Species('Lannister')

lannister_male_infant = pd.AgeGroup('Lannister', 'Infant', 'male', 0, 3)
lannister_female_infant = pd.AgeGroup('Lannister', 'Infant', 'female', 0, 3)

lannister_male_adolescent = pd.AgeGroup('Lannister', 'Adolescent', 'male', 4, 12)
lannister_female_adolescent = pd.AgeGroup('Lannister', 'Adolescent', 'female', 4, 12)

lannister_male_adult = pd.AgeGroup('Lannister', 'Adult', 'male', 13, 50)
lannister_female_adult = pd.AgeGroup('Lannister', 'Adult', 'female', 13, 50)

# White walker
white_walker = pd.Species('White Walker')

# Mortality
#================================================================================

# For humans
disease = pd.Mortality('Disease')
accident = pd.Mortality('Accident')
white_walker_death = pd.Mortality('White Walker')
white_walker_death.add_as_species(white_walker, [(0, 0), (0.1, 0.1), (0.5, 0.8), (1, 0.9)])

# Stark soldiers
lannister_male_soldier = pd.Mortality('Lannister Male')
lannister_male_soldier.add_as_species(lannister_male_adult, [(0, 0), (0.2, 0.05), (1, 0.1)])
lannister_female_soldier = pd.Mortality('Lannister Female')
lannister_female_soldier.add_as_species(lannister_female_adult, [(0, 0), (0.2, 0.01), (1, 0.05)])

# Lannister soldiers
stark_male_soldier = pd.Mortality('Stark Male')
stark_male_soldier.add_as_species(stark_male_adult, [(0, 0), (0.2, 0.05), (1, 0.1)])
stark_female_soldier = pd.Mortality('Stark Female')
stark_female_soldier.add_as_species(stark_male_adult, [(0, 0), (0.2, 0.05), (1, 0.1)])

# Carrying capacity
#================================================================================

stark_k_data = some_random_k(shape, 1.)
stark_k = pd.CarryingCapacity('Stark Habitat')
lannister_k_data = some_random_k(shape, 1.2)
lannister_k = pd.CarryingCapacity('Lannister Habitat')
white_walker_k = pd.CarryingCapacity('White Walker Habitat')
white_walker_k_data = some_random_k(shape, 2.)

# Fecundity
#================================================================================
male_fecundity = pd.Fecundity('males', multiplies=False)
female_fecundity = pd.Fecundity('females')

# Tests
#================================================================================

def no_species():
    """Should raise an exception in the solver error checker"""
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=1, left=0) as domain:
        pd.solvers.discrete_explicit(domain, 0, 1).execute()

def incorrect_ages():
    """Should raise an exception in the solver error checker"""
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=1, left=0) as domain:
        stark_female_infant = pd.AgeGroup('Stark', 'Infant', 'female', 0, 3)
        stark_female_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'female', 3, 12)  # Age range overlaps infants
        domain.add_population(stark_female_infant, 1, 0)
        domain.add_population(stark_female_adolescent, 1, 0)

        pd.solvers.discrete_explicit(domain, 0, 1).execute()

def incorrect_species():
    """Should raise an exception in the solver error checker"""
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=1, left=0) as domain:
        domain.add_population(stark_female_adolescent, 1, 0)
        # Add white walker as a form of mortality, but purposefully do not add white walkers to the domain
        domain.add_mortality(stark_female_adolescent, white_walker_death, time=0)

        pd.solvers.discrete_explicit(domain, 0, 1).execute()

def single_species():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # Population should remain consistent at 10,000 for time steps 1 and 2
        for i in range(1, 3):
            tot_pop = summary.total_population(domain, 'stark', i).sum()
            if not np.isclose(tot_pop, 10000):
                raise Exception('The population has changed to {} at time {}'.format(tot_pop, i))

def single_species_emigration():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)

        # Add population at time 1
        domain.add_population(starks, 10000., 1, distribute_by_habitat=True)

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # Population should 20,000 at time step 1
        tot_pop = summary.total_population(domain, 'stark', 1).sum()
        if not np.isclose(tot_pop, 20000):
            raise Exception('The population should be {}, got {}'.format(20000, tot_pop))

def single_species_mvp():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        stark_with_mvp = pd.Species('stark', minimum_viable_population=100, minimum_viable_area=shape[0] * shape[1])

        domain.add_carrying_capacity(stark_with_mvp, stark_k, 0, shape[0] * shape[1] * 100.)
        domain.add_population(stark_with_mvp, 50., 0)

        pd.solvers.discrete_explicit(domain, 0, 10).execute()

        # Population should be reduced
        tp = []
        for i in range(0, 11):
            tp.append(summary.total_population(domain, 'stark', i).sum())
        if not any([t == 0 for t in tp]):
            raise Exception('The population did not succumb to MVP (populations: {})'.format(tp))

def single_species_random_k():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        rand_k = pd.CarryingCapacity('Stark Habitat')
        rand_k_data = some_random_k(shape, 1.)
        rand_k.random('normal', **{'scale': 10})
        domain.add_carrying_capacity(starks, rand_k, 0, rand_k_data, distribute=False)
        domain.add_population(starks, rand_k_data, 0, distribute=False)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # Population is at k, so any k above the previous year should result in:
        #   1. density dependent mortality
        #   2. a reciprocal reduction in population
        prv_pop = rand_k_data
        for i in range(1, 3):
            # Make sure k is variable
            cc = summary.total_carrying_capacity(domain, 'stark', i)
            if np.all(cc == rand_k_data):
                raise Exception('The k is unchanged')
            tot_pop = summary.total_population(domain, 'stark', i)
            dd_mort = summary.total_mortality(domain, 'stark', i, mortality_name='Density Dependent')
            if not np.allclose(-dd_mort, (tot_pop - prv_pop), atol=1E-04):
                raise Exception('DD mortality not correct at time {}'.format(i))
            prv_pop = tot_pop

def single_species_sex():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)
        domain.add_population(pd.Sex('Stark', 'male'), 5000., 0, distribute_by_habitat=True)
        domain.add_population(pd.Sex('Stark', 'female'), 5000., 0, distribute_by_habitat=True)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # Population should remain consistent at 10,000 for time steps 1 and 2
        for i in range(1, 3):
            tot_pop = summary.total_population(domain, 'stark', i).sum()
            if not np.isclose(tot_pop, 20000):
                raise Exception('The total population has changed to {} at time {}'.format(tot_pop, i))
            tot_pop = summary.total_population(domain, 'stark', i, 'male').sum()
            if not np.isclose(tot_pop, 10000):
                raise Exception('The male population has changed to {} at time {}'.format(tot_pop, i))
            tot_pop = summary.total_population(domain, 'stark', i, 'female').sum()
            if not np.isclose(tot_pop, 10000):
                raise Exception('The female population has changed to {} at time {}'.format(tot_pop, i))
            cc = summary.total_carrying_capacity(domain, 'stark', i, 'male').sum()
            if not np.isclose(cc, (stark_k_data).sum()):
                raise Exception('The male carrying capacity has changed to {} at time {}'.format(cc, i))
            cc = summary.total_carrying_capacity(domain, 'stark', i, 'male').sum()
            if not np.isclose(cc, (stark_k_data).sum()):
                raise Exception('The female carrying capacity has changed to {} at time {}'.format(cc, i))

def single_species_fecundity():
    # F:M ratio allows offspring
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 1000., distribute=False)
        males = pd.Sex('Stark', 'male')
        females = pd.Sex('Stark', 'female')
        domain.add_population(males, 50, 0)
        domain.add_population(females, 50, 0)
        domain.add_fecundity(males, male_fecundity, 0, 1.1)
        domain.add_fecundity(females, female_fecundity, 0, 1.1)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # A constant set of offspring should result
        prv_pop = 50.
        for i in range(1, 3):
            tot_pop = summary.total_population(domain, 'stark', i, 'female').sum()
            tot_off = summary.total_offspring(domain, 'stark', i, 'female').sum()
            if not np.isclose(tot_off, prv_pop * 1.1):
                raise Exception(
                    'The offspring is {} and should be {} at time {}'.format(tot_off, prv_pop * 1.1, i)
                )
            if not np.isclose(tot_pop, prv_pop + (prv_pop * 1.1) / 2):
                raise Exception(
                    'The population is {} and should be {} at time {}'.format(
                        tot_pop, prv_pop + (prv_pop * 1.1) / 2, i
                    )
                )
            prv_pop = tot_pop

    if os.path.isfile('seven_kingdoms.popdyn'):
        os.remove('seven_kingdoms.popdyn')
    else:
        shutil.rmtree('seven_kingdoms.popdyn')

    # F:M ratio disallows offspring
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 60)

        stark_male = pd.Sex('Stark', 'male')
        stark_female = pd.Sex('Stark', 'female')

        domain.add_population(stark_male, 20., 0)
        domain.add_population(stark_female, 40., 0)

        test_female_fecundity_1 = pd.Fecundity('females', fecundity_lookup=[(0., 1.), (2., 0.)])
        # Explore density-dependence
        test_female_fecundity_2 = pd.Fecundity('females',
                                               density_fecundity_threshold=0.9,
                                               density_fecundity_max=1.1,
                                               fecundity_lookup=[(0., 1.), (1.5, 1. / 1.1)])
        domain.add_fecundity(stark_male, male_fecundity, 0, 1.0)
        domain.add_fecundity(stark_female, test_female_fecundity_1, 0, 1.1)

        domain.add_fecundity(stark_female, test_female_fecundity_2, 2, 1.1)

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        male_pop = summary.total_offspring(domain, 'stark', 1, 'male').sum()
        if male_pop != 0:
            raise Exception('Offspring is {} when should be 0'.format(male_pop))

        female_pop = summary.total_population(domain, 'stark', 2, 'female').sum()
        # Fecundity should be 0.5
        if not np.isclose(female_pop, 40 + (40 * 0.5 / 2)):
            raise Exception(
                'The population is {} and should be {}'.format(female_pop, 40 + (40 * 0.5 / 2))
            )

def single_species_dispersion():
    for _iter, method in enumerate(['density-based dispersion', 'distance propagation']):
        with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
            starks = pd.Species('Stark')

            starks.add_dispersal(method, (10,))

            # Avoid density-dependent mortality
            domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data + 1000., distribute=False)

            seed_population = np.zeros(shape=shape, dtype='float32')
            seed_population[(np.random.randint(0, shape[0], 10), np.random.randint(0, shape[1], 10))] = 100
            domain.add_population(starks, seed_population, 0, distribute=False)

            try:
                pd.solvers.discrete_explicit(domain, 0, 2).execute()
            except NotImplementedError:
                domain.file.close()
                if os.path.isfile('seven_kingdoms.popdyn'):
                    os.remove('seven_kingdoms.popdyn')
                else:
                    shutil.rmtree('seven_kingdoms.popdyn')
                continue

            prv_pop = seed_population
            for i in range(1, 3):
                _pop = summary.total_population(domain, 'stark', i)
                if not np.isclose(_pop.sum(), prv_pop.sum()):
                    raise Exception('The population changed from {} to {} at time {}'.format(prv_pop.sum(), _pop.sum() ,i))
                if np.allclose(_pop, prv_pop):
                    raise Exception('The population did not disperse using {}'.format(method))
                prv_pop = _pop

        if _iter == 0:
            if os.path.isfile('seven_kingdoms.popdyn'):
                os.remove('seven_kingdoms.popdyn')
            else:
                shutil.rmtree('seven_kingdoms.popdyn')

def single_species_agegroups():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:
        pops = [250., 250., 600., 600., 3000., 3000.]
        tot_pop = np.sum(pops)

        domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data + tot_pop, distribute=False)
        domain.add_population(stark_male_infant, 250., 0, distribute_by_habitat=True)
        domain.add_population(stark_female_infant, 250., 0, distribute_by_habitat=True)
        domain.add_population(stark_male_adolescent, 600., 0, distribute_by_habitat=True)
        domain.add_population(stark_female_adolescent, 600., 0, distribute_by_habitat=True)
        domain.add_population(stark_male_adult, 3000., 0, distribute_by_habitat=True)
        domain.add_population(stark_female_adult, 3000., 0, distribute_by_habitat=True)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        for i in range(1, 3):
            _pop = summary.total_population(domain, 'stark', i).sum()
            # Subtract old age mortality for the total population
            if not np.isclose(_pop, tot_pop - 3000. / (50 - 13. + 1) * 2 * i):
                raise Exception('The population should be {}, not {} at time {}'.format(
                    tot_pop - 3000. / (50 - 13. + 1) * 2, _pop, i)
                )
            # Check propagation
            _pop = summary.total_population(domain, 'stark', i, 'male', group=stark_male_infant.group_name).sum()
            if not np.isclose(_pop, 250 - ((250 / 4.) * i)):
                raise Exception('{} population should be {}, but is {} at time {}'.format(
                    stark_male_infant.group_name, 250 - ((250 / 4.) * i), _pop, i
                ))
            _pop = summary.total_population(domain, 'stark', i, 'male', group=stark_male_adolescent.group_name).sum()
            if not np.isclose(_pop, 600 - ((600 / 9.) * i) + ((250 / 4.) * i)):
                raise Exception('{} population should be {}, but is {} at time {}'.format(
                    stark_male_adolescent.group_name, 600 - ((600 / 9.) * i) + ((250 / 4.) * i), _pop, i
                ))
            _pop = summary.total_population(domain, 'stark', i, 'male', group=stark_male_adult.group_name).sum()
            if not np.isclose(_pop, 3000. - (3000. / (50 - 13. + 1) * i) + ((600 / 9.) * i)):
                raise Exception('{} population should be {}, but is {} at time {}'.format(
                    stark_male_adult.group_name, 3000. - (3000. / (50 - 13. + 1) * i) + ((600 / 9.) * i), _pop, i
                ))

def single_species_mask():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=shape, top=shape[0], left=0) as domain:

        starks = pd.Species('Stark')
        stark_male_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'male', 4, 12)
        stark_male_adult = pd.AgeGroup('Stark', 'Adult', 'male', 13, 50)
        stark_male_adolescent.add_dispersal('masked density-based dispersion', (10,))

        domain.add_carrying_capacity(starks, stark_k, 0, stark_k_data, distribute=False)
        domain.add_population(stark_male_adolescent, np.random.random(shape), 0, distribute=False)
        domain.add_population(stark_male_adult, np.random.random(shape), 0, distribute=False)
        domain.add_mask(starks, 0, np.random.random(shape), distribute=False)  # Should be inherited

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        for i in range(1, 2):
            _pop = summary.total_population(domain, 'stark', i).sum() + summary.total_mortality(domain, 'stark', i).sum()
            if not np.isclose(_pop, summary.total_population(domain, 'stark', 0).sum()):
                raise Exception('Population should be {}, got {}'.format(
                    summary.total_population(domain, 'stark', 0).sum(), _pop)
                )
            if np.allclose(summary.total_population(domain, 'stark', i) + summary.total_mortality(domain, 'stark', i),
                           summary.total_population(domain, 'stark', 0)):
                raise Exception('No dispersal occurred')

def single_species_mortality():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        # Test old age propagation here
        stark_male_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'male', 4, 12, live_past_max=True)
        stark_female_adolescent = pd.AgeGroup('Stark', 'Adolescent', 'female', 4, 12, live_past_max=True)

        domain.add_carrying_capacity(starks, stark_k, 0, 1000., distribute=False)
        domain.add_population(stark_male_adolescent, 10., 0)
        domain.add_population(stark_female_adolescent, 10., 0)
        domain.add_mortality(stark_male_adolescent, disease, 0, 0.2)
        domain.add_mortality(stark_female_adolescent, disease, 0, 0.2)
        domain.add_mortality(stark_female_adolescent, accident, 0, 0.1)
        domain.add_mortality(stark_male_adolescent, accident, 0, 0.1)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        for sex in ['male', 'female']:
            prv_pop = 10
            for i in range(1, 3):
                tot_pop = summary.total_population(domain, 'stark', i, sex, group='adolescent')
                mort1 = summary.total_mortality(domain, 'stark', i, sex, mortality_name='Accident')
                mort2 = summary.total_mortality(domain, 'stark', i, sex, mortality_name='Disease')
                if not np.isclose(tot_pop.sum(), prv_pop - (mort1.sum() + mort2.sum())):
                    raise Exception('{} adolescent mortality not correct at time {}'.format(sex, i))
                prv_pop = tot_pop

def single_species_recipient():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        recipient = pd.Species('recipient')

        test_mortality = pd.Mortality('test recipient')
        test_mortality.add_recipient_species(recipient)

        domain.add_carrying_capacity(starks, stark_k, 0, 10000, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)
        domain.add_mortality(starks, test_mortality, 0, 0.1)
        domain.add_population(recipient, 0, 10)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        tot_pop = summary.total_population(domain, 'recipient', 1).sum()
        if tot_pop != 10000 * .1:
            raise Exception('Recipient population should be {}, got {} at time 1'.format(
                10000 * .1, tot_pop)
            )
        tot_pop = summary.total_population(domain, 'recipient', 2).sum()
        if tot_pop != (10000 - (10000 * .1)) * .1:
            raise Exception('Recipient population should be {}, got {} at time 2'.format(
                (10000 - (10000 * .1)) * .1, tot_pop)
            )

def reverse_conversion_at_birth():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        recipient_males = pd.Sex('recipient', 'male')
        recipient_females = pd.Sex('recipient', 'female')
        test_mortality = pd.Mortality('test recipient')
        test_mortality.add_recipient_species(recipient_females)

        domain.add_population(recipient_males, 10, 0)
        domain.add_population(recipient_females, 10, 0)
        domain.add_fecundity(recipient_males, pd.Fecundity('test'), 0, .05)
        domain.add_fecundity(recipient_females, pd.Fecundity('test'), 0, .05)

        domain.add_carrying_capacity(starks, stark_k, 0, 10000, distribute=False)
        domain.add_population(pd.Sex('Stark', 'male'), 0.5, 0, distribute_by_habitat=True)
        domain.add_population(pd.Sex('Stark', 'female'), 0.5, 0, distribute_by_habitat=True)
        domain.add_mortality(starks, test_mortality, 0, 0.1)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        # There should be 2 starks at time 1
        tot_pop = summary.total_population(domain, 'stark', 1).sum()
        if tot_pop != 2:
            raise Exception('No offspring from recipient species in Starks at time 1')

def species_as_mortality():
    # White Walker mortality lookup: [(0, 0), (0.1, 0.1), (0.5, 0.8), (1, 0.9)]
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 100)
        domain.add_population(starks, 100, 0)
        domain.add_mortality(starks, white_walker_death, 0)

        # 30% density, mortality should be 0.45
        domain.add_carrying_capacity(white_walker, white_walker_k, 0, 100)
        domain.add_population(white_walker, 30., 0)

        # 50% density, mortality should be 0.8
        domain.add_carrying_capacity(white_walker, white_walker_k, 2, 60)

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        tot_pop = summary.total_population(domain, 'stark', 1).sum()
        if not np.isclose(tot_pop, 100 - (100 * 0.45)):
            raise Exception('Population is {}, but should be {} at time step 1'.format(tot_pop, 100 - (100 * 0.45)))
        tot_pop = summary.total_population(domain, 'stark', 2).sum()
        if not np.isclose(tot_pop, 55. - (55 * 0.8)):
            raise Exception('Population is {}, but should be {} at time step 1'.format(tot_pop, 55. - (55 * 0.8)))

def species_as_carrying_capacity():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 100)
        domain.add_carrying_capacity(lannister, lannister_k, 0, 100)
        domain.add_population(starks, 20, 0)  # Density should be 20%, increasing lannister habitat by 20%
        domain.add_population(lannister, 100, 0)

        stark_as_k = pd.CarryingCapacity('Starks')
        stark_as_k.add_as_species(starks, [(0., 1), (0.2, 1.2), (0.5, 0.9), (1., 0.1)])

        domain.add_carrying_capacity(lannister, stark_as_k, 0)

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        cc = summary.total_carrying_capacity(domain, 'lannister', 1).sum()
        if not np.isclose(cc, 120):
            raise Exception('Carrying capacity should be {}, but it is {}'.format(120, cc))

def circular_species():
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 100)
        domain.add_carrying_capacity(lannister, lannister_k, 0, 100)
        domain.add_population(starks, 20, 0)  # Density should be 20%, increasing lannister habitat by 20%
        domain.add_population(lannister, 100, 0)  # Density should be 100%, decreasing stark habitat by 90%

        stark_as_k = pd.CarryingCapacity('Starks')
        stark_as_k.add_as_species(starks, [(0., 1), (0.2, 1.2), (0.5, 0.9), (1., 0.1)])

        lannister_as_k = pd.CarryingCapacity('Lannister')
        lannister_as_k.add_as_species(lannister, [(0., 1), (0.2, 1.2), (0.5, 0.9), (1., 0.1)])

        domain.add_carrying_capacity(lannister, stark_as_k, 0)
        domain.add_carrying_capacity(starks, lannister_as_k, 0)

        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        if not np.isclose(summary.total_carrying_capacity(domain, 'lannister', 1).sum(), 120):
            raise Exception('Lannister carrying capacity should be {}, but it is {}'.format(120, cc))
        if not np.isclose(summary.total_carrying_capacity(domain, 'stark', 1).sum(), 10):
            raise Exception('Stark carrying capacity should be {}, but it is {}'.format(10, cc))

def max_age():
    with pd.Domain('seven_kingdoms.popdyn', csx=1, csy=1, shape=(1, 1), top=1, left=0) as domain:
        gp = pd.AgeGroup('Stark', 'Infant', 'male', 0, 3)
        gp.live_past_max = True
        domain.add_carrying_capacity(gp, stark_k, 0, 100)
        domain.add_population(gp, 20, 0)

        pd.solvers.discrete_explicit(domain, 0, 10).execute()

        # The age range of the group should have changed
        if domain.species['stark']['male']['infant'].max_age != 13:
            raise Exception('The max age is {}, and it should be 13'.format(domain.species['stark']['male']['infant'].max_age))

def global_n_interspecies():
    # White Walker mortality lookup: [(0, 0), (0.1, 0.1), (0.5, 0.8), (1, 0.9)]
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        domain.add_carrying_capacity(starks, stark_k, 0, 100)
        domain.add_population(starks, 100, 0)

        white_walker_death = pd.Mortality('White Walker')
        white_walker_death.add_as_species(white_walker, [(0, 0), (1/3., 0.1), (0.35714, 0.5), (1, 1)],
                                          population_type='global ratio')

        domain.add_mortality(starks, white_walker_death, 0)

        # should be 1/3, making mortality 0.1
        domain.add_carrying_capacity(white_walker, white_walker_k, 0, 100)
        domain.add_population(white_walker, 50., 0)

        pd.solvers.discrete_explicit(domain, 0, 2, global_density=True).execute()

        tot_pop = summary.total_population(domain, 'stark', 1).sum()
        if not np.isclose(tot_pop, 90):
            raise Exception('Population is {}, but should be {} at time step 1'.format(tot_pop, 90))
        tot_pop = summary.total_population(domain, 'stark', 2).sum()
        if not np.isclose(tot_pop, 90 * .5):
            raise Exception('Population is {}, but should be {} at time step 2'.format(tot_pop, 90 * .5))

def rate_based_mortality():
    # Pasted from recipient species
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        recipient = pd.Species('recipient')

        time_based_mortality = pd.Mortality('time based')
        time_based_mortality.add_time_based_mortality([0.2, 0.5, 1.])

        domain.add_mortality(recipient, time_based_mortality, 0)

        test_mortality = pd.Mortality('test recipient')
        test_mortality.add_recipient_species(recipient)

        domain.add_carrying_capacity(starks, stark_k, 0, 10000, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)
        domain.add_mortality(starks, test_mortality, 0, 0.1)
        pd.solvers.discrete_explicit(domain, 0, 5).execute()

        tot_shouldbe = [1000, 900, 810, 729, 656.10004]
        mort_shouldbe = [0, 200, 450, 810, 0]

        tot_pop, mort = [], []
        for year in range(1, 6):
            tot_pop.append(summary.total_population(domain, 'recipient', year).sum())
            mort.append(np.squeeze(summary.total_mortality(domain, 'recipient', year, mortality_name='time based')))

        if not np.allclose(tot_pop, tot_shouldbe) or not np.allclose(mort, mort_shouldbe):
            raise Exception('Population should be {}, got {}\nMortality should be {}, got {}'.format(
                tot_shouldbe, tot_pop, mort_shouldbe, mort)
            )

def conversion_and_transmission():
    # Pasted from recipient species
    with pd.Domain('seven_kingdoms.popdyn', csx=1., csy=1., shape=(1, 1), top=shape[0], left=0) as domain:
        recipient = pd.Species('recipient')

        test_mortality = pd.Mortality('test recipient')
        test_mortality.add_recipient_species(recipient)

        test_file = os.path.join(os.path.dirname(__file__), 'direct_transmission.csv')

        starks.add_disease(test_file, **{
            'direct_transmission': True,
            'environmental_transmission': True,
            'E_data': {
                1: np.array([[0.1]]),
                2: np.array([[0.1]])
            }
        })

        domain.add_carrying_capacity(starks, stark_k, 0, 10000, distribute=False)
        domain.add_population(starks, 10000., 0, distribute_by_habitat=True)
        domain.add_mortality(starks, test_mortality, 0, 0.1)
        domain.add_population(recipient, 0, 10)
        pd.solvers.discrete_explicit(domain, 0, 2).execute()

        Pi = 10 / 10000
        FOI = 0.1 * Pi
        env = 0.1 * 0.1
        conv = 10000 * (FOI + env)
        tot_pop = summary.total_population(domain, 'recipient', 1).sum()
        if not np.isclose(tot_pop, conv):
            raise Exception('Recipient population should be {}, got {} at time 1'.format(
                conv, tot_pop)
            )

if __name__ == '__main__':
    antitests = [no_species, incorrect_ages, incorrect_species]

    tests = [single_species, single_species_emigration, single_species_mvp, single_species_random_k,
             single_species_sex, single_species_fecundity, single_species_dispersion, single_species_agegroups,
             single_species_mask, single_species_mortality, single_species_recipient, species_as_mortality,
             species_as_carrying_capacity, circular_species, max_age, rate_based_mortality, global_n_interspecies,
             conversion_and_transmission, reverse_conversion_at_birth]

    error_check = 0
    now = time.time()
    for test in antitests:
        try:
            test()
            print("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n"
                  "{}: FAIL:\nNo Exception raised\n"
                  "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n".format(test.__name__))
        except Exception as e:
            if (isinstance(e, pd.solvers.SolverError) or
                isinstance(e, pd.PopdynError)):
                error_check += 1
                print("--------------------------------------\n"
                      "{}: PASS\n"
                      "--------------------------------------\n".format(test.__name__))
            else:
                print("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n"
                      "{}: FAIL:\n{}\n"
                      "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n".format(test.__name__, e))
        if os.path.isfile('seven_kingdoms.popdyn'):
            os.remove('seven_kingdoms.popdyn')
        else:
            shutil.rmtree('seven_kingdoms.popdyn')

    function_check = 0
    for test in tests:
        try:
            test()
            function_check += 1
            print("--------------------------------------\n"
                  "{}: PASS\n"
                  "--------------------------------------\n".format(test.__name__))
        except Exception as e:
            print("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n"
                  "{}: FAIL:\n{}\n"
                  "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n".format(test.__name__, e))
        if os.path.isfile('seven_kingdoms.popdyn'):
            os.remove('seven_kingdoms.popdyn')
        else:
            shutil.rmtree('seven_kingdoms.popdyn')

    stars = np.array([[' '] * 50 for i in range(10)])
    stars[(np.random.randint(0, 5, 40), np.random.randint(0, 49, 40))] = '*'
    stars[(np.random.randint(5, 10, 10), np.random.randint(0, 49, 10))] = '*'
    stars[:, -1] = '\n'
    stars = ''.join(stars.ravel())

    print("{}"
          "\n\n\n                Tests Completed in {} seconds\n\n"
          "".format(stars, time.time() - now))
    print("          {} of {} exception checks passed".format(error_check, len(antitests)))
    print("          {} of {} solver checks passed\n\n\n"
          "___wv_____wwWww____vW___|____wWv___|_|____w____wwvWWv".format(function_check, len(tests)))

max_stars_repo_name: ieg-tech/popdyn

max_issues_repo_name: ieg-tech/popdyn

max_forks_repo_name: ieg-tech/popdyn

reason: import numpy,from scipy


text:
from setuptools import find_packages
from numpy.distutils.core import setup, Extension
from json import load as json_load

ext1 = Extension(
        name = 'samos.lib.gaussian_density',
        sources = ['samos/lib/gaussian_density.f90'],
    )
ext2 = Extension(
        name = 'samos.lib.mdutils',
        sources = ['samos/lib/mdutils.f90'],
    )
ext3 = Extension(
        name = 'samos.lib.rdf',
        sources = ['samos/lib/rdf.f90'],
    )

if __name__ == '__main__':
    with open('setup.json', 'r') as info:
        kwargs = json_load(info)
    setup(
        include_package_data=True,
        packages=find_packages(),
        package_data = {'': ['*.f90']},
        ext_modules = [ext1, ext2, ext3],
        **kwargs
    )

max_stars_repo_name: ahzeeshan/samos

max_issues_repo_name: ahzeeshan/samos

max_forks_repo_name: ahzeeshan/samos

reason: from numpy


text:
# !/usr/bin/env python
# Created by "Thieu" at 14:51, 17/03/2020 ----------%
#       Email: nguyenthieu2102@gmail.com            %
#       Github: https://github.com/thieu1995        %
# --------------------------------------------------%

import numpy as np
from math import gamma
from copy import deepcopy
from mealpy.optimizer import Optimizer

class BaseHHO(Optimizer):
    """
    The original version of: Harris Hawks Optimization (HHO)

    Links:
        1. https://doi.org/10.1016/j.future.2019.02.028

    Examples
    ~~~~~~~~
    >>> import numpy as np
    >>> from mealpy.swarm_based.HHO import BaseHHO
    >>>
    >>> def fitness_function(solution):
    >>>     return np.sum(solution**2)
    >>>
    >>> problem_dict1 = {
    >>>     "fit_func": fitness_function,
    >>>     "lb": [-10, -15, -4, -2, -8],
    >>>     "ub": [10, 15, 12, 8, 20],
    >>>     "minmax": "min",
    >>> }
    >>>
    >>> epoch = 1000
    >>> pop_size = 50
    >>> model = BaseHHO(problem_dict1, epoch, pop_size)
    >>> best_position, best_fitness = model.solve()
    >>> print(f"Solution: {best_position}, Fitness: {best_fitness}")

    References
    ~~~~~~~~~~
    [1] Heidari, A.A., Mirjalili, S., Faris, H., Aljarah, I., Mafarja, M. and Chen, H., 2019.
    Harris hawks optimization: Algorithm and applications. Future generation computer systems, 97, pp.849-872.
    """

    def __init__(self, problem, epoch=10000, pop_size=100, **kwargs):
        """
        Args:
            problem (dict): The problem dictionary
            epoch (int): maximum number of iterations, default = 10000
            pop_size (int): number of population size, default = 100
        """
        super().__init__(problem, kwargs)
        self.epoch = self.validator.check_int("epoch", epoch, [1, 100000])
        self.pop_size = self.validator.check_int("pop_size", pop_size, [10, 10000])

        self.nfe_per_epoch = 1.5 * self.pop_size
        self.sort_flag = False

    def evolve(self, epoch):
        """
        The main operations (equations) of algorithm. Inherit from Optimizer class

        Args:
            epoch (int): The current iteration
        """
        pop_new = []
        for idx in range(0, self.pop_size):
            # -1 < E0 < 1
            E0 = 2 * np.random.uniform() - 1
            # factor to show the decreasing energy of rabbit
            E = 2 * E0 * (1 - (epoch + 1) * 1.0 / self.epoch)
            J = 2 * (1 - np.random.uniform())

            # -------- Exploration phase Eq. (1) in paper -------------------
            if np.abs(E) >= 1:
                # Harris' hawks perch randomly based on 2 strategy:
                if np.random.rand() >= 0.5:  # perch based on other family members
                    X_rand = deepcopy(self.pop[np.random.randint(0, self.pop_size)][self.ID_POS])
                    pos_new = X_rand - np.random.uniform() * np.abs(X_rand - 2 * np.random.uniform() * self.pop[idx][self.ID_POS])

                else:  # perch on a random tall tree (random site inside group's home range)
                    X_m = np.mean([x[self.ID_POS] for x in self.pop])
                    pos_new = (self.g_best[self.ID_POS] - X_m) - np.random.uniform() * \
                              (self.problem.lb + np.random.uniform() * (self.problem.ub - self.problem.lb))
                pos_new = self.amend_position(pos_new, self.problem.lb, self.problem.ub)
                pop_new.append([pos_new, None])
            # -------- Exploitation phase -------------------
            else:
                # Attacking the rabbit using 4 strategies regarding the behavior of the rabbit
                # phase 1: ----- surprise pounce (seven kills) ----------
                # surprise pounce (seven kills): multiple, short rapid dives by different hawks
                if (np.random.rand() >= 0.5):
                    delta_X = self.g_best[self.ID_POS] - self.pop[idx][self.ID_POS]
                    if np.abs(E) >= 0.5:  # Hard besiege Eq. (6) in paper
                        pos_new = delta_X - E * np.abs(J * self.g_best[self.ID_POS] - self.pop[idx][self.ID_POS])
                    else:  # Soft besiege Eq. (4) in paper
                        pos_new = self.g_best[self.ID_POS] - E * np.abs(delta_X)
                    pos_new = self.amend_position(pos_new, self.problem.lb, self.problem.ub)
                    pop_new.append([pos_new, None])
                else:
                    xichma = np.power((gamma(1 + 1.5) * np.sin(np.pi * 1.5 / 2.0)) /
                                      (gamma((1 + 1.5) * 1.5 * np.power(2, (1.5 - 1) / 2)) / 2.0), 1.0 / 1.5)
                    LF_D = 0.01 * np.random.uniform() * xichma / np.power(np.abs(np.random.uniform()), 1.0 / 1.5)
                    if np.abs(E) >= 0.5:  # Soft besiege Eq. (10) in paper
                        Y = self.g_best[self.ID_POS] - E * np.abs(J * self.g_best[self.ID_POS] - self.pop[idx][self.ID_POS])
                    else:  # Hard besiege Eq. (11) in paper
                        X_m = np.mean([x[self.ID_POS] for x in self.pop])
                        Y = self.g_best[self.ID_POS] - E * np.abs(J * self.g_best[self.ID_POS] - X_m)
                    pos_Y = self.amend_position(Y, self.problem.lb, self.problem.ub)
                    target_Y = self.get_target_wrapper(pos_Y)
                    Z = Y + np.random.uniform(self.problem.lb, self.problem.ub) * LF_D
                    pos_Z = self.amend_position(Z, self.problem.lb, self.problem.ub)
                    target_Z = self.get_target_wrapper(pos_Z)
                    if self.compare_agent([pos_Y, target_Y], self.pop[idx]):
                        pop_new.append([pos_Y, target_Y])
                        continue
                    if self.compare_agent([pos_Z, target_Z], self.pop[idx]):
                        pop_new.append([pos_Z, target_Z])
                        continue
                    pop_new.append(deepcopy(self.pop[idx]))
        self.pop = self.update_target_wrapper_population(pop_new)

max_stars_repo_name: rishavpramanik/mealpy

max_issues_repo_name: rishavpramanik/mealpy

max_forks_repo_name: rishavpramanik/mealpy

reason: import numpy


text:
import os
from random import shuffle

import cv2
import numpy as np
from sklearn.model_selection import train_test_split
from tensorflow.python.keras.applications import ResNet50
from tensorflow.python.keras.layers import Dense
from tensorflow.python.keras.models import Sequential
from tqdm import tqdm  # Helps in vis

TEST_SIZE = 0.5
RANDOM_STATE = 2018
BATCH_SIZE = 64
NO_EPOCHS = 5
NUM_CLASSES = 2
TRAIN_DIR = '/Users/administrator/workspace/ai_data/dog_vs_cat/Kaggle/train'
TEST_DIR = '/Users/administrator/workspace/ai_data/dog_vs_cat/test'
IMG_SIZE = 224

SHORT_LIST_TRAIN = os.listdir(TRAIN_DIR)[0:10000]  # using a subset of data as resouces as limited.
SHORT_LIST_TEST = os.listdir(TEST_DIR)

def label_img(img):
    word_label = img.split('.')[-3]
    # conversion to one-hot array [cat,dog]
    #                            [much cat, no dog]
    if word_label == 'cat':
        return [1, 0]
    #                             [no cat, very doggo]
    elif word_label == 'dog':
        return [0, 1]

def create_train_data():
    training_data = []
    for img in tqdm(SHORT_LIST_TRAIN):
        label = label_img(img)
        path = os.path.join(TRAIN_DIR, img)
        img = cv2.imread(path, cv2.IMREAD_COLOR)
        img = cv2.resize(img, (IMG_SIZE, IMG_SIZE))
        training_data.append([np.array(img), np.array(label)])
    shuffle(training_data)
    return training_data

def process_test_data():
    testing_data = []
    for img in tqdm(os.listdir(TEST_DIR)):
        path = os.path.join(TEST_DIR, img)
        img_num = img.split('.')[0]
        img = cv2.imread(path, cv2.IMREAD_COLOR)
        img = cv2.resize(img, (IMG_SIZE, IMG_SIZE))
        testing_data.append([np.array(img), img_num])
    shuffle(testing_data)
    return testing_data

train = create_train_data()
X = np.array([i[0] for i in train]).reshape(-1, IMG_SIZE, IMG_SIZE, 3)
y = np.array([i[1] for i in train])

my_new_model = Sequential()
my_new_model.add(ResNet50(include_top=False, pooling='max'))
my_new_model.add(Dense(NUM_CLASSES, activation='softmax'))

# Say not to train first layer (ResNet) model. It is already trained
my_new_model.layers[0].trainable = True

my_new_model.compile(optimizer='sgd', loss='categorical_crossentropy', metrics=['accuracy'])

my_new_model.summary()

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=TEST_SIZE, random_state=RANDOM_STATE)

train_model = my_new_model.fit(X_train, y_train,
                               batch_size=BATCH_SIZE,
                               epochs=NO_EPOCHS,
                               verbose=1,
                               validation_data=(X_val, y_val))

my_new_model.save('weights.h5')

score = my_new_model.evaluate(X_val, y_val, verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1])


text:
from data import *
from utils.augmentations import SSDAugmentation
from layers.modules import MultiBoxLoss
from ssd import build_ssd
import os
import re
import random
import sys
import time
import warnings
import torch
from torch.autograd import Variable
import torch.nn as nn
import torch.optim as optim
import torch.backends.cudnn as cudnn
import torch.nn.init as init
import torch.utils.data as data
import numpy as np
import argparse
import torch.multiprocessing as mp
import torch.utils.data.distributed
import torch.distributed as dist
from torch.cuda.amp import autocast, GradScaler

from tensorboardX import SummaryWriter

import sys
cur_dir = os.path.dirname(os.path.abspath(__file__))
sys.path.append(cur_dir + "/../../tools/utils/")
from metric import MetricCollector

def str2bool(v):
    return v.lower() in ("yes", "true", "t", "1")

parser = argparse.ArgumentParser(
    description='Single Shot MultiBox Detector Training With Pytorch')
train_set = parser.add_mutually_exclusive_group()
parser.add_argument('--dataset', default='VOC', choices=['VOC', 'COCO'],
                    type=str, help='VOC or COCO')
parser.add_argument('--dataset_root', default=VOC_ROOT,
                    help='Dataset root directory path')
parser.add_argument('--basenet', default='vgg16_reducedfc.pth',
                    help='Pretrained base model')
parser.add_argument('--batch_size', default=32, type=int,
                    help='Batch size for training')
parser.add_argument('--resume', default=None, type=str,
                    help='Checkpoint state_dict file to resume training from')
parser.add_argument('--start_iter', default=0, type=int,
                    help='Resume training at this iter')
parser.add_argument('--num_workers', default=4, type=int,
                    help='Number of workers used in dataloading')
parser.add_argument('--device', default='mlu', type=str,
                    help='Use MLU or GPU to train model, defaultly be MLU')
parser.add_argument('--lr', '--learning-rate', default=1e-3, type=float,
                    help='initial learning rate')
parser.add_argument('--momentum', default=0.9, type=float,
                    help='Momentum value for optim')
parser.add_argument('--weight_decay', default=5e-4, type=float,
                    help='Weight decay for SGD')
parser.add_argument('--gamma', default=0.1, type=float,
                    help='Gamma update for SGD')
parser.add_argument('--visdom', default=False, type=str2bool,
                    help='Use visdom for loss visualization')
parser.add_argument('--save_folder', default='weights/',
                    help='Directory for saving checkpoint models')
parser.add_argument('--seed', default=42, type=int)
parser.add_argument('--iters', type=int, default=1000, metavar='N',
                    help='train how many iterations totaly')
parser.add_argument('--world-size', default=-1, type=int,
                    help='number of nodes for distributed training')
parser.add_argument('--rank', default=-1, type=int,
                    help='node rank for distributed training')
parser.add_argument('--dist_url', default='tcp://127.0.0.10:28500', type=str,
                    help='url used to set up distributed training')
parser.add_argument('--device_id', default=None, type=int,
                    help='Use specified device for training, useless in '
                         'multiprocessing distributed training')
parser.add_argument('--multiprocessing-distributed', action='store_true',
                    help='Use multi-processing distributed training to launch '
                         'N processes per node, which has N GPUs. This is the '
                         'fastest way to use PyTorch for either single node or '
                         'multi node data parallel training')
parser.add_argument('--dist-backend', default='nccl', type=str,
                    help='distributed backend')
parser.add_argument('--pyamp', action='store_true', default=False,
                    help='use pytorch amp for mixed precision training')
args = parser.parse_args()

if args.device == 'mlu':
    import torch_mlu
    import torch_mlu.core.mlu_model as ct
    import torch_mlu.core.mlu_quantize as qt
    print("Use MLU Deivce ......")

if torch.cuda.is_available():
    if args.device == 'gpu':
        # There are no such APIs like torch.mlu.FloatTensor, in order to
        # maintains the consistensy before we support torch.mlu.FloatTensor,
        # we temporarily disable the torch.cuda.FloatTensor.
        #torch.set_default_tensor_type('torch.cuda.FloatTensor')
        torch.set_default_tensor_type('torch.FloatTensor')
    if args.device != 'gpu':
        print("WARNING: It looks like you have a CUDA device, but aren't " +
              "using CUDA.\nRun with --cuda for optimal training speed.")
        torch.set_default_tensor_type('torch.FloatTensor')
else:
    torch.set_default_tensor_type('torch.FloatTensor')

def load_amp_scaler_state(scaler, base_file):
    other, ext = os.path.splitext(base_file)
    if ext == '.pkl' or '.pth':
        ckpt = torch.load(base_file, map_location=lambda storage, loc: storage)
        if isinstance(ckpt, dict) and 'amp' in ckpt:
            print('Loading the scaler state of AMP...')
            scaler.load_state_dict(ckpt['amp'])
            print('Finished!')

def main():
    print('Training SSD on:', args.dataset)
    print('Using the specified args:')
    print(args)
    f = open("param_train.txt", "w")
    f.write(str(args)+"\n")
    f.close()

    if args.device_id is not None:
        warnings.warn('You have chosen a specific device. This will completely '
                      'disable data parallelism.')

    if args.dist_url == "env://" and args.world_size == -1:
        args.world_size = int(os.environ["WORLD_SIZE"])
    args.distributed = args.multiprocessing_distributed or args.world_size > 1

    ndevs_per_node = ct.device_count() if args.device == 'mlu' else torch.cuda.device_count()
    if args.multiprocessing_distributed:
        args.world_size = ndevs_per_node * args.world_size
        mp.spawn(train, nprocs=ndevs_per_node, args=(ndevs_per_node, args))
    else:
        train(args.device_id, ndevs_per_node, args)

def train(dev_id, ndevs_per_node, args):
    batch_time = AverageMeter('Time' , ':6.3f')
    losses     = AverageMeter('Loss' , ':.4e' )
    if args.seed is not None:
        random.seed(args.seed)
        np.random.seed(args.seed)
        torch.manual_seed(args.seed)
        if torch.cuda.is_available():
            torch.cuda.manual_seed(args.seed)
            if os.getenv('BENCHMARK_LOG') is None:
                torch.backends.cudnn.deterministic = True
                torch.backends.cudnn.benchmark = False
            else:
                torch.backends.cudnn.benchmark = True

    args.device_id = dev_id
    if args.device_id is not None:
        print("Use Device: {} for training".format(args.device_id))

    if args.device == 'mlu':
        ct.set_device(args.device_id)
    elif args.device == 'gpu':
        torch.cuda.set_device(args.device_id)
    if args.distributed:
        if args.dist_url == "env://" and args.rank == -1:
            args.rank = int(os.environ["RANK"])
        if args.multiprocessing_distributed:
            args.rank = args.rank * ndevs_per_node + dev_id
        dist.init_process_group(backend=args.dist_backend, init_method=args.dist_url,
            world_size=args.world_size, rank=args.rank)

    if args.dataset == 'COCO':
        if args.dataset_root == VOC_ROOT:
            if not os.path.exists(COCO_ROOT):
                parser.error('Must specify dataset_root if specifying dataset')
            print("WARNING: Using default COCO dataset_root because " +
                  "--dataset_root was not specified.")
            args.dataset_root = COCO_ROOT
        cfg = coco
        dataset = COCODetection(root=args.dataset_root,
                                transform=SSDAugmentation(cfg['min_dim'],
                                                          MEANS))
    elif args.dataset == 'VOC':
        if args.dataset_root == COCO_ROOT:
            parser.error('Must specify dataset if specifying dataset_root')
        cfg = voc
        dataset = VOCDetection(root=args.dataset_root,
                               transform=SSDAugmentation(cfg['min_dim'],
                                                         MEANS))

    if args.distributed:
        args.batch_size = int(args.batch_size / ndevs_per_node)
        args.num_workers = int((args.num_workers + ndevs_per_node - 1) / ndevs_per_node)
        train_sampler = torch.utils.data.distributed.DistributedSampler(dataset)
    else:
        train_sampler = None

    if args.visdom:
        import visdom
        viz = visdom.Visdom()

    ssd_net = build_ssd('train', cfg['min_dim'], cfg['num_classes'])
    net = ssd_net

    if args.resume:
        print('Resuming training, loading {}...'.format(args.resume))
        ssd_net.load_weights(args.resume)
        if args.pyamp:
            scaler = GradScaler()
            load_amp_scaler_state(scaler, args.resume)
    else:
        vgg_weights = torch.load(args.save_folder + args.basenet)
        print('Loading base network...')
        ssd_net.vgg.load_state_dict(vgg_weights)

    if not args.resume:
        print('Initializing weights...')
        # initialize newly added layers' weights with xavier method
        ssd_net.extras.apply(weights_init)
        ssd_net.loc.apply(weights_init)
        ssd_net.conf.apply(weights_init)

    if args.device == 'gpu':
        net = net.cuda()
        net.priors = net.priors.cuda()
    elif args.device == 'mlu':
        iters_per_epoch = len(dataset) // args.batch_size
        if train_sampler is not None:
            iters_per_epoch = iters_per_epoch // ndevs_per_node
        print("iters_per_epoch is: ", iters_per_epoch)
        net.to(ct.mlu_device())
        net.priors = net.priors.to('mlu')
    else:
        print('using CPU, this will be slow')

    if args.device != 'cpu' and args.distributed:
        net = torch.nn.parallel.DistributedDataParallel(net, device_ids=[args.device_id])

    optimizer = optim.SGD(net.parameters(), lr=args.lr, momentum=args.momentum,
                          weight_decay=args.weight_decay)
    criterion = MultiBoxLoss(cfg['num_classes'], 0.5, True, 0, True, 3, 0.5,
                             False, args.device)
    if args.device == 'mlu':
        ct.to(optimizer, torch.device('mlu'))

    net.train()
    # loss counters
    loc_loss = 0
    conf_loss = 0
    epoch = 0
    print('Loading the dataset...')

    epoch_size = len(dataset) // args.batch_size

    step_index = 0

    if args.visdom:
        vis_title = 'SSD.PyTorch on ' + dataset.name
        vis_legend = ['Loc Loss', 'Conf Loss', 'Total Loss']
        iter_plot = create_vis_plot('Iteration', 'Loss', vis_title, vis_legend)
        epoch_plot = create_vis_plot('Epoch', 'Loss', vis_title, vis_legend)

    data_loader = data.DataLoader(dataset, args.batch_size,
                                  num_workers=args.num_workers,
                                  sampler=train_sampler,
                                  shuffle=(train_sampler is None),
                                  collate_fn=detection_collate,
                                  pin_memory=True)

    #writer = SummaryWriter('runs/fp')
    # create batch iterator
    #batch_iterator = iter(data_loader)
    iteration = args.start_iter
    cfg['max_iter'] = iteration + args.iters

    # for internal benchmark test
    metric_collector = MetricCollector(
        enable_only_benchmark=True,
        record_elapsed_time=True,
        record_hardware_time=True if args.device == 'mlu' else False)
    metric_collector.place()

    while iteration < cfg['max_iter']:
        if args.visdom and iteration != 0 and (iteration % epoch_size == 0):
            update_vis_plot(epoch, loc_loss, conf_loss, epoch_plot, None,
                            'append', epoch_size)
            # reset epoch loss counters
            loc_loss = 0
            conf_loss = 0
        if args.distributed:
            train_sampler.set_epoch(epoch)
        t1 = time.time()
        for _, (images, targets) in enumerate(data_loader):
            iteration += 1
            if iteration in cfg['lr_steps']:
                step_index += 1
            adjust_learning_rate(optimizer, args.gamma, iteration, step_index)

            if args.device == 'gpu':
                images = images.cuda()
                targets = [ann.cuda() for ann in targets]
            elif args.device == 'mlu':
                images = images.to("mlu")
                targets = [ann.to("mlu") for ann in targets]
            else:
                images = images
                targets = [ann for ann in targets]
            # forward
            t0 = time.time()

            with autocast(enabled=args.pyamp):
                out = net(images)
                loss_l, loss_c = criterion(out, targets)
                loss = loss_l + loss_c

            optimizer.zero_grad()
            metric_collector.record()
            metric_collector.place()
            losses.update(loss.item(), images.size(0))
            # backprop
            if args.pyamp:
                scaler.scale(loss).backward()
                scaler.step(optimizer)
                scaler.update()
            else:
                loss.backward()
                optimizer.step()
            loc_loss += loss_l.item()
            conf_loss += loss_c.item()
            batch_time.update(time.time() - t1)
            t1 = time.time()

            #writer.add_scalar('train/loss', loss.item(), iteration)
            if iteration % 1 == 0:
                print('proc id: ' + str(args.device_id) + ' || iter ' +
                      repr(iteration) + ' || Loss: %.4f ||' % (loss.item()) +
                      ' timer: %.4f sec.' % (t1 - t0))

            if iteration == cfg['max_iter']:
                if (args.distributed == False) or (args.rank == 0):
                    if args.distributed:
                        state_dict = net.module.state_dict()
                    else:
                        state_dict = net.state_dict()
                    if args.pyamp:
                        state_dict = {"state_dict":state_dict, "amp":scaler.state_dict()}
                    if args.device == 'mlu':
                        save_pre = "mlu_weights"
                    elif args.device == 'gpu':
                        save_pre = "gpu_weights"
                    else:
                        save_pre = "cpu_weights"
                    print('Saving state, iter:', iteration)
                    if not os.path.exists(save_pre):
                        os.makedirs(save_pre)
                    torch.save(state_dict, save_pre + '/ssd300_COCO_' +
                        repr(iteration) + '.pth')
                break
        epoch += 1
    if ((args.distributed == False) or (args.rank == 0)) and os.getenv('AVG_LOG'):
        with open(os.getenv('AVG_LOG'), 'a') as train_avg:
            train_avg.write('net:SSD-VGG16, iter:{}, cards:{}, avg_loss:{}, avg_time:{}, '.
            format(-1 if args.iters == 120000 else args.iters, args.world_size if args.distributed else 1,
            losses.avg, batch_time.avg))
    metric_collector.insert_metrics(
        net = "SSD-VGG16",
        batch_size = args.batch_size,
        precision = "amp" if args.pyamp else "fp32",
        cards = args.world_size if args.distributed else 1,
        DPF_mode = "ddp " if args.distributed else "single")
    if ((args.distributed == False) or (args.rank == 0)):
        metric_collector.dump()

def adjust_learning_rate(optimizer, gamma, iteration, step):
    """Sets the learning rate to the initial LR decayed by 10 at every
        specified step
    # Adapted from PyTorch Imagenet example:
    # https://github.com/pytorch/examples/blob/master/imagenet/main.py
    """
    if iteration < 1000:
        lr = args.lr * (iteration // 258 + 1) / (1000 // 258 + 1)
    else:
        lr = args.lr * (gamma ** (step))
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr

def xavier(param):
    init.xavier_uniform(param)

def weights_init(m):
    if isinstance(m, nn.Conv2d):
        xavier(m.weight.data)
        m.bias.data.zero_()

def create_vis_plot(_xlabel, _ylabel, _title, _legend):
    return viz.line(
        X=torch.zeros((1,)).cpu(),
        Y=torch.zeros((1, 3)).cpu(),
        opts=dict(
            xlabel=_xlabel,
            ylabel=_ylabel,
            title=_title,
            legend=_legend
        )
    )

def update_vis_plot(iteration, loc, conf, window1, window2, update_type,
                    epoch_size=1):
    viz.line(
        X=torch.ones((1, 3)).cpu() * iteration,
        Y=torch.Tensor([loc, conf, loc + conf]).unsqueeze(0).cpu() / epoch_size,
        win=window1,
        update=update_type
    )
    # initialize epoch plot on first iteration
    if iteration == 0:
        viz.line(
            X=torch.zeros((1, 3)).cpu(),
            Y=torch.Tensor([loc, conf, loc + conf]).unsqueeze(0).cpu(),
            win=window2,
            update=True
        )

class AverageMeter(object):
    """Computes and stores the average and current value"""
    def __init__(self, name, fmt=':f'):
        self.name = name
        self.fmt = fmt
        self.reset()

    def reset(self):
        self.val = 0
        self.avg = 0
        self.sum = 0
        self.count = 0

    def update(self, val, n=1):
        self.val = val
        self.sum += val * n
        self.count += n
        self.avg = self.sum / self.count

    def __str__(self):
        fmtstr = '{name} {val' + self.fmt + '} ({avg' + self.fmt + '})'
        return fmtstr.format(**self.__dict__)

if __name__ == '__main__':
    main()


text:
import pandas as pd
from pandas.io.json import json_normalize
import json
import numpy as np

dataset_file = open('Dataset_SimpleQA.json','r')
dataset_decode = json.load(dataset_file)
tab = json_normalize(dataset_decode)

distribution = pd.read_csv("relation_distribution.csv",sep='|').head(50)

for rel in distribution['rel']:
	tab_questions = tab[tab.relation == rel]['question']
	question = ""
	for q in tab_questions:
		question += q + " "
	wordlist = question.split()
	word_unique_list = np.unique(np.array(wordlist))
	term_count = pd.DataFrame(columns = ['word','frequency'])
	i = 0
	for w in word_unique_list:
		frames = pd.DataFrame({'word':w,'frequency': wordlist.count(w)},index=[i])
		term_count = pd.concat([term_count, frames],axis=0)
		i =+ 1
	filepath= "Term_Count/" + rel + ".csv"
	term_count.to_csv(filepath,sep = '|',index=False,encoding='utf-8')

max_stars_repo_name: AskNowQA/VANiLLa

max_issues_repo_name: AskNowQA/VANiLLa

max_forks_repo_name: AskNowQA/VANiLLa


text:
import numpy as np

def rbf_kernel_distance(x, y, gamma=0.1) -> float:
    """ RBF 距离

    Args:
        x ([type]): [description]
        y ([type]): [description]
        gamma (int, optional): [description]. Defaults to 1.

    Returns:
        float: [description]
    """
    euclid_dis = np.sum((x-y)**2)
    d = 2 * gamma ** 2
    return 1 / np.exp(euclid_dis / d)

def euclid_distance(x, y):
    """计算欧几里得距离

    Args:
        x ([type]): [description]
        y ([type]): [description]

    Returns:
        [type]: [description]
    """
    return np.sum((x-y)**2)

def cal_distance_matrix(X):
    """生成距离矩阵

    Args:
        X (numpy.ndarray): 一组样本点

    Returns:
        numpy.ndarray: 样本点构成的距离矩阵
    """
    X = np.array(X)
    S = np.zeros((len(X), len(X)))
    for i in range(len(X)):
        for j in range(i+1, len(X)):
            S[i][j] = 1.0 * euclid_distance(X[i], X[j])
            # S[i][j] = 1.0 * rbf_kernel_distance(X[i], X[j])
            S[j][i] = S[i][j]
    return S

def cal_laplacian_matrix(adjacentMatrix):

    # compute the Degree Matrix: D=sum(A)
    degreeMatrix = np.sum(adjacentMatrix, axis=1)

    # compute the Laplacian Matrix: L=D-A
    laplacianMatrix = np.diag(degreeMatrix) - adjacentMatrix

    # normailze
    # D^(-1/2) L D^(-1/2)
    sqrtDegreeMatrix = np.diag(1.0 / (degreeMatrix ** (0.5)))
    return np.dot(np.dot(sqrtDegreeMatrix, laplacianMatrix), sqrtDegreeMatrix)


text:
from utils import symetric_random , true_positive , false_positive , precision , false_negative , recall , F1score , logical_and
import numpy as np

actual =   np.array([ 0 , 1 , 1 , 0 , 1 , 0 , 0 , 0])
estimated = np.array([0 , 1 , 0 , 1 , 1 , 1 , 0 , 1])

def test_F2Score():
    print(true_positive(actual,estimated) == 2)

    print(false_positive(actual,estimated) == 3)

    print(false_negative(actual,estimated) == 1)

    expected_prec =  2/(2+3)
    print(precision(actual,estimated) == expected_prec)

    expected_recall = 2/(2+1)
    print(recall(actual,estimated) == expected_recall)

    print(F1score(actual,estimated) == 2*expected_recall*expected_prec/(expected_recall+expected_prec))

def test_logical_and():
    print(logical_and(0,0) == 0)
    print(logical_and(1,0) == 0)
    print(logical_and(0,1) == 0)
    print(logical_and(1,1) == 1)

def test_symetric_random():
    e = 0.5
    i = 0
    while i < 100:
        r = symetric_random(e)
        if r <= e and r >= -e:
            i += 1
        else:
            print (False)
            return

    return True

# main
test_F2Score()
test_logical_and()
test_symetric_random()


text:
import copy
import os
import random

import numpy as np

sep = os.sep

####################### GLOBAL PARAMETERS ##################################################
############################################################################################
Params = {
    'num_channels': 1,
    'num_classes': 2,
    'batch_size': 4,
    'epochs': 251,
    'learning_rate': 0.001,
    'patch_shape': (388, 388),
    'patch_offset': (150, 150),
    'expand_patch_by': (184, 184),
    'use_gpu': True,
    'distribute': True,
    'shuffle': True,
    'log_frequency': 5,
    'validation_frequency': 1,
    'mode': 'train',
    'parallel_trained': False,
}

dparm_1_100_1 = lambda x: np.random.choice(np.arange(1, 101, 1), 2)
dparm_1_1 = lambda x: [1, 1]
d_parm_weighted = lambda x: [x['Params']['cls_weights'][0], x['Params']['cls_weights'][1]]
##############################################################################################

DRIVE = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'DRIVE' + sep + 'images',
        'mask': 'data' + sep + 'DRIVE' + sep + 'mask',
        'truth': 'data' + sep + 'DRIVE' + sep + 'manual',
        'splits_json': 'data' + sep + 'DRIVE' + sep + 'splits'
    },

    'Funcs': {
        'truth_getter': lambda file_name: file_name.split('_')[0] + '_manual1.gif',
        'mask_getter': lambda file_name: file_name.split('_')[0] + '_mask.gif'
    }
}

DRIVE_1_100_1 = copy.deepcopy(DRIVE)
DRIVE_1_100_1['Dirs']['logs'] = 'logs' + sep + 'DRIVE' + sep + 'UNET_1_100_1'
DRIVE_1_100_1['Funcs']['dparm'] = dparm_1_100_1

DRIVE_1_1 = copy.deepcopy(DRIVE)
DRIVE_1_1['Dirs']['logs'] = 'logs' + sep + 'DRIVE' + sep + 'UNET_1_1'
DRIVE_1_1['Funcs']['dparm'] = dparm_1_1

DRIVE_WEIGHTED = copy.deepcopy(DRIVE)
DRIVE_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'DRIVE' + sep + 'UNET_WEIGHTED'
DRIVE_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# --------------------------------------------------------------------------------------------

WIDE = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'AV-WIDE' + sep + 'images',
        'truth': 'data' + sep + 'AV-WIDE' + sep + 'manual',
        'splits_json': 'data' + sep + 'AV-WIDE' + sep + 'splits'
    },

    'Funcs': {
        'truth_getter': lambda file_name: file_name.split('.')[0] + '_vessels.png',
        'mask_getter': None
    }
}

WIDE_1_100_1 = copy.deepcopy(WIDE)
WIDE_1_100_1['Dirs']['logs'] = 'logs' + sep + 'AV_WIDE' + sep + 'UNET_1_100_1'
WIDE_1_100_1['Funcs']['dparm'] = dparm_1_100_1

WIDE_1_1 = copy.deepcopy(WIDE)
WIDE_1_1['Dirs']['logs'] = 'logs' + sep + 'AV_WIDE' + sep + 'UNET_1_1'
WIDE_1_1['Funcs']['dparm'] = dparm_1_1

WIDE_WEIGHTED = copy.deepcopy(WIDE)
WIDE_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'AV_WIDE' + sep + 'UNET_WEIGHTED'
WIDE_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# ---------------------------------------------------------------------------------------------

STARE = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'STARE' + sep + 'stare-images',
        'truth': 'data' + sep + 'STARE' + sep + 'labels-ah',
        'splits_json': 'data' + sep + 'STARE' + sep + 'splits'
    },

    'Funcs': {
        'truth_getter': lambda file_name: file_name.split('.')[0] + '.ah.pgm',
        'mask_getter': None
    }
}

STARE_1_100_1 = copy.deepcopy(STARE)
STARE_1_100_1['Dirs']['logs'] = 'logs' + sep + 'STARE' + sep + 'UNET_1_100_1'
STARE_1_100_1['Funcs']['dparm'] = dparm_1_100_1

STARE_1_1 = copy.deepcopy(STARE)
STARE_1_1['Dirs']['logs'] = 'logs' + sep + 'STARE' + sep + 'UNET_1_1'
STARE_1_1['Funcs']['dparm'] = dparm_1_1

STARE_WEIGHTED = copy.deepcopy(STARE)
STARE_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'STARE' + sep + 'UNET_WEIGHTED'
STARE_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# ------------------------------------------------------------------------------------------------

CHASEDB = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'CHASEDB' + sep + 'images',
        'truth': 'data' + sep + 'CHASEDB' + sep + 'manual',
        'splits_json': 'data' + sep + 'CHASEDB' + sep + 'splits'
    },

    'Funcs': {
        'truth_getter': lambda file_name: file_name.split('.')[0] + '_1stHO.png',
        'mask_getter': None
    }
}

CHASEDB_1_100_1 = copy.deepcopy(CHASEDB)
CHASEDB_1_100_1['Dirs']['logs'] = 'logs' + sep + 'CHASEDB' + sep + 'UNET_1_100_1'
CHASEDB_1_100_1['Funcs']['dparm'] = dparm_1_100_1

CHASEDB_1_1 = copy.deepcopy(CHASEDB)
CHASEDB_1_1['Dirs']['logs'] = 'logs' + sep + 'CHASEDB' + sep + 'UNET_1_1'
CHASEDB_1_1['Funcs']['dparm'] = dparm_1_1

CHASEDB_WEIGHTED = copy.deepcopy(CHASEDB)
CHASEDB_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'CHASEDB' + sep + 'UNET_WEIGHTED'
CHASEDB_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# -------------------------------------------------------------------------------------------------

VEVIO_MOSAICS = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'VEVIO' + sep + 'mosaics',
        'mask': 'data' + sep + 'VEVIO' + sep + 'mosaics_masks',
        'truth': 'data' + sep + 'VEVIO' + sep + 'mosaics_manual_01_bw',
        'splits_json': 'data' + sep + 'VEVIO' + sep + 'splits_mosaics'
    },

    'Funcs': {
        'truth_getter': lambda file_name: 'bw_' + file_name.split('.')[0] + '_black.' + file_name.split('.')[1],
        'mask_getter': lambda file_name: 'mask_' + file_name
    }
}

VEVIO_MOSAICS_1_100_1 = copy.deepcopy(VEVIO_MOSAICS)
VEVIO_MOSAICS_1_100_1['Dirs']['logs'] = 'logs' + sep + 'VEVIO_MOSAICS' + sep + 'UNET_1_100_1'
VEVIO_MOSAICS_1_100_1['Funcs']['dparm'] = dparm_1_100_1

VEVIO_MOSAICS_1_1 = copy.deepcopy(VEVIO_MOSAICS)
VEVIO_MOSAICS_1_1['Dirs']['logs'] = 'logs' + sep + 'VEVIO_MOSAICS' + sep + 'UNET_1_1'
VEVIO_MOSAICS_1_1['Funcs']['dparm'] = dparm_1_1

VEVIO_MOSAICS_WEIGHTED = copy.deepcopy(VEVIO_MOSAICS)
VEVIO_MOSAICS_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'VEVIO_MOSAICS' + sep + 'UNET_WEIGHTED'
VEVIO_MOSAICS_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# ---------------------------------------------------------------------------------------------------------

VEVIO_FRAMES = {
    'Params': Params,
    'Dirs': {
        'image': 'data' + sep + 'VEVIO' + sep + 'frames',
        'mask': 'data' + sep + 'VEVIO' + sep + 'frames_masks',
        'truth': 'data' + sep + 'VEVIO' + sep + 'frames_manual_01_bw',
        'splits_json': 'data' + sep + 'VEVIO' + sep + 'splits_frames'
    },

    'Funcs': {
        'truth_getter': lambda file_name: 'bw_' + file_name.split('.')[0] + '_black.' + file_name.split('.')[1],
        'mask_getter': lambda file_name: 'mask_' + file_name
    }
}

VEVIO_FRAMES_1_100_1 = copy.deepcopy(VEVIO_FRAMES)
VEVIO_FRAMES_1_100_1['Dirs']['logs'] = 'logs' + sep + 'VEVIO_FRAMES' + sep + 'UNET_1_100_1'
VEVIO_FRAMES_1_100_1['Funcs']['dparm'] = dparm_1_100_1

VEVIO_FRAMES_1_1 = copy.deepcopy(VEVIO_FRAMES)
VEVIO_FRAMES_1_1['Dirs']['logs'] = 'logs' + sep + 'VEVIO_FRAMES' + sep + 'UNET_1_1'
VEVIO_FRAMES_1_1['Funcs']['dparm'] = dparm_1_1

VEVIO_FRAMES_WEIGHTED = copy.deepcopy(VEVIO_FRAMES)
VEVIO_FRAMES_WEIGHTED['Dirs']['logs'] = 'logs' + sep + 'VEVIO_FRAMES' + sep + 'UNET_WEIGHTED'
VEVIO_FRAMES_WEIGHTED['Funcs']['dparm'] = d_parm_weighted
# -------------------------------------------------------------------------------------------------------------

max_issues_repo_name: LucasLee-ff/deepdyn

max_forks_repo_name: aksish/ature


text:
import matplotlib
import __main__ as main
if hasattr(main, '__file__'):
    matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
import astropy.units as u
from scipy.interpolate import interp1d, interp2d, RectBivariateSpline, griddata

matplotlib.rcParams.update({'font.size': 14})

tmp = np.load('47UMac_comps.npz')
allhs = tmp['allhs']
allcs = tmp['allcs']
ts = np.linspace(0,2391,100)
allcs = allcs[:100]

wfirstcontr = np.genfromtxt('WFIRST_pred_imaging.txt')
contr = wfirstcontr[:,1]
angsep = wfirstcontr[:,0] #l/D
angsep = (angsep * (575.0*u.nm)/(2.37*u.m)*u.rad).decompose().to(u.mas).value #mas
wfirstc = interp1d(angsep,contr,bounds_error = False, fill_value = 'extrapolate')

minangsep=150
maxangsep=450
WAbins0 = np.arange(minangsep,maxangsep+1,1)
WAbins = np.hstack((0, WAbins0, np.inf))
dMagbins0 = np.arange(0,26.1,0.1)
dMagbins = np.hstack((dMagbins0,np.inf))

WAc,dMagc = np.meshgrid(WAbins0[:-1]+np.diff(WAbins0)/2.0,dMagbins0[:-1]+np.diff(dMagbins0)/2.0)
WAc = WAc.T
dMagc = dMagc.T

fig,ax = plt.subplots(1,2,figsize=(10,4))

cs = ax[0].contourf(WAc,dMagc,np.log10(allhs[0]),np.arange(-6,-2.7,0.2))
ax[0].set_xlabel('Separation (mas)')
ax[0].set_ylabel('$\Delta$mag')
cbar = fig.colorbar(cs,ax=ax[0])
cbar.ax.set_title('log$_{10}(f_{\Delta\\mathrm{mag},\\alpha})$',pad=15)

ax[0].set_xlim([150,350])
ax[0].set_ylim([19,26])

ax[0].plot(angsep,-2.5*np.log10(contr),'r')

csplt = ax[1].plot(ts/ts.max(),allcs)
ax[1].set_xlabel('Fraction of Orbital Period')
ax[1].set_ylabel('Completeness')
ax[1].set_ylim([0.61,0.68])
ax[1].set_xlim([0,1])

plt.subplots_adjust(bottom=0.15, right=0.985, top=0.9,left=0.075)

def drawFrame(j):
    print(j)

    #for c in cs.collections:
    #    c.remove()
    #csplt[0].remove()

    ax[0].cla()
    ax[1].cla()
    cs = ax[0].contourf(WAc,dMagc,np.log10(allhs[j]),np.arange(-6,-2.7,0.2))
    ax[0].set_xlabel('Separation (mas)')
    ax[0].set_ylabel('$\Delta$mag')
    ax[0].set_xlim([150,350])
    ax[0].set_ylim([19,26])
    ax[0].plot(angsep,-2.5*np.log10(contr),'r')

    csplt = ax[1].plot(ts[:j+1]/ts.max(),allcs[:j+1])
    ax[1].set_xlabel('Fraction of Orbital Period')
    ax[1].set_ylabel('Completeness')
    ax[1].set_ylim([0.61,0.68])
    ax[1].set_xlim([0,1])

    return cs,csplt

fps = 15

if __name__ == '__main__':
    anim = animation.FuncAnimation(fig, drawFrame, frames=100, interval=1./fps*1000, blit=False)
    Writer = animation.writers['ffmpeg']
    writer = Writer(fps=fps)
    anim.save('47UMa_c_comps_anim.mp4', writer=writer)

reason: import numpy,from scipy,import astropy