text:
from datetime import (
    datetime,
    timedelta,
)

import numpy as np
import pytest

from pandas.errors import UnsupportedFunctionCall

from pandas import (
    DataFrame,
    DatetimeIndex,
    MultiIndex,
    Series,
    Timedelta,
    Timestamp,
    date_range,
    period_range,
    to_datetime,
    to_timedelta,
)
import pandas._testing as tm
from pandas.api.indexers import BaseIndexer
from pandas.core.window import Rolling

def test_doc_string():

    df = DataFrame({"B": [0, 1, 2, np.nan, 4]})
    df
    df.rolling(2).sum()
    df.rolling(2, min_periods=1).sum()

def test_constructor(frame_or_series):
    # GH 12669

    c = frame_or_series(range(5)).rolling

    # valid
    c(0)
    c(window=2)
    c(window=2, min_periods=1)
    c(window=2, min_periods=1, center=True)
    c(window=2, min_periods=1, center=False)

    # GH 13383

    msg = "window must be an integer 0 or greater"

    with pytest.raises(ValueError, match=msg):
        c(-1)

@pytest.mark.parametrize("w", [2.0, "foo", np.array([2])])
def test_invalid_constructor(frame_or_series, w):
    # not valid

    c = frame_or_series(range(5)).rolling

    msg = (
        "window must be an integer|"
        "passed window foo is not compatible with a datetimelike index"
    )
    with pytest.raises(ValueError, match=msg):
        c(window=w)

    msg = "min_periods must be an integer"
    with pytest.raises(ValueError, match=msg):
        c(window=2, min_periods=w)

    msg = "center must be a boolean"
    with pytest.raises(ValueError, match=msg):
        c(window=2, min_periods=1, center=w)

@pytest.mark.parametrize("window", [timedelta(days=3), Timedelta(days=3)])
def test_constructor_with_timedelta_window(window):
    # GH 15440
    n = 10
    df = DataFrame(
        {"value": np.arange(n)}, index=date_range("2015-12-24", periods=n, freq="D")
    )
    expected_data = np.append([0.0, 1.0], np.arange(3.0, 27.0, 3))

    result = df.rolling(window=window).sum()
    expected = DataFrame(
        {"value": expected_data},
        index=date_range("2015-12-24", periods=n, freq="D"),
    )
    tm.assert_frame_equal(result, expected)
    expected = df.rolling("3D").sum()
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize("window", [timedelta(days=3), Timedelta(days=3), "3D"])
def test_constructor_timedelta_window_and_minperiods(window, raw):
    # GH 15305
    n = 10
    df = DataFrame(
        {"value": np.arange(n)}, index=date_range("2017-08-08", periods=n, freq="D")
    )
    expected = DataFrame(
        {"value": np.append([np.NaN, 1.0], np.arange(3.0, 27.0, 3))},
        index=date_range("2017-08-08", periods=n, freq="D"),
    )
    result_roll_sum = df.rolling(window=window, min_periods=2).sum()
    result_roll_generic = df.rolling(window=window, min_periods=2).apply(sum, raw=raw)
    tm.assert_frame_equal(result_roll_sum, expected)
    tm.assert_frame_equal(result_roll_generic, expected)

@pytest.mark.parametrize("method", ["std", "mean", "sum", "max", "min", "var"])
def test_numpy_compat(method):
    # see gh-12811
    r = Rolling(Series([2, 4, 6]), window=2)

    msg = "numpy operations are not valid with window objects"

    with pytest.raises(UnsupportedFunctionCall, match=msg):
        getattr(r, method)(1, 2, 3)
    with pytest.raises(UnsupportedFunctionCall, match=msg):
        getattr(r, method)(dtype=np.float64)

def test_closed_fixed(closed, arithmetic_win_operators):
    # GH 34315
    func_name = arithmetic_win_operators
    df_fixed = DataFrame({"A": [0, 1, 2, 3, 4]})
    df_time = DataFrame({"A": [0, 1, 2, 3, 4]}, index=date_range("2020", periods=5))

    result = getattr(df_fixed.rolling(2, closed=closed, min_periods=1), func_name)()
    expected = getattr(df_time.rolling("2D", closed=closed), func_name)().reset_index(
        drop=True
    )

    tm.assert_frame_equal(result, expected)

def test_closed_fixed_binary_col():
    # GH 34315
    data = [0, 1, 1, 0, 0, 1, 0, 1]
    df = DataFrame(
        {"binary_col": data},
        index=date_range(start="2020-01-01", freq="min", periods=len(data)),
    )

    rolling = df.rolling(window=len(df), closed="left", min_periods=1)
    result = rolling.mean()
    expected = DataFrame(
        [np.nan, 0, 0.5, 2 / 3, 0.5, 0.4, 0.5, 0.428571],
        columns=["binary_col"],
        index=date_range(start="2020-01-01", freq="min", periods=len(data)),
    )
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize("closed", ["neither", "left"])
def test_closed_empty(closed, arithmetic_win_operators):
    # GH 26005
    func_name = arithmetic_win_operators
    ser = Series(data=np.arange(5), index=date_range("2000", periods=5, freq="2D"))
    roll = ser.rolling("1D", closed=closed)

    result = getattr(roll, func_name)()
    expected = Series([np.nan] * 5, index=ser.index)
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize("func", ["min", "max"])
def test_closed_one_entry(func):
    # GH24718
    ser = Series(data=[2], index=date_range("2000", periods=1))
    result = getattr(ser.rolling("10D", closed="left"), func)()
    tm.assert_series_equal(result, Series([np.nan], index=ser.index))

@pytest.mark.parametrize("func", ["min", "max"])
def test_closed_one_entry_groupby(func):
    # GH24718
    ser = DataFrame(
        data={"A": [1, 1, 2], "B": [3, 2, 1]}, index=date_range("2000", periods=3)
    )
    result = getattr(
        ser.groupby("A", sort=False)["B"].rolling("10D", closed="left"), func
    )()
    exp_idx = MultiIndex.from_arrays(arrays=[[1, 1, 2], ser.index], names=("A", None))
    expected = Series(data=[np.nan, 3, np.nan], index=exp_idx, name="B")
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize("input_dtype", ["int", "float"])
@pytest.mark.parametrize(
    "func,closed,expected",
    [
        ("min", "right", [0.0, 0, 0, 1, 2, 3, 4, 5, 6, 7]),
        ("min", "both", [0.0, 0, 0, 0, 1, 2, 3, 4, 5, 6]),
        ("min", "neither", [np.nan, 0, 0, 1, 2, 3, 4, 5, 6, 7]),
        ("min", "left", [np.nan, 0, 0, 0, 1, 2, 3, 4, 5, 6]),
        ("max", "right", [0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),
        ("max", "both", [0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),
        ("max", "neither", [np.nan, 0, 1, 2, 3, 4, 5, 6, 7, 8]),
        ("max", "left", [np.nan, 0, 1, 2, 3, 4, 5, 6, 7, 8]),
    ],
)
def test_closed_min_max_datetime(input_dtype, func, closed, expected):
    # see gh-21704
    ser = Series(
        data=np.arange(10).astype(input_dtype), index=date_range("2000", periods=10)
    )

    result = getattr(ser.rolling("3D", closed=closed), func)()
    expected = Series(expected, index=ser.index)
    tm.assert_series_equal(result, expected)

def test_closed_uneven():
    # see gh-21704
    ser = Series(data=np.arange(10), index=date_range("2000", periods=10))

    # uneven
    ser = ser.drop(index=ser.index[[1, 5]])
    result = ser.rolling("3D", closed="left").min()
    expected = Series([np.nan, 0, 0, 2, 3, 4, 6, 6], index=ser.index)
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize(
    "func,closed,expected",
    [
        ("min", "right", [np.nan, 0, 0, 1, 2, 3, 4, 5, np.nan, np.nan]),
        ("min", "both", [np.nan, 0, 0, 0, 1, 2, 3, 4, 5, np.nan]),
        ("min", "neither", [np.nan, np.nan, 0, 1, 2, 3, 4, 5, np.nan, np.nan]),
        ("min", "left", [np.nan, np.nan, 0, 0, 1, 2, 3, 4, 5, np.nan]),
        ("max", "right", [np.nan, 1, 2, 3, 4, 5, 6, 6, np.nan, np.nan]),
        ("max", "both", [np.nan, 1, 2, 3, 4, 5, 6, 6, 6, np.nan]),
        ("max", "neither", [np.nan, np.nan, 1, 2, 3, 4, 5, 6, np.nan, np.nan]),
        ("max", "left", [np.nan, np.nan, 1, 2, 3, 4, 5, 6, 6, np.nan]),
    ],
)
def test_closed_min_max_minp(func, closed, expected):
    # see gh-21704
    ser = Series(data=np.arange(10), index=date_range("2000", periods=10))
    ser[ser.index[-3:]] = np.nan
    result = getattr(ser.rolling("3D", min_periods=2, closed=closed), func)()
    expected = Series(expected, index=ser.index)
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize(
    "closed,expected",
    [
        ("right", [0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8]),
        ("both", [0, 0.5, 1, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5]),
        ("neither", [np.nan, 0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5]),
        ("left", [np.nan, 0, 0.5, 1, 2, 3, 4, 5, 6, 7]),
    ],
)
def test_closed_median_quantile(closed, expected):
    # GH 26005
    ser = Series(data=np.arange(10), index=date_range("2000", periods=10))
    roll = ser.rolling("3D", closed=closed)
    expected = Series(expected, index=ser.index)

    result = roll.median()
    tm.assert_series_equal(result, expected)

    result = roll.quantile(0.5)
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize("roller", ["1s", 1])
def tests_empty_df_rolling(roller):
    # GH 15819 Verifies that datetime and integer rolling windows can be
    # applied to empty DataFrames
    expected = DataFrame()
    result = DataFrame().rolling(roller).sum()
    tm.assert_frame_equal(result, expected)

    # Verifies that datetime and integer rolling windows can be applied to
    # empty DataFrames with datetime index
    expected = DataFrame(index=DatetimeIndex([]))
    result = DataFrame(index=DatetimeIndex([])).rolling(roller).sum()
    tm.assert_frame_equal(result, expected)

def test_empty_window_median_quantile():
    # GH 26005
    expected = Series([np.nan, np.nan, np.nan])
    roll = Series(np.arange(3)).rolling(0)

    result = roll.median()
    tm.assert_series_equal(result, expected)

    result = roll.quantile(0.1)
    tm.assert_series_equal(result, expected)

def test_missing_minp_zero():
    # https://github.com/pandas-dev/pandas/pull/18921
    # minp=0
    x = Series([np.nan])
    result = x.rolling(1, min_periods=0).sum()
    expected = Series([0.0])
    tm.assert_series_equal(result, expected)

    # minp=1
    result = x.rolling(1, min_periods=1).sum()
    expected = Series([np.nan])
    tm.assert_series_equal(result, expected)

def test_missing_minp_zero_variable():
    # https://github.com/pandas-dev/pandas/pull/18921
    x = Series(
        [np.nan] * 4,
        index=DatetimeIndex(["2017-01-01", "2017-01-04", "2017-01-06", "2017-01-07"]),
    )
    result = x.rolling(Timedelta("2d"), min_periods=0).sum()
    expected = Series(0.0, index=x.index)
    tm.assert_series_equal(result, expected)

def test_multi_index_names():

    # GH 16789, 16825
    cols = MultiIndex.from_product([["A", "B"], ["C", "D", "E"]], names=["1", "2"])
    df = DataFrame(np.ones((10, 6)), columns=cols)
    result = df.rolling(3).cov()

    tm.assert_index_equal(result.columns, df.columns)
    assert result.index.names == [None, "1", "2"]

def test_rolling_axis_sum(axis_frame):
    # see gh-23372.
    df = DataFrame(np.ones((10, 20)))
    axis = df._get_axis_number(axis_frame)

    if axis == 0:
        expected = DataFrame({i: [np.nan] * 2 + [3.0] * 8 for i in range(20)})
    else:
        # axis == 1
        expected = DataFrame([[np.nan] * 2 + [3.0] * 18] * 10)

    result = df.rolling(3, axis=axis_frame).sum()
    tm.assert_frame_equal(result, expected)

def test_rolling_axis_count(axis_frame):
    # see gh-26055
    df = DataFrame({"x": range(3), "y": range(3)})

    axis = df._get_axis_number(axis_frame)

    if axis in [0, "index"]:
        expected = DataFrame({"x": [1.0, 2.0, 2.0], "y": [1.0, 2.0, 2.0]})
    else:
        expected = DataFrame({"x": [1.0, 1.0, 1.0], "y": [2.0, 2.0, 2.0]})

    result = df.rolling(2, axis=axis_frame, min_periods=0).count()
    tm.assert_frame_equal(result, expected)

def test_readonly_array():
    # GH-27766
    arr = np.array([1, 3, np.nan, 3, 5])
    arr.setflags(write=False)
    result = Series(arr).rolling(2).mean()
    expected = Series([np.nan, 2, np.nan, np.nan, 4])
    tm.assert_series_equal(result, expected)

def test_rolling_datetime(axis_frame, tz_naive_fixture):
    # GH-28192
    tz = tz_naive_fixture
    df = DataFrame(
        {i: [1] * 2 for i in date_range("2019-8-01", "2019-08-03", freq="D", tz=tz)}
    )
    if axis_frame in [0, "index"]:
        result = df.T.rolling("2D", axis=axis_frame).sum().T
    else:
        result = df.rolling("2D", axis=axis_frame).sum()
    expected = DataFrame(
        {
            **{
                i: [1.0] * 2
                for i in date_range("2019-8-01", periods=1, freq="D", tz=tz)
            },
            **{
                i: [2.0] * 2
                for i in date_range("2019-8-02", "2019-8-03", freq="D", tz=tz)
            },
        }
    )
    tm.assert_frame_equal(result, expected)

def test_rolling_window_as_string(using_array_manager):
    # see gh-22590
    date_today = datetime.now()
    days = date_range(date_today, date_today + timedelta(365), freq="D")

    npr = np.random.RandomState(seed=421)

    data = npr.randint(1, high=100, size=len(days))
    df = DataFrame({"DateCol": days, "metric": data})

    df.set_index("DateCol", inplace=True)
    result = df.rolling(window="21D", min_periods=2, closed="left")["metric"].agg("max")

    expData = (
        [np.nan] * 2
        + [88.0] * 16
        + [97.0] * 9
        + [98.0]
        + [99.0] * 21
        + [95.0] * 16
        + [93.0] * 5
        + [89.0] * 5
        + [96.0] * 21
        + [94.0] * 14
        + [90.0] * 13
        + [88.0] * 2
        + [90.0] * 9
        + [96.0] * 21
        + [95.0] * 6
        + [91.0]
        + [87.0] * 6
        + [92.0] * 21
        + [83.0] * 2
        + [86.0] * 10
        + [87.0] * 5
        + [98.0] * 21
        + [97.0] * 14
        + [93.0] * 7
        + [87.0] * 4
        + [86.0] * 4
        + [95.0] * 21
        + [85.0] * 14
        + [83.0] * 2
        + [76.0] * 5
        + [81.0] * 2
        + [98.0] * 21
        + [95.0] * 14
        + [91.0] * 7
        + [86.0]
        + [93.0] * 3
        + [95.0] * 20
    )

    index = days.rename("DateCol")
    if not using_array_manager:
        # INFO(ArrayManager) preserves the frequence of the index
        index = index._with_freq(None)
    expected = Series(expData, index=index, name="metric")
    tm.assert_series_equal(result, expected)

def test_min_periods1():
    # GH#6795
    df = DataFrame([0, 1, 2, 1, 0], columns=["a"])
    result = df["a"].rolling(3, center=True, min_periods=1).max()
    expected = Series([1.0, 2.0, 2.0, 2.0, 1.0], name="a")
    tm.assert_series_equal(result, expected)

def test_rolling_count_with_min_periods(frame_or_series):
    # GH 26996
    result = frame_or_series(range(5)).rolling(3, min_periods=3).count()
    expected = frame_or_series([np.nan, np.nan, 3.0, 3.0, 3.0])
    tm.assert_equal(result, expected)

def test_rolling_count_default_min_periods_with_null_values(frame_or_series):
    # GH 26996
    values = [1, 2, 3, np.nan, 4, 5, 6]
    expected_counts = [1.0, 2.0, 3.0, 2.0, 2.0, 2.0, 3.0]

    # GH 31302
    with tm.assert_produces_warning(FutureWarning, check_stacklevel=False):
        result = frame_or_series(values).rolling(3).count()
    expected = frame_or_series(expected_counts)
    tm.assert_equal(result, expected)

@pytest.mark.parametrize(
    "df,expected,window,min_periods",
    [
        (
            DataFrame({"A": [1, 2, 3], "B": [4, 5, 6]}),
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [1, 2], "B": [4, 5]}, [0, 1]),
                ({"A": [1, 2, 3], "B": [4, 5, 6]}, [0, 1, 2]),
            ],
            3,
            None,
        ),
        (
            DataFrame({"A": [1, 2, 3], "B": [4, 5, 6]}),
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [1, 2], "B": [4, 5]}, [0, 1]),
                ({"A": [2, 3], "B": [5, 6]}, [1, 2]),
            ],
            2,
            1,
        ),
        (
            DataFrame({"A": [1, 2, 3], "B": [4, 5, 6]}),
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [1, 2], "B": [4, 5]}, [0, 1]),
                ({"A": [2, 3], "B": [5, 6]}, [1, 2]),
            ],
            2,
            2,
        ),
        (
            DataFrame({"A": [1, 2, 3], "B": [4, 5, 6]}),
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [2], "B": [5]}, [1]),
                ({"A": [3], "B": [6]}, [2]),
            ],
            1,
            1,
        ),
        (
            DataFrame({"A": [1, 2, 3], "B": [4, 5, 6]}),
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [2], "B": [5]}, [1]),
                ({"A": [3], "B": [6]}, [2]),
            ],
            1,
            0,
        ),
        (DataFrame({"A": [1], "B": [4]}), [], 2, None),
        (DataFrame({"A": [1], "B": [4]}), [], 2, 1),
        (DataFrame(), [({}, [])], 2, None),
        (
            DataFrame({"A": [1, np.nan, 3], "B": [np.nan, 5, 6]}),
            [
                ({"A": [1.0], "B": [np.nan]}, [0]),
                ({"A": [1, np.nan], "B": [np.nan, 5]}, [0, 1]),
                ({"A": [1, np.nan, 3], "B": [np.nan, 5, 6]}, [0, 1, 2]),
            ],
            3,
            2,
        ),
    ],
)
def test_iter_rolling_dataframe(df, expected, window, min_periods):
    # GH 11704
    expected = [DataFrame(values, index=index) for (values, index) in expected]

    for (expected, actual) in zip(
        expected, df.rolling(window, min_periods=min_periods)
    ):
        tm.assert_frame_equal(actual, expected)

@pytest.mark.parametrize(
    "expected,window",
    [
        (
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [1, 2], "B": [4, 5]}, [0, 1]),
                ({"A": [2, 3], "B": [5, 6]}, [1, 2]),
            ],
            "2D",
        ),
        (
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [1, 2], "B": [4, 5]}, [0, 1]),
                ({"A": [1, 2, 3], "B": [4, 5, 6]}, [0, 1, 2]),
            ],
            "3D",
        ),
        (
            [
                ({"A": [1], "B": [4]}, [0]),
                ({"A": [2], "B": [5]}, [1]),
                ({"A": [3], "B": [6]}, [2]),
            ],
            "1D",
        ),
    ],
)
def test_iter_rolling_on_dataframe(expected, window):
    # GH 11704
    df = DataFrame(
        {
            "A": [1, 2, 3, 4, 5],
            "B": [4, 5, 6, 7, 8],
            "C": date_range(start="2016-01-01", periods=5, freq="D"),
        }
    )

    expected = [DataFrame(values, index=index) for (values, index) in expected]
    for (expected, actual) in zip(expected, df.rolling(window, on="C")):
        tm.assert_frame_equal(actual, expected)

@pytest.mark.parametrize(
    "ser,expected,window, min_periods",
    [
        (
            Series([1, 2, 3]),
            [([1], [0]), ([1, 2], [0, 1]), ([1, 2, 3], [0, 1, 2])],
            3,
            None,
        ),
        (
            Series([1, 2, 3]),
            [([1], [0]), ([1, 2], [0, 1]), ([1, 2, 3], [0, 1, 2])],
            3,
            1,
        ),
        (Series([1, 2, 3]), [([1], [0]), ([1, 2], [0, 1]), ([2, 3], [1, 2])], 2, 1),
        (Series([1, 2, 3]), [([1], [0]), ([1, 2], [0, 1]), ([2, 3], [1, 2])], 2, 2),
        (Series([1, 2, 3]), [([1], [0]), ([2], [1]), ([3], [2])], 1, 0),
        (Series([1, 2, 3]), [([1], [0]), ([2], [1]), ([3], [2])], 1, 1),
        (Series([1, 2]), [([1], [0]), ([1, 2], [0, 1])], 2, 0),
        (Series([], dtype="int64"), [], 2, 1),
    ],
)
def test_iter_rolling_series(ser, expected, window, min_periods):
    # GH 11704
    expected = [Series(values, index=index) for (values, index) in expected]

    for (expected, actual) in zip(
        expected, ser.rolling(window, min_periods=min_periods)
    ):
        tm.assert_series_equal(actual, expected)

@pytest.mark.parametrize(
    "expected,expected_index,window",
    [
        (
            [[0], [1], [2], [3], [4]],
            [
                date_range("2020-01-01", periods=1, freq="D"),
                date_range("2020-01-02", periods=1, freq="D"),
                date_range("2020-01-03", periods=1, freq="D"),
                date_range("2020-01-04", periods=1, freq="D"),
                date_range("2020-01-05", periods=1, freq="D"),
            ],
            "1D",
        ),
        (
            [[0], [0, 1], [1, 2], [2, 3], [3, 4]],
            [
                date_range("2020-01-01", periods=1, freq="D"),
                date_range("2020-01-01", periods=2, freq="D"),
                date_range("2020-01-02", periods=2, freq="D"),
                date_range("2020-01-03", periods=2, freq="D"),
                date_range("2020-01-04", periods=2, freq="D"),
            ],
            "2D",
        ),
        (
            [[0], [0, 1], [0, 1, 2], [1, 2, 3], [2, 3, 4]],
            [
                date_range("2020-01-01", periods=1, freq="D"),
                date_range("2020-01-01", periods=2, freq="D"),
                date_range("2020-01-01", periods=3, freq="D"),
                date_range("2020-01-02", periods=3, freq="D"),
                date_range("2020-01-03", periods=3, freq="D"),
            ],
            "3D",
        ),
    ],
)
def test_iter_rolling_datetime(expected, expected_index, window):
    # GH 11704
    ser = Series(range(5), index=date_range(start="2020-01-01", periods=5, freq="D"))

    expected = [
        Series(values, index=idx) for (values, idx) in zip(expected, expected_index)
    ]

    for (expected, actual) in zip(expected, ser.rolling(window)):
        tm.assert_series_equal(actual, expected)

@pytest.mark.parametrize(
    "grouping,_index",
    [
        (
            {"level": 0},
            MultiIndex.from_tuples(
                [(0, 0), (0, 0), (1, 1), (1, 1), (1, 1)], names=[None, None]
            ),
        ),
        (
            {"by": "X"},
            MultiIndex.from_tuples(
                [(0, 0), (1, 0), (2, 1), (3, 1), (4, 1)], names=["X", None]
            ),
        ),
    ],
)
def test_rolling_positional_argument(grouping, _index, raw):
    # GH 34605

    def scaled_sum(*args):
        if len(args) < 2:
            raise ValueError("The function needs two arguments")
        array, scale = args
        return array.sum() / scale

    df = DataFrame(data={"X": range(5)}, index=[0, 0, 1, 1, 1])

    expected = DataFrame(data={"X": [0.0, 0.5, 1.0, 1.5, 2.0]}, index=_index)
    result = df.groupby(**grouping).rolling(1).apply(scaled_sum, raw=raw, args=(2,))
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize("add", [0.0, 2.0])
def test_rolling_numerical_accuracy_kahan_mean(add):
    # GH: 36031 implementing kahan summation
    df = DataFrame(
        {"A": [3002399751580331.0 + add, -0.0, -0.0]},
        index=[
            Timestamp("19700101 09:00:00"),
            Timestamp("19700101 09:00:03"),
            Timestamp("19700101 09:00:06"),
        ],
    )
    result = (
        df.resample("1s").ffill().rolling("3s", closed="left", min_periods=3).mean()
    )
    dates = date_range("19700101 09:00:00", periods=7, freq="S")
    expected = DataFrame(
        {
            "A": [
                np.nan,
                np.nan,
                np.nan,
                3002399751580330.5,
                2001599834386887.25,
                1000799917193443.625,
                0.0,
            ]
        },
        index=dates,
    )
    tm.assert_frame_equal(result, expected)

def test_rolling_numerical_accuracy_kahan_sum():
    # GH: 13254
    df = DataFrame([2.186, -1.647, 0.0, 0.0, 0.0, 0.0], columns=["x"])
    result = df["x"].rolling(3).sum()
    expected = Series([np.nan, np.nan, 0.539, -1.647, 0.0, 0.0], name="x")
    tm.assert_series_equal(result, expected)

def test_rolling_numerical_accuracy_jump():
    # GH: 32761
    index = date_range(start="2020-01-01", end="2020-01-02", freq="60s").append(
        DatetimeIndex(["2020-01-03"])
    )
    data = np.random.rand(len(index))

    df = DataFrame({"data": data}, index=index)
    result = df.rolling("60s").mean()
    tm.assert_frame_equal(result, df[["data"]])

def test_rolling_numerical_accuracy_small_values():
    # GH: 10319
    s = Series(
        data=[0.00012456, 0.0003, -0.0, -0.0],
        index=date_range("1999-02-03", "1999-02-06"),
    )
    result = s.rolling(1).mean()
    tm.assert_series_equal(result, s)

def test_rolling_numerical_too_large_numbers():
    # GH: 11645
    dates = date_range("2015-01-01", periods=10, freq="D")
    ds = Series(data=range(10), index=dates, dtype=np.float64)
    ds[2] = -9e33
    result = ds.rolling(5).mean()
    expected = Series(
        [np.nan, np.nan, np.nan, np.nan, -1.8e33, -1.8e33, -1.8e33, 5.0, 6.0, 7.0],
        index=dates,
    )
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize(
    ("func", "value"),
    [("sum", 2.0), ("max", 1.0), ("min", 1.0), ("mean", 1.0), ("median", 1.0)],
)
def test_rolling_mixed_dtypes_axis_1(func, value):
    # GH: 20649
    df = DataFrame(1, index=[1, 2], columns=["a", "b", "c"])
    df["c"] = 1.0
    result = getattr(df.rolling(window=2, min_periods=1, axis=1), func)()
    expected = DataFrame(
        {"a": [1.0, 1.0], "b": [value, value], "c": [value, value]}, index=[1, 2]
    )
    tm.assert_frame_equal(result, expected)

def test_rolling_axis_one_with_nan():
    # GH: 35596
    df = DataFrame(
        [
            [0, 1, 2, 4, np.nan, np.nan, np.nan],
            [0, 1, 2, np.nan, np.nan, np.nan, np.nan],
            [0, 2, 2, np.nan, 2, np.nan, 1],
        ]
    )
    result = df.rolling(window=7, min_periods=1, axis="columns").sum()
    expected = DataFrame(
        [
            [0.0, 1.0, 3.0, 7.0, 7.0, 7.0, 7.0],
            [0.0, 1.0, 3.0, 3.0, 3.0, 3.0, 3.0],
            [0.0, 2.0, 4.0, 4.0, 6.0, 6.0, 7.0],
        ]
    )
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize(
    "value",
    ["test", to_datetime("2019-12-31"), to_timedelta("1 days 06:05:01.00003")],
)
def test_rolling_axis_1_non_numeric_dtypes(value):
    # GH: 20649
    df = DataFrame({"a": [1, 2]})
    df["b"] = value
    result = df.rolling(window=2, min_periods=1, axis=1).sum()
    expected = DataFrame({"a": [1.0, 2.0]})
    tm.assert_frame_equal(result, expected)

def test_rolling_on_df_transposed():
    # GH: 32724
    df = DataFrame({"A": [1, None], "B": [4, 5], "C": [7, 8]})
    expected = DataFrame({"A": [1.0, np.nan], "B": [5.0, 5.0], "C": [11.0, 13.0]})
    result = df.rolling(min_periods=1, window=2, axis=1).sum()
    tm.assert_frame_equal(result, expected)

    result = df.T.rolling(min_periods=1, window=2).sum().T
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize(
    ("index", "window"),
    [
        (
            period_range(start="2020-01-01 08:00", end="2020-01-01 08:08", freq="T"),
            "2T",
        ),
        (
            period_range(start="2020-01-01 08:00", end="2020-01-01 12:00", freq="30T"),
            "1h",
        ),
    ],
)
@pytest.mark.parametrize(
    ("func", "values"),
    [
        ("min", [np.nan, 0, 0, 1, 2, 3, 4, 5, 6]),
        ("max", [np.nan, 0, 1, 2, 3, 4, 5, 6, 7]),
        ("sum", [np.nan, 0, 1, 3, 5, 7, 9, 11, 13]),
    ],
)
def test_rolling_period_index(index, window, func, values):
    # GH: 34225
    ds = Series([0, 1, 2, 3, 4, 5, 6, 7, 8], index=index)
    result = getattr(ds.rolling(window, closed="left"), func)()
    expected = Series(values, index=index)
    tm.assert_series_equal(result, expected)

def test_rolling_sem(frame_or_series):
    # GH: 26476
    obj = frame_or_series([0, 1, 2])
    result = obj.rolling(2, min_periods=1).sem()
    if isinstance(result, DataFrame):
        result = Series(result[0].values)
    expected = Series([np.nan] + [0.707107] * 2)
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize(
    ("func", "third_value", "values"),
    [
        ("var", 1, [5e33, 0, 0.5, 0.5, 2, 0]),
        ("std", 1, [7.071068e16, 0, 0.7071068, 0.7071068, 1.414214, 0]),
        ("var", 2, [5e33, 0.5, 0, 0.5, 2, 0]),
        ("std", 2, [7.071068e16, 0.7071068, 0, 0.7071068, 1.414214, 0]),
    ],
)
def test_rolling_var_numerical_issues(func, third_value, values):
    # GH: 37051
    ds = Series([99999999999999999, 1, third_value, 2, 3, 1, 1])
    result = getattr(ds.rolling(2), func)()
    expected = Series([np.nan] + values)
    tm.assert_series_equal(result, expected)

def test_timeoffset_as_window_parameter_for_corr():
    # GH: 28266
    exp = DataFrame(
        {
            "B": [
                np.nan,
                np.nan,
                0.9999999999999998,
                -1.0,
                1.0,
                -0.3273268353539892,
                0.9999999999999998,
                1.0,
                0.9999999999999998,
                1.0,
            ],
            "A": [
                np.nan,
                np.nan,
                -1.0,
                1.0000000000000002,
                -0.3273268353539892,
                0.9999999999999966,
                1.0,
                1.0000000000000002,
                1.0,
                1.0000000000000002,
            ],
        },
        index=MultiIndex.from_tuples(
            [
                (Timestamp("20130101 09:00:00"), "B"),
                (Timestamp("20130101 09:00:00"), "A"),
                (Timestamp("20130102 09:00:02"), "B"),
                (Timestamp("20130102 09:00:02"), "A"),
                (Timestamp("20130103 09:00:03"), "B"),
                (Timestamp("20130103 09:00:03"), "A"),
                (Timestamp("20130105 09:00:05"), "B"),
                (Timestamp("20130105 09:00:05"), "A"),
                (Timestamp("20130106 09:00:06"), "B"),
                (Timestamp("20130106 09:00:06"), "A"),
            ]
        ),
    )

    df = DataFrame(
        {"B": [0, 1, 2, 4, 3], "A": [7, 4, 6, 9, 3]},
        index=[
            Timestamp("20130101 09:00:00"),
            Timestamp("20130102 09:00:02"),
            Timestamp("20130103 09:00:03"),
            Timestamp("20130105 09:00:05"),
            Timestamp("20130106 09:00:06"),
        ],
    )

    res = df.rolling(window="3d").corr()

    tm.assert_frame_equal(exp, res)

@pytest.mark.parametrize("method", ["var", "sum", "mean", "skew", "kurt", "min", "max"])
def test_rolling_decreasing_indices(method):
    """
    Make sure that decreasing indices give the same results as increasing indices.

    GH 36933
    """
    df = DataFrame({"values": np.arange(-15, 10) ** 2})
    df_reverse = DataFrame({"values": df["values"][::-1]}, index=df.index[::-1])

    increasing = getattr(df.rolling(window=5), method)()
    decreasing = getattr(df_reverse.rolling(window=5), method)()

    assert np.abs(decreasing.values[::-1][:-4] - increasing.values[4:]).max() < 1e-12

@pytest.mark.parametrize(
    "method,expected",
    [
        (
            "var",
            [
                float("nan"),
                43.0,
                float("nan"),
                136.333333,
                43.5,
                94.966667,
                182.0,
                318.0,
            ],
        ),
        ("mean", [float("nan"), 7.5, float("nan"), 21.5, 6.0, 9.166667, 13.0, 17.5]),
        ("sum", [float("nan"), 30.0, float("nan"), 86.0, 30.0, 55.0, 91.0, 140.0]),
        (
            "skew",
            [
                float("nan"),
                0.709296,
                float("nan"),
                0.407073,
                0.984656,
                0.919184,
                0.874674,
                0.842418,
            ],
        ),
        (
            "kurt",
            [
                float("nan"),
                -0.5916711736073559,
                float("nan"),
                -1.0028993131317954,
                -0.06103844629409494,
                -0.254143227116194,
                -0.37362637362637585,
                -0.45439658241367054,
            ],
        ),
    ],
)
def test_rolling_non_monotonic(method, expected):
    """
    Make sure the (rare) branch of non-monotonic indices is covered by a test.

    output from 1.1.3 is assumed to be the expected output. Output of sum/mean has
    manually been verified.

    GH 36933.
    """
    # Based on an example found in computation.rst
    use_expanding = [True, False, True, False, True, True, True, True]
    df = DataFrame({"values": np.arange(len(use_expanding)) ** 2})

    class CustomIndexer(BaseIndexer):
        def get_window_bounds(self, num_values, min_periods, center, closed):
            start = np.empty(num_values, dtype=np.int64)
            end = np.empty(num_values, dtype=np.int64)
            for i in range(num_values):
                if self.use_expanding[i]:
                    start[i] = 0
                    end[i] = i + 1
                else:
                    start[i] = i
                    end[i] = i + self.window_size
            return start, end

    indexer = CustomIndexer(window_size=4, use_expanding=use_expanding)

    result = getattr(df.rolling(indexer), method)()
    expected = DataFrame({"values": expected})
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize(
    ("index", "window"),
    [([0, 1, 2, 3, 4], 2), (date_range("2001-01-01", freq="D", periods=5), "2D")],
)
def test_rolling_corr_timedelta_index(index, window):
    # GH: 31286
    x = Series([1, 2, 3, 4, 5], index=index)
    y = x.copy()
    x[0:2] = 0.0
    result = x.rolling(window).corr(y)
    expected = Series([np.nan, np.nan, 1, 1, 1], index=index)
    tm.assert_almost_equal(result, expected)

def test_groupby_rolling_nan_included():
    # GH 35542
    data = {"group": ["g1", np.nan, "g1", "g2", np.nan], "B": [0, 1, 2, 3, 4]}
    df = DataFrame(data)
    result = df.groupby("group", dropna=False).rolling(1, min_periods=1).mean()
    expected = DataFrame(
        {"B": [0.0, 2.0, 3.0, 1.0, 4.0]},
        # GH-38057 from_tuples puts the NaNs in the codes, result expects them
        # to be in the levels, at the moment
        # index=MultiIndex.from_tuples(
        #     [("g1", 0), ("g1", 2), ("g2", 3), (np.nan, 1), (np.nan, 4)],
        #     names=["group", None],
        # ),
        index=MultiIndex(
            [["g1", "g2", np.nan], [0, 1, 2, 3, 4]],
            [[0, 0, 1, 2, 2], [0, 2, 3, 1, 4]],
            names=["group", None],
        ),
    )
    tm.assert_frame_equal(result, expected)

@pytest.mark.parametrize("method", ["skew", "kurt"])
def test_rolling_skew_kurt_numerical_stability(method):
    # GH#6929
    ser = Series(np.random.rand(10))
    ser_copy = ser.copy()
    expected = getattr(ser.rolling(3), method)()
    tm.assert_series_equal(ser, ser_copy)
    ser = ser + 50000
    result = getattr(ser.rolling(3), method)()
    tm.assert_series_equal(result, expected)

@pytest.mark.parametrize(
    ("method", "values"),
    [
        ("skew", [2.0, 0.854563, 0.0, 1.999984]),
        ("kurt", [4.0, -1.289256, -1.2, 3.999946]),
    ],
)
def test_rolling_skew_kurt_large_value_range(method, values):
    # GH: 37557
    s = Series([3000000, 1, 1, 2, 3, 4, 999])
    result = getattr(s.rolling(4), method)()
    expected = Series([np.nan] * 3 + values)
    tm.assert_series_equal(result, expected)

def test_invalid_method():
    with pytest.raises(ValueError, match="method must be 'table' or 'single"):
        Series(range(1)).rolling(1, method="foo")

@pytest.mark.parametrize("window", [1, "1d"])
def test_rolling_descending_date_order_with_offset(window, frame_or_series):
    # GH#40002
    idx = date_range(start="2020-01-01", end="2020-01-03", freq="1d")
    obj = frame_or_series(range(1, 4), index=idx)
    result = obj.rolling("1d", closed="left").sum()
    expected = frame_or_series([np.nan, 1, 2], index=idx)
    tm.assert_equal(result, expected)

    result = obj.iloc[::-1].rolling("1d", closed="left").sum()
    idx = date_range(start="2020-01-03", end="2020-01-01", freq="-1d")
    expected = frame_or_series([np.nan, 3, 2], index=idx)
    tm.assert_equal(result, expected)

def test_rolling_var_floating_artifact_precision():
    # GH 37051
    s = Series([7, 5, 5, 5])
    result = s.rolling(3).var()
    expected = Series([np.nan, np.nan, 4 / 3, 0])
    tm.assert_series_equal(result, expected, atol=1.0e-15, rtol=1.0e-15)

def test_rolling_std_small_values():
    # GH 37051
    s = Series(
        [
            0.00000054,
            0.00000053,
            0.00000054,
        ]
    )
    result = s.rolling(2).std()
    expected = Series([np.nan, 7.071068e-9, 7.071068e-9])
    tm.assert_series_equal(result, expected, atol=1.0e-15, rtol=1.0e-15)

max_stars_repo_name: mdengler/pandas

max_issues_repo_name: mdengler/pandas

max_forks_repo_name: mdengler/pandas


text:
################## HEADER DO NOT EDIT ##################
import os
import __main__

### Get the name of program importing this package ###
if hasattr(__main__,"__file__"):
    main_file = os.path.basename(__main__.__file__)
else:
    main_file = "ipython"

### This checks if we are just doing documentation ###
if not main_file in ["sphinx-build", "__main__.py"]:
    from dolfin import *
    from dolfin_adjoint import *
########################################################

### Additional import statements ###
import numpy as np
import math
import os

### Declare Unique name
name = "KE_entrainment"

### Set default keyword argument values ###
keyword_defaults = {
    "ke_location": "rotor", # location of measurement hub, rotor, tip (only rotor works for now)
    }

### Define objective function
def objective(solver, inflow_angle = 0.0, first_call=False, **kwargs):
    '''
    The "KE_entrainment" objective function computed the vertical kinetic
    entrainment behind a single turbine

    Keyword arguments:
        ke_location: location of measurement, hub, rotor, tip (only rotor works for now)
    '''

    ### Extract keyword arguments
    ke_location = kwargs.pop("ke_location")

    solver.fprint("Using Kinetic Energy Entrainment Functional")
    turb_id = solver.opt_turb_id[0]
    HH = solver.problem.farm.HH[turb_id]
    R = solver.problem.farm.RD[turb_id]/2.0

    # mark cells in area of interest
    if first_call:
        solver.KE_objective_markers = MeshFunction("size_t", solver.problem.dom.mesh, solver.problem.dom.mesh.topology().dim())
        # solver.KE_objective_markers = MeshFunction("size_t", solver.problem.dom.mesh, solver.problem.dom.mesh.topology().dim()-1)

        solver.KE_objective_markers.set_all(0)
        x0 = min(solver.problem.farm.x)
        x1 = solver.problem.dom.x_range[1]
        y0 = min(solver.problem.farm.y)-3.0*max(solver.problem.farm.RD)/2.0
        y1 = max(solver.problem.farm.y)+3.0*max(solver.problem.farm.RD)/2.0
        # if ke_location =="tip":
        #     z0 = min(solver.problem.farm.z)
        #     z1 = max(solver.problem.farm.z)+max(solver.problem.farm.RD)/2.0
        # elif ke_location == "hub":
        #     z0 = min(solver.problem.farm.z)-solver.problem.dom.mesh.hmin()*1.0
        #     z1 = max(solver.problem.farm.z)+solver.problem.dom.mesh.hmin()*1.0
        # elif ke_location == "rotor":
        z0_in = min(solver.problem.farm.z) + R - solver.problem.dom.mesh.hmin()*1.0
        z1_in = max(solver.problem.farm.z) + R + solver.problem.dom.mesh.hmin()*1.0
        z0_out = min(solver.problem.farm.z) - R - solver.problem.dom.mesh.hmin()*1.0
        z1_out = max(solver.problem.farm.z) - R + solver.problem.dom.mesh.hmin()*1.0

        fluxIn =  CompiledSubDomain("x[0]>x0 && x[0]<x1 && x[1]>y0 && x[1]<y1  && x[2]>z0 && x[2]<z1",x0=x0,x1=x1,y0=y0,y1=y1,z0=z0_in,z1=z1_in)
        fluxOut = CompiledSubDomain("x[0]>x0 && x[0]<x1 && x[1]>y0 && x[1]<y1  && x[2]>z0 && x[2]<z1",x0=x0,x1=x1,y0=y0,y1=y1,z0=z0_out,z1=z1_out)
        # AOI  = CompiledSubDomain("x[0]>x0 && x[0]<x1 && x[1]>y0 && x[1]<y1  && x[2]>z0 && x[2]<z1",x0=x0,x1=x1,y0=y0,y1=y1,z0=z0,z1=z1)
        # AOI.mark(solver.KE_objective_markers,1)
        fluxIn.mark(solver.KE_objective_markers,1)
        fluxOut.mark(solver.KE_objective_markers,2)
        File(solver.params.folder+"mesh/KE_objective_markers.pvd")<<solver.KE_objective_markers

        if solver.save_objective:
            folder_string = solver.params.folder+"data/"
            if not os.path.exists(folder_string): os.makedirs(folder_string)
            f = open(folder_string+"KE_entrainment_data.txt",'w')
            header = str("Time    "+"Entrainment\n")
            f.write(header)
            f.close()

    dx_KE = Measure('dx', subdomain_data=solver.KE_objective_markers)
    # ds_KE = Measure('ds', subdomain_data=solver.KE_objective_markers)

    J = assemble((solver.problem.vertKE*dx_KE(1) - solver.problem.vertKE*dx_KE(2)))

    if solver.save_objective:
        folder_string = solver.params.folder+"data/"
        f = open(folder_string+"KE_entrainment_data.txt",'a')
        out_data = [solver.simTime]
        out_data.extend([float(J)])
        np.savetxt(f,[out_data])
        f.close()

    return J

ext: py

max_stars_repo_name: jefalon/WindSE

max_issues_repo_name: jefalon/WindSE

max_forks_repo_name: jefalon/WindSE

reason: import numpy


text:
import numpy as np
from keras.models import Model
from keras.layers import Input, Add, Activation, Dense, Conv3D, ZeroPadding3D, MaxPooling3D, BatchNormalization, AveragePooling3D, Flatten
from keras.utils import plot_model
from keras.applications.resnet50 import ResNet50

# define a function to add bottleneck
def add_bottle_neck(x, stride_val, filters_1, filters_2, regress_indentity):
    x_identity = x
    x_layer = Conv3D(filters_1, kernel_size=(1, 1, 1), padding='valid', strides=stride_val)(x)
    x_layer = BatchNormalization(axis=-1)(x_layer)
    x_layer = Activation(activation='relu')(x_layer)
    x_layer = Conv3D(filters_1, kernel_size=(3, 3, 1), padding='same', strides=(1, 1, 1))(x_layer)
    x_layer = BatchNormalization(axis=-1)(x_layer)
    x_layer = Activation(activation='relu')(x_layer)
    x_layer = Conv3D(filters_2, kernel_size=(1, 1, 1), padding='valid', strides=(1, 1, 1))(x_layer)
    x_layer = BatchNormalization(axis=-1)(x_layer)
    if (regress_indentity == True):
        shortcut_path = Conv3D(filters_2, kernel_size=(1, 1, 1), padding='valid', strides=stride_val)(x_identity)
        shortcut_path = BatchNormalization(axis=-1)(shortcut_path)
        final_layer = Add()([x_layer, shortcut_path])
    else:
        final_layer = Add()([x_layer, x_identity])

    final_layer = Activation(activation='relu')(final_layer)
    return final_layer

# define model input
input_x = Input(shape=(224, 224, 8, 3))

# define 3d ResNet50
X = ZeroPadding3D((3, 3, 0))(input_x)
X = Conv3D(64, kernel_size=(7, 7, 1), padding='valid', strides=(2, 2, 1))(X)
X = BatchNormalization(axis=-1)(X)
X = Activation(activation='relu')(X)
X = MaxPooling3D(pool_size=(3, 3, 1), strides=(2, 2, 1),padding='valid')(X)
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=64, filters_2=256, regress_indentity=True) # conv2_1
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=64, filters_2=256, regress_indentity=False) # conv2_2
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=64, filters_2=256, regress_indentity=False) # conv2_3
X = add_bottle_neck(X, stride_val=(2, 2, 1), filters_1=128, filters_2=512, regress_indentity=True) # conv3_1
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=128, filters_2=512, regress_indentity=False) # conv3_2
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=128, filters_2=512, regress_indentity=False) # conv3_3
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=128, filters_2=512, regress_indentity=False) # conv3_4
X = add_bottle_neck(X, stride_val=(2, 2, 1), filters_1=256, filters_2=1024, regress_indentity=True) # conv4_1
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=256, filters_2=1024, regress_indentity=False) # conv4_2
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=256, filters_2=1024, regress_indentity=False) # conv4_3
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=256, filters_2=1024, regress_indentity=False) # conv4_4
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=256, filters_2=1024, regress_indentity=False) # conv4_5
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=256, filters_2=1024, regress_indentity=False) # conv4_6
X = add_bottle_neck(X, stride_val=(2, 2, 1), filters_1=512, filters_2=2048, regress_indentity=True) # conv5_1
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=512, filters_2=2048, regress_indentity=False) # conv5_2
X = add_bottle_neck(X, stride_val=(1, 1, 1), filters_1=512, filters_2=2048, regress_indentity=False) # conv5_3
X = AveragePooling3D(pool_size=(7,7,8), strides=None, padding='valid')(X) # Modify for glimpse clouds
X = Flatten()(X)
X = Dense(1000, activation='softmax')(X)

# Create model
resnet50_3d_model = Model(inputs = input_x, outputs = X, name='ResNet50_3D')

# model summary and plot
print(resnet50_3d_model.summary())
plot_model(resnet50_3d_model, to_file='model.png')

# Load pretrained weights for a 2D resnet trained on imagenet

model = ResNet50(weights='imagenet')
print("Model Loaded")

# read weights
weights_2d = model.get_weights()

# Declare array containing weight inxedes which have to be expanded
index_array = [0, 6, 12, 18, 20, 30, 36, 42, 48, 54, 60, 66, 72, 78, 80, 90, 96, 102, 108, 114, 120, 126, 132, 138, 144, 150, 156, 158, 168, 174, 180, 186, 192, 198, 204, 210, 216, 222, 228, 234, 240, 246, 252, 258, 264, 270, 272, 282, 288, 294, 300, 306, 312]

weights_3d = []
for i in range(0, len(weights_2d)):
    if i in index_array:
        #N = weights_2d[i].shape[1]
        N = 1
        weights_3d.append(np.repeat(weights_2d[i][:,:, np.newaxis,:,:], N, axis=2))
        # Normalize weights for H, W, and D to keep activation values same as the 2D network
        weights_3d[i][:, :, :, 1, 1] /= (N+1)
    else:
        weights_3d.append(weights_2d[i])
    print(weights_3d[i].shape)

print("Weight dimensions inflated")

# Load weights into the model
resnet50_3d_model.set_weights(weights_3d)
print("weights loaded into the model")


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Apr  4 14:45:37 2022

@author: xavier.mouy
"""

import argparse
import os
import numpy as np
import tensorflow as tf
import ketos.data_handling.database_interface as dbi
from ketos.neural_networks.resnet import ResNetInterface
from ketos.data_handling.data_feeding import BatchGenerator

## INPUT ARGUMENTS ############################################################

# command line arguments
parser = argparse.ArgumentParser(description="Input arguments")
parser.add_argument("--db_file", type=str, action='store', required=False)
parser.add_argument("--out_dir", type=str, action='store', required=False)
parser.add_argument("--recipe_file", type=str, action='store', required=False)
parser.add_argument("--spec_config_file", type=str, action='store', required=False)
parser.add_argument("--checkpoints_dir", type=str, action='store', required=False)
parser.add_argument("--logs_dir", type=str, action='store', required=False)
parser.add_argument("--batch_size", type=int, action='store', required=False)
parser.add_argument("--n_epochs", type=int, action='store', required=False)
args = parser.parse_args()

# default values
params = dict()
params['db_file'] = r'C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\ketos_databases\spectro-5s_fft-0.128_step-0.064_fmin-0_fmax-800\database.h5'
params['out_dir'] = r'C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\ketos_models'
params['recipe_file'] = r'C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\recipe.json'
params['spec_config_file'] = r'C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\ketos_databases\spectro-5s_fft-0.128_step-0.064_fmin-0_fmax-800\spec_config.json'
params['checkpoints_dir'] = r"C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\checkpoints"
params['logs_dir'] = r"C:\Users\xavier.mouy\Documents\GitHub\Ketos\minke\20220401T140308\logs"
params['batch_size'] = 20 #128
params['shuffle'] = True
params['refresh_on_epoch_end_train'] = True
params['refresh_on_epoch_end_test'] = False
params['n_epochs'] = 30

# parse cli arguments if any
args = parser.parse_args()
if args.db_file:
    params['db_file'] = args.db_file
if args.out_dir:
    params['out_dir'] = args.out_dir
if args.recipe_file:
    params['recipe_file'] = args.recipe_file
if args.spec_config_file:
    params['spec_config_file'] = args.spec_config_file
if args.checkpoints_dir:
    params['checkpoints_dir'] = args.checkpoints_dir
if args.logs_dir:
    params['logs_dir'] = args.logs_dir
if args.batch_size:
    params['batch_size'] = args.batch_size
if args.n_epochs:
    params['n_epochs'] = args.n_epochs

## ############################################################################

# init random seed to keep same results for each run
np.random.seed(1000)
tf.random.set_seed(2000)

# create output folder if it doesn't exist yet
if os.path.isdir(params['out_dir']) is False:
    os.mkdir(params['out_dir'])

# creating data feed
db = dbi.open_file(params['db_file'], 'r')
train_data = dbi.open_table(db, "/train/data")
val_data = dbi.open_table(db, "/test/data")

train_generator = BatchGenerator(batch_size=params['batch_size'],
                                 data_table=train_data,
                                 output_transform_func=ResNetInterface.transform_batch,
                                 shuffle=params['shuffle'],
                                 refresh_on_epoch_end=params['refresh_on_epoch_end_train']
                                 )

val_generator = BatchGenerator(batch_size=params['batch_size'],
                               data_table=val_data,
                               output_transform_func=ResNetInterface.transform_batch,
                               shuffle=True,
                               refresh_on_epoch_end=params['refresh_on_epoch_end_test']
                               )

# Creating neural net structure (from recipe file)
resnet = ResNetInterface.build_from_recipe_file(params['recipe_file'])
resnet.train_generator = train_generator
resnet.val_generator = val_generator
resnet.checkpoint_dir = params['checkpoints_dir']
resnet.log_dir = params['logs_dir']

# start training the network
resnet.train_loop(n_epochs=params['n_epochs'], log_csv=True, verbose=True)

db.close()
resnet.save_model(os.path.join(params['out_dir'],'ketos_model.kt'), audio_repr_file=params['spec_config_file'])
resnet.save_model(os.path.join(params['out_dir'],'ketos_model.ktpb'), audio_repr_file=params['spec_config_file'])

max_stars_repo_name: xaviermouy/minke-whale-dataset

max_issues_repo_name: xaviermouy/minke-whale-dataset

max_forks_repo_name: xaviermouy/minke-whale-dataset


text:
# Copyright (c) 2013-2015 Siphon Contributors.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""
================
NCSS and CartoPy
================

Use Siphon to query the NetCDF Subset Service (NCSS) and plot on a map.

This example uses Siphon's NCSS class to provide temperature data
for contouring a basic map using CartoPy.
"""
from datetime import datetime

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from netCDF4 import num2date
import numpy as np

from siphon.catalog import TDSCatalog

###########################################
# First we construct a `TDSCatalog` instance pointing to our dataset of interest, in
# this case TDS' "Best" virtual dataset for the GFS global 0.25 degree collection of
# GRIB files. This will give us a good resolution for our map. This catalog contains a
# single dataset.
best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
                      'Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
print(list(best_gfs.datasets))

###########################################
# We pull out this dataset and get the NCSS access point
best_ds = best_gfs.datasets[0]
ncss = best_ds.subset()

###########################################
# We can then use the `ncss` object to create a new query object, which
# facilitates asking for data from the server.
query = ncss.query()

###########################################
# We construct a query asking for data corresponding to a latitude and longitude box where 43
# lat is the northern extent, 35 lat is the southern extent, -111 long is the western extent
# and -100 is the eastern extent. We request the data for the current time.
#
# We also ask for NetCDF version 4 data, for the variable 'temperature_surface'. This request
# will return all surface temperatures for points in our bounding box for a single time,
# nearest to that requested. Note the string representation of the query is a properly encoded
# query string.
query.lonlat_box(north=43, south=35, east=-100, west=-111).time(datetime.utcnow())
query.accept('netcdf4')
query.variables('Temperature_surface')

###########################################
# We now request data from the server using this query. The `NCSS` class handles parsing
# this NetCDF data (using the `netCDF4` module). If we print out the variable names, we see
# our requested variable, as well as the coordinate variables (needed to properly reference
# the data).
data = ncss.get_data(query)
print(list(data.variables))

###########################################
# We'll pull out the useful variables for temperature, latitude, and longitude, and time
# (which is the time, in hours since the forecast run).
temp_var = data.variables['Temperature_surface']

# Time variables can be renamed in GRIB collections. Best to just pull it out of the
# coordinates attribute on temperature
time_name = temp_var.coordinates.split()[1]
time_var = data.variables[time_name]
lat_var = data.variables['lat']
lon_var = data.variables['lon']

###########################################
# Now we make our data suitable for plotting.

# Get the actual data values and remove any size 1 dimensions
temp_vals = temp_var[:].squeeze()
lat_vals = lat_var[:].squeeze()
lon_vals = lon_var[:].squeeze()

# Convert the number of hours since the reference time to an actual date
time_val = num2date(time_var[:].squeeze(), time_var.units)

# Convert temps to Fahrenheit from Kelvin
temp_vals = temp_vals * 1.8 - 459.67

# Combine 1D latitude and longitudes into a 2D grid of locations
lon_2d, lat_2d = np.meshgrid(lon_vals, lat_vals)

###########################################
# Now we can plot these up using matplotlib and cartopy.

# Create a new figure
fig = plt.figure(figsize=(15, 12))

# Add the map and set the extent
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-100., -111., 35, 43])

# Add state boundaries to plot
ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=2)

# Contour temperature at each lat/long
cf = ax.contourf(lon_2d, lat_2d, temp_vals, 200, transform=ccrs.PlateCarree(), zorder=0,
                 cmap='coolwarm')

# Plot a colorbar to show temperature and reduce the size of it
plt.colorbar(cf, ax=ax, fraction=0.032)

# Make a title with the time value
ax.set_title(u'Temperature forecast (\u00b0F) for {0:%d %B %Y %H:%MZ}'.format(time_val),
             fontsize=20)

# Plot markers for each lat/long to show grid points for 0.25 deg GFS
ax.plot(lon_2d.flatten(), lat_2d.flatten(), marker='o', color='black', markersize=2,
        alpha=0.3, transform=ccrs.Geodetic(), zorder=2, linestyle='none')

max_stars_repo_name: DanielWatkins/siphon

max_issues_repo_name: DanielWatkins/siphon

max_forks_repo_name: DanielWatkins/siphon


text:
import random
import numpy as np
import datetime

a = random.randint(3,10)
print(a)
b = random.random()
print("%.02f"%b)
print(round(b,2))
c = np.random.randn(3)
print("%.2f"%c[2])
print("{1:.2f}".format(c[1],c[2]))

max_stars_repo_name: LeeBeral/python

max_issues_repo_name: LeeBeral/python

max_forks_repo_name: LeeBeral/python


text:
import numpy as np
import os
from numba import jit
from collections import OrderedDict
import scipy.sparse

class OrderedDefaultListDict(OrderedDict): #name according to default
    def __missing__(self, key):
        self[key] = value = [] #change to whatever default you want
        return value

def scalability_classes(A, method):
    """returns a dictionary with the scalability classes:
the keys are unique couples (k_out, k_in), the objects are list of nodes with the exact (k_out, k_in) of the key.
    :param A: adjacency matrix
    :param method str: denote the method
    :return: scalability classes
    :rtype: dict
    """
    if method == 'dcm_rd':
        k_out = out_degree(A)
        k_in = in_degree(A)
        d = OrderedDefaultListDict()
        tup = tuple([k_out[0], k_in[0]])
        d[tup] = [0]
        n = len(k_out)
        for i in range(1, n):
            # visit each couple (in, out) and add new classes to the dict
            tup = tuple([k_out[i], k_in[i]])
            d[tup].append(i)
        return d

def rd2full(x, d, method):
    """converts a reduced array to full form
    :param x numpy.ndarray: reduced array
    :param d dict: scalability classes
    :param method str: denote the method
    :return: the full extent array
    :rtype: numpy.ndarray
    """
    if method == 'dcm_rd':
        return rd2full_dcm_rd(x, d)

def rd2full_dcm_rd(x, d):
    val = list(d.values())
    n = 0  # dimension of the full solution
    m = len(val)
    for i in range(0, m):  # todo: usare un sum invece del ciclo
        n += len(val[i])
    # allocate full solution
    y1 = np.zeros(n, dtype=x.dtype)
    y2 = np.zeros(n, dtype=x.dtype)
    for i in range(0, m):
        y1[val[i]] = x[i]
        y2[val[i]] = x[i+m]
    return np.hstack((y1, y2))

def setup(A, method):
    """Setup function for the solver.

    Takes in input adjacency matrix and method string
    and returns the parameters array and the initial point

    :param A: adjacency matrix
    :param method str: denote the method
    :return: list of method parameters and solver initial point
    :rtype: list
    """
    if method == 'cm':
        # A should be symmetric!!!
        par = out_degree(A)
        L = A.sum()
        # starting point
        v0 = par/np.sqrt(L)

        return [par, v0]

    if method == 'dcm':
        k_out = out_degree(A)
        k_in = in_degree(A)
        par = np.concatenate((k_out, k_in))
        L = A.sum()
        # starting point
        x = k_out/np.sqrt(L)
        y = k_in/np.sqrt(L)
        v0 = np.concatenate((x, y))

        return [par, v0]

    if method == 'dcm_rd':
        d = scalability_classes(A, method='dcm_rd')
        keys = list(d)
        k_out = np.array([x[0] for x in keys])
        k_in = np.array([x[1] for x in keys])
        c = np.array([len(d[(kout,kin)]) for kout,kin in zip(k_out, k_in)])
        par = np.concatenate((k_out, k_in, c))
        # starting point
        L = A.sum()
        x = k_out/np.sqrt(L)
        y = k_in/np.sqrt(L)
        v0 = np.concatenate((x, y))

        return [par, v0]

    if method == 'rdcm':
        k_out_nr = non_reciprocated_out_degree(A)
        k_in_nr = non_reciprocated_in_degree(A)
        k_r = reciprocated_degree(A)
        par = np.concatenate((k_out_nr, k_in_nr, k_r))
        L = A.sum()
        # starting point
        x = k_out_nr/np.sqrt(L)
        y = k_in_nr/np.sqrt(L)
        z = k_r/np.sqrt(L)
        v0 = np.concatenate((x, y, z))

        return [par, v0]

    if method == 'decm':
        k_out = out_degree(A)
        k_in = in_degree(A)
        s_out = out_strength(A)
        s_in = in_strength(A)
        par = np.concatenate((k_out, k_in, s_out, s_in))
        n = len(k_in)
        L = int(k_in.sum())
        W = int(s_in.sum())
        # starting point
        a_out = k_out/L
        a_in = k_in/L
        b_out = s_out/W
        b_in = s_in/W

        """
        # still to decide the right initial point
        a_out = k_out/np.sqrt(L)
        a_in = k_in/np.sqrt(L)
        b_out = s_out/W
        b_in = s_in/W
        a_out = 0.5*np.ones(len(k_in))
        a_in = 0.5*np.ones(len(k_in))
        b_out = 0.9*np.ones(len(k_in))
        b_in = 0.9*np.ones(len(k_in))
        """
        v0 = np.concatenate((a_out, a_in, b_out, b_in))

        return [par, v0]

@jit(nopython=True)
def iterative_fun_cm(v, par):
    """Return the next iterative step for the Configuration Model.

    :param numpy.ndarray v: old iteration step
    :param numpy.ndarray par: constant parameters of the cm function
    :return: next iteration step
    :rtype: numpy.ndarray

    """

    # problem dimension
    n = len(v)
    x = v
    k = par
    # calculate the denominators
    xd = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                xd[i] += x[j]/(1 + x[i]*x[j])

    # calculate final solutions xx and yy
    xx = k/xd

    return xx

@jit(nopython=True)
def iterative_fun_dcm(v, par):
    """Return the next iterative step for the Directed Configuration Model.

    :param numpy.ndarray v: old iteration step
    :param numpy.ndarray par: constant parameters of the cm function
    :return: next iteration step
    :rtype: numpy.ndarray
    """

    n = int(len(v)/2)
    x = v[0:n]
    y = v[n:2*n]
    k_out = par[0:n]
    k_in = par[n:2*n]
    # calculate the denominators
    xd = np.zeros(n)
    yd = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                xd[i] += y[j]/(1 + x[i]*y[j])
                yd[i] += x[j]/(1 + y[i]*x[j])

    # calculate final solutions xx and yy
    xx = k_out/xd
    yy = k_in/yd

    return np.concatenate((xx, yy))

@jit(nopython=True)
def iterative_fun_dcm_rd(v, par):
    """Return the next iterative step for the Directed Configuration Model Reduced version.

    :param numpy.ndarray v: old iteration step
    :param numpy.ndarray par: constant parameters of the cm function
    :return: next iteration step
    :rtype: numpy.ndarray
    """

    n = int(len(v)/2)
    x = v[0:n]
    y = v[n:2*n]
    k_out = par[0:n]
    k_in = par[n:2*n]
    c = par[2*n:3*n]

    xd = np.zeros(n)
    yd = np.zeros(n)

    for i in range(0, n):
        for j in range(0, n):
            if j != i:
                xd[i] += c[j]*y[j]/(1 + x[i]*y[j])
                yd[i] += c[j]*x[j]/(1 + y[i]*x[j])
            else:
                xd[i] += (c[i] - 1)*y[i]/(1 + x[i]*y[i])
                yd[i] += (c[i] - 1)*x[i]/(1 + x[i]*y[i])

    # calculate final solutions xx and yy
    xx = k_out/xd
    yy = k_in/yd

    return np.concatenate((xx, yy))

@jit(nopython=True)
def iterative_fun_rdcm(v, par):
    """Return the next iterative step for the Reciprocated Directed Configuration Model Reduced version.

    :param numpy.ndarray v: old iteration step
    :param numpy.ndarray par: constant parameters of the cm function
    :return: next iteration step
    :rtype: numpy.ndarray
    """

    # problem dimension
    n = int(len(v)/3)
    x = v[0:n]
    y = v[n:2*n]
    z = v[2*n:3*n]
    k_out_nr = par[0:n]
    k_in_nr = par[n:2*n]
    k_r = par[2*n:3*n]
    # calculate the denominators
    xd = np.zeros(n)
    yd = np.zeros(n)
    zd = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                den = 1 + x[i]*y[j] + x[j]*y[i] \
                    + z[i]*z[j]
                xd[i] += y[j]/den
                yd[i] += x[j]/den
                zd[i] += z[j]/den

    # calculate final solutions xx and yy
    xx = k_out_nr/xd
    yy = k_in_nr/yd
    zz = k_r/zd

    return np.concatenate((xx, yy, zz))

@jit(nopython=True)
def iterative_fun_decm(v, par):
    """Return the next iterative step for the Directed Enhanced Configuration Model Reduced version.

    :param numpy.ndarray v: old iteration step
    :param numpy.ndarray par: constant parameters of the cm function
    :return: next iteration step
    :rtype: numpy.ndarray
    """

   # problem dimension
    n = int(len(v)/4)
    a_out = v[0:n]
    a_in = v[n:2*n]
    b_out = v[2*n:3*n]
    b_in = v[3*n:4*n]
    k_out = par[0:n]
    k_in = par[n:2*n]
    s_out = par[2*n:3*n]
    s_in = par[3*n:4*n]
    # calculate the denominators
    a_out_d = np.zeros(n)
    a_in_d = np.zeros(n)
    b_out_d = np.zeros(n)
    b_in_d = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                """
                print((1 - b_in[j]*b_out[i] + a_in[j]*a_out[i]*b_in[j]*b_out[i]))
                print((1 - b_out[j]*b_in[i] + a_out[j]*a_in[i]*b_out[j]*b_in[i]))
                print(1 - b_in[j]*b_out[i])
                print(1 - b_in[i]*b_out[j])
                """
                a_out_d[i] += a_in[j]*b_in[j]*b_out[i] \
                        /(1 - b_in[j]*b_out[i] \
                        + a_in[j]*a_out[i]*b_in[j]*b_out[i])
                a_in_d[i] += a_out[j]*b_out[j]*b_in[i] \
                        /(1 - b_out[j]*b_in[i] \
                        + a_out[j]*a_in[i]*b_out[j]*b_in[i])
                b_out_d[i] += (a_in[j]*b_in[j]*a_out[i] - b_in[j]) \
                        /(1 - b_in[j]*b_out[i] \
                        + a_in[j]*a_out[i]*b_in[j]*b_out[i]) \
                        + b_in[j]/(1 - b_in[j]*b_out[i])
                b_in_d[i] += (a_out[j]*b_out[j]*a_in[i] - b_out[j]) \
                        /(1 - b_in[i]*b_out[j] \
                        + a_in[i]*a_out[j]*b_in[i]*b_out[j]) \
                        + b_out[j]/(1 - b_in[i]*b_out[j])

        """
        if a_out_d[i] == 0:
            a_out_d[i] = 1
        if a_in_d[i] == 0:
            a_in_d[i] = 1
        if b_out_d[i] == 0:
            b_out_d[i] = 1
        if b_in_d[i] == 0:
            b_in_d[i] = 1
        """
    """
    print(a_out_d)
    print(a_in_d)
    print(b_out_d)
    print(b_in_d)
    """

    # calculate final solutions
    aa_out = k_out/a_out_d
    aa_in = k_in/a_in_d
    bb_out = s_out/b_out_d
    bb_in = s_in/b_in_d

    return np.concatenate((aa_out, aa_in, bb_out, bb_in))

def iterative_solver(A, x0 = None, max_steps = 300, eps = 0.01, method = 'dcm', alfa=1, verbose = False):
    """Return the next iterative step for the Directed Enhanced Configuration Model Reduced version.

    :param numpy.ndarray A: adjacency matrix
    :param int max_steps: maximum number of steps allowed
    :param float eps: solver precision
    :param str method: method the solver implements
    :param float alfa: relax parameter
    :param bool verbose: if True, prints convergence information while running
    :return: model array solution, number of steps for convergence and difference between last two steps
    :rtype: list
    """

    # function choice
    f_dict = {
            'cm' : iterative_fun_cm,
            'dcm' : iterative_fun_dcm,
            'dcm_rd': iterative_fun_dcm_rd,
            'rdcm': iterative_fun_rdcm,
            'decm': iterative_fun_decm
            }
    iterative_fun = f_dict[method]

    # initial setup
    par, v = setup(A, method)
    if x0 != None:
        v = x0

    # verbose
    if verbose == True:
        print('\nProblem parameters = \n{}'.format(par))
        print('\nInitial point  = \n{}'.format(v))

    # iteration steps
    step = 0
    diff = eps + 1
    while diff > eps and step < max_steps:
        # iterative step
        f_step = iterative_fun(v, par)
        dv = f_step - v
        alfa = alfa_choice(v, dv , method=method)
        vv = v + alfa*dv
        old_v = v  # save old step
        # convergence step
        diff = np.linalg.norm(v - vv)/np.linalg.norm(v)  # 2-norm
        # set next step
        v = vv
        step += 1
        # verbose
        if verbose == True:
            print('\n\nstep = {}'.format(step))
            print('\nalfa = {}'.format(alfa))
            print('\nsol = {}'.format(v))
            print('\ndiff = {}'.format(diff))
            # expectation error
            # err = np.nan_to_num(old_v*par/v)
            # print('\nexpectation = {}'.format(err))

    # output
    sol = v

    return sol, step, diff

def out_degree(a):
    #todo : for out_degree and in_degree...check np.int32 is always returned
    """returns matrix A out degrees.

    :param a numpy.ndarray: adjacency matrix
    :return: out degree sequence
    :rtype: numpy.ndarray
    """
    # if the matrix is a numpy array
    if type(a) == np.ndarray:
        return np.sum(a > 0, 1)
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        return np.sum(a > 0, 1).A1

def in_degree(a):
    """returns matrix A in degrees.

    :param a np.ndarray: adjacency matrix
    :return: in degree sequence
    :rtype: numpy.ndarray
    """
    # if the matrix is a numpy array
    if type(a) == np.ndarray:
        return np.sum(a > 0, 0)
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        return np.sum(a > 0, 0).A1

def out_strength(a):
    """returns matrix A out strengths.

    :param a np.ndarray: adjacency matrix
    :return: out strengths sequence
    :rtype: numpy.ndarray
    """

    # if the matrix is a numpy array
    if type(a) == np.ndarray:
        return np.sum(a, 1)
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        return np.sum(a, 1).A1

def in_strength(a):
    """returns matrix A in strengths

    :param a np.ndarray: adjacency matrix
    :return: in stregths sequence
    :rtype: numpy.ndarray
    """

    # if the matrix is a numpy array
    if type(a) == np.ndarray:
        return np.sum(a, 0)
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        return np.sum(a, 0).A1

def non_reciprocated_out_degree(a):
    """returns matrix A non reciprocated in degree

    :param a np.ndarray: adjacency matrix
    :return: out degree non reciprocatedsequence
    :rtype: numpy.ndarray
    """

    if type(a) == np.ndarray:
        s = a.shape
        one = np.ones(shape=s)
        return np.diagonal(a@(one - a))
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        s = a.shape
        one = np.ones(shape=s)
        o = scipy.sparse.csr_matrix(one)
        return (a.dot(o - a)).diagonal()

def non_reciprocated_in_degree(a):
    """returns matrix A non reciprocated in degree

    :param a np.ndarray: adjacency matrix
    :return: in degree non reciprocated sequence
    :rtype: numpy.ndarray
    """

    if type(a) == np.ndarray:
        s = a.shape
        one = np.ones(shape=s)
        return np.diagonal(a.transpose()@(one - a.transpose()))
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        s = a.shape
        one = np.ones(shape=s)
        o = scipy.sparse.csr_matrix(one)
        return (a.transpose().dot(o - a.transpose())).diagonal()

def reciprocated_degree(a):
    """returns matrix A reciprocated degree

    :param a np.ndarray: adjacency matrix
    :return: degree reciprocated sequence
    :rtype: numpy.ndarray
    """

    if type(a) == np.ndarray:
        return np.diagonal(a@a)
    # if the matrix is a scipy sparse matrix
    elif type(a) in [scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix]:
        return (a.dot(a)).diagonal()

def dyads_count(a):
    """Counts number of dyads

    :param a np.ndarray: adjacency matrix
    :return: dyads count
    :rtype: int
    """

    at = a.transpose()
    tmp = a + at
    if isinstance(a, np.ndarray):
        return int(len(tmp[tmp == 2]))
    if isinstance(a, (scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix)):
        return int(tmp[tmp == 2].shape[1])

def singles_count(a):
    """Counts number of singles

    :param a np.ndarray: adjacency matrix
    :return: singles count
    :rtype: int
    """

    at = a.transpose()
    tmp = a + at
    if isinstance(a, np.ndarray):
        return int(len(tmp[tmp == 1])/2)
    if isinstance(a, (scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix)):
        return int(tmp[tmp == 1].shape[1]/2)

def zeros_count(a):
    """Counts number of zeros

    :param a np.ndarray: adjacency matrix
    :return: zeros count
    :rtype: int
    """

    n = a.shape[0]
    at = a.transpose()
    tmp = a + at
    if isinstance(a, np.ndarray):
        return int((n*(n-1) - np.count_nonzero(tmp)))
    if isinstance(a, (scipy.sparse.csr.csr_matrix, scipy.sparse.coo.coo_matrix)):
        return int((n*(n-1) - tmp.count_nonzero()))

def expected_degree(sol, method, directed=None, d=None):
    """returns expected degree after ERGM method

    Parameters
    ----------

    :param sol: :class:`~numpy.ndarray`
        Solution of the ERGM problem
    :param method: :class:`~string`
        String stands for the requested method
    :param directed: :class:`~string`
        If the method is for directed network.
        Accepted values are "In" and "Out"
    :param d: :class:`~dict`
        Scalability map for reduced method

    Returns
    -------

    :return k: :class:`~numpy.ndarray`
        array of expected out-degrees
    """
    # undericted methods
    if method == 'cm':
        return expected_degree_cm(sol)
    # directed methods
    if directed == 'Out':
        if method == 'dcm':
            return expected_out_degree_dcm(sol)
        if method == 'dcm_rd':
            sol_full = rd2full(sol, d, 'dcm_rd')
            return expected_out_degree_dcm(sol_full)
        if method == 'decm':
            return expected_out_degree_decm(sol)
    if directed == 'In':
        if method == 'dcm':
            return expected_in_degree_dcm(sol)
        if method == 'dcm_rd':
            sol_full = rd2full(sol, d, 'dcm_rd')
            return expected_in_degree_dcm(sol_full)
        if method == 'decm':
            return expected_in_degree_decm(sol)

def expected_out_degree(sol, method, d=None):
    # TODO: controllare che funzioni
    """returns expected out degree after ERGM method, on undirected networks
    returns just the expected degree

    Parameters
    ----------

    :param sol: :class:`~numpy.ndarray`
        Solution of the ERGM problem

    Returns
    -------

    :return k: :class:`~numpy.ndarray`
        array of expected out-degrees

    """
    if method == 'dcm':
        return expected_out_degree_dcm(sol)

    if method == 'dcm_rd':
        # cardinality of reduced equivalent classes
        c = [len(d[key]) for key in d.keys()]
        k = expected_out_degree_dcm_rd(sol, c)
        # convert rd to full array
        m = len(d)
        d_vals = list(d.values())
        n = np.array([len(d[x]) for x in d]).sum()
        y = np.zeros(n, dtype=k.dtype)
        for i in range(m):
            y[d_vals[i]] = k[i]

        return y

@jit(nopython=True)
def expected_degree_cm(sol):
    n = len(sol)
    a = sol

    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += a[j]*a[i] / (1 + a[j]*a[i])

    return k

@jit(nopython=True)
def expected_out_degree_dcm(sol):
    n = int(sol.size / 2)
    a_out = sol[0:n]
    a_in = sol[n:]

    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += a_in[j]*a_out[i] / (1 + a_in[j]*a_out[i])

    return k

@jit(nopython=True)
def expected_out_degree_dcm_rd(sol, c):
    n = int(sol.size/2)
    a_out = sol[0:n]
    a_in = sol[n:]
    k = np.zeros(n)  # allocate k

    for i in range(n):
        for j in range(n):
            if j != i:
                k[i] += c[j]*a_in[j]*a_out[i]/(1 + a_in[j]*a_out[i])
            else:
                k[i] += (c[i] - 1)*a_in[i]*a_out[i]/(1 + a_in[i]*a_out[i])
    return k

def expected_in_degree(sol, method, d=None):
    """returns expected in degree after ERGM method, on undirected networks
    returns just the expected degree

    Parameters
    ----------

    :param sol: :class:`~numpy.ndarray`
        Solution of the ERGM problem

    Returns
    -------

    :return k: :class:`~numpy.ndarray`
        array of expected in-degrees
    """
    if method == 'dcm':
        return expected_in_degree_dcm(sol)

    if method == 'dcm_rd':
        # cardinality of scalability classes
        c = [len(d[key]) for key in d.keys()]
        # expected in degree by class
        k = expected_in_degree_dcm_rd(sol, c)
        # convert rd to full array
        m = len(d)
        d_vals = list(d.values())
        n = np.array([len(d[x]) for x in d]).sum()
        y = np.zeros(n, dtype=k.dtype)
        for i in range(m):
            y[d_vals[i]] = k[i]

        return y

@jit(nopython=True)
def expected_in_degree_dcm(sol):
    n = int(sol.size/2)
    a_out = sol[0:n]
    a_in = sol[n:]
    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += a_in[i]*a_out[j]/(1 + a_in[i]*a_out[j])

    return k

@jit(nopython=True)
def expected_in_degree_dcm_rd(sol, c):
    n = int(sol.size/2)
    a_out = sol[0:n]
    a_in = sol[n:]
    k = np.zeros(n)  # allocate k

    for i in range(n):
        for j in range(n):
            if j != i:
                k[i] += c[j]*a_out[j]*a_in[i]/(1 + a_out[j]*a_in[i])
            else:
                k[i] += (c[i] - 1)*a_out[i]*a_in[i]/(1 + a_out[i]*a_in[i])

    return k

@jit(nopython=True)
def expected_non_reciprocated_in_degree_rdcm(sol):
    n = int(sol.size/3)
    x = sol[0:n]
    y = sol[n:2*n]
    z = sol[2*n:3*n]
    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += x[j]*y[i]/(1 + x[i]*y[j] + x[j]*y[i] + z[i]*z[j])

    return k

@jit(nopython=True)
def expected_non_reciprocated_out_degree_rdcm(sol):
    n = int(sol.size/3)
    x = sol[0:n]
    y = sol[n:2*n]
    z = sol[2*n:3*n]
    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += x[i]*y[j]/(1 + x[i]*y[j] + x[j]*y[i] + z[i]*z[j])

    return k

@jit(nopython=True)
def expected_reciprocated_degree_rdcm(sol):
    n = int(sol.size/3)
    x = sol[0:n]
    y = sol[n:2*n]
    z = sol[2*n:3*n]
    k = np.zeros(n)  # allocate k
    for i in range(n):
        for j in range(n):
            if i != j:
                k[i] += z[i]*z[j]/(1 + x[i]*y[j] + x[j]*y[i] + z[i]*z[j])

    return k

@jit(nopython=True)
def expected_out_degree_decm(sol):
    # problem dimension
    n = int(len(sol)/4)
    a_out = sol[0:n]
    a_in = sol[n:2*n]
    b_out = sol[2*n:3*n]
    b_in = sol[3*n:4*n]
    # calculate the denominators
    a_out_d = np.zeros(n)
    a_in_d = np.zeros(n)
    b_out_d = np.zeros(n)
    b_in_d = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                a_out_d[i] += a_in[j]*b_in[j]*b_out[i]/(1 \
                        - b_in[j]*b_out[i] \
                        + a_in[j]*a_out[i]*b_in[j]*b_out[i])

    # calculate final solutions xx and yy
    aa_out = a_out*a_out_d

    return aa_out

@jit(nopython=True)
def expected_in_degree_decm(sol):
    # problem dimension
    n = int(len(sol)/4)
    a_out = sol[0:n]
    a_in = sol[n:2*n]
    b_out = sol[2*n:3*n]
    b_in = sol[3*n:4*n]
    # calculate the denominators
    a_out_d = np.zeros(n)
    a_in_d = np.zeros(n)
    b_out_d = np.zeros(n)
    b_in_d = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                a_in_d[i] += a_out[j]*b_out[j]*b_in[i]/(1 \
                        - b_out[j]*b_in[i] \
                        + a_out[j]*a_in[i]*b_out[j]*b_in[i])

    # calculate final solutions xx and yy
    aa_in = a_in*a_in_d

    return aa_in

def expected_strength(sol, method, directed=None, d=None):
    """returns expected strengths after ERGM method

    Parameters
    ----------

    :param sol: :class:`~numpy.ndarray`
        Solution of the ERGM problem
    :param method: :class:`~string`
        String stands for the requested method
    :param directed: :class:`~string`
        If the method is for directed network.
        Accepted values are "In" and "Out"
    :param d: :class:`~dict`
        Scalability map for reduced method

    Returns
    -------

    :return k: :class:`~numpy.ndarray`
        array of expected strengths
    """
    # undericted methods
    # directed methods
    if directed == 'Out':
        if method == 'decm':
            return expected_out_strength_decm(sol)
    if directed == 'In':
        if method == 'decm':
            return expected_in_strength_decm(sol)

@jit(nopython=True)
def expected_out_strength_decm(sol):
    # problem dimension
    n = int(len(sol)/4)
    a_out = sol[0:n]
    a_in = sol[n:2*n]
    b_out = sol[2*n:3*n]
    b_in = sol[3*n:4*n]
    # calculate the denominators
    a_out_d = np.zeros(n)
    a_in_d = np.zeros(n)
    b_out_d = np.zeros(n)
    b_in_d = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                b_out_d[i] += (a_in[j]*b_in[j]*a_out[i] - b_in[j]) \
                        /(1 - b_in[j]*b_out[i] \
                        + a_in[j]*a_out[i]*b_in[j]*b_out[i]) \
                        + b_in[j]/(1 - b_in[j]*b_out[i])

    # calculate final solutions xx and yy
    bb_out = b_out*b_out_d

    return bb_out

@jit(nopython=True)
def expected_in_strength_decm(sol):
    # problem dimension
    n = int(len(sol)/4)
    a_out = sol[0:n]
    a_in = sol[n:2*n]
    b_out = sol[2*n:3*n]
    b_in = sol[3*n:4*n]
    # calculate the denominators
    a_out_d = np.zeros(n)
    a_in_d = np.zeros(n)
    b_out_d = np.zeros(n)
    b_in_d = np.zeros(n)

    for i in range(n):
        for j in range(n):
            if j != i:
                b_in_d[i] += (a_out[j]*b_in[i]*a_out[j] - b_in[i]) \
                        /(1 - b_in[i]*b_out[j] \
                        + a_in[i]*a_out[j]*b_in[i]*b_out[j]) \
                        + b_out[j]/(1 - b_in[i]*b_out[j])

    # calculate final solutions xx and yy
    bb_in = b_in*b_in_d

    return bb_in

def expected_dyads(sol, method, A=None, t="dyads"):
    """
    Computes the expected number of dyads on the ERGM ensemble
    :param sol: np.ndarray, problem's solution
    :param method: string, the available ERGM methods:
        'dcm':
        'dcm_rd':
    :param d: ordered Dict, contains the info about the reduced system
    :return:
    """
    if method == 'dcm':
        if t == 'dyads':
            return expected_dyads_dcm(sol)
        if t == 'singles':
            return expected_singles_dcm(sol)
        if t == 'zeros':
            return expected_zeros_dcm(sol)

    if method == 'dcm_rd':
        d = scalability_classes(A, 'dcm_rd')
        sol_full = rd2full(sol, d, 'dcm_rd')

        if t == 'dyads':
            return expected_dyads_dcm(sol_full)
        if t == 'singles':
            return expected_singles_dcm(sol_full)
        if t == 'zeros':
            return expected_zeros_dcm(sol_full)

        #TODO check the following commented code and dcm_rd method for dyads
        """
        # cardinality of scalability classes
        c = [len(d[key]) for key in d.keys()]
        # expected in degree by class
        ed = expected_dyads_dcm_rd(sol, c)
        # convert rd to full array
        m = len(d)
        d_vals = list(d.values())
        n = np.array([len(d[x]) for x in d]).sum()
        y = np.zeros(n, dtype=ed.dtype)
        for i in range(m):
            y[d_vals[i]] = ed[i]
        return y
        return y
        """

def std_dyads(sol, method, A=None, t="dyads"):
    """
    Computes the standard deviation of the number of dyads on the ERGM ensemble
    :param sol: np.ndarray, problem's solution
    :param method: string, the available ERGM methods:
        'dcm':
        'dcm_rd':
    :param d: ordered Dict, contains the info about the reduced system
    :return:
    """
    if method == 'dcm':
        if t == 'dyads':
            return std_dyads_dcm(sol)
        if t == 'singles':
            return std_singles_dcm(sol)
        if t == 'zeros':
            return std_zeros_dcm(sol)

    if method == 'dcm_rd':
        d = scalability_classes(A, 'dcm_rd')
        sol_full = rd2full(sol, d, 'dcm_rd')

        if t == 'dyads':
            return std_dyads_dcm(sol_full)
        if t == 'singles':
            return std_singles_dcm(sol_full)
        if t == 'zeros':
            return std_zeros_dcm(sol_full)

@jit(nopython=True)
def expected_dyads_dcm(sol):
    """ compute the expected number of reciprocated links
    """
    # edges
    n = int(len(sol)/2)
    y = sol[:n]
    x = sol[n:]
    er = 0
    for i in range(n):
        temp = 0
        for j in range(n):
            temp += x[j]*y[j]/((1 + x[i]*y[j])*(1 + y[i]*x[j]))
        # i != j should not be accounted
        temp -= x[i]*y[i]/((1 + x[i]*y[i])*(1 + y[i]*x[i]))
        er += x[i]*y[i]*temp
    return er

@jit(nopython=True)
def expected_singles_dcm(sol):
    """ compute the expected number of non reciprocated links
    """
    # edges
    n = int(len(sol)/2)
    y = sol[:n]
    x = sol[n:]
    er = 0
    for i in range(n):
        temp = 0
        for j in range(n):
            temp += y[j]*x[i]/((1 + x[i]*y[j])*(1 + y[i]*x[j]))
        # i != j should not be accounted
        temp -= x[i]*y[i]/((1 + x[i]*y[i])*(1 + y[i]*x[i]))
        er += temp
    return er

@jit(nopython=True)
def expected_zeros_dcm(sol):
    """ compute the expected number of non present links (number of couples not connected)
    """
    # edges
    n = int(len(sol)/2)
    y = sol[:n]
    x = sol[n:]
    er = 0
    for i in range(n):
        temp = 0
        for j in range(n):
            temp += 1/((1 + x[i]*y[j])*(1 + y[i]*x[j]))
        # i != j should not be accounted
        temp -= 1/((1 + x[i]*y[i])*(1 + y[i]*x[i]))
        er += temp
    return er

@jit(nopython=True)
def std_dyads_dcm(sol):
    """ compute the expected number of reciprocated links
    """
    # edges
    n = int(len(sol)/2)
    x = sol[:n]
    y = sol[n:]
    temp = 0
    for i in range(n):
        for j in range(n):
            if i != j:
                pij = x[i]*y[j]/(1 + x[i]*y[j])
                pji = x[j]*y[i]/(1 + x[j]*y[i])
                temp += 2*pij*pji*(1 - pij*pji)
    return np.sqrt(temp)

@jit(nopython=True)
def std_singles_dcm(sol):
    """ compute the expected number of reciprocated links
    """
    # edges
    n = int(len(sol)/2)
    x = sol[:n]
    y = sol[n:]
    temp = 0
    for i in range(n):
        for j in range(n):
            if i != j:
                pij = x[i]*y[j]/(1 + x[i]*y[j])
                pji = x[j]*y[i]/(1 + x[j]*y[i])
                temp += pij*(1 - pji)*(1 - pij*(1 - pji) - pji*(1 - pij))
    return np.sqrt(temp)

@jit(nopython=True)
def std_zeros_dcm(sol):
    """ compute the expected number of zeros couples
    """
    # edges
    n = int(len(sol)/2)
    x = sol[:n]
    y = sol[n:]
    temp = 0
    for i in range(n):
        for j in range(n):
            if i != j:
                pij = x[i]*y[j]/(1 + x[i]*y[j])
                pji = x[j]*y[i]/(1 + x[j]*y[i])
                temp += 2*(1 - pij)*(1 - pji)*(1 - (1 - pij)*(1 - pji))
    return np.sqrt(temp)

@jit(nopython=True)
def expected_dyads_dcm_rd(sol, c):
    #TODO: doesn't work
    n = int(len(sol)/2)
    y = sol[:n] #TODO: tmeporary fix from an old notation
    x = sol[n:]
    er = 0
    for i in range(n):
        temp = 0
        for j in range(n):
            if i != j:
                temp += c[j]*x[j]*y[j]/((1 + x[i]*y[j])*(1 + x[j]*y[i]))
            else:
                temp += (c[j] - 1)*x[j]*y[j] / \
                        ((1 + x[i]*y[j])*(1 + x[j]*y[i]))
        er += c[i]*x[i]*y[i]*temp
    return er

def ensemble_sampler(sol, m, method, sample_dir='.', start=0, seed=None):
    """ samples m adjacency matrices in diretory sampler_dir, after the method ergm solution given by sol
    :param sol numpy.ndarray:
    :param m int:
    :param method str:
    :param sample_dir str:
    :param start int:
    :param seed int:
    :return:
    """
    np.random.seed(seed)
    # if sample_dir doesn't exists, creates it
    try:
        os. mkdir(sample_dir)
    except FileExistsError:
        pass

    if method == 'cm':
        n = len(sol)
        x = sol
        for k in range(start, start + m):
            r = np.random.random((n, n))
            a = np.outer(x, x)/(np.ones((n, n)) + np.outer(x, x))
            np.fill_diagonal(a, 0)
            c = np.zeros((n, n))
            c[a.__gt__(r)] = 1
            del r, a
            sparse_matrix = scipy.sparse.coo_matrix(c)
            del c

            outfile = sample_dir + '/' + 'cm_graph_{}.npz'.format(k)
            scipy.sparse.save_npz(outfile, sparse_matrix)

    if method == 'dcm':
        n = int(len(sol)/2)
        x = sol[:n]
        y = sol[n:]
        # sampling
        for k in range(start, start + m):
            r = np.random.random((n, n))
            a = np.outer(x, y)/(np.ones((n, n)) + np.outer(x, y))
            np.fill_diagonal(a, 0)
            c = np.zeros((n, n))
            c[a.__gt__(r)] = 1
            del r, a
            sparse_matrix = scipy.sparse.coo_matrix(c)
            del c

            outfile = sample_dir + '/' + 'dcm_graph_{}.npz'.format(k)
            scipy.sparse.save_npz(outfile, sparse_matrix)

        return 0

def probability_matrix(sol, method):
    """
    :param sol numpy.ndarray:
    :param method str:
    :return:
    :rtype: numpy.ndarray
    """
    if method == 'dcm':
        n = int(len(sol)/2)
        x = sol[:n]
        y = sol[n:]
        p = np.outer(x, y)/(np.ones((n, n)) + np.outer(x, y))
        np.fill_diagonal(p, 0)
        return p

def alfa_choice(x, dx, eps=1e-2, alfa=0.1, method='decm'):
    """
    :param dv numpy.ndarray:
    :param eps float:
    :param alfa float:
    :param method str:
    :return:
    :rtype: float

    """
    if method == 'decm':
        alfa0 = (eps-1)*x/dx
        for a in alfa0:
            if a>=0:
                alfa = min(alfa, a)
    else:
        alfa = 1

    return alfa

def nearest_neighbour_degree(A, method='undirected'):
    """
    """
    if method == 'undirected':
        knn = nearest_neighbour_degree_undirected(A)

    if method == 'outout':
        knn = nearest_neighbour_degree_outout(A)

    if method == 'inin':
        knn = nearest_neighbour_degree_inin(A)

    return knn

def nearest_neighbour_degree_undirected(A):
    """
    """
    # S = np.maximum( A, A.transpose() )
    S = A
    del A
    rows, cols = S.nonzero()
    S[cols,rows] = S[rows, cols]
    k = out_degree(S)
    v = S.dot(k)
    knn = v/k

    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

def nearest_neighbour_degree_inin(A):
    """
    """
    k = in_degree(A)
    v = A.transpose().dot(k)
    knn = v/k

    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

def nearest_neighbour_degree_outout(A):
    """
    """
    k = out_degree(A)
    v = A.dot(k)
    knn = v/k

    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

def expected_nearest_neighbour_degree_undirected(A, k):
    """
    """
    # S = np.maximum( A, A.transpose() )
    S = A
    del A
    rows, cols = S.nonzero()
    S[cols,rows] = S[rows, cols]
    v = S.dot(k)
    knn = v/k

    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

def expected_dcm_nearest_neighbour_degree_inin(sol):
    """
    """
    # expectation calculation
    n = int(len(sol)/2)
    x = sol[:n]
    y = sol[n:]
    k = expected_in_degree_dcm(sol)
    # v = A.dot(k)
    knn = np.array([ np.sum(y[i]*np.delete(x, i)*np.delete(k, i)/(np.ones(n-1) + y[i]*np.delete(x, i)))/k[i] for i in range(n)])
    print(knn)

    # dict making
    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

def expected_dcm_nearest_neighbour_degree_outout(sol):
    """
    """
    # expectation calculation
    n = int(len(sol)/2)
    x = sol[:n]
    y = sol[n:]
    k = expected_out_degree_dcm(sol)
    # v = A.dot(k)
    knn = np.array([ np.sum(x[i]*np.delete(y, i)*np.delete(k, i)/(np.ones(n-1) + x[i]*np.delete(y, i)))/k[i] for i in range(n)])

    # dict making
    b, c = np.unique(k, return_inverse = True)
    knn_unique = [ knn[np.where(c == i)].sum()/len(np.where(c== i)[0]) for i in range(len(b))]
    k_unique = np.unique(k)

    # return knn
    # return dict with degrees as keys and knn average value as value
    return dict(zip(k_unique, knn_unique))

reason: import numpy,import scipy,from numba