text:
using Test
using CLIMA.MoistThermodynamics

using CLIMA.PlanetParameters
using LinearAlgebra

@testset "moist thermodynamics - correctness" begin
  FT = Float64
  # ideal gas law
  @test air_pressure(FT(1), FT(1), PhasePartition(FT(1))) === FT(R_v)
  @test air_pressure(FT(1), FT(2), PhasePartition(FT(1), FT(0.5), FT(0))) === FT(R_v)
  @test air_pressure(FT(1), FT(1)) === FT(R_d)
  @test air_pressure(FT(1), FT(2)) === 2*FT(R_d)
  @test air_density(FT(1), FT(1)) === 1/FT(R_d)
  @test air_density(FT(1), FT(2)) === 2/FT(R_d)

  # gas constants and heat capacities
  @test gas_constant_air(PhasePartition(FT(0))) === FT(R_d)
  @test gas_constant_air(PhasePartition(FT(1))) === FT(R_v)
  @test gas_constant_air(PhasePartition(FT(0.5), FT(0.5))) ≈ FT(R_d)/2
  @test gas_constant_air(FT) == FT(R_d)

  @test cp_m(PhasePartition(FT(0))) === FT(cp_d)
  @test cp_m(PhasePartition(FT(1))) === FT(cp_v)
  @test cp_m(PhasePartition(FT(1), FT(1))) === FT(cp_l)
  @test cp_m(PhasePartition(FT(1), FT(0), FT(1))) === FT(cp_i)
  @test cp_m(FT) == FT(cp_d)

  @test cv_m(PhasePartition(FT(0))) === FT(cp_d - R_d)
  @test cv_m(PhasePartition(FT(1))) === FT(cp_v - R_v)
  @test cv_m(PhasePartition(FT(1), FT(1))) === FT(cv_l)
  @test cv_m(PhasePartition(FT(1), FT(0), FT(1))) === FT(cv_i)
  @test cv_m(FT) == FT(cv_d)

  # speed of sound
  @test soundspeed_air(T_0 + 20, PhasePartition(FT(0))) == sqrt(cp_d/cv_d * R_d * (T_0 + 20))
  @test soundspeed_air(T_0 + 100, PhasePartition(FT(1))) == sqrt(cp_v/cv_v * R_v * (T_0 + 100))

  # specific latent heats
  @test latent_heat_vapor(FT(T_0))  ≈ LH_v0
  @test latent_heat_fusion(FT(T_0)) ≈ LH_f0
  @test latent_heat_sublim(FT(T_0)) ≈ LH_s0

  # saturation vapor pressure and specific humidity
  p=FT(1.e5); q_tot=FT(0.23); ρ=FT(1.);
  ρ_v_triple = press_triple / R_v / T_triple;
  @test saturation_vapor_pressure(FT(T_triple), Liquid()) ≈ press_triple
  @test saturation_vapor_pressure(FT(T_triple), Ice()) ≈ press_triple

  @test q_vap_saturation(FT(T_triple), ρ, PhasePartition(FT(0))) == ρ_v_triple / ρ
  @test q_vap_saturation(FT(T_triple), ρ, PhasePartition(q_tot,q_tot)) == ρ_v_triple / ρ

  @test q_vap_saturation_generic(FT(T_triple), ρ; phase=Liquid()) == ρ_v_triple / ρ
  @test q_vap_saturation_generic(FT(T_triple), ρ; phase=Ice()) == ρ_v_triple / ρ
  @test q_vap_saturation_generic.(FT(T_triple-20), ρ; phase=Liquid()) >=
        q_vap_saturation_generic.(FT(T_triple-20), ρ; phase=Ice())

  @test saturation_excess(FT(T_triple), ρ, PhasePartition(q_tot)) ==  q_tot - ρ_v_triple/ρ
  @test saturation_excess(FT(T_triple), ρ, PhasePartition(q_tot/1000)) == 0.0

  # energy functions and inverse (temperature)
  T=FT(300); e_kin=FT(11); e_pot=FT(13);
  @test air_temperature(FT(cv_d*(T-T_0))) === FT(T)
  @test air_temperature(FT(cv_d*(T-T_0)), PhasePartition(FT(0))) === FT(T)

  @test air_temperature(cv_m(PhasePartition(FT(0)))*(T-T_0), PhasePartition(FT(0))) === FT(T)
  @test air_temperature(cv_m(PhasePartition(FT(q_tot)))*(T-T_0) + q_tot*e_int_v0, PhasePartition(q_tot)) ≈ FT(T)

  @test total_energy(FT(e_kin),FT(e_pot), FT(T_0)) === FT(e_kin) + FT(e_pot)
  @test total_energy(FT(e_kin),FT(e_pot), FT(T)) ≈ FT(e_kin) + FT(e_pot) + cv_d*(T-T_0)
  @test total_energy(FT(0),FT(0), FT(T_0), PhasePartition(q_tot)) ≈ q_tot*e_int_v0

  # phase partitioning in equilibrium
  q_liq = FT(0.1);
  T = FT(T_icenuc-10); ρ = FT(1.0); q_tot = FT(0.21)
  @test liquid_fraction_equil(T) === FT(0)
  @test liquid_fraction_equil(T, PhasePartition(q_tot,q_liq,q_liq)) === FT(0.5)
  q = PhasePartition_equil(T, ρ, q_tot)
  @test q.liq ≈ FT(0)
  @test 0 < q.ice <= q_tot

  T = FT(T_freeze+10); ρ = FT(0.1); q_tot = FT(0.60)
  @test liquid_fraction_equil(T) === FT(1)
  @test liquid_fraction_equil(T, PhasePartition(q_tot,q_liq,q_liq/2)) === FT(2/3)
  q = PhasePartition_equil(T, ρ, q_tot)
  @test 0 < q.liq <= q_tot
  @test q.ice ≈ 0

  # saturation adjustment in equilibrium (i.e., given the thermodynamic
  # variables E_int, p, q_tot, compute the temperature and partitioning of the phases
  q_tot = FT(0); ρ = FT(1)
  @test saturation_adjustment(internal_energy_sat(300.0, ρ, q_tot), ρ, q_tot) ≈ 300.0

  q_tot = FT(0.21); ρ = FT(0.1)
  @test saturation_adjustment(internal_energy_sat(200.0, ρ, q_tot), ρ, q_tot) ≈ 200.0
  q = PhasePartition_equil(T, ρ, q_tot)
  @test q.tot - q.liq - q.ice ≈ q_vap_saturation(T, ρ)

  q_tot = FT(0.78); ρ = FT(1)
  @test saturation_adjustment(internal_energy_sat(300.0, ρ, q_tot), ρ, q_tot) ≈ 300.0
  q = PhasePartition_equil(T, ρ, q_tot)
  @test q.tot - q.liq - q.ice ≈ q_vap_saturation(T, ρ)

  ρ = FT(1)
  ρu = FT[1,2,3]
  ρe = FT(1100)
  e_pot = FT(93)
  @test internal_energy(ρ, ρe, ρu, e_pot) ≈ 1000.0

  # potential temperatures
  T = FT(300)
  @test liquid_ice_pottemp_given_pressure(T, FT(MSLP)) === T
  @test liquid_ice_pottemp_given_pressure(T, FT(MSLP)/10) ≈ T * 10^(R_d/cp_d)
  @test liquid_ice_pottemp_given_pressure(T, FT(MSLP)/10, PhasePartition(FT(1))) ≈ T * 10^(R_v/cp_v)

  # dry potential temperatures. FIXME: add correctness tests
  T = FT(300); p=FT(1.e5); q_tot=FT(0.23)
  @test dry_pottemp_given_pressure(T, p, PhasePartition(q_tot)) isa typeof(p)
  @test air_temperature_from_liquid_ice_pottemp_given_pressure(
    dry_pottemp_given_pressure(T, p, PhasePartition(q_tot)), p, PhasePartition(q_tot)) ≈ T

  # Exner function. FIXME: add correctness tests
  p=FT(1.e5); q_tot=FT(0.23)
  @test exner_given_pressure(p, PhasePartition(q_tot)) isa typeof(p)
end

@testset "moist thermodynamics - constructor consistency" begin
  # Make sure `ThermodynamicState` arguments are returned unchanged

  # FIXME/TODO: Increasing `n` to 40 results in broken tests
  # due to failure in saturation adjustment. This needs
  # to be improved/made more robust.
  n = 35
  FT = Float64
  θ_liq_ice = collect(range(250,400,length=n))
  T = collect(range(250,400,length=n))
  ρ = collect(range(0.5,2.0,length=n))
  q_tot_dry = zeros(n)
  q_tot_moist = collect(range(0.0,0.2,length=n))
  q_liq = collect(range(0.0,0.1,length=n))
  q_ice = collect(range(0.0,0.05,length=n))

  e_int = collect(range(10_000.0,200_000.0,length=n))
  p = collect(range(10_000.0,200_000.0,length=n))
  q_pt_dry = PhasePartition.(q_tot_dry, q_liq .* 0, q_ice .* 0)
  q_pt_moist = PhasePartition.(q_tot_moist, q_liq, q_ice)

  # air_temperature_from_liquid_ice_pottemp_given_pressure-liquid_ice_pottemp inverse
  θ_liq_ice = liquid_ice_pottemp_given_pressure.(T, p, q_pt_dry)
  @test all(air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice, p, q_pt_dry) .≈ T)

  # liquid_ice_pottemp-air_temperature_from_liquid_ice_pottemp_given_pressure inverse
  T = air_temperature_from_liquid_ice_pottemp_given_pressure.(θ_liq_ice, p, q_pt_dry)
  @test all(liquid_ice_pottemp_given_pressure.(T, p, q_pt_dry) .≈ θ_liq_ice)

  # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice)
  T_non_linear = air_temperature_from_liquid_ice_pottemp_non_linear.(θ_liq_ice, ρ, q_pt_moist)
  e_int_ = internal_energy.(T_non_linear, q_pt_moist)
  ts = PhaseNonEquil.(e_int_, q_pt_moist, ρ)
  @test all(T_non_linear .≈ air_temperature.(ts))
  @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))

  # PhaseDry
  ts = PhaseDry.(e_int, ρ)
  @test all(internal_energy.(ts) .≈ e_int)
  @test all(air_density.(ts) .≈ ρ)

  # PhaseEquil: Dry case
  ts = PhaseEquil.(e_int, q_tot_dry, ρ)
  @test all(internal_energy.(ts) .≈ e_int)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot_dry)
  @test all(air_density.(ts) .≈ ρ)

  # PhaseEquil: Moist case
  ts = PhaseEquil.(e_int, q_tot_moist, ρ)
  @test all(internal_energy.(ts) .≈ e_int)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot_moist)
  @test all(air_density.(ts) .≈ ρ)

  # PhaseNonEquil: Dry case
  ts = PhaseNonEquil.(e_int, q_pt_dry, ρ)
  @test all(internal_energy.(ts) .≈ e_int)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt_dry, :tot))
  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt_dry, :liq))
  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt_dry, :ice))
  @test all(air_density.(ts) .≈ ρ)

  # PhaseNonEquil: Moist case
  ts = PhaseNonEquil.(e_int, q_pt_moist, ρ)
  @test all(internal_energy.(ts) .≈ e_int)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt_moist, :tot))
  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt_moist, :liq))
  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt_moist, :ice))
  @test all(air_density.(ts) .≈ ρ)

  # LiquidIcePotTempSHumEquil: Dry case
  ts = LiquidIcePotTempSHumEquil.(θ_liq_ice, q_tot_dry, ρ)
  @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
  @test all(air_density.(ts) .≈ ρ)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot_dry)

  # LiquidIcePotTempSHumEquil: Moist case
  ts = LiquidIcePotTempSHumEquil.(θ_liq_ice, q_tot_moist, ρ)
  @test all(abs.(liquid_ice_pottemp.(ts) - θ_liq_ice) .< FT(1e-8))
  @test maximum(abs.(liquid_ice_pottemp.(ts) - θ_liq_ice)) ./ eps(FT) < sqrt(1/(eps(FT)))
  @test all(air_density.(ts) .≈ ρ)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ q_tot_moist)

  # LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice::FT, q_pt::PhasePartition{FT}, p::FT): Moist case
  ts = LiquidIcePotTempSHumNonEquil_given_pressure.(θ_liq_ice, q_pt_moist, p)
  @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
  @test all(air_pressure.(ts) .≈ p)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt_moist,:tot))
  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt_moist,:liq))
  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt_moist,:ice))

  # LiquidIcePotTempSHumNonEquil(θ_liq_ice::FT, q_pt::PhasePartition{FT}, ρ::FT): Moist case
  ts = LiquidIcePotTempSHumNonEquil.(θ_liq_ice, q_pt_moist, ρ)
  @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts))
  @test all(air_density.(ts) .≈ ρ)
  @test all(getproperty.(PhasePartition.(ts),:tot) .≈ getproperty.(q_pt_moist,:tot))
  @test all(getproperty.(PhasePartition.(ts),:liq) .≈ getproperty.(q_pt_moist,:liq))
  @test all(getproperty.(PhasePartition.(ts),:ice) .≈ getproperty.(q_pt_moist,:ice))

end

@testset "moist thermodynamics - type-stability" begin

  FT = Float32
  ρ = FT(1.0)
  p = FT(1000.0*100)

  ρu = FT[1,2,3]
  ρe = FT(1000)
  e_kin = FT(5.0)
  e_pot = FT(100.0)
  e_int = FT(2.0)
  @test internal_energy(ρ, ρe, ρu, e_pot) isa typeof(e_int)

  q_tot = FT(0.01)
  q_ice = FT(0.001)
  q_liq = FT(0.001)
  q_pt = PhasePartition(q_tot, q_liq, q_ice)
  θ_liq_ice = FT(300.0)
  ts_eq = PhaseEquil(e_int, q_tot, ρ)
  ts_dry = PhaseDry(e_int, ρ)
  ts_T = TemperatureSHumEquil(air_temperature(ts_dry), q_pt.tot, air_pressure(ts_dry))
  ts_neq = PhaseNonEquil(e_int, PhasePartition(q_tot, q_liq, q_ice), ρ)
  ts_θ_liq_ice_eq = LiquidIcePotTempSHumEquil(θ_liq_ice, q_tot, ρ)
  ts_θ_liq_ice_neq = LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice, q_pt, p)
  ts_θ_liq_ice_neq_ρ = LiquidIcePotTempSHumNonEquil(θ_liq_ice, q_pt, ρ)
  for ts in (ts_eq, ts_dry, ts_T, ts_neq, ts_θ_liq_ice_eq, ts_θ_liq_ice_neq, ts_θ_liq_ice_neq_ρ)
    @test linearized_air_pressure(e_kin, e_pot, ts) isa typeof(e_int)
    @test soundspeed_air(ts) isa typeof(e_int)
    @test gas_constant_air(ts) isa typeof(e_int)
    @test relative_humidity(ts) isa typeof(e_int)
    @test air_pressure(ts) isa typeof(e_int)
    @test air_density(ts) isa typeof(e_int)
    @test cp_m(ts) isa typeof(e_int)
    @test cv_m(ts) isa typeof(e_int)
    @test eltype(gas_constants(ts)) == typeof(e_int)
    @test air_temperature(ts) isa typeof(e_int)
    @test internal_energy_sat(ts) isa typeof(e_int)
    @test internal_energy(ts) isa typeof(e_int)
    @test latent_heat_vapor(ts) isa typeof(e_int)
    @test latent_heat_sublim(ts) isa typeof(e_int)
    @test latent_heat_fusion(ts) isa typeof(e_int)
    @test q_vap_saturation(ts) isa typeof(e_int)
    @test saturation_excess(ts) isa typeof(e_int)
    @test liquid_fraction_equil(ts) isa typeof(e_int)
    @test liquid_fraction_nonequil(ts) isa typeof(e_int)
    @test PhasePartition(ts) isa PhasePartition{typeof(e_int)}
    @test liquid_ice_pottemp(ts) isa typeof(e_int)
    @test dry_pottemp(ts) isa typeof(e_int)
    @test exner(ts) isa typeof(e_int)
    @test liquid_ice_pottemp_sat(ts) isa typeof(e_int)
    @test specific_volume(ts) isa typeof(e_int)
    @test virtual_pottemp(ts) isa typeof(e_int)
  end
end

@testset "moist thermodynamics - dry limit" begin

  e_int_range = -6.5e4:10000:4.5e4
  ρ_range = 0.1:0.5:2.0

  FT = eltype(ρ_range)
  e_kin = FT(10)
  e_pot = FT(100.0)
  # PhasePartition test is noisy, so do this only once:
  ts_dry = PhaseDry(first(e_int_range), first(ρ_range))
  ts_eq  = PhaseEquil(first(e_int_range), typeof(first(ρ_range))(0), first(ρ_range))
  @test PhasePartition(ts_eq).tot                  ≈ PhasePartition(ts_dry).tot
  @test PhasePartition(ts_eq).liq                  ≈ PhasePartition(ts_dry).liq
  @test PhasePartition(ts_eq).ice                  ≈ PhasePartition(ts_dry).ice

  for e_int in e_int_range
    for ρ in ρ_range
      ts_dry = PhaseDry(e_int, ρ)
      ts_eq  = PhaseEquil(e_int, typeof(ρ)(0), ρ)

      @test linearized_air_pressure(FT(0), e_pot, ts_dry) ≈ air_pressure(ts_dry)
      @test linearized_air_pressure(e_kin, e_pot, ts_eq)  ≈ linearized_air_pressure(e_kin, e_pot, ts_dry)
      @test gas_constant_air(ts_eq)                       ≈ gas_constant_air(ts_dry)
      @test relative_humidity(ts_eq)                      ≈ relative_humidity(ts_dry)
      @test air_pressure(ts_eq)                           ≈ air_pressure(ts_dry)
      @test air_density(ts_eq)                            ≈ air_density(ts_dry)
      @test specific_volume(ts_eq)                        ≈ specific_volume(ts_dry)
      @test cp_m(ts_eq)                                   ≈ cp_m(ts_dry)
      @test cv_m(ts_eq)                                   ≈ cv_m(ts_dry)
      @test air_temperature(ts_eq)                        ≈ air_temperature(ts_dry)
      @test internal_energy(ts_eq)                        ≈ internal_energy(ts_dry)
      @test internal_energy_sat(ts_eq)                    ≈ internal_energy_sat(ts_dry)
      @test soundspeed_air(ts_eq)                         ≈ soundspeed_air(ts_dry)
      @test latent_heat_vapor(ts_eq)                      ≈ latent_heat_vapor(ts_dry)
      @test latent_heat_sublim(ts_eq)                     ≈ latent_heat_sublim(ts_dry)
      @test latent_heat_fusion(ts_eq)                     ≈ latent_heat_fusion(ts_dry)
      @test q_vap_saturation(ts_eq)                       ≈ q_vap_saturation(ts_dry)
      @test saturation_excess(ts_eq)                      ≈ saturation_excess(ts_dry)
      @test liquid_fraction_equil(ts_eq)                  ≈ liquid_fraction_equil(ts_dry)
      @test liquid_fraction_nonequil(ts_eq)               ≈ liquid_fraction_nonequil(ts_dry)
      @test liquid_ice_pottemp(ts_eq)                     ≈ liquid_ice_pottemp(ts_dry)
      @test dry_pottemp(ts_eq)                            ≈ dry_pottemp(ts_dry)
      @test virtual_pottemp(ts_eq)                        ≈ virtual_pottemp(ts_dry)
      @test liquid_ice_pottemp_sat(ts_eq)                 ≈ liquid_ice_pottemp_sat(ts_dry)
      @test exner(ts_eq)                                  ≈ exner(ts_dry)
      @test saturation_vapor_pressure(ts_eq, Ice())       ≈ saturation_vapor_pressure(ts_dry, Ice())
      @test saturation_vapor_pressure(ts_eq, Liquid())    ≈ saturation_vapor_pressure(ts_dry, Liquid())
      @test all(gas_constants(ts_eq)               .≈ gas_constants(ts_dry))
    end
  end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using SossMLJ
using Test

import DataFrames
import Documenter
import MLJBase
import MLJModelInterface
import Soss

include("examples-list.jl")

@testset "SossMLJ.jl" begin
    @testset "Unit tests" begin
        @testset "types.jl" begin
            @testset "default transform" begin
                df = DataFrames.DataFrame()
                df[:a] = [1.0, 2.0, 3.0]
                df[:b] = [4.0, 5.0, 6.0]
                @test SossMLJ.default_transform(df) == (X = [1.0 4.0; 2.0 5.0; 3.0 6.0],)
                @test SossMLJ.default_transform(df).X isa Matrix{Float64}
            end
            @testset "Constructors for the `SossMLJModel` type" begin
                @testset "default response" begin
                    m = Soss.@model begin
                        y ~ Normal()
                        return y
                    end
                    hyperparams = (;)
                    transform = () -> ()
                    infer = Soss.dynamicHMC
                    model = SossMLJModel(;
                        model = m,
                        hyperparams = hyperparams,
                        infer = infer,
                        transform = transform,
                    )
                    @test model.hyperparams == hyperparams
                    @test model.transform == transform
                    @test model.model == m
                    @test model.infer == infer
                    @test model.response isa Symbol
                    @test model.response == :y
                    @test MLJModelInterface.clean!(model) == ""
                end
                @testset "specify response" begin
                    m = Soss.@model begin
                        y ~ Normal()
                        return y
                    end
                    hyperparams = (;)
                    transform = () -> ()
                    infer = Soss.dynamicHMC
                    model = SossMLJModel(;
                        model = m,
                        hyperparams = hyperparams,
                        infer = infer,
                        response = :y,
                        transform = transform,
                    )
                    @test model.hyperparams == hyperparams
                    @test model.transform == transform
                    @test model.model == m
                    @test model.infer == infer
                    @test model.response isa Symbol
                    @test model.response == :y
                end
            end
        end
    end

    @testset "Examples" begin
        for example in EXAMPLES
            @testset "Run example: $(example[1])" begin
                example_file = joinpath(EXAMPLESROOT, "example-$(example[2]).jl")
                extra_example_tests = joinpath(TESTROOT, "extra-example-tests", "$(example[2]).jl")
                @info("Running $(example_file)")
                include(example_file)
                if isfile(extra_example_tests)
                    @info("Running $(extra_example_tests)")
                    include(extra_example_tests)
                end
            end
        end
    end

    @testset "Doctests" begin
        Documenter.doctest(SossMLJ)
    end
end

ext: jl


text:
"""
    count_total_bits(n::Int64)

Calculates the total number of bits from 1 to `n` for a given integer `n` by combing a
while-loop for building the binary subsets with a range for-loop for counting the bits.

# Arguments
- `n::Int64`: Integer value which has to be evaluated.

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.count_total_bits(17)
35
```
"""
function count_total_bits(n::Int64)

    i = 0
    result = 0

    while n > (1 << i)

        k =0
        change = 1 << i

        for _ in 1:n+1
            result += k

            if change == 1
                # Flipping the bit
                if k == 0
                    k = 1
                else
                    k = 0
                end
                change = 1 << i

            else
                change -= 1
            end
        end
        i += 1

    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using HorseML.Regression

@testset "Regression" begin
    mse = LossFunction.mse

    #generating data
    x = 5 .+ 25 .* rand(20)
    t = 170 .- 108 .* exp.(-0.2 .* x) .+ 4 .* rand(20)
    x1 = 23 .* (t ./ 100).^2 .+ 2 .* rand(20)
    x = hcat(x, x1)

    #Linear Regression
    @testset "Linear Regression" begin
        model = LinearRegression()
        @test_throws DimensionMismatch Regression.fit!(model, x', t)
        @test_nowarn Regression.fit!(model, x, t)
        @test size(model(x)) == size(t)
    end

    #Ridge Regression
    @testset "Ridge Regression" begin
        model = Ridge()
        @test_nowarn Regression.fit!(model, x, t)
        @test size(model(x)) == size(t)
    end

    #Lasso Regression
    @testset "Lasso Regression" begin
        model = Lasso()
        @test_nowarn Regression.fit!(model, x, t)
        @test_nowarn size(model(x)) == size(t)
    end

    #Polynomial Regression
    @testset "Polynomial Regression" begin
        model = LinearRegression()
        @test_nowarn x1 = make_design_matrix(x, dims = 2)
        @test_nowarn Regression.fit!(model, x1, t)
        @test size(model(x1)) == size(t)
    end
end


text:
module XRayLabTool
#=
XRayLabTool  Material property when interacted with x-rays.

Usage: result = Refrac([Chemical Formulas],[Energy],[MassDensity])
       result = SubRefrac(Chemical Formula, [Energy], MassDensity)

Input:
    1) Chemical formula. Can be either a vector (using Refrac for multiple formulas) or a single string (using SubRefrac for a single formula). The formula is case sensitive (e.g. CO for Carbon Monoxide vs Co for Cobalt);
    2) Energy (0.03KeV~30KeV). Have to be a vector, even one 1-element.
    3) List of mass densities (g/cm^3). Can be either a vector of numbers (using Refrac for multiple formulas) or a single number (using SubRefrac for a single formula).

Output:
    Can be an array of dictionaries or a single dictionary (if input is single string and using SubRefrac) with the following elements:
    1) Chemical formula;
    2) Molecular weight;
    3) Number of electrons per molecule;
    4) Mass density (g/cm^3);
    5) Electron density (1/A^3);
    6) Real part of the atomic scattering factor (f1)
    7) Imaginary part of the atomic scattering factor (f2);
    8) X-ray energy (KeV);
    9) Corresponding X-ray wavelength (A);
    10) Dispersion;
    11) Absorption;
    12) Critical angle (degree);
    13) Attenuation length (cm);
    14) Real part of scattering length density (SLD) (A^-2);
    15) Imaginary part of SLD (A^-2).

Example 1: >> result=Refrac(["H2O","Si3N4"],[8.04778],[1,3.44])
           Output: a Dict array
               result["Si3N4"] =
               Dict{String, Any}(
                    "Critical Angle" => [0.2707377852086883],
                    "reSLD" => [2.955439361638965e-5],
                    "MW" => 140.283,
                    "Electron Density" => 1.0337188334951493,
                    "Wavelength" => [1.5406009069367983],
                    "f2" => [1.0482169679625088],
                    "Formula" => "Si3N4",
                    "Number Of Electrons" => 70.0,
                    "f1" => [71.0208511434114],
                    "Attenuation Length" => [0.0074403349864984964],
                    "Energy" => [8.04778],
                    "Dispersion" => [1.116406825816022e-5],
                    "imSLD" => [4.362017121420607e-7],
                    "Density" => 3.44,
                    "Absorption" => [1.6477366282283455e-7])
                result["H2O"] =
                Dict{String, Any}(
                    "Critical Angle" => [0.1532479196629315],
                    "reSLD" => [9.469204261542127e-6],
                    "MW" => 18.015,
                    "Electron Density" => 0.3342848731612545,
                    "Wavelength" => [1.5406009069367983],
                    "f2" => [0.033709251629429504],
                    "Formula" => "H2O",
                    "Number Of Electrons" => 10.0,
                    "f1" => [10.052289064465802],
                    "Attenuation Length" => [0.10220737349289194],
                    "Energy" => [8.04778],
                    "Dispersion" => [3.576958610569931e-6],
                    "imSLD" => [3.175394053391685e-8],
                    "Density" => 1.0,
                    "Absorption" => [1.1994939371370336e-8])

Example 2: >>result=SubRefrac("SiO2",Vector(8:0.5:10),2.33)
          Output: Dict
           result =
           Dict{String, Any} with 15 entries:
              "Critical Angle"      => [0.224095, 0.210836, 0.199074, 0.188608, 0.179106]
              "reSLD"               => [2.00086e-5, 1.9994e-5, 1.99843e-5, 1.99867e-5, 1.99707e-5]
              "MW"                  => 60.083
              "Electron Density"    => 0.70061
              "Wavelength"          => [1.5498, 1.45864, 1.3776, 1.3051, 1.23984]
              "f2"                  => [0.396933, 0.352482, 0.314986, 0.283079, 0.255704]
              "Formula"             => "SiO2"
              "Number Of Electrons" => 30.0
              "f1"                  => [30.4039, 30.3817, 30.367, 30.3706, 30.3463]
              "Attenuation Length"  => [0.0123506, 0.0147774, 0.0175093, 0.0205652, 0.0239651]
              "Energy"              => [8.0, 8.5, 9.0, 9.5, 10.0]
              ⋮                     => ⋮

For more information about X-ray interactions with matter, go to
   http://www.cxro.lbl.gov
   http://www.nist.gov/

Atomic scattering factor table is taken from the above two websites.

This is translated from Zhang Jiang's MATLAB script
=#

using BenchmarkTools
using CSV
using DataFrames
using Interpolations
using LinearAlgebra
using PCHIPInterpolation
using PeriodicTable
using SpecialFunctions
using Unitful
using Base.Threads

export Refrac, SubRefrac

function Refrac(formulaList, energy, massDensityList)

    if !isa(formulaList, Vector{String})
        println("Need a styring for chemical formula input argument.")
        return
    end
    if !isa(energy, Vector) | isempty(energy)
        println("Invalid x-ray energy.")
        return
    end
    if count(!iszero, energy .< 0.03) | count(!iszero, energy .> 30) > 0
        println("Energy is out of range 0.03KeV~30KeV.")
        return
    end
    if !isa(massDensityList, Vector)
        println("Invalid mass density.")
        return
    end
    if length(formulaList) != length(massDensityList)
        println("Input arguements do not match.")
        return
    end

    # --- sort energy list and
    energy = sort!(energy)  # sort energy from min to max

    result = Dict()

    @threads for ii = 1:length(formulaList)
        result[formulaList[ii]] =
            SubRefrac(formulaList[ii], energy, massDensityList[ii])
    end
    return result
end

function SubRefrac(formulaStr, energy, massDensity)
    # Some physical constants
    thompson = 2.8179403227e-15 # m
    speedoflight = 299792458 # m/sec
    plank = 6.626068e-34 # m²kg/sec
    elementcharge = 1.60217646e-19 # Coulombs
    avogadro = 6.02214199e23 # mole⁻¹

    # convert energy to wavelength
    if !isa(energy, Vector)
        println("Invalid x-ray energy.")
    end
    wavelength = (speedoflight * plank / elementcharge) ./ (energy * 1000.0)

    # determine elements and number of atoms in the formula
    nElements = 0
    formulaElement = []
    nAtoms = []

    formulaChar = split(formulaStr, "")

    for iFormula = 1:length(formulaStr)
        if count(!iszero, formulaChar[iFormula] <= "Z") &
           count(!iszero, formulaChar[iFormula] >= "A") > 0
            nElements = nElements + 1
            push!(formulaElement, formulaChar[iFormula])
            push!(nAtoms, "0")
        elseif count(!iszero, formulaChar[iFormula] <= "z") &
               count(!iszero, formulaChar[iFormula] >= "a") &
               (
                   (
                       count(!iszero, string(formulaChar[iFormula-1]) <= "Z") &
                       count(!iszero, string(formulaChar[iFormula-1]) >= "A")
                   ) | (
                       count(!iszero, string(formulaChar[iFormula-1]) <= "z") &
                       count(!iszero, string(formulaChar[iFormula-1]) >= "a")
                   )
               ) > 0
            formulaElement[nElements] =
                join([formulaElement[nElements], formulaChar[iFormula]])
        elseif (
            count(!iszero, formulaChar[iFormula] <= "9") &
            count(!iszero, formulaChar[iFormula] >= "0")
        ) | count(!iszero, formulaChar[iFormula] == ".") > 0
            nAtoms[nElements] = join([nAtoms[nElements], formulaChar[iFormula]])
        else
            println("Invalid chemical formula.")
            return
        end
    end

    for iElements = 1:nElements
        nAtoms[iElements] = parse(Float64, nAtoms[iElements])
        if abs.(nAtoms[iElements]) .- 0.0 < eps()
            nAtoms[iElements] = 1.0
        end
    end

    # read f1 and f2 from tables
    f1f2Table = []
    for iElements = 1:nElements
        fname = join([lowercase(formulaElement[iElements]), ".nff"])
        file = normpath(joinpath(@__DIR__, "AtomicScatteringFactor", fname))
        # println(file)
        try
            table = CSV.File(file) |> DataFrame
            push!(f1f2Table, table)
        catch
            prtstr = formulaElement[iElements]
            println("Element $prtstr is NOT in the table list.")
        end
    end

    # determine the atomic number and atomic weight
    atomicNumber = []
    atomicWeight = []
    for iElements = 1:nElements
        AN = findall(
            x -> x == formulaElement[iElements],
            [elements[iAtomicnum].symbol for iAtomicnum = 1:length(elements)],
        )
        push!(atomicNumber, AN[1])
        push!(atomicWeight, elements[atomicNumber[iElements]].atomic_mass)
    end

    # determine molecular weight and number of electrons
    molecularWeight = 0
    numberOfElectrons = 0
    for iElements = 1:nElements
        molecularWeight =
            molecularWeight +
            nAtoms[iElements] * ustrip(atomicWeight[iElements])
        numberOfElectrons =
            numberOfElectrons + atomicNumber[iElements] * nAtoms[iElements]
    end

    # interpolate to get f1 and f2 for given energies
    interpf1 = []
    interpf2 = []
    for iElements = 1:nElements
        # use PCHIPInterpolation
        itp1 = Interpolator(f1f2Table[iElements].E, f1f2Table[iElements].f1)
        itp2 = Interpolator(f1f2Table[iElements].E, f1f2Table[iElements].f2)
        # use Interpolations
        # itp1 = interpolate((f1f2Table[iElements].E,), f1f2Table[iElements].f1, Gridded(Linear()))
        # itp2 = interpolate((f1f2Table[iElements].E,), f1f2Table[iElements].f2, Gridded(Linear()))
        push!(interpf1, [itp1(E) for E in energy * 1000])
        push!(interpf2, [itp2(E) for E in energy * 1000])
    end

    # calculate dispersion, absorption, critical angle and attenuation length
    Dispersion = zeros(length(energy))
    Absorption = zeros(length(energy))
    f1 = zeros(length(energy))
    f2 = zeros(length(energy))
    CriticalAngle = zeros(length(energy))
    ElectronDensity =
        1e6 * massDensity / molecularWeight * avogadro * numberOfElectrons /
        1e30
    for iElements = 1:nElements
        Dispersion =
            Dispersion +
            wavelength .^ 2 / (2 * pi) *
            thompson *
            avogadro *
            massDensity *
            1e6 / molecularWeight * nAtoms[iElements] .* interpf1[iElements]
        Absorption =
            Absorption +
            wavelength .^ 2 / (2 * pi) *
            thompson *
            avogadro *
            massDensity *
            1e6 / molecularWeight * nAtoms[iElements] .* interpf2[iElements]
        f1 = f1 + nAtoms[iElements] .* interpf1[iElements]
        f2 = f2 + nAtoms[iElements] .* interpf2[iElements]
    end
    CriticalAngle = sqrt.((2 * Dispersion)) * 180 / pi
    AttLength = wavelength ./ Absorption / (4 * pi) * 1e2
    reSLD = Dispersion * 2 * pi ./ wavelength .^ 2 / 1e20
    imSLD = Absorption * 2 * pi ./ wavelength .^ 2 / 1e20
    Xresult = Dict(
        "Formula" => formulaStr,
        "MW" => molecularWeight,
        "Number Of Electrons" => numberOfElectrons,
        "Density" => massDensity,
        "Electron Density" => ElectronDensity,
        "Energy" => energy,
        "Wavelength" => wavelength * 1e10,
        "Dispersion" => Dispersion,
        "Absorption" => Absorption,
        "f1" => f1,
        "f2" => f2,
        "Critical Angle" => CriticalAngle,
        "Attenuation Length" => AttLength,
        "reSLD" => reSLD,
        "imSLD" => imSLD,
    )
    return Xresult
end
end


text:
# Reverse interval constraint propagation using Symbolics

# import Base: ∩
# import Base: ∪

#

# @register a ∩ b
# @register a ∪ b

# Possibly should replace all calls to `rev` with calls to the actual
# reverse functions instead for speed

remove_constant(s) = (s isa Num || s isa Sym) ? s : Variable(:_)
remove_parameters(s, params) = (any(x -> isequal(x, s), params)) ? Variable(:_) : s

function rev(eq::Assignment, params)

    # @show eq, params

    vars = tuple(args(eq)...)
    return_vars = remove_constant.(tuple(lhs(eq), vars...))
    # return_vars = remove_constant.(tuple(lhs(eq), vars...))

    return_vars = remove_parameters.(return_vars, Ref(params))

    reverse = rev(op(eq), lhs(eq), vars...)

    return Assignment(Symbolics.MakeTuple(return_vars), reverse)

end

@register rev(a::Any, b, c, d)
@register rev(a::Any, b, c)

# difference between reverse mode AD and reverse propagation:
# reverse mode AD introduces *new* variables
# reverse propagation can use the *same* variables

# x ~ x ∩ (z - y)

# reverse ops from IntervalContractors:

const binary_functions = Dict(
                    :+     => :plus_rev,
                    :-     => :minus_rev,
                    :*     => :mul_rev,
                    :/     => :div_rev,
                    :^     => :power_rev,
                    );

for (f, f_rev) in binary_functions
    @eval rev(::typeof($f), z::Real, x::Real, y::Real) = $f_rev(z, x, y)
end

const unary_functions = [:sqrt, :abs,
            :exp, :exp2, :exp10, :expm1,
            :log, :log2, :log10, :log1p,
            :sin, :cos, :tan,
            :asin, :acos, :atan,
            :sinh, :cosh, :tanh,
            :asinh, :acosh, :atanh,
            :inv, :sign, :max, :min];

for f in unary_functions
    f_rev = Symbol(f, :_rev)
    @eval rev(::typeof($f), z::Real, x::Real) = $f_rev(z, x)
end

"Generate code (as Symbolics.Assignment) for forward--backward (HC4Revise) contractor"
function forward_backward_code(ex, vars, params=[])

    final_var = make_variable(:value)  # to record the output of running the forward interval function
    constraint_var = make_variable(:constraint)   # symbolic constraint variable

    forward_code, last = cse_equations(ex)

    # @show forward_code

    # @show constraint_var, final_var

    constraint_code = [Assignment(final_var, last),
                        Assignment(last, last ∩ constraint_var)]

    reverse_code = rev.(reverse(forward_code), Ref(params))

    code = forward_code ∪ constraint_code ∪ reverse_code

    return code, final_var, constraint_var
end

# code, final, constraint = forward_backward_code(x^2 + a*y^2, [x, y], [a])

# code

Symbolics.toexpr(t::Tuple) = Symbolics.toexpr(Symbolics.MakeTuple(t))

# vars = @variables x, y

# ex = x^2 + a*y^2
# code, final, constraint = forward_backward_code(ex, vars, [a])

# code
# final
# constraint

# dump(toexpr.(code))

# function forward_backward_expr(vars, ex)
#     symbolic_code, final, constraint = forward_backward_code(vars, ex)

#     code = toexpr.(symbolic_code)

#     return_tuple = toexpr(vars)
#     # push!(code.args, :(return $return_tuple))

#     return code, final, return_tuple
# end

"Build Julia code for forward_backward contractor"
function forward_backward_expr(ex, vars, params=[])

    symbolic_code, final_var, constraint_var = forward_backward_code(ex, vars, params)

    # @show symbolic_code

    code = toexpr.(symbolic_code)
    all_code = Expr(:block, code...)

    return all_code, final_var, constraint_var
end

# forward_backward_expr(ex, vars, [a])

function forward_backward_contractor(ex, vars, params=[])

    code, final_var, constraint_var = forward_backward_expr(ex, vars, params)

    input_vars = toexpr(Symbolics.MakeTuple(vars))
    final = toexpr(final_var)

    constraint = toexpr(constraint_var)

    if !isempty(params)
        params_tuple = toexpr(Symbolics.MakeTuple(params))

        function_code =
            quote
                ($input_vars, $constraint, $params_tuple) -> begin
                    $code
                    return $input_vars, $(final)
                end
            end

    else

        function_code =
            quote
                ($input_vars, $constraint) -> begin
                    $code
                    return $input_vars, $(final)
                end
            end

    end

    return eval(function_code)

end

# ex = x^2 + a * y^2
# C = forward_backward_contractor(ex, vars, [a])

# const CC2 = forward_backward_contractor(ex, vars, [a])

# CC2((-10..10, -10..10), 0..1, 5)
# @btime CC2((-10..10, -10..10), 0..1, 6)
# @btime CC2((-10..10, -10..10), 0..1, 7)

# using BenchmarkTools

# @btime CC((-10..10, -10..10), 0..1)

# @code_native C((-10..10, -10..10), 0..1)

# CC(IntervalBox(-10..10, 2), 0..1)

# vars = @variables x, y
# ex = 3x^2 + 4x^2 * y

# code, final, return_tuple = forward_backward_code(vars, 3x^2 + 4x^2 * y)

# expr = forward_backward_expr(vars, ex)

# C = forward_backward_contractor(vars, ex)

# using IntervalArithmetic

# C( (1..2, 3..4) )

# ex = x^2 + y^2

# C( (-10..10, -10..10) )


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: codestar_notifications
using AWS.Compat
using AWS.UUIDs

"""
    create_notification_rule(detail_type, event_type_ids, name, resource, targets)
    create_notification_rule(detail_type, event_type_ids, name, resource, targets, params::Dict{String,<:Any})

Creates a notification rule for a resource. The rule specifies the events you want
notifications about and the targets (such as SNS topics) where you want to receive them.

# Arguments
- `detail_type`: The level of detail to include in the notifications for this resource.
  BASIC will include only the contents of the event as it would appear in AWS CloudWatch.
  FULL will include any supplemental information provided by AWS CodeStar Notifications
  and/or the service for the resource for which the notification is created.
- `event_type_ids`: A list of event types associated with this notification rule. For a
  list of allowed events, see EventTypeSummary.
- `name`: The name for the notification rule. Notifictaion rule names must be unique in
  your AWS account.
- `resource`: The Amazon Resource Name (ARN) of the resource to associate with the
  notification rule. Supported resources include pipelines in AWS CodePipeline, repositories
  in AWS CodeCommit, and build projects in AWS CodeBuild.
- `targets`: A list of Amazon Resource Names (ARNs) of SNS topics to associate with the
  notification rule.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`: A unique, client-generated idempotency token that, when provided
  in a request, ensures the request cannot be repeated with a changed parameter. If a request
  with the same parameters is received and a token is included, the request returns
  information about the initial request that used that token.  The AWS SDKs prepopulate
  client request tokens. If you are using an AWS SDK, an idempotency token is created for
  you.
- `"Status"`: The status of the notification rule. The default value is ENABLED. If the
  status is set to DISABLED, notifications aren't sent for the notification rule.
- `"Tags"`: A list of tags to apply to this notification rule. Key names cannot start with
  \"aws\".
"""
create_notification_rule(DetailType, EventTypeIds, Name, Resource, Targets; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/createNotificationRule", Dict{String, Any}("DetailType"=>DetailType, "EventTypeIds"=>EventTypeIds, "Name"=>Name, "Resource"=>Resource, "Targets"=>Targets, "ClientRequestToken"=>string(uuid4())); aws_config=aws_config)
create_notification_rule(DetailType, EventTypeIds, Name, Resource, Targets, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/createNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DetailType"=>DetailType, "EventTypeIds"=>EventTypeIds, "Name"=>Name, "Resource"=>Resource, "Targets"=>Targets, "ClientRequestToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_notification_rule(arn)
    delete_notification_rule(arn, params::Dict{String,<:Any})

Deletes a notification rule for a resource.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule you want to delete.

"""
delete_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
delete_notification_rule(Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), params)); aws_config=aws_config)

"""
    delete_target(target_address)
    delete_target(target_address, params::Dict{String,<:Any})

Deletes a specified target for notifications.

# Arguments
- `target_address`: The Amazon Resource Name (ARN) of the SNS topic to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ForceUnsubscribeAll"`: A Boolean value that can be used to delete all associations with
  this SNS topic. The default value is FALSE. If set to TRUE, all associations between that
  target and every notification rule in your AWS account are deleted.
"""
delete_target(TargetAddress; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteTarget", Dict{String, Any}("TargetAddress"=>TargetAddress); aws_config=aws_config)
delete_target(TargetAddress, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteTarget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("TargetAddress"=>TargetAddress), params)); aws_config=aws_config)

"""
    describe_notification_rule(arn)
    describe_notification_rule(arn, params::Dict{String,<:Any})

Returns information about a specified notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.

"""
describe_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/describeNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
describe_notification_rule(Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/describeNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), params)); aws_config=aws_config)

"""
    list_event_types()
    list_event_types(params::Dict{String,<:Any})

Returns information about the event types available for configuring notifications.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  default number is 50. The maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
list_event_types(; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listEventTypes"; aws_config=aws_config)
list_event_types(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listEventTypes", params; aws_config=aws_config)

"""
    list_notification_rules()
    list_notification_rules(params::Dict{String,<:Any})

Returns a list of the notification rules for an AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type. For
  valid values, see ListNotificationRulesFilter.  A filter with the same name can appear more
  than once when used with OR statements. Filters with different names should be applied with
  AND statements.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
list_notification_rules(; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listNotificationRules"; aws_config=aws_config)
list_notification_rules(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listNotificationRules", params; aws_config=aws_config)

"""
    list_tags_for_resource(arn)
    list_tags_for_resource(arn, params::Dict{String,<:Any})

Returns a list of the tags associated with a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) for the notification rule.

"""
list_tags_for_resource(Arn; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTagsForResource", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
list_tags_for_resource(Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), params)); aws_config=aws_config)

"""
    list_targets()
    list_targets(params::Dict{String,<:Any})

Returns a list of the notification rule targets for an AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type. Valid
  filters include target type, target address, and target status.  A filter with the same
  name can appear more than once when used with OR statements. Filters with different names
  should be applied with AND statements.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
list_targets(; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTargets"; aws_config=aws_config)
list_targets(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTargets", params; aws_config=aws_config)

"""
    subscribe(arn, target)
    subscribe(arn, target, params::Dict{String,<:Any})

Creates an association between a notification rule and an SNS topic so that the associated
target can receive notifications when the events described in the rule are triggered.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule for which you want to
  create the association.
- `target`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`: An enumeration token that, when provided in a request, returns
  the next batch of the results.
"""
subscribe(Arn, Target; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/subscribe", Dict{String, Any}("Arn"=>Arn, "Target"=>Target); aws_config=aws_config)
subscribe(Arn, Target, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/subscribe", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "Target"=>Target), params)); aws_config=aws_config)

"""
    tag_resource(arn, tags)
    tag_resource(arn, tags, params::Dict{String,<:Any})

Associates a set of provided tags with a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule to tag.
- `tags`: The list of tags to associate with the resource. Tag key names cannot start with
  \"aws\".

"""
tag_resource(Arn, Tags; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/tagResource", Dict{String, Any}("Arn"=>Arn, "Tags"=>Tags); aws_config=aws_config)
tag_resource(Arn, Tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/tagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "Tags"=>Tags), params)); aws_config=aws_config)

"""
    unsubscribe(arn, target_address)
    unsubscribe(arn, target_address, params::Dict{String,<:Any})

Removes an association between a notification rule and an Amazon SNS topic so that
subscribers to that topic stop receiving notifications when the events described in the
rule are triggered.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.
- `target_address`: The ARN of the SNS topic to unsubscribe from the notification rule.

"""
unsubscribe(Arn, TargetAddress; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/unsubscribe", Dict{String, Any}("Arn"=>Arn, "TargetAddress"=>TargetAddress); aws_config=aws_config)
unsubscribe(Arn, TargetAddress, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/unsubscribe", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "TargetAddress"=>TargetAddress), params)); aws_config=aws_config)

"""
    untag_resource(arn, tag_keys)
    untag_resource(arn, tag_keys, params::Dict{String,<:Any})

Removes the association between one or more provided tags and a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule from which to remove the
  tags.
- `tag_keys`: The key names of the tags to remove.

"""
untag_resource(Arn, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/untagResource", Dict{String, Any}("Arn"=>Arn, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(Arn, TagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/untagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "TagKeys"=>TagKeys), params)); aws_config=aws_config)

"""
    update_notification_rule(arn)
    update_notification_rule(arn, params::Dict{String,<:Any})

Updates a notification rule for a resource. You can change the events that trigger the
notification rule, the status of the rule, and the targets that receive the notifications.
To add or remove tags for a notification rule, you must use TagResource and UntagResource.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DetailType"`: The level of detail to include in the notifications for this resource.
  BASIC will include only the contents of the event as it would appear in AWS CloudWatch.
  FULL will include any supplemental information provided by AWS CodeStar Notifications
  and/or the service for the resource for which the notification is created.
- `"EventTypeIds"`: A list of event types associated with this notification rule.
- `"Name"`: The name of the notification rule.
- `"Status"`: The status of the notification rule. Valid statuses include enabled (sending
  notifications) or disabled (not sending notifications).
- `"Targets"`: The address and type of the targets to receive notifications from this
  notification rule.
"""
update_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/updateNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
update_notification_rule(Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = codestar_notifications("POST", "/updateNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), params)); aws_config=aws_config)


text:
using FewShotAnomalyDetection
using Adapt
using Flux
using MLDataPattern
using JLD2
using FileIO
using FluxExtensions
using ADatasets
using StatsBase

folderpath = "D:/dev/julia/"
# folderpath = "/home/bimjan/dev/julia/

"""
    gridSearch(f, parameters...)
Maps `f` to product of `parameters`.
"""
gridSearch(f, parameters...) = map(p -> printAndRun(f, p), Base.product(parameters...))

function printAndRun(f, p)
    println(p)
    (p, f(p))
end

function createSVAEWithMem(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, memorySize, k, labelCount, β, α = 0.1, T = Float64)
    encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, hiddenDim, numLayers - 1, nonlinearity, "", layerType))
    decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))

    svae = SVAEbase(encoder, decoder, hiddenDim, latentDim, T)
    # train!, classify, trainOnLatent! = augmentModelWithMemory((x) -> zfromx(svae, x), memorySize, latentDim, k, labelCount, α, T)

    function learnRepresentation!(data, labels)
        # return wloss(svae, data, β, (x, y) -> mmd_imq(x, y, 1))
        return loss(svae, data, β)
    end

    function learnAnomaly!(data, labels)
        # trainOnLatent!(zfromx(svae, data), labels)
        return wloss(svae, data, β, (x, y) -> FewShotAnomalyDetection.mmd_imq(x, y, 1))
        # return loss(svae, data, β)
    end

    return svae, learnRepresentation!, learnAnomaly!
end

function runExperiment(datasetName, train, test, createModel, anomalyCounts, batchSize = 100, numBatches = 10000)
    (model, learnRepresentation!, learnAnomaly!) = createModel()
    opt = Flux.Optimise.ADAM(Flux.params(model), 1e-3)
    cb = Flux.throttle(() -> println("$datasetName : $(learnAnomaly!(train[1], zero(train[2]) .+ 2))"), 5)
    Flux.train!(learnAnomaly!, RandomBatches((train[1], zero(train[2]) .+ 2), batchSize, numBatches), opt, cb = cb)
    # FluxExtensions.learn(learnRepresentation!, opt, RandomBatches((train[1], train[2] .- 1), batchSize, numBatches), ()->(), 100)

    learnRepresentation!(train[1], train[2] .- 1)

    results = []
    anomalies = train[1][:, train[2] .- 1 .== 1] # TODO needs to be shuffled!!!
    for ac in anomalyCounts
        if ac <= size(anomalies, 2)
            l = learnAnomaly!(anomalies[:, ac], [1])
        else
            break;
        end

        values, probScore = classify(test[1])
        values = Flux.Tracker.data(values)
        probScore = Flux.Tracker.data(probScore)

        rocData = roc(test[2] .- 1, values)
        f1 = f1score(rocData)
        # auc = EvalCurves.auc(EvalCurves.roccurve(probScore, test[2] .- 1)...)
        auc = pyauc(test[2] .- 1, probScore)
        push!(results, (ac, f1, auc, values, probScore, rstrn, rstst, knnauc, knnprec, knnrecall))
    end
    return results
end

outputFolder = folderpath * "OSL/experiments/NewImplementationTests/"
mkpath(outputFolder)

# datasets = ["breast-cancer-wisconsin", "sonar", "wall-following-robot", "waveform-1"]
datasets = ["abalone"]
difficulties = ["easy"]
const dataPath = folderpath * "data/loda/public/datasets/numerical"
batchSize = 100
iterations = 10000

loadData(datasetName, difficulty) =  ADatasets.makeset(ADatasets.loaddataset(datasetName, difficulty, dataPath)..., 0.8, "low")

if length(ARGS) != 0
    datasets = [ARGS[1]]
    difficulties = ["easy"]
end

for (dn, df) in zip(datasets, difficulties)
    train, test, clusterdness = loadData(dn, df)

    println("$dn")
    println("Running svae...")

    evaluateOneConfig = p -> runExperiment(dn, train, test, () -> createSVAEWithMem(size(train[1], 1), p...), 1:5, batchSize, iterations)
    results = gridSearch(evaluateOneConfig, [32], [8], [3], ["relu"], ["Dense"], [16 32 64 128 256], [16 32], [1], [0.01 0.05 0.1])
    results = reshape(results, length(results), 1)
    save(outputFolder * dn * "-svae.jld2", "results", results)
end


text:
#
# Describes Restricted Two-body orbits through Cartesian coordinates and Orbital Elements.
#

"""
Abstract type for all four conic sections.
"""
abstract type AbstractConic end

"""
Type for orbits in the circular conic section.
"""
struct Circular <: AbstractConic end

"""
Type for orbits in the elliptical conic section.
"""
struct Elliptical <: AbstractConic end

"""
Type for orbits in the parabolic conic section.
"""
struct Parabolic <: AbstractConic end

"""
Type for orbits in the hyperbolic conic section.
"""
struct Hyperbolic <: AbstractConic end

"""
Type for invalid orbits (orbits with NaN fields)
"""
struct Invalid <: AbstractConic end

"""
Type representing large bodies in space. Currently, the following
solar system bodies are supported:

Sun, Mercury, Venus, Earth, Moon (Luna), Mars, Jupiter,
Saturn, Uranus, Neptune, Pluto.
"""
struct RestrictedTwoBodySystem{F, LU, TU} <: AbstractSystem{F, LU, TU}
    μ::F
    name::String

    function RestrictedTwoBodySystem(μ::Real, name=""; lengthunit = u"km", timeunit = u"s")
        T = typeof(μ)
        if !(T <: AbstractFloat)
            @warn "Non-float type $(string(T)) provided. Defaulting to Float64."
            T = Float64
        end
        return new{T, typeof(lengthunit), typeof(timeunit)}(T(μ), name)
    end

    function RestrictedTwoBodySystem(μ::MassParameter, name="")
        T = typeof(ustrip(μ))
        if !(T <: AbstractFloat)
            @warn "Non-float type $(string(T)) provided. Defaulting to Float64."
            T = Float64
        end

        M = unit(μ)
        lengthaxis = findfirst(T -> T isa Unitful.Dimension{:Length}, collect(typeof(typeof(M).parameters[2]).parameters[1]))
        timeaxis   = findfirst(T -> T isa Unitful.Dimension{:Time}, collect(typeof(typeof(M).parameters[2]).parameters[1]))

        lengthunit = (Unitful.FreeUnits{(typeof(M).parameters[1][lengthaxis], ), Unitful.𝐋^3})()^(1//3)
        timeunit   = (Unitful.FreeUnits{(typeof(M).parameters[1][timeaxis],), Unitful.𝐓^2})()^(-1//2)

        return new{T, typeof(lengthunit), typeof(timeunit)}(T(ustrip(lengthunit^3 / timeunit^2, μ)), name)
    end
end

"""
Conversions between `RestrictedTwoBodySystem` instances.
"""
function Base.convert(::Type{RestrictedTwoBodySystem{F, LU, TU}}, sys::RestrictedTwoBodySystem) where {F, LU, TU}
    return RestrictedTwoBodySystem(ustrip(LU()^3 / TU()^2, mass_parameter(sys)), sys.name; lengthunit = LU(), timeunit = TU())
end

"""
Custom display for `RestrictedTwoBodySystem` instances.
"""
function Base.show(io::IO, body::RestrictedTwoBodySystem{F,LU,TU}) where {F,LU,TU}
    println(io, "  Restricted Two-body System ", body.name == "" ? "" : string("(",body.name,")"), ":")
    println(io, "")
    println(io, "    μ = ", body.μ, " ", string(LU()^3 / TU()^2))
end

"""
Struct for storing `Keplerian` states.
"""
struct KeplerianState{F, LU, TU, AU} <: AbstractState{F, LU, TU, Inertial} where {AU<:Unitful.Units{U, NoDims, nothing} where U}
    t::F
    e::F
    a::F
    i::F
    Ω::F
    ω::F
    ν::F

    function KeplerianState(e::E, a::A, i::I, Ω::O, ω::W, ν::V, epoch::Real=0; lengthunit = u"km", timeunit = u"s", angularunit = u"rad") where {
        E <: Real, A <: Real, I <: Real, O <: Real, W <: Real, V <: Real
    }
        F = promote_type(E, A, I, O, W, V, typeof(epoch))
        if !(F <: AbstractFloat)
            @warn "Promoted type ($(string(F)) is not a float: defaulting to Float64."
            F = Float64
        end
        return new{F, typeof(lengthunit), typeof(u"s"), typeof(angularunit)}(epoch, e, a, i, Ω, ω, ν)
    end

    function KeplerianState(e::Real, a::Unitful.Length, i::DimensionlessQuantity,
                            Ω::DimensionlessQuantity, ω::DimensionlessQuantity,
                            ν::DimensionlessQuantity, epoch::Unitful.Time = 0.0u"s")

        F = promote_type(typeof(e), typeof(ustrip(a)), typeof(ustrip(i)), typeof(ustrip(Ω)), typeof(ustrip(ω)), typeof(ustrip(ν)), typeof(ustrip(epoch)))
        if !(F <: AbstractFloat)
            @warn "Non-float parameters provided. Defaulting to Float64."
            F = Float64
        end
        lengthunit  = unit(a)
        angularunit = u"rad" ∈ unit.((i, Ω, ω, ν)) ? u"rad" : unit(i)
        timeunit = unit(epoch)
        return KeplerianState(e, ustrip(lengthunit, a), ustrip(angularunit, i),
                              ustrip(angularunit, Ω), ustrip(angularunit, ω), ustrip(angularunit, ν), ustrip(timeunit, epoch);
                              lengthunit = lengthunit, timeunit = timeunit, angularunit = angularunit)
    end

end

"""
Returns `Unitful` Restricted Two-body orbit eccentricity.
"""
eccentricity(kep::KeplerianState) = kep.e

"""
Returns `Unitful` Restricted Two-body orbit semimajor axis.
"""
semimajor_axis(kep::KeplerianState) = kep.a * lengthunit(kep)

"""
Returns `Unitful` Restricted Two-body orbit inclination.
"""
inclination(kep::KeplerianState) = kep.i * angularunit(kep)

"""
Returns `Unitful` Restricted Two-body orbit right ascension of the ascending node (RAAN).
"""
RAAN(kep::KeplerianState) = kep.Ω * angularunit(kep)

"""
Returns `Unitful` Restricted Two-body orbit argument of periapsis.
"""
argument_of_periapsis(kep::KeplerianState) = kep.ω * angularunit(kep)

"""
Returns `Unitful` Restricted Two-body orbit true anomoly.
"""
true_anomoly(kep::KeplerianState) = kep.ν * angularunit(kep)

"""
Returns the timepoint associated with a `KeplerianState`.
"""
epoch(kep::KeplerianState) = kep.t * timeunit(kep)

"""
Returns the dimmensionless angular unit associated with the Keplerian state.
"""
angularunit(::KeplerianState{F, LU, TU, AU}) where {F, LU, TU, AU} = AU()

"""
An orbital state within the Restricted Two-body Problem.
"""
const R2BPState{F,LU,TU} = Union{CartesianState{F, LU, TU, <:Union{Inertial, BarycentricInertial, ECI, HCI}}, KeplerianState{F, LU, TU, <:Unitful.DimensionlessUnits}} where {F,LU,TU}

"""
A structure which contains __all__ relevant values for Restricted Two-body Problem orbits.
"""
struct RestrictedTwoBodyOrbit{C<:AbstractConic, F, LU, TU, T<:R2BPState{F,LU,TU}} <: AbstractOrbit{F, LU, TU}
    state::T
    system::RestrictedTwoBodySystem{F,LU,TU}

    function RestrictedTwoBodyOrbit(r, v, body::RestrictedTwoBodySystem{T}, epoch=0) where {T}
        F = promote_type(eltype(ustrip.(r)), eltype(ustrip.(v)), T, typeof(ustrip(epoch)))
        if !(F <: AbstractFloat)
            @warn "Promoted type $(string(F)) is not of type float. Defaulting to Float64."
            F = Float64
        end
        state = CartesianState(r, v, epoch, Inertial)
        LU = lengthunit(state)
        TU   = timeunit(state)
        newbody = RestrictedTwoBodySystem(F(ustrip(LU^3 / TU^2, mass_parameter(body))), body.name; lengthunit = LU, timeunit = TU)
        return new{conic(eccentricity(state.r, state.v, newbody.μ)), F, typeof(LU), typeof(TU), typeof(state)}(state, newbody)
    end
    RestrictedTwoBodyOrbit(r, v, μ::Number, epoch=0) = RestrictedTwoBodyOrbit(r, v, RestrictedTwoBodySystem(μ), epoch)

    function RestrictedTwoBodyOrbit(e, a, i, Ω, ω, ν, body::RestrictedTwoBodySystem{T}, epoch=0) where {T}
        state = KeplerianState(e, a, i, Ω, ω, ν, epoch)
        F = typeof(state).parameters[1]
        LU = lengthunit(state)
        TU   = timeunit(state)
        newbody = RestrictedTwoBodySystem(F(ustrip(LU^3 / TU^2, mass_parameter(body))), body.name; lengthunit = LU, timeunit = TU)
        return new{conic(state.e), F, typeof(LU), typeof(TU), typeof(state)}(state, newbody)
    end
    RestrictedTwoBodyOrbit(e, a, i, Ω, ω, ν, μ::Number, epoch=0) = RestrictedTwoBodyOrbit(e, a, i, Ω, ω, ν, RestrictedTwoBodySystem(μ), epoch)

    function RestrictedTwoBodyOrbit(state::CartesianState, sys::RestrictedTwoBodySystem)
        F = promote_type(eltype(state), eltype(sys))
        LU = lengthunit(state)
        TU = timeunit(state)
        e = eccentricity(position_vector(state), velocity_vector(state), mass_parameter(sys))
        newstate = convert(CartesianState{F, typeof(LU), typeof(TU)}, state)
        newsys   = convert(RestrictedTwoBodySystem{F, typeof(LU), typeof(TU)}, sys)
        return new{conic(e), F, typeof(LU), typeof(TU), typeof(newstate)}(newstate, newsys)
    end

    function RestrictedTwoBodyOrbit(state::KeplerianState, sys::RestrictedTwoBodySystem)
        F = promote_type(eltype(state), eltype(sys))
        LU = lengthunit(state)
        TU = timeunit(state)
        AU = angularunit(state)
        e = eccentricity(state)
        newstate = convert(KeplerianState{F, LU, TU, AU}, state)
        newsys   = convert(RestrictedTwoBodySystem{F, LU, TU}, sys)
        return new{conic(e), F, typeof(LU), typeof(TU), typeof(newstate)}(newstate, newsys)
    end

end

"""
Convert between `eltype`, `lengthunit`, and `timeunit` values for `KeplerianStates`.
"""
function Base.convert(::Type{KeplerianState{F,LU,TU,AU}}, kep::KeplerianState) where {F,LU,TU,AU}
    e = F(eccentricity(kep))
    a = F(ustrip(LU(), semimajor_axis(kep)))
    i = F(ustrip(AU(), inclination(kep)))
    Ω = F(ustrip(AU(), RAAN(kep)))
    ω = F(ustrip(AU(), argument_of_periapsis(kep)))
    ν = F(ustrip(AU(), true_anomoly(kep)))
    t = F(ustrip(TU(), epoch(kep)))
    return KeplerianState(e, a, i, Ω, ω, ν, t; lengthunit = LU(), timeunit = TU(), angularunit = AU())
end

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.LengthFreeUnits)(state::KeplerianState{F,LU,TU}) where {F,LU,TU} = convert(KeplerianState{F, typeof(u), TU}, state)

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.TimeFreeUnits)(state::KeplerianState{F,LU,TU}) where {F,LU,TU} = convert(KeplerianState{F, LU, typeof(u)}, state)

"""
`RestrictedTwoBodyOrbit` instanes are likely the most commonly used feature
of this package. To accomadate this, `Orbit` is an alias for common
`RestrictedTwoBodyOrbit` constructers.
"""
Orbit(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, body::RestrictedTwoBodySystem, t = 0) = CartesianOrbit(r, v, body, t)

"""
`RestrictedTwoBodyOrbit` instanes are likely the most commonly used feature
of this package. To accomadate this, `Orbit` is an alias for common
`RestrictedTwoBodyOrbit` constructers.
"""
Orbit(r::AbstractVector{<:Unitful.Length}, v::AbstractVector{<:Unitful.Velocity}, body::RestrictedTwoBodySystem, t = 0u"s") = CartesianOrbit(r, v, body, t)

"""
`RestrictedTwoBodyOrbit` instanes are likely the most commonly used feature
of this package. To accomadate this, `Orbit` is an alias for common
`RestrictedTwoBodyOrbit` constructers.
"""
Orbit(e::Real, a::Real, i::Real, Ω::Real, ω::Real, ν::Real, body::RestrictedTwoBodySystem, t = 0) = KeplerianOrbit(e, a, i, Ω, ω, ν, body, t)

"""
`RestrictedTwoBodyOrbit` instanes are likely the most commonly used feature
of this package. To accomadate this, `Orbit` is an alias for common
`RestrictedTwoBodyOrbit` constructers.
"""
Orbit(e::Real, a::Unitful.Length, i::Unitful.DimensionlessQuantity, Ω::Unitful.DimensionlessQuantity, ω::Unitful.DimensionlessQuantity, ν::Unitful.DimensionlessQuantity, body::RestrictedTwoBodySystem, t = 0u"s") = KeplerianOrbit(e, a, i, Ω, ω, ν, body, t)

"""
Print `KeplerianState` instances to `io`.
"""
function Base.show(io::IO, orbit::KeplerianState{F,LU,TU,AU}) where {F,LU,TU,AU}

    println(io, "  Keplerian State:")
    println(io, "")
    println(io, "    e:  ",
                orbit.e)
    println(io, "    a:  ",
                orbit.a, " ", string(LU()))
    println(io, "    i:  ",
                orbit.i, AU == u"rad" ? " " : "", string(AU()))
    println(io, "    Ω:  ",
                orbit.Ω, AU == u"rad" ? " " : "", string(AU()))
    println(io, "    ω:  ",
                orbit.ω, AU == u"rad" ? " " : "", string(AU()))
    println(io, "    ν:  ",
                orbit.ν, AU == u"rad" ? " " : "", string(AU()))

end

"""
Print `RestrictedTwoBodyOrbit` instances to `io`.
"""
function Base.show(io::IO, orbit::RestrictedTwoBodyOrbit{C, F, T}) where {C, F, T}
    print(io, string(C), " Restricted Two-body Orbit")
    println(io, " (", string(F), "):")
    println(io, "")
    show(io, orbit.state)
    println(io, "")
    show(io, orbit.system)
end

"""
An alias for `RestrictedTwoBodyOrbit`.
"""
const R2BPOrbit = RestrictedTwoBodyOrbit

"""
An alias for `RestrictedTwoBodySystem`.
"""
const R2BPSystem = RestrictedTwoBodySystem

"""
An alias for `RestrictedTwoBodyOrbit` instances with `KeplerianState` values.
"""
const KeplerianOrbit{C,F,LU,TU} = RestrictedTwoBodyOrbit{C, F, LU, TU, K} where {K <: KeplerianState{F, LU, TU, AU} where AU}

"""
An alias for `RestrictedTwoBodyOrbit` instances with `CartesianState` values.
"""
const CartesianOrbit{C,F,LU,TU} = RestrictedTwoBodyOrbit{C, F, LU, TU, R} where {R <: CartesianState{F, LU, TU, <:Union{Inertial, BarycentricInertial, ECI, HCI}}}

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `KeplerianState` values.
"""
KeplerianOrbit(e, a, i, Ω, ω, ν, system, epoch) = RestrictedTwoBodyOrbit(e, a, i, Ω, ω, ν, system, epoch)

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `CartesianState` values.
"""
CartesianOrbit(r, v, system, epoch) = RestrictedTwoBodyOrbit(r, v, system, epoch)

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `KeplerianState` values.
"""
KeplerianOrbit(orbit::CartesianOrbit) = KeplerianOrbit(keplerian(position_vector(orbit), velocity_vector(orbit), orbit.system)..., orbit.system, epoch(orbit.state))

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `CartesianState` values.
"""
CartesianOrbit(orbit::KeplerianOrbit) = CartesianOrbit(cartesian(eccentricity(orbit), semimajor_axis(orbit), inclination(orbit), RAAN(orbit), argument_of_periapsis(orbit), true_anomoly(orbit), mass_parameter(orbit.system))..., orbit.system, epoch(orbit.state))

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `KeplerianState` values.
"""
KeplerianOrbit(state::KeplerianState, system) = RestrictedTwoBodyOrbit(state, system)

"""
Alias for a `RestrictedTwoBodyOrbit` constructor with `CartesianState` values.
"""
CartesianOrbit(state::CartesianOrbit, system) = RestrictedTwoBodyOrbit(state, system)

"""
Prints a `Trajectory` instance to `io`.
"""
Base.show(io::IO, traj::Trajectory{<:CartesianOrbit}) = println(io, "Cartesian Two-body trajectory with ", length(traj), " steps")

"""
Prints a `Trajectory` instance to `io`.
"""
Base.show(io::IO, traj::Trajectory{<:KeplerianOrbit}) = println(io, "Keplerian Two-body trajectory with ", length(traj), " steps")

Base.string(::Type{Circular}) = "Circular"
Base.string(::Type{Elliptical}) = "Elliptical"
Base.string(::Type{Hyperbolic}) = "Hyperbolic"
Base.string(::Type{Parabolic}) = "Parabolic"
Base.string(::Type{Invalid}) = "Invalid"


text:
type DataSet <: FeatureSpace
  clusters::Dict{Any, Cluster}
  vector_sum::FeatureVector
  mdata::Any
  DataSet() = new(Dict{Any,Cluster}(),FeatureVector())
end

# returns the Cluster indexed by [key]
function getindex(ds::DataSet, key)
  return ds.clusters[key]
end

# maps a [key] to a Cluster [c]
function setindex!(ds::DataSet, c::Cluster, key)
  if haskey(ds, key)
    subtract!(ds.vector_sum,ds[key].vector_sum)
  end
  add!(ds.vector_sum,c.vector_sum)
  ds.clusters[key] = c
end

function length(ds::DataSet)
  size = 0
  for clust in ds.clusters
    size += length(clust[2])
  end
  return size
end

# returns all keys in the DataSet
function keys(ds::DataSet)
  return Base.keys(ds.clusters)
end

# returns all Clusters in the DataSet
function values(ds::DataSet)
  return Base.values(ds.clusters)
end

# check to see if the DataSet is empty.
function isempty(ds::DataSet)
  return Base.isempty(ds.clusters)
end

function haskey(ds::DataSet, key)
  return Base.haskey(ds.clusters, key)
end

# prints out the dataset cleanly
function display(ds::DataSet)
  print_with_color(:white,string(typeof(ds)))
  k = length(ds.clusters)

  print_with_color(:white," with $k Clusters")
  if k > 0
    print_with_color(:white,":\n")
    if k > 5
      k = 5
    end

    for key in collect(keys(ds))[1:k]
      l = length(ds[key])
      if length(key) > 40
        key = key[1:40]
        key = key*"…"
      end
      print_with_color(:white,"    $key with $l FeatureVectors\n")
    end
    if length(ds) > 5
      print_with_color(:white,"\t\t\t⋮\n")
    end
  end
end

function show(io::IO, ds::DataSet)
  print(io,string(typeof(ds)))
  k = length(ds.clusters)
  print(io," with $k Clusters")
end


text:
using Plots: display
using Serialization

demo_config = TrainConfig(
    # metadata
    name = "demo",
    data_path = "data/demo/demo_data.jls",
    # train configuration
    num_epochs = 200,
    p_init_fn = initial_params,
    patience = 10,
    logevery = 3,
    # model configuration
    ode = CTRNN1Canonical(3, 6, 2; σ_u = identity, σ = relu),
    link = AffineCanonical(3, 3),
    mod = ISSCTRNN1(3, 6; enable = true, is_fixed = false, ϵ = 1e-3),
    feed = NoFeedthrough(),
    # DynamicsTrainer arguments
    solver = BS3(),
    adaptive = true,
    sensealg = InterpolatingAdjoint(checkpointing = true, autojacvec = ZygoteVJP()),
    dt = deserialize("data/demo/demo_metadata.jls")["dt"],
    reltol = 1e-4,
    abstol = 1e-6,
    initial_condition = solve_ss,
    loss = Flux.Losses.mse,
    optim = ADAMW(1e-3),
    wdecay = WeightDecay(0.0),
    sched = Exp(λ = 5e0, γ = 1.0),
    validation_split = 0.16,
    test_split = 0.0,
    train_split = 0.84,
    batch_size = 8,
    time_dropout = 0.5,
    time_horizon = 500,
    time_horizon_increment = 50,
    time_horizon_max = 500,
    callback_tasks = [
        generate_netlist_task(
            model_name = "demo",
            num_ports = 2,
            input_config = [("", "V", 1), ("", "V", 2)],
            output_config = [("", "I", 1), ("", "I", 2)],
            time_scale = dataset =
                deserialize("data/demo/demo_metadata.jls")["t_scale_factor"],
            input_scale = deserialize("data/demo/demo_metadata.jls")["u_scale"],
            input_bias = deserialize("data/demo/demo_metadata.jls")["u_bias"],
            output_scale = deserialize("data/demo/demo_metadata.jls")["y_scale"],
            output_bias = deserialize("data/demo/demo_metadata.jls")["y_bias"],
        ),
    ],
)

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
abstract type EntryType end
struct KeepEntry <: EntryType end
struct DropEntry <: EntryType end
struct NewEntry <: EntryType end

struct DefaultBranch end

struct Package
    name::String
    uuid::UUIDs.UUID
end

mutable struct DepInfo
    package::Package
    latest_version::Union{VersionNumber,Nothing}
    version_spec::Union{VersionSpec,Nothing}
    version_verbatim::Union{String,Nothing}

    function DepInfo(
        package::Package;
        latest_version::Union{VersionNumber,Nothing}=nothing,
        version_spec::Union{VersionSpec,Nothing}=nothing,
        version_verbatim::Union{String,Nothing}=nothing,
    )
        return new(package, latest_version, version_spec, version_verbatim)
    end
end

function Base.in(p::Package, s::Set{DepInfo})
    for i in s
        if i.package == p
            return true
        end
    end

    return false
end

const DEFAULT_REGISTRIES = Pkg.RegistrySpec[Pkg.RegistrySpec(;
    name="General",
    uuid="23338594-aafe-5451-b93e-139f81909106",
    url="https://github.com/JuliaRegistries/General.git",
)]

Base.@kwdef struct Options
    entry_type::EntryType = KeepEntry()
    registries::Vector{Pkg.RegistrySpec} = DEFAULT_REGISTRIES
    use_existing_registries::Bool = false
    depot::String = DEPOT_PATH[1]
    subdirs::Vector{String} = [""]
    master_branch::Union{DefaultBranch,AbstractString} = DefaultBranch()
    bump_compat_containing_equality_specifier::Bool = true
    pr_title_prefix::String = ""
    include_jll::Bool = false
    unsub_from_prs::Bool = false
    cc_user::Bool = false
    bump_version::Bool = false
    include_yanked::Bool = false
end


text:
# Randomization option slows down processing
# but improves quality of solution considerably
# Would be better to do randomization in place
function cddual{T<:Real}(X::AbstractMatrix{T},
                         Y::AbstractVector{T};
                         C::Real = 1.0,
                         norm::Integer = 2,
                         randomized::Bool = true,
                         maxpasses::Integer = 2)
    # l: # of samples
    # n: # of features
    n, l = size(X)
    alpha = zeros(Float64, l)
    w = zeros(Float64, n)

    # Set U and D
    #  * L1-SVM: U = C, D[i] = 0
    #  * L2-SVM: U = Inf, D[i] = 1 / (2C)
    U = 0.0
    if norm == 1
        U = C
        D = zeros(Float64, l)
    elseif norm == 2
        U = Inf
        D = fill(1.0 / (2.0 * C), l)
    else
        throw(ArgumentError("Only L1-SVM and L2-SVM are supported"))
    end

    # Set Qbar
    Qbar = Array{Float64,1}(l)
    for i in 1:l
        Qbar[i] = D[i] + dot(X[:, i], X[:, i])
    end

    # Loop over examples
    converged = false
    pass = 0

    while !converged
        # Assess convergence
        pass += 1
        if pass == maxpasses
            converged = true
        end

        # Choose order of observations to process
        if randomized
            indices = randperm(l)
        else
            indices = 1:l
        end

        # Process all observations
        for i in indices
            g = Y[i] * dot(w, X[:, i]) - 1.0 + D[i] * alpha[i]

            if alpha[i] == 0.0
                pg = min(g, 0.0)
            elseif alpha[i] == U
                pg = max(g, 0.0)
            else
                pg = g
            end

            if abs(pg) > 0.0
                alphabar = alpha[i]
                alpha[i] = min(max(alpha[i] - g / Qbar[i], 0.0), U)
                for j in 1:n
                    w[j] = w[j] + (alpha[i] - alphabar) * Y[i] * X[j, i]
                end
            end
        end
    end

    return SVMFit(w, pass, converged)
end


text:
################################################################################
#                         Euler Angles <=> Quaternions
################################################################################

for rot_seq in rot_seq_array
    for k = 1:samples
        # Sample a vector.
        v = randn(3)

        # Sample three angles form a uniform distribution [-pi,pi].
        θx       = -pi + 2*pi*rand()
        θy       = -pi + 2*pi*rand()
        θz       = -pi + 2*pi*rand()
        eulerang = EulerAngles(θx, θy, θz, rot_seq)

        # Rotate the vector using a DCM (which was already tested).
        D      = angle2rot(eulerang)
        rv_dcm = D*v

        # Rotate the vector using a quaternion.
        q1    = angle2rot(Quaternion,eulerang)
        rv_q1 = vect(inv(q1)*v*q1)

        q2    = angle2rot(Quaternion,θx,θy,θz,rot_seq)
        rv_q2 = vect(inv(q2)*v*q2)

        # Compare.
        @test norm(rv_q1-rv_dcm) < 1e-10
        @test norm(rv_q2-rv_dcm) < 1e-10
        @test rv_q1 ≈ rv_q2

        # Compute the Euler angles using the Quaternion.
        eulerang_quat = quat2angle(q1, eulerang.rot_seq)

        # Create the quaternion using those Euler angles and compare to the
        # original.
        q3 = angle2rot(Quaternion,eulerang_quat)

        @test q3.q0 ≈ q1.q0 atol=1e-10
        @test q3.q1 ≈ q1.q1 atol=1e-10
        @test q3.q2 ≈ q1.q2 atol=1e-10
        @test q3.q3 ≈ q1.q3 atol=1e-10
    end
end

for k = 1:samples
    # Sample three angles form a uniform distribution [-0.0001,0.0001].
    eulerang = EulerAngles(-0.0001 + 0.0002*rand(),
                           -0.0001 + 0.0002*rand(),
                           -0.0001 + 0.0002*rand(),
                           :XYZ)

    # Get the error between the exact rotation and the small angle
    # approximation.
    error = angle2rot(Quaternion,eulerang) -
            smallangle2rot(Quaternion,eulerang.a1, eulerang.a2, eulerang.a3)

    # If everything is fine, the norm of the matrix error should be small.
    @test norm(error) < 1e-7
end


text:
using DSGE, ModelConstructors, Test
using HDF5, Random
path = dirname(@__FILE__)
writing_output = false
Random.seed!(1793)

## Assess post-estimation calculation of mode
#-----------------------------------------------------------------
# csminwel
#-----------------------------------------------------------------
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Load data
file = "$path/../reference/optimize.h5"
x0   = h5read(file, "params")
data = h5read(file, "data")'

# For regenerating test file
params_test = deepcopy(x0)
data_test   = Matrix{Float64}(data)

minimizer  = h5read(file, "minimizer")
minimum    = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
DSGE.update!(m, x0)
n_iterations = 3
m <= Setting(:optimization_iterations, n_iterations)
m <= Setting(:optimization_ftol, 1e-14)
m <= Setting(:optimization_xtol, 1e-32)
m <= Setting(:optimization_gtol, 1e-8)
m <= Setting(:optimization_attempts, 0)
m <= Setting(:optimization_step_size, .01)

modal_minimizer, modal_out, modal_H, modal_hessian =
    DSGE.calculate_mode(m, data_test, vec(params_test), :csminwel, mle = false,
                        verbose = :none, get_all_results = true,
                        save_results = false)

file = "$path/../reference/estimate_util.h5"
exp_hessian = h5read(file, "hessian")

@testset "Check optimize minimizers are the same [csminwel]" begin
    @test minimizer ≈ modal_out.minimizer atol=5e-4
    @test @test_matrix_approx_eq H_expected modal_H
    # this works when ran in REPL but breaks in Test mode in Julia 1.1, 1.3, and 1.5
    @test_broken @test_matrix_approx_eq exp_hessian modal_hessian
end


text:
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using JSON
using CMSSPs
using MCTS
using Logging
using TOML

struct MCTSParams
    depth::Int64
    exploration::Float64
    init_N::Int64
    n_iters::Int64
end

function parse_mcts_params(filename::String)
    params_key = TOML.parsefile(filename)
    return MCTSParams(params_key["DEPTH"], params_key["C"], params_key["INITN"], params_key["NITERS"])
end

global_logger(SimpleLogger(stderr, Logging.Warn))

rng = MersenneTwister(3456)

param_files = ["./paramsets/scale-1.toml","./paramsets/simtime-1.toml","./paramsets/cost-1.toml"]

scale_file = param_files[1]
simtime_file = param_files[2]
cost_file = param_files[3]

# Load flight policy and create an alias for estimate value
infhor_fn = "./policies/dreamr-uf-params111.jld2"
flight_policy = load_localapproxvi_policy_from_jld2(infhor_fn)
estimate_value(dmcts::DREAMRMCTSType, state::DREAMRMCTSState, depth::Int64) = estimate_value_dreamr(flight_policy, dmcts, state, depth)

init_Q(dmcts::DREAMRMCTSType, state::DREAMRMCTSState, a::DREAMRActionType) = init_q_dreamr(flight_policy, dmcts, state, a)

# ep_file_prefix = "/scratch/shushman/HitchhikingDrones/set-2-hard/set-2-100-to-1000"
mcts_param_file = ARGS[1]
outfn = ARGS[2]
num_eps = parse(Int64, ARGS[3])
ep_file_prefix = ARGS[4]
log_out_file_prefix = ARGS[5]
to_log = parse(Bool,ARGS[6])

# ep_file_prefix = ARGS[1]
# num_eps = parse(Int64, ARGS[2])
# mcts_param_file = ARGS[3]
# outfn = ARGS[4]

params = parse_params(scale_file=scale_file, simtime_file=simtime_file, cost_file=cost_file)
mcts_params = parse_mcts_params(mcts_param_file)

costs = Vector{Float64}(undef, 0)
steps = Vector{Int64}(undef, 0)
mode_switches = Vector{Int64}(undef, 0)
unsuccessful = Vector{Int64}(undef, 0)

for iter=1:num_eps

    episode_dict = Dict

    ep_file = string(ep_file_prefix,"-",iter,".json")

    @show ep_file

    open(ep_file,"r") do f
        episode_dict = JSON.parse(f)
    end

    context_set = get_dreamr_episode_context(episode_dict["epochs"])

    start_pos = Point(episode_dict["start_pos"][1], episode_dict["start_pos"][2])
    goal_pos = Point(episode_dict["goal_pos"][1], episode_dict["goal_pos"][2])

    log_soln_dict = Dict()
    if to_log
        log_fn = string(log_out_file_prefix,"-",iter,"-output.json")
        log_soln_dict["start_pos"] = episode_dict["start_pos"]
        log_soln_dict["goal_pos"] = episode_dict["goal_pos"]
        log_soln_dict["epochs"] = Dict(0 => Dict("drone-info"=>Dict("pos"=>[start_pos.x,start_pos.y], "on_car"=>"")))
    end

    start_state = DREAMRStateType{MultiRotorUAVState}(FLIGHT, get_state_at_rest(MultiRotorUAVState, start_pos))

    # @show start_state
    # @show goal_pos
    # readline()

    # Create CMSSP
    cmssp = create_dreamr_cmssp(MultiRotorUAVState, MultiRotorUAVAction, context_set, goal_pos, params)

    dmcts = DREAMRMCTSType(cmssp)

    # Create solver and planner
    solver = DPWSolver(depth=mcts_params.depth,
                       exploration_constant = mcts_params.exploration,
                       init_N = mcts_params.init_N,
                       n_iterations = mcts_params.n_iters,
                       estimate_value=estimate_value,
                       # init_Q = init_Q,
                       rng=rng)
    planner = solve(solver, dmcts)

    curr_dmcts_state = DREAMRMCTSState(start_state)

    timesteps = 0
    cost = 0
    num_mode_switches = 0
    successful = false

    while true

        # @debug curr_dmcts_state
        a = action(planner, curr_dmcts_state)
        # @debug a

        # readline()

        (next_state, reward, _, timeout) = simulate_cmssp!(dmcts.cmssp, curr_dmcts_state.cmssp_state, a, timesteps, rng)

        if next_state.mode != curr_dmcts_state.cmssp_state.mode
            num_mode_switches += 1
        end

        cost += -1.0*reward

        curr_dmcts_state = DREAMRMCTSState(next_state, cmssp.curr_context_set.car_id, 0)

        timesteps += 1
        # @debug timesteps

        if to_log
            log_soln_dict["epochs"][timesteps] = Dict("drone-info"=>Dict("pos"=>[next_state.continuous.x, next_state.continuous.y],
                                                  "on_car"=>cmssp.curr_context_set.car_id))
        end

        if isterminal(dmcts, curr_dmcts_state)
            successful = true
            break
        end

        if timeout
            @debug "Timeout!"
            break
        end
    end

    if successful
        energy_cost = cost - params.cost_params.TIME_COEFFICIENT*params.time_params.MDP_TIMESTEP*timesteps
        @show energy_cost
        push!(costs, energy_cost)
        push!(steps, timesteps)
        push!(mode_switches, num_mode_switches)

        log_soln_dict["num_epochs"] = timesteps+1
        log_soln_dict["success"] = false
        if to_log
            open(log_fn,"w") do f
                JSON.print(f, log_soln_dict, 2)
            end
        end
    else
        @info "unsuccessful"
        push!(unsuccessful, iter)

        log_soln_dict["num_epochs"] = episode_dict["num_epochs"]
        log_soln_dict["success"] = true
    end
end

results_dict = Dict("costs"=>costs, "steps"=>steps, "mode_switches"=>mode_switches, "failed_iters"=>unsuccessful)

open(outfn, "w") do f
    JSON.print(f, results_dict, 2)
end


text:
##########################################################################
# Timing for list of methods. loop_over works like in simulate.
function timing(network, simulation_params; loop_over::Symbol=:precoding_methods)
    println("-- timing on $network.")
    Lumberjack.info("Starting timing.", @compat Dict(:network => network, :simulation_params => simulation_params))

    # Get assignment method if needed
    assignment_method() = nothing
    if loop_over == :precoding_methods
        if haskey(simulation_params, "assignment_methods")
            if length(simulation_params["assignment_methods"]) > 1
                Lumberjack.warn("Looping over precoding methods: will only use first assignment method provided.")
            end
            assignment_method(channel, network) = simulation_params["assignment_methods"][1](channel, network)
        else
            assignment_method(channel, network) = IDCellAssignment!(channel, network)
        end
    end

    # Set initial aux params
    set_initial_aux_params!(simulation_params, network)

    # Ensure that we are not storing all intermediate iterations.
    set_aux_precoding_param!(network, :final_iteration, "output_protocol")

    # Get channel realization to work on
    draw_user_drop!(network)
    channel = draw_channel(network)

    if loop_over == :precoding_methods
        # Perform assignment
        assignment_method(channel, network)

        # Make sure things are JITed
        for method in simulation_params["precoding_methods"]
            method(channel, network)
        end

        # Run performance test
        for method in simulation_params["precoding_methods"]
            println("--- Testing performance of ", stringify_method(method))
            @time for i = 1:simulation_params["Ntest"]; method(channel, network); end
        end
    elseif loop_over == :assignment_methods
        # Make sure things are JITed
        for method in simulation_params["assignment_methods"]
            method(channel, network)
        end

        # Run performance test
        for method in simulation_params["assignment_methods"]
            println("--- Testing performance of ", stringify_method(method))
            @time for i = 1:simulation_params["Ntest"]; method(channel, network); end
        end
    end
end


text:
#transform tests

# test Parseval's theorem for wavefunctions
function parsevaltest(L,N)
    X,K,dX,dK = makearrays(L,N)
    T = makeT(X,K)

    ψ = randn(N...) + im*randn(N...)

    n1 = sum(abs2.(ψ))*prod(dX)
    ϕ = T.Txk*ψ

    n2 = sum(abs2.(ϕ))*prod(dK)
    return n1,n2
end

L = (88.,)
N = (128,)
# X,K,dX,dK = makearrays(L,N)
# T = makeT(X,K) # TODO: norm problem in makeT for 2,3 dims

n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

## debug offending test
L = (30.,20.)
N = (30,34)

# X,K,dX,dK = makearrays(L,N)
# T = makeT(X,K) # TODO: norm problem in makeT for 2,3 dims

# ψ = randn(N...) + im*randn(N...)

# n1 = sum(abs2.(ψ))*prod(dX)
# ϕ = T.Txk*ψ

# n2 = sum(abs2.(ϕ))*prod(dK)

n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

L = (30.,20.,10)
N = (30,34,24)
n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

L = (30.,20.,10,50.)
N = (30,34,24,22)
n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2


text:
# Large example based on simulated data
# Create LOSO adjusted phenotypes based on GenSel output
# Run GWAS based on these
#
# Download dataset at https://datadryad.org/
# Extract:
# 800K genotypes on 80K indivs: gen11.1.bed/bim/fam
# 31K genotypes on 80K indivs: gen12.1.bed/bim/fam
# Simulated phenotype: phen_8.1.txt
#
# These two functions will typically be run in separate chunks as they will typically have different CPU/MEM requirements
#
using MANA
# create adjusted phenotypes
function make_loso_phenotypes()
  dir = "/path/to/files"
  chr=26
  interval=10000000
  phenotypesfilename= dir .* "phen_8.1.txt"
  markersfileprefix=dir .* "gen12.1"
  alphasfilename=dir .* "other_13.txt"
  outputprefix=dir .* "adjusted_phenotype.chr26"
  REGION=0
  CODE="LOSO"
  ycreator = build_ycreator()
  add_alphas(ycreator, alphasfilename) # marker effects
  add_phenotypes(ycreator, phenotypesfilename) # phenotype
  add_markers_gensel(ycreator, markersfileprefix) # marker genotypes
  sort_ycreator(ycreator)
  initiate_ycreator(ycreator, interval, chr, outputprefix, REGION, CODE)  # set up creator object
  create_yk_loso(ycreator) # create adjusted phenotypes based on LOSO and the interval given
end
#
make_loso_phenotype()
#
# run GWAS
function main()
	dir = "/path/to/files/" # May need to be edited to directory of interest
	snpfilename = dir .* "gen11.1"
  phenotypesprefix = dir .* "adjusted_phenotype.chr26"
  outputfilename = dir .* "out.simulation_example"
  printchain = false
  weightsfilename= dir .* "weights.simulation.txt" # file of IDs and weights (default = 1)
  interval = 10000000
  setsize = 1000
  type = "dominance"
  # Build Model
  model = build_gwas()
  add_snps(model, snpfilename)
  setup_phenotypes(model, phenotypesprefix, weightsfilename)
  initiate_model(model, interval, setsize, outputfilename, printchain, type)
  # Run Model
  run_model(model)
  # Summarise
  concatenate_summaries(model)
end
#
main()

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test Base.Cartesian.exprresolve(:(1 + 3)) == 4
ex = Base.Cartesian.exprresolve(:(if 5 > 4; :x; else :y; end))
@test ex.args[2] == QuoteNode(:x)

@test Base.Cartesian.lreplace!("val_col", Base.Cartesian.LReplace{String}(:col, "col", 1)) == "val_1"
@test Base.setindex(CartesianIndex(1,5,4),3,2) == CartesianIndex(1, 3, 4)

# test conversions for CartesianIndex

@testset "CartesianIndex Conversions" begin
    @test convert(Int, CartesianIndex(42)) === 42
    @test convert(Float64, CartesianIndex(42)) === 42.0
    @test convert(Tuple, CartesianIndex(42, 1)) === (42, 1)

    # can't convert higher-dimensional indices to Int
    @test_throws MethodError convert(Int, CartesianIndex(42, 1))
end

@testset "CartesianIndices overflow" begin
    I = CartesianIndices((1:typemax(Int),))
    i = last(I)
    @test iterate(I, i) === nothing

    I = CartesianIndices((1:(typemax(Int)-1),))
    i = CartesianIndex(typemax(Int))
    @test iterate(I, i) === nothing

    I = CartesianIndices((1:typemax(Int), 1:typemax(Int)))
    i = last(I)
    @test iterate(I, i) === nothing

    i = CartesianIndex(typemax(Int), 1)
    @test iterate(I, i) === (CartesianIndex(1, 2), CartesianIndex(1,2))

    # reverse cartesian indices
    I = CartesianIndices((typemin(Int):(typemin(Int)+3),))
    i = last(I)
    @test iterate(I, i) === nothing
end

@testset "CartesianIndices iteration" begin
    I = CartesianIndices((2:4, 0:1, 1:1, 3:5))
    indices = Vector{eltype(I)}()
    for i in I
        push!(indices, i)
    end
    @test length(I) == length(indices)
    @test vec(I) == indices

    empty!(indices)
    I = Iterators.reverse(I)
    for i in I
        push!(indices, i)
    end
    @test length(I) == length(indices)
    @test vec(collect(I)) == indices

    # test invalid state
    I = CartesianIndices((2:4, 3:5))
    @test iterate(I, CartesianIndex(typemax(Int), 3))[1] == CartesianIndex(2,4)
    @test iterate(I, CartesianIndex(typemax(Int), 4))[1] == CartesianIndex(2,5)
    @test iterate(I, CartesianIndex(typemax(Int), 5))    === nothing

    @test iterate(I, CartesianIndex(3, typemax(Int)))[1] == CartesianIndex(4,typemax(Int))
    @test iterate(I, CartesianIndex(4, typemax(Int)))    === nothing
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module Models

export Domain, Physics, Controls, Economics, ClimateModelParameters, ClimateModel

include("domain.jl")
include("physics.jl")
include("controls.jl")
include("economics.jl")

"""
    ClimateModelParameters(name, domain::Domain, economics::Economics, physics::Physics)

Create a named instance of the MARGO climate model parameters, which include
economic input parameters (`economics`), physical climate parameters (`physics`),
and climate control policies (`controls`) on some spatial-temporal grid (`domain`).

Use these to construct a [`ClimateModel`](@ref), which also contains the optimized
controls.
"""
mutable struct ClimateModelParameters
    name::String
    domain::Domain
    economics::Economics
    physics::Physics
end

mutable struct ClimateModel
    name::String
    domain::Domain
    economics::Economics
    physics::Physics
    controls::Controls
end

"""
    ClimateModel(params::ClimateModelParameters[, controls::Controls])

Create an instance of an extremely idealized multi-control climate model. The
returned object contains the [`ClimateModelParameters`](@ref), and will contain
the optimized [`Controls`](@ref). These can be computed using
[`optimize_controls!`](@ref).
"""
ClimateModel(params::ClimateModelParameters, controls::Controls) = ClimateModel(
    params.name,
    params.domain,
    params.economics,
    params.physics,
    controls
)
function ClimateModel(params::ClimateModelParameters)
    dom = params.domain
    t = collect(dom.initial_year:dom.dt:dom.final_year)
    return ClimateModel(
        params,
        Controls(
            zeros(size(t)),
            zeros(size(t)),
            zeros(size(t)),
            zeros(size(t))
            )
        )
end

end


text:
int main() {
    if (false);
    return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler