text:
# pylint: disable=I0011,W0613,W0201,W0212,E1101,E1103

from __future__ import absolute_import, division, print_function

from itertools import product

import numpy as np

from ..component_link import ComponentLink
from ..data import ComponentID, DerivedComponent, Data, Component
from ..coordinates import Coordinates
from ..data_collection import DataCollection
from ..link_manager import (LinkManager, accessible_links, discover_links,
                            find_dependents, is_convertible_to_single_pixel_cid,
                            equivalent_pixel_cids)
from ..link_helpers import LinkSame

comp = Component(data=np.array([1, 2, 3]))

def example_components(self, add_derived=True):
    """ Link Topology

           --- c1---c3--\
    data --|             --c5,c6   (c7,c8 disconnected)
           --- c2---c4--/
    """
    self.data = Data()
    c1 = ComponentID('c1')
    c2 = ComponentID('c2')
    c3 = ComponentID('c3')
    c4 = ComponentID('c4')
    c5 = ComponentID('c5')
    c6 = ComponentID('c6')
    c7 = ComponentID('c7')
    c8 = ComponentID('c8')

    dummy_using = lambda x, y: (x, y)
    self.cs = [c1, c2, c3, c4, c5, c6, c7, c8]
    self.links = [ComponentLink([c1], c3, lambda x:x),
                  ComponentLink([c2], c4, lambda x:x),
                  ComponentLink([c3], c1, lambda x:x),
                  ComponentLink([c4], c2, lambda x:x),
                  ComponentLink([c3, c4], c5, dummy_using),
                  ComponentLink([c3, c4], c6, dummy_using)]

    self.data.add_component(comp, c1)
    self.data.add_component(comp, c2)
    if add_derived:
        for i in [0, 1, 4, 5]:
            dc = DerivedComponent(self.data, self.links[i])
            self.data.add_component(dc, dc.link.get_to_id())

    self.primary = [c1, c2]
    self.direct = [c3, c4]
    self.derived = []
    self.externally_derived = [c5, c6]
    self.inaccessible = [c7, c8]

class TestAccessibleLinks(object):

    def setup_method(self, method):
        self.cs = [ComponentID("%i" % i) for i in range(10)]

    def test_returned_if_available(self):
        cids = self.cs[0:5]
        links = [ComponentLink([self.cs[0]], self.cs[1])]
        assert links[0] in accessible_links(cids, links)

    def test_returned_if_reachable(self):
        cids = self.cs[0:5]
        links = [ComponentLink([self.cs[0]], self.cs[6])]
        assert links[0] in accessible_links(cids, links)

    def test_not_returned_if_not_reachable(self):
        cids = self.cs[0:5]
        links = [ComponentLink([self.cs[6]], self.cs[7])]
        assert not links[0] in accessible_links(cids, links)

class TestDiscoverLinks(object):

    def setup_method(self, method):
        example_components(self)

    def test_correct_discover(self):
        """discover_links finds the correct links"""
        links = discover_links(self.data, self.links)

        for i in self.inaccessible:
            assert i not in links

        for d in self.direct:
            assert d in links

        for d in self.derived:
            assert d in links

        for p in self.primary:
            assert p not in links

    def test_links_point_to_proper_ids(self):
        """ Dictionary values are ComponentLinks which
        point to the keys """
        links = discover_links(self.data, self.links)
        for cid in links:
            assert cid == links[cid].get_to_id()

    def test_shortest_path(self):
        """ Shortcircuit c5 to c1, yielding 2 ways to get to c5.
        Ensure that the shortest path is chosen """
        self.links.append(ComponentLink([self.cs[0]], self.cs[4]))
        links = discover_links(self.data, self.links)

        assert links[self.cs[4]] is self.links[-1]

class TestFindDependents(object):

    def setup_method(self, method):
        example_components(self)

    def test_propagated(self):
        to_remove = self.links[0]
        result = find_dependents(self.data, to_remove)
        expected = set([self.cs[2], self.cs[4], self.cs[5]])
        assert expected == result

    def test_basic(self):
        to_remove = self.links[4]
        result = find_dependents(self.data, to_remove)
        expected = set([self.cs[4]])
        assert expected == result

class TestLinkManager(object):

    def test_add_links(self):
        id1 = ComponentID('id1')
        id2 = ComponentID('id2')
        id3 = ComponentID('id3')
        lm = LinkManager()
        using = lambda x, y: 0
        link = ComponentLink([id1, id2], id3, using)
        lm.add_link(link)
        links = lm.links
        assert links == [link]

    def test_remove_link(self):
        id1 = ComponentID('id1')
        id2 = ComponentID('id2')
        id3 = ComponentID('id3')
        lm = LinkManager()
        using = lambda x, y: 0
        link = ComponentLink([id1, id2], id3, using)
        lm.add_link(link)
        lm.remove_link(link)
        links = lm.links
        assert links == []

    def test_setup(self):
        example_components(self, add_derived=False)
        expected = set()
        assert set(self.data.derived_components) == expected

    def test_update_externally_derivable_components_adds_correctly(self):
        example_components(self, add_derived=False)
        lm = LinkManager()
        list(map(lm.add_link, self.links))

        lm.update_externally_derivable_components(self.data)
        derived = set(self.data.externally_derivable_components)
        expected = set(self.externally_derived + self.direct)
        assert derived == expected

    def test_update_externally_derivable_components_removes_correctly(self):
        # add all but last link to manager
        example_components(self, add_derived=False)
        lm = LinkManager()
        list(map(lm.add_link, self.links[:-1]))
        lm.update_externally_derivable_components(self.data)

        # manually add last link as derived component
        dc = DerivedComponent(self.data, self.links[-1])
        self.data._externally_derivable_components.update({dc.link.get_to_id(): dc})
        removed = set([dc.link.get_to_id()])
        assert dc.link.get_to_id() in self.data.externally_derivable_components

        # this link should be removed upon update_components
        lm.update_externally_derivable_components(self.data)
        derived = set(self.data.externally_derivable_components)
        expected = set(self.direct + self.externally_derived) - removed
        assert derived == expected

    def test_derived_links_correctwith_mergers(self):
        """When the link manager merges components, links that depend on the
        merged components remain functional"""

        d1 = Data(x=[[1, 2], [3, 4]])
        d2 = Data(u=[[5, 6], [7, 8]])

        dc = DataCollection([d1, d2])

        # link world coordinates...
        dc.add_link(LinkSame(
            d1.world_component_ids[0], d2.world_component_ids[0]))
        dc.add_link(LinkSame(
            d1.world_component_ids[1], d2.world_component_ids[1]))

        # and then retrieve pixel coordinates
        np.testing.assert_array_equal(
            d2[d1.pixel_component_ids[0]], [[0, 0], [1, 1]])
        np.testing.assert_array_equal(
            d1[d2.pixel_component_ids[1]], [[0, 1], [0, 1]])

    def test_binary_links_correct_with_mergers(self):
        """Regression test. BinaryComponentLinks should work after mergers"""

        d1 = Data(x=[1, 2, 3], y=[2, 3, 4])
        d2 = Data(u=[2, 3, 4], v=[3, 4, 5])

        z = d1.id['x'] + d1.id['y']
        d1.add_component_link(z, 'z')

        dc = DataCollection([d1, d2])
        dc.add_link(LinkSame(d2.id['u'], d1.id['x']))

        np.testing.assert_array_equal(d1['z'], [3, 5, 7])

    def test_complex_links_correct_with_mergers(self):
        """Regression test. multi-level links should work after mergers"""

        d1 = Data(x=[1, 2, 3], y=[2, 3, 4])
        d2 = Data(u=[2, 3, 4], v=[3, 4, 5])
        x = d1.id['x']

        z = d1.id['x'] + d1.id['y'] + 5
        d1.add_component_link(z, 'z')

        dc = DataCollection([d1, d2])
        dc.add_link(LinkSame(d2.id['u'], d1.id['x']))

        # NOTE: the behavior tested here is not desirable anymore, so the
        #       following assert is no longer True
        # assert x not in d1.components

        np.testing.assert_array_equal(d1['z'], [8, 10, 12])

    def test_remove_data_removes_links(self):

        d1 = Data(x=[[1, 2], [3, 4]], label='image')
        d2 = Data(a=[1, 2, 3], b=[4, 5, 6], label='catalog')

        dc = DataCollection([d1, d2])

        assert len(dc.links) == 6

        dc.add_link(LinkSame(d1.id['x'], d2.id['a']))

        assert len(dc.links) == 7

        # Removing dataset should remove related links
        dc.remove(d1)
        assert len(dc.links) == 2

    def test_remove_component_removes_links(self):

        d1 = Data(x=[[1, 2], [3, 4]], label='image')
        d2 = Data(a=[1, 2, 3], b=[4, 5, 6], label='catalog')

        dc = DataCollection([d1, d2])

        assert len(dc.links) == 6

        assert len(d1.components) == 5
        assert len(d2.components) == 4

        assert len(d1.externally_derivable_components) == 0
        assert len(d2.externally_derivable_components) == 0

        dc.add_link(LinkSame(d1.id['x'], d2.id['a']))

        assert len(dc.links) == 7

        assert len(d1.components) == 5
        assert len(d2.components) == 4

        assert len(d1.externally_derivable_components) == 1
        assert len(d2.externally_derivable_components) == 1

        assert d1.id['x'] in d2.externally_derivable_components
        assert d2.id['a'] in d1.externally_derivable_components

        # Removing component a from d2 should remove related links and
        # derived components.
        d2.remove_component(d2.id['a'])

        assert len(dc.links) == 6

        assert len(d1.components) == 5
        assert len(d2.components) == 3

        assert len(d1.externally_derivable_components) == 0
        assert len(d2.externally_derivable_components) == 0

def test_is_convertible_to_single_pixel_cid():

    # This tests the function is_convertible_to_single_pixel_cid, which gives
    # for a given dataset the pixel component ID that can be uniquely
    # transformed into the requested component ID.

    # Set up a coordinate object which has an independent first axis and
    # has the second and third axes depend on each other. The transformation
    # itself is irrelevant since for this function to work we only care about
    # whether or not an axis is independent.

    class CustomCoordinates(Coordinates):
        def dependent_axes(self, axis):
            if axis == 0:
                return (0,)
            else:
                return (1, 2)

    data1 = Data()
    data1.coords = CustomCoordinates()
    data1['x'] = np.ones((4, 3, 4))
    px1, py1, pz1 = data1.pixel_component_ids
    wx1, wy1, wz1 = data1.world_component_ids
    data1['a'] = px1 * 2
    data1['b'] = wx1 * 2
    data1['c'] = wy1 * 2
    data1['d'] = wx1 * 2 + px1
    data1['e'] = wx1 * 2 + wy1

    # Pixel component IDs should just be returned directly
    for cid in data1.pixel_component_ids:
        assert is_convertible_to_single_pixel_cid(data1, cid) is cid

    # Only the first world component should return a valid pixel component
    # ID since the two other world components are interlinked
    assert is_convertible_to_single_pixel_cid(data1, wx1) is px1
    assert is_convertible_to_single_pixel_cid(data1, wy1) is None
    assert is_convertible_to_single_pixel_cid(data1, wz1) is None

    # a and b are ultimately linked to the first pixel coordinate, whereas c
    # depends on the second world coordinate which is interlinked with the third
    # Finally, d is ok because it really only depends on px1
    assert is_convertible_to_single_pixel_cid(data1, data1.id['a']) is px1
    assert is_convertible_to_single_pixel_cid(data1, data1.id['b']) is px1
    assert is_convertible_to_single_pixel_cid(data1, data1.id['c']) is None
    assert is_convertible_to_single_pixel_cid(data1, data1.id['d']) is px1
    assert is_convertible_to_single_pixel_cid(data1, data1.id['e']) is None

    # We now create a second dataset and set up links
    data2 = Data(y=np.ones((4, 5, 6, 7)), z=np.zeros((4, 5, 6, 7)))
    dc = DataCollection([data1, data2])
    dc.add_link(ComponentLink([data1.id['a'], px1], data2.id['y'], using=lambda x: 2 * x))
    dc.add_link(ComponentLink([wy1], data2.id['z'], using=lambda x: 2 * x))

    assert is_convertible_to_single_pixel_cid(data1, data2.id['y']) is px1
    assert is_convertible_to_single_pixel_cid(data1, data2.id['z']) is None

def test_equivalent_pixel_cids():

    # This tests the equivalent_pixel_cids function which checks whether all
    # pixel component IDs in a dataset are equivalent to pixel component IDs
    # in another dataset.

    data1 = Data(x=np.ones((2, 4, 3)))
    data2 = Data(y=np.ones((4, 3, 2)))
    data3 = Data(z=np.ones((2, 3)))

    dc = DataCollection([data1, data2, data3])

    # Checking data against itself should give the normal axis order

    assert equivalent_pixel_cids(data1, data1) == [0, 1, 2]
    assert equivalent_pixel_cids(data3, data3) == [0, 1]

    # Checking data against unlinked data should give None

    for d1, d2 in product(dc, dc):
        if d1 is not d2:
            assert equivalent_pixel_cids(d1, d2) is None

    # Add one set of links, which shouldn't change anything

    dc.add_link(LinkSame(data1.pixel_component_ids[0], data2.pixel_component_ids[2]))
    dc.add_link(LinkSame(data1.pixel_component_ids[0], data3.pixel_component_ids[0]))

    for d1, d2 in product(dc, dc):
        if d1 is not d2:
            assert equivalent_pixel_cids(d1, d2) is None

    # Add links between a second set of axes

    dc.add_link(LinkSame(data1.pixel_component_ids[2], data2.pixel_component_ids[1]))
    dc.add_link(LinkSame(data1.pixel_component_ids[2], data3.pixel_component_ids[1]))

    # At this point, data3 has all its axes contained in data1 and data2

    assert equivalent_pixel_cids(data1, data2) is None
    assert equivalent_pixel_cids(data2, data1) is None
    assert equivalent_pixel_cids(data1, data3) == [0, 2]
    assert equivalent_pixel_cids(data3, data1) is None
    assert equivalent_pixel_cids(data2, data3) == [2, 1]
    assert equivalent_pixel_cids(data3, data2) is None

    # Finally we can link the third set of axes

    dc.add_link(LinkSame(data1.pixel_component_ids[1], data2.pixel_component_ids[0]))

    # At this point, data1 and data2 should now be linked. Note that in cases
    # where the target has more dimensions than the reference, we always get None

    assert equivalent_pixel_cids(data1, data2) == [1, 2, 0]
    assert equivalent_pixel_cids(data2, data1) == [2, 0, 1]
    assert equivalent_pixel_cids(data1, data3) == [0, 2]
    assert equivalent_pixel_cids(data3, data1) is None
    assert equivalent_pixel_cids(data2, data3) == [2, 1]
    assert equivalent_pixel_cids(data3, data2) is None

max_stars_repo_name: tiagopereira/glue

max_issues_repo_name: tiagopereira/glue

max_forks_repo_name: tiagopereira/glue


text:
# import packages used
import numpy as np
import tools_Exercise_1_6 as tools
import scipy.optimize as optimize
import scipy.interpolate as interpolate
import time as time

def setup():
    class par: pass
    # Model parameters
    par.beta = 0.999
    par.B = 0.33
    par.upsillon=20
    par.Lt = 1
    par.W = 20
    par.G= 0.99
    par.chi = 30
    par.xi1 = 0
    par.xi2 = 0
    par.D = 0.005
    par.d = 0.005
    par.varphi = 0
    par.kappa1 = 1
    par.kappa2 = 8
    par.Upsillon = 0.51*par.upsillon
    par.Z = 75000
    par.gamma1 = 0.052
    par.tests = 0.01
    par.varsigma = 13
    par.varrho = 2
    par.t = 1.8
    par.phi1 = 0.2*0.37
    par.phi2 = 0.2*0.33
    par.phi3 = 0.2*0.3
    par.sigma = 0.001
    par.varrho = 0.4
    par.alpha=0.3
    par.rho=5
    par.g=40
    par.mu = 2
    par.H = 4.7
    par.d = 0.00025
    par.nd = 0.001
    par.cap = 0.02
    par.reinf = 0.000288
    # Shock parameters
    par.num_M = 7
    par.M_max = 0.3
    par.num_shocks = 8

    # Convergens settings
    par.max_iter = 10000   # maximum number of iterations
    par.tol = 10e-2

    # Simulation parameters
    par.simN = 700
    par.I_ini = 0.01
    par.Q_ini = 0.001
    par.R_ini = 0.001
    par.D_ini = 0.001
    par.lw_ini = 1

    # Setup grid
    setup_grids(par)
    return par

def setup_grids(par):

    #Grid of disease parameters
    par.grid_I = tools.nonlinspace(1.0e-10,par.M_max,par.num_M,1.5) # non-linear spaced points: like np.linspace with unequal spacing
    par.grid_Q = tools.nonlinspace(1.0e-10,par.M_max,par.num_M,1.5) # non-linear spaced points: like np.linspace with unequal spacing
    par.grid_R = tools.nonlinspace(1.0e-10,0.7,par.num_M,1.5) # non-linear spaced points: like np.linspace with unequal spacing
    par.grid_D = tools.nonlinspace(1.0e-10,0.05,par.num_M,1.5)
    par.grid_lw = tools.nonlinspace(1.0e-10,1,100,1)
    #Gauss-Hermite
   # x,w = tools.gauss_hermite(par.num_shocks)
   # par.eps = np.exp(par.sigma*np.sqrt(2)*x)
   #  par.eps_w = w/np.sqrt(np.pi)
    return par
def solve_cons_inf(par):

    # Initalize
    class sol: pass
    sol.V = np.ones([par.num_M, par.num_M, par.num_M, par.num_M])*1e-5
    sol.lw = np.zeros([par.num_M, par.num_M, par.num_M, par.num_M])

    sol.it = 0   #Number of iteration
    sol.delta = 1000.0 #Different between V+ and V
    sol.S=[]
    sol.lo=[]
    sol.s=[]
    sol.wi=[]
    sol.Y=[]
    sol.i=[]
    sol.l=[]
    sol.gamma2=[]
    sol.gamma3=[]
    sol.I_plus=[]
    sol.Q_plus=[]
    sol.R_plus=[]
    sol.D_plus=[]
    sol.p=[]
    sol.pi=[]
    prcompo = 0
        #precomp
    for  I in (par.grid_I):
        for Q in (par.grid_Q):
            for R in (par.grid_R):
                for D in (par.grid_D):
                    for lw in (par.grid_lw):
                        if lw+Q+D > 1:
                            break
                        S=(1-I-Q-R)
                        lo=(1 - lw - Q - par.D*R)
                        s=min(max((lw-(1-D)*R)*(1-I/(S+I)),0),1)
                        wi=min(max((lw-(1-D)*R)*(I/(S+I)),0),1)
                        Y=max(par.H*np.log(par.upsillon*lw+par.Upsillon*lo)-(par.chi*I)**2 - par.varphi*R, 1.0e-8)
                        #w=(lw+Q+lo*par.G)*par.W
                        w=par.W
                        #print(Y)
                        l=((par.Z*par.phi2*I*max(1-R-Q,0.1)/(par.alpha*par.varsigma))**(1/(par.alpha-1)))*Y
                        if l<0:
                            p=0
                            l=0
                        elif l>1:
                            l=1
                            p=((1-par.alpha)*par.varsigma)*Y**-par.alpha
                        else:
                            p=((1-par.alpha)*par.varsigma*l**(par.alpha) * Y**(-par.alpha))
                        if p*Y>w+par.g:
                            p=(w+par.g)/Y
                        #print(p)
                        #print(l)
                        gamma2=np.array(par.sigma + (par.t*par.tests/((1 + I*par.rho)**par.mu)))
                        gamma3=np.array(par.gamma1 * (1+ par.kappa1/(1+Q**(1/par.kappa2))))
                        if Q+I > par.cap:
                            sol.I_plus.append(max(min((1-par.gamma1-gamma2)*I + par.phi1*s*wi + par.phi2*S*I*l*l + par.phi3*S*I-(I/(I+Q))*par.nd*(I+Q-par.cap) -(I/(I+Q))*par.cap*par.d,1),1.0e-9))
                            sol.Q_plus.append(max(min((1- gamma3)*Q + gamma2*I-(Q/(I+Q))*par.nd*(I+Q-par.cap)-(Q/(I+Q))*par.cap*par.d,1),1.0e-9))
                            sol.D_plus.append(D+par.nd*(I+Q-par.cap)+par.cap*par.d)
                        else:
                            sol.I_plus.append(max(min((1-par.gamma1-gamma2)*I + par.phi1*s*wi + par.phi2*S*I*l*l + par.phi3*S*I-par.d*I,1),1.0e-9))
                            sol.Q_plus.append(max(min((1- gamma3)*Q + gamma2*I-par.d*Q,1),1.0e-9))
                            sol.D_plus.append(D+par.d*(I+Q))
                        sol.R_plus.append(max(min(R*(1-par.reinf) + par.gamma1*I + gamma3*Q,1),1.0e-9))
                        sol.pi.append(Y*p - (lw+Q)*par.W - lo*par.G*par.W - par.xi1*I**2 - par.xi2*par.d*R)
                        #print(Y*p - (lw+Q)*par.W - lo*par.G*par.W - (par.xi1*I)**2 - par.xi2*par.d*R)
                        #print(par.W+par.g-Y*p)
                        prcompo +=1

    #points=np.meshgrid(par.grid_I, par.grid_Q, par.grid_R, copy=False, indexing='xy')
    points = (par.grid_I, par.grid_Q, par.grid_R, par.grid_D)
    #print(np.shape(points))
    #print(max(sol.I_plus))
    #print(min(sol.I_plus))
    #print(max(sol.Q_plus))
    #print(min(sol.Q_plus))
    #print(max(sol.R_plus))
    #print(min(sol.R_plus))
    point = np.transpose(np.array([sol.I_plus, sol.Q_plus, sol.R_plus, sol.D_plus]))
    while (sol.delta >= par.tol and sol.it < par.max_iter):
        V_next = sol.V.copy()
        V_plus = interpolate.interpn(points, V_next, point, method='linear', bounds_error=False, fill_value=None)
        print(np.sum(V_plus))
        ind = -1
        # find V
        Ih = -1
        Qh = -1
        Rh = -1
        Dh = -1

        for I in (par.grid_I):
            Ih +=1
            for Q in (par.grid_Q):
                Qh +=1
                for R in (par.grid_R):
                    Rh +=1
                    for D in (par.grid_D):
                        Dh +=1
                        for lw in (par.grid_lw):
                            if lw+Q+D > 1:
                                break
                            ind += 1
                            V_guess =sol.pi[ind] + par.beta*V_plus[ind]
                            if V_guess > sol.V[Ih, Qh, Rh, Dh]:
                                sol.V[Ih, Qh, Rh, Dh]=V_guess
                                sol.lw[Ih, Qh, Rh, Dh]=lw
                    Dh=-1
                Rh=-1
            Qh=-1

        # opdate delta and it

        sol.it += 1
        c_new = np.ravel(sol.V)
        c_old = np.ravel(V_next)
        #sol.delta = max(abs(sol.V - V_next))
        sol.delta = max(abs(c_new - c_old))
        print(sol.delta)
    return(sol)
def simu(par, sol):
    class simu: pass

    simu.S=np.zeros([par.simN])
    simu.lo=np.zeros([par.simN])
    simu.s=np.zeros([par.simN])
    simu.wi=np.zeros([par.simN])
    simu.Y=np.zeros([par.simN])
    simu.l=np.zeros([par.simN])
    simu.p=np.zeros([par.simN])
    simu.gamma2=np.zeros([par.simN])
    simu.gamma3=np.zeros([par.simN])
    simu.pi=np.zeros([par.simN])
    simu.util=np.zeros([par.simN])
    simu.c=np.zeros([par.simN])
    simu.I=np.zeros([par.simN+1])
    simu.Q=np.zeros([par.simN+1])
    simu.R=np.zeros([par.simN+1])
    simu.D=np.zeros([par.simN+1])
    simu.w=np.zeros([par.simN])
    simu.Pos=np.zeros([par.simN])
    simu.I[0]=(par.I_ini)
    simu.Q[0]=(par.Q_ini)
    simu.R[0]=(par.R_ini)
    simu.D[0]=(par.D_ini)
    simu.lw =np.zeros([par.simN])
    simu.lw[0] = 1
    ite=0
    points = (par.grid_I, par.grid_Q, par.grid_R, par.grid_D)
    while ite < par.simN:
    #Start of simulation.
        simu.lw[ite]=min(interpolate.interpn(points, sol.lw, ([simu.I[ite], simu.Q[ite], simu.R[ite], simu.D[ite]]), method='linear', bounds_error=False, fill_value=None), 1-simu.Q[ite]-simu.R[ite]*par.D)

        simu.lw[ite]=max(min(simu.lw[ite], 1-simu.Q[ite]-simu.D[ite]), 0)

        if ite == 0:

            simu.lw[ite]=1

        simu.lw[ite]=max(simu.lw[ite],0)
        print(simu.lw[ite])

        simu.S[ite]=(1-simu.I[ite]-simu.Q[ite]-simu.R[ite])

        simu.lo[ite]=(1 - simu.lw[ite] - simu.Q[ite] - par.D*simu.R[ite])

        simu.s[ite]=(max((simu.lw[ite]-(1-par.D)*simu.R[ite])*(1-simu.I[ite]/(simu.S[ite]+simu.I[ite])),1.0e-9))

        simu.wi[ite]=(max((simu.lw[ite]-(1-par.D)*simu.R[ite])*(simu.I[ite]/(simu.S[ite]+simu.I[ite])),1.0e-9))

        simu.Y[ite]=(max(par.H*np.log(par.upsillon*simu.lw[ite]+par.Upsillon*simu.lo[ite])-(par.chi*simu.I[ite])**2 - par.varphi*simu.R[ite], 1.0e-9))

        simu.l[ite]=(par.Z*par.phi2*simu.I[ite]*(1-simu.R[ite]-simu.Q[ite])/(par.alpha*par.varsigma))**(1/(par.alpha-1))*simu.Y[ite]

        simu.w[ite]=(simu.lo[ite]*par.G+simu.lw[ite]+simu.Q[ite])*par.W

        if  simu.l[ite] < 0:
            simu.l[ite]=0
            simu.p[ite] = 0

        elif simu.l[ite]>1:
            simu.l[ite]=1
            simu.p[ite]=((1-par.alpha)*par.varsigma)/(simu.Y[ite]**par.alpha)
        else:
            simu.p[ite]=(1-par.alpha)*par.varsigma*simu.l[ite]**(par.alpha) * simu.Y[ite]**(-par.alpha)

        if simu.p[ite]*simu.Y[ite]>simu.w[ite]+par.g:
            simu.p[ite]=(simu.w[ite]+par.g)/simu.Y[ite]

        simu.gamma2[ite]=(np.array(par.sigma + (par.t*par.tests/(1 + simu.I[ite]*par.rho)**par.mu)))

        simu.gamma3[ite]=(np.array(par.gamma1 * (1+ par.kappa1/(1+simu.Q[ite]**(1/par.kappa2)))))

        simu.pi[ite]=(simu.Y[ite]*simu.p[ite] -(simu.lw[ite]+simu.Q[ite])*par.W - simu.lo[ite]*par.G*par.W - (par.xi1*simu.I[ite])**2 - par.xi2*par.d*simu.R[ite])

        simu.util[ite]=(par.varsigma*simu.l[ite]**par.alpha*simu.Y[ite]**(1-par.alpha)+simu.w[ite]+par.g-simu.p[ite]*simu.Y[ite]-par.Z*par.phi2*simu.I[ite]*simu.l[ite]*(1-simu.R[ite]-simu.Q[ite])- par.Z*par.phi3*simu.I[ite]*(1-simu.R[ite]-simu.Q[ite]))

        simu.c[ite]=simu.w[ite]+par.g-simu.p[ite]*simu.Y[ite]

        simu.Pos[ite]=simu.gamma2[ite]*simu.I[ite]/(par.tests)*100

        if simu.Q[ite]+simu.I[ite] > par.cap:
            simu.I[ite+1]=(max(min((1-par.gamma1-simu.gamma2[ite])*simu.I[ite] + par.phi1*simu.s[ite]*simu.wi[ite] + par.phi2*simu.S[ite]*simu.I[ite]*simu.l[ite]*simu.l[ite] + par.phi3*simu.S[ite]*simu.I[ite] -(simu.I[ite]/(simu.I[ite]+simu.Q[ite]))*par.cap*par.d -(simu.I[ite]/(simu.I[ite]+simu.Q[ite]))*par.nd*(simu.Q[ite]+simu.I[ite]-par.cap),1),1.0e-9))
            simu.Q[ite+1]=(max(min((1- simu.gamma3[ite])*simu.Q[ite] + simu.gamma2[ite]*simu.I[ite]-(simu.Q[ite]/(simu.I[ite]+simu.Q[ite]))*par.cap*par.d - (simu.Q[ite]/(simu.I[ite]+simu.Q[ite]))*par.nd*(simu.Q[ite]+simu.I[ite]-par.cap),1),1.0e-9))
            simu.D[ite+1]=simu.D[ite] +  par.cap*par.d +par.nd*(simu.I[ite]+simu.Q[ite]-par.cap)
        else:
            simu.I[ite+1]=(max(min((1-par.gamma1-simu.gamma2[ite])*simu.I[ite] + par.phi1*simu.s[ite]*simu.wi[ite] + par.phi2*simu.S[ite]*simu.I[ite]*simu.l[ite]*simu.l[ite] + par.phi3*simu.S[ite]*simu.I[ite]- par.d*simu.I[ite],1),1.0e-9))
            simu.Q[ite+1]=(max(min((1- simu.gamma3[ite])*simu.Q[ite] + simu.gamma2[ite]*simu.I[ite] - par.d*simu.Q[ite],1),1.0e-9))
            simu.D[ite+1]=simu.D[ite]+par.d*(simu.I[ite]+simu.Q[ite])
        simu.R[ite+1]=(max(min(simu.R[ite]*(1-par.reinf) + par.gamma1*simu.I[ite] + simu.gamma3[ite]*simu.Q[ite],1),1.0e-9))
        ite+=1

    simu.grid = np.linspace(0,ite,ite)
    simu.I = simu.I[0:ite]
    simu.Q = simu.Q[0:ite]
    simu.R = simu.R[0:ite]
    simu.D = simu.D[0:ite]
    simu.GDP = simu.p*simu.Y
    return(simu)

ext: py

max_stars_repo_name: Madskmeier/MastersThesis

max_issues_repo_name: Madskmeier/MastersThesis

max_forks_repo_name: Madskmeier/MastersThesis

reason: import numpy,import scipy


text:
"""
The :mod:`constants` module provides a number of constants used by the :class:`SLRealizer` worker class and its descendants.
"""
from astropy.coordinates import SkyCoord

def return_zeropoint():
    """
    returns zeropoint for the SDSS magnitude system
    """
    return 22.5

def return_obs_RA_DEC():
    """
    Return OM10's lensed system's position
    """
    return SkyCoord('03h 32m 30s', '10d 00m 24s')

def get_filter_AB_offset():
    return 0.0

def get_flux_err():
    return 0.01 # unit : percentage

def get_flux_err_std():
    return 0.005 # unit : percentage

def get_first_moment_err():
    return 0.01 # unit : percentage

def get_first_moment_err_std():
    return 0.005

def get_second_moment_err():
    return 0.15 # changed from 0.05

def get_second_moment_err_std():
    return 0.03 # changed from 0.01

def get_x_min():
    return -5.0

def get_y_min():
    return -5.0

def get_x_max():
    return 5.0

def get_y_max():
    return 5.0

def get_ellipticity_cut():
    return 0.4

def get_SDSS_pixel_arcsec_conversion():
    return 0.396

max_stars_repo_name: LSSTDESC/SLRealizer

max_issues_repo_name: LSSTDESC/SLRealizer

max_forks_repo_name: LSSTDESC/SLRealizer

reason: from astropy


text:
"""
Interpolates MCD12C1 data to lidar trajectory
"""

import os, sys, subprocess
from   datetime               import date, datetime, timedelta
from   dateutil.parser        import parse as isoparser
from   dateutil.relativedelta import relativedelta
from   pyhdf.SD               import SD, HDF4Error
import numpy                  as     np
from   glob                   import glob
from   scipy.interpolate      import RegularGridInterpolator
from   netCDF4                import Dataset
from   py_leo_vlidort.copyvar                import _copyVar

nkernels = 3
SDS = { 'Majority_Land_Cover_Type_1'                 : 'IGBPlc'}

LANDCOVER = {"water"                              : 0,
             "evergreen needleleaf forest"        : 1,
             "evergreen broadleaf forest"         : 2,
             "deciduous needleleaf forest"        : 3,
             "deciduous broadleaf forest"         : 4,
             "mixed forests"                      : 5,
             "closed shrubland"                   : 6,
             "open shrublands"                    : 7,
             "woody savannas"                     : 8,
             "savannas"                           : 9,
             "grasslands"                         : 10,
             "permanent wetlands"                 : 11,
             "croplands"                          : 12,
             "urban and built-up"                 : 13,
             "cropland/natural vegetation mosaic" : 14,
             "snow and ice"                       : 15,
             "barren or sparsely vegetated"       : 16}

#from Maignan et al. RSE, 2009 Table 1
COEF = {0 : -999,
        1 : 4.98,
        2 : 7.36,
        3 : 5.98,
        4 : 6.87,
        5 : 4.95,
        6 : 4.71,
        7 : 5.99,
        8 : 6.16,
        9 : 6.66,
        10: 6.11,
        11: -999,
        12: 7.58,
        13: -999,
        14: 6.86,
        15: 7.07,
        16: 7.29}

class MCD12C(object):
    def __init__(self,datadir,verbose=False):
        self.verbose = verbose
        self.dlon = 0.05
        self.dlat = 0.05

        self.nEW = 360./self.dlon
        self.nNS = 180./self.dlat

        self.lon = np.linspace(-180 + 0.5*self.dlon,180 - 0.5*self.dlon,self.nEW)
        self.lat = np.linspace(-90 + 0.5*self.dlat,90 - 0.5*self.dlat,self.nNS)

        self.fileName = datadir+'/MCD12C1.A%year001.051.hdf'
        self.SDS      = SDS

    def readFile(self,tyme):

        hfile = SD(self.fileName.replace('%year',str(tyme.year)))

        for sds in SDS:
            v = hfile.select(sds).get()
            a = hfile.select(sds).attributes()

            v = v.astype('float')
            fill_value = a['_FillValue']
            v[np.abs(v - fill_value)/fill_value < 0.01] = -999.
            self.__dict__[SDS[sds]] = np.flipud(v)

    #----
    def writenc(self,nctrj,outFile,verbose=False):
        """
        Write a NetCDF file with sampled MCD12C1 Land Cover Type
        and BPDF coefficient on lidar trajectory
        """

        # Dimensions
        # ------------------
        ntime = len(self.tyme)

        # Open NC file
        # ------------
        nc = Dataset(outFile,'w',format='NETCDF4_CLASSIC')

        # Set global attributes
        # ---------------------
        nc.title = "MCD12C1 Trajectory Sampler"
        nc.institution = 'NASA/Goddard Space Flight Center'
        nc.source = 'Global Model and Assimilation Office'
        nc.history = 'Created from MCD12C1 v051 collections by mcd12c_sampler.py'
        nc.references = 'n/a'
        nc.comment = 'This file contains IGBP Land Cover and BPDF coefficient from Table 1 in Maignan et al. 2009 sampled on a satellite trajectory'
        nc.contact = 'Patricia Castellanos <patricia.castellanos@nasa.gov>'
        nc.Conventions = 'CF'

        for lc in LANDCOVER:
            nc.__dict__['TYPE'+str(LANDCOVER[lc])] = lc

        # Create dimensions
        # -----------------
        t = nc.createDimension('time',ntime)
        s = nc.createDimension('ls',19)
        x = nc.createDimension('x',1)
        y = nc.createDimension('y',1)

        # Save lon/lat
        # --------------------------
        _copyVar(nctrj,nc,u'time', dtype='i4',zlib=False,verbose=verbose)
        _copyVar(nctrj,nc,u'x',dtype='f4',zlib=False,verbose=verbose)
        _copyVar(nctrj,nc,u'y',dtype='f4',zlib=False,verbose=verbose)
        _copyVar(nctrj,nc,u'trjLon',dtype='f4',zlib=False,verbose=verbose)
        _copyVar(nctrj,nc,u'trjLat',dtype='f4',zlib=False,verbose=verbose)
        _copyVar(nctrj,nc,u'isotime', dtype='S1',zlib=False,verbose=verbose)

        # Write land cover to file
        #---------------------------------------
        dim = ('time',)
        for sds in SDS:
          this = nc.createVariable(SDS[sds],'i2',dim)

          this.long_name = sds
          this.missing_value = -999
          this.units = 'none'
          this[:] = self.__dict__['trj'+SDS[sds]]

        this = nc.createVariable('BPDFcoef','f4',dim)
        this.long_name = 'BPDF coefficient from Maignan et al.'
        this.units = 'none'
        this[:] = self.BPDFcoef
        nc.close()

    def sample(self,inFile,outFile,Verbose=False):
        # Open lidar sampled file
        nctrj = Dataset(inFile)

        # Read in tymes
        iso = nctrj.variables['isotime'][:]
        tyme = []
        for isotime in iso:
            tyme.append(isoparser(''.join(isotime)))

        tyme = np.array(tyme)

        # Read in locations
        trjLat = nctrj.variables['trjLat'][:]
        trjLon = nctrj.variables['trjLon'][:]

        # Round dates to year
        dtyme = np.array([datetime(d.year,1,1) for d in tyme])
        utyme = np.sort(np.unique(dtyme))

        for sds in SDS:
            self.__dict__['trj'+SDS[sds]] = np.zeros(len(tyme))

        self.tyme = tyme
        self.trjLon = trjLon
        self.trjLat = trjLat
        for ut in utyme:
            if Verbose:
                print 'Working on '+ str(ut.date())
            self.readFile(ut)

            Ityme = dtyme == ut

            lat = trjLat[Ityme]
            lon = trjLon[Ityme]
            pts = []
            for LAT,LON in zip(lat,lon): pts.append([LAT,LON])

            for sds in SDS:
                interpFunc = RegularGridInterpolator((self.lat, self.lon), self.__dict__[SDS[sds]],
                                method='nearest',bounds_error=False,fill_value=None)
                self.__dict__['trj'+SDS[sds]][Ityme] = interpFunc(pts)

        self.BPDFcoef = np.array([COEF[lc] for lc in self.trjIGBPlc])
        self.writenc(nctrj,outFile,verbose=Verbose)
        nctrj.close()

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps

reason: import numpy


text:
"""Code for reformatting supervised learning targets."""
from operator import mul

import numpy as np
from six.moves import reduce
import theano.sparse
if theano.sparse.enable_sparse:
    scipy_available = True
    import scipy.sparse
else:
    scipy_available = False
from theano import tensor, config
from pylearn2.utils.exc import reraise_as

class OneHotFormatter(object):
    """
    A target formatter that transforms labels from integers in both single
    and batch mode.

    Parameters
    ----------
    max_labels : int
        The number of possible classes/labels. This means that all labels
        should be < max_labels. Example: For MNIST there are 10 numbers
        and hence max_labels = 10.
    dtype : dtype, optional
        The desired dtype for the converted one-hot vectors. Defaults to
        `config.floatX` if not given.
    """
    def __init__(self, max_labels, dtype=None):
        """
        Initializes the formatter given the number of max labels.
        """
        try:
            np.empty(max_labels)
        except (ValueError, TypeError):
            reraise_as(ValueError("%s got bad max_labels argument '%s'" %
                                  (self.__class__.__name__, str(max_labels))))
        self._max_labels = max_labels
        if dtype is None:
            self._dtype = config.floatX
        else:
            try:
                np.dtype(dtype)
            except TypeError:
                reraise_as(TypeError("%s got bad dtype identifier %s" %
                                     (self.__class__.__name__, str(dtype))))
            self._dtype = dtype

    def format(self, targets, mode='stack', sparse=False):
        """
        Formats a given array of target labels into a one-hot
        vector. If labels appear multiple times, their value
        in the one-hot vector is incremented.

        Parameters
        ----------
        targets : ndarray
            A 1D array of targets, or a batch (2D array) where
            each row is a list of targets.
        mode : string
            The way in which to convert the labels to arrays. Takes
            three different options:

              - "concatenate" : concatenates the one-hot vectors from
                multiple labels
              - "stack" : returns a matrix where each row is the
                one-hot vector of a label
              - "merge" : merges the one-hot vectors together to
                form a vector where the elements are
                the result of an indicator function
                NB: As the result of an indicator function
                the result is the same in case a label
                is duplicated in the input.
        sparse : bool
            If true then the return value is sparse matrix. Note that
            if sparse is True, then mode cannot be 'stack' because
            sparse matrices need to be 2D

        Returns
        -------
        one_hot : a NumPy array (can be 1D-3D depending on settings)
            where normally the first axis are the different batch items,
            the second axis the labels, the third axis the one_hot
            vectors. Can be dense or sparse.
        """
        if mode not in ('concatenate', 'stack', 'merge'):
            raise ValueError("%s got bad mode argument '%s'" %
                             (self.__class__.__name__, str(self._max_labels)))
        elif mode == 'stack' and sparse:
            raise ValueError("Sparse matrices need to be 2D, hence they"
                             "cannot be stacked")
        if targets.ndim > 2:
            raise ValueError("Targets needs to be 1D or 2D, but received %d "
                             "dimensions" % targets.ndim)
        if 'int' not in str(targets.dtype):
            raise TypeError("need an integer array for targets")
        if sparse:
            if not scipy_available:
                raise RuntimeError("The converting of indices to a sparse "
                                   "one-hot vector requires scipy to be "
                                   "installed")
            if mode == 'concatenate':
                one_hot = scipy.sparse.csr_matrix(
                    (np.ones(targets.size, dtype=self._dtype),
                     (targets.flatten() + np.arange(targets.size)
                      * self._max_labels)
                     % (self._max_labels * targets.shape[1]),
                     np.arange(targets.shape[0] + 1) * targets.shape[1]),
                    (targets.shape[0], self._max_labels * targets.shape[1])
                )
            elif mode == 'merge':
                one_hot = scipy.sparse.csr_matrix(
                    (np.ones(targets.size), targets.flatten(),
                     np.arange(targets.shape[0] + 1) * targets.shape[1]),
                    (targets.shape[0], self._max_labels)
                )
        else:
            one_hot = np.zeros(targets.shape + (self._max_labels,),
                               dtype=self._dtype)
            shape = (np.prod(one_hot.shape[:-1]), one_hot.shape[-1])
            one_hot.reshape(shape)[np.arange(shape[0]), targets.flatten()] = 1
            if mode == 'concatenate':
                shape = one_hot.shape[-3:-2] + (reduce(mul,
                                                       one_hot.shape[-2:], 1),)
                one_hot = one_hot.reshape(shape)
            elif mode == 'merge':
                one_hot = np.minimum(one_hot.sum(axis=one_hot.ndim - 2), 1)
        return one_hot

    def theano_expr(self, targets, mode='stack', sparse=False):
        """
        Return the one-hot transformation as a symbolic expression.
        If labels appear multiple times, their value in the one-hot
        vector is incremented.

        Parameters
        ----------
        targets : tensor_like, 1- or 2-dimensional, integer dtype
            A symbolic tensor representing labels as integers
            between 0 and `max_labels` - 1, `max_labels` supplied
            at formatter construction.
        mode : string
            The way in which to convert the labels to arrays. Takes
            three different options:

              - "concatenate" : concatenates the one-hot vectors from
                multiple labels
              - "stack" : returns a matrix where each row is the
                one-hot vector of a label
              - "merge" : merges the one-hot vectors together to
                form a vector where the elements are
                the result of an indicator function
                NB: As the result of an indicator function
                the result is the same in case a label
                is duplicated in the input.
        sparse : bool
            If true then the return value is sparse matrix. Note that
            if sparse is True, then mode cannot be 'stack' because
            sparse matrices need to be 2D

        Returns
        -------
        one_hot : TensorVariable, 1, 2 or 3-dimensional, sparse or dense
            A symbolic tensor representing a one-hot encoding of the
            supplied labels.
        """
        if mode not in ('concatenate', 'stack', 'merge'):
            raise ValueError("%s got bad mode argument '%s'" %
                             (self.__class__.__name__, str(self._max_labels)))
        elif mode == 'stack' and sparse:
            raise ValueError("Sparse matrices need to be 2D, hence they"
                             "cannot be stacked")
        squeeze_required = False
        if targets.ndim != 2:
            if targets.ndim == 1:
                squeeze_required = True
                targets = targets.dimshuffle('x', 0)
            else:
                raise ValueError("targets tensor must be 1 or 2-dimensional")
        if 'int' not in str(targets.dtype):
            raise TypeError("need an integer tensor for targets")
        if sparse:
            if mode == 'concatenate':
                one_hot = theano.sparse.CSR(
                    tensor.ones_like(targets, dtype=self._dtype).flatten(),
                    (targets.flatten() + tensor.arange(targets.size) *
                     self._max_labels) % (self._max_labels * targets.shape[1]),
                    tensor.arange(targets.shape[0] + 1) * targets.shape[1],
                    tensor.stack(targets.shape[0],
                                 self._max_labels * targets.shape[1])
                )
            else:
                one_hot = theano.sparse.CSR(
                    tensor.ones_like(targets, dtype=self._dtype).flatten(),
                    targets.flatten(),
                    tensor.arange(targets.shape[0] + 1) * targets.shape[1],
                    tensor.stack(targets.shape[0], self._max_labels)
                )
        else:
            if mode == 'concatenate':
                one_hot = tensor.zeros((targets.shape[0] * targets.shape[1],
                                        self._max_labels), dtype=self._dtype)
                one_hot = tensor.set_subtensor(
                    one_hot[tensor.arange(targets.size),
                            targets.flatten()], 1)
                one_hot = one_hot.reshape(
                    (targets.shape[0], targets.shape[1] * self._max_labels)
                )
            elif mode == 'merge':
                one_hot = tensor.zeros((targets.shape[0], self._max_labels),
                                       dtype=self._dtype)
                one_hot = tensor.set_subtensor(
                    one_hot[tensor.arange(targets.size) % targets.shape[0],
                            targets.T.flatten()], 1)
            else:
                one_hot = tensor.zeros((targets.shape[0], targets.shape[1],
                                        self._max_labels), dtype=self._dtype)
                one_hot = tensor.set_subtensor(one_hot[
                    tensor.arange(targets.shape[0]).reshape((targets.shape[0],
                                                             1)),
                    tensor.arange(targets.shape[1]),
                    targets
                ], 1)
            if squeeze_required:
                if one_hot.ndim == 2:
                    one_hot = one_hot.reshape((one_hot.shape[1],))
                if one_hot.ndim == 3:
                    one_hot = one_hot.reshape((one_hot.shape[1],
                                               one_hot.shape[2]))
        return one_hot

def convert_to_one_hot(integer_vector, dtype=None, max_labels=None,
                       mode='stack', sparse=False):
    """
    Formats a given array of target labels into a one-hot
    vector.

    Parameters
    ----------
    max_labels : int, optional
        The number of possible classes/labels. This means that
        all labels should be < max_labels. Example: For MNIST
        there are 10 numbers and hence max_labels = 10. If not
        given it defaults to max(integer_vector) + 1.
    dtype : dtype, optional
        The desired dtype for the converted one-hot vectors.
        Defaults to config.floatX if not given.
    integer_vector : ndarray
        A 1D array of targets, or a batch (2D array) where
        each row is a list of targets.
    mode : string
        The way in which to convert the labels to arrays. Takes
        three different options:

          - "concatenate" : concatenates the one-hot vectors from
            multiple labels
          - "stack" : returns a matrix where each row is the
            one-hot vector of a label
          - "merge" : merges the one-hot vectors together to
            form a vector where the elements are
            the result of an indicator function
    sparse : bool
        If true then the return value is sparse matrix. Note that
        if sparse is True, then mode cannot be 'stack' because
        sparse matrices need to be 2D

    Returns
    -------
    one_hot : NumPy array
       Can be 1D-3D depending on settings. Normally, the first axis are
       the different batch items, the second axis the labels, the third
       axis the one_hot vectors. Can be dense or sparse.
    """
    if dtype is None:
        dtype = config.floatX
    if isinstance(integer_vector, list):
        integer_vector = np.array(integer_vector)
    assert np.min(integer_vector) >= 0
    assert integer_vector.ndim <= 2
    if max_labels is None:
        max_labels = max(integer_vector) + 1
    return OneHotFormatter(max_labels, dtype=dtype).format(
        integer_vector, mode=mode, sparse=sparse
    )

def _validate_labels(labels, ndim):
    """
    Validate that the passed label is in a right data type, and convert
    it into the desired shape.

    Parameters
    ----------
    labels : array_like, 1-dimensional (or 2-dimensional (nlabels, 1))
        The integer labels to use to construct the one hot matrix.
    ndim : int
        Number of dimensions the label have.

    Returns
    -------
    labels : ndarray, (nlabels, ) or (nlabels, )
        The resulting label vector.
    """
    labels = np.asarray(labels)
    if labels.dtype.kind not in ('u', 'i'):
        raise ValueError("labels must have int or uint dtype")
    if ndim == 1 and labels.ndim != 1:
        if labels.ndim == 2 and labels.shape[1] == 1:
            labels = labels.squeeze()
        else:
            raise ValueError("labels must be 1-dimensional")
    elif ndim == 2 and labels.ndim != 2:
        raise ValueError("labels must be 2-dimensional, no ragged "
                         "lists-of-lists")
    return labels

def compressed_one_hot(labels, dtype=None, out=None, simplify_binary=True,
                       mode='stack', sparse=False):
    """
    Construct a one-hot matrix from a vector of integer labels, but
    only including columns corresponding to integer labels that
    actually appear.

    Parameters
    ----------
    labels : array_like, 1-dimensional (or 2-dimensional (nlabels, 1))
        The integer labels to use to construct the one hot matrix.

    dtype : str or dtype object, optional
        The dtype you wish the returned array to have. Defaults
        to `labels.dtype` if not provided.

    out : ndarray, optional
        An array to use in lieu of allocating one. Must be the
        right shape, i.e. same first dimension as `labels` and
        second dimension greater than or equal to the number of
        unique values in `labels`.

    simplify_binary : bool, optional
        If `True`, if there are only two distinct labels, return
        an `(nlabels, 1)` matrix with 0 denoting the lesser integer
        label and 1 denoting the greater, instead of a redundant
        `(nlabels, 2)` matrix.
    mode : string
        The way in which to convert the labels to arrays. Takes
        three different options:

            - "concatenate" : concatenates the one-hot vectors from
              multiple labels
            - "stack" : returns a matrix where each row is the
              one-hot vector of a label
            - "merge" : merges the one-hot vectors together to
              form a vector where the elements are
              the result of an indicator function
              NB: As the result of an indicator function
              the result is the same in case a label
              is duplicated in the input.
    sparse : bool
        If true then the return value is sparse matrix. Note that
        if sparse is True, then mode cannot be 'stack' because
        sparse matrices need to be 2D

    Returns
    -------
    out : ndarray, (nlabels, max_label + 1) or (nlabels, 1)
        The resulting one-hot matrix.

    uniq : ndarray, 1-dimensional
        The array of unique values in `labels` in the order
        in which the corresponding columns appear in `out`.
    """
    labels = _validate_labels(labels, ndim=1)
    labels_ = labels.copy()
    uniq = np.unique(labels_)
    for i, e in enumerate(uniq):
        labels_[labels_ == e] = i
    if simplify_binary and len(uniq) == 2:
        return labels_.reshape((labels_.shape[0], 1)), uniq
    else:
        return OneHotFormatter(len(uniq), dtype=dtype).format(
            labels_, mode=mode, sparse=sparse), uniq

reason: import numpy,import scipy,import theano,from theano


text:
from __future__ import absolute_import
# Positive Control Enricher class

from math import sqrt
import math
import scipy.stats as s
import array as a

from scipy.optimize import fminbound
from scipy.special import gammaln as gamln
from numpy import *
from numpy import log, pi, log10, e, log1p, exp
import numpy as np

from .MultiSiteCommon import MultiSiteCommon
from .MixtureEstimationMethods import MixtureEstimationMethods

class PositiveControlEnricher(MultiSiteCommon):

    def __init__(self, gbmModel, sequence, rawKinetics):

        MultiSiteCommon.__init__(self, gbmModel, sequence, rawKinetics)
        self.fwd_model = np.genfromtxt("/home/UNIXHOME/obanerjee/initial_lr_model_weights_fwd.csv", delimiter=',')
        self.rev_model = np.genfromtxt("/home/UNIXHOME/obanerjee/initial_lr_model_weights_rev.csv", delimiter=',')
        self.fwd_model = np.squeeze(np.asarray(self.fwd_model))
        self.rev_model = np.squeeze(np.asarray(self.rev_model))

    def fn(self, l):
        if l == "A":
            return 1
        if l == "C":
            return 2
        if l == "G":
            return 3
        return 4

    def tStatisticDenominator(self, mu0, tErr):

        em = 0.01 + 0.03 * mu0 + 0.06 * mu0 ** 1.7
        den = sqrt(em ** 2 + tErr ** 2)
        return den

    def applyLRmodel(self, kinetics, pos, unmodIPDs, modifIPDs, model, up, down, context):
        """ Test out LDA model """

        res = np.zeros((up + down + 1, 7))
        ind = 0

        # range from -down to +up
        for offset in range(-down, (up + 1)):
            a = pos + offset
            tmp = np.squeeze(np.asarray([kinetics[a]["tMean"], kinetics[a]["tErr"], kinetics[a]["coverage"], unmodIPDs[offset + down], modifIPDs[offset + down]]))

            # get t-statistics corresponding to mu0 and mu1:
            den = self.tStatisticDenominator(tmp[3], tmp[1])
            k1 = -(tmp[0] - tmp[3]) / den
            k2 = -(tmp[0] - tmp[4]) / den
            tmp = np.append(np.log(tmp + 0.01), k1)
            tmp = np.append(tmp, k2)

            res[ind, ] = tmp
            ind += 1

        # collected features for prediction:
        apply = np.hstack(res.transpose())
        apply = np.squeeze(np.asarray(apply))
        apply[isnan(apply)] = 0

        # include context:
        context = context[(pos - 10):(pos + 11)]
        del context[11]
        context = np.array(map(self.fn, context))
        apply = np.concatenate([apply, context])

        # calculate logistic regression score:
        z = sum(np.multiply(apply, model[1:])) + model[0]
        score = -z - np.log(1 + np.exp(-z))

        return score

    def callLRstrand(self, kinetics, strand, model, up, down):

        tmp = [d for d in kinetics if d["strand"] == strand]
        tmp.sort(key=lambda x: x["tpl"])

        modSeq = a.array('c')
        modSeq.fromstring(self.sequence)

        L = len(tmp)

        for pos in range(down, (L - up)):

            if tmp[pos]["base"] == 'C':

                # Get H1 means:
                modSeq[pos] = 'I'
                modifIPDs = self.getContextMeans(pos - 10, pos + 10, modSeq)

                # Get H0 means:
                modSeq[pos] = 'C'
                unmodIPDs = self.getContextMeans(pos - 10, pos + 10, modSeq)

                tmp[pos]["Ca5C"] = self.applyLRmodel(tmp, pos, unmodIPDs, modifIPDs, model, up, down, modSeq)

        return tmp

    def callEnricherFunction(self, kinetics, up=10, down=10):

        # Compute all the required mean ipds under all possible composite hypotheses
        self.computeContextMeans()

        fwd = self.callLRstrand(kinetics, 0, self.fwd_model, up, down)
        rev = self.callLRstrand(kinetics, 1, self.rev_model, up, down)
        res = fwd + rev
        res.sort(key=lambda x: x["tpl"])
        return res

    def scoreMods(self, modCalls):
        """
        For each modification in the best scoring configuration, score a config excluding the current mod against the winning config
        use this value as the Qmod for the deleted modification
        """

        qvModCalls = dict()

        modSeq = a.array('c')
        modSeq.fromstring(self.sequence)

        modCalls = [i for i in range(len(modSeq)) if modSeq[i] == 'C']

        for pos in modCalls:

            # now apply the modification at this position:
            modSeq[pos] = 'K'
            modifIPDs = self.getContextMeans(pos - self.post, pos + self.pre, modSeq)

            # using canonical base map, try to get H0 means:
            modSeq[pos] = canonicalBaseMap[call]
            unmodIPDs = self.getContextMeans(pos - self.post, pos + self.pre, modSeq)

            # try to collect related statistics:  tMean, tErr, tStatistic
            tmp = self.applyLRmodel(kinetics, pos, unmodIPDs, modifIPDs, modSeq, up, down)

            # now try to use these vectors to make a basic decision:
            basicDecision[pos] = {'score': tmp}

        return basicDecision

reason: import numpy,from numpy,import scipy,from scipy


text:
import argparse
from pathlib import Path
from typing import Dict, List, Optional, Any

import albumentations as albu
import cv2
import numpy as np
import torch
import torch.nn.parallel
import torch.utils.data
import torch.utils.data.distributed
import yaml
from albumentations.core.serialization import from_dict
from iglovikov_helper_functions.config_parsing.utils import object_from_dict
from iglovikov_helper_functions.dl.pytorch.utils import state_dict_from_disk, tensor_from_rgb_image
from iglovikov_helper_functions.utils.image_utils import load_rgb, pad_to_size, unpad_from_size
from torch.utils.data import Dataset
from torch.utils.data.distributed import DistributedSampler
from tqdm import tqdm

def get_args():
    parser = argparse.ArgumentParser()
    arg = parser.add_argument
    arg("-i", "--input_path", type=Path, help="Path with images.", required=True)
    arg("-c", "--config_path", type=Path, help="Path to config.", required=True)
    arg("-o", "--output_path", type=Path, help="Path to save masks.", required=True)
    arg("-b", "--batch_size", type=int, help="batch_size", default=1)
    arg("-j", "--num_workers", type=int, help="num_workers", default=12)
    arg("-w", "--weight_path", type=str, help="Path to weights.", required=True)
    arg("--world_size", default=-1, type=int, help="number of nodes for distributed training")
    arg("--local_rank", default=-1, type=int, help="node rank for distributed training")
    arg("--fp16", action="store_true", help="Use fp6")
    return parser.parse_args()

class InferenceDataset(Dataset):
    def __init__(self, file_paths: List[Path], transform: albu.Compose) -> None:
        self.file_paths = file_paths
        self.transform = transform

    def __len__(self) -> int:
        return len(self.file_paths)

    def __getitem__(self, idx: int) -> Optional[Dict[str, Any]]:
        image_path = self.file_paths[idx]

        image = load_rgb(image_path)
        height, width = image.shape[:2]

        image = self.transform(image=image)["image"]
        pad_dict = pad_to_size((max(image.shape[:2]), max(image.shape[:2])), image)

        return {
            "torched_image": tensor_from_rgb_image(pad_dict["image"]),
            "image_path": str(image_path),
            "pads": pad_dict["pads"],
            "original_width": width,
            "original_height": height,
        }

def main():
    args = get_args()
    torch.distributed.init_process_group(backend="nccl")

    with open(args.config_path) as f:
        hparams = yaml.load(f, Loader=yaml.SafeLoader)

    hparams.update(
        {
            "local_rank": args.local_rank,
            "fp16": args.fp16,
        }
    )

    output_mask_path = args.output_path / "masks"
    output_mask_path.mkdir(parents=True, exist_ok=True)
    hparams["output_mask_path"] = output_mask_path

    device = torch.device("cuda", args.local_rank)

    model = object_from_dict(hparams["model"])
    model = model.to(device)

    if args.fp16:
        model = model.half()

    corrections: Dict[str, str] = {"model.": ""}
    state_dict = state_dict_from_disk(file_path=args.weight_path, rename_in_layers=corrections)
    model.load_state_dict(state_dict)

    model = torch.nn.parallel.DistributedDataParallel(
        model, device_ids=[args.local_rank], output_device=args.local_rank
    )

    file_paths = []

    for regexp in ["*.jpg", "*.png", "*.jpeg", "*.JPG"]:
        file_paths += sorted([x for x in tqdm(args.input_path.rglob(regexp))])

    dataset = InferenceDataset(file_paths, transform=from_dict(hparams["test_aug"]))

    sampler = DistributedSampler(dataset, shuffle=False)

    dataloader = torch.utils.data.DataLoader(
        dataset,
        batch_size=args.batch_size,
        num_workers=args.num_workers,
        pin_memory=True,
        shuffle=False,
        drop_last=False,
        sampler=sampler,
    )

    predict(dataloader, model, hparams, device)

def predict(dataloader, model, hparams, device):
    model.eval()

    if hparams["local_rank"] == 0:
        loader = tqdm(dataloader)
    else:
        loader = dataloader

    with torch.no_grad():
        for batch in loader:
            torched_images = batch["torched_image"]  # images that are rescaled and padded

            if hparams["fp16"]:
                torched_images = torched_images.half()

            image_paths = batch["image_path"]
            pads = batch["pads"]
            heights = batch["original_height"]
            widths = batch["original_width"]

            batch_size = torched_images.shape[0]

            predictions = model(torched_images.to(device))

            for batch_id in range(batch_size):
                file_id = Path(image_paths[batch_id]).stem

                mask = (predictions[batch_id][0].cpu().numpy() > 0).astype(np.uint8) * 255
                mask = unpad_from_size(pads, image=mask)["image"]
                mask = cv2.resize(
                    mask, (widths[batch_id].item(), heights[batch_id].item()), interpolation=cv2.INTER_NEAREST
                )

                cv2.imwrite(str(hparams["output_mask_path"] / f"{file_id}.png"), mask)

if __name__ == "__main__":
    main()


text:
# License: MIT
# Author: Karl Stelzner

import numpy as np

def z_where_to_bboxes(z_where, conf):
    bbox = np.zeros(list(z_where.shape[:2]) + [4])
    width, height = conf.patch_width, conf.patch_height
    bbox[..., 0] = z_where[..., 5] - 5
    bbox[..., 1] = z_where[..., 2] - 5

    bbox[..., 2] = bbox[..., 0] + height * z_where[..., 4] + 10
    bbox[..., 3] = bbox[..., 1] + width * z_where[..., 0] + 10

    return bbox

def scene_intersection_over_union(n_pred, n_real, bboxes_pred, bboxes_real):
    num_images = bboxes_pred.shape[0]

    score = 0.0
    num_boxes = np.sum(n_real)
    for i in range(num_images):
        matched = [False] * n_real[i]
        for j in range(min(n_real[i], n_pred[i])):
            best_iou = -1
            best_match = -1
            for k in range(n_real[i]):
                if matched[k]:
                    continue
                cur_iou = intersection_over_union(bboxes_pred[i, j],
                                                  bboxes_real[i, k])
                if cur_iou > best_iou:
                    best_iou = cur_iou
                    best_match = k

            if best_match >= 0:
                score += best_iou
                matched[best_match] = True

    return score / num_boxes

def intersection_over_union(boxA, boxB):
    # determine the (x, y)-coordinates of the intersection rectangle
    xA = max(boxA[0], boxB[0])
    yA = max(boxA[1], boxB[1])
    xB = min(boxA[2], boxB[2])
    yB = min(boxA[3], boxB[3])

    # compute the area of intersection rectangle
    interArea = max(0, xB - xA + 1) * max(0, yB - yA + 1)

    # compute the area of both the prediction and ground-truth
    # rectangles
    boxAArea = (boxA[2] - boxA[0] + 1) * (boxA[3] - boxA[1] + 1)
    boxBArea = (boxB[2] - boxB[0] + 1) * (boxB[3] - boxB[1] + 1)

    # compute the intersection over union by taking the intersection
    # area and dividing it by the sum of prediction + ground-truth
    # areas - the interesection area
    iou = interArea / float(boxAArea + boxBArea - interArea)

    # return the intersection over union value
    return iou

max_stars_repo_name: dalessioluca/reproduce-sp-air

max_issues_repo_name: dalessioluca/reproduce-sp-air

max_forks_repo_name: dalessioluca/reproduce-sp-air