text:
# -*- coding: utf-8 -*-

# Copyright (c) 2016-2020 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.

import numpy as np
import pandapower as pp
import pytest

import pandas as pd
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)

def test_convenience_create_functions():
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 110.)
    b2 = pp.create_bus(net, 110.)
    b3 = pp.create_bus(net, 20)
    pp.create_ext_grid(net, b1)
    pp.create_line_from_parameters(net, b1, b2, length_km=20., r_ohm_per_km=0.0487,
                                   x_ohm_per_km=0.1382301, c_nf_per_km=160., max_i_ka=0.664)

    l0 = pp.create_load_from_cosphi(net, b2, 10, 0.95, "ind", name="load")
    pp.runpp(net, init="flat")
    assert net.load.p_mw.at[l0] == 9.5
    assert net.load.q_mvar.at[l0] > 0
    assert np.sqrt(net.load.p_mw.at[l0] ** 2 + net.load.q_mvar.at[l0] ** 2) == 10
    assert np.isclose(net.res_bus.vm_pu.at[b2], 0.99990833838)
    assert net.load.name.at[l0] == "load"

    sh0 = pp.create_shunt_as_capacitor(net, b2, 10, loss_factor=0.01, name="shunt")
    pp.runpp(net, init="flat")
    assert np.isclose(net.res_shunt.q_mvar.at[sh0], -10.043934174)
    assert np.isclose(net.res_shunt.p_mw.at[sh0], 0.10043933665)
    assert np.isclose(net.res_bus.vm_pu.at[b2], 1.0021942964)
    assert net.shunt.name.at[sh0] == "shunt"

    sg0 = pp.create_sgen_from_cosphi(net, b2, 5, 0.95, "cap", name="sgen")
    pp.runpp(net, init="flat")
    assert np.sqrt(net.sgen.p_mw.at[sg0] ** 2 + net.sgen.q_mvar.at[sg0] ** 2) == 5
    assert net.sgen.p_mw.at[sg0] == 4.75
    assert net.sgen.q_mvar.at[sg0] > 0
    assert np.isclose(net.res_bus.vm_pu.at[b2], 1.0029376578)
    assert net.sgen.name.at[sg0] == "sgen"

    tol = 1e-6
    base_z = 110 ** 2 / 100
    sind = pp.create_series_reactor_as_impedance(net, b1, b2, r_ohm=100, x_ohm=200, sn_mva=100)
    assert net.impedance.at[sind, 'rft_pu'] - 100 / base_z < tol
    assert net.impedance.at[sind, 'xft_pu'] - 200 / base_z < tol

    tid = pp.create_transformer_from_parameters(net, hv_bus=b2, lv_bus=b3, sn_mva=0.1, vn_hv_kv=110,
                                                vn_lv_kv=20, vkr_percent=5, vk_percent=20,
                                                pfe_kw=1, i0_percent=1)
    pp.create_load(net, b3, 0.1)
    assert net.trafo.at[tid, 'df'] == 1
    pp.runpp(net)
    tr_l = net.res_trafo.at[tid, 'loading_percent']
    net.trafo.at[tid, 'df'] = 2
    pp.runpp(net)
    tr_l_2 = net.res_trafo.at[tid, 'loading_percent']
    assert tr_l == tr_l_2 * 2
    net.trafo.at[tid, 'df'] = 0
    with pytest.raises(UserWarning):
        pp.runpp(net)

def test_nonexistent_bus():
    from functools import partial
    net = pp.create_empty_network()
    create_functions = [partial(pp.create_load, net=net, p_mw=0, q_mvar=0, bus=0, index=0),
                        partial(pp.create_sgen, net=net, p_mw=0, q_mvar=0, bus=0, index=0),
                        partial(pp.create_dcline, net, from_bus=0, to_bus=1, p_mw=0.1,
                                loss_percent=0, loss_mw=0.01, vm_from_pu=1., vm_to_pu=1., index=0),
                        partial(pp.create_gen, net=net, p_mw=0, bus=0, index=0),
                        partial(pp.create_ward, net, 0, 0, 0, 0, 0, index=0),
                        partial(pp.create_xward, net, 0, 0, 0, 0, 0, 1, 1, 1, index=0),
                        partial(pp.create_shunt, net=net, q_mvar=0, bus=0, index=0),
                        partial(pp.create_ext_grid, net=net, bus=1, index=0),
                        partial(pp.create_line, net=net, from_bus=0, to_bus=1, length_km=1.,
                                std_type="NAYY 4x50 SE", index=0),
                        partial(pp.create_line_from_parameters, net=net, from_bus=0, to_bus=1,
                                length_km=1., r_ohm_per_km=0.1, x_ohm_per_km=0.1, max_i_ka=0.4,
                                c_nf_per_km=10, index=1),
                        partial(pp.create_transformer, net=net, hv_bus=0, lv_bus=1,
                                std_type="63 MVA 110/20 kV", index=0),
                        partial(pp.create_transformer3w, net=net, hv_bus=0, lv_bus=1, mv_bus=2,
                                std_type="63/25/38 MVA 110/20/10 kV", index=0),
                        partial(pp.create_transformer3w_from_parameters, net=net, hv_bus=0,
                                lv_bus=1, mv_bus=2, i0_percent=0.89, pfe_kw=3.5,
                                vn_hv_kv=110, vn_lv_kv=10, vn_mv_kv=20, sn_hv_mva=63,
                                sn_lv_mva=38, sn_mv_mva=25, vk_hv_percent=10.4,
                                vk_lv_percent=10.4, vk_mv_percent=10.4, vkr_hv_percent=0.28,
                                vkr_lv_percent=0.35, vkr_mv_percent=0.32, index=1),
                        partial(pp.create_transformer_from_parameters, net=net, hv_bus=0, lv_bus=1,
                                sn_mva=60, vn_hv_kv=20., vn_lv_kv=0.4, vk_percent=10,
                                vkr_percent=0.1, pfe_kw=0, i0_percent=0, index=1),
                        partial(pp.create_impedance, net=net, from_bus=0, to_bus=1,
                                rft_pu=0.1, xft_pu=0.1, sn_mva=0.6, index=0),
                        partial(pp.create_switch, net, bus=0, element=1, et="b", index=0)]
    for func in create_functions:
        with pytest.raises(Exception):  # exception has to be raised since bus doesn't exist
            func()
    pp.create_bus(net, 0.4)
    pp.create_bus(net, 0.4)
    pp.create_bus(net, 0.4)
    for func in create_functions:
        func()  # buses exist, element can be created
        with pytest.raises(Exception):  # exception is raised because index already exists
            func()

def test_tap_phase_shifter_default():
    expected_default = False
    net = pp.create_empty_network()
    pp.create_bus(net, 110)
    pp.create_bus(net, 20)
    data = pp.load_std_type(net, "25 MVA 110/20 kV", "trafo")
    if "tap_phase_shifter" in data:
        del data["tap_phase_shifter"]
    pp.create_std_type(net, data, "without_tap_shifter_info", "trafo")
    pp.create_transformer_from_parameters(net, 0, 1, 25e3, 110, 20, 0.4, 12, 20, 0.07)
    pp.create_transformer(net, 0, 1, "without_tap_shifter_info")
    assert (net.trafo.tap_phase_shifter == expected_default).all()

def test_create_line_conductance():
    net = pp.create_empty_network()
    pp.create_bus(net, 20)
    pp.create_bus(net, 20)
    pp.create_std_type(net, {'c_nf_per_km': 210, 'max_i_ka': 0.142, 'q_mm2': 50,
                             'r_ohm_per_km': 0.642, 'type': 'cs', 'x_ohm_per_km': 0.083,
                             "g_us_per_km": 1}, "test_conductance")

    l = pp.create_line(net, 0, 1, 1., "test_conductance")
    assert net.line.g_us_per_km.at[l] == 1

def test_create_buses():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_buses(net, 3, 110)
    # with geodata
    b2 = pp.create_buses(net, 3, 110, geodata=(10, 20))
    # with geodata as array
    geodata = np.array([[10, 20], [20, 30], [30, 40]])
    b3 = pp.create_buses(net, 3, 110, geodata=geodata)

    assert len(net.bus) == 9
    assert len(net.bus_geodata) == 6

    for i in b2:
        assert net.bus_geodata.at[i, 'x'] == 10
        assert net.bus_geodata.at[i, 'y'] == 20

    assert (net.bus_geodata.loc[b3, ['x', 'y']].values == geodata).all()

    # no way of creating buses with not matching shape
    with pytest.raises(ValueError):
        pp.create_buses(net, 2, 110, geodata=geodata)

def test_create_lines():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines(net, [b1, b1], [b2, b2], 4, std_type="48-AL1/8-ST1A 10.0")
    assert len(net.line) == 2
    assert len(net.line_geodata) == 0
    assert sum(net.line.std_type == "48-AL1/8-ST1A 10.0") == 2
    assert len(set(net.line.r_ohm_per_km)) == 1

    # with geodata
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines(net, [b1, b1], [b2, b2], [1.5, 3], std_type="48-AL1/8-ST1A 10.0",
                         geodata=[[(1,1),(2,2),(3,3)], [(1,1),(1,2)]])

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert net.line_geodata.at[l[0], "coords"] == [(1,1),(2,2),(3,3)]
    assert net.line_geodata.at[l[1], "coords"] == [(1,1),(1,2)]

    # setting params as single value
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines(net, [b1, b1], [b2, b2], length_km=5, df=0.8, in_service=False,
                        geodata=[(10, 10), (20, 20)], parallel=1, max_loading_percent=90,
                        name="test", std_type="48-AL1/8-ST1A 10.0")

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert net.line.length_km.at[l[0]] == 5
    assert net.line.length_km.at[l[1]] == 5
    assert net.line.at[l[0], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line.at[l[1], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line_geodata.at[l[0], "coords"] == [(10,10), (20,20)]
    assert net.line_geodata.at[l[1], "coords"] == [(10,10), (20,20)]
    assert net.line.at[l[0], "name"] == "test"
    assert net.line.at[l[1], "name"] == "test"
    assert net.line.at[l[0], "max_loading_percent"] == 90
    assert net.line.at[l[1], "max_loading_percent"] == 90
    assert net.line.at[l[0], "parallel"] == 1
    assert net.line.at[l[1], "parallel"] == 1

    # setting params as array
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines(net, [b1, b1], [b2, b2], length_km=[1, 5], df=[0.8, 0.7],
                        in_service=[True, False],
                        geodata=[[(10, 10), (20, 20)], [(100, 10), (200, 20)]], parallel=[2, 1],
                        max_loading_percent=[80, 90], name=["test1", "test2"],
                        std_type="48-AL1/8-ST1A 10.0")

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert net.line.at[l[0], "length_km"] == 1
    assert net.line.at[l[1], "length_km"] == 5
    assert net.line.at[l[0], "in_service"] == True  # is actually <class 'numpy.bool_'>
    assert net.line.at[l[1], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line_geodata.at[l[0], "coords"] == [(10,10), (20,20)]
    assert net.line_geodata.at[l[1], "coords"] == [(100,10), (200,20)]
    assert net.line.at[l[0], "name"] == "test1"
    assert net.line.at[l[1], "name"] == "test2"
    assert net.line.at[l[0], "max_loading_percent"] == 80
    assert net.line.at[l[1], "max_loading_percent"] == 90
    assert net.line.at[l[0], "parallel"] == 2
    assert net.line.at[l[1], "parallel"] == 1

def test_create_lines_from_parameters():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], length_km=[10., 5.], x_ohm_per_km=[1., 1.],
                                        r_ohm_per_km=[0.2, 0.2], c_nf_per_km=[0, 0], max_i_ka=[100, 100])
    assert len(net.line) == 2
    assert len(net.line_geodata) == 0
    assert len(net.line.x_ohm_per_km) == 2
    assert len(net.line.r_ohm_per_km) == 2
    assert len(net.line.c_nf_per_km) == 2
    assert len(net.line.max_i_ka) == 2
    assert len(net.line.df) == 2

    # with geodata
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], length_km=[10., 5.], x_ohm_per_km=[1., 1.],
                                        r_ohm_per_km=[0.2, 0.2], c_nf_per_km=[0, 0], max_i_ka=[100, 100],
                                        geodata=[[(1,1),(2,2),(3,3)], [(1,1),(1,2)]])

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert net.line_geodata.at[l[0], "coords"] == [(1,1),(2,2),(3,3)]
    assert net.line_geodata.at[l[1], "coords"] == [(1,1),(1,2)]

    # setting params as single value
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], length_km=5, x_ohm_per_km=1,
                                        r_ohm_per_km=0.2, c_nf_per_km=0, max_i_ka=100, df=0.8, in_service=False,
                                        geodata=[(10, 10), (20, 20)], parallel=1, max_loading_percent=90, name='test',
                                        r0_ohm_per_km=0.1, g0_us_per_km=0., c0_nf_per_km=0.,
                                        temperature_degree_celsius=20, alpha=0.04)

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert all(net.line['length_km'].values == 5)
    assert all(net.line['x_ohm_per_km'].values == 1)
    assert all(net.line['r_ohm_per_km'].values == 0.2)
    assert all(net.line["r0_ohm_per_km"].values == 0.1)
    assert all(net.line["g0_us_per_km"].values == 0)
    assert all(net.line["c0_nf_per_km"].values == 0)
    assert net.line.at[l[0], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line.at[l[1], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line_geodata.at[l[0], "coords"] == [(10,10), (20,20)]
    assert net.line_geodata.at[l[1], "coords"] == [(10,10), (20,20)]
    assert all(net.line["name"].values == "test")
    assert all(net.line["max_loading_percent"].values == 90)
    assert all(net.line["parallel"].values == 1)
    assert all(net.line["temperature_degree_celsius"].values == 20.)
    assert all(net.line["alpha"].values == 0.04)

    # setting params as array
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    l = pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], length_km=[1, 5], r_ohm_per_km=[1, 2],
                                        x_ohm_per_km=[0.3, 0.5], c_nf_per_km=[0., 0.1], r0_ohm_per_km=[0.1, 0.15],
                                        x0_ohm_per_km=[0.2, 0.25], g0_us_per_km=[0., 0.], c0_nf_per_km=[0., 0.],
                                        df=[0.8, 0.7], in_service=[True, False],
                                        geodata=[[(10, 10), (20, 20)], [(100, 10), (200, 20)]], parallel=[2, 1],
                                        max_loading_percent=[80, 90], name=["test1", "test2"], max_i_ka=[100, 200])

    assert len(net.line) == 2
    assert len(net.line_geodata) == 2
    assert net.line.at[l[0], "length_km"] == 1
    assert net.line.at[l[1], "length_km"] == 5
    assert net.line.at[l[0], "r_ohm_per_km"] == 1
    assert net.line.at[l[1], "r_ohm_per_km"] == 2
    assert net.line.at[l[0], "x_ohm_per_km"] == 0.3
    assert net.line.at[l[1], "x_ohm_per_km"] == 0.5
    assert net.line.at[l[0], "c_nf_per_km"] == 0.
    assert net.line.at[l[1], "c_nf_per_km"] == 0.1
    assert net.line.at[l[0], "r0_ohm_per_km"] == 0.1
    assert net.line.at[l[1], "r0_ohm_per_km"] == 0.15
    assert net.line.at[l[0], "x0_ohm_per_km"] == 0.2
    assert net.line.at[l[1], "x0_ohm_per_km"] == 0.25
    assert all(net.line["g0_us_per_km"].values == 0)
    assert all(net.line["c0_nf_per_km"].values == 0)
    assert net.line.at[l[0], "in_service"] == True  # is actually <class 'numpy.bool_'>
    assert net.line.at[l[1], "in_service"] == False  # is actually <class 'numpy.bool_'>
    assert net.line_geodata.at[l[0], "coords"] == [(10,10), (20,20)]
    assert net.line_geodata.at[l[1], "coords"] == [(100,10), (200,20)]
    assert net.line.at[l[0], "name"] == "test1"
    assert net.line.at[l[1], "name"] == "test2"
    assert net.line.at[l[0], "max_loading_percent"] == 80
    assert net.line.at[l[1], "max_loading_percent"] == 90
    assert net.line.at[l[0], "parallel"] == 2
    assert net.line.at[l[1], "parallel"] == 1
    assert net.line.at[l[0], "max_i_ka"] == 100
    assert net.line.at[l[1], "max_i_ka"] == 200

def test_create_lines_raise_except():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], length_km=[10., 5.], x_ohm_per_km=[1., 1.],
                                    r_ohm_per_km=[0.2, 0.2], c_nf_per_km=[0, 0], max_i_ka=[100, 100])

    with pytest.raises(UserWarning, match="Lines trying to attach .*"):
        pp.create_lines_from_parameters(net, [b1, 2], [2, b2], length_km=[10., 5.], x_ohm_per_km=[1., 1.],
                                        r_ohm_per_km=[0.2, 0.2], c_nf_per_km=[0, 0], max_i_ka=[100, 100])

    with pytest.raises(UserWarning, match="Lines with indexes .*"):
        pp.create_lines_from_parameters(net, [b1, b1], [b2, b2], index=[0, 0], length_km=[10., 5.], x_ohm_per_km=[1., 1.],
                                        r_ohm_per_km=[0.2, 0.2], c_nf_per_km=[0, 0], max_i_ka=[100, 100])

def test_create_line_alpha_temperature():
    net=pp.create_empty_network()
    b = pp.create_buses(net, 5, 110)

    l1=pp.create_line(net,0,1, 10, "48-AL1/8-ST1A 10.0")
    l2=pp.create_line(net,1,2, 10, "48-AL1/8-ST1A 10.0", alpha=4.03e-3, temperature_degree_celsius=80)
    l3=pp.create_line(net,2,3, 10, "48-AL1/8-ST1A 10.0")
    l4=pp.create_line_from_parameters(net, 3,4,10, 1,1,1,100)
    l5=pp.create_line_from_parameters(net, 3,4,10, 1,1,1,100, alpha=4.03e-3)

    assert 'alpha' in net.line.columns
    assert all(net.line.loc[[l2,l3,l5], 'alpha'] == 4.03e-3)
    assert all(net.line.loc[[l1,l4], 'alpha'].isnull())
    assert net.line.loc[l2, 'temperature_degree_celsius'] == 80
    assert all(net.line.loc[[l1,l3,l4,l5], 'temperature_degree_celsius'].isnull())

def test_create_transformers_from_parameters():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 15)
    b2 = pp.create_bus(net, 0.4)
    t = pp.create_transformers_from_parameters(net, [b1, b1], [b2, b2], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                               sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                               i0_percent=0.3, foo=2)
    assert len(net.trafo) == 2
    assert len(net.trafo.vk_percent) == 2
    assert len(net.trafo.vkr_percent) == 2
    assert len(net.trafo.pfe_kw) == 2
    assert len(net.trafo.i0_percent) == 2
    assert len(net.trafo.df) == 2
    assert len(net.trafo.foo) == 2

    # setting params as single value
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 15)
    b2 = pp.create_bus(net, 0.4)
    pp.create_transformers_from_parameters(net, hv_buses=[b1, b1], lv_buses=[b2, b2], vn_hv_kv=15., vn_lv_kv=0.45,
                                               sn_mva= 0.5, vk_percent=1., vkr_percent=0.3, pfe_kw=0.2, i0_percent=0.3,
                                               vk0_percent=0.4, vkr0_percent=1.7,  mag0_rx=0.4, mag0_percent=0.3,
                                               tap_neutral=0., vector_group='Dyn', si0_hv_partial=0.1,
                                               max_loading_percent=80)
    assert len(net.trafo) == 2
    assert all(net.trafo.hv_bus == 0)
    assert all(net.trafo.lv_bus == 1)
    assert all(net.trafo.sn_mva == 0.5)
    assert all(net.trafo.vn_hv_kv == 15.)
    assert all(net.trafo.vn_lv_kv == 0.45)
    assert all(net.trafo.vk_percent == 1.)
    assert all(net.trafo.vkr_percent == 0.3)
    assert all(net.trafo.pfe_kw == 0.2)
    assert all(net.trafo.i0_percent == 0.3)
    assert all(net.trafo.vk0_percent == 0.4)
    assert all(net.trafo.mag0_rx == 0.4)
    assert all(net.trafo.mag0_percent == 0.3)
    assert all(net.trafo.tap_neutral == 0.)
    assert all(net.trafo.tap_pos == 0.)
    assert all(net.trafo.vector_group == 'Dyn')
    assert all(net.trafo.max_loading_percent == 80.)
    assert all(net.trafo.si0_hv_partial == 0.1)

    # setting params as array
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    t = pp.create_transformers_from_parameters(net, hv_buses=[b1, b1], lv_buses=[b2, b2], vn_hv_kv=[15., 15.], sn_mva=[0.6, 0.6],
                                               vn_lv_kv=[0.45, 0.45], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3],
                                               pfe_kw=[0.2, 0.2], i0_percent=[0.3, 0.3],  vk0_percent=[0.4, 0.4],
                                               mag0_rx=[0.4, 0.4], mag0_percent=[0.3, 0.3], tap_neutral=[0., 1.],
                                               tap_pos=[-1, 4])
    print(net.trafo.tap_neutral)
    assert len(net.trafo) == 2
    assert all(net.trafo.hv_bus == 0)
    assert all(net.trafo.lv_bus == 1)
    assert all(net.trafo.vn_hv_kv == 15.)
    assert all(net.trafo.vn_lv_kv == 0.45)
    assert all(net.trafo.sn_mva == 0.6)
    assert all(net.trafo.vk_percent == 1.)
    assert all(net.trafo.vkr_percent == 0.3)
    assert all(net.trafo.pfe_kw == 0.2)
    assert all(net.trafo.i0_percent == 0.3)
    assert all(net.trafo.vk0_percent == 0.4)
    assert all(net.trafo.mag0_rx == 0.4)
    assert all(net.trafo.mag0_percent == 0.3)
    assert net.trafo.tap_neutral.at[t[0]] == 0
    assert net.trafo.tap_neutral.at[t[1]] == 1
    assert net.trafo.tap_pos.at[t[0]] == -1
    assert net.trafo.tap_pos.at[t[1]] == 4

def test_create_transformers_raise_except():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    pp.create_transformers_from_parameters(net, [b1, b1], [b2, b2], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                           sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                           i0_percent=0.3, foo=2)

    with pytest.raises(UserWarning, match="A trafo with index 1 already exists"):
        pp.create_transformers_from_parameters(net, [b1, b1], [b2, b2], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                               sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                               i0_percent=0.3, index=[2, 1])
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 10)
    pp.create_transformers_from_parameters(net, [b1, b1], [b2, b2], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                           sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                           i0_percent=0.3, foo=2)
    with pytest.raises(UserWarning, match="Transformer trying to attach to non existing buses \[2\]"):
        pp.create_transformers_from_parameters(net, [b1, 2], [b2, b2], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                               sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                               i0_percent=0.3, foo=2)

    with pytest.raises(UserWarning, match="Transformer trying to attach to non existing buses \[3\]"):
        pp.create_transformers_from_parameters(net, [b1, b1], [b2, 3], vn_hv_kv=[15., 15.], vn_lv_kv=[0.45, 0.45],
                                                sn_mva=[0.5, 0.7], vk_percent=[1., 1.], vkr_percent=[0.3, 0.3], pfe_kw=0.2,
                                                i0_percent=0.3, foo=2)

def test_create_transformers3w_from_parameters():

    # setting params as single value
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 15)
    b2 = pp.create_bus(net, 0.4)
    b3 = pp.create_bus(net, 0.9)
    pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, b3], lv_buses=[b2, b2],
                                               vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                               sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                               vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                               pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                               mag0_rx=0.4, mag0_percent=0.3)
    assert len(net.trafo3w) == 2
    assert all(net.trafo3w.hv_bus == 0)
    assert all(net.trafo3w.lv_bus == 1)
    assert all(net.trafo3w.mv_bus == 2)
    assert all(net.trafo3w.sn_hv_mva == 0.6)
    assert all(net.trafo3w.sn_mv_mva == 0.5)
    assert all(net.trafo3w.sn_lv_mva == 0.4)
    assert all(net.trafo3w.vn_hv_kv == 15.)
    assert all(net.trafo3w.vn_mv_kv == 0.9)
    assert all(net.trafo3w.vn_lv_kv == 0.45)
    assert all(net.trafo3w.vk_hv_percent == 1.)
    assert all(net.trafo3w.vk_mv_percent == 1.)
    assert all(net.trafo3w.vk_lv_percent == 1.)
    assert all(net.trafo3w.vkr_hv_percent == 0.3)
    assert all(net.trafo3w.vkr_mv_percent == 0.3)
    assert all(net.trafo3w.vkr_lv_percent == 0.3)
    assert all(net.trafo3w.pfe_kw == 0.2)
    assert all(net.trafo3w.i0_percent == 0.3)
    assert all(net.trafo3w.mag0_rx == 0.4)
    assert all(net.trafo3w.mag0_percent == 0.3)
    assert all(net.trafo3w.tap_neutral == 0.)
    assert all(net.trafo3w.tap_pos == 0.)

    # setting params as array
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 10)
    b2 = pp.create_bus(net, 0.4)
    b3 = pp.create_bus(net, 0.9)
    pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, b3], lv_buses=[b2, b2],
                                                 vn_hv_kv=[15., 14.5], vn_mv_kv=[0.9, 0.7], vn_lv_kv=[0.45, 0.5],
                                                 sn_hv_mva= [0.6, 0.7], sn_mv_mva=[0.5, 0.4], sn_lv_mva=[0.4, 0.3],
                                                 vk_hv_percent=[1., 1.], vk_mv_percent=[1., 1.], vk_lv_percent=[1.,1.],
                                                 vkr_hv_percent=[0.3, 0.3], vkr_mv_percent=[0.3, 0.3],
                                                 vkr_lv_percent=[0.3, 0.3], pfe_kw=[0.2, 0.1], i0_percent=[0.3, 0.2],
                                                 tap_neutral=[0., 5.], tap_pos=[1, 2], in_service=[True, False],
                                                 custom_arg=['foo', 'bar']
                                                 )
    assert len(net.trafo3w) == 2
    assert all(net.trafo3w.hv_bus == 0)
    assert all(net.trafo3w.lv_bus == 1)
    assert all(net.trafo3w.mv_bus == 2)
    assert all(net.trafo3w.sn_hv_mva == [0.6, 0.7])
    assert all(net.trafo3w.sn_mv_mva == [0.5, 0.4])
    assert all(net.trafo3w.sn_lv_mva == [0.4, 0.3])
    assert all(net.trafo3w.vn_hv_kv == [15., 14.5])
    assert all(net.trafo3w.vn_mv_kv == [0.9, 0.7])
    assert all(net.trafo3w.vn_lv_kv == [0.45, 0.5])
    assert all(net.trafo3w.vk_hv_percent == 1.)
    assert all(net.trafo3w.vk_mv_percent == 1.)
    assert all(net.trafo3w.vk_lv_percent == 1.)
    assert all(net.trafo3w.vkr_hv_percent == 0.3)
    assert all(net.trafo3w.vkr_mv_percent == 0.3)
    assert all(net.trafo3w.vkr_lv_percent == 0.3)
    assert all(net.trafo3w.pfe_kw == [0.2, 0.1])
    assert all(net.trafo3w.i0_percent == [0.3, 0.2])
    assert all(net.trafo3w.tap_neutral == [0., 5.])
    assert all(net.trafo3w.tap_pos == [1, 2])
    assert all(net.trafo3w.in_service == [True, False])
    assert all(net.trafo3w.custom_arg == ['foo', 'bar'])

def test_create_transformers3w_raise_except():
    # standard
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 15)
    b2 = pp.create_bus(net, 0.4)
    b3 = pp.create_bus(net, 0.9)
    pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, b3], lv_buses=[b2, b2],
                                             vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                             sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                             vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                             pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                             mag0_rx=0.4, mag0_percent=0.3)

    with pytest.raises(UserWarning, match="A three winding transformer with index 1 already exists"):
        pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, b3], lv_buses=[b2, b2],
                                                 vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                                 sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                                 vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                                 pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                                 mag0_rx=0.4, mag0_percent=0.3, index=[2, 1])
    net = pp.create_empty_network()
    b1 = pp.create_bus(net, 15)
    b2 = pp.create_bus(net, 0.4)
    b3 = pp.create_bus(net, 0.9)
    with pytest.raises(UserWarning, match="Transformer trying to attach to non existing buses \{6\}"):
        pp.create_transformers3w_from_parameters(net, hv_buses=[6, b1], mv_buses=[b3, b3], lv_buses=[b2, b2],
                                                 vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                                 sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                                 vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                                 pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                                 mag0_rx=0.4, mag0_percent=0.3, index=[0, 1])

    with pytest.raises(UserWarning, match="Transformer trying to attach to non existing buses \{3\}"):
        pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, 3], lv_buses=[b2, b2],
                                                 vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                                 sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                                 vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                                 pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                                 mag0_rx=0.4, mag0_percent=0.3)
    with pytest.raises(UserWarning, match="Transformer trying to attach to non existing buses \{3, 4\}"):
        pp.create_transformers3w_from_parameters(net, hv_buses=[b1, b1], mv_buses=[b3, b3], lv_buses=[4, 3],
                                                vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                                sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                                vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                                pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.,
                                                mag0_rx=0.4, mag0_percent=0.3)

def test_create_switches():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 15)
    b4 = pp.create_bus(net, 15)
    l1 = pp.create_line(net, b1, b2,  length_km=1, std_type="48-AL1/8-ST1A 10.0")
    t1 = pp.create_transformer(net, b2, b3, std_type='160 MVA 380/110 kV')

    sw = pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, t1, b4], et=["l", "t", "b"], z_ohm=0., foo='aaa')

    assert(net.switch.bus.at[0] == b1)
    assert(net.switch.bus.at[1] == b2)
    assert(net.switch.bus.at[2] == b3)
    assert(net.switch.element.at[sw[0]] == l1)
    assert(net.switch.element.at[sw[1]] == t1)
    assert(net.switch.element.at[sw[2]] == b4)
    assert(net.switch.et.at[0] == "l")
    assert(net.switch.et.at[1] == "t")
    assert(net.switch.et.at[2] == "b")
    assert(net.switch.z_ohm.at[0] == 0)
    assert(net.switch.z_ohm.at[1] == 0)
    assert(net.switch.z_ohm.at[2] == 0)
    assert(net.switch.foo.at[0] == 'aaa')
    assert(net.switch.foo.at[1] == 'aaa')
    assert(net.switch.foo.at[2] == 'aaa')

def test_create_switches_raise_except():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 15)
    b4 = pp.create_bus(net, 15)
    b5 = pp.create_bus(net, 0.9)
    b6 = pp.create_bus(net, 0.4)
    l1 = pp.create_line(net, b1, b2,  length_km=1, std_type="48-AL1/8-ST1A 10.0")
    t1 = pp.create_transformer(net, b2, b3, std_type='160 MVA 380/110 kV')
    t3w1 = pp.create_transformer3w_from_parameters(
                                                net,
                                                hv_bus=b4, mv_bus=b5, lv_bus=b6,
                                                vn_hv_kv=15., vn_mv_kv=0.9, vn_lv_kv=0.45, sn_hv_mva= 0.6, sn_mv_mva=0.5,
                                                sn_lv_mva=0.4, vk_hv_percent=1., vk_mv_percent=1., vk_lv_percent=1.,
                                                vkr_hv_percent=0.3, vkr_mv_percent=0.3, vkr_lv_percent=0.3,
                                                pfe_kw=0.2, i0_percent=0.3, tap_neutral=0.
                                            )
    sw = pp.create_switch(net, bus=b1, element=l1, et="l", z_ohm=0.)
    with pytest.raises(UserWarning, match="A switch with index 0 already exists"):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, t1, b4], et=["l", "t", "b"], z_ohm=0.,
                           index=[sw, 1, 2])

    with pytest.raises(UserWarning, match="Buses \{6\} do not exist"):
        pp.create_switches(net, buses=[6, b2, b3], elements=[l1, t1, b4], et=["l", "t", "b"], z_ohm=0.)

    with pytest.raises(UserWarning, match="Line 1 does not exist"):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[1, t1, b4], et=["l", "t", "b"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Line %s not connected to bus %s" % (l1, b3)):
        pp.create_switches(net, buses=[b3, b2, b3], elements=[l1, t1, b4], et=["l", "t", "b"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Trafo 1 does not exist"):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, 1, b4], et=["l", "t", "b"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Trafo %s not connected to bus %s" % (t1, b1)):
        pp.create_switches(net, buses=[b1, b1, b3], elements=[l1, t1, b4], et=["l", "t", "b"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Unknown bus index 6"):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, t1, 6], et=["l", "t", "b"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Trafo3w 1 does not exist"):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, t1, 1], et=["l", "t", "t3"], z_ohm=0.)
    with pytest.raises(UserWarning, match="Trafo3w %s not connected to bus %s" % (t3w1, b3)):
        pp.create_switches(net, buses=[b1, b2, b3], elements=[l1, t1, t3w1], et=["l", "t", "t3"], z_ohm=0.)

def test_create_loads():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)
    pp.create_loads(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                    max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                    )

    assert(net.load.bus.at[0] == b1)
    assert(net.load.bus.at[1] == b2)
    assert(net.load.bus.at[2] == b3)
    assert(net.load.p_mw.at[0] == 0)
    assert(net.load.p_mw.at[1] == 0)
    assert(net.load.p_mw.at[2] == 1)
    assert(net.load.q_mvar.at[0] == 0)
    assert(net.load.q_mvar.at[1] == 0)
    assert(net.load.q_mvar.at[2] == 0)
    assert(net.load.controllable.at[0] == True)
    assert(net.load.controllable.at[1] == False)
    assert(net.load.controllable.at[2] == False)
    assert(all(net.load.max_p_mw.values == 0.2))
    assert(all(net.load.min_p_mw.values == [0, 0.1, 0]))
    assert(all(net.load.max_q_mvar.values == 0.2))
    assert(all(net.load.min_q_mvar.values == [0, 0.1, 0]))

def test_create_loads_raise_except():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)

    with pytest.raises(UserWarning, match="Cannot attach to buses \{3, 4, 5\}, they does not exist"):
        pp.create_loads(net, buses=[3, 4, 5], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        )

    l = pp.create_loads(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        )
    with pytest.raises(UserWarning, match="Loads with the ids \[0 1 2\] already exists"):
        pp.create_loads(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        index = l)

def test_create_sgens():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)
    pp.create_sgens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                    max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                    k=1.3, rx=0.4, current_source=True)

    assert(net.sgen.bus.at[0] == b1)
    assert(net.sgen.bus.at[1] == b2)
    assert(net.sgen.bus.at[2] == b3)
    assert(net.sgen.p_mw.at[0] == 0)
    assert(net.sgen.p_mw.at[1] == 0)
    assert(net.sgen.p_mw.at[2] == 1)
    assert(net.sgen.q_mvar.at[0] == 0)
    assert(net.sgen.q_mvar.at[1] == 0)
    assert(net.sgen.q_mvar.at[2] == 0)
    assert(net.sgen.controllable.at[0] == True)
    assert(net.sgen.controllable.at[1] == False)
    assert(net.sgen.controllable.at[2] == False)
    assert(all(net.sgen.max_p_mw.values == 0.2))
    assert(all(net.sgen.min_p_mw.values == [0, 0.1, 0]))
    assert(all(net.sgen.max_q_mvar.values == 0.2))
    assert(all(net.sgen.min_q_mvar.values == [0, 0.1, 0]))
    assert(all(net.sgen.k.values == 1.3))
    assert(all(net.sgen.rx.values == 0.4))
    assert(all(net.sgen.current_source))

def test_create_sgens_raise_except():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)

    with pytest.raises(UserWarning, match="Cannot attach to buses \{3, 4, 5\}, they does not exist"):
        pp.create_sgens(net, buses=[3, 4, 5], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        k=1.3, rx=0.4, current_source=True)

    sg = pp.create_sgens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        k=1.3, rx=0.4, current_source=True)
    with pytest.raises(UserWarning, match="Sgens with the ids \[0 1 2\] already exists"):
        pp.create_sgens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], q_mwar=0., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        k=1.3, rx=0.4, current_source=True, index=sg)

def test_create_gens():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)
    pp.create_gens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], vm_pu=1., controllable=[True, False, False],
                   max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                   min_vm_pu=0.85, max_vm_pu=1.15, vn_kv=0.4, xdss_pu=0.1, rdss_pu=0.1, cos_phi=1.
                )
    assert(net.gen.bus.at[0] == b1)
    assert(net.gen.bus.at[1] == b2)
    assert(net.gen.bus.at[2] == b3)
    assert(net.gen.p_mw.at[0] == 0)
    assert(net.gen.p_mw.at[1] == 0)
    assert(net.gen.p_mw.at[2] == 1)
    assert(net.gen.controllable.at[0] == True)
    assert(net.gen.controllable.at[1] == False)
    assert(net.gen.controllable.at[2] == False)
    assert(all(net.gen.max_p_mw.values == 0.2))
    assert(all(net.gen.min_p_mw.values == [0, 0.1, 0]))
    assert(all(net.gen.max_q_mvar.values == 0.2))
    assert(all(net.gen.min_q_mvar.values == [0, 0.1, 0]))
    assert(all(net.gen.min_vm_pu.values == 0.85))
    assert(all(net.gen.max_vm_pu.values == 1.15))
    assert(all(net.gen.vn_kv.values == 0.4))
    assert(all(net.gen.xdss_pu.values == 0.1))
    assert(all(net.gen.rdss_pu.values == 0.1))
    assert(all(net.gen.cos_phi.values == 1.))

def test_create_gens_raise_except():
    net = pp.create_empty_network()
    # standard
    b1 = pp.create_bus(net, 110)
    b2 = pp.create_bus(net, 110)
    b3 = pp.create_bus(net, 110)

    with pytest.raises(UserWarning, match="Cannot attach to buses \{3, 4, 5\}, they does not exist"):
        pp.create_gens(net, buses=[3, 4, 5], p_mw=[0, 0, 1], vm_pu=1., controllable=[True, False, False],
                       max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                       min_vm_pu=0.85, max_vm_pu=1.15, vn_kv=0.4, xdss_pu=0.1, rdss_pu=0.1, cos_phi=1.
                       )

    g = pp.create_gens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], vm_pu=1., controllable=[True, False, False],
                        max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                        min_vm_pu=0.85, max_vm_pu=1.15, vn_kv=0.4, xdss_pu=0.1, rdss_pu=0.1, cos_phi=1.
                        )
    with pytest.raises(UserWarning, match="gens with the ids \[0 1 2\] already exists"):
        pp.create_gens(net, buses=[b1, b2, b3], p_mw=[0, 0, 1], vm_pu=1., controllable=[True, False, False],
                       max_p_mw=0.2, min_p_mw=[0, 0.1, 0],  max_q_mvar=0.2, min_q_mvar=[0, 0.1, 0],
                       min_vm_pu=0.85, max_vm_pu=1.15, vn_kv=0.4, xdss_pu=0.1, rdss_pu=0.1, cos_phi=1.,
                       index=g)

if __name__ == '__main__':
    pytest.main(["test_create.py"])

max_stars_repo_name: jurasofish/pandapower

max_issues_repo_name: jurasofish/pandapower

max_forks_repo_name: jurasofish/pandapower


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

from spcdist.scipy_1d import EntmaxGaussian1D

if __name__ == "__main__":

    ts = np.linspace(-3, 3, num=1000).reshape(-1, 1)

    for alpha in (2, 1.5, 0.9):
        mbg = EntmaxGaussian1D(mu=0, sigma_sq=1, alpha=alpha)
        plt.plot(ts, mbg.pdf(ts), label=fr"$\alpha={alpha}$")
    plt.legend()
    plt.show()

ext: py

reason: import numpy


text:
#%%

import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

#%%
# draw cifar 10 model
for type in ['plain','residual']:
    df_train = pd.DataFrame(columns=['iter','model','error'])
    dataset = 'train'
    for model in ['res20','res32','res44','res56']:
        file_path = '../logs/cifar10/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 50
        else:
            temp['iter'] = temp['iter'] * 354

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_train = df_train.append(temp)

    model = 'res110'
    for i in range(1,6):
        file_path = '../logs/cifar10/residual/'+model+'_residual_'+dataset+'_err'+str(i)+'.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 50
        else:
            temp['iter'] = temp['iter'] * 354

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_train = df_train.append(temp)

    # 过滤非整epoch
    df_train = df_train[df_train['iter'] % 2000 ==0]
    df_train['dataset'] = 'train'
    df_train['error'] = df_train['error']*100

    df_test = pd.DataFrame(columns=['iter','model','error'])
    dataset = 'test'
    for model in ['res20','res32','res44','res56']:
        file_path = '../logs/cifar10/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 50
        else:
            temp['iter'] = temp['iter'] * 354

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_test = df_test.append(temp)

    model = 'res110'
    for i in range(1,6):
        file_path = '../logs/cifar10/residual/'+model+'_residual_'+dataset+'_err'+str(i)+'.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 50
        else:
            temp['iter'] = temp['iter'] * 354

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_test = df_test.append(temp)

    # 过滤非整epoch
    df_test = df_test[df_test['iter'] % 4 ==0]
    df_test['error'] = df_test['error']*100
    df_test['dataset'] = 'val'
    # sns.lineplot(x='iter',y='error',hue='model',style='model',data=df_test, dashes=[0.1,0.])
    # df = df_train.append(df_test)
    df = df_test.append(df_train)

    plt.figure(figsize=(8, 5))
    if type == 'plain':
        plt.ylim(0, 70)
    else:
        plt.ylim(0, 30)
    sns.lineplot(x='iter',y='error',hue='model',style='dataset',data=df)

# %%
# draw tiny_imagenet model
# model_list = ['res18','res34']
model_list = ['res18','res34','res50','res101']
for type in ['plain','residual']:
    df_train = pd.DataFrame(columns=['iter','model','error'])
    dataset = 'train'
    for model in model_list:
        file_path = '../logs/tiny_imagenet/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 35
            pass
        else:
            temp['iter'] = temp['iter'] * 250

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_train = df_train.append(temp)

    if 'res101' in model_list and type == 'residual':
        model = 'res152'
        file_path = '../logs/tiny_imagenet/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 35
            pass
        else:
            temp['iter'] = temp['iter'] * 250

        temp['model'] = model
        df_train = df_train.append(temp)

    # 过滤非整epoch
    df_train = df_train[df_train['iter'] % 100 ==0]
    df_train['dataset'] = 'train'
    df_train['error'] = df_train['error']*100

    df_test = pd.DataFrame(columns=['iter','model','error'])
    dataset = 'test'
    for model in model_list:
        file_path = '../logs/tiny_imagenet/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 35
        else:
            temp['iter'] = temp['iter'] * 250

        if type == 'residual':
            temp['model'] = model
        else:
            temp['model'] = type + model.replace('res','')
        df_test = df_test.append(temp)

    if 'res101' in model_list and type == 'residual':
        model = 'res152'
        file_path = '../logs/tiny_imagenet/'+type+'/'+model+'_'+type+'_'+dataset+'_err.csv'
        temp = pd.read_csv(file_path,names=['iter','unknown','error'])
        temp = temp.drop(columns=['unknown'],axis=1)
        if dataset == 'train':
            temp['iter'] = temp['iter'] * 35
        else:
            temp['iter'] = temp['iter'] * 250

        temp['model'] = model
        df_test = df_test.append(temp)

    # 过滤非整epoch
    df_test = df_test[df_test['iter'] % 4 ==0]
    df_test['error'] = df_test['error']*100
    df_test['dataset'] = 'val'

    df = df_test.append(df_train)

    plt.figure(figsize=(8, 5))
    if type == 'plain':
        plt.ylim(0, 100)
    else:
        plt.ylim(0, 100)
    sns.lineplot(x='iter',y='error',hue='model',style='dataset',data=df)
    plt.show()

# %%


text:
import os
import sys
import time
import logging
import argparse
import numpy as np
from tqdm import tqdm

import torch
import torch.utils.data

from lib.averageMeter import AverageMeters
from lib.logger import colorlogger
from lib.timer import Timers
from lib.averageMeter import AverageMeters
from lib.torch_utils import adjust_learning_rate

from modeling.build_model import Pose2Seg
from datasets.CocoDatasetInfo import CocoDatasetInfo, annToMask
from test import test

NAME = "release_base"

# Set `LOG_DIR` and `SNAPSHOT_DIR`
def setup_logdir():
    timestamp = time.strftime("%Y-%m-%d_%H:%M:%S", time.localtime())
    LOGDIR = os.path.join('logs', '%s_%s'%(NAME, timestamp))
    SNAPSHOTDIR = os.path.join('snapshot', '%s_%s'%(NAME, timestamp))
    if not os.path.exists(LOGDIR):
        os.makedirs(LOGDIR)
    if not os.path.exists(SNAPSHOTDIR):
        os.makedirs(SNAPSHOTDIR)
    return LOGDIR, SNAPSHOTDIR
LOGDIR, SNAPSHOTDIR = setup_logdir()

# Set logging
logger = colorlogger(log_dir=LOGDIR, log_name='train_logs.txt')

# Set Global Timer
timers = Timers()

# Set Global AverageMeter
averMeters = AverageMeters()

def train(model, dataloader, optimizer, epoch, iteration):
    # switch to train mode
    model.train()

    averMeters.clear()
    end = time.time()
    for i, inputs in enumerate(dataloader):
        averMeters['data_time'].update(time.time() - end)
        iteration += 1

        lr = adjust_learning_rate(optimizer, iteration, BASE_LR=0.0002,
                         WARM_UP_FACTOR=1.0/3, WARM_UP_ITERS=1000,
                         STEPS=(0, 14150*15, 14150*20), GAMMA=0.1)

        # forward
        outputs = model(**inputs)

        # loss
        loss = outputs

        # backward
        averMeters['loss'].update(loss.data.item())
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        # measure elapsed time
        averMeters['batch_time'].update(time.time() - end)
        end = time.time()

        if i % 10 == 0:
            logger.info('Epoch: [{0}][{1}/{2}]\t'
                  'Lr: [{3}]\t'
                  'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                  'Data {data_time.val:.3f} ({data_time.avg:.3f})\t'
                  'loss {loss.val:.5f} ({loss.avg:.5f})\t'
                  .format(
                      epoch, i, len(dataloader), lr,
                      batch_time=averMeters['batch_time'], data_time=averMeters['data_time'],
                      loss=averMeters['loss'])
                 )

        if i % 10000 == 0:
            torch.save(model.state_dict(), os.path.join(SNAPSHOTDIR, '%d_%d.pkl'%(epoch,i)))
            torch.save(model.state_dict(), os.path.join(SNAPSHOTDIR, 'last.pkl'))

    return iteration

class Dataset():
    def __init__(self):
        ImageRoot = './data/coco2017/train2017'
        AnnoFile = './data/coco2017/annotations/person_keypoints_train2017_pose2seg.json'
        self.datainfos = CocoDatasetInfo(ImageRoot, AnnoFile, onlyperson=True, loadimg=True)

    def __len__(self):
        return len(self.datainfos)

    def __getitem__(self, idx):
        rawdata = self.datainfos[idx]
        img = rawdata['data']
        image_id = rawdata['id']

        height, width = img.shape[0:2]
        gt_kpts = np.float32(rawdata['gt_keypoints']).transpose(0, 2, 1) # (N, 17, 3)
        gt_segms = rawdata['segms']
        gt_masks = np.array([annToMask(segm, height, width) for segm in gt_segms])

        return {'img': img, 'kpts': gt_kpts, 'masks': gt_masks}

    def collate_fn(self, batch):
        batchimgs = [data['img'] for data in batch]
        batchkpts = [data['kpts'] for data in batch]
        batchmasks = [data['masks'] for data in batch]
        return {'batchimgs': batchimgs, 'batchkpts': batchkpts, 'batchmasks':batchmasks}

if __name__=='__main__':
    logger.info('===========> loading model <===========')
    model = Pose2Seg().cuda()
    #model.init("")
    model.train()

    logger.info('===========> loading data <===========')
    datasetTrain = Dataset()
    dataloaderTrain = torch.utils.data.DataLoader(datasetTrain, batch_size=4, shuffle=True,
                                                   num_workers=4, pin_memory=False,
                                                   collate_fn=datasetTrain.collate_fn)

    logger.info('===========> set optimizer <===========')
    ''' set your optimizer like this. Normally is Adam/SGD. '''
    #optimizer = torch.optim.SGD(model.parameters(), 0.0002, momentum=0.9, weight_decay=0.0005)
    optimizer = torch.optim.Adam(model.parameters(), 0.0002, weight_decay=0.0000)

    iteration = 0
    epoch = 0
    try:
        while iteration < 14150*25:
            logger.info('===========>   training    <===========')
            iteration = train(model, dataloaderTrain, optimizer, epoch, iteration)
            epoch += 1

            logger.info('===========>   testing    <===========')
            test(model, dataset='cocoVal', logger=logger.info)
            test(model, dataset='OCHumanVal', logger=logger.info)

    except (KeyboardInterrupt):
        logger.info('Save ckpt on exception ...')
        torch.save(model.state_dict(), os.path.join(SNAPSHOTDIR, 'interrupt_%d_%d.pkl'%(epoch,iteration)))
        logger.info('Save ckpt done.')


text:
import copy
import os
import sys
import urllib.request

from sure_tosca_client import Configuration, ApiClient, NodeTemplate
from sure_tosca_client.api import default_api
import networkx as nx
import matplotlib.pyplot as plt

class ToscaHelper:

    def __init__(self, sure_tosca_base_url, tosca_template_path):
        self.sure_tosca_base_url = sure_tosca_base_url
        self.tosca_template_path = tosca_template_path
        self.tosca_client = self.init_sure_tosca_client(sure_tosca_base_url)
        self.doc_id = self.upload_tosca_template(tosca_template_path)

    def upload_tosca_template(self, file_path):
        file_id = self.tosca_client.upload_tosca_template(file_path)
        return file_id

    def init_sure_tosca_client(self,sure_tosca_base_path):
        configuration = Configuration()
        configuration.host = sure_tosca_base_path
        api_client = ApiClient(configuration=configuration)
        api = default_api.DefaultApi(api_client=api_client)
        return api

    def get_interface_types(target):
        interface_types = []
        for interface in target.node_template.interfaces:
            interface_types.append(interface)

        return interface_types

    def get_application_nodes(self):
        return self.tosca_client.get_node_templates(self.doc_id, type_name='tosca.nodes.QC.Application')

    def get_deployment_node_pipeline(self):
        nodes_to_deploy = self.get_application_nodes()
        G = nx.DiGraph()
        sorted_nodes = []
        for node in nodes_to_deploy:
            related_nodes = self.tosca_client.get_related_nodes(self.doc_id,node.name)
            for related_node in related_nodes:
                G.add_edge(node.name, related_node.name)
            #     # We need to deploy the docker orchestrator on the VMs not the topology.
            #     # But the topology is directly connected to the orchestrator not the VMs.
            #     # So we explicitly get the VMs
            #     # I don't like this solution but I can't think of something better.
            #     if related_node.node_template.type == 'tosca.nodes.QC.VM.topology':
            #         vms = self.get_vms()
            #         related_node = vms
            #     pair = (related_node, node)
            #     nodes_pairs.append(pair)
        sorted_graph = sorted(G.in_degree, key=lambda x: x[1], reverse=True)
        for node_tuple in sorted_graph:
            node_name = node_tuple[0]
            for node in nodes_to_deploy:
                if node.name == node_name:
                    sorted_nodes.append(node)
        return sorted_nodes

    @classmethod
    def service_is_up(cls, url):
        code = None
        try:
            code = urllib.request.urlopen(url).getcode()
        except Exception as e:
            exc_type, exc_obj, exc_tb = sys.exc_info()
            fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
            if not e.reason and not e.reason.errno and e.code:
                return False
            else:
                return True

        return True

    def get_vms(self):
        return self.tosca_client.get_node_templates(self.doc_id, type_name='tosca.nodes.QC.VM.Compute')

    def set_node(self, updated_node, tosca_template_dict):
        node_templates = tosca_template_dict['topology_template']['node_templates']
        for node_name in node_templates:
            if node_name == updated_node.name:
                node_templates[node_name] = updated_node.node_template.to_dict()
                return tosca_template_dict

def get_interface_types(node):
    interface_type_names = []
    if node.node_template.interfaces:
        for interface in node.node_template.interfaces:
            interface_type_names.append(interface)
        return interface_type_names

max_stars_repo_name: skoulouzis/DRIP

max_issues_repo_name: skoulouzis/DRIP

max_forks_repo_name: skoulouzis/CONF

reason: import networkx


text:
"""
Base class for models of periodic sources
"""

import numpy as np
from .findpeaks import peakdetect

class PeriodicModeler(object):
    """Base class for periodic modeling"""
    def __init__(self, *args, **kwargs):
        raise NotImplementedError("PeriodicModeler")

    def fit(self, t, y, dy=None, filts=None):
        """Provide data for the fit"""
        raise NotImplementedError()

    def predict(self, t, filts=None, period=None):
        """Predict the best-fit model at t for the given period"""
        raise NotImplementedError()

    def score(self, period):
        """Compute the score for a period or array of periods"""
        raise NotImplementedError()

    def period_search(self, pmin, pmax, resolution=1e4, nperiods=1,
                      return_scores=False, linspace=False):
        """Find the best period for the model"""
        # TODO: implement using the score() function provided by subclasses

        if linspace:
            periods = np.linspace(pmin, pmax, resolution)
        else:
            periods = np.logspace(np.log10(pmin),
                                np.log10(pmax),
                                resolution)

        scores = self.score(periods)
        maxes,mins = peakdetect(scores, periods, lookahead=5)

        maxes = np.array(maxes)
        inds = np.argsort(maxes[:,1])

        pks,hts = maxes[inds,0][-nperiods:][::-1],maxes[inds,1][-nperiods:][::-1]

        if len(pks)==1:
            pks = pks[0]
            hts = hts[0]

        if return_scores:
            return pks, hts
        else:
            return pks

max_stars_repo_name: stephtdouglas/periodogram

max_issues_repo_name: stephtdouglas/periodogram

max_forks_repo_name: stephtdouglas/periodogram


text:
"""Tests for the UnitQuaternion class.
"""

import numpy as np
import pytest

from walle.core import Quaternion, UnitQuaternion

class TestUnitQuaternion(object):
  def axis_angle_vector(self, deg):
    theta = np.deg2rad(deg)
    unit_vec = np.array([0, 0, 1])
    return unit_vec, theta

  def rotation_vector(self, deg):
    unit_vec, theta = self.axis_angle_vector(deg)
    return theta * unit_vec

  def quat_from_axang(self, deg):
    unit_vec, theta = self.axis_angle_vector(deg)
    s = np.cos(theta / 2)
    v = unit_vec * np.sin(theta / 2)
    return UnitQuaternion(s, v)

  def test_inv(self):
    """Tests the inverse of the unit quaternion.
    """
    quat = self.quat_from_axang(90)
    expected = quat.conjugate()
    actual = quat.inv()
    assert actual == expected

  def test_rotation_equiv(self):
    """Tests that Rodrigues' formula is equivalent to conjugation.
    """
    unit_quat = self.quat_from_axang(90)
    quat = Quaternion(unit_quat.array)
    pos_vec = np.array([1, 2, 3])
    actual = Quaternion(unit_quat * pos_vec)
    expected = quat * pos_vec * quat.inv()
    assert actual == expected

  def test_rotation_i2j(self):
    """Tests the rotation of unit vector `i` by 90 degrees.
    """
    quat = self.quat_from_axang(90)
    unit_i = np.array([1, 0, 0])
    actual = (quat * unit_i)
    expected = np.array([0, 1, 0])
    assert np.allclose(actual, expected)

  def test_rotation_i2negi(self):
    """Tests the rotation of unit vector `i` by 180 degrees.
    """
    quat = self.quat_from_axang(180)
    unit_i = np.array([1, 0, 0])
    actual = (quat * unit_i)
    expected = np.array([-1, 0, 0])
    assert np.allclose(actual, expected)

  def test_slerp_zero(self):
    """Tests slerp with `t=0`.
    """
    quat1 = self.quat_from_axang(90)
    quat2 = self.quat_from_axang(180)
    expected = quat1
    actual = quat1.slerp(quat2, 0)
    assert actual == expected

  def test_slerp_one(self):
    """Tests slerp with `t=1`.
    """
    quat1 = self.quat_from_axang(90)
    quat2 = self.quat_from_axang(180)
    expected = quat2
    actual = quat1.slerp(quat2, 1.)
    assert actual == expected

  def test_slerp_rand(self):
    """Tests slerp with a random `t`.

    Ref: https://en.wikipedia.org/wiki/Slerp#Geometric_Slerp
    """
    quat1 = self.quat_from_axang(90)
    quat2 = self.quat_from_axang(180)
    t = np.random.random()
    omega = quat1.angle(quat2)
    coeff1 = np.sin((1 - t) * omega) / np.sin(omega)
    coeff2 = np.sin(t * omega) / np.sin(omega)
    expected = UnitQuaternion((coeff1 * quat1.array) + (coeff2 * quat2.array))
    actual = quat1.slerp(quat2, t)
    assert actual == expected

  def test_rotation_n_times(self):
    """Tests the rotation of unit vector `i` by 90 degrees 4 times.

    We should essentially recover `i`.
    """
    quat = self.quat_from_axang(90)
    unit_i = np.array([1, 0, 0])
    actual = (quat**4 * unit_i)
    expected = unit_i
    assert np.allclose(actual, expected)

  def test_quaternion_to_axang(self):
    """Tests quaternion to axis-angle conversion.
    """
    quat = self.quat_from_axang(90)
    axis_actual, angle_actual = quat.axis_angle
    axis_expected, angle_expected = self.axis_angle_vector(90)
    assert np.isclose(angle_actual, angle_expected) and np.allclose(axis_actual, axis_expected)

  def test_quaternion_to_rot_vec(self):
    """Tests quaternion to rotation vector conversion.
    """
    quat = self.quat_from_axang(90)
    expected = self.rotation_vector(90)
    actual = quat.rot_vec
    assert np.allclose(actual, expected)

max_stars_repo_name: kaixin-bai/walle

max_issues_repo_name: kaixin-bai/walle

max_forks_repo_name: kaixin-bai/walle


text:
import numpy as np
import tensorflow as tf
from tensorflow.contrib.framework import add_arg_scope

from tfsnippet.ops import (assert_rank, assert_scalar_equal, flatten_to_ndims,
                           unflatten_from_ndims)
from tfsnippet.utils import (validate_positive_int_arg, ParamSpec,
                             is_tensor_object, assert_deps,
                             get_shape,
                             add_name_and_scope_arg_doc, model_variable,
                             maybe_check_numerics, maybe_add_histogram, InputSpec,
                             get_static_shape, validate_enum_arg, validate_int_tuple_arg)

__all__ = [
    'conv1d', 'deconv1d', 'validate_conv1d_input', 'get_deconv_output_length', 'batch_norm_1d'
]

@add_arg_scope
def batch_norm_1d(input, channels_last=True, training=False, name=None,
                  scope=None):
    """
    Apply batch normalization on 1D convolutional layer.

    Args:
        input (tf.Tensor): The input tensor.
        channels_last (bool): Whether or not the channel dimension is at last?
        training (bool or tf.Tensor): Whether or not the model is under
            training stage?

    Returns:
        tf.Tensor: The normalized tensor.
    """
    with tf.variable_scope(scope, default_name=name or 'batch_norm_1d'):
        input, s1, s2 = flatten_to_ndims(input, ndims=3)
        output = tf.layers.batch_normalization(
            input,
            axis=-1 if channels_last else -2,
            training=training,
            name='norm'
        )
        output = unflatten_from_ndims(output, s1, s2)
        return output

def validate_conv1d_input(input, channels_last, arg_name='input'):
    """
    Validate the input for 1-d convolution.
    Args:
        input: The input tensor, must be at least 3-d.
        channels_last (bool): Whether or not the last dimension is the
            channels dimension? (i.e., the data format is (batch, length, channels))
        arg_name (str): Name of the input argument.
    Returns:
        (tf.Tensor, int, str): The validated input tensor, the number of input
            channels, and the data format.
    """
    if channels_last:
        channel_axis = -1
        input_spec = InputSpec(shape=('...', '?', '?', '*'))
        data_format = "NWC"
    else:
        channel_axis = -2
        input_spec = InputSpec(shape=('...', '?', '*', '?'))
        data_format = "NCW"
    input = input_spec.validate(arg_name, input)
    input_shape = get_static_shape(input)
    in_channels = input_shape[channel_axis]

    return input, in_channels, data_format

def get_deconv_output_length(input_length, kernel_size, strides, padding):
    """
    Get the output length of deconvolution at a specific dimension.
    Args:
        input_length: Input tensor length.
        kernel_size: The size of the kernel.
        strides: The stride of convolution.
        padding: One of {"same", "valid"}, case in-sensitive
    Returns:
        int: The output length of deconvolution.
    """
    padding = validate_enum_arg(
        'padding', str(padding).upper(), ['SAME', 'VALID'])
    output_length = input_length * strides
    if padding == 'VALID':
        output_length += max(kernel_size - strides, 0)
    return output_length

@add_arg_scope
@add_name_and_scope_arg_doc
def conv1d(input,
           out_channels,
           kernel_size,
           strides=1,
           dilations=1,
           padding='same',
           channels_last=True,
           activation_fn=None,
           normalizer_fn=None,
           gated=False,
           gate_sigmoid_bias=2.,
           kernel=None,
           kernel_mask=None,
           kernel_initializer=None,
           kernel_regularizer=None,
           kernel_constraint=None,
           use_bias=None,
           bias=None,
           bias_initializer=tf.zeros_initializer(),
           bias_regularizer=None,
           bias_constraint=None,
           trainable=True,
           name=None,
           scope=None):
    """
    1D convolutional layer.
    Args:
        input (Tensor): The input tensor, at least 3-d.
        out_channels (int): The channel numbers of the output.
        kernel_size (int or tuple(int,)): Kernel size over spatial dimensions.
        strides (int): Strides over spatial dimensions.
        dilations (int): The dilation factor over spatial dimensions.
        padding: One of {"valid", "same"}, case in-sensitive.
        channels_last (bool): Whether or not the channel axis is the last
            axis in `input`? (i.e., the data format is "NWC")
        activation_fn: The activation function.
        normalizer_fn: The normalizer function.
        gated (bool): Whether or not to use gate on output?
            `output = activation_fn(output) * sigmoid(gate)`.
        gate_sigmoid_bias (Tensor): The bias added to `gate` before applying
            the `sigmoid` activation.
        kernel (Tensor): Instead of creating a new variable, use this tensor.
        kernel_mask (Tensor): If specified, multiply this mask onto `kernel`,
            i.e., the actual kernel to use will be `kernel * kernel_mask`.
        kernel_initializer: The initializer for `kernel`.
            Would be ``default_kernel_initializer(...)`` if not specified.
        kernel_regularizer: The regularizer for `kernel`.
        kernel_constraint: The constraint for `kernel`.
        use_bias (bool or None): Whether or not to use `bias`?
            If :obj:`True`, will always use bias.
            If :obj:`None`, will use bias only if `normalizer_fn` is not given.
            If :obj:`False`, will never use bias.
            Default is :obj:`None`.
        bias (Tensor): Instead of creating a new variable, use this tensor.
        bias_initializer: The initializer for `bias`.
        bias_regularizer: The regularizer for `bias`.
        bias_constraint: The constraint for `bias`.
        trainable (bool): Whether or not the parameters are trainable?
    Returns:
        tf.Tensor: The output tensor.
    """
    if not channels_last:
        raise ValueError('Currently only channels_last=True is supported.')
    input, in_channels, data_format = \
        validate_conv1d_input(input, channels_last)
    out_channels = validate_positive_int_arg('out_channels', out_channels)
    dtype = input.dtype.base_dtype
    if gated:
        out_channels *= 2

    # check functional arguments
    padding = validate_enum_arg(
        'padding', str(padding).upper(), ['VALID', 'SAME'])
    dilations = validate_positive_int_arg('dilations', dilations)
    strides = validate_positive_int_arg('strides', strides)

    if dilations > 1 and not channels_last:
        raise ValueError('`channels_last` == False is incompatible with '
                         '`dilations` > 1.')

    if strides > 1 and dilations > 1:
        raise ValueError('`strides` > 1 is incompatible with `dilations` > 1.')

    if use_bias is None:
        use_bias = normalizer_fn is None

    # get the specification of outputs and parameters
    kernel_size = validate_int_tuple_arg('kernel_size', kernel_size)
    kernel_shape = kernel_size + (in_channels, out_channels)
    bias_shape = (out_channels,)

    # validate the parameters
    if kernel is not None:
        kernel_spec = ParamSpec(shape=kernel_shape, dtype=dtype)
        kernel = kernel_spec.validate('kernel', kernel)
    if kernel_mask is not None:
        kernel_mask_spec = InputSpec(dtype=dtype)
        kernel_mask = kernel_mask_spec.validate('kernel_mask', kernel_mask)
    if kernel_initializer is None:
        kernel_initializer = tf.glorot_normal_initializer()
    if bias is not None:
        bias_spec = ParamSpec(shape=bias_shape, dtype=dtype)
        bias = bias_spec.validate('bias', bias)

    # the main part of the conv1d layer
    with tf.variable_scope(scope, default_name=name or 'conv1d'):
        c_axis = -1 if channels_last else -2

        # create the variables
        if kernel is None:
            kernel = model_variable(
                'kernel',
                shape=kernel_shape,
                dtype=dtype,
                initializer=kernel_initializer,
                regularizer=kernel_regularizer,
                constraint=kernel_constraint,
                trainable=trainable
            )

        if kernel_mask is not None:
            kernel = kernel * kernel_mask

        maybe_add_histogram(kernel, 'kernel')
        kernel = maybe_check_numerics(kernel, 'kernel')

        if use_bias and bias is None:
            bias = model_variable(
                'bias',
                shape=bias_shape,
                initializer=bias_initializer,
                regularizer=bias_regularizer,
                constraint=bias_constraint,
                trainable=trainable
            )
            maybe_add_histogram(bias, 'bias')
            bias = maybe_check_numerics(bias, 'bias')

        # flatten to 3d
        output, s1, s2 = flatten_to_ndims(input, 3)

        # do convolution
        if dilations > 1:
            output = tf.nn.convolution(
                input=output,
                filter=kernel,
                dilation_rate=(dilations,),
                padding=padding,
                data_format=data_format
            )
        else:
            output = tf.nn.conv1d(
                value=output,
                filters=kernel,
                stride=strides,
                padding=padding,
                data_format=data_format
            )

        # add bias
        if use_bias:
            output = tf.add(output, bias)

        # apply the normalization function if specified
        if normalizer_fn is not None:
            output = normalizer_fn(output)

        # split into halves if gated
        if gated:
            output, gate = tf.split(output, 2, axis=c_axis)

        # apply the activation function if specified
        if activation_fn is not None:
            output = activation_fn(output)

        # apply the gate if required
        if gated:
            if gate_sigmoid_bias is None:
                gate_sigmoid_bias = model_variable(
                    'gate_sigmoid_bias',
                    shape=bias_shape,
                    initializer=bias_initializer,
                    regularizer=bias_regularizer,
                    constraint=bias_constraint,
                    trainable=trainable
                )
                maybe_add_histogram(gate_sigmoid_bias, 'gate_sigmoid_bias')
                gate_sigmoid_bias = maybe_check_numerics(gate_sigmoid_bias, 'gate_sigmoid_bias')
            output = output * tf.sigmoid(gate + gate_sigmoid_bias, name='gate')

        # unflatten back to original shape
        output = unflatten_from_ndims(output, s1, s2)

        maybe_add_histogram(output, 'output')
        output = maybe_check_numerics(output, 'output')
    return output

@add_arg_scope
@add_name_and_scope_arg_doc
def deconv1d(input,
             out_channels,
             kernel_size,
             strides=1,
             padding='same',
             channels_last=True,
             output_shape=None,
             activation_fn=None,
             normalizer_fn=None,
             gated=False,
             gate_sigmoid_bias=2.,
             kernel=None,
             kernel_initializer=None,
             kernel_regularizer=None,
             kernel_constraint=None,
             use_bias=None,
             bias=None,
             bias_initializer=tf.zeros_initializer(),
             bias_regularizer=None,
             bias_constraint=None,
             trainable=True,
             name=None,
             scope=None):
    """
    1D deconvolutional layer.
    Args:
        input (Tensor): The input tensor, at least 3-d.
        out_channels (int): The channel numbers of the deconvolution output.
        kernel_size (int or tuple(int,)): Kernel size over spatial dimensions.
        strides (int): Strides over spatial dimensions.
        padding: One of {"valid", "same"}, case in-sensitive.
        channels_last (bool): Whether or not the channel axis is the last
            axis in `input`? (i.e., the data format is "NWC")
        output_shape: If specified, use this as the shape of the
            deconvolution output; otherwise compute the size of each dimension
            by::
                output_size = input_size * strides
                if padding == 'valid':
                    output_size += max(kernel_size - strides, 0)
        activation_fn: The activation function.
        normalizer_fn: The normalizer function.
        gated (bool): Whether or not to use gate on output?
            `output = activation_fn(output) * sigmoid(gate)`.
        gate_sigmoid_bias (Tensor): The bias added to `gate` before applying
            the `sigmoid` activation.
        kernel (Tensor): Instead of creating a new variable, use this tensor.
        kernel_initializer: The initializer for `kernel`.
            Would be ``default_kernel_initializer(...)`` if not specified.
        kernel_regularizer: The regularizer for `kernel`.
        kernel_constraint: The constraint for `kernel`.
        use_bias (bool or None): Whether or not to use `bias`?
            If :obj:`True`, will always use bias.
            If :obj:`None`, will use bias only if `normalizer_fn` is not given.
            If :obj:`False`, will never use bias.
            Default is :obj:`None`.
        bias (Tensor): Instead of creating a new variable, use this tensor.
        bias_initializer: The initializer for `bias`.
        bias_regularizer: The regularizer for `bias`.
        bias_constraint: The constraint for `bias`.
        trainable (bool): Whether or not the parameters are trainable?
    Returns:
        tf.Tensor: The output tensor.
    """
    if not channels_last:
        raise ValueError('Currently only channels_last=True is supported.')
    input, in_channels, data_format = \
        validate_conv1d_input(input, channels_last)
    out_channels = validate_positive_int_arg('out_channels', out_channels)
    dtype = input.dtype.base_dtype
    if gated:
        out_channels *= 2

    # check functional arguments
    padding = validate_enum_arg(
        'padding', str(padding).upper(), ['VALID', 'SAME'])
    strides = validate_positive_int_arg('strides', strides)

    if use_bias is None:
        use_bias = normalizer_fn is None

    # get the specification of outputs and parameters
    kernel_size = validate_int_tuple_arg('kernel_size', kernel_size)
    kernel_shape = kernel_size + (out_channels, in_channels)
    bias_shape = (out_channels,)

    given_w = None
    given_output_shape = output_shape

    if is_tensor_object(given_output_shape):
        given_output_shape = tf.convert_to_tensor(given_output_shape)
    elif given_output_shape is not None:
        given_w = given_output_shape

    # validate the parameters
    if kernel is not None:
        kernel_spec = ParamSpec(shape=kernel_shape, dtype=dtype)
        kernel = kernel_spec.validate('kernel', kernel)
    if kernel_initializer is None:
        kernel_initializer = tf.glorot_normal_initializer()
    if bias is not None:
        bias_spec = ParamSpec(shape=bias_shape, dtype=dtype)
        bias = bias_spec.validate('bias', bias)

    # the main part of the conv2d layer
    with tf.variable_scope(scope, default_name=name or 'deconv1d'):
        with tf.name_scope('output_shape'):
            # detect the input shape and axis arrangements
            input_shape = get_static_shape(input)
            if channels_last:
                c_axis, w_axis = -1, -2
            else:
                c_axis, w_axis = -2, -1

            output_shape = [None, None, None]
            output_shape[c_axis] = out_channels
            if given_output_shape is None:
                if input_shape[w_axis] is not None:
                    output_shape[w_axis] = get_deconv_output_length(
                        input_shape[w_axis], kernel_shape[0], strides[0],
                        padding
                    )
            else:
                if not is_tensor_object(given_output_shape):
                    output_shape[w_axis] = given_w

            # infer the batch shape in 3-d
            batch_shape = input_shape[:-2]
            if None not in batch_shape:
                output_shape[0] = int(np.prod(batch_shape))

            # now the static output shape is ready
            output_static_shape = tf.TensorShape(output_shape)

            # prepare for the dynamic batch shape
            if output_shape[0] is None:
                output_shape[0] = tf.reduce_prod(get_shape(input)[:-2])

            # prepare for the dynamic spatial dimensions
            if output_shape[w_axis] is None:
                if given_output_shape is None:
                    input_shape = get_shape(input)
                    if output_shape[w_axis] is None:
                        output_shape[w_axis] = get_deconv_output_length(
                            input_shape[w_axis], kernel_shape[0],
                            strides[0], padding
                        )
                else:
                    assert(is_tensor_object(given_output_shape))
                    with assert_deps([
                        assert_rank(given_output_shape, 1),
                        assert_scalar_equal(
                            tf.size(given_output_shape), 1)
                    ]):
                        output_shape[w_axis] = given_output_shape[0]

            # compose the final dynamic shape
            if any(is_tensor_object(s) for s in output_shape):
                output_shape = tf.stack(output_shape)
            else:
                output_shape = tuple(output_shape)

        # create the variables
        if kernel is None:
            kernel = model_variable(
                'kernel',
                shape=kernel_shape,
                dtype=dtype,
                initializer=kernel_initializer,
                regularizer=kernel_regularizer,
                constraint=kernel_constraint,
                trainable=trainable
            )

        maybe_add_histogram(kernel, 'kernel')
        kernel = maybe_check_numerics(kernel, 'kernel')

        if use_bias and bias is None:
            bias = model_variable(
                'bias',
                shape=bias_shape,
                initializer=bias_initializer,
                regularizer=bias_regularizer,
                constraint=bias_constraint,
                trainable=trainable
            )
            maybe_add_histogram(bias, 'bias')
            bias = maybe_check_numerics(bias, 'bias')

        # flatten to 3d
        output, s1, s2 = flatten_to_ndims(input, 3)

        # do convolution or deconvolution
        output = tf.contrib.nn.conv1d_transpose(
            value=output,
            filter=kernel,
            output_shape=output_shape,
            stride=strides,
            padding=padding,
            data_format=data_format
        )
        if output_static_shape is not None:
            output.set_shape(output_static_shape)

        # add bias
        if use_bias:
            output = tf.add(output, bias)

        # apply the normalization function if specified
        if normalizer_fn is not None:
            output = normalizer_fn(output)

        # split into halves if gated
        if gated:
            output, gate = tf.split(output, 2, axis=c_axis)

        # apply the activation function if specified
        if activation_fn is not None:
            output = activation_fn(output)

        # apply the gate if required
        if gated:
            if gate_sigmoid_bias is None:
                gate_sigmoid_bias = model_variable(
                    'gate_sigmoid_bias',
                    shape=bias_shape,
                    initializer=bias_initializer,
                    regularizer=bias_regularizer,
                    constraint=bias_constraint,
                    trainable=trainable
                )
                maybe_add_histogram(gate_sigmoid_bias, 'gate_sigmoid_bias')
                gate_sigmoid_bias = maybe_check_numerics(gate_sigmoid_bias, 'gate_sigmoid_bias')
            output = output * tf.sigmoid(gate + gate_sigmoid_bias, name='gate')

        # unflatten back to original shape
        output = unflatten_from_ndims(output, s1, s2)

        maybe_add_histogram(output, 'output')
        output = maybe_check_numerics(output, 'output')

    return output

max_stars_repo_name: hhhercules/InterFusion

max_issues_repo_name: hhhercules/InterFusion

max_forks_repo_name: hhhercules/InterFusion


text:
import numpy as np
import pyomo.environ as pyEnv

path = "/home/tito/cplex/bin/x86-64_linux/cplex"
#            Cliente 1 Cliente 2
#Fabrica 1    c11        c12
#Fabrica 2    c21        c22
#Fabrica 3    c31        c32

#Custo de transporte da fabrica i pro cliente j
cij = [[162,247],
       [117,193],
       [131,185]]

capacidades = [1000,1500,1200]
demandas = [2300,1400]

m = len(capacidades)
n = len(demandas)

##-------------------------DECLARACAO DO MODELO E SEUS PARAMETROS--------------------##
#Modelo
modelo = pyEnv.ConcreteModel()
#Indice para as fabricas
modelo.I = pyEnv.RangeSet(m)
#Indice para os clientes
modelo.J = pyEnv.RangeSet(n)

#Variaveis de decisao xij
modelo.x=pyEnv.Var(modelo.I,modelo.J, within=pyEnv.NonNegativeReals)

#Custo de transporte da fabrica i pro cliente j
modelo.c = pyEnv.Param(modelo.I, modelo.J,initialize=lambda modelo, i, j: cij[i-1][j-1])
#Capacidade de cada fabrica
modelo.b = pyEnv.Param(modelo.I, initialize = lambda modelo,i: capacidades[i-1])
#Demanda de cada cliente
modelo.d = pyEnv.Param(modelo.J, initialize = lambda modelo,j: demandas[j-1])

##-------------------------DECLARACAO DA FUNCAO OBJETIVO E RESTRICOES--------------------##

def func_objetivo(modelo):
    return sum(modelo.x[i,j] * modelo.c[i,j] for i in modelo.I for j in modelo.J)

modelo.objetivo = pyEnv.Objective(rule=func_objetivo,sense=pyEnv.minimize)

##------------------------------------------------------##
#Cada fabrica nao exceder sua capacidade

def rule_rest1(modelo,i):
    return sum(modelo.x[i,j] for j in modelo.J) <= modelo.b[i]

modelo.rest1 = pyEnv.Constraint(modelo.I,rule=rule_rest1)

##------------------------------------------------------##
#Cada cliente ter sua demanda atendida

def rule_rest2(modelo,j):
    return sum(modelo.x[i,j] for i in modelo.I) >= modelo.d[j]

modelo.rest2 = pyEnv.Constraint(modelo.J,rule=rule_rest2)

##-------------------------------RESOLUCAO DO MODELO-----------------------------------##
solver = pyEnv.SolverFactory('cplex',executable=path)
resultado = solver.solve(modelo,tee = True)
print('\n-------------------------\n')
#modelo.pprint()
modelo.objetivo()
print(resultado)

##-------------------------PRINT DAS VARIAVEIS DE DECISAO ESCOLHIDAS--------------------##
##(cliente,facilidade)
lista = list(modelo.x.keys())
for i in lista:
    print(i,'--', modelo.x[i]())

max_stars_repo_name: hpaucar/autonomous-system-repo

max_issues_repo_name: hpaucar/autonomous-system-repo

max_forks_repo_name: hpaucar/autonomous-system-repo

reason: import numpy,import pyomo


text:
# Copyright 2021 The NetKet Authors - All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from functools import partial
import warnings
from textwrap import dedent

import jax

import netket.jax as nkjax

@partial(jax.jit, static_argnums=(0, 4, 5))
def _choose_jacobian_mode(apply_fun, pars, model_state, samples, mode, holomorphic):
    homogeneous_vars = nkjax.tree_ishomogeneous(pars)

    if holomorphic is True:
        if not homogeneous_vars:
            warnings.warn(
                dedent(
                    """The ansatz has non homogeneous variables, which might not behave well with the
                       holomorhic implemnetation.
                       Use `holomorphic=False` or mode='complex' for more accurate results but
                       lower performance.
                    """
                )
            )
        mode = "holomorphic"
    else:
        leaf_iscomplex = nkjax.tree_leaf_iscomplex(pars)
        complex_output = nkjax.is_complex(
            jax.eval_shape(
                apply_fun,
                {"params": pars, **model_state},
                samples.reshape(-1, samples.shape[-1]),
            )
        )

        if complex_output:
            if leaf_iscomplex:
                if holomorphic is None:
                    warnings.warn(
                        dedent(
                            """
                                Complex-to-Complex model detected. Defaulting to `holomorphic=False` for
                                the implementation of QGTJacobianDense.
                                If your model is holomorphic, specify `holomorphic=True` to use a more
                                performant implementation.
                                To suppress this warning specify `holomorphic`.
                                """
                        ),
                        UserWarning,
                    )
                mode = "complex"
            else:
                mode = "complex"
        else:
            mode = "real"

    if mode == "real":
        return 0
    elif mode == "complex":
        return 1
    elif mode == "holomorphic":
        return 2
    else:
        raise ValueError(f"unknown mode {mode}")

def choose_jacobian_mode(afun, pars, state, samples, *, mode, holomorphic):
    """
    Select an implementation of Jacobian
    """
    i = _choose_jacobian_mode(afun, pars, state, samples, mode, holomorphic).item()
    if i == 0:
        return "real"
    elif i == 1:
        return "complex"
    elif i == 2:
        return "holomorphic"
    else:
        raise ValueError(f"unknown mode {i}")

max_stars_repo_name: gpescia/MyNetKet

max_issues_repo_name: gpescia/MyNetKet

max_forks_repo_name: gpescia/MyNetKet

reason: import jax


text:
# coding: utf-8

from astropy import convolution as astro_conv
from scipy.stats import vonmises, expon
from scipy import signal
from scipy import ndimage
from scipy.special import iv
from math import pi#, pow #this doesn't support elementwise operations
from numpy import zeros_like, linspace, roll, exp, sqrt, where
import psrchive as psr

# TODO make it possible to set nbin or get it from data
nbin = 1024

def sum_of_gaussians(sigma, amp, centre, nbin):
    ''' Produces a sum of Gaussians on a periodic domain.

        sigma, amp, centre are lists of equal length
        nbin is a positive integer

        returns two arrays, one of x values and one of the model'''
    kappa = []
    mu = []
    for i in range(len(sigma)):
        kappa.append(1.0 / (sigma[i]**2) * 2.35428 * 2.35428 * nbin**2 / (4 * pi**2))
        mu.append(centre[i]/nbin)

    x = linspace(0.0, 1.0, nbin)
    xbins = linspace(1, nbin, num=nbin)

    data = zeros_like(xbins)

    for i in range(len(sigma)):
        model = vonmises.pdf(x*2*pi-pi, kappa[i], loc=mu[i]*2*pi - pi, scale=1)
        model *= amp[i] / max(model)
        data += model

    return xbins, data

def filteredGaussians(sigma, amp, centre, nbin, _filter):
    '''
    params:'''
    xbins = linspace(1,nbin,nbin)
    return signal.convolve( sum_of_gaussians(sigma, amp, centre, nbin), _filter, mode='full') / sum(_filter)

def truncateScreen(xbins, _filter, cutoff):
    return where(xbins < cutoff, _filter, [0]*len(xbins))

def infiniteScreenFilter(xbins, tau):
    '''
    produces a PBF for a infinite screen halfway between obs and source
    g_{thin}=\exp(-t/\tau_{s}) (Rusul et al. 2012, original Cronyn 1970, Williamson 1972)
    Based on Cronyn, W. M. 1970, Science 1968, 1453
    params:
    '''
    return expon.pdf(xbins, scale=tau)

def infiniteScreenFilterTruncated(xbins, tau, cutoff):
    return truncateScreen(xbins, infiniteScreenFilter(xbins, tau), cutoff)

def thickScreenNear(xbins, tau):
    '''
    TODO: Is this screan near observer or source? Scalings are different, but are they ultimately just hidden away in tau?
    produces a PBF for a infinite screen halfway between obs and source
    g_{thin}=g_{thick}=(\frac{\pi\tau_{s}}{4t^3})^{1/2}\exp(\frac{-\pi^2\tau_{s}}{16t})
    (Rusul et al. 2012, original Cronyn 1970, Williamson 1972)
    params:
    '''
    xzero = xbins[0]
    xbins[0] = 1.0
    filter_g  = sqrt(pi * tau / 4. / pow(xbins, 3)) * exp(-pow(pi, 2) * tau / 16. / xbins )
    filter_g[0] = 0.0
    xbins[0] = xzero
    return filter_g

def extendedScreen(xbins, tau):
    '''
    A screen spreading uniformly between the observer and the source.
    g_{extend}=(\frac{\pi^5\tau_{s}^3}{8t^5})^{1/2}\exp(\frac{-\pi^2\tau_{s}}{4t}) (Rusul et al. 2012, originally, Williamson 1972)
    '''
    xzero = xbins[0]
    xbins[0] = 1.0
    filter_g = sqrt(pow(pi, 5) * pow(tau, 3) / 8. / pow(xbins, 5)) * exp(-pow(pi, 2) * tau / 4. / xbins )
    filter_g[0] = 0.0
    xbins[0] = xzero
    return filter_g

def anisotropicScreen(xbins, tau_x, tau_y):
    '''
    A screen with different scattering properties in two dimensions:
    g_{anisotropic} = \frac{1}{\sqrt{\tau_x \tau_y}} \exp\left(-t/2\left(\frac{1}{\tau_x}+\frac{1}{\tau_y}\right)\right) I\left(0, t/2\left(\frac{1}{\tau_x}-\frac{1}{\tau_y}\right)\right)
    where I is the modified Bessel function of the first kind.

    From Marisa Geyer's talk at PWG. She calculated it using Jim Cordes' eq. 9 from his 2001 paper "Anomalous radio-wave scattering..."
    '''
    exp_ = exp(-xbins/2. * (1. / tau_x + 1. / tau_y))
    Iv = iv(0, xbins/2. * ( 1. / tau_x - 1. / tau_y))
    scale =  1. / sqrt(tau_x * tau_y)
    return scale * exp_ * Iv

def anisotropicScreenTruncated(xbins, tau_x, tau_y, cutoff):
    '''
    '''
    return truncateScreen(xbins, anisotropicScreen(xbins, tau_x, tau_y), cutoff)

def scatterProfileReturnAsIs(profile, filter_, mode, boundary):
    '''
    This function returns a profile convolved with the filter.
    It will return the result exactly as returned by the convolving function.
    This may be different than what is expected, e.g., for signal.convolve
    the length of output will be different than for the input
    '''
    if mode == 'signal':
        filtered = signal.convolve( data, filter_, mode='full') / sum(filter_)
    elif mode == 'astropy':
        filtered = astro_conv.convolve_fft(profile, roll(filter_,
            len(profile)/2), boundary=boundary)
    else:
        raise ValueError("Mode " + mode + " is not supported. Only 'astropy' and 'signal' are available")
    return filtered

def scatterProfile(profile, filter_, mode='astropy', boundary='wrap'):
    '''
    This function returns profile convolved with filter.
    It ensures that the returned profile is the same length as input data
    '''
    return scatterProfileReturnAsIs(profile, filter_, mode, boundary)[0:nbin] #range needed if mode is signal

def getTemplateFromData(input, smooth=False):
    ''' params:
    input: File used as a template
    smooth: Apply Savitzky-Golay filter to the data
    '''
    ar = psr.Archive_load(input)
    ar.remove_baseline()
    ar.pscrunch()
    ar.fscrunch()
    template = ar.get_data().squeeze()
    if smooth:
        from scipy.signal import savgol_filter
        template_sg = savgol_filter(template, 17, 5)
        # TODO:
        # J2051 specific 550 and 180 below is the off-pulse region I want to blank out
        # for simulations
        template_sg[:180] = 0.0
        template_sg[550:] = 0.0
        # ensure smooth join with the zeroed part
        template = savgol_filter(template_sg, 17, 5)

    return template

def scattered_gaussian_residual(params, data=None, errs=None):
    #prms = np.array([param.value for param in params.itervalues()])
    return (data - scattered_gaussian(params, nbin)) / errs

max_stars_repo_name: sosl/scattsim

max_issues_repo_name: sosl/scattsim

max_forks_repo_name: sosl/scattsim

reason: from numpy,from scipy,from astropy