text:
import copy
import re
from datetime import datetime
from os import listdir, remove
from os.path import exists, join
from subprocess import CalledProcessError
from unittest.mock import patch

import f90nml
import numpy as np
import pandas as pd
import pytest
from openscm_units import unit_registry

from pymagicc import MAGICC6, MAGICC7, zero_emissions
from pymagicc.core import MAGICCBase, _clean_value, config
from pymagicc.io import MAGICCData, read_cfg_file

from .test_io import MAGICC6_DIR

@pytest.fixture(scope="function")
def magicc_base():
    yield MAGICCBase()

def write_config(p):
    emis_key = "file_emissionscenario" if p.version == 6 else "FILE_EMISSCEN"
    outpath = join(p.run_dir, "MAGTUNE_SIMPLE.CFG")
    f90nml.write({"nml_allcfgs": {emis_key: "RCP26.SCEN"}}, outpath, force=True)

    # Write years config.
    outpath_years = join(p.run_dir, "MAGCFG_NMLYEARS.CFG")
    f90nml.write(
        {"nml_years": {"startyear": 1765, "endyear": 2100, "stepsperyear": 12}},
        outpath_years,
        force=True,
    )

def test_not_initalise():
    p = MAGICC6()
    assert p.root_dir is None
    assert p.run_dir is None
    assert p.out_dir is None

def test_initalise_and_clean(package):
    # fixture package has already been initialised
    assert exists(package.run_dir)
    assert exists(join(package.run_dir, "MAGCFG_USER.CFG"))
    assert exists(package.out_dir)

def test_run_failure(package):
    if exists(join(package.run_dir, "HISTRCP_CO2I_EMIS.IN")):
        remove(join(package.run_dir, "HISTRCP_CO2I_EMIS.IN"))

    if exists(join(package.run_dir, "HISTSSP_CO2I_EMIS.IN")):
        remove(join(package.run_dir, "HISTSSP_CO2I_EMIS.IN"))

    with pytest.raises(CalledProcessError):
        package.run()

    assert package.config is None

@pytest.mark.parametrize("tonly", (["junk"], ["junk", "junkier"]))
def test_run_no_output(package, tonly):
    error_msg = re.escape("No output found for only={}".format(tonly))
    with pytest.raises(ValueError, match=error_msg):
        package.run(only=tonly)

def test_run_failure_no_magcfg_user(package):
    if exists(join(package.run_dir, "MAGCFG_USER.CFG")):
        remove(join(package.run_dir, "MAGCFG_USER.CFG"))

    with pytest.raises(FileNotFoundError):
        package.run()

    assert package.config is None

@pytest.mark.parametrize(
    "invalid_config",
    [
        {"nml_allcfgs": {"file_tuningmodel_1": "junk"}},
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "pymagicc",
                "file_tuningmodel_2": "pymagicc",
            }
        },
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "pymagicc",
                "file_tuningmodel_2": "c4mip_default",
            }
        },
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "c4mip_default",
                "file_tuningmodel_2": "pymagicc",
            }
        },
    ],
)
def test_run_failure_no_pymagicc_in_magcfg_user(package, invalid_config):
    f90nml.write(invalid_config, join(package.run_dir, "MAGCFG_USER.CFG"), force=True)

    error_msg = re.escape(
        "PYMAGICC is not the only tuning model that will be used by "
        "`MAGCFG_USER.CFG`: your run is likely to fail/do odd things"
    )
    with pytest.raises(ValueError, match=error_msg):
        package.run()

    assert package.config is None

@pytest.mark.parametrize(
    "invalid_config",
    [
        {"nml_allcfgs": {"file_emisscen": "SSP1", "file_emisscen_2": "SSP2"}},
        {
            "nml_allcfgs": {
                "file_emisscen": "SSP1",
                "file_emisscen_2": "",
                "file_emisscen_3": "SSP2",
            }
        },
        {
            "nml_allcfgs": {
                "file_emisscen": "",
                "file_emisscen_2": "SSP1",
                "file_emisscen_3": "SSP2",
            }
        },
    ],
)
def test_run_failure_confusing_emissions_scenarios(package, invalid_config):
    f90nml.write(invalid_config, join(package.run_dir, "MAGCFG_USER.CFG"), force=True)

    error_msg = re.escape(
        "You have more than one `FILE_EMISSCEN_X` flag set. Using more than one "
        "emissions scenario is hard to debug and unnecessary with Pymagicc's "
        "Dataframe scenario input. Please combine all your scenarios into one "
        "Dataframe with Pymagicc and Pandas, then feed this single Dataframe into "
        "Pymagicc's run API."
    )
    with pytest.raises(ValueError, match=error_msg):
        package.run()

    assert package.config is None

@pytest.mark.parametrize(
    "invalid_config",
    [
        {"nml_allcfgs": {"file_emisscen": "SSP1", "file_emisscen_2": "SSP2"}},
        {
            "nml_allcfgs": {
                "file_emisscen": "SSP1",
                "file_emisscen_2": "",
                "file_emisscen_3": "SSP2",
            }
        },
        {
            "nml_allcfgs": {
                "file_emisscen": "",
                "file_emisscen_2": "SSP1",
                "file_emisscen_3": "SSP2",
            }
        },
    ],
)
def test_check_config_non_strict(package_non_strict, invalid_config):
    f90nml.write(
        invalid_config, join(package_non_strict.run_dir, "MAGCFG_USER.CFG"), force=True
    )

    error_msg = re.escape(
        "You have more than one `FILE_EMISSCEN_X` flag set. Using more than one "
        "emissions scenario is hard to debug and unnecessary with Pymagicc's "
        "Dataframe scenario input. Please combine all your scenarios into one "
        "Dataframe with Pymagicc and Pandas, then feed this single Dataframe into "
        "Pymagicc's run API."
    )
    with pytest.warns(UserWarning, match=error_msg):
        package_non_strict.check_config()

@pytest.mark.parametrize(
    "invalid_config",
    [
        {"nml_allcfgs": {"file_tuningmodel_1": "junk"}},
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "pymagicc",
                "file_tuningmodel_2": "pymagicc",
            }
        },
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "pymagicc",
                "file_tuningmodel_2": "c4mip_default",
            }
        },
        {
            "nml_allcfgs": {
                "file_tuningmodel_1": "c4mip_default",
                "file_tuningmodel_2": "pymagicc",
            }
        },
    ],
)
def test_check_config_non_strict_no_pymagicc(package_non_strict, invalid_config):
    f90nml.write(
        invalid_config, join(package_non_strict.run_dir, "MAGCFG_USER.CFG"), force=True
    )

    error_msg = re.escape(
        "PYMAGICC is not the only tuning model that will be used by "
        "`MAGCFG_USER.CFG`: your run is likely to fail/do odd things"
    )
    with pytest.warns(UserWarning, match=error_msg):
        package_non_strict.check_config()

def test_run_success(package):
    write_config(package)
    results = package.run(out_parameters=1)

    assert isinstance(results, MAGICCData)
    assert len(results["variable"].unique()) > 1
    assert "Surface Temperature" in results["variable"].values

    assert (results["climate_model"] == "MAGICC{}".format(package.version)).all()
    # running with preset information should result in unspecified everywhere
    # only running from an IamDataFrame instance should automatically fill model
    # and scenario columns
    assert (results["model"] == "unspecified").all()
    assert (results["scenario"] == "unspecified").all()

    assert len(package.config.keys()) != 0

def test_run_with_magiccdata(package, temp_dir):
    tmodel = "IMAGE"
    tscenario = "RCP26"
    scen = MAGICCData(
        join(MAGICC6_DIR, "RCP26.SCEN"),
        columns={"model": [tmodel], "scenario": [tscenario]},
    )

    results = package.run(scen, only=["Surface Temperature"])

    assert len(results["variable"].unique()) == 1
    assert "Surface Temperature" in results["variable"].values

    assert (results["climate_model"] == "MAGICC{}".format(package.version)).all()
    assert (results["model"] == tmodel).all()
    assert (results["scenario"] == tscenario).all()

def test_run_success_binary(package):
    results = package.run(out_ascii_binary="BINARY", out_keydata_2=True)

    assert isinstance(results, MAGICCData)
    assert len(results["variable"].unique()) > 1
    assert "Surface Temperature" in results["variable"].values

    assert len(package.config.keys()) != 0

def test_run_success_update_config(package):
    package.set_output_variables(
        keydata_2=True, parameters=1, write_ascii=False, write_binary=True
    )
    results = package.run()

    assert isinstance(results, MAGICCData)
    assert len(results["variable"].unique()) > 1
    assert "Surface Temperature" in results["variable"].values

    assert len(package.config.keys()) != 0

def test_run_ambiguous_config(package):
    error_msg = re.escape(
        "The following configuration keys clash because configs are case "
        "insensitive: ['pf_apply', 'PF_APPLY'], "
        "['OUT_TempERATURE', 'out_temperature']"
    )
    with pytest.raises(ValueError, match=error_msg):
        package.run(
            pf_apply=1, PF_APPLY=0, OUT_TempERATURE=0, out_temperature=1, out_forcing=1,
        )

def test_run_only(package):
    write_config(package)
    results = package.run(only=["Surface Temperature"])

    assert len(results["variable"].unique()) == 1
    assert "Surface Temperature" in results["variable"].values

def test_run_rewritten_scen_file(package, temp_dir):
    starting_scen = join(MAGICC6_DIR, "RCP26.SCEN")
    written_scen = join(package.run_dir, "RCP26.SCEN7")

    cols = {"model": ["IMAGE"], "scenario": ["RCP26"], "climate_model": ["MAGICC6"]}
    mdata_initial = MAGICCData(starting_scen, columns=cols)

    mdata_initial.write(written_scen, magicc_version=7)

    mdata_written = MAGICCData(written_scen, columns=cols)
    results = package.run(mdata_written, only=["Surface Temperature"])

    assert len(results["variable"].unique()) == 1
    assert "Surface Temperature" in results["variable"].values

def test_override_config():
    config["EXECUTABLE_6"] = "/tmp/magicc"
    magicc = MAGICC6()

    # Stop this override impacting other tests
    del config.overrides["EXECUTABLE_6"]
    assert magicc.executable == "/tmp/magicc"

def test_dont_create_dir():
    magicc = MAGICC6()
    # Dir isn't created yet
    assert magicc.root_dir is None
    magicc.create_copy()
    root_dir = magicc.root_dir
    assert exists(root_dir)
    magicc.remove_temp_copy()
    assert not exists(root_dir)
    assert magicc.root_dir is None

@pytest.mark.parametrize("strict", (True, False))
def test_overwrite_config(strict):
    magicc = MAGICC7(strict=strict)
    if magicc.executable is None or not exists(magicc.original_dir):
        pytest.skip("MAGICC7 unavailable")

    magicc.create_copy()
    cfgs = magicc.update_config(filename="MAGCFG_USER.CFG")["nml_allcfgs"]
    assert cfgs["file_tuningmodel_1"] == "PYMAGICC"
    assert cfgs["file_tuningmodel_2"] == "USER"
    if strict:
        assert cfgs["file_emisscen_2"] == "NONE"
    else:
        assert "file_emisscen_2" not in cfgs

def test_clean_value_simple():
    assert "SF6" == _clean_value("SF6                 ")

    assert 1970 == _clean_value(1970)
    assert 2012.123 == _clean_value(2012.123)

def test_clean_value_nulls():
    in_str = [
        "SF6                 ",
        "SO2F2               ",
        "\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000",
        "\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000",
    ]
    expected = ["SF6", "SO2F2", "", ""]
    out_str = _clean_value(in_str)

    assert len(out_str) == len(expected)
    for o, e in zip(out_str, expected):
        assert o == e

def test_incorrect_subdir():
    config["EXECUTABLE_6"] = "/tmp/magicc"
    magicc = MAGICC6()
    try:
        with pytest.raises(FileNotFoundError):
            magicc.create_copy()
    finally:
        del config.overrides["EXECUTABLE_6"]
        magicc.remove_temp_copy()

def test_root_dir():
    with MAGICC6() as magicc:
        m2 = MAGICC6(root_dir=magicc.root_dir)

        assert m2.root_dir == magicc.root_dir

        # Does nothing
        m2.remove_temp_copy()
        # Can be called many times
        m2.remove_temp_copy()

        assert m2.root_dir is not None

def test_no_root_dir():
    assert not exists("/tmp/magicc/")
    magicc = MAGICC6(root_dir="/tmp/magicc/")

    with pytest.raises(FileNotFoundError):
        magicc.run()

@patch.object(MAGICCBase, "update_config")
@patch.object(MAGICCBase, "set_years")
def test_diagnose_ecs_config_setup(mock_set_years, mock_update_config, magicc_base):
    magicc_base._diagnose_ecs_config_setup()
    mock_set_years.assert_called_with(startyear=1750, endyear=4200)
    mock_update_config.assert_called_with(
        FILE_CO2_CONC="ABRUPT2XCO2_CO2_CONC.IN",
        CO2_SWITCHFROMCONC2EMIS_YEAR=30000,
        RF_TOTAL_RUNMODUS="CO2",
        RF_TOTAL_CONSTANTAFTERYR=2000,
    )

@patch.object(MAGICCBase, "update_config")
@patch.object(MAGICCBase, "set_years")
def test_diagnose_tcr_tcre_config_setup(
    mock_set_years, mock_update_config, magicc_base
):
    magicc_base._diagnose_tcr_tcre_config_setup()
    mock_set_years.assert_called_with(startyear=1750, endyear=2020)
    mock_update_config.assert_called_with(
        FILE_CO2_CONC="1PCTCO2_CO2_CONC.IN",
        CO2_SWITCHFROMCONC2EMIS_YEAR=30000,
        RF_TOTAL_RUNMODUS="CO2",
        RF_TOTAL_CONSTANTAFTERYR=3000,
        OUT_INVERSEEMIS=1,
    )

def _get_mock_diagnosis_results_as_magiccdata(
    fake_time, fake_co2_concs, fake_rf, fake_temp, fake_inverse_emms=None
):
    fake_regions = ["World"] * len(fake_time)

    mock_co2_conc = pd.DataFrame(
        {
            "time": fake_time,
            "unit": ["ppm"] * len(fake_time),
            "variable": ["Atmospheric Concentrations|CO2"] * len(fake_time),
            "value": fake_co2_concs,
            "region": fake_regions,
            "model": ["N/A"] * len(fake_time),
            "scenario": ["1%/yr_co2"] * len(fake_time),
        }
    )
    mock_rf = pd.DataFrame(
        {
            "time": fake_time,
            "unit": ["W / m^2"] * len(fake_time),
            "variable": ["Radiative Forcing"] * len(fake_time),
            "value": fake_rf,
            "region": fake_regions,
            "model": ["N/A"] * len(fake_time),
            "scenario": ["1%/yr_co2"] * len(fake_time),
        }
    )
    mock_temp = pd.DataFrame(
        {
            "time": fake_time,
            "unit": ["K"] * len(fake_time),
            "variable": ["Surface Temperature"] * len(fake_time),
            "value": fake_temp,
            "region": fake_regions,
            "model": ["N/A"] * len(fake_time),
            "scenario": ["1%/yr_co2"] * len(fake_time),
        }
    )

    to_concat = [mock_co2_conc, mock_rf, mock_temp]

    if fake_inverse_emms is not None:
        mock_inverse_emms = pd.DataFrame(
            {
                "time": fake_time,
                "unit": ["GtC/yr"] * len(fake_time),
                "variable": ["Inverse Emissions|CO2|MAGICC Fossil and Industrial"]
                * len(fake_time),
                "value": fake_inverse_emms,
                "region": fake_regions,
                "model": ["N/A"] * len(fake_time),
                "scenario": ["1%/yr_co2"] * len(fake_time),
            }
        )

        to_concat.append(mock_inverse_emms)

    mock_results = MAGICCData(pd.concat(to_concat))

    return mock_results

@pytest.fixture
def valid_ecs_diagnosis_results():
    startyear = 1700
    endyear = 4000
    spin_up_time = 50
    ecs_start_time = startyear + spin_up_time
    ecs_yr = endyear
    fake_PI_conc = 278.0

    fake_time = np.arange(startyear, endyear + 1)
    fake_concs = np.concatenate(
        (
            fake_PI_conc * np.ones(spin_up_time),
            fake_PI_conc * 2 * np.ones(endyear - startyear - spin_up_time + 1),
        )
    )
    fake_rf = 2.0 * np.log(fake_concs / fake_PI_conc)
    fake_temp = np.log(fake_rf + 1.0) + fake_time / fake_time[1400]

    mock_results = _get_mock_diagnosis_results_as_magiccdata(
        fake_time, fake_concs, fake_rf, fake_temp
    )

    yield {
        "mock_results": mock_results,
        "ecs_time": datetime(ecs_yr, 1, 1),
        "ecs_start_time": datetime(ecs_start_time, 1, 1),
    }

@pytest.fixture
def valid_tcr_tcre_diagnosis_results():
    startyear = 1700
    endyear = 1960
    spin_up_time = 50
    tcre_start_yr = startyear + spin_up_time
    tcr_yr = tcre_start_yr + 70

    fake_PI_conc = 278.0

    fake_time = np.arange(startyear, endyear + 1)
    fake_concs = np.concatenate(
        (
            fake_PI_conc * np.ones(spin_up_time),
            fake_PI_conc * 1.01 ** (np.arange(endyear - startyear - spin_up_time + 1)),
        )
    )
    fake_inverse_emms = np.linspace(0, 60, endyear - startyear + 1)
    fake_rf = 2.0 * np.log(fake_concs / fake_PI_conc)
    fake_temp = np.log(fake_rf + 1.0) + fake_time / fake_time[len(fake_time) // 3]

    mock_results = _get_mock_diagnosis_results_as_magiccdata(
        fake_time, fake_concs, fake_rf, fake_temp, fake_inverse_emms=fake_inverse_emms
    )

    yield {
        "mock_results": mock_results,
        "tcr_time": datetime(tcr_yr, 1, 1),
        "tcr_start_time": datetime(tcre_start_yr, 1, 1),
    }

@patch.object(MAGICCBase, "_check_ecs_temp")
@patch.object(MAGICCBase, "_check_ecs_total_RF")
@patch.object(MAGICCBase, "_get_ecs_ecs_start_yr_from_CO2_concs")
def test_get_ecs_from_diagnosis_results(
    mock_get_ecs_ecs_start_yr_from_CO2_concs,
    mock_check_ecs_total_RF,
    mock_check_ecs_temp,
    valid_ecs_diagnosis_results,
    magicc_base,
):
    test_ecs_time = valid_ecs_diagnosis_results["ecs_time"]
    test_ecs_start_time = valid_ecs_diagnosis_results["ecs_start_time"]
    test_results_df = valid_ecs_diagnosis_results["mock_results"]

    mock_get_ecs_ecs_start_yr_from_CO2_concs.return_value = (
        test_ecs_time,
        test_ecs_start_time,
    )

    expected_ecs = (
        test_results_df.filter(variable="Surface Temperature", time=test_ecs_time)
        .timeseries()
        .squeeze()
    ) * unit_registry("K")

    actual_ecs = magicc_base.get_ecs_from_diagnosis_results(test_results_df)
    assert actual_ecs == expected_ecs

    assert mock_get_ecs_ecs_start_yr_from_CO2_concs.call_count == 1
    assert mock_check_ecs_total_RF.call_count == 1
    assert mock_check_ecs_temp.call_count == 1

@patch.object(MAGICCBase, "_check_tcr_tcre_temp")
@patch.object(MAGICCBase, "_check_tcr_tcre_total_RF")
@patch.object(MAGICCBase, "_get_tcr_tcr_start_yr_from_CO2_concs")
def test_get_tcr_tcre_from_diagnosis_results(
    mock_get_tcr_tcr_start_yr_from_CO2_concs,
    mock_check_tcr_tcre_total_RF,
    mock_check_tcr_tcre_temp,
    valid_tcr_tcre_diagnosis_results,
    magicc_base,
):
    test_tcr_time = valid_tcr_tcre_diagnosis_results["tcr_time"]
    test_tcre_start_time = valid_tcr_tcre_diagnosis_results["tcr_start_time"]
    test_results_df = valid_tcr_tcre_diagnosis_results["mock_results"]

    mock_get_tcr_tcr_start_yr_from_CO2_concs.return_value = [
        test_tcr_time,
        test_tcre_start_time,
    ]

    expected_tcr = (
        test_results_df.filter(variable="Surface Temperature", time=test_tcr_time)
        .timeseries()
        .squeeze()
    ) * unit_registry("K")
    expected_tcre_cumulative_co2 = test_results_df.filter(
        variable="Inverse Emissions|CO2|MAGICC Fossil and Industrial",
        year=range(test_tcre_start_time.year, test_tcr_time.year),
    ).values.sum() * unit_registry("GtC")
    expected_tcre = expected_tcr / expected_tcre_cumulative_co2

    (actual_tcr, actual_tcre,) = magicc_base.get_tcr_tcre_from_diagnosis_results(
        test_results_df
    )
    assert actual_tcr == expected_tcr
    assert actual_tcre == expected_tcre

    assert mock_get_tcr_tcr_start_yr_from_CO2_concs.call_count == 1
    assert mock_check_tcr_tcre_total_RF.call_count == 1
    assert mock_check_tcr_tcre_temp.call_count == 1

def _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
    base_data, method_to_run, regexp_to_match, *args, test_target="other"
):
    test_time = base_data["time"].values
    times_to_break = [
        test_time[3],
        test_time[15],
        test_time[115],
        test_time[-100],
        test_time[-1],
    ]
    if test_target != "temperature":
        times_to_break.append(test_time[0])
    for time_to_break in times_to_break:
        broken_data = base_data.timeseries()
        col_to_adjust = broken_data.columns == time_to_break
        if test_target == "temperature":
            broken_data.iloc[:, col_to_adjust] -= 0.1
        elif test_target == "rf-tcr":
            broken_data.iloc[:, col_to_adjust] -= 0.1
        else:
            broken_data.iloc[:, col_to_adjust] *= 1.01
            broken_data.iloc[:, col_to_adjust] += 0.01
            broken_data.iloc[-1, col_to_adjust] += 0.1
        with pytest.raises(ValueError, match=regexp_to_match):
            method_to_run(MAGICCData(broken_data), *args)

def test_get_ecs_ecs_start_yr_from_CO2_concs(valid_ecs_diagnosis_results, magicc_base):
    test_results_df = valid_ecs_diagnosis_results["mock_results"]
    test_CO2_data = test_results_df.filter(variable="Atmospheric Concentrations|CO2")

    (
        actual_ecs_yr,
        actual_ecs_start_yr,
    ) = magicc_base._get_ecs_ecs_start_yr_from_CO2_concs(test_CO2_data)
    assert actual_ecs_yr == valid_ecs_diagnosis_results["ecs_time"]
    assert actual_ecs_start_yr == valid_ecs_diagnosis_results["ecs_start_time"]

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_CO2_data,
        magicc_base._get_ecs_ecs_start_yr_from_CO2_concs,
        r"The ECS CO2 concs look wrong.*",
    )

def test_get_tcr_tcr_start_yr_from_CO2_concs(
    valid_tcr_tcre_diagnosis_results, magicc_base
):
    test_results_df = valid_tcr_tcre_diagnosis_results["mock_results"]
    test_CO2_data = test_results_df.filter(variable="Atmospheric Concentrations|CO2")

    (
        actual_tcr_yr,
        actual_tcr_start_yr,
    ) = magicc_base._get_tcr_tcr_start_yr_from_CO2_concs(test_CO2_data)
    assert actual_tcr_yr == valid_tcr_tcre_diagnosis_results["tcr_time"]
    assert actual_tcr_start_yr == valid_tcr_tcre_diagnosis_results["tcr_start_time"]

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_CO2_data,
        magicc_base._get_tcr_tcr_start_yr_from_CO2_concs,
        r"The TCR/TCRE CO2 concs look wrong.*",
    )

def test_check_ecs_total_RF(valid_ecs_diagnosis_results, magicc_base):
    test_results_df = valid_ecs_diagnosis_results["mock_results"]
    test_RF_data = test_results_df.filter(variable="Radiative Forcing")
    magicc_base._check_ecs_total_RF(
        test_RF_data, valid_ecs_diagnosis_results["ecs_start_time"],
    )

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_RF_data,
        magicc_base._check_ecs_total_RF,
        r"The ECS total radiative forcing looks wrong.*",
        valid_ecs_diagnosis_results["ecs_start_time"],
    )

def test_check_tcr_tcre_total_RF(valid_tcr_tcre_diagnosis_results, magicc_base):
    test_results_df = valid_tcr_tcre_diagnosis_results["mock_results"]
    test_RF_data = test_results_df.filter(variable="Radiative Forcing")
    magicc_base._check_tcr_tcre_total_RF(
        test_RF_data, valid_tcr_tcre_diagnosis_results["tcr_time"],
    )

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_RF_data,
        magicc_base._check_tcr_tcre_total_RF,
        r"The TCR/TCRE total radiative forcing looks wrong.*",
        valid_tcr_tcre_diagnosis_results["tcr_time"],
        test_target="rf-tcr",
    )

def test_check_ecs_temp(valid_ecs_diagnosis_results, magicc_base):
    test_results_df = valid_ecs_diagnosis_results["mock_results"]
    test_temp_data = test_results_df.filter(variable="Surface Temperature")

    magicc_base._check_ecs_temp(test_temp_data)

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_temp_data,
        magicc_base._check_ecs_temp,
        r"The ECS surface temperature looks wrong, it decreases",
        test_target="temperature",
    )

def test_check_tcr_tcre_temp(valid_tcr_tcre_diagnosis_results, magicc_base):
    test_results_df = valid_tcr_tcre_diagnosis_results["mock_results"]
    test_temp_data = test_results_df.filter(variable="Surface Temperature")

    magicc_base._check_tcr_tcre_temp(test_temp_data)

    _assert_bad_tcr_ecs_tcre_diagnosis_values_caught(
        test_temp_data,
        magicc_base._check_tcr_tcre_temp,
        r"The TCR/TCRE surface temperature looks wrong, it decreases",
        test_target="temperature",
    )

# integration test (i.e. actually runs magicc) hence slow
@pytest.mark.slow
def test_integration_diagnose_tcr_ecs_tcre(package):
    actual_result = package.diagnose_tcr_ecs_tcre()
    assert isinstance(actual_result, dict)
    assert "tcr" in actual_result
    assert "ecs" in actual_result
    assert "tcre" in actual_result
    assert actual_result["tcr"] < actual_result["ecs"]
    if isinstance(package, MAGICC6):
        # MAGICC6 shipped with pymagicc should be stable
        np.testing.assert_allclose(
            actual_result["tcr"].to("K").magnitude, 1.9733976, rtol=1e-5
        )
        np.testing.assert_allclose(
            actual_result["ecs"].to("K").magnitude, 2.98326, rtol=1e-5
        )
        np.testing.assert_allclose(
            actual_result["tcre"].to("K / TtC").magnitude, 2.28698, rtol=1e-5
        )

    if isinstance(package, MAGICC7):
        # see how stable this is, can delete the test later if it's overly restrictive
        np.testing.assert_allclose(
            actual_result["tcr"].to("K").magnitude, 1.982697, rtol=1e-5
        )
        np.testing.assert_allclose(
            actual_result["ecs"].to("K").magnitude, 2.9948422, rtol=1e-5
        )
        np.testing.assert_allclose(
            actual_result["tcre"].to("K / TtC").magnitude, 2.3189736, rtol=1e-5
        )

@patch.object(MAGICCBase, "_diagnose_ecs_config_setup")
@patch.object(MAGICCBase, "run")
@patch.object(MAGICCBase, "get_ecs_from_diagnosis_results")
def test_diagnose_ecs(
    mock_get_ecs_from_results,
    mock_run,
    mock_diagnose_ecs_setup,
    valid_ecs_diagnosis_results,
    magicc_base,
):
    mock_ecs_val = 3.1 * unit_registry("K")
    mock_run_results = valid_ecs_diagnosis_results["mock_results"]

    mock_run.return_value = mock_run_results
    mock_get_ecs_from_results.return_value = mock_ecs_val

    assert magicc_base.diagnose_ecs()["ecs"] == mock_ecs_val
    assert mock_diagnose_ecs_setup.call_count == 1
    mock_run.assert_called_with(
        only=[
            "Atmospheric Concentrations|CO2",
            "Radiative Forcing",
            "Surface Temperature",
        ],
        scenario=None,
    )
    assert mock_get_ecs_from_results.call_count == 1

    assert magicc_base.diagnose_ecs()["ecs"] == mock_ecs_val
    assert mock_diagnose_ecs_setup.call_count == 2
    assert mock_get_ecs_from_results.call_count == 2

    full_results = magicc_base.diagnose_ecs()
    assert isinstance(full_results, dict)

@patch.object(MAGICCBase, "_diagnose_tcr_tcre_config_setup")
@patch.object(MAGICCBase, "run")
@patch.object(MAGICCBase, "get_tcr_tcre_from_diagnosis_results")
def test_diagnose_tcr_tcre(
    mock_get_tcr_tcre_from_results,
    mock_run,
    mock_diagnose_tcr_tcre_setup,
    valid_tcr_tcre_diagnosis_results,
    magicc_base,
):
    mock_tcr_val = 1.8 * unit_registry("K")
    mock_tcre_val = 2.5 * unit_registry("K / GtC")
    mock_run_results = valid_tcr_tcre_diagnosis_results["mock_results"]

    mock_run.return_value = mock_run_results
    mock_get_tcr_tcre_from_results.return_value = [
        mock_tcr_val,
        mock_tcre_val,
    ]

    assert magicc_base.diagnose_tcr_tcre()["tcr"] == mock_tcr_val
    assert mock_diagnose_tcr_tcre_setup.call_count == 1
    mock_run.assert_called_with(
        only=[
            "Atmospheric Concentrations|CO2",
            "INVERSEEMIS",
            "Radiative Forcing",
            "Surface Temperature",
        ],
        scenario=None,
    )
    assert mock_get_tcr_tcre_from_results.call_count == 1

    assert magicc_base.diagnose_tcr_tcre()["tcre"] == mock_tcre_val
    assert mock_diagnose_tcr_tcre_setup.call_count == 2
    assert mock_get_tcr_tcre_from_results.call_count == 2

    full_results = magicc_base.diagnose_tcr_tcre()
    assert isinstance(full_results, dict)

def test_missing_config(config_override):
    with MAGICC6():
        pass
    config_override("EXECUTABLE_6", "")
    with pytest.raises(FileNotFoundError):
        with MAGICC6():
            pass

    config_override("EXECUTABLE_6", "/invalid/path")
    with pytest.raises(FileNotFoundError):
        with MAGICC6():
            pass

    config_override("EXECUTABLE_7", "")
    with pytest.raises(FileNotFoundError):
        with MAGICC7():
            pass

def test_read_parameters(package):
    with type(package)() as magicc:
        # parameters don't exist
        with pytest.raises(FileNotFoundError):
            magicc.read_parameters()

        # Don't read config if it doesn't exist
        magicc.set_config(out_parameters=0)
        magicc.run()

        assert magicc.config is None

    with type(package)() as magicc:
        magicc.run()
        assert isinstance(magicc.config, dict)
        assert "allcfgs" in magicc.config

def test_updates_namelist(package):
    write_config(package)

    fname = join(package.run_dir, "MAGTUNE_SIMPLE.CFG")
    raw_conf = f90nml.read(fname)
    assert "test_value" not in raw_conf["nml_allcfgs"]

    package.update_config("MAGTUNE_SIMPLE.CFG", test_value=1.2)

    updated_conf = f90nml.read(fname)
    assert "test_value" in updated_conf["nml_allcfgs"]

def test_updates_namelist_missing(package):
    fname = join(package.run_dir, "MAGTUNE_NOTEXISTS.CFG")

    assert not exists(fname)

    package.update_config("MAGTUNE_NOTEXISTS.CFG", test_value=1.2)

    updated_conf = f90nml.read(fname)
    assert "test_value" in updated_conf["nml_allcfgs"]

def test_update_config_conflicting_keys_error(package):
    error_msg = re.escape(
        "The following configuration keys clash because configs are case "
        "insensitive: ['core_climatesensitivity', 'CORE_CLIMATESENSITIVITY']"
    )
    with pytest.raises(ValueError, match=error_msg):
        package.update_config(core_climatesensitivity=3, CORE_CLIMATESENSITIVITY=3)

def test_set_config_conflicting_keys_error(package):
    error_msg = re.escape(
        "The following configuration keys clash because configs are case "
        "insensitive: ['core_climatesensitivity', 'CORE_CLIMATESENSITIVITY']"
    )
    with pytest.raises(ValueError, match=error_msg):
        package.set_config(core_climatesensitivity=3, CORE_CLIMATESENSITIVITY=2)

def test_ascii_output(package):
    fname = join(package.run_dir, "MAGTUNE_PYMAGICC.CFG")

    package.set_output_variables(write_ascii=True, write_binary=True)
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_ASCII_BINARY"] == "BOTH"

    package.set_output_variables(write_ascii=False, write_binary=True)
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_ASCII_BINARY"] == "BINARY"

    package.set_output_variables()
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_ASCII_BINARY"] == "ASCII"

    with pytest.raises(AssertionError):
        package.set_output_variables(write_ascii=False, write_binary=False)

def test_output_variables(package):
    fname = join(package.run_dir, "MAGTUNE_PYMAGICC.CFG")

    package.set_output_variables()
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_TEMPERATURE"] == 0

    package.set_output_variables(temperature=True)
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_TEMPERATURE"] == 1

    # Even accepts invalid variable names
    package.set_output_variables(this_doesnt_exist=False)
    raw_conf = f90nml.read(fname)
    assert raw_conf["nml_allcfgs"]["OUT_THIS_DOESNT_EXIST"] == 0

def test_persistant_state(package):
    test_ecs = 4.3
    package.update_config(CORE_CLIMATESENSITIVITY=test_ecs)

    actual_results = package.run(out_parameters=1)
    assert (
        actual_results.metadata["parameters"]["allcfgs"]["core_climatesensitivity"]
        == test_ecs
    )

# TODO: move to integration tests folder
@pytest.mark.parametrize("rf_method", ["OLBL", "IPCCTAR"])
def test_persistant_state_integration(package, rf_method):
    test_ecs = 1.75
    package.update_config(CORE_CLIMATESENSITIVITY=test_ecs)
    if package.version == 7:
        # MAGICC6 doesn't have this flag so we don't need to adjust anything (yes, we
        # run the test twice, small price to pay for clarity and scalability)
        package.update_config(CORE_CO2CH4N2O_RFMETHOD=rf_method)
        if rf_method == "OLBL":
            pytest.xfail(
                "MAGICC7's ECS only follows its definition with IPCCTAR forcing"
            )

    actual_results = package.diagnose_tcr_ecs_tcre()
    np.testing.assert_allclose(actual_results["ecs"], test_ecs, rtol=1e-02)

# TODO: move to integration tests folder
@pytest.mark.parametrize(
    "test_filename, relevant_config, outputs_to_check, time_check_min, time_check_max",
    [
        (
            "HISTRCP_N2OI_EMIS.IN",
            {
                "file_n2oi_emis": "test_filename",
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            ["Emissions|N2O|MAGICC Fossil and Industrial"],
            1,
            1999,
        ),
        (
            "HISTRCP_CH4_CONC.IN",
            {
                "file_ch4_conc": "test_filename",
                "out_concentrations": 1,
                "ch4_switchfromconc2emis_year": 2010,
            },
            ["Atmospheric Concentrations|CH4"],
            1,
            1999,
        ),
        (
            "RCP26.SCEN",
            {
                "file_emissionscenario": "test_filename",
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            [
                "Emissions|CO2|MAGICC Fossil and Industrial",
                "Emissions|CO2|MAGICC AFOLU",
                "Emissions|CH4|MAGICC Fossil and Industrial",
                "Emissions|BC|MAGICC Fossil and Industrial",
                "Emissions|SOx|MAGICC AFOLU",
                "Emissions|HFC32",
            ],
            2030,
            20000,
        ),
        (
            "RCP26.SCEN7",
            {
                "file_emissionscenario": "test_filename",
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            [
                "Emissions|CO2|MAGICC Fossil and Industrial",
                "Emissions|CO2|MAGICC AFOLU",
                "Emissions|CH4|MAGICC Fossil and Industrial",
                "Emissions|BC|MAGICC Fossil and Industrial",
                "Emissions|SOx|MAGICC AFOLU",
                "Emissions|HFC32",
            ],
            2030,
            20000,
        ),
        (
            "SRESA2.SCEN",  # all other SRES scenarios have been removed from MAGICC7's run directory
            {
                "file_emisscen": "test_filename",  # use MAGICC7 flag here, should still pass
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            [
                "Emissions|CO2|MAGICC Fossil and Industrial",
                "Emissions|CO2|MAGICC AFOLU",
                "Emissions|CH4|MAGICC Fossil and Industrial",
                "Emissions|SOx|MAGICC AFOLU",
                "Emissions|HFC32",
            ],
            2030,
            20000,
        ),
        (
            "RCPODS_WMO2006_Emissions_A1.prn",
            {
                "file_mhalo_emis": "test_filename",
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            ["Emissions|CFC11", "Emissions|CCl4", "Emissions|CH3Cl"],
            1,
            20000,
        ),
        (
            "RCPODS_WMO2006_MixingRatios_A1.prn",
            {
                "file_mhalo_conc": "test_filename",
                "out_concentrations": 1,
                "scen_histadjust_0no1scale2shift": 0,
                "mhalo_switch_conc2emis_yr": 20000,
            },
            [
                "Atmospheric Concentrations|CFC11",
                "Atmospheric Concentrations|CCl4",
                "Atmospheric Concentrations|CH3Cl",
            ],
            1,
            20000,
        ),
        (
            "HISTRCP85_SOLAR_RF.IN",
            {"file_solar_rf": "test_filename", "out_forcing": 1},
            ["Radiative Forcing|Solar"],
            1,
            20000,
        ),
    ],
)
def test_pymagicc_writing_has_an_effect(
    package,
    test_filename,
    relevant_config,
    outputs_to_check,
    time_check_min,
    time_check_max,
):
    if ("SCEN" in test_filename) and (package.version == 7):
        # special undocumented flags!!!
        relevant_config["fgas_adjstfutremis2past_0no1scale"] = 0
        relevant_config["mhalo_adjstfutremis2past_0no1scale"] = 0

    iter_dict = copy.deepcopy(relevant_config)
    for key, value in iter_dict.items():
        if value == "test_filename":
            relevant_config[key] = test_filename
        # Handle adjustment to `.prn` handling in MAGICC
        if key == "file_mhalo_emis" and package.version == 7:
            relevant_config["mhalo_prnfile_emis"] = relevant_config.pop(key)
            relevant_config["mhalo_take_prnfile"] = 1
        if key == "file_mhalo_conc" and package.version == 7:
            relevant_config["mhalo_prnfile_conc"] = relevant_config.pop(key)
            relevant_config["mhalo_take_prnfile"] = 1

    package.set_config(conflict="ignore", **relevant_config)

    if (package.version == 6) and test_filename.endswith("SCEN7"):
        error_msg = re.compile("MAGICC6 cannot run SCEN7 files")
        with pytest.raises(ValueError, match=error_msg):
            package.run(only=outputs_to_check)
        return

    if ("SRES" in test_filename) and (package.version == 7):
        # MAGICC7 cannot run SRES SCEN files
        with pytest.raises(CalledProcessError):
            package.run(only=outputs_to_check)
        return

    if ".prn" in test_filename and package.version == 7:
        # MAGICC7's prn handling is not working
        with pytest.raises(CalledProcessError):
            package.run(only=outputs_to_check)
        return

    initial_results = package.run(only=outputs_to_check)

    ttweak_factor = 0.9

    mdata = MAGICCData(
        join(package.run_dir, test_filename),
        columns={"model": ["unspecified"], "scenario": ["unspecified"]},
    )
    mdata = mdata * ttweak_factor
    mdata.write(join(package.run_dir, test_filename), package.version)

    tweaked_results = package.run(only=outputs_to_check)

    for output_to_check in outputs_to_check:
        result = (
            tweaked_results.filter(
                variable=output_to_check, year=range(time_check_min, time_check_max + 1)
            )
            .timeseries()
            .values
        )
        initial = (
            initial_results.filter(
                variable=output_to_check, year=range(time_check_min, time_check_max + 1)
            )
            .timeseries()
            .values
        )
        expected = ttweak_factor * initial

        abstol = np.max([result, expected]) * 10 ** -3
        np.testing.assert_allclose(result, expected, rtol=1e-5, atol=abstol)

# TODO: move to integration tests folder
@pytest.mark.parametrize(
    "test_filename, relevant_config, outputs_to_check",
    [
        (
            "RCP26.SCEN",
            {
                "file_emisscen": "test_filename",
                "out_emissions": 1,
                "scen_histadjust_0no1scale2shift": 0,
            },
            [("Emissions|C2F6", "World", 2050, "kt C2F6 / yr", 0.4712)],
        )
    ],
)
def test_pymagicc_writing_compatibility_203(
    package, test_filename, relevant_config, outputs_to_check
):
    if ("SCEN" in test_filename) and (package.version == 7):
        # special undocumented flags!!!
        relevant_config["fgas_adjstfutremis2past_0no1scale"] = 0
        relevant_config["mhalo_adjstfutremis2past_0no1scale"] = 0

    for key, value in relevant_config.items():
        if value == "test_filename":
            relevant_config[key] = test_filename

    package.set_config(conflict="ignore", **relevant_config)
    results = package.run(
        out_emissions=1, out_ascii_binary="ASCII"  # ensure units included
    )

    for output_to_check in outputs_to_check:
        expected = output_to_check[-1]
        result = (
            results.filter(
                variable=output_to_check[0],
                region=output_to_check[1],
                year=output_to_check[2],
                unit=output_to_check[3],
            )
            .timeseries()
            .squeeze()
        )

        assert expected == result

def test_zero_run(package):
    package.set_zero_config()
    vars_to_check = ["Surface Temperature", "Radiative Forcing"]
    results = package.run(only=vars_to_check, endyear=2500)
    for var in vars_to_check:
        np.testing.assert_allclose(results.filter(variable=var).timeseries().values, 0)

def test_external_forcing_only_run(package):
    time = zero_emissions["time"]

    forcing_external = 2.0 * np.arange(0, len(time)) / len(time)
    forcing_ext = MAGICCData(
        forcing_external,
        index=time,
        columns={
            "scenario": ["idealised"],
            "model": ["unspecified"],
            "climate_model": ["unspecified"],
            "variable": ["Radiative Forcing|Extra"],
            "unit": ["W / m^2"],
            "todo": ["SET"],
            "region": ["World"],
        },
    )
    forcing_ext_filename = "CUSTOM_EXTRA_RF.IN"
    forcing_ext.metadata = {"header": "External radiative forcing file for testing"}
    forcing_ext.write(join(package.run_dir, forcing_ext_filename), package.version)

    results = package.run(
        rf_extra_read=1,
        file_extra_rf=forcing_ext_filename,
        rf_total_runmodus="QEXTRA",
        endyear=max(time).year,
        rf_initialization_method="ZEROSTARTSHIFT",  # this is default but just in case
        rf_total_constantafteryr=5000,
    )

    # MAGICC's weird last year business means that last result is just constant from previous
    # year and is not treated properly
    # TODO: add this in docs
    validation_output = (
        results.filter(variable="Radiative Forcing", region="World")
        .timeseries()
        .values.squeeze()[:-1]
    )
    validation_input = forcing_external[:-1]
    np.testing.assert_allclose(validation_input, validation_output, rtol=1e-5)
    temperature_global = (
        results.filter(variable="Surface Temperature", region="World")
        .timeseries()
        .values.squeeze()
    )
    assert (temperature_global[1:] - temperature_global[:-1] >= 0).all()

def test_co2_emissions_only(package):
    df = zero_emissions.timeseries()

    time = zero_emissions["time"]
    emms_fossil_co2 = np.linspace(0, 30, len(time))
    df.loc[
        (
            df.index.get_level_values("variable")
            == "Emissions|CO2|MAGICC Fossil and Industrial"
        ),
        :,
    ] = emms_fossil_co2

    scen = MAGICCData(df)
    scen.metadata["header"] = "Test CO2 emissions only file"
    results = package.run(
        scen,
        endyear=max(scen["time"]).year,
        rf_total_constantafteryr=5000,
        rf_total_runmodus="CO2",
        co2_switchfromconc2emis_year=min(scen["time"]).year - 1,
        out_emissions=1,
        only=["Surface Temperature", "Emissions|CO2|MAGICC Fossil and Industrial"],
    )

    output_co2 = (
        results.filter(variable="Em*CO2*Fossil*", region="World")
        .timeseries()
        .values.squeeze()
    )
    assert not (output_co2 == 0).all()
    np.testing.assert_allclose(output_co2, emms_fossil_co2, rtol=0.0005)

    temperature_global = (
        results.filter(variable="Surface Temperature", region="World")
        .timeseries()
        .values.squeeze()
    )
    assert (temperature_global[1:] - temperature_global[:-1] >= 0).all()

@pytest.mark.parametrize("emms_co2_level", [0, 5])
def test_co2_emms_other_rf_run(package, emms_co2_level):
    package.set_zero_config()

    df = zero_emissions.timeseries()

    time = zero_emissions["time"]
    emms_fossil_co2 = np.zeros(len(time))
    emms_fossil_co2[20:] = emms_co2_level

    df.loc[
        (
            df.index.get_level_values("variable")
            == "Emissions|CO2|MAGICC Fossil and Industrial"
        ),
        :,
    ] = emms_fossil_co2

    scen = MAGICCData(df)
    scen.metadata["header"] = "Test CO2 emissions with other RF file"

    forcing_external = 2.0 * np.arange(0, len(time)) / len(time)
    forcing_ext = MAGICCData(
        forcing_external,
        index=time,
        columns={
            "scenario": ["idealised"],
            "model": ["unspecified"],
            "climate_model": ["unspecified"],
            "variable": ["Radiative Forcing|Extra"],
            "unit": ["W / m^2"],
            "todo": ["SET"],
            "region": ["World"],
        },
    )
    forcing_ext.metadata = {"header": "External radiative forcing file for testing"}
    forcing_ext_filename = "CUSTOM_EXTRA_RF.IN"
    forcing_ext.write(join(package.run_dir, forcing_ext_filename), package.version)

    # TODO: fix endyear so it takes from scenario input by default
    results = package.run(
        scen,
        endyear=max(time).year,
        rf_extra_read=1,  # fix writing of 'True'
        file_extra_rf=forcing_ext_filename,
        rf_total_runmodus="ALL",
        rf_initialization_method="ZEROSTARTSHIFT",
        rf_total_constantafteryr=5000,
        co2_switchfromconc2emis_year=min(scen["time"]).year,
        out_emissions=1,
        only=["Radiative Forcing", "Emissions|CO2|MAGICC Fossil and Industrial"],
    )

    np.testing.assert_allclose(
        results.filter(variable="Em*CO2*Fossil*", region="World")
        .timeseries()
        .values.squeeze(),
        emms_fossil_co2,
    )
    # CO2 temperature feedbacks mean that you get a CO2 outgassing, hence CO2 forcing. As a
    # result radiative forcing values don't match exactly. Numerical precision adds to this.
    ext_rf_output_vals = (
        results.filter(variable="Radiative Forcing", region="World")
        .timeseries()
        .values.squeeze()
    )
    zero_rows = (forcing_external == 0) & (ext_rf_output_vals == 0)

    greater_equal_rows = ext_rf_output_vals >= forcing_external
    close_rows_denominator = forcing_external
    close_rows_denominator[zero_rows] = 10 ** -10  # avoid divide by zero
    close_rows = (
        np.abs(ext_rf_output_vals - forcing_external) / close_rows_denominator
        <= 10 ** -3
    )
    matching_rows = greater_equal_rows | close_rows
    assert matching_rows.all()

@patch("pymagicc.core.listdir")
def test_get_output_filenames(mock_listdir):
    mock_listdir.return_value = [
        "DAT_SLR_SEMIEMPI_RATE.OUT",
        "DAT_SLR_SEMIEMPI_RATE.BINOUT",
        "TEMP_OCEANLAYERS.BINOUT",
        "TEMP_OCEANLAYERS.OUT",
        "DAT_SLR_AIS_SMB.OUT",
        "EXTRA.OTHER",
        "PARAMETERS.OUT",
    ]

    m = MAGICC6()
    obs = sorted(m._get_output_filenames())
    exp = sorted(
        [
            "DAT_SLR_SEMIEMPI_RATE.BINOUT",
            "DAT_SLR_AIS_SMB.OUT",
            "TEMP_OCEANLAYERS.OUT",
            "EXTRA.OTHER",
        ]
    )
    np.testing.assert_array_equal(obs, exp)

def test_default_config(package):
    package.default_config
    if package.version == 6:
        expected_config = {
            "file_emissionscenario": "RCP26.SCEN",
            "file_tuningmodel": "PYMAGICC",
        }
    else:
        expected_config = {
            "file_emisscen_2": "NONE",
            "file_emisscen_3": "NONE",
            "file_emisscen_4": "NONE",
            "file_emisscen_5": "NONE",
            "file_emisscen_6": "NONE",
            "file_emisscen_7": "NONE",
            "file_emisscen_8": "NONE",
            "file_tuningmodel_1": "PYMAGICC",
            "file_tuningmodel_2": "USER",
            "file_tuningmodel_3": "USER",
            "file_tuningmodel_4": "USER",
            "file_tuningmodel_5": "USER",
            "file_tuningmodel_6": "USER",
            "file_tuningmodel_7": "USER",
            "file_tuningmodel_8": "USER",
            "file_tuningmodel_9": "USER",
            "file_tuningmodel_10": "USER",
        }

    cfg = read_cfg_file(join(package.run_dir, "MAGCFG_USER.CFG"))
    for key, expected in expected_config.items():
        assert cfg["nml_allcfgs"][key] == expected

def test_out_forcing(package):
    with pytest.warns(UserWarning, match="Not reading file: DAT_VOLCANIC_RF.OUT"):
        run_kwargs = {"out_forcing": True, "out_ascii_binary": "ASCII"}
        if package.version == 7:
            run_kwargs["out_forcing_subannual"] = True

        res = package.run(**run_kwargs)

    # The results should include sub-annual timeseries by default
    volc_monthly = MAGICCData(join(package.out_dir, "DAT_VOLCANIC_RF.OUT"))
    assert (volc_monthly["time"][1] - volc_monthly["time"][0]).days == 30
    volc = res.filter(variable="Radiative Forcing|Volcanic")
    assert (volc["time"][1] - volc["time"][0]).days == 365

    # make sure annual series also read sensibly
    if package.version == 6:
        exp = 289.2
        assert (
            res.filter(variable="*Conc*CO2", year=1876, region="World").values.squeeze()
            == exp
        )

def test_format_config(package):
    inp = {
        "out_temperature": True,
        "out_allowdynamicvars": False,
        "out_keydata1_vars": ["DAT_SURF_TEMP"],
        "out_dynamic_vars": ["DAT_SURF_TEMP"],
        "out_zero_temp_period": [1990, 2000],
    }
    exp = {
        "out_temperature": 1,
        "out_allowdynamicvars": 0,
        "out_keydata1_vars": ["DAT_SURF_TEMP"],
        "out_dynamic_vars": ["DAT_SURF_TEMP"],
        "out_zero_temp_period": [1990, 2000],
    }

    res = package._convert_out_config_flags_to_integers(inp)
    assert exp == res

@pytest.mark.slow
def test_limit_output(package):
    # Check that we can run the model and only output a single variable
    package.set_output_variables(write_binary=True, write_ascii=False)
    if package.version == 6:
        # MAGICC6 doesn't have this capability
        with pytest.raises(CalledProcessError):
            package.run(out_dynamic_vars=["DAT_SURFACE_TEMP"])
        return

    res = package.run(out_dynamic_vars=["DAT_SURFACE_TEMP"])
    assert listdir(package.out_dir) == ["DAT_SURFACE_TEMP.BINOUT"]
    assert res["variable"].unique() == ["Surface Temperature"]

@pytest.mark.slow
def test_stderr_debug(package):
    if package.version == 6:
        # MAGICC6 doesn't have this capability
        with pytest.raises(ValueError):
            res = package.run(debug=True, only=["Surface Temperature"])
        return

    res = package.run(debug=True, only=["Surface Temperature"])

    assert "stderr" in res.metadata
    assert "<DEBUG>" in res.metadata["stderr"]

@pytest.mark.slow
def test_stderr_verbose(package):
    if package.version == 6:
        # MAGICC6 doesn't have this capability
        with pytest.raises(ValueError):
            res = package.run(debug=True, only=["Surface Temperature"])
        return

    res = package.run(debug="verbose", only=["Surface Temperature"])

    assert "stderr" in res.metadata

    assert "<DEBUG>" not in res.metadata["stderr"]
    assert "<INFO>" in res.metadata["stderr"]

@pytest.mark.slow
def test_stderr_accessible_on_failure(package):
    raised = False
    try:
        package.run(invalid_parameter=True, verbose=True)
    except CalledProcessError as e:
        stderr = e.stderr.decode("ascii")
        assert stderr
        raised = True
    finally:
        assert raised

@pytest.mark.slow
@pytest.mark.parametrize(
    "level,raises",
    [("WARNING", True), ("ERROR", True), ("FATAL", True), ("INFO", False)],
)
def test_stderr_warning_raises_warning(mocker, level, raises):

    # Run magicc, but replaces the error message
    def run(*args, **kwargs):
        import subprocess

        r = subprocess.run(*args, **kwargs)
        r.stderr = level.encode("ascii")
        return r

    mock_run = mocker.patch("pymagicc.core.subprocess").run
    mock_run.side_effect = run

    try:
        with MAGICC7() as m:
            if raises:
                with pytest.warns(
                    UserWarning, match=r"magicc logged a {} message*".format(level)
                ):
                    m.run(out_dynamic_vars=["DAT_SURFACE_TEMP"])
            else:
                with pytest.warns(None) as record:
                    m.run(out_dynamic_vars=["DAT_SURFACE_TEMP"])
                assert len(record) == 0
    except FileNotFoundError:
        pytest.skip("MAGICC7 not installed")

@pytest.mark.slow
def test_empty_output(package):
    package.strict = False
    package.set_output_variables()

    if package.version == 7:
        res = package.run(out_dynamic_vars=[""])
    else:
        res = package.run()

    assert len(res) == 0

@pytest.mark.slow
def test_empty_output_strict(package):
    package.strict = True
    package.set_output_variables()
    with pytest.raises(ValueError, match="No output found. Check configuration"):
        if package.version == 7:
            package.run(out_dynamic_vars=[""])
        else:
            package.run()

@pytest.mark.slow
def test_empty_output_with_params(package):
    package.strict = False
    package.set_output_variables(parameters=True)
    if package.version == 7:
        res = package.run(out_dynamic_vars=[""])
    else:
        res = package.run()

    assert len(res) == 0
    assert len(res.metadata["parameters"])

@pytest.mark.slow
def test_failure_message(package, capsys):
    package.strict = False
    emisscen_key = "file_emissionscenario" if package.version == 6 else "file_emisscen"
    remove(join(package.run_dir, package.default_config["nml_allcfgs"][emisscen_key]))
    with pytest.raises(CalledProcessError):
        package.run()

    assert "stderr:\n" in capsys.readouterr().out

max_stars_repo_name: gabriel-abrahao/pymagicc

max_issues_repo_name: gabriel-abrahao/pymagicc

max_forks_repo_name: gabriel-abrahao/pymagicc


text:
"""Test the measurement package."""

from typing import Dict

import numpy as np
import pytest
import xarray as xr

from weldx.constants import Q_, U_
from weldx.core import MathematicalExpression
from weldx.measurement import (
    Error,
    MeasurementChain,
    MeasurementEquipment,
    Signal,
    SignalSource,
    SignalTransformation,
)

from ._helpers import get_test_name

class TestMeasurementChain:
    """Test the `MeasurementChain` class."""

    # helper functions -----------------------------------------------------------------

    @staticmethod
    def _default_source_kwargs(kwargs: Dict = None) -> Dict:
        """Update a dict with default keyword arguments to create a `SignalSource`."""
        default_kwargs = dict(
            name="source", output_signal=Signal("analog", "V"), error=Error(0.01)
        )

        if kwargs is not None:
            default_kwargs.update(kwargs)

        return default_kwargs

    @classmethod
    def _default_init_kwargs(
        cls, kwargs: Dict = None, source_kwargs: Dict = None
    ) -> Dict:
        """Return a dictionary of keyword arguments required by the `__init__` method.

        Parameters
        ----------
        kwargs :
            A dictionary containing some key word arguments that should replace the
            default ones.
        source_kwargs :
            A dictionary of key word arguments that should replace the arguments used
            for the default source.

        Returns
        -------
        Dict :
            Dictionary with keyword arguments for the `__init__` method

        """
        source_kwargs = cls._default_source_kwargs(source_kwargs)

        default_kwargs = dict(
            name="name",
            source=SignalSource(**source_kwargs),
            signal_data=[1, 3, 5],
        )
        if kwargs is not None:
            default_kwargs.update(kwargs)

        return default_kwargs

    @staticmethod
    def _default_transformation(kwargs: Dict = None) -> SignalTransformation:
        """Return a default `SignalTransformation`.

        Use the kwargs parameter to modify the default values.

        """
        default_kwargs = dict(
            name="transformation",
            error=Error(0.1),
            func=MathematicalExpression("a*x", parameters={"a": Q_(1, "1/V")}),
            type_transformation="AD",
        )
        if kwargs is not None:
            default_kwargs.update(kwargs)

        return SignalTransformation(**default_kwargs)

    @classmethod
    def _default_add_transformation_kwargs(cls, kwargs: Dict = None) -> Dict:
        """Update a dict with default keyword arguments to call `add_transformation`."""
        default_kwargs = dict(
            transformation=cls._default_transformation(),
            error=Error(0.02),
            output_signal_type="digital",
            output_signal_unit="",
        )
        if kwargs is not None:
            default_kwargs.update(kwargs)

        return default_kwargs

    # test_init ------------------------------------------------------------------------

    @staticmethod
    @pytest.mark.parametrize(
        "kwargs, source_kwargs",
        [
            ({}, {}),
            (dict(signal_data=None), dict(output_signal=Signal("analog", "V", [1]))),
        ],
    )
    def test_init(kwargs: Dict, source_kwargs: Dict):
        """Test the `__init__` method of the `MeasurementChain`.

        Parameters
        ----------
        kwargs:
            A dictionary with keyword arguments that are passed to the `__init__`
            method. Missing arguments are added.
        source_kwargs :
            A dictionary with keyword arguments that are used to construct the
            `SignalSource` that is passed to the `__init__` method. Missing arguments
            are added.

        """
        kwargs = TestMeasurementChain._default_init_kwargs(kwargs, source_kwargs)
        MeasurementChain(**kwargs)

    # test_init_exceptions -------------------------------------------------------------

    @staticmethod
    @pytest.mark.parametrize(
        "kwargs, source_kwargs,  exception_type, test_name",
        [({}, {"output_signal": Signal("analog", "V", [1])}, KeyError, "# 2x data")],
        ids=get_test_name,
    )
    def test_init_exceptions(
        kwargs: Dict, source_kwargs: Dict, exception_type, test_name: str
    ):
        """Test the exceptions of the `__init__` method.

        Parameters
        ----------
        kwargs :
            A dictionary with keyword arguments that are passed to the `__init__`
            method. Missing arguments are added.
        source_kwargs :
            A dictionary with keyword arguments that are used to construct the
            `SignalSource` that is passed to the `__init__` method. Missing arguments
            are added.
        exception_type :
            The expected exception type
        test_name :
            Name of the test

        """
        kwargs = TestMeasurementChain._default_init_kwargs(kwargs, source_kwargs)
        with pytest.raises(exception_type):
            MeasurementChain(**kwargs)

    # test_from_equipment --------------------------------------------------------------

    @pytest.mark.parametrize(
        "num_sources, source_name, exception",
        [
            (0, None, ValueError),
            (1, None, None),
            (2, "source_1", None),
            (2, "wrong name", KeyError),
            (2, None, ValueError),
        ],
    )
    def test_from_equipment(
        self, num_sources: int, source_name: str, exception: Exception
    ):
        """Test the `from_equipment` factory and its exceptions.

        Parameters
        ----------
        num_sources :
            Number of sources of the generated equipment
        source_name :
            Corresponding parameter of `from_equipment`
        exception :
            Expected exception

        """
        sources = [
            SignalSource(**self._default_source_kwargs({"name": f"source_{i}"}))
            for i in range(num_sources)
        ]
        equipment = MeasurementEquipment("Equipment", sources=sources)

        if exception is not None:
            with pytest.raises(exception):
                MeasurementChain.from_equipment(
                    name="name", equipment=equipment, source_name=source_name
                )
        else:
            MeasurementChain.from_equipment(
                name="name", equipment=equipment, source_name=source_name
            )

    # test_add_transformations ---------------------------------------------------------

    @pytest.mark.parametrize(
        "tf_kwargs, exp_signal_type, exp_signal_unit",
        [
            ({}, "digital", U_("")),
            (dict(type_transformation="AA"), "analog", U_("")),
            (dict(type_transformation=None), "analog", U_("")),
            (dict(func=None), "digital", U_("V")),
        ],
    )
    def test_add_transformation(self, tf_kwargs, exp_signal_type, exp_signal_unit):
        """Test the `add_transformation` method of the `MeasurementChain`.

        Parameters
        ----------
        tf_kwargs:
            A dictionary with keyword arguments that are used to construct the
            `SignalTransformation` that is passed to the `add_transformation` method.
            Missing arguments are added.
        exp_signal_type :
            The expected signal type after the transformation
        exp_signal_unit :
            The expected unit after the transformation

        """
        mc = MeasurementChain(**self._default_init_kwargs())

        mc.add_transformation(self._default_transformation(tf_kwargs))

        signal = mc.output_signal
        assert signal.signal_type == exp_signal_type
        assert U_(signal.units) == exp_signal_unit

    # test_add_transformation_exceptions -----------------------------------------------

    @pytest.mark.parametrize(
        "tf_kwargs, input_signal_source, exception_type, test_name",
        [
            (dict(type_transformation="DA"), None, ValueError, "# inv. signal type #1"),
            (dict(type_transformation="DD"), None, ValueError, "# inv. signal type #2"),
            ({}, "not found", KeyError, "# invalid input signal source"),
            (dict(name="source"), None, KeyError, "# Name does already exist"),
            (
                dict(func=MathematicalExpression("x+a", parameters={"a": Q_(1, "A")})),
                None,
                ValueError,
                "# incompatible function",
            ),
        ],
        ids=get_test_name,
    )
    def test_add_transformation_exceptions(
        self, tf_kwargs: Dict, input_signal_source: str, exception_type, test_name: str
    ):
        """Test the exceptions of the `add_transformation` method.

        Parameters
        ----------
        tf_kwargs:
            A dictionary with keyword arguments that are used to construct the
            `SignalTransformation` that is passed to the `add_transformation` method.
            Missing arguments are added.
        input_signal_source :
            The value of the corresponding parameter of 'add_transformation'
        exception_type :
            The expected exception type
        test_name :
            Name of the test

        """
        mc = MeasurementChain(**self._default_init_kwargs())

        tf = self._default_transformation(tf_kwargs)

        with pytest.raises(exception_type):
            mc.add_transformation(tf, input_signal_source=input_signal_source)

    # test_add_transformation_from_equipment -------------------------------------------

    @pytest.mark.parametrize(
        "num_transformations, transformation_name, exception",
        [
            (0, None, ValueError),
            (1, None, None),
            (2, "transformation_1", None),
            (2, "wrong name", KeyError),
            (2, None, ValueError),
        ],
    )
    def test_add_transformation_from_equipment(
        self, num_transformations: int, transformation_name: str, exception
    ):
        """Test `add_transformation_from_equipment` and its exceptions.

        Parameters
        ----------
        num_transformations :
            Number of transformations of the generated equipment
        transformation_name :
            Corresponding parameter of `add_transformation_from_equipment`
        exception :
            Expected exception

        """
        mc = MeasurementChain(**self._default_init_kwargs())
        transformations = [
            self._default_transformation({"name": f"transformation_{i}"})
            for i in range(num_transformations)
        ]
        equipment = MeasurementEquipment(name="name", transformations=transformations)

        if exception is not None:
            with pytest.raises(exception):
                mc.add_transformation_from_equipment(
                    equipment=equipment, transformation_name=transformation_name
                )
        else:
            mc.add_transformation_from_equipment(
                equipment=equipment, transformation_name=transformation_name
            )

    # test_add_signal_data -------------------------------------------------------------

    @pytest.mark.parametrize(
        "kwargs",
        [
            dict(data=xr.DataArray([2, 3])),
            dict(signal_source="source"),
        ],
    )
    def test_add_signal_data(self, kwargs):
        """Test the `add_signal_data` method of the `MeasurementChain`.

        Parameters
        ----------
        kwargs:
            A dictionary with keyword arguments that are passed to the
            `add_signal_data` method. If no name is in the kwargs, a default one is
            added.

        """
        mc = MeasurementChain(**self._default_init_kwargs({"signal_data": None}))
        mc.add_transformation(self._default_transformation())

        full_kwargs = dict(data=xr.DataArray([1, 2]))
        full_kwargs.update(kwargs)

        mc.add_signal_data(**full_kwargs)

    # test_add_signal_data_exceptions --------------------------------------------------

    @pytest.mark.parametrize(
        "kwargs,  exception_type, test_name",
        [
            (dict(signal_source="what"), KeyError, "# invalid signal source"),
            (dict(signal_source="source"), KeyError, "# already has data #1"),
            (dict(signal_source="transformation"), KeyError, "# already has data #2"),
        ],
        ids=get_test_name,
    )
    def test_add_signal_data_exceptions(
        self, kwargs: Dict, exception_type, test_name: str
    ):
        """Test the exceptions of the `add_signal_data` method.

        Parameters
        ----------
        kwargs :
            A dictionary with keyword arguments that are passed to the `add_signal_data`
            method. Missing arguments are added.
        exception_type :
            The expected exception type
        test_name :
            Name of the test

        """
        mc = MeasurementChain(**self._default_init_kwargs())
        mc.add_transformation(self._default_transformation(), data=[1, 2, 3])
        mc.add_transformation(
            self._default_transformation(
                dict(name="transformation 2", type_transformation="DA")
            )
        )

        full_kwargs = dict(data=xr.DataArray([1, 2]))
        full_kwargs.update(kwargs)

        with pytest.raises(exception_type):
            mc.add_signal_data(**full_kwargs)

    # test_get_equipment ---------------------------------------------------------------

    @pytest.mark.parametrize(
        "signal_source, exception",
        [
            ("source", None),
            ("transformation_1", None),
            ("transformation_2", None),
            ("transformation_3", KeyError),
        ],
    )
    def test_get_equipment(self, signal_source, exception):
        """Test the `get_equipment` function and their exceptions.

        Parameters
        ----------
        signal_source :
            Corresponding function parameter
        exception :
            Expected exception

        """
        src_eq = MeasurementEquipment(
            "Source Eq", sources=[SignalSource(**self._default_source_kwargs())]
        )
        tf_eq = MeasurementEquipment(
            "Transformation_eq",
            transformations=[
                self._default_transformation({"name": "transformation_1"})
            ],
        )

        mc = MeasurementChain.from_equipment("Chain", src_eq)
        mc.add_transformation_from_equipment(tf_eq)
        mc.create_transformation("transformation_2", None, output_signal_unit="A")

        if exception is not None:
            with pytest.raises(exception):
                mc.get_equipment(signal_source=signal_source)
        else:
            mc.get_equipment(signal_source=signal_source)

    # test_get_signal_data -------------------------------------------------------------

    def test_get_signal_data(self):
        """Test the `get_signal_data` method.

        This test assures that the returned data is identical to the one passed
        to the
        measurement chain and that a key error is raised if the requested data is
        not
        present.

        """
        data = xr.DataArray([1, 2, 3])

        mc = MeasurementChain(**self._default_init_kwargs())
        mc.add_transformation(self._default_transformation(), data=data)
        mc.create_transformation("transformation_2", None, output_signal_unit="A")

        assert np.all(mc.get_signal_data("transformation") == data)

        # no data
        with pytest.raises(KeyError):
            mc.get_signal_data("transformation_2")

        # source not present
        with pytest.raises(KeyError):
            mc.get_signal_data("not found")

    # test_get_transformation ----------------------------------------------------------

    def test_get_transformation(self):
        """Test the `get_transformation` method."""
        mc = MeasurementChain(**self._default_init_kwargs())
        mc.add_transformation(self._default_transformation())

        transformation = mc.get_transformation("transformation")

        assert transformation == self._default_transformation()

    # test_get_transformation_exception ------------------------------------------------

    def test_get_transformation_exception(self):
        """Test that a `KeyError` is raised if the transformation does not exist."""
        mc = MeasurementChain(**self._default_init_kwargs())
        mc.add_transformation(self._default_transformation())

        with pytest.raises(KeyError):
            mc.get_transformation("not found")

ext: py

max_stars_repo_name: CagtayFabry/weldx

max_issues_repo_name: BAMWelDX/weldx

max_forks_repo_name: CagtayFabry/weldx

reason: import numpy


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-

from __future__ import division

import numpy as np
import matplotlib.pyplot as plt
from data import DataTrain
import os

class Shape(object):
        def __init__(self, landmarks, name):
            self.xy = np.array(landmarks, dtype=np.int)
            self.nr_landmarks =  self.xy.shape[0]/2
            self.name = name
            self.centroid = self.get_centroid()

        def load_shape(self, shape):
            pass

        def translate_shape(self, direction):

            T = np.identity(3)
            T[:2, 2] = direction[:2]

            return T

        def get_centroid(self):
            centroid = self.xy.mean(0)

            return centroid

        def align_shape(self, mean_shape):
            """
            Procrustes fit
            """
            n, m = mean_shape.xy.shape
            ny, my = self.xy.shape
            reflection = 'best'
            scaling = True

            X0 = mean_shape.xy - mean_shape.centroid
            Y0 = self.xy - self.centroid

            ssX = (X0**2.).sum()
            ssY = (Y0**2.).sum()

            # centred Frobenius norm
            normX = np.sqrt(ssX)
            normY = np.sqrt(ssY)

            # scale to equal (unit) norm
            X0 /= normX
            Y0 /= normY

            if my < m:
                Y0 = np.concatenate((Y0, np.zeros(n, m-my)),0)

            # optimum rotation matrix of Y
            A = np.dot(X0.T, Y0)
            U,s,Vt = np.linalg.svd(A,full_matrices=False)
            V = Vt.T
            T = np.dot(V, U.T)

            if reflection is not 'best':

                # does the current solution use a reflection?
                have_reflection = np.linalg.det(T) < 0

                # if that's not what was specified, force another reflection
                if reflection != have_reflection:
                    V[:,-1] *= -1
                    s[-1] *= -1
                    T = np.dot(V, U.T)

            traceTA = s.sum()

            if scaling:

                # optimum scaling of Y
                b = traceTA * normX / normY

                # standarised distance between X and b*Y*T + c
                d = 1 - traceTA**2

                # transformed coords
                Z = normX*traceTA*np.dot(Y0, T) + mean_shape.centroid

            else:
                b = 1
                d = 1 + ssY/ssX - 2 * traceTA * normY / normX
                Z = normY*np.dot(Y0, T) + mean_shape.centroid

            # transformation matrix
            if my < m:
                T = T[:my,:]
            c = mean_shape.centroid - b*np.dot(self.centroid, T)

            #transformation values
            tform = {'rotation':T, 'scale':b, 'translation':c}

            self.xy = Z.copy()
            #print d

        def plot_shape(self):
            return plt.plot(self.xy[:,0], self.xy[:,1], 'o')

def main(filename_xy, path_images):
    my_shapes = []
    my_data = DataTrain(filename_xy, path_images, 200, ids=False)

    for row in my_data.data.iterrows():
        my_shapes.append(Shape(row[1][:].reshape(45, 2),'lero'))
    for shape in my_shapes:
        shape.align_shape(my_shapes[0])

    for shape in my_shapes:
        shape.plot_shape()
    plt.axis('equal')
    plt.show()

if __name__ == '__main__':
    #TODO this should be a test
    package_directory = os.path.dirname(os.path.abspath(__file__))
    main(os.path.join(package_directory, "data/use.txt"), os.path.join(package_directory,"images/"))

max_stars_repo_name: celiacintas/aampy

max_issues_repo_name: celiacintas/aampy

max_forks_repo_name: celiacintas/aampy


text:
#%%
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import community
import scipy.ndimage as nd
import collections
from scipy.stats import entropy

#import warnings; warnings.simplefilter('ignore')

#%%
KERNEL = np.zeros(9, dtype=np.uint8).reshape((3,3))
#print(KERNEL)

#%%
# use of convolution, courtesy of salt-die
class generator_state:
    seed = nx.to_numpy_matrix(nx.erdos_renyi_graph(33, 0.01, seed = 1, directed=True))

    def next_state(self):
        seed = self.seed
        seed = np.pad(seed, (1, 1), 'constant')
        neighbor_count = nd.convolve(seed, KERNEL, mode="constant")
        self.seed = np.where(((seed == 0) & (neighbor_count > 0) & (neighbor_count <= 4)) |
                                 ((seed == 1) & (neighbor_count > 4)), 1, 0)
        return(self.seed)

    def run(self):
        while True:
            for _ in range(100):
                self.next_state()
            return(self.seed)

#%%
for i in range(0, 512, 1):
    print(i)

    kernel_seed = f'{i:09b}'
    KERNEL = np.array(list(kernel_seed), dtype=np.uint8).reshape((3,3))

    gen = generator_state()
    gen_g1 = nx.from_numpy_matrix(gen.run())

    degree_sequence = sorted([d for n, d in gen_g1.degree()], reverse=True)
    degreeCount = collections.Counter(degree_sequence)

    degreeCount_init = collections.Counter(degree_sequence)
    deg_init, cnt_init = zip(*degreeCount_init.items())
    cnt_avg_init_0 = cnt_init[-1]/sum(cnt_init)

    del degreeCount[0] # exclude graphs with zero edges

    if deg_init[-1] == 0 and cnt_avg_init_0 > 0.10: # exclude graphs with edge-less nodes over 10%
        ent = 100
    else:
        if len(degreeCount) > 0:
            deg, cnt = zip(*degreeCount.items())
            cnt_avg = [cnt[i]/sum(cnt) for i in range(len(cnt))]
            ent = entropy(cnt_avg)
        else:
            ent = 100

#    print(ent)

    if ent < 2:
        print('hit')
        fig, ax = plt.subplots()
        plt.bar(deg_init, cnt_init, width=0.80, color='b')
        plt.title("Degree Histogram")
        plt.ylabel("Count")
        plt.xlabel("Degree")
        ax.set_xticks([d + 0.4 for d in deg_init])
        ax.set_xticklabels(deg_init)
        # draw graph in inset
        plt.axes([0.4, 0.4, 0.5, 0.5])
        Gcc = gen_g1.subgraph(sorted(nx.connected_components(gen_g1), key=len, reverse=True)[0])
        pos = nx.spring_layout(gen_g1)
        plt.axis('off')
        nx.draw_networkx_nodes(gen_g1, pos, node_size=20)
        nx.draw_networkx_edges(gen_g1, pos, alpha=0.4)
        plt.savefig('%s_degree_histogram.png' % i)

max_stars_repo_name: kazuyamagiwa/graphnet-automata

max_issues_repo_name: kazuyamagiwa/graphnet-automata

max_forks_repo_name: kazuyamagiwa/graphnet-automata

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


text:
import torch
from torch.utils.data.sampler import Sampler
from collections import defaultdict
import numpy as np

def safe_random_choice(input_list, size):
    replace = len(input_list) < size
    return np.random.choice(input_list, size=size, replace=replace)

class MPerClassSampler(Sampler):
    """
    Give m samples per class.

    Args:
        labels (np.ndarray):
            Ground truth of datasets
        m (int):
            M samples per class
        batch_size (int):
            Batch size must be an interger multiple of m
    """
    def __init__(
        self,
        labels,
        m,
        batch_size,
    ):
        # raise NotImplementedError()
        self.labels = labels
        self.m = m
        self.batch_size = batch_size

        self.init_indices()

    def init_indices(self):
        # construct indices
        self.labels_to_indices = defaultdict(list)
        for i, label in enumerate(self.labels):
            self.labels_to_indices[label].append(i)
        for k, v in self.labels_to_indices.items():
            self.labels_to_indices[k] = np.array(v).astype('int')
        self.labels_set = list(self.labels_to_indices.keys())
        self.list_size = self.m * len(self.labels_set)

        # assert
        assert self.list_size >= self.batch_size
        assert self.batch_size % self.m == 0
        self.list_size -= self.list_size % self.batch_size
        self.num_iters = self.list_size // self.batch_size

    def __len__(self):
        return self.list_size

    def __iter__(self):
        idx_list = [0] * self.list_size
        i = 0
        for _ in range(self.num_iters):
            np.random.shuffle(self.labels_set)
            curr_labels_set = self.labels_set[:self.batch_size // self.m]
            for label in curr_labels_set:
                curr_indices = self.labels_to_indices[label]
                idx_list[i:i+self.m] = safe_random_choice(curr_indices, size=self.m)
                i += self.m
        return iter(idx_list)


text:
"""Data loading utilities."""

import glob
import os

import numpy as np
import portalocker
import torch

from utils.vocabulary import OpenAIVocab, Vocab

class LMOrderedIterator:
    def __init__(self, data, bsz, bptt, device='cpu', ext_len=None):
        """
            data -- LongTensor -- the LongTensor is strictly ordered
        """
        self.bsz = bsz
        self.bptt = bptt
        self.ext_len = ext_len if ext_len is not None else 0

        self.device = device

        # Work out how cleanly we can divide the dataset into bsz parts.
        self.n_step = data.size(0) // bsz

        # Trim off any extra elements that wouldn't cleanly fit (remainders).
        data = data.narrow(0, 0, self.n_step * bsz)

        # Evenly divide the data across the bsz batches.
        self.data = data.view(bsz, -1).t().contiguous().to(device)

        # Number of mini-batches
        self.n_batch = (self.n_step + self.bptt - 1) // self.bptt

    def get_batch(self, i, bptt=None):
        if bptt is None: bptt = self.bptt
        seq_len = min(bptt, self.data.size(0) - 1 - i)

        end_idx = i + seq_len
        beg_idx = max(0, i - self.ext_len)

        data = self.data[beg_idx:end_idx]
        target = self.data[i+1:i+1+seq_len]

        return data, target, seq_len

    def get_fixlen_iter(self, start=0):
        for i in range(start, self.data.size(0) - 1, self.bptt):
            yield self.get_batch(i)

    def get_varlen_iter(self, start=0, std=5, min_len=5, max_deviation=3):
        max_len = self.bptt + max_deviation * std
        i = start
        while True:
            bptt = self.bptt if np.random.random() < 0.95 else self.bptt / 2.
            bptt = min(max_len, max(min_len, int(np.random.normal(bptt, std))))
            data, target, seq_len = self.get_batch(i, bptt)
            i += seq_len
            yield data, target, seq_len
            if i >= self.data.size(0) - 2:
                break

    def __iter__(self):
        """Wrapper for get_fixlen_iter."""
        return self.get_fixlen_iter()

class LMShuffledIterator:
    def __init__(self, data, bsz, bptt, device='cpu', ext_len=None, shuffle=False):
        """
            data -- list[LongTensor] -- there is no order among the LongTensors
        """
        self.data = data

        self.bsz = bsz
        self.bptt = bptt
        self.ext_len = ext_len if ext_len is not None else 0

        self.device = device
        self.shuffle = shuffle

    def get_sent_stream(self):
        # index iterator
        epoch_indices = np.random.permutation(len(self.data)) if self.shuffle \
            else np.array(range(len(self.data)))

        # sentence iterator
        for idx in epoch_indices:
            yield self.data[idx]

    def stream_iterator(self, sent_stream):
        # streams for each data in the batch
        streams = [None] * self.bsz

        data = torch.LongTensor(self.bptt, self.bsz)
        target = torch.LongTensor(self.bptt, self.bsz)

        n_retain = 0

        while True:
            # data   : [n_retain+bptt x bsz]
            # target : [bptt x bsz]
            data[n_retain:].fill_(-1)
            target.fill_(-1)

            valid_batch = True

            for i in range(self.bsz):
                n_filled = 0
                try:
                    while n_filled < self.bptt:
                        if streams[i] is None or len(streams[i]) <= 1:
                            streams[i] = next(sent_stream)
                        # number of new tokens to fill in
                        n_new = min(len(streams[i]) - 1, self.bptt - n_filled)
                        # first n_retain tokens are retained from last batch
                        data[n_retain+n_filled:n_retain+n_filled+n_new, i] = \
                            streams[i][:n_new]
                        target[n_filled:n_filled+n_new, i] = \
                            streams[i][1:n_new+1]
                        streams[i] = streams[i][n_new:]
                        n_filled += n_new
                except StopIteration:
                    valid_batch = False
                    break

            if not valid_batch:
                return

            data = data.to(self.device)
            target = target.to(self.device)

            yield data, target, self.bptt

            n_retain = min(data.size(0), self.ext_len)
            if n_retain > 0:
                data[:n_retain] = data[-n_retain:]
            data.resize_(n_retain + self.bptt, data.size(1))

    def __iter__(self):
        # sent_stream is an iterator
        sent_stream = self.get_sent_stream()

        for batch in self.stream_iterator(sent_stream):
            yield batch

class LMMultiFileIterator(LMShuffledIterator):
    def __init__(self, paths, vocab, bsz, bptt, device='cpu', ext_len=None,
        shuffle=False):

        self.paths = paths
        self.vocab = vocab

        self.bsz = bsz
        self.bptt = bptt
        self.ext_len = ext_len if ext_len is not None else 0

        self.device = device
        self.shuffle = shuffle

    def get_sent_stream(self, path):
        sents = self.vocab.encode_file(path, add_double_eos=True)
        if self.shuffle:
            np.random.shuffle(sents)
        # Create virtual sentences for wikipedia data.
        # TODO: Load a few files at a time, then chunk?
        if type(sents) == torch.Tensor:
            return iter(sents.split(len(sents) // self.bsz))
        return iter(sents)

    def __iter__(self):
        if self.shuffle:
            np.random.shuffle(self.paths)

        for path in self.paths:
            # sent_stream is an iterator
            sent_stream = self.get_sent_stream(path)
            for batch in self.stream_iterator(sent_stream):
                yield batch

class Corpus:
    def __init__(self, path, dataset, use_bpe, *args, **kwargs):
        self.dataset = dataset
        if use_bpe:
            self.vocab = OpenAIVocab(kwargs['max_size'], kwargs.get('vocab_file'))
        else:
            self.vocab = Vocab(*args, **kwargs)

        if self.dataset in ['ptb', 'wt2', 'enwik8', 'text8']:
            self.vocab.count_file(os.path.join(path, 'train.txt'))
            self.vocab.count_file(os.path.join(path, 'valid.txt'))
            self.vocab.count_file(os.path.join(path, 'test.txt'))
        elif self.dataset == 'wt103' or self.dataset == 'wt2':
            self.vocab.count_file(os.path.join(path, 'train.txt'))
        elif self.dataset == 'wt103-normal':
            self.vocab.count_file(os.path.join(path, 'wiki.train.tokens'))
        elif self.dataset == 'lm1b':
            train_path_pattern = os.path.join(
                path, '1-billion-word-language-modeling-benchmark-r13output',
                'training-monolingual.tokenized.shuffled', 'news.en-*')
            train_paths = glob.glob(train_path_pattern)
        elif self.dataset == 'wiki':
            file_path_pattern = os.path.join(path, '*/wiki_*.txt')
            file_paths = glob.glob(file_path_pattern)
            assert file_paths, f'Nothing found at {file_path_pattern}'

        # the vocab will load from file when build_vocab() is called
        self.vocab.build_vocab()

        if self.dataset in ['ptb', 'wt2', 'wt103']:
            self.train = self.vocab.encode_file(
                os.path.join(path, 'train.txt'), ordered=True)
            self.valid = self.vocab.encode_file(
                os.path.join(path, 'valid.txt'), ordered=True)
            self.test = self.vocab.encode_file(
                os.path.join(path, 'test.txt'), ordered=True)
        elif self.dataset in ['enwik8', 'text8']:
            self.train = self.vocab.encode_file(
                os.path.join(path, 'train.txt'), ordered=True, add_eos=False)
            self.valid = self.vocab.encode_file(
                os.path.join(path, 'valid.txt'), ordered=True, add_eos=False)
            self.test = self.vocab.encode_file(
                os.path.join(path, 'test.txt'), ordered=True, add_eos=False)
        elif self.dataset == 'lm1b':
            self.train = train_paths
            self.valid = self.vocab.encode_file(
                os.path.join(path, 'valid.txt'), ordered=False, add_double_eos=True)
            self.test = self.vocab.encode_file(
                os.path.join(path, 'test.txt'), ordered=False, add_double_eos=True)
        elif self.dataset == 'wiki':
            valid_path = sorted(file_paths)[42]
            test_path = sorted(file_paths)[1337]
            self.valid = self.vocab.encode_file(valid_path, ordered=True)
            self.test = self.vocab.encode_file(test_path, ordered=True)
            self.train = list(set(file_paths) - set((valid_path, test_path)))
        elif self.dataset in ['wt103-normal']:
            self.train = self.vocab.encode_file(
                os.path.join(path, 'wiki.train.tokens'), ordered=True, add_eos=False)
            self.valid = self.vocab.encode_file(
                os.path.join(path, 'wiki.valid.tokens'), ordered=True, add_eos=False)
            self.test = self.vocab.encode_file(
                os.path.join(path, 'wiki.test.tokens'), ordered=True, add_eos=False)

    def get_dist_iterator(self, split: str, *args, rank: int = 0, max_rank: int = 1, **kwargs):
        """Get an iterator that only operates on rank//max_rank independent subset of the data."""
        data = self.__getattribute__(split)
        subset = list(chunk(data, max_rank))[rank]
        if self.dataset in ['lm1b', 'wiki']:
            if split == 'train':
                return LMMultiFileIterator(subset, self.vocab, *args, **kwargs)

        return LMOrderedIterator(subset, *args, **kwargs)

    def get_iterator(self, split: str, *args, **kwargs):
        """Get an iterator over the corpus.

        Each next() returns (data, target, seq_length).
        data and target have shape (bptt, bsz) and seq_length is a scalar.
        """
        data = self.__getattribute__(split)
        if self.dataset in ['ptb', 'wt2', 'wt103', 'enwik8', 'text8', 'wt103-normal']:
            return LMOrderedIterator(data, *args, **kwargs)
        if self.dataset == 'lm1b':
            if split in ['valid', 'test']:
                return LMShuffledIterator(data, *args, **kwargs)

            kwargs['shuffle'] = True
            return LMMultiFileIterator(data, self.vocab, *args, **kwargs)
        if self.dataset == 'wiki':
            if split == 'train':
                return LMMultiFileIterator(data, self.vocab, *args, **kwargs)
            return LMOrderedIterator(data, *args, **kwargs)

def get_lm_corpus(datadir: str, dataset: str, use_bpe=False, max_size=None) -> Corpus:
    """Factory method for Corpus.

    Arguments:
        datadir: Where does the data live?
        dataset: eg 'wt103' which tells the Corpus how to parse the data.
    """
    cache_filepath = os.path.join(datadir, 'cache.pt.bpe' if use_bpe else 'cache.pt')
    # Don't cache dataset for wiki, it's just a file list.
    if os.path.exists(cache_filepath) and dataset != 'wiki':
        print('Loading cached dataset...')
        corpus = torch.load(cache_filepath)
    else:
        print('Producing dataset {}...'.format(dataset))
        kwargs = {'max_size': max_size}
        if dataset in ['wt103', 'wt2', 'wt103-normal']:
            kwargs['special'] = ['<eos>']
            kwargs['lower_case'] = False
        elif dataset == 'ptb':
            kwargs['special'] = ['<eos>']
            kwargs['lower_case'] = True
        elif dataset == 'lm1b':
            kwargs['special'] = []
            kwargs['lower_case'] = False
            kwargs['vocab_file'] = os.path.join(datadir, '1b_word_vocab.txt')
        elif dataset in ['enwik8', 'text8']:
            pass

        corpus = Corpus(datadir, dataset, use_bpe, **kwargs)
        with portalocker.Lock(cache_filepath, timeout=60) as _:
            torch.save(corpus, cache_filepath)

    return corpus

def chunk(a: list, n: int):
    """Split `a` into `n` chunks, with the last bucket taking the remaining.

    https://stackoverflow.com/a/2135920
    """
    k, m = divmod(len(a), n)
    return (a[i * k + min(i, m):(i + 1) * k + min(i + 1, m)] for i in range(n))

def main():
    import argparse
    parser = argparse.ArgumentParser(description='unit test')
    parser.add_argument('--datadir', type=str, default='../data/text8',
                        help='location of the data corpus')
    parser.add_argument('--dataset', type=str, default='text8',
                        choices=['ptb', 'wt2', 'wt103', 'lm1b', 'enwik8', 'text8', 'wt103-normal', 'wiki'],
                        help='dataset name')
    parser.add_argument('--debug', action='store_true')
    args = parser.parse_args()

    corpus = get_lm_corpus(args.datadir, args.dataset, use_bpe=True)
    print(f'Vocab size : {len(corpus.vocab)}')
    #tr_iter = corpus.get_iterator('train', 16, 150, 'cpu', ext_len=0)
    # Loop through all the data to force caching.
    for split in ('train', 'valid', 'test'):
        tr_iter = corpus.get_dist_iterator(split, rank=0, max_rank=1, bsz=16, bptt=150, device='cpu', ext_len=0)
        for data, target, seq_len in tr_iter:
            if args.debug:
                print(data.shape, target.shape, seq_len, len(list(tr_iter)))
                break

if __name__ == '__main__':
    main()

max_stars_repo_name: JasonBenn/duet

max_issues_repo_name: JasonBenn/duet

max_forks_repo_name: JasonBenn/duet


text:
import pandas as pd

import numpy as np

from sklearn.ensemble import GradientBoostingRegressor,BaggingRegressor,RandomForestRegressor
from sklearn.model_selection import train_test_split

from sklearn import metrics

# plot
import matplotlib.pyplot as plt

pd.set_option('display.max_columns', None)
DATA_PATH ="C:/Users/user/Desktop/Data/Data.xlsx"
full_dataframe = pd.read_excel(DATA_PATH, sep=',')

X = full_dataframe[["game_positive","game_negative","game_owners","game_initialprice","game_discount"]]
y = full_dataframe["game_price"]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=43)

gbrt = GradientBoostingRegressor()
# GradientBoostingRegressor 선언
gbrt.fit(X_train,y_train)
# 나의 데이터를 GradientBoostingRegressor에 적용시킨다.
gbrt_y_pred = gbrt.predict(X_test)
# GradientBoostingRegressor 알고리즘을 통해 Y를 예측한다.

# print("Feature importances")
# print(gbrt.feature_importances_)
# GradientBoostingRegressor 메소드중 feature 중요도를 나타내는 것이 있었다.

# print('Gradientboostingregressor Mean Absolute Error:', metrics.mean_absolute_error(y,gbrt_y_pred))
# print('Gradientboostingregressor Mean Squared Error:', metrics.mean_squared_error(y,gbrt_y_pred))
# print('Gradientboostingregressor Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y,gbrt_y_pred)))
# print('Gradientboostingregressor Accuracy:', metrics.r2_score(y,gbrt_y_pred))
# print("R-squared for Train: %.2f" %gbrt.score(X,y))
# print("R-squared for Test: %.2f" %gbrt.score(X,y))

# 예측 전의 Y 값과 예측 후의 Y값
gbrt_df = pd.DataFrame({'Actual':y_test, 'Predicted':gbrt_y_pred})

# Gradientboosting regressor 알고리즘을 사용해서 가격예측을 해야한다.
# 예측되어진 값에서 예측전의 값을 빼면 두 값의 편차가 생긴다.
# 이 편차값을 예측되어진 값에서 뺀다면 원래 초기가격에 맞춰지지 않을까 생각을 하였다.
# 이 때 두가지 가정이 주어진다
# 적용 전의 값 - 적용 후의 값 : 양수 or 음수
# 양수일 때는 적용전의 값이 클때이며
# 음수일 때는 적용 후의 값이 더 클때이다.
# 양수일때와 음수일때의 계산을 다르게 해야할 필요를 느꼈다.

# 예측전의 값과 예측후의 값의 차이를 양수와 음수로 비교해서 , 실제 가격에 맞춘다.
# 예측후의 값에서 예측 전의 값을 뺏을 때 양수일 경우
gbrt_df.loc[gbrt_df['Predicted']-gbrt_df['Actual']>0,"Final_Predict"]=gbrt_df['Predicted']-(gbrt_df['Predicted']-gbrt_df['Actual'])
# 예측후의 값에서 예측 전의 값을 뺏을 때 음수일 경우
gbrt_df.loc[gbrt_df['Predicted']-gbrt_df['Actual']<0,"Final_Predict"]=gbrt_df['Predicted']-(gbrt_df['Predicted']-gbrt_df['Actual'])
print(gbrt_df)

# https://scikit-learn.org/stable/auto_examples/ensemble/plot_voting_regressor.html#sphx-glr-auto-examples-ensemble-plot-voting-regressor-py

plt.figure()
plt.plot(gbrt_df['Actual'],'gd', label='Actual_price')
plt.plot(gbrt_df['Predicted'],'b^', label='Predict_price')
plt.plot(gbrt_df['Final_Predict'],'ys', label='Final_Predict_price')
plt.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False)

plt.ylabel('Prieces')
plt.xlabel('Feature_S')
plt.legend(loc="best")
plt.show()


text:
#!/usr/bin/env python3

import boost_histogram as bh
import numpy as np
import matplotlib.pyplot as plt

# Create 2d-histogram with two axes with 10 equidistant bins from -3 to 3
h = bh.Histogram(
    bh.axis.Regular(10, -3, 3, metadata="x"), bh.axis.Regular(10, -3, 3, metadata="y")
)

# Generate some Numpy arrays with data to fill into histogram,
# in this case normal distributed random numbers in x and y
x_data = np.random.randn(1000)
y_data = 0.5 * np.random.randn(1000)

# Fill histogram with numpy arrays, this is very fast
h.fill(x_data, y_data)

# Get representations of the bin edges as Numpy arrays
x = h.axes[0].edges
y = h.axes[1].edges

# Creates a view of the counts (no copy involved)
count_matrix = h.view()

# Draw the count matrix
plt.pcolor(x, y, count_matrix.T)
plt.xlabel(h.axes[0].metadata)
plt.ylabel(h.axes[1].metadata)
plt.savefig("simple_numpy.png")

max_stars_repo_name: HDembinski/boost-histogram

max_issues_repo_name: HDembinski/boost-histogram

max_forks_repo_name: HDembinski/boost-histogram