<commit_before># -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2021, GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.
import os
import unittest
import numpy as np
from numpy.testing import assert_allclose as aac
import pandas
from openquake.hazardlib import read_input, valid
from openquake.hazardlib.cross_correlation import BakerJayaram2008
from openquake.hazardlib.calc.filters import IntegrationDistance
from openquake.hazardlib.calc.cond_spectra import get_cs_out, cond_spectra

OVERWRITE_EXPECTED = False

CWD = os.path.dirname(__file__)
SOURCES_XML = os.path.join(CWD, 'data', 'sm01.xml')
GSIM_XML = os.path.join(CWD, 'data', 'lt02.xml')
PARAM = dict(source_model_file=SOURCES_XML,
             gsim_logic_tree_file=GSIM_XML,
             sites=[(0, -0.8)],
             reference_vs30_value=600,
             reference_depth_to_2pt5km_per_sec=5,
             reference_depth_to_1pt0km_per_sec=100,
             maximum_distance=IntegrationDistance.new('200'),
             rupture_mesh_spacing=5.,
             width_of_mfd_bin=1.,
             investigation_time=1,
             truncation_level=3,
             cross_correl=BakerJayaram2008(),
             imtls={"SA(0.05)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.075)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.1)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.15)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.2)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.25)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.3)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.5)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.75)": valid.logscale(0.005, 2.13, 45),
                    "SA(1.0)": valid.logscale(0.005, 2.13, 45),
                    "SA(2.0)": valid.logscale(0.005, 2.13, 45)})
imti = 4  # corresponds to SA(0.2)
poes = [0.000404]


# useful while debugging
def plot(df, imts):
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(1, 2)
    axs[0].plot(df.period[:11], df.cs_exp[:11], 'x-')
    axs[0].plot(df.period[11:], df.cs_exp[11:], 'x-')
    axs[1].plot(df.period[:11], df.cs_std[:11], 'x-')
    axs[1].plot(df.period[11:], df.cs_std[11:], 'x-')
    axs[0].grid(which='both')
    axs[1].grid(which='both')
    axs[0].set_xscale('log')
    axs[1].set_xscale('log')
    axs[0].set_yscale('log')
    axs[0].set_ylim([1e-3, 10])
    axs[0].set_xlabel('Mean spectrum, Period[s]')
    axs[1].set_xlabel('Std spectrum, Period[s]')
    plt.show()


# used to create the expected file the first time
def outdic_to_dframe(outdic, imts, n, p):
    """
    :param outdic: a dictionary g_ -> array
    :param imts: M intensity measure types
    :param rlzs: R realization indices
    :param n: an index in the range 0..N-1 where N is the number of sites
    :param p: an index in the range 0..P-1 where P is the number of IMLs
    """
    dic = dict(rlz_id=[], period=[], cs_exp=[], cs_std=[])
    for r, c in outdic.items():
        for m, imt in enumerate(imts):
            dic['rlz_id'].append(r)
            dic['period'].append(imt.period)
            dic['cs_exp'].append(np.exp(c[m, n, 1, p] / c[m, n, 0, p]))
            dic['cs_std'].append(np.sqrt(c[m, n, 2, p] / c[m, n, 0, p]))
    return pandas.DataFrame(dic)


class CondSpectraTestCase(unittest.TestCase):

    def test_1_rlz(self):
        # test with one GMPE, 1 TRT, checking additivity
        inp = read_input(
            PARAM, gsim_logic_tree_file=os.path.join(CWD, 'data', 'lt01.xml'))
        [cmaker] = inp.cmakerdict.values()
        [src_group] = inp.groups
        [ctx] = cmaker.from_srcs(src_group, inp.sitecol)
        assert len(ctx) == 100
        ctx1 = ctx[:50]
        ctx2 = ctx[50:]

        # The hazard for the target IMT and poe
        cmaker.poes = [0.000404]
        imls = [0.394359437]

        mom1 = get_cs_out(cmaker, ctx1, imti, imls)[0]
        mom2 = get_cs_out(cmaker, ctx2, imti, imls)[0]
        mom = get_cs_out(cmaker, ctx, imti, imls)[0]
        aac(mom1 + mom2, mom)

        spectra, s_sigma = cond_spectra(
            cmaker, src_group, inp.sitecol, 'SA(0.2)', imls)
        aac(spectra.flatten(), [0.19236242, 0.23961989, 0.27838065, 0.35216192,
                                0.39435944, 0.36501786, 0.34676928, 0.23458421,
                                0.15669297, 0.11154595, 0.0409729], atol=1E-5)
        aac(s_sigma.flatten(), [0.327456, 0.368969, 0.388289, 0.270122,
                                0.006058, 0.235236, 0.319312, 0.463179,
                                0.556208, 0.596132, 0.70371], atol=1E-5)

    def test_2_rlzs(self):
        # test with two GMPEs, 1 TRT
        inp = read_input(PARAM)
        [cmaker] = inp.cmakerdict.values()
        [src_group] = inp.groups<commit_msg>Raised tolerance<commit_after># -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2021, GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.
import os
import unittest
import numpy as np
from numpy.testing import assert_allclose as aac
import pandas
from openquake.hazardlib import read_input, valid
from openquake.hazardlib.cross_correlation import BakerJayaram2008
from openquake.hazardlib.calc.filters import IntegrationDistance
from openquake.hazardlib.calc.cond_spectra import get_cs_out, cond_spectra

OVERWRITE_EXPECTED = False

CWD = os.path.dirname(__file__)
SOURCES_XML = os.path.join(CWD, 'data', 'sm01.xml')
GSIM_XML = os.path.join(CWD, 'data', 'lt02.xml')
PARAM = dict(source_model_file=SOURCES_XML,
             gsim_logic_tree_file=GSIM_XML,
             sites=[(0, -0.8)],
             reference_vs30_value=600,
             reference_depth_to_2pt5km_per_sec=5,
             reference_depth_to_1pt0km_per_sec=100,
             maximum_distance=IntegrationDistance.new('200'),
             rupture_mesh_spacing=5.,
             width_of_mfd_bin=1.,
             investigation_time=1,
             truncation_level=3,
             cross_correl=BakerJayaram2008(),
             imtls={"SA(0.05)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.075)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.1)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.15)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.2)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.25)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.3)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.5)": valid.logscale(0.005, 2.13, 45),
                    "SA(0.75)": valid.logscale(0.005, 2.13, 45),
                    "SA(1.0)": valid.logscale(0.005, 2.13, 45),
                    "SA(2.0)": valid.logscale(0.005, 2.13, 45)})
imti = 4  # corresponds to SA(0.2)
poes = [0.000404]


# useful while debugging
def plot(df, imts):
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(1, 2)
    axs[0].plot(df.period[:11], df.cs_exp[:11], 'x-')
    axs[0].plot(df.period[11:], df.cs_exp[11:], 'x-')
    axs[1].plot(df.period[:11], df.cs_std[:11], 'x-')
    axs[1].plot(df.period[11:], df.cs_std[11:], 'x-')
    axs[0].grid(which='both')
    axs[1].grid(which='both')
    axs[0].set_xscale('log')
    axs[1].set_xscale('log')
    axs[0].set_yscale('log')
    axs[0].set_ylim([1e-3, 10])
    axs[0].set_xlabel('Mean spectrum, Period[s]')
    axs[1].set_xlabel('Std spectrum, Period[s]')
    plt.show()


# used to create the expected file the first time
def outdic_to_dframe(outdic, imts, n, p):
    """
    :param outdic: a dictionary g_ -> array
    :param imts: M intensity measure types
    :param rlzs: R realization indices
    :param n: an index in the range 0..N-1 where N is the number of sites
    :param p: an index in the range 0..P-1 where P is the number of IMLs
    """
    dic = dict(rlz_id=[], period=[], cs_exp=[], cs_std=[])
    for r, c in outdic.items():
        for m, imt in enumerate(imts):
            dic['rlz_id'].append(r)
            dic['period'].append(imt.period)
            dic['cs_exp'].append(np.exp(c[m, n, 1, p] / c[m, n, 0, p]))
            dic['cs_std'].append(np.sqrt(c[m, n, 2, p] / c[m, n, 0, p]))
    return pandas.DataFrame(dic)


class CondSpectraTestCase(unittest.TestCase):

    def test_1_rlz(self):
        # test with one GMPE, 1 TRT, checking additivity
        inp = read_input(
            PARAM, gsim_logic_tree_file=os.path.join(CWD, 'data', 'lt01.xml'))
        [cmaker] = inp.cmakerdict.values()
        [src_group] = inp.groups
        [ctx] = cmaker.from_srcs(src_group, inp.sitecol)
        assert len(ctx) == 100
        ctx1 = ctx[:50]
        ctx2 = ctx[50:]

        # The hazard for the target IMT and poe
        cmaker.poes = [0.000404]
        imls = [0.394359437]

        mom1 = get_cs_out(cmaker, ctx1, imti, imls)[0]
        mom2 = get_cs_out(cmaker, ctx2, imti, imls)[0]
        mom = get_cs_out(cmaker, ctx, imti, imls)[0]
        aac(mom1 + mom2, mom)

        spectra, s_sigma = cond_spectra(
            cmaker, src_group, inp.sitecol, 'SA(0.2)', imls)
        aac(spectra.flatten(), [0.19236242, 0.23961989, 0.27838065, 0.35216192,
                                0.39435944, 0.36501786, 0.34676928, 0.23458421,
                                0.15669297, 0.11154595, 0.0409729], atol=2E-5)
        aac(s_sigma.flatten(), [0.327456, 0.368969, 0.388289, 0.270122,
                                0.006058, 0.235236, 0.319312, 0.463179,
                                0.556208, 0.596132, 0.70371], atol=1E-5)

    def test_2_rlzs(self):
        # test with two GMPEs, 1 TRT
        inp = read_input(PARAM)
        [cmaker] = inp.cmakerdict.values()
        [src_group] = inp.groups
{'fraction_non_alphanumeric': 0.14104088539284734, 'fraction_numerical': 0.0812219884752584, 'mean_word_length': 2.5, 'pattern_counts': {'":': 22, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '366686', 'n_tokens_mistral': 4546, 'n_tokens_neox': 3992, 'n_words': 936}