text:
using StochasticGroundMotionSimulation
using Test
using ForwardDiff
using ForwardDiff: Dual
using FastGaussQuadrature
using QuadGK
using LinearAlgebra

@testset "StochasticGroundMotionSimulation.jl" begin

    @testset "Performance" begin
        Ti = [0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 5.0, 7.5, 10.0]
        m = 4.0 + π
        r = 500.0 + π

        src = SourceParameters(100.0)
        geo = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5])
        ane = AnelasticAttenuationParameters(200.0, 0.4)
        sat = NearSourceSaturationParameters(:BT15)
        path = PathParameters(geo, sat, ane)
        site = SiteParameters(0.039)

        fas = FourierParameters(src, path, site)
        rvt = RandomVibrationParameters(:DK80)

        sdof = Oscillator(1.0)
        # @code_warntype Oscillator(1.0)
        # @code_warntype period(sdof)
        # @code_warntype transfer(1.0, sdof)

        # @code_warntype rvt_response_spectral_ordinate(Ti[1], m, r, fas, rvt)
        # @code_warntype rvt_response_spectrum(Ti, m, r, fas, rvt)
        # @time Sai = rvt_response_spectrum(Ti, m, r, fas, rvt)
        # @profile Sai = rvt_response_spectrum(Ti, m, r, fas, rvt)

    end

    @testset "Source" begin

        m = 6.0
        Δσ = 100.0
        β = 3.5

        # @code_warntype magnitude_to_moment(m)
        #
        # @code_warntype corner_frequency_brune(m, Δσ)
        # @code_warntype corner_frequency_brune(m, Δσ, β)
        # @code_warntype corner_frequency_atkinson_silva_2000(m)

        srcf = SourceParameters(Δσ)
        srcd = SourceParameters(Dual{Float64}(Δσ))

        # @code_warntype corner_frequency(m, srcf)
        # @code_warntype corner_frequency(m, srcd)
        # @code_warntype corner_frequency(Dual(m), srcf)
        # @code_warntype corner_frequency(Dual(m), srcd)

        T = StochasticGroundMotionSimulation.get_parametric_type(srcf)
        @test T == Float64
        T = StochasticGroundMotionSimulation.get_parametric_type(srcd)
        @test T <: Dual

        faf, fbf, fεf = corner_frequency(m, srcf)
        fad, fbd, fεd = corner_frequency(m, srcd)
        # @time faf, fbf, fεf = corner_frequency(m, srcf)
        # @time fad, fbd, fεd = corner_frequency(m, srcd)

        @test faf == fad.value

        srcf = SourceParameters(Δσ, :Atkinson_Silva_2000)
        srcd = SourceParameters(Dual{Float64}(Δσ), :Atkinson_Silva_2000)

        faf, fbf, fεf = corner_frequency(m, srcf)
        fad, fbd, fεd = corner_frequency(m, srcd)
        # @time faf, fbf, fεf = corner_frequency(m, srcf)
        # @time fad, fbd, fεd = corner_frequency(m, srcd)

        @test faf == fad.value
        @test fbf == fbd.value
        @test fεf == fεd.value

        srcd = SourceParameters(Dual(100.0), 3.5, 2.75)
        srcf = SourceParameters(100.0, 3.5, 2.75)
        @test srcd.Δσ.value == srcf.Δσ

        @test StochasticGroundMotionSimulation.magnitude_to_moment(6.0) ≈ exp10(25.05)

    end

    @testset "Path" begin

        Rrefi = [1.0, 50.0, Inf]
        γi = [1.0, 0.5]

        @testset "Geometric Spreading Constructors" begin
            # test floating point instantiation
            @test typeof(GeometricSpreadingParameters(Rrefi, γi)) <: GeometricSpreadingParameters
            # test Dual instantiation
            @test typeof(GeometricSpreadingParameters(Rrefi, [0.5], [Dual{Float64}(1.0)], BitVector([1, 0]), :Piecewise)) <: GeometricSpreadingParameters

            geo = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5])
            geoc = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5], :CY14)
            @test geo.model == :Piecewise
            @test geoc.model == :CY14
            geod = GeometricSpreadingParameters([1.0, 50.0, Inf], [Dual(1.0), Dual(0.5)])
            @test geo.γconi[1] == geod.γvari[1].value

        end

        geof = GeometricSpreadingParameters(Rrefi, γi)
        geod = GeometricSpreadingParameters(Rrefi, [0.5], [Dual{Float64}(1.0)], BitVector([1, 0]), :Piecewise)

        @testset "Geometric Spreading Types" begin
            T = StochasticGroundMotionSimulation.get_parametric_type(geof)
            @test T == Float64
            T = StochasticGroundMotionSimulation.get_parametric_type(geod)
            @test T <: Dual
        end

        @testset "Near Source Saturation Constructors" begin
            # test standard instantiation with known models
            @test typeof(NearSourceSaturationParameters(:BT15)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(:YA15)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(:CY14)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(:CY14mod)) <: NearSourceSaturationParameters

            @test typeof(NearSourceSaturationParameters(1, :BT15)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters([5.5, 7.0, Inf], [4.0, 6.0], :ConstantConstrained)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters([5.5, 7.0, Inf], [Dual(4.0), Dual(6.0)], :ConstantVariable)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters([5.5, 7.0, Inf], [4.0, 6.0])) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters([5.5, 7.0, Inf], [Dual(4.0), Dual(6.0)])) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(5.0)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(Dual(5.0))) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(5.0, 2)) <: NearSourceSaturationParameters
            @test typeof(NearSourceSaturationParameters(Dual(5.0), 2)) <: NearSourceSaturationParameters

        end

        sat = NearSourceSaturationParameters(:BT15)
        satd = NearSourceSaturationParameters(Dual{Float64}(5.0))

        @testset "Near Source Saturation Types" begin
            T = StochasticGroundMotionSimulation.get_parametric_type(sat)
            @test T == Float64
        end

        @testset "Anelastic Attenuation Constructors" begin
            # test floating point instantiation
            @test typeof(AnelasticAttenuationParameters(200.0)) <: AnelasticAttenuationParameters
            # test Dual instantiation
            @test typeof(AnelasticAttenuationParameters(Dual{Float64}(200.0))) <: AnelasticAttenuationParameters

            ane = AnelasticAttenuationParameters(200.0, 0.5, 3.5)
            @test ane.rmetric == :Rps
            ane = AnelasticAttenuationParameters(200.0, 0.5, 3.5, :Rrup)
            @test ane.rmetric == :Rrup

            # test the segmented versions
            # scalar inputs (internally mapped to vectors)
            @test typeof(AnelasticAttenuationParameters(200.0, 0.5, :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters(200.0, Dual{Float64}(0.5), :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters(Dual{Float64}(200.0), 0.5, :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters(Dual{Float64}(200.0), Dual{Float64}(0.5), :Rrup)) <: AnelasticAttenuationParameters
            # vector inputs
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)], :Rrup)) <: AnelasticAttenuationParameters

            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [0.5, 0.5])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [0.5, 0.5], :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)], [Dual{Float64}(0.5), Dual{Float64}(0.5)])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)], [Dual{Float64}(0.5), Dual{Float64}(0.5)], :Rrup)) <: AnelasticAttenuationParameters

            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [Dual{Float64}(0.5), Dual{Float64}(0.5)])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [Dual{Float64}(0.5), Dual{Float64}(0.5)], :Rrup)) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)], [0.5, 0.5])) <: AnelasticAttenuationParameters
            @test typeof(AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual{Float64}(200.0), Dual{Float64}(200.0)], [0.5, 0.5], :Rrup)) <: AnelasticAttenuationParameters

        end

        Q0 = 200.0
        anef = AnelasticAttenuationParameters(Q0)
        aned = AnelasticAttenuationParameters(Dual{Float64}(Q0))

        @testset "Anelastic Attenuation Types" begin
            T = StochasticGroundMotionSimulation.get_parametric_type(anef)
            @test T == Float64
            T = StochasticGroundMotionSimulation.get_parametric_type(aned)
            @test T <: Dual

            # artificial test, but for ensuring complete coverage
            anet = AnelasticAttenuationParameters([Dual{Float64}(0.0), Dual{Float64}(Inf)], [Dual{Float64}(200.0)], Vector{Dual{Float64,Float64,0}}(), zeros(Dual{Float64,Float64,0}, 1), zeros(Dual{Float64,Float64,0}, 1), [3.5], BitVector(ones(1)), BitVector(ones(1)), :Rps)
            T = StochasticGroundMotionSimulation.get_parametric_type(anet)
            @test T <: Dual
        end

        @testset "Path Constructors" begin
            # test floating point components
            @test typeof(PathParameters(geof, sat, anef)) <: PathParameters
            # test Dual components
            @test typeof(PathParameters(geod, sat, aned)) <: PathParameters

            path = PathParameters(geof, anef)
            @test path.saturation.model == :None

            path = PathParameters(geof, satd, anef)
            @test StochasticGroundMotionSimulation.get_parametric_type(path) <: Dual

        end

        pathf = PathParameters(geof, sat, anef)
        pathd = PathParameters(geod, sat, aned)

        @testset "Path Types" begin
            T = StochasticGroundMotionSimulation.get_parametric_type(pathf)
            @test T == Float64
            T = StochasticGroundMotionSimulation.get_parametric_type(pathd)
            @test T <: Dual

        end

        r = 10.0
        m = 6.0

        @testset "Geometric Spreading Functionality" begin

            geo_cy14 = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5], :CY14)
            geo_cy14mod = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5], :CY14mod)
            geo_null = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5], :Null)

            sat = NearSourceSaturationParameters(:BT15)
            r_ps = equivalent_point_source_distance(r, m, sat)

            gr_cy14 = geometric_spreading(r_ps, m, geo_cy14, sat)
            gr_cy14mod = geometric_spreading(r_ps, m, geo_cy14mod, sat)
            @test gr_cy14mod < gr_cy14
            @test isnan(geometric_spreading(r_ps, m, geo_null, sat))

            grp = geometric_spreading(r_ps, pathf)
            fasf = FourierParameters(SourceParameters(50.0), pathf)
            grf = geometric_spreading(r_ps, fasf)
            @test grp == grf

            geop = GeometricSpreadingParameters([1.0, Inf], [1.0], Vector{Float64}(), BitVector(undef, 0), :Piecewise)
            @test StochasticGroundMotionSimulation.geometric_spreading_piecewise(r_ps, geop) == 1.0
            @test StochasticGroundMotionSimulation.geometric_spreading_piecewise(Dual(r_ps), geop) == 1.0
            fasp = FourierParameters(SourceParameters(50.0), PathParameters(geop, sat, anef))
            @test StochasticGroundMotionSimulation.geometric_spreading_piecewise(r_ps, fasp) == 1.0
            @test StochasticGroundMotionSimulation.geometric_spreading_piecewise(Dual(r_ps), fasp) == 1.0

            geo_cy14d = GeometricSpreadingParameters([1.0, 50.0, Inf], [Dual(1.0), Dual(1.0)], :CY14)
            sat = NearSourceSaturationParameters(:None)
            r_ps = equivalent_point_source_distance(1.0, -5.0, sat)
            fas_cy14d = FourierParameters(SourceParameters(50.0), PathParameters(geo_cy14d, sat, anef))

            @test StochasticGroundMotionSimulation.geometric_spreading_cy14(r_ps, fas_cy14d).value ≈ 1.0

            geo_cy14d = GeometricSpreadingParameters([1.0, 50.0, Inf], [Dual(1.0), Dual(1.0)], :CY14mod)
            sat = NearSourceSaturationParameters(:None)
            r_ps = equivalent_point_source_distance(1.0, -5.0, sat)
            fas_cy14d = FourierParameters(SourceParameters(50.0), PathParameters(geo_cy14d, sat, anef))

            @test StochasticGroundMotionSimulation.geometric_spreading_cy14mod(r_ps, -5.0, fas_cy14d).value ≈ 1.0

            # @code_warntype equivalent_point_source_distance(r, m, pathf)
            # @code_warntype equivalent_point_source_distance(r, m, pathd)
            r_psf = equivalent_point_source_distance(r, m, pathf)
            r_psd = equivalent_point_source_distance(r, m, pathd)

            # @code_warntype StochasticGroundMotionSimulation.geometric_spreading_piecewise(r, geof)
            # @code_warntype StochasticGroundMotionSimulation.geometric_spreading_piecewise(r, geod)
            grf = StochasticGroundMotionSimulation.geometric_spreading_piecewise(r, geof)
            grd = StochasticGroundMotionSimulation.geometric_spreading_piecewise(r, geod)
            @test grf == grd.value

            # @code_warntype StochasticGroundMotionSimulation.geometric_spreading_cy14(r, geof)
            # @code_warntype StochasticGroundMotionSimulation.geometric_spreading_cy14(r, geod)
            grf = StochasticGroundMotionSimulation.geometric_spreading_cy14(r, geof)
            grd = StochasticGroundMotionSimulation.geometric_spreading_cy14(r, geod)
            @test grf == grd.value

            # @code_warntype geometric_spreading(r, geof)
            # @code_warntype geometric_spreading(r, geod)

            grf = geometric_spreading(r, m, geof, sat)
            grd = geometric_spreading(r, m, geod, sat)
            @test grf == grd.value
        end

        @testset "Near Source Saturation Functionality" begin

            # @code_warntype near_source_saturation(m, pathf.saturation)
            # @code_warntype near_source_saturation(m, pathd.saturation)
            hf = near_source_saturation(m, pathf.saturation)
            hd = near_source_saturation(m, pathd.saturation)
            if StochasticGroundMotionSimulation.get_parametric_type(pathd.saturation) <: Dual
                @test hf == hd.value
            else
                @test hf == hd
            end

            src = SourceParameters(100.0)
            geo = GeometricSpreadingParameters([1.0, Inf], [1.0])
            ane = AnelasticAttenuationParameters(200.0, 0.5)
            sat_ya = NearSourceSaturationParameters(:YA15)
            sat_cy = NearSourceSaturationParameters(:CY14)
            sat_none = NearSourceSaturationParameters(:None)
            sat_con = NearSourceSaturationParameters(5.0)
            sat_var = NearSourceSaturationParameters(Dual(5.0))
            sat_null = NearSourceSaturationParameters(:Null)

            path_ya = PathParameters(geo, sat_ya, ane)
            path_cy = PathParameters(geo, sat_cy, ane)
            path_none = PathParameters(geo, sat_none, ane)
            path_con = PathParameters(geo, sat_con, ane)
            path_var = PathParameters(geo, sat_var, ane)
            path_null = PathParameters(geo, sat_null, ane)

            h_ya = near_source_saturation(5.0, sat_ya)
            h_cy = near_source_saturation(5.0, sat_cy)
            h_none = near_source_saturation(5.0, sat_none)
            h_con = near_source_saturation(5.0, sat_con)
            h_var = near_source_saturation(5.0, sat_var)
            h_null = near_source_saturation(5.0, sat_null)

            @test h_con == h_var.value
        end

        f = 1.0
        r = 100.0

        @testset "Anelastic Attenuation Functionality" begin
            # @code_warntype anelastic_attenuation(f, r, anef)
            # @code_warntype anelastic_attenuation(f, r, aned)
            qrf = anelastic_attenuation(f, r, anef)
            qrd = anelastic_attenuation(f, r, aned)
            @test qrf ≈ qrd.value

            fas = FourierParameters(SourceParameters(100.0), pathf)
            @test fas.site.model == :Unit

            q_r = anelastic_attenuation(f, r, fas)
            @test qrf ≈ q_r
        end

        @testset "Anelastic Attenuation Segmentation" begin
            # test the vector descriptions of anelastic attenuation
            ane_vec = AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [0.4, 0.4])
            ane_con = AnelasticAttenuationParameters(200.0, 0.4)

            @test anelastic_attenuation(5.0, 50.0, ane_vec) ≈ anelastic_attenuation(5.0, 50.0, ane_con)
            @test anelastic_attenuation(5.0, 150.0, ane_vec) ≈ anelastic_attenuation(5.0, 150.0, ane_con)

            ane_vecd = AnelasticAttenuationParameters([0.0, 80.0, Inf], [Dual(200.0), Dual(200.0)], [0.4, 0.4])
            ane_cond = AnelasticAttenuationParameters(Dual(200.0), 0.4)

            @test anelastic_attenuation(5.0, 50.0, ane_vecd) ≈ anelastic_attenuation(5.0, 50.0, ane_cond)
            @test anelastic_attenuation(5.0, 150.0, ane_vecd) ≈ anelastic_attenuation(5.0, 150.0, ane_cond)

            ane_vecd = AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0, 200.0], [Dual(0.4), Dual(0.4)])
            ane_cond = AnelasticAttenuationParameters(200.0, Dual(0.4))

            @test anelastic_attenuation(5.0, 50.0, ane_vecd) ≈ anelastic_attenuation(5.0, 50.0, ane_cond)
            @test anelastic_attenuation(5.0, 150.0, ane_vecd) ≈ anelastic_attenuation(5.0, 150.0, ane_cond)

            ane_vecd = AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0], [Dual(200.0)], [0.4], [Dual(0.4)], 3.5 * ones(2), BitVector([1, 0]), BitVector([0, 1]), :Rrup)
            ane_cond = AnelasticAttenuationParameters(200.0, 0.4)

            @test anelastic_attenuation(5.0, 50.0, ane_vecd) ≈ anelastic_attenuation(5.0, 50.0, ane_cond)
            @test anelastic_attenuation(5.0, 150.0, ane_vecd) ≈ anelastic_attenuation(5.0, 150.0, ane_cond)

            ane_vecd = AnelasticAttenuationParameters([0.0, 80.0, Inf], [200.0], [Dual(200.0)], [0.4], [Dual(0.4)], 3.5 * ones(2), BitVector([0, 1]), BitVector([1, 0]), :Rrup)
            ane_cond = AnelasticAttenuationParameters(200.0, 0.4)

            @test anelastic_attenuation(5.0, 50.0, ane_vecd) ≈ anelastic_attenuation(5.0, 50.0, ane_cond)
            @test anelastic_attenuation(5.0, 150.0, ane_vecd) ≈ anelastic_attenuation(5.0, 150.0, ane_cond)

            ane_inf = AnelasticAttenuationParameters([0.0, 80.0, Inf], [Inf], [Dual{Float64}(200.0)], [0.0], [Dual{Float64}(0.5)], 3.5 * ones(2), BitVector([0, 1]), BitVector([0, 1]), :Rrup)
            @test anelastic_attenuation(5.0, 50.0, ane_inf) == 1.0
        end
    end

    @testset "Site" begin

        κ0f = 0.039
        κ0d = Dual{Float64}(κ0f)
        ζ0f = 0.039
        ζ0d = Dual{Float64}(ζ0f)
        ηf = 0.75
        ηd = Dual{Float64}(ηf)

        @testset "Site Constructors" begin
            @test typeof(SiteParameters(κ0f)) <: SiteParameters
            @test typeof(SiteParameters(κ0f, :AlAtik2021_cy14_760)) <: SiteParameters
            @test typeof(SiteParameters(κ0f, :Boore2016)) <: SiteParameters
            @test typeof(SiteParameters(κ0f, :Unit)) <: SiteParameters
            @test typeof(SiteParameters(κ0f, :NaN)) <: SiteParameters

            @test typeof(SiteParameters(κ0d)) <: SiteParameters
            @test typeof(SiteParameters(κ0d, :AlAtik2021_cy14_760)) <: SiteParameters
            @test typeof(SiteParameters(κ0d, :Boore2016)) <: SiteParameters
            @test typeof(SiteParameters(κ0d, :Unit)) <: SiteParameters
            @test typeof(SiteParameters(κ0d, :NaN)) <: SiteParameters

            @test typeof(SiteParameters(ζ0f, ηf)) <: SiteParameters
            @test typeof(SiteParameters(ζ0d, ηf)) <: SiteParameters
            @test typeof(SiteParameters(ζ0f, ηd)) <: SiteParameters
            @test typeof(SiteParameters(ζ0d, ηd)) <: SiteParameters
            @test typeof(SiteParameters(ζ0f, ηf, :AlAtik2021_cy14_760)) <: SiteParameters
            @test typeof(SiteParameters(ζ0d, ηf, :AlAtik2021_cy14_760)) <: SiteParameters
            @test typeof(SiteParameters(ζ0f, ηd, :AlAtik2021_cy14_760)) <: SiteParameters
            @test typeof(SiteParameters(ζ0d, ηd, :AlAtik2021_cy14_760)) <: SiteParameters

        end

        site0f = SiteParameters(κ0f)
        siteAf = SiteParameters(κ0f, :AlAtik2021_cy14_760)
        siteBf = SiteParameters(κ0f, :Boore2016)
        siteUf = SiteParameters(κ0f, :Unit)
        siteNf = SiteParameters(κ0f, :NaN)

        site0d = SiteParameters(κ0d)
        siteAd = SiteParameters(κ0d, :AlAtik2021_cy14_760)
        siteBd = SiteParameters(κ0d, :Boore2016)
        siteUd = SiteParameters(κ0d, :Unit)
        siteNd = SiteParameters(κ0d, :NaN)

        f = 0.05

        @testset "Site Amplification" begin
            # @code_warntype site_amplification(f, site0f)
            # @code_warntype site_amplification(f, site0d)
            Sff = site_amplification(f, site0f)
            Sfd = site_amplification(f, site0d)
            @test Sff == Sfd

            Af0f = site_amplification(f, site0f)
            Af1f = site_amplification(f, siteAf)
            Af2f = site_amplification(f, siteBf)
            Af3f = site_amplification(f, siteUf)
            Af4f = site_amplification(f, siteNf)

            @test Af0f == Af1f
            @test Af3f == 1.0
            @test Af2f < Af1f
            @test isnan(Af4f)

            Af0d = site_amplification(f, site0d)
            Af1d = site_amplification(f, siteAd)
            Af2d = site_amplification(f, siteBd)
            Af3d = site_amplification(f, siteUd)
            Af4d = site_amplification(f, siteNd)

            @test Af0d == Af1d
            @test Af3d == 1.0
            @test Af2d < Af1d
            @test isnan(Af4d)

            @test Af0f == Af0d

            f0 = 80.0
            f1 = 100.0
            Af0 = site_amplification(f0, siteBf)
            Af1 = site_amplification(f1, siteBf)

            @test Af0 == Af1

            ft = 10.0
            Af620 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_ask14_620))
            Af760 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_ask14_760))
            Af1100 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_ask14_1100))

            @test Af620 > Af760
            @test Af760 > Af1100

            Af620 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_bssa14_620))
            Af760 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_bssa14_760))
            Af1100 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_bssa14_1100))

            @test Af620 > Af760
            @test Af760 > Af1100

            Af620 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cb14_620))
            Af760 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cb14_760))
            Af1100 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cb14_1100))

            @test Af620 > Af760
            @test Af760 > Af1100

            Af620 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cy14_620))
            Af760 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cy14_760))
            Af1100 = site_amplification(ft, SiteParameters(0.039, :AlAtik2021_cy14_1100))

            @test Af620 > Af760
            @test Af760 > Af1100
        end

        @testset "Impedance Functions" begin
            @test StochasticGroundMotionSimulation.boore_2016_generic_amplification(0.015) == 1.01
            @test isnan(StochasticGroundMotionSimulation.boore_2016_generic_amplification(NaN))

            numf = length(StochasticGroundMotionSimulation.fii_b16_760)
            for i = 1:numf
                fi = StochasticGroundMotionSimulation.fii_b16_760[i]
                @test StochasticGroundMotionSimulation.boore_2016_generic_amplification(fi) == StochasticGroundMotionSimulation.Aii_b16_760[i]
            end

            numf = length(StochasticGroundMotionSimulation.fii_aa21_cy14_760)
            for i = 1:numf
                fi = StochasticGroundMotionSimulation.fii_aa21_cy14_760[i]
                @test StochasticGroundMotionSimulation.itp_aa21_cy14_760(log(fi)) == StochasticGroundMotionSimulation.Aii_aa21_cy14_760[i]
            end
        end

        @testset "Kappa Filter" begin
            f = 10.0
            # @code_warntype kappa_filter(f, siteAf)
            # @code_warntype kappa_filter(f, siteAd)
            Kff = kappa_filter(f, siteAf)
            Kfd = kappa_filter(f, siteAd)
            @test Kff == Kfd.value
        end

        @testset "Zeta Filter" begin
            f = 10.0
            # @code_warntype kappa_filter(f, siteAf)
            # @code_warntype kappa_filter(f, siteAd)
            siteAzf = SiteParameters(ζ0f, ηf)
            siteAzd = SiteParameters(ζ0d, ηd)
            Kff = kappa_filter(f, siteAzf)
            Kfd = kappa_filter(f, siteAzd)
            @test Kff == Kfd.value

            # fi = [1.0, 5.0, 10.0, 15.0, 20.0, 30.0]
            # [map(f -> kappa_filter(f, SiteParameters(0.04)), fi) map(f -> kappa_filter(f, SiteParameters(0.04, 0.1)), fi) ]

        end

    end

    @testset "Oscillator" begin
        ζ = 0.05
        f_n = 1.0
        sdof = Oscillator(f_n, ζ)

        @test f_n ≈ 1.0 / period(sdof)

        @test transfer(0.5, sdof)^2 ≈ StochasticGroundMotionSimulation.squared_transfer(0.5, sdof)

        fi = [0.5, 1.0, 2.0]

        # @code_warntype transfer(fi, sdof)

        Hfi = transfer(fi, sdof)
        tfi = 2 * fi
        transfer!(Hfi, tfi, sdof)
        @test Hfi ≈ transfer(tfi, sdof)

        StochasticGroundMotionSimulation.squared_transfer!(Hfi, fi, sdof)
        @test Hfi ≈ transfer(fi, sdof) .^ 2

    end

    @testset "Duration" begin

        src = SourceParameters(100.0)
        geo = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5])
        ane = AnelasticAttenuationParameters(200.0, 0.4)
        sat = NearSourceSaturationParameters(:BT15)
        path = PathParameters(geo, sat, ane)
        site = SiteParameters(0.039)

        fas = FourierParameters(src, path, site)

        # Boore & Thompson 2014
        m = 6.0
        fa, fb, ε = corner_frequency(m, src)
        Ds = 1.0 / fa

        Δσf = 100.0
        Δσd = Dual{Float64}(Δσf)
        srcf = SourceParameters(Δσf)
        srcd = SourceParameters(Δσd)
        fasf = FourierParameters(srcf, path, site)
        fasd = FourierParameters(srcd, path, site)

        faf, fbf, εf = corner_frequency(m, srcf)
        # @code_warntype StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, srcf)
        Dsf = StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, srcf)
        @test Dsf ≈ 1.0 / faf
        @test isnan(fbf)
        @test isnan(εf)
        fad, fbd, εd = corner_frequency(m, srcd)
        # @code_warntype StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, srcd)
        Dsd = StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, srcd)
        @test Dsd ≈ 1.0 / fad
        @test isnan(fbd)
        @test isnan(εd)

        # @code_warntype StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, fasf)
        @test StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, fasf) ≈ Ds
        @test StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, fasd) ≈ Ds

        m = 6.0
        r = 7.0
        fa, fb, ε = corner_frequency(m, fasf)
        Dur = 1.0 / fa + 2.4
        @test StochasticGroundMotionSimulation.boore_thompson_2014(m, r, fasf) ≈ Dur
        fa, fb, ε = corner_frequency(m, fasd)
        Dur = 1.0 / fa + 2.4
        @test StochasticGroundMotionSimulation.boore_thompson_2014(m, r, fasd) ≈ Dur

        srcAS = SourceParameters(Δσf, :Atkinson_Silva_2000)
        fa, fb, ε = corner_frequency(m, srcAS)
        Ds = 0.5 * (1.0 / fa + 1.0 / fb)
        @test StochasticGroundMotionSimulation.boore_thompson_2014(m, 0.0, srcAS) ≈ Ds
        @test isnan(fb) == false
        @test isnan(ε) == false
        @test fa < fb

        # test gradient of BT14 duration model w.r.t. magnitude
        h = 0.05
        m1 = 8.0
        m2 = m1 + h
        r_ps1 = 1.0 + near_source_saturation(m1, fasf)
        r_ps2 = 1.0 + near_source_saturation(m2, fasf)
        Dex1 = StochasticGroundMotionSimulation.boore_thompson_2014(m1, r_ps1, fasf)
        Dex2 = StochasticGroundMotionSimulation.boore_thompson_2014(m2, r_ps2, fasf)
        fdg = log(Dex2 / Dex1) / h

        d(x) = log(StochasticGroundMotionSimulation.boore_thompson_2014(x[1], 1.0 + near_source_saturation(x[1], fasf), fasf))
        gd(x) = ForwardDiff.gradient(d, x)
        adg = gd([8.0])[1]

        @test fdg ≈ adg atol = 1e-2

        rvt = RandomVibrationParameters(:BT14)
        # @code_warntype excitation_duration(m, r, fasf, rvt)
        # @code_warntype excitation_duration(m, r, fasd, rvt)
        Dexf = excitation_duration(m, r, fasf, rvt)
        Dexd = excitation_duration(m, r, fasd, rvt)
        @test Dexf == Dexd.value

        c11 = [8.4312e-01, -2.8671e-02, 2.0, 1.7316e+00, 1.1695e+00, 2.1671e+00, 9.6224e-01]
        c11f = StochasticGroundMotionSimulation.StochasticGroundMotionSimulation.boore_thompson_2012_coefs(1, 1)
        @test c11f[1] == c11[1]
        @test all(isapprox.(c11, c11f))

        # @code_warntype StochasticGroundMotionSimulation.StochasticGroundMotionSimulation.boore_thompson_2012_coefs(1, 1)

        m = 8.0
        r = 1.0
        Dex = StochasticGroundMotionSimulation.boore_thompson_2014(m, r, srcf)
        # get the oscillator period
        sdof = Oscillator(100.0)
        T_n = period(sdof)
        ζ = sdof.ζ_n
        # define the η parameter as T_n/Dex
        η = T_n / Dex

        # @time c = StochasticGroundMotionSimulation.StochasticGroundMotionSimulation.boore_thompson_2012_coefs(1, 1)
        # @time StochasticGroundMotionSimulation.boore_thompson_2012_base(η, c, ζ)
        #
        # @time StochasticGroundMotionSimulation.boore_thompson_2012(m, r, srcf, sdof, rvt)
        # @code_warntype StochasticGroundMotionSimulation.boore_thompson_2012(m, r, srcf, sdof, rvt)
        #
        # @time StochasticGroundMotionSimulation.boore_thompson_2012(m, r, srcd, sdof, rvt)
        # @code_warntype StochasticGroundMotionSimulation.boore_thompson_2012(m, r, srcd, sdof, rvt)

        sdof = Oscillator(1.0)

        Drms, Dex, Dratio = StochasticGroundMotionSimulation.boore_thompson_2012(m, r, fas, sdof, rvt)
        Dex0 = StochasticGroundMotionSimulation.boore_thompson_2014(m, r, fas)

        @test Dex == Dex0

        # @code_warntype rms_duration(m, r, srcf, path, sdof, rvt)
        # @code_warntype rms_duration(m, r, srcd, path, sdof, rvt)
        # @code_warntype rms_duration(m, r, fasf, sdof, rvt)
        # @code_warntype rms_duration(m, r, fasd, sdof, rvt)

        # @time Drmsf, Dexf, Dratiof = rms_duration(m, r, fasf, sdof, rvt)
        # @time Drmsd, Dexd, Dratiod = rms_duration(m, r, fasd, sdof, rvt)
        Drmsf, Dexf, Dratiof = rms_duration(m, r, fasf, sdof, rvt)
        Drmsd, Dexd, Dratiod = rms_duration(m, r, fasd, sdof, rvt)

        @test Drmsf == Drmsd.value
        @test Dexf == Dexd.value
        @test Dratiof == Dratiod.value

        # Boore & Thompson 2014
        m = 6.0
        fa, fb, ε = corner_frequency(m, src)
        Ds = 1.0 / fa

        Dex270 = Ds + 34.2
        Dex300 = Dex270 + 0.156 * 30.0
        @test Dex270 ≈ StochasticGroundMotionSimulation.boore_thompson_2014(m, 270.0, src)
        @test Dex300 ≈ StochasticGroundMotionSimulation.boore_thompson_2014(m, 300.0, src)

        @test isnan(StochasticGroundMotionSimulation.boore_thompson_2014(m, -1.0, src))
        @test isnan(StochasticGroundMotionSimulation.boore_thompson_2014(m, -1.0, srcd))
        @test isnan(StochasticGroundMotionSimulation.boore_thompson_2014(m, Dual(-1.0), src))

        @test isnan(excitation_duration(m, -1.0, src, rvt))
        @test isnan(excitation_duration(m, -1.0, srcd, rvt))
        @test isnan(excitation_duration(m, Dual(-1.0), src, rvt))

        c = StochasticGroundMotionSimulation.StochasticGroundMotionSimulation.boore_thompson_2012_coefs(1, 1, region = :ENA)
        idx = 1
        @test all(isapprox(c, StochasticGroundMotionSimulation.coefs_ena_bt12[idx, 3:9]))

        d1a, d2a, d3a = StochasticGroundMotionSimulation.boore_thompson_2012(6.1234, 2.0, src, sdof, rvt)
        d1b, d2b, d3b = StochasticGroundMotionSimulation.boore_thompson_2012(6.1234, 2.1234, src, sdof, rvt)
        d1c, d2c, d3c = StochasticGroundMotionSimulation.boore_thompson_2012(6.0, 2.1234, src, sdof, rvt)
        d1d, d2d, d3d = StochasticGroundMotionSimulation.boore_thompson_2012(6.0, 2.0, src, sdof, rvt)
        @test d1a < d1b
        @test d2a < d2b
        @test d3a > d3b
        @test d1a > d1c
        @test d1d < d1a

        c = StochasticGroundMotionSimulation.boore_thompson_2015_coefs(1, 1, region = :ENA)
        idx = 1
        @test all(isapprox(c, StochasticGroundMotionSimulation.coefs_ena_bt15[idx, 3:9]))

        d1a, d2a, d3a = StochasticGroundMotionSimulation.boore_thompson_2015(6.1234, 2.0, src, sdof, rvt)
        d1b, d2b, d3b = StochasticGroundMotionSimulation.boore_thompson_2015(6.1234, 2.1234, src, sdof, rvt)
        d1c, d2c, d3c = StochasticGroundMotionSimulation.boore_thompson_2015(6.0, 2.1234, src, sdof, rvt)
        d1d, d2d, d3d = StochasticGroundMotionSimulation.boore_thompson_2015(6.0, 2.0, src, sdof, rvt)
        @test d1a < d1b
        @test d2a < d2b
        @test d3a > d3b
        @test d1a > d1c
        @test d1d < d1a

        d1as, d2as, d3as = StochasticGroundMotionSimulation.boore_thompson_2015(6.1234, 2.0, src, sdof, rvt)
        d1af, d2af, d3af = StochasticGroundMotionSimulation.boore_thompson_2015(6.1234, 2.0, fas, sdof, rvt)
        @test d1as == d1af
        @test d2as == d2af
        @test d3as == d3af

        rvt = RandomVibrationParameters(:PS, :PS, :PS, :PS)
        d1, d2, d3 = rms_duration(6.0, 1.0, fas, sdof, rvt)
        @test isnan(d1)
        @test isnan(d2)
        @test isnan(d3)

        srcf = SourceParameters(100.0)
        srcd = SourceParameters(Dual(100.0))
        rvt = RandomVibrationParameters(:PS, :PS, :PS, :PS)

        @test isnan(excitation_duration(6.0, 10.0, srcf, rvt))
        @test isnan(excitation_duration(6.0, 10.0, srcd, rvt))
        @test isnan(excitation_duration(6.0, Dual(10.0), srcf, rvt))

        D50 = excitation_duration(6.0, 50.0, srcf, RandomVibrationParameters())
        D150 = excitation_duration(6.0, 150.0, srcf, RandomVibrationParameters())
        @test D150 > D50

    end

    @testset "Fourier" begin

        Δσf = 100.0
        γ1f = 1.0
        γ2f = 0.5
        Q0f = 200.0
        ηf = 0.4
        κ0f = 0.039

        Δσd = Dual{Float64}(Δσf)
        γ1d = Dual{Float64}(γ1f)
        γ2d = Dual{Float64}(γ2f)
        Q0d = Dual{Float64}(Q0f)
        ηd = Dual{Float64}(ηf)
        κ0d = Dual{Float64}(κ0f)

        srcf = SourceParameters(Δσf)
        srcd = SourceParameters(Δσd)

        Rrefi = [1.0, 50.0, Inf]
        geof = GeometricSpreadingParameters(Rrefi, [γ1f, γ2f])
        geod = GeometricSpreadingParameters(Rrefi, [γ1d, γ2d])
        geom = GeometricSpreadingParameters(Rrefi, [γ1f], [γ2d], BitVector([0, 1]), :Piecewise)

        anef = AnelasticAttenuationParameters(Q0f, ηf)
        aned = AnelasticAttenuationParameters(Q0d, ηd)
        anem = AnelasticAttenuationParameters(Q0f, ηd)

        sat = NearSourceSaturationParameters(:BT15)

        pathf = PathParameters(geof, sat, anef)
        pathd = PathParameters(geod, sat, aned)
        pathm = PathParameters(geom, sat, anem)

        sitef = SiteParameters(κ0f)
        sited = SiteParameters(κ0d)

        fasf = FourierParameters(srcf, pathf, sitef)
        fasd = FourierParameters(srcd, pathd, sited)
        fasm = FourierParameters(srcf, pathm, sited)

        Tf = StochasticGroundMotionSimulation.get_parametric_type(fasf)
        Td = StochasticGroundMotionSimulation.get_parametric_type(fasd)
        Tm = StochasticGroundMotionSimulation.get_parametric_type(fasm)
        @test Tf == Float64
        @test Td <: Dual
        @test Tm <: Dual
        @test Td == Tm

        # @code_warntype fourier_constant(srcf)
        Cfs = fourier_constant(srcf)
        Cff = fourier_constant(fasf)
        @test Cfs == Cff

        # @code_warntype fourier_constant(srcd)
        Cfsd = fourier_constant(srcd)
        Cffd = fourier_constant(fasd)
        @test Cfsd == Cffd

        f = 0.001
        m = 6.0

        # @code_warntype fourier_source_shape(f, m, srcf)
        # @code_warntype fourier_source_shape(f, m, srcd)
        Affs = fourier_source_shape(f, m, srcf)
        Afff = fourier_source_shape(f, m, fasf)
        Afds = fourier_source_shape(f, m, srcd)
        Afdf = fourier_source_shape(f, m, fasd)
        @test Affs == Afds.value
        @test Afff == Afdf.value
        @test Affs == Afff
        @test Afds == Afdf

        @test Afff ≈ 1.0 atol = 1e-3

        fa, fb, ε = corner_frequency(m, srcf)
        # @code_warntype fourier_source_shape(f, fa, fb, ε, srcf.model)
        Afc = fourier_source_shape(f, fa, fb, ε, srcf.model)
        @test Afc ≈ 1.0 atol = 1e-3

        fa, fb, ε = corner_frequency(m, srcd)
        # @code_warntype fourier_source_shape(f, fa, fb, ε, srcd.model)
        Afcd = fourier_source_shape(f, fa, fb, ε, srcd.model)
        @test Afcd ≈ 1.0 atol = 1e-3

        src_a = SourceParameters(100.0, :Atkinson_Silva_2000)
        Af_a = fourier_source_shape(f, m, src_a)
        @test Af_a ≈ 1.0 atol = 1e-3
        src_n = SourceParameters(100.0, :Null)
        Af_n = fourier_source_shape(f, m, src_n)
        src_b = SourceParameters(100.0)
        Af_b = fourier_source_shape(f, m, src_b)
        @test Af_n == Af_b

        fa, fb, ε = corner_frequency(m, src_a)
        Af_a = fourier_source_shape(f, fa, fb, ε, src_a.model)
        @test Af_a ≈ 1.0 atol = 1e-3
        fa, fb, ε = corner_frequency(m, src_b)
        Af_n = fourier_source_shape(f, fa, fb, ε, src_n.model)
        @test Af_n ≈ Af_b

        # @code_warntype fourier_source(f, m, srcf)
        # @code_warntype fourier_source(f, m, srcd)
        Afs = fourier_source(f, m, srcf)
        Aff = fourier_source(f, m, fasf)
        @test Afs == Aff
        # @time fourier_source(f, m, srcd)
        # @time fourier_source(f, m, fasd)

        f = 10.0
        r = 100.0
        m = 6.0

        ane = AnelasticAttenuationParameters(200.0, 0.5, :Rrup)
        Pfr = fourier_path(f, r, m, geof, sat, ane)
        path = PathParameters(geof, sat, ane)
        Pfp = fourier_path(f, r, m, path)
        fas = FourierParameters(SourceParameters(100.0), path)
        Pff = fourier_path(f, r, m, fas)
        @test Pfr == Pfp
        @test Pfr == Pff

        # @code_warntype fourier_path(f, r, geof, anef)
        # @code_warntype fourier_path(f, r, geod, aned)
        # @code_warntype fourier_path(f, r, geom, anef)
        # @code_warntype fourier_path(f, r, pathf)
        # @code_warntype fourier_path(f, r, pathd)
        # @code_warntype fourier_path(f, r, pathm)
        # @code_warntype fourier_path(f, r, fasf)
        # @code_warntype fourier_path(f, r, fasd)
        # @code_warntype fourier_path(f, r, fasm)

        Pf = fourier_path(f, r, fasf)
        Pd = fourier_path(f, r, fasd)
        Pm = fourier_path(f, r, fasm)
        @test Pf == Pd.value
        @test Pd == Pm

        f = 10.0
        # @code_warntype fourier_attenuation(f, r, anef, sitef)
        # @code_warntype fourier_attenuation(f, r, aned, sited)
        # @code_warntype fourier_attenuation(f, r, anef, sited)
        # @code_warntype fourier_attenuation(f, r, pathf, sitef)
        # @code_warntype fourier_attenuation(f, r, pathd, sited)
        # @code_warntype fourier_attenuation(f, r, pathm, sited)
        # @code_warntype fourier_attenuation(f, r, fasf)
        # @code_warntype fourier_attenuation(f, r, fasd)
        # @code_warntype fourier_attenuation(f, r, fasm)

        Qf = fourier_attenuation(f, r, fasf)
        Qd = fourier_attenuation(f, r, fasd)
        Qm = fourier_attenuation(f, r, fasm)
        @test Qf == Qd.value
        @test Qd == Qm

        @test fourier_attenuation(-1.0, r, fasf) == 1.0

        # @code_warntype fourier_site(f, sitef)
        # @code_warntype fourier_site(f, sited)
        # @time fourier_site(f, sitef)
        # @time fourier_site(f, sited)

        Sf = fourier_site(f, fasf)
        Sd = fourier_site(f, fasd)
        Sm = fourier_site(f, fasm)
        @test Sf == Sd.value
        @test Sd == Sm

        f = 1.0
        m = 6.0
        r = 10.0
        r_psf = equivalent_point_source_distance(r, m, fasf)
        r_psd = equivalent_point_source_distance(r, m, fasd)
        r_psm = equivalent_point_source_distance(r, m, fasm)

        # @code_warntype fourier_spectral_ordinate(f, m, r_psf, fasf)
        # @code_warntype fourier_spectral_ordinate(f, m, r_psd, fasd)
        # @code_warntype fourier_spectral_ordinate(f, m, r_psm, fasm)

        Af = fourier_spectral_ordinate(f, m, r_psf, fasf)
        Ad = fourier_spectral_ordinate(f, m, r_psd, fasd)
        Am = fourier_spectral_ordinate(f, m, r_psm, fasm)
        @test Af == Ad.value
        @test Ad == Am

        fi = [0.01, 0.1, 1.0, 10.0, 100.0]

        # @code_warntype fourier_spectrum(fi, m, r_psf, fasf)
        # @code_warntype fourier_spectrum(fi, m, r_psf, fasd)
        # @code_warntype fourier_spectrum(fi, m, r_psd, fasd)
        # @code_warntype fourier_spectrum(fi, m, r_psd, fasm)

        Afif = fourier_spectrum(fi, m, r_psf, fasf)
        Afid = fourier_spectrum(fi, m, r_psf, fasd)
        Afim = fourier_spectrum(fi, m, r_psd, fasm)
        for i = 1:length(fi)
            @test Afif[i] == Afid[i].value
        end
        @test all(isapprox.(Afid, Afim))

        fourier_spectrum!(Afif, fi, m, r_psf, fasf)
        fourier_spectrum!(Afid, fi, m, r_psf, fasd)
        fourier_spectrum!(Afim, fi, m, r_psd, fasm)
        for i = 1:length(fi)
            @test Afif[i] == Afid[i].value
        end
        @test all(isapprox.(Afid, Afim))

        StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afif, fi, m, r_psf, fasf)
        StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afid, fi, m, r_psf, fasd)
        StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afim, fi, m, r_psd, fasm)
        for i = 1:length(fi)
            @test Afif[i] == Afid[i].value
        end
        @test all(isapprox.(Afid, Afim))

        ane = AnelasticAttenuationParameters(200.0, 0.0, :Rrup)
        path = PathParameters(geof, sat, ane)
        fas = FourierParameters(SourceParameters(100.0), path)

        m = Dual(6.0)
        r_rup = 10.0
        r_ps = equivalent_point_source_distance(r_rup, m, fas)

        Afid = fourier_spectrum(Vector{Float64}(), m, r_ps, fas)
        @test eltype(Afid) <: Dual
        Afid = fourier_spectrum(fi, m, r_ps, fas)

        sqAfid = StochasticGroundMotionSimulation.squared_fourier_spectrum(fi, m, r_ps, fas)
        @test any(isapprox.(sqAfid, Afid .^ 2))

        # @code_warntype fourier_spectrum!(Afif, fi, m, r_psf, fasf)
        # @code_warntype fourier_spectrum!(Afid, fi, m, r_psf, fasd)
        # @code_warntype fourier_spectrum!(Afim, fi, m, r_psd, fasm)

        Afid = fourier_spectrum(Vector{Float64}(), m, r_ps, fas)
        fourier_spectrum!(Afid, Vector{Float64}(), m, r_ps, fas)
        StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afid, Vector{Float64}(), m, r_ps, fas)

        Afid = fourier_spectrum(fi, m, r_ps, fas)
        fourier_spectrum!(Afid, fi, m, r_ps, fas)
        StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afid, fi, m, r_ps, fas)

        # @code_warntype StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afif, fi, m, r_psf, fasf)
        # @code_warntype StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afid, fi, m, r_psf, fasd)
        # @code_warntype StochasticGroundMotionSimulation.squared_fourier_spectrum!(Afim, fi, m, r_psd, fasm)

        # @code_warntype combined_kappa_frequency(r_psf, fasf)
        # @code_warntype combined_kappa_frequency(r_psd, fasd)
        fkf = combined_kappa_frequency(r_psf, 0.5, fasf)
        fkd = combined_kappa_frequency(r_psd, 0.5, fasd)
        @test fkf == fkd.value

        fkf = combined_kappa_frequency(r_psf, 0.5, fasf)
        fkfd = combined_kappa_frequency(Dual(r_psf), 0.5, fasf)
        @test fkf == fkfd.value

        fkf0 = combined_kappa_frequency(r_psf, 0.5, fas)
        fkf1 = combined_kappa_frequency(r_psf, 0.5, fasf)
        @test fkf0 > fkf1

    end

    @testset "RVT" begin

        @testset "Integration" begin

            Δσf = 100.0
            γ1f = 1.158
            γ2f = 0.5
            Q0f = 212.5
            ηf = 0.65
            κ0f = 0.038

            Δσd = Dual{Float64}(Δσf)
            γ1d = Dual{Float64}(γ1f)
            γ2d = Dual{Float64}(γ2f)
            Q0d = Dual{Float64}(Q0f)
            ηd = Dual{Float64}(ηf)
            κ0d = Dual{Float64}(κ0f)

            srcf = SourceParameters(Δσf)
            srcd = SourceParameters(Δσd)

            Rrefi = [1.0, 50.0, Inf]
            geof = GeometricSpreadingParameters(Rrefi, [γ1f, γ2f])
            geod = GeometricSpreadingParameters(Rrefi, [γ1d, γ2d])
            geom = GeometricSpreadingParameters(Rrefi, [γ1f], [γ2d], BitVector([0, 1]), :Piecewise)

            anef = AnelasticAttenuationParameters(Q0f, ηf)
            aned = AnelasticAttenuationParameters(Q0d, ηd)
            anem = AnelasticAttenuationParameters(Q0f, ηd)

            sat = NearSourceSaturationParameters(:BT15)

            pathf = PathParameters(geof, sat, anef)
            pathd = PathParameters(geod, sat, aned)
            pathm = PathParameters(geom, sat, anem)

            sitef = SiteParameters(κ0f, :Boore2016)
            sited = SiteParameters(κ0d, :Boore2016)

            fasf = FourierParameters(srcf, pathf, sitef)
            fasd = FourierParameters(srcd, pathd, sited)
            fasm = FourierParameters(srcf, pathm, sited)

            m = 5.5
            r = 10.88
            r_psf = equivalent_point_source_distance(r, m, fasf)
            r_psd = equivalent_point_source_distance(r, m, fasd)
            r_psm = equivalent_point_source_distance(r, m, fasm)

            sdof = Oscillator(100.0)

            # Boore comparison (assume his are cgs units)
            ps2db(f) = (1.0 / (2π * sdof.f_n))^2 * 1e4

            dbm0_integrand(f) = StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof) * ps2db(f)
            dbm0ln_integrand(lnf) = exp(lnf) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof) * ps2db(exp(lnf))

            dbm1_integrand(f) = (2π * f) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof) * ps2db(f)
            dbm1ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf)) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof) * ps2db(exp(lnf))

            dbm2_integrand(f) = (2π * f)^2 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof) * ps2db(f)
            dbm2ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf))^2 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof) * ps2db(exp(lnf))

            dbm4_integrand(f) = (2π * f)^4 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof) * ps2db(f)
            dbm4ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf))^4 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof) * ps2db(exp(lnf))

            @time igk = 2 * quadgk(dbm0_integrand, 0.0, Inf)[1]
            @time igk = 2 * quadgk(dbm0_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = 2 * StochasticGroundMotionSimulation.gauss_intervals(dbm0ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)
            @test igk ≈ iglelnm rtol = 1e-4

            @time igk = 2 * quadgk(dbm1_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = 2 * StochasticGroundMotionSimulation.gauss_intervals(dbm1ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)
            @test igk ≈ iglelnm rtol = 1e-4

            @time igk = 2 * quadgk(dbm2_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = 2 * StochasticGroundMotionSimulation.gauss_intervals(dbm2ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)
            @test igk ≈ iglelnm rtol = 1e-3

            @time igk = 2 * quadgk(dbm4_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = 2 * StochasticGroundMotionSimulation.gauss_intervals(dbm4ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)
            @test igk ≈ iglelnm rtol = 1e-3

            m0_integrand(f) = StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof)
            m0ln_integrand(lnf) = exp(lnf) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof)

            m1_integrand(f) = (2π * f) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof)
            m1ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf)) * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof)

            m2_integrand(f) = (2π * f)^2 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof)
            m2ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf))^2 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof)

            m4_integrand(f) = (2π * f)^4 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(f, m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(f, sdof)
            m4ln_integrand(lnf) = exp(lnf) * (2π * exp(lnf))^4 * StochasticGroundMotionSimulation.squared_fourier_spectral_ordinate(exp(lnf), m, r_psf, fasf) * StochasticGroundMotionSimulation.squared_transfer(exp(lnf), sdof)

            @time igk = quadgk(m0_integrand, exp(-7.0), exp(7.0))[1]
            @time igle = StochasticGroundMotionSimulation.gauss_interval(m0_integrand, 2000, 0.0, 300.0)
            @time igleln = StochasticGroundMotionSimulation.gauss_interval(m0ln_integrand, 750, -7.0, 7.0)

            @time iglelnm = StochasticGroundMotionSimulation.gauss_intervals(m0ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)

            # @test igk ≈ igle rtol=1e-2
            @test igk ≈ igleln rtol = 1e-4
            @test igk ≈ iglelnm rtol = 1e-4

            lnfi = log.([1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, sdof.f_n])
            sort!(lnfi)

            @time igk = quadgk(m1_integrand, exp(-7.0), exp(7.0))[1]
            @time igle = StochasticGroundMotionSimulation.gauss_interval(m1_integrand, 1500, 0.0, 300.0)
            @time igleln = StochasticGroundMotionSimulation.gauss_interval(m1ln_integrand, 750, -7.0, 7.0)
            @time iglelnm = StochasticGroundMotionSimulation.gauss_intervals(m1ln_integrand, 250, -7.0, log(sdof.f_n), 7.0)

            @time iglelnm = StochasticGroundMotionSimulation.gauss_intervals(m1ln_integrand, 30, lnfi...)

            @time itr = StochasticGroundMotionSimulation.trapezoidal(m1ln_integrand, 60, lnfi...)

            # @test igk ≈ igle rtol=1e-2
            @test igk ≈ iglelnm rtol = 1e-4
            @test igk ≈ itr rtol = 1e-3

            @time igk = quadgk(m2_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = StochasticGroundMotionSimulation.gauss_intervals(m2ln_integrand, 30, lnfi...)
            @time itr = StochasticGroundMotionSimulation.trapezoidal(m2ln_integrand, 60, lnfi...)

            @test igk ≈ iglelnm rtol = 1e-4
            @test igk ≈ itr rtol = 1e-3

            @time igk = quadgk(m4_integrand, exp(-7.0), exp(7.0))[1]
            @time iglelnm = StochasticGroundMotionSimulation.gauss_intervals(m4ln_integrand, 30, lnfi...)
            @time itr = StochasticGroundMotionSimulation.trapezoidal(m4ln_integrand, 60, lnfi...)

            @test igk ≈ iglelnm rtol = 1e-4
            @test igk ≈ itr rtol = 1e-3

            integrand(x) = sin(x)

            intervals = 101
            x_min = 0.0
            x_max = 2π
            xx = collect(range(x_min, stop = x_max, length = intervals))
            yy = integrand.(xx)

            isr = StochasticGroundMotionSimulation.simpsons_rule(xx, yy)
            ist = StochasticGroundMotionSimulation.trapezoidal_rule(xx, yy)

            igk = quadgk(integrand, x_min, x_max)[1]

            @test isr ≈ igk atol = 10 * eps()
            @test ist ≈ igk atol = 10 * eps()

            m = 6.0
            r = 100.0
            src = SourceParameters(100.0)
            geo = GeometricSpreadingParameters([1.0, 50.0, Inf], [1.0, 0.5])
            sat = NearSourceSaturationParameters(:BT15)
            ane = AnelasticAttenuationParameters(200.0, 0.4, :Rrup)
            path = PathParameters(geo, sat, ane)
            site = SiteParameters(0.039)
            fas = FourierParameters(src, path, site)
            sdof = Oscillator(1.0)

            integrand(f) = StochasticGroundMotionSimulation.squared_transfer(f, sdof) * fourier_spectral_ordinate(f, m, r, fas)^2

            intervals = 101
            x_min = sdof.f_n / 1.1
            x_max = sdof.f_n * 1.1
            xx = collect(range(x_min, stop = x_max, length = intervals))
            yy = integrand.(xx)

            isr = StochasticGroundMotionSimulation.simpsons_rule(xx, yy)
            igk = quadgk(integrand, x_min, x_max)[1]

            @test isr ≈ igk rtol = 1e-6
            @test isr ≈ igk atol = 1e-6

            integrand(f) = StochasticGroundMotionSimulation.squared_transfer(f, sdof) * fourier_spectral_ordinate(f, m, r, fas)^2

            intervals = 101
            x_min = 100.0
            x_max = 200.0
            xx = collect(range(x_min, stop = x_max, length = intervals))
            yy = integrand.(xx)

            isr = StochasticGroundMotionSimulation.simpsons_rule(xx, yy)
            igk = quadgk(integrand, x_min, x_max)[1]

            @test isr ≈ igk rtol = 1e-3
            @test isr ≈ igk atol = 1e-6

            integrand(f) = (2π * f)^4 * StochasticGroundMotionSimulation.squared_transfer(f, sdof) * fourier_spectral_ordinate(f, m, r, fas)^2

            intervals = 101
            x_min = 100.0
            x_max = 200.0
            xx = collect(range(x_min, stop = x_max, length = intervals))
            yy = integrand.(xx)

            isr = StochasticGroundMotionSimulation.simpsons_rule(xx, yy)
            igk = quadgk(integrand, x_min, x_max)[1]

            @test isr ≈ igk rtol = 1e-3
            @test isr ≈ igk atol = 1e-6

            x_min = 300.0
            x_max = 500.0
            xx = collect(range(x_min, stop = x_max, length = intervals))
            yy = integrand.(xx)

            isr = StochasticGroundMotionSimulation.simpsons_rule(xx, yy)
            igk = quadgk(integrand, x_min, x_max)[1]

            @test isr ≈ igk rtol = 1e-3
            @test isr ≈ igk atol = 1e-6

        end

        @testset "Spectral Moments" begin

            Δσf = 100.0
            γ1f = 1.0
            γ2f = 0.5
            Q0f = 200.0
            ηf = 0.4
            κ0f = 0.039

            Δσd = Dual{Float64}(Δσf)
            γ1d = Dual{Float64}(γ1f)
            γ2d = Dual{Float64}(γ2f)
            Q0d = Dual{Float64}(Q0f)
            ηd = Dual{Float64}(ηf)
            κ0d = Dual{Float64}(κ0f)

            srcf = SourceParameters(Δσf)
            srcd = SourceParameters(Δσd)

            Rrefi = [1.0, 50.0, Inf]
            geof = GeometricSpreadingParameters(Rrefi, [γ1f, γ2f])
            geod = GeometricSpreadingParameters(Rrefi, [γ1d, γ2d])
            geom = GeometricSpreadingParameters(Rrefi, [γ1f], [γ2d], BitVector([0, 1]), :Piecewise)

            anef = AnelasticAttenuationParameters(Q0f, ηf)
            aned = AnelasticAttenuationParameters(Q0d, ηd)
            anem = AnelasticAttenuationParameters(Q0f, ηd)

            sat = NearSourceSaturationParameters(:BT15)

            pathf = PathParameters(geof, sat, anef)
            pathd = PathParameters(geod, sat, aned)
            pathm = PathParameters(geom, sat, anem)

            sitef = SiteParameters(κ0f)
            sited = SiteParameters(κ0d)

            fasf = FourierParameters(srcf, pathf, sitef)
            fasd = FourierParameters(srcd, pathd, sited)
            fasm = FourierParameters(srcf, pathm, sited)

            m = 6.0
            r = 10.0
            r_psf = equivalent_point_source_distance(r, m, fasf)
            r_psd = equivalent_point_source_distance(r, m, fasd)
            r_psm = equivalent_point_source_distance(r, m, fasm)

            sdof = Oscillator(0.1)

            order = 0
            m0f = spectral_moment(order, m, r_psf, fasf, sdof)
            m0d = spectral_moment(order, m, r_psd, fasd, sdof)
            m0m = spectral_moment(order, m, r_psm, fasm, sdof)
            @test m0f ≈ m0d.value
            @test m0d ≈ m0m

            # @code_warntype spectral_moment(order, m, r_psf, fasf, sdof)
            # @code_warntype spectral_moment(order, m, r_psd, fasd, sdof)
            # @code_warntype spectral_moment(order, m, r_psm, fasm, sdof)

            # @time spectral_moments([0, 1, 2, 4], m, r_psf, fasf, sdof)
            # @time spectral_moments([0, 1, 2, 4], m, r_psd, fasd, sdof)
            # @time spectral_moments([0, 1, 2, 4], m, r_psm, fasm, sdof)

            # @code_warntype spectral_moments([0, 1, 2, 4], m, r_psf, fasf, sdof)
            # @code_warntype spectral_moments([0, 1, 2, 4], m, r_psd, fasd, sdof)
            # @code_warntype spectral_moments([0, 1, 2, 4], m, r_psm, fasm, sdof)

            smi = spectral_moments([0, 1, 2, 4], m, r_psf, fasf, sdof)
            smigk = StochasticGroundMotionSimulation.spectral_moments_gk([0, 1, 2, 4], m, r_psf, fasf, sdof)

            @test all(isapprox.(smi, smigk, rtol = 1e-3))

            sdof = Oscillator(1 / 3)
            smi = spectral_moments([0, 1, 2, 4], m, r_psf, fasf, sdof)
            smigk = StochasticGroundMotionSimulation.spectral_moments_gk([0, 1, 2, 4], m, r_psf, fasf, sdof)

            @test all(isapprox.(smi, smigk, rtol = 1e-3))

            m = Dual(6.0)
            smdi = spectral_moments([0, 1, 2, 4], m, r_psf, fasf, sdof)
            smfi = spectral_moments([0, 1, 2, 4], m.value, r_psf, fasf, sdof)
            @test any(isapprox.(smdi, smfi))

            m = Dual(6.0)
            smdi = spectral_moment(0, m, r_psf, fasf, sdof)
            smfi = spectral_moment(0, m.value, r_psf, fasf, sdof)
            @test smdi ≈ smfi

            rvt = RandomVibrationParameters()
            Ti = [0.01, 0.1, 1.0]
            m = Dual(6.0)
            Sadi = rvt_response_spectrum(Ti, m, 10.0, fasf, rvt)
            Safi = rvt_response_spectrum(Ti, m.value, 10.0, fasf, rvt)
            for i = 1:length(Ti)
                @test Sadi[i].value ≈ Safi[i]
            end

            rvt_response_spectrum!(Sadi, Ti, m, 10.0, fasf, rvt)
            rvt_response_spectrum!(Safi, Ti, m.value, 10.0, fasf, rvt)
            for i = 1:length(Ti)
                @test Sadi[i].value ≈ Safi[i]
            end

        end

        @testset "Peak Factor" begin

            Δσf = 100.0
            γ1f = 1.0
            γ2f = 0.5
            Q0f = 200.0
            ηf = 0.4
            κ0f = 0.039

            Δσd = Dual{Float64}(Δσf)
            γ1d = Dual{Float64}(γ1f)
            γ2d = Dual{Float64}(γ2f)
            Q0d = Dual{Float64}(Q0f)
            ηd = Dual{Float64}(ηf)
            κ0d = Dual{Float64}(κ0f)

            srcf = SourceParameters(Δσf)
            srcd = SourceParameters(Δσd)

            Rrefi = [1.0, 50.0, Inf]
            geof = GeometricSpreadingParameters(Rrefi, [γ1f, γ2f])
            geod = GeometricSpreadingParameters(Rrefi, [γ1d, γ2d])
            geom = GeometricSpreadingParameters(Rrefi, [γ1f], [γ2d], BitVector([0, 1]), :Piecewise)

            anef = AnelasticAttenuationParameters(Q0f, ηf)
            aned = AnelasticAttenuationParameters(Q0d, ηd)
            anem = AnelasticAttenuationParameters(Q0f, ηd)

            sat = NearSourceSaturationParameters(:BT15)

            pathf = PathParameters(geof, sat, anef)
            pathd = PathParameters(geod, sat, aned)
            pathm = PathParameters(geom, sat, anem)

            sitef = SiteParameters(κ0f)
            sited = SiteParameters(κ0d)

            fasf = FourierParameters(srcf, pathf, sitef)
            fasd = FourierParameters(srcd, pathd, sited)
            fasm = FourierParameters(srcf, pathm, sited)

            m = 7.0
            r = 1.0

            r_psf = equivalent_point_source_distance(r, m, fasf)
            r_psd = equivalent_point_source_distance(r, m, fasd)
            r_psm = equivalent_point_source_distance(r, m, fasm)

            sdof = Oscillator(1.0)

            @time pfps = StochasticGroundMotionSimulation.peak_factor_dk80(m, r_psf, fasf, sdof)
            @time pfpsn = StochasticGroundMotionSimulation.peak_factor_dk80(m, r_psf, fasf, sdof, nodes = 30)
            @time pfgk = StochasticGroundMotionSimulation.peak_factor_dk80_gk(m, r_psf, fasf, sdof)

            @test pfps ≈ pfgk rtol = 1e-6
            @test pfpsn ≈ pfgk rtol = 1e-6

            @time pfps = StochasticGroundMotionSimulation.peak_factor_cl56(m, r_psf, fasf, sdof)
            @time pfpsn = StochasticGroundMotionSimulation.peak_factor_cl56(m, r_psf, fasf, sdof, nodes = 40)
            @time pfgk = StochasticGroundMotionSimulation.peak_factor_cl56_gk(m, r_psf, fasf, sdof)

            @test pfps ≈ pfgk rtol = 1e-5
            @test pfpsn ≈ pfgk rtol = 1e-5

            @test StochasticGroundMotionSimulation.vanmarcke_cdf(-1.0, 10.0, 0.5) == 0.0

            rvt = RandomVibrationParameters(:CL56)
            pf = peak_factor(6.0, 10.0, fasf, sdof, rvt)
            rvt = RandomVibrationParameters(:DK80)
            pf = peak_factor(6.0, 10.0, fasf, sdof, rvt)
            rvt = RandomVibrationParameters(:PS)
            pf = peak_factor(6.0, 10.0, fasf, sdof, rvt)
            @test isnan(pf)

            m = 6.0
            r_rup = 10.0
            r_ps = equivalent_point_source_distance(r_rup, m, fasf)
            Dex = excitation_duration(m, r_ps, fasf, rvt)
            m0 = spectral_moment(0, m, r_ps, fasf, sdof)
            rvt = RandomVibrationParameters(:PS)
            pf = peak_factor(m, r_ps, Dex, m0, fasf, sdof, rvt)
            @test isnan(pf)
            pf = peak_factor(Dual(m), r_ps, Dex, m0, fasf, sdof, rvt)
            @test isnan(pf)
            pf = peak_factor(m, r_ps, Dex, Dual(m0), fasf, sdof, rvt)
            @test isnan(pf)
            rvt = RandomVibrationParameters(:CL56)
            pf = peak_factor(m, r_ps, Dex, m0, fasf, sdof, rvt)
            @test isnan(pf) == false
            @test pf ≈ peak_factor(m, r_ps, fasf, sdof, RandomVibrationParameters(:CL56))

            pfi = StochasticGroundMotionSimulation.peak_factor_integrand_cl56(0.0, 10.0, 10.0)
            @test pfi ≈ 1.0
            pfi = StochasticGroundMotionSimulation.peak_factor_integrand_cl56(Inf, 10.0, 10.0)
            @test pfi ≈ 0.0

            pf0 = StochasticGroundMotionSimulation.peak_factor_cl56(10.0, 10.0)
            pf1 = StochasticGroundMotionSimulation.peak_factor_cl56(10.0, 10.0, nodes = 50)
            @test pf0 ≈ pf1

        end

        @testset "Response Spectra" begin

            Δσf = 100.0
            γ1f = 1.0
            γ2f = 0.5
            Q0f = 200.0
            ηf = 0.4
            κ0f = 0.039

            Δσd = Dual{Float64}(Δσf)
            γ1d = Dual{Float64}(γ1f)
            γ2d = Dual{Float64}(γ2f)
            Q0d = Dual{Float64}(Q0f)
            ηd = Dual{Float64}(ηf)
            κ0d = Dual{Float64}(κ0f)

            srcf = SourceParameters(Δσf)
            srcd = SourceParameters(Δσd)

            Rrefi = [1.0, 50.0, Inf]
            geof = GeometricSpreadingParameters(Rrefi, [γ1f, γ2f])
            geod = GeometricSpreadingParameters(Rrefi, [γ1d, γ2d])
            geom = GeometricSpreadingParameters(Rrefi, [γ1f], [γ2d], BitVector([0, 1]), :Piecewise)

            anef = AnelasticAttenuationParameters(Q0f, ηf)
            aned = AnelasticAttenuationParameters(Q0d, ηd)
            anem = AnelasticAttenuationParameters(Q0f, ηd)

            sat = NearSourceSaturationParameters(:BT15)

            pathf = PathParameters(geof, sat, anef)
            pathd = PathParameters(geod, sat, aned)
            pathm = PathParameters(geom, sat, anem)

            sitef = SiteParameters(κ0f)
            sited = SiteParameters(κ0d)

            fasf = FourierParameters(srcf, pathf, sitef)
            fasd = FourierParameters(srcd, pathd, sited)
            fasm = FourierParameters(srcf, pathm, sited)

            m = 7.0
            r = 1.0

            r_psf = equivalent_point_source_distance(r, m, fasf)
            r_psd = equivalent_point_source_distance(r, m, fasd)
            r_psm = equivalent_point_source_distance(r, m, fasm)

            Ti = [0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 5.0, 7.5, 10.0]

            rvt = RandomVibrationParameters()

            Saif = rvt_response_spectrum(Ti, m, r_psf, fasf, rvt)
            Said = rvt_response_spectrum(Ti, m, r_psd, fasd, rvt)
            Saim = rvt_response_spectrum(Ti, m, r_psm, fasm, rvt)
            for i = 1:length(Ti)
                @test Saif[i] ≈ Said[i].value
            end
            @test all(isapprox.(Said, Saim))

            # @code_warntype rvt_response_spectrum(Ti, m, r_psf, fasf, rvt)
            # @code_warntype rvt_response_spectrum(Ti, m, r_psd, fasd, rvt)
            # @code_warntype rvt_response_spectrum(Ti, m, r_psm, fasm, rvt)

        end
    end
end


text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractThreadLocal{T} end

Abstract type for thread-local values of type `T`.

The value for the current thread is accessed via
`getindex(::AbstractThreadLocal)` and `setindex(::AbstractThreadLocal, x).

To access both regular and thread-local values in a unified manner, use
the function [`getlocalvalue`](@ref).

To get the all values across all threads, use the function
[`getallvalues`](@ref).

Default implementation is [`ThreadLocal`](@ref).
"""
abstract type AbstractThreadLocal{T} end
export AbstractThreadLocal

@deprecate isdefined_local(x) isassigned(x)

"""
    getlocalvalue(x::Any) = x
    getlocalvalue(x::ThreadLocal) = x[]

Access plain values and thread-local values in a unified fashion.
"""
function getlocalvalue end
export getlocalvalue

"""
    getallvalues(v::AbstractThreadLocal{T})::AbstractVector{T}

Access the all values (one for each thread) of a thread-local value as a
vector. Can only be called in single-threaded code sections.
"""
function getallvalues end
export getallvalues

function _protect_from_resize(V::Vector)
    reshape(V, length(V), 1)
    V
end

"""
    ThreadLocal{T} <: AbstractThreadLocal{T}

Represents a thread-local value. See [`AbstractThreadLocal`](@ref) for
the API.

Constructors:

```julia
ThreadLocal{T}() where {T}
ThreadLocal(value::T) where {T}
ThreadLocal{T}(f::Base.Callable) where {T}
```

Examples:

```julia
tlvalue = ThreadLocal(0)
@onthreads allthreads() tlvalue[] = Base.Threads.threadid()
getallvalues(tlvalue) == allthreads()
```

```julia
rand_value_on_each_thread = ThreadLocal{Float64}(rand)
all(x -> 0 < x < 1, getallvalues(rand_value_on_each_thread))
```
"""
struct ThreadLocal{T} <: AbstractThreadLocal{T}
    value::Vector{T}

    ThreadLocal{T}(::UndefInitializer) where {T} =
        new{T}(_protect_from_resize(Vector{T}(undef, nthreads())))

    ThreadLocal{T}(value::T) where {T} =
        new{T}(_protect_from_resize([deepcopy(value) for i in 1:nthreads()]))
end

export ThreadLocal

ThreadLocal{T}() where {T} = ThreadLocal{T}(T)

ThreadLocal(value::T) where {T} = ThreadLocal{T}(value)

@inline Base.getindex(x::ThreadLocal) = x.value[threadid()]

@inline Base.setindex!(x::ThreadLocal, y) = x.value[threadid()] = y

Base.eltype(x::Type{<:ThreadLocal{T}}) where {T} = @isdefined(T) ? T : Any

@inline Base.get(x::ThreadLocal) = x[]

@inline Base.get(default::Base.Callable, x::ThreadLocal{T}) where {T} = isassigned(x) ? x[] : convert(T, default())

@inline Base.get(x::ThreadLocal, default) = get(() -> default, x)

function Base.get!(default::Base.Callable, x::ThreadLocal{T}) where {T}
    if isassigned(x)
        x[]
    else
        x[] = convert(T, default())
    end
end

Base.get!(x::ThreadLocal, default) = get!(() -> default, x)

@inline getlocalvalue(x) = x
@inline getlocalvalue(x::ThreadLocal) = x[]

Base.isassigned(x::ThreadLocal) = isassigned(x.value, threadid())

function getallvalues(x::ThreadLocal)
    @static if VERSION >= v"1.3.0-alpha.0"
        x.value
    else
        if !Base.Threads.in_threaded_loop[]
            x.value
        else
            throw(InvalidStateException("Can not access thread local values across threads in multi-threaded code sections"))
        end
    end
end

ext: jl


text:
module Flags

export Flag, FlagParser, BoolFlag, parse, formathelp, printhelp

struct Flag
    name::String
    default::Any
    usage::String
end

Flag(
    name::String,
    default=nothing,
    usage=""
) = Flag(name, default, usage)

function BoolFlag(name::String, usage::String)::Flag
    Flag(name, false, usage)
end

function formathelp(flag::Flag)::String
    if length(flag.name) == 1
        f = "-$(flag.name)"
    else
        f = "--$(flag.name)"
    end

    if flag.default === nothing
        d = ""
    else
        d = "=$(flag.default)"
    end

    """  $f$d: $(flag.usage)"""
end

struct FlagParser
    name::String
    flags::Vector{Flag}

    FlagParser(name="", flags=Vector{Flag}()) = new(name, flags)
end

function printhelp(parser::FlagParser)
    if length(parser.name) > 0
        name = parser.name
    else
        name = PROGRAM_FILE
    end

    println("Usage of $name:")
    for flag in parser.flags
        println(formathelp(flag))
    end
end

function parse(parser::FlagParser, args::Vector{String})

end

end # module


text:
abstract type AbstractGPPrior{K<:Kernel,Tmean<:PriorMean} end# <: AbstractGPs.AbstractGPModel end

kernel(gp::AbstractGPPrior) = gp.kernel
setkernel!(gp::AbstractGPPrior, k::Kernel) = gp.kernel = k
Distributions.mean(gp::AbstractGPPrior) = gp.μ₀
Distributions.mean(gp::AbstractGPPrior, X::AbstractVector) = gp.μ₀(X)
setmean!(gp::AbstractGPPrior, μ₀::PriorMean) = gp.μ₀ = μ₀
Distributions.cov(gp::AbstractGPPrior) = gp.K

mutable struct GPPrior{K<:Kernel,Tmean<:PriorMean} <: AbstractGPPrior{K,Tmean}
    kernel::K
    μ₀::Tmean
end

mutable struct TPrior{K<:Kernel,Tmean<:PriorMean,T} <: AbstractGPPrior{K,Tmean}
    kernel::K
    μ₀::Tmean
    ν::T # Number of degrees of freedom
    l²::T # Expectation of ||L^{-1}(f-μ⁰)||₂² (L2 norm)
    χ::T  # Expectation of σ
end


text:
function print_startup_message()

    println("")

    s = """

      ███████ ███████ ██████
      ██      ██      ██
      ███████ ███████ ██████
           ██      ██ ██
      ███████ ███████ ██

      """
    println(s)
end

function print_output(training_setup::Traindata_Settings,num_entries, num_trouble)
    print_startup_message()
    header = (["Setting", "value"])
    data =  [ "number meshes"                length(training_setup.meshes) ;
              "different polydegrees"        length(training_setup.polydegrees);
              "# Entries in Traindataset"    num_entries;
              "good Cells"                   num_entries - num_trouble;
              "Troubled Cells"               num_trouble;];
    pretty_table(data; header = header)
end


text:
#=

  Euler #44 in Julia.

  Problem 44:
  """
  Pentagonal numbers are generated by the formula, P(n)=n(3n−1)/2.
  The first ten pentagonal numbers are:

  1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...

  It can be seen that P(4) + P(7) = 22 + 70 = 92 = P(8). However,
  their difference,  70 − 22 = 48, is not pentagonal.

  Find the pair of pentagonal numbers, P(j) and P(k), for which their sum
  and difference is pentagonal and D = |P(k) − P(j)| is minimised; what
  is the value of D?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.04571280s
function euler44a()
    s = (1:2500).|>n->floor(Int,n*(3*n-1) / 2)
    t = Set(s)
    d = 10_000_000
    for j in s
        for k in s
            a = j+k
            b = abs(j-k)
            if j < k &&
                a < d &&
                a in t &&
                b < d &&
                b in t

                d = b
                break
            end
        end
    end

    return d
end

# 0.04486643s
function euler44b()
    n = 2500
    s = (1:2500).|>n->floor(Int,n*(3*n-1) / 2)
    t = Set(s)
    d = 10_000_000

    for j in 1:n
        for k = 1:n
            if s[j] >= s[k]
                continue;
            end
            a = s[j]+s[k]
            if a >= d
                continue
            end
            b = abs(s[j]-s[k])
            if b >= d
                continue
            end
            if a in t && b in t
                d = b
                break
            end
        end
    end

    return d
end

# run_euler(euler44a)
run_euler(euler44b)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"""
    PlanetParameters

Module containing physical constants and parameters characterizing the planet.
"""
module PlanetParameters
using ..ParametersType

# Physical constants
@exportparameter gas_constant      8.3144598     "Universal gas constant (J/mol/K)"
@exportparameter light_speed       2.99792458e8  "Speed of light in vacuum (m/s)"
@exportparameter h_Planck          6.626e-34     "Planck constant (m^2 kg/s)"
@exportparameter k_Boltzmann       1.381e-23     "Boltzmann constant (m^2 kg/s^2/K)"
@exportparameter Stefan            5.670e-8      "Stefan-Boltzmann constant (W/m^2/K^4)"
@exportparameter astro_unit        1.4959787e11  "Astronomical unit (m)"
@exportparameter k_Karman          0.4           "Von Karman constant (1)"

# Properties of dry air
@exportparameter molmass_dryair    28.97e-3      "Molecular weight dry air (kg/mol)"
@exportparameter R_d               gas_constant/
                                   molmass_dryair "Gas constant dry air (J/kg/K)"
@exportparameter kappa_d           2//7          "Adiabatic exponent dry air"
@exportparameter cp_d              R_d/kappa_d   "Isobaric specific heat dry air"
@exportparameter cv_d              cp_d - R_d    "Isochoric specific heat dry air"

# Properties of water
@exportparameter ρ_cloud_liq       1e3           "Density of liquid water (kg/m^3)"
@exportparameter ρ_cloud_ice       916.7         "Density of ice water (kg/m^3)"
@exportparameter molmass_water     18.01528e-3   "Molecular weight (kg/mol)"
@exportparameter molmass_ratio     molmass_dryair/
                                   molmass_water "Molar mass ratio dry air/water"
@exportparameter R_v               gas_constant/
                                   molmass_water "Gas constant water vapor (J/kg/K)"
@exportparameter cp_v              1859          "Isobaric specific heat vapor (J/kg/K)"
@exportparameter cp_l              4181          "Isobaric specific heat liquid (J/kg/K)"
@exportparameter cp_i              2100          "Isobaric specific heat ice (J/kg/K)"
@exportparameter cv_v              cp_v - R_v    "Isochoric specific heat vapor (J/kg/K)"
@exportparameter cv_l              cp_l          "Isochoric specific heat liquid (J/kg/K)"
@exportparameter cv_i              cp_i          "Isochoric specific heat ice (J/kg/K)"
@exportparameter T_freeze          273.15        "Freezing point temperature (K)"
@exportparameter T_min             150.0         "Minimum temperature guess in saturation adjustment (K)"
@exportparameter T_max             1000.0        "Maximum temperature guess in saturation adjustment (K)"
@exportparameter T_icenuc          233.00        "Homogeneous nucleation temperature (K)"
@exportparameter T_triple          273.16        "Triple point temperature (K)"
@exportparameter T_0               T_triple      "Reference temperature (K)"
@exportparameter LH_v0             2.5008e6      "Latent heat vaporization at T_0 (J/kg)"
@exportparameter LH_s0             2.8344e6      "Latent heat sublimation at T_0 (J/kg)"
@exportparameter LH_f0             LH_s0 - LH_v0  "Latent heat of fusion at T_0 (J/kg)"
@exportparameter e_int_v0          LH_v0 - R_v*T_0 "Specific internal energy of
                                                    vapor at T_0 (J/kg)"
@exportparameter e_int_i0          LH_f0         "Specific internal energy of
                                                  ice at T_0 (J/kg)"
@exportparameter press_triple      611.657       "Triple point vapor pressure (Pa)"

# Properties of sea water
@exportparameter ρ_ocean           1.035e3       "Reference density sea water (kg/m^3)"
@exportparameter cp_ocean          3989.25       "Specific heat sea water (J/kg/K)"

# Planetary parameters
@exportparameter planet_radius     6.371e6       "Mean planetary radius (m)"
@exportparameter day               86400         "Length of day (s)"
@exportparameter Omega             7.2921159e-5  "Ang. velocity planetary rotation (1/s)"
@exportparameter grav              9.81          "Gravitational acceleration (m/s^2)"
@exportparameter year_anom         365.26*day    "Length of anomalistic year (s)"
@exportparameter orbit_semimaj     1*astro_unit  "Length of semimajor orbital axis (m)"
@exportparameter TSI               1362          "Total solar irradiance (W/m^2)"
@exportparameter MSLP              1.01325e5     "Mean sea level pressure (Pa)"

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# de Kort, S. R., Dickinson, A., and Clayton, N. S. (2005), Retrospective
# cognition by food-caching western scrub-jays, Learning and Motivation
# 36:159-176

function _summarize(::Experiment{:deKort05}, results)
    datainspectpeanut = @where(results, :action .== "inspect", :foodtype .== "peanut")
    datainspectwaxworm = @where(results, :action .== "inspect", :foodtype .== "other")
    data = @transform datainspectwaxworm counts = datainspectwaxworm.counts ./
                                (1e-16 .+ datainspectpeanut.counts .+ datainspectwaxworm.counts)
    combine(groupby(data, [:group, :action, :set, :RI]),
       d -> DataFrame(μ = mean(d.counts), sem = sem(d.counts),
                      firstinspection = mean(d.firstinspection),
                      n = length(d.counts)))
end
function statistical_tests(exp::Experiment{:deKort05}, data)
    tests = Test[]
    id = "inspect"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :action .== "inspect"),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = [:RI, :set, :foodtype],
                                     id = :id, data = :counts),
                      exp.tests))
    d = @where(data, :action .== "inspect", :set .== 8)
    sort!(d, [:id, :foodtype, :RI])
    push!(tests, Test("first inspect trial 8",
                      chisqtest(firstinspect_dist(@where(d, :group .== "degrade")),
                                firstinspect_dist(@where(d, :group .== "ripen"))),
                      exp.tests))
    TestCollection(tests)
end
function run!(::Experiment{:deKort05}, models)
    results = DataFrame(group = Union{String,Missing}[],
                        foodtype = Union{String,Missing}[],
                        id = Int[], set = Int[], trial = Union{String,Missing}[],
                        RI = Union{Int,Missing}[],
                        action = String[], counts = Int[],
                        firstinspection = Union{Missing,Int}[])
    id = 0
    for group in ("degrade", "ripen")
        for i in 1:(group == "degrade" ? 5 : 7)
            id += 1
            m = models[id]
            for set in 1:8
                for RI in (4u"hr", 28u"hr")
                    Clayton01_trial!(m, results, group, id, set, Peanut,
                                     Waxworm, RI,
                                     degrade = isdegrade(group, Waxworm, RI),
                                     pilfer = false)
                end
            end
        end
    end
    results
end


text:
using DelimitedFiles
using FFTW

function derive(data, dim, dx)
    out = Array{Complex,2}(undef,size(data,1), size(data,2))
    dim2 = dim+1
    if dim2 == 3
        dim2 = 1
    end
    for i = 1:size(data,dim)
        for j = 1:size(data,dim2)
            #println(j,'\t',i)
            if (dim == 1)
                if (i == size(data,dim))
                    out[i,j] = (data[1, j] - data[i, j])/dx
                else
                    out[i,j] = (data[i + 1, j] - data[i, j])/dx
                end
            else
                if (i == size(data,dim))
                    out[j,i] = (data[j,1] - data[j,i])/dx
                else
                    out[j,i] = (data[j,i + 1] - data[j,i])/dx
                end
            end
        end
    end
    return out
end

# We are calculating the energy to check <Psi|H|Psi>
function calculate_energy(wfc, H_k, H_r, Ax, Ay, g, xDim, yDim, dx, dy)
    hbar = 1.05457148e-34
    omega = 0.6
    omegaX = 6.283

    # Creating momentum and conjugate wavefunctions
    wfc_k = fft(wfc)
    wfc_c = conj(wfc)

    # Finding the momentum and real-space energy terms
    energy_k = wfc_c.*ifft((H_k) .* wfc_k)
    energy_r = wfc_c.* (H_r).* wfc
    energy_i = wfc_c.* (0.5*g*abs2.(wfc)).* wfc

    energy_l = wfc_c.*(im*hbar*(Ax.*derive(wfc,1,dx) + Ay.*derive(wfc,2,dy)))

    # Integrating over all space
    energy_if = 0
    energy_lf = 0
    energy_kf = 0
    energy_rf = 0
    for i = 1:xDim*yDim
        energy_if += real(energy_i[i])*dx*dy
        energy_rf += real(energy_r[i])*dx*dy
        energy_lf += real(energy_l[i])*dx*dy
        energy_kf += real(energy_k[i])*dx*dy
    end

    println("Kinetic energy:", "\t\t\t", energy_kf)
    println("Potential energy:", "\t\t", energy_rf)
    println("Internal energy:", "\t\t", energy_if)
    println("Angular Momentum energy:", '\t', energy_lf)
    println("Total energy:", "\t\t\t", energy_kf+energy_if+energy_rf+energy_lf)

end

# Read in param.cfg file
function calculate(param_file::String, data_dir::String)
    parameters = Dict()

    for line in readlines(data_dir*param_file)
        if line != "[Params]"
            tmp = split(line,"=")
            parameters[tmp[1]] = tmp[2]
        end
    end

    xDim = parse(Int64, parameters["xDim"])
    yDim = parse(Int64, parameters["yDim"])
    dx = parse(Float64, parameters["dx"])
    dy = parse(Float64, parameters["dy"])

    omega = parse(Float64, parameters["omega"])
    omegaX = parse(Float64, parameters["omegaX"])

    g = parse(Float64, parameters["gDenConst"])

    Ax = readdlm(data_dir*"Ax_0")
    Ay = readdlm(data_dir*"Ay_0")
    K = readdlm(data_dir*"K_0")
    V = readdlm(data_dir*"V_0")

    Ax = reshape(Ax, xDim, yDim)
    Ay = reshape(Ay, xDim, yDim)
    K = reshape(K, xDim, yDim)
    V = reshape(V, xDim, yDim)

    start = 0
    ende = parse(Int64, parameters["esteps"])
    endg = parse(Int64, parameters["gsteps"])
    incr = parse(Int64, parameters["printSteps"])

    # Ground State Evolution
    println("Starting imaginary time energy calculation")
    for i = start:incr:endg
        wfc = readdlm(data_dir*"wfc_0_const_"*string(i)) +
              readdlm(data_dir*"wfc_0_consti_"*string(i))*im
        println(data_dir*"wfc_0_const_"*string(i), '\t',
                data_dir*"wfc_0_consti_"*string(i))
        wfc = reshape(wfc, xDim, yDim)
        calculate_energy(wfc, K, V, Ax, Ay, g, xDim, yDim, dx, dy)
        println()
    end

    println()

    # Ground State Evolution
    println("Starting real time energy calculation")
    for i = start:incr:ende
        wfc = readdlm(data_dir*"wfc_ev_"*string(i)) +
              readdlm(data_dir*"wfc_evi_"*string(i))*im
        println(data_dir*"wfc_0_const_"*string(i), '\t',
                data_dir*"wfc_0_consti_"*string(i))
        wfc = reshape(wfc, xDim, yDim)
        calculate_energy(wfc, K, V, Ax, Ay, g, xDim, yDim, dx, dy)
        println()
    end

end

calculate("Params.dat", "../data/")

max_stars_repo_name: hifabian/GPUE

max_issues_repo_name: hifabian/GPUE

max_forks_repo_name: hifabian/GPUE


text:
using OffsetArrays

coefftype(ring::Type) = eltype(ring)
modulus(ring) = Nemo.modulus(ring)
degree(ring) = Nemo.degree(ring)

struct RingSampler{Ring} <: Random.Sampler{Ring}
    ring::Ring
    coeff_distribution::Any #DiscreteUnivariateDistribution
end
Random.gentype(r::RingSampler) = Any
ring(r::RingSampler) = r.ring

function sample_ring_array(rng, ℛ, coeff_distribution)
    [coefftype(ℛ)(rand(rng, coeff_distribution)) for _ in 1:degree(modulus(ℛ))]
end

function Random.rand(rng::Random.AbstractRNG, r::RingSampler)
    ℛ = ring(r)
    ℛ(OffsetArray(
        sample_ring_array(rng, ℛ, r.coeff_distribution),
        0:degree(modulus(ℛ))-1))
end

struct Zero; end


text:
export SO2,
        matrix,
        generator,
        transformation

using StaticArrays
using LinearAlgebra

mutable struct SO2{T<:AbstractFloat}
    DoF::Int8
    tangent::T
    generator::SMatrix{2, 2, T}

    function SO2(tangent::AbstractFloat)
        generator = SMatrix{2, 2, typeof(tangent)}(0, -1, 1, 0)
        new{typeof(tangent)}(Int8(1), tangent, generator)
    end
end

function generator(group::SO2)
    return group.generator
end

function matrix(group::SO2)
    return group.tangent * SMatrix{2, 2, typeof(group.tangent)}(0, -1, 1, 0)
end

function transformation(group::SO2)
    # println(typeof(group.tangent))
    return exp.(group.tangent * group.generator)
end


text:
using HomotopyContinuation.DoubleDouble
const DD = HomotopyContinuation.DoubleDouble

@testset "DoubleF64" begin
    @test DoubleF64(Float32(2.0)) == DoubleF64(2.0)
    @test DoubleF64(Float16(2.0)) == DoubleF64(2.0)
    @test DoubleF64(2) == DoubleF64(2.0)
    @test DoubleF64(big(π)) isa DoubleF64
    @test DoubleF64(π) == DoubleF64(big(π))

    @test DD.hi(DoubleF64(2.0)) == 2.0
    @test DD.lo(DoubleF64(2.0)) == 0.0

    @test isbitstype(DoubleF64)
    @test isbits(DoubleF64(3))

    @test zero(DoubleF64(2.0)) == zero(DoubleF64)

    @test convert(Integer, DoubleF64(2.0)) isa Int64
    @test 32 ≤ length(string(DoubleF64(rand()))) ≤ 40

    @test convert(DoubleF64, 2.0) isa DoubleF64
    @test convert(DoubleF64, Float32(2.0)) isa DoubleF64
    @test convert(DoubleF64, Float16(2.0)) isa DoubleF64
    @test convert(DoubleF64, 2) isa DoubleF64
    @test convert(DoubleF64, π) isa DoubleF64

    @test promote_type(DoubleF64, Float64) == DoubleF64
    @test promote_type(DoubleF64, Float32) == DoubleF64
    @test promote_type(DoubleF64, Float16) == DoubleF64

    @test DD.wide_sub(rand(), rand()) isa DoubleF64
    @test DD.wide_add(rand(), rand()) isa DoubleF64
    @test DD.wide_div(rand(), rand()) isa DoubleF64
    @test DD.wide_mul(rand(), rand()) isa DoubleF64

    x = DoubleF64(2.321)
    @test x isa DoubleF64
    @test x / 3.2 isa DoubleF64

    @test x + 2 isa DoubleF64
    @test 2 + x isa DoubleF64

    @test im * DoubleF64(rand()) isa Complex{DoubleF64}
    @test im * DoubleF64(rand()) isa ComplexDF64
    @test x < 2 * x
    @test x ≤ x
    @test convert(Float64, x) ≤ x
    @test x ≤ 2 * convert(Float64, x)
    @test DoubleF64(2.0) == 2.0
    @test 2.0 == DoubleF64(2.0)

    @test isnan(DoubleF64(NaN))
    @test isinf(DoubleF64(Inf))

    @test floor(DoubleF64(3.2)) == 3.0
    @test floor(Int, DoubleF64(3.2)) == 3
    @test ceil(DoubleF64(3.2)) == 4.0
    @test ceil(Int, DoubleF64(3.2)) == 4
    @test trunc(DoubleF64(3.2)) == 3.0
    @test trunc(Int, DoubleF64(3.2)) == 3
    @test isinteger(DoubleF64(3.2)) == false
    @test isinteger(DoubleF64(3.0)) == true

    @test DoubleF64(0.0)^0.14 == 1.0

    for k = 1:10
        x = DoubleF64(rand()) * 20 - 10
        y = DoubleF64(rand()) * 20 - 10
        @test x * y ≈ big(x) * big(y) atol = 1e-29
        @test x + y ≈ big(x) + big(y) atol = 1e-30
        @test x - y ≈ big(x) - big(y) atol = 1e-30
        @test x / y ≈ big(x) / big(y) atol = 1e-26
        u = rand()
        @test x / u ≈ big(x) / big(u) atol = 1e-26
        @test x^5 ≈ BigFloat(x)^5 atol = 1e-25

        @test exp(x) ≈ exp(big(x)) atol = 1e-26
        @test log(abs(x)) ≈ log(big(abs(x))) atol = 1e-29
        @test sin(x) ≈ sin(big(x)) atol = 1e-30
        @test cos(x) ≈ cos(big(x)) atol = 1e-30
        @test atan(x) ≈ atan(big(x)) atol = 1e-30
        @test atan(y, x) ≈ atan(big(y), big(x)) atol = 1e-30
        @test tan(x) ≈ tan(big(x)) atol = 1e-21
        sin_x, cos_x = DD.sincos(x)
        @test sin_x ≈ sin(big(x)) atol = 1e-30
        @test cos_x ≈ cos(big(x)) atol = 1e-30

        z = DoubleF64(rand() * 2 - 1)
        @test asin(z) ≈ asin(big(z)) atol = 1e-30
        @test acos(z) ≈ acos(big(z)) atol = 1e-30

        n, p, d = Base.decompose(x)
        @test x ≈ BigFloat(n) * BigFloat(2)^p / BigFloat(d)
    end

    @test eps(DoubleF64(1.0)) == 4.93038065763132e-32 # 2^-104
    @test eps(DoubleF64) == 4.93038065763132e-32 # 2^-104
    @test eps(DoubleF64) == 4.93038065763132e-32 # 2^-104
end


text:
"""
@inbounds variant of bloom_v0

benchmark shows no significant effect.
I tried also a replacement of Base._nthbyte declared @propagate_inbounds,
which did not change results. Maybe LLVM does an excellent job and eliminates
bounds checks within loops by index limit analysis?!

"""
function bloom_v0i end

function bloom_v0i(t::SearchSequence)
    n = sizeof(t)
    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j - 1
        end
    end
    return t,bloom_mask,skip,tlast
end

function bloom_v0i(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_v0i(s,bloom_v0(t),i)
end

function bloom_v0i(s::SearchSequence, p::Tuple,i::Integer,sv::MaybeVector=nothing)
    @inbounds begin
    ( t,bloom_mask,skip,tlast) = p
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)

    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end

    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

#=
    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j - 1
        end
    end
=#
    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i -= 1
    while i <= w
        if DOSTATS loops += 1 end
        if _nthbyte(s,i+n) == tlast
            # check candidate
            j = 0
            while j < n - 1
                if _nthbyte(s,i+j+1) != _nthbyte(t,j+1)
                    break
                end
                j += 1
            end

            # match found
            if j == n - 1
                if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                return i+1
            end

            # no match, try to rule out the next character
            if DOSTATS bloomtests += 1 end
            if i < w && bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                if DOSTATS bloomskips += 1 end
                i += n
            else
                i += skip
            end
        elseif i < w
            if DOSTATS bloomtests += 1 end
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                if DOSTATS bloomskips += 1 end
                i += n
            end
        end
        i += 1
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
    0
end
end


text:
using StaticArrays: SMatrix

using ..Continuations: AbstractContinuationProblem, AbstractContinuationSolver,
    AbstractContinuationCache,
    as, SweepSetup, ContinuationSweep, ContinuationSolution,
    ContinuationCache, ContinuationOptions, ContinuationSolver,
    residual_jacobian!

using ..BifurcationsBase: timekind, Continuous, Discrete
import ..BifurcationsBase: TimeKind

using ..BifurcationsBase
using ..BifurcationsBase:
    AbstractSpecialPoint, SpecialPoint, SpecialPointInterval,
    BifurcationSweep, BifurcationSolution, BifurcationSolver,
    BifurcationProblem, BifurcationCache,
    eigvals_prototype, allocate_sweep!, check_sweep_length, record!
import ..BifurcationsBase: analyze!, re_analyze!

abstract type Codim1Problem{skind, tkind} <: BifurcationProblem{skind, tkind} end

module PointTypes
@enum(
    PointType,
    none,
    # initial,
    # simple_bifurcation,
    saddle_node,
    hopf,
    period_doubling,
)
end  # module
using .PointTypes: PointType

const AbstractCodim1SpecialPoint{tkind <: TimeKind} =
    AbstractSpecialPoint{tkind, PointType}

const Codim1Sweep{tkind <: TimeKind, ptType <: PointType} =
    BifurcationSweep{tkind, ptType}
const Codim1Solution{S <: ContinuationSolution, W <: Codim1Sweep} =
    BifurcationSolution{S, W}
# The constraints `<: TimeKind` and `<: ContinuationSolution` are
# important.  For example, the compiler didn't infer that
# `Codim1Sweep` was more specific than `BifurcationSweep` when `tkind`
# type parameter for `Codim1Sweep` was unconstrained.

BifurcationsBase.point_type_type(::Codim1Problem) = PointType

BifurcationsBase.regular_point(::Type{PointType}) = PointTypes.none

BifurcationsBase.eigvals_prototype(prob::Codim1Problem,
                                   cache::ContinuationCache) =
    cache.u[1:end - 1] .+ 1im
# TODO: improve it for SVector

mutable struct Codim1Cache{P, C <: ContinuationCache{P},
                           JType, eType,
                           } <: BifurcationCache{P}
    super::C
    J::JType
    prev_det::Float64
    eigvals::eType
    point_type::PointType
end

function Codim1Cache(super::C,
                     J::JType,
                     eigvals::eType,
                     point_type = PointTypes.none,
                     ) where {P, C <: ContinuationCache{P},
                              JType, eType}
    return Codim1Cache{P, C, JType, eType}(
        super, J, NaN, eigvals, point_type)
end
# TODO: Remove this constructor after removing the type parameter `P`.

Codim1Cache(prob::Codim1Problem, super::ContinuationCache) =
    Codim1Cache(
        super,
        ds_jacobian(super),
        copy(eigvals_prototype(prob, super)),
    )

BifurcationsBase.BifurcationCache(prob::Codim1Problem,
                                  super::ContinuationCache) =
    Codim1Cache(prob, super)

const Codim1Solver{
        R <: ContinuationSolver,
        P <: Codim1Problem,
        C <: Codim1Cache,
        S <: Codim1Solution,
        } =
    BifurcationSolver{R, P, C, S}

function re_analyze!(solver::Codim1Solver, u::AbstractVector)
    residual_jacobian!(as(solver.cache, ContinuationCache), u)
    analyze!(solver.cache, solver.opts)

    # Suppress special point recording:
    # It's a bit ugly hack... (communicate by sharing!)
    solver.cache.point_type = PointTypes.none  # TODO: FIX!

    record!(solver.sol, solver.cache)
end

function analyze!(cache::Codim1Cache, opts)
    cache.J = J = ds_jacobian(cache)
    eigvals = ds_eigvals(timekind(cache), J)
    cache.point_type = guess_point_type(timekind(cache), cache,
                                        J, eigvals, opts)
    cache.eigvals = eigvals
end

ds_jacobian(solver) = ds_jacobian(as(solver, ContinuationSolver).cache)
ds_jacobian(cache::Codim1Cache) = ds_jacobian(as(cache, ContinuationCache))
ds_jacobian(cache::ContinuationCache) = ds_jacobian(cache.J)
ds_jacobian(HJ::AbstractArray) = @view HJ[:, 1:end-1]
ds_jacobian(HJ::SMatrix) = HJ[:, 1:end-1]
# TOOD: optimize it for StaticArrays using generated functions


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

function stage_update!(Q, Qstages, Rstages, Qhat, Qtt, RKA_explicit, RKA_implicit, dt,
                       ::Val{is}, ::Val{split_nonlinear_linear}, slow_δ,
                       slow_dQ) where {is, split_nonlinear_linear}
  @inbounds @loop for i = (1:length(Q);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    Qhat_i = Q[i]
    Qstages_is_i = -zero(eltype(Q))

    if slow_δ !== nothing
      Rstages[is-1][i] += slow_δ * slow_dQ[i]
    end

    @unroll for js = 1:is-1
      if split_nonlinear_linear
        rkcoeff = RKA_implicit[is, js] / RKA_implicit[is, is]
      else
        rkcoeff = (RKA_implicit[is, js] - RKA_explicit[is, js]) / RKA_implicit[is, is]
      end
      commonterm = rkcoeff * Qstages[js][i]
      Qhat_i += commonterm + dt * RKA_explicit[is, js] * Rstages[js][i]
      Qstages_is_i -= commonterm
    end
    Qstages[is][i] = Qstages_is_i
    Qhat[i] = Qhat_i
    Qtt[i] = Qhat_i
  end
end

function solution_update!(Q, Rstages, RKB, dt, ::Val{nstages}, slow_δ,
                          slow_dQ, slow_scaling) where nstages
  @inbounds @loop for i = (1:length(Q);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    if slow_δ !== nothing
      Rstages[nstages][i] += slow_δ * slow_dQ[i]
    end
    if slow_scaling !== nothing
      slow_dQ[i] *= slow_scaling
    end

    @unroll for is = 1:nstages
      Q[i] += RKB[is] * dt * Rstages[is][i]
    end
  end
end