text:
module InitialGuessModule
export computeDensityGuessSAD, InitialGuess, ZeroGuess, SADGuess
using ..AtomModule

abstract type InitialGuess end
type ZeroGuess_ <: InitialGuess end
ZeroGuess = ZeroGuess_()
type SADGuess_ <: InitialGuess end
SADGuess = SADGuess_()

function computeDensityMatrixGuess(mat::Matrix,dim::Int64)
  if (size(mat) != (dim,dim))
    error("Guess matrix has wrong dimensions: $(size(mat)) (expected $((dim,dim)))")
  end
  return mat
end

function computeDensityMatrixGuess(::ZeroGuess_,dim::Int64)
  return zeros(dim,dim)
end

type DensityGuessSAD
  α
  β
end

# These are constructed by running a unrestricted SCF calculation on the single atom with the corresponding method and basis set
densityGuessSADHartreeFock = DensityGuessSAD(
#alpha spin density
Dict("STO-3G" => Dict(
	Element("H") => [ 1.000000 ],
	Element("O") => [ 1.065739   -0.264689    0.000000    0.000000    0.000000;
                         -0.264689    1.065739    0.000000    0.000000    0.000000;
                          0.000000    0.000000    1.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    1.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000 ]
	),
     "def2-SVP" => Dict(
        Element("H") => [ 0.345487    0.294542    0.000000    0.000000    0.000000;
                          0.294542    0.251108    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000 ],
	Element("O") => [ 1.083840    0.188236    0.143345    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.188236    0.470190    0.284077    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.143345    0.284077    0.172502    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.500580    0.000000    0.000000    0.308265   -0.000000   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.473992    0.000000    0.000000    0.315137    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.500580    0.000000   -0.000000    0.308265    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.308265    0.000000    0.000000    0.189835    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000   -0.000000    0.315137   -0.000000    0.000000    0.209522   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000   -0.000000    0.000000    0.308265    0.000000   -0.000000    0.189835    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 ]
	)
),
#beta spin density
Dict("STO-3G" => Dict(
	Element("H") => [ 0.000000 ],
	Element("O") => [ 1.065739   -0.264689    0.000000    0.000000    0.000000;
                         -0.264689    1.065739    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000 ]
	),
     "def2-SVP" => Dict(
        Element("H") => [ 0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000 ],
	Element("O") => [ 1.075394    0.156592    0.157593    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.156592    0.367558    0.298568    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.157593    0.298568    0.243442    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000   -0.000000   -0.000000    0.000000   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000   -0.000000    0.425423    0.000000   -0.000000    0.324962   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000   -0.000000    0.000000    0.000000   -0.000000    0.000000   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000   -0.000000   -0.000000    0.000000   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000   -0.000000    0.324962    0.000000   -0.000000    0.248224   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000   -0.000000   -0.000000    0.000000   -0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000;
                          0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000 ]
	)
)
)

function computeDensityGuessSAD(methodName, basisSetName, geometry)
  if methodName=="HF" || methodName=="HartreeFock"
    αDensities = densityGuessSADHartreeFock.α[basisSetName]
    βDensities = densityGuessSADHartreeFock.β[basisSetName]
    αResult = cat([1,2],[αDensities[atom.element] for atom in geometry.atoms]...)
    βResult = cat([1,2],[βDensities[atom.element] for atom in geometry.atoms]...)
    return (αResult, βResult)
  else
    error("No SAD guess implemented for method "*methodName)
  end
end
end #module


text:
module TestExamples2DEuler

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler" begin
  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5])
  end

  @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fv.jl"),
      l2   = [0.026440292358506527, 0.013245905852168414, 0.013245905852168479, 0.03912520302609374],
      linf = [0.042130817806361964, 0.022685499230187034, 0.022685499230187922, 0.06999771202145322])
  end

  @trixi_testset "elixir_euler_density_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0010600778457964775, 0.00010600778457634275, 0.00021201556915872665, 2.650194614399671e-5],
      linf = [0.006614198043413566, 0.0006614198043973507, 0.001322839608837334, 0.000165354951256802],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.3653424742011065e-6, 2.1388875092652046e-6, 2.1388875093792834e-6, 6.010896863407163e-6],
      linf = [1.4080465934984687e-5, 1.7579850582816192e-5, 1.757985059525069e-5, 5.95689353266593e-5])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
      linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
  end

  @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481471610306124, 0.027694280613944234, 0.027697905866996532, 0.12932052501462554],
      linf = [0.31052098400669004, 0.3481295959664616, 0.34807152194137336, 1.1044947556170719],
      maxiters = 10,
      surface_flux = flux_kennedy_gruber,
      volume_flux = flux_kennedy_gruber)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481122603050542, 0.027662840593087695, 0.027665658732350273, 0.12927455860656786],
      linf = [0.3110089578739834, 0.34888111987218107, 0.3488278669826813, 1.1056349046774305],
      maxiters = 10,
      surface_flux = flux_chandrashekar,
      volume_flux = flux_chandrashekar)
  end

  @trixi_testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.05380378209319121, 0.04696482040740498, 0.046967440141927635, 0.19686385916646665],
      linf = [0.1852693811520424, 0.24028641582658528, 0.2326441454434102, 0.6874069967414047])
  end

  @trixi_testset "elixir_euler_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
      l2   = [0.1399704878290245, 0.11580885255188444, 0.115808885313645, 0.3391969791404581],
      linf = [1.4473441276982646, 1.3355435317505784, 1.335543531751034, 1.8203831482774346],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"),
      l2   = [0.470793846459279, 0.20772743006225688, 0.20772649286847356, 0.7011946441199758],
      linf = [1.459980543578883, 0.976674580177999, 0.9759214032734075, 2.960406757135902],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
      l2   = [0.4667685563568254, 0.20611874363505156, 0.20588110990842964, 0.6993986008683889],
      linf = [1.4753582207406488, 0.9904209156187197, 0.9841399353882521, 2.9346283218997806],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl with mortars" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
      l2   = [0.016398689755743465, 0.03074685501652015, 0.03078660152836722, 0.15176702545248866],
      linf = [0.27612812614937343, 0.5598351620373528, 0.5601415330040513, 2.48372589562021],
      refinement_patches=(
          (type="box", coordinates_min=(-0.25, -0.25), coordinates_max=(0.25, 0.25)),
          (type="box", coordinates_min=(-0.125, -0.125), coordinates_max=(0.125, 0.125)),),
      initial_refinement_level = 4,
      maxiters = 5)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_cnn.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_cnn.jl"),
      l2   = [0.47425390570612547, 0.20990014031106777, 0.2103622732820483, 0.7022475294762982],
      linf = [1.6366874508175182, 0.9734578279971777, 0.9412921347930632, 2.949655726363399],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_pure_fv.jl"),
      l2   = [0.39957047631960346, 0.21006912294983154, 0.21006903549932, 0.6280328163981136],
      linf = [2.20417889887697, 1.5487238480003327, 1.5486788679247812, 2.4656795949035857],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
      l2   = [0.6777351803934385, 0.28143110338521693, 0.28143035507881886, 0.7211949078634373],
      linf = [2.881942488461436, 1.8037294505391677, 1.803677581960154, 3.037915614619938],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.4819231474652106, 0.1655453843075888, 0.16554538430758944, 0.6182451347074112],
      linf = [2.4849541392059202, 1.2816207838889486, 1.2816207838888944, 6.4743208586529875],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl"),
      l2   = [8.29516192e-02, 9.18870398e-02, 9.18870398e-02, 4.37449170e-01],
      linf = [1.23402442e+00, 1.63254552e+00, 1.63254552e+00, 4.21977821e+00],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_positivity.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"),
      l2   = [0.48401054859424775, 0.16613328471131006, 0.16613328471131006, 0.6181793382171535],
      linf = [2.4903327485761544, 1.2898210694161085, 1.2898210694161072, 6.4723873993158385],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_blob_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_mortar.jl"),
      l2   = [0.22206706482833133, 0.6296849068200803, 0.24351238064525166, 2.9118777552679345],
      linf = [9.584083314003589, 27.72635365994101, 9.343044146443756, 101.53747277352629],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blob_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
      l2   = [0.20150296718390992, 1.1816675092370863, 0.10129652123574924, 5.230357225596295],
      linf = [14.174958287930844, 71.16876942972455, 7.2799739019809415, 291.67608863630244],
      tspan = (0.0, 0.12))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      l2   = [0.055689867638859906, 0.032985369421826734, 0.052243738479375024, 0.08009259073134713],
      linf = [0.24051443859674326, 0.16611572070980668, 0.123559477474734, 0.2694484247345663],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr.jl"),
      l2   = [0.05569140599245201, 0.033108757924402134, 0.05223526059785845, 0.08007239985228143],
      linf = [0.2540392117514094, 0.17432468822204936, 0.12323401271048683, 0.26897998287356195],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl"),
      l2   = [0.0009932162527962933, 0.004946122604416242, 0.004863174253792848, 0.005019630958200973],
      linf = [0.008603757244273202, 0.02090493316272357, 0.01712529102069524, 0.028966987799588395],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_colliding_flow.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"),
      l2   = [7.23705116e-03,   4.48861301e-02,   8.47273157e-07,   6.62667920e-01],
      linf = [1.93754293e-01,   5.52253475e-01,   4.98055004e-05,   1.50879674e+01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_colliding_flow_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow_amr.jl"),
      l2   = [6.76960244e-03,   3.21735698e-02,   2.63316434e-07,   6.78416591e-01],
      linf = [2.50119529e-01,   4.07493507e-01,   9.96993324e-05,   2.23204482e+01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_astro_jet_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_astro_jet_amr.jl"),
      l2   = [0.011338493575859507, 10.097552800226566, 0.0038670992882799088, 4031.828021246965],
      linf = [3.316670548834158, 2992.6286682433883, 7.919515633718346, 1.1914747265545786e6],
      tspan = (0.0, 1.0e-7))
  end

  @trixi_testset "elixir_euler_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [3.6343141788303172e-6, 0.0032111379945554378, 0.0032111482803927763, 0.004545715899533244],
      linf = [7.901851921288117e-5, 0.030561510550305093, 0.03050266851613559, 0.042876690674344076])
  end

  @trixi_testset "elixir_euler_vortex_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar.jl"),
      l2   = [2.1203067597905396e-6, 2.792922844344066e-5, 3.7593422075434476e-5, 8.813644239472934e-5],
      linf = [5.932046017209647e-5, 0.0007491265252463908, 0.0008165690091972433, 0.0022122638482677814])
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1201025902895562e-6, 2.8056524989293713e-5, 3.759500393896144e-5, 8.84137216637438e-5],
      linf = [5.934103282212444e-5, 0.0007552316670187409, 0.0008152449093751235, 0.0022069874183294758])
  end

  @trixi_testset "elixir_euler_vortex_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_shockcapturing.jl"),
      l2   = [3.7998676433961623e-6, 5.563687203373343e-5, 5.561628017766063e-5, 0.00015707075214260475],
      linf = [8.522682141076654e-5, 0.0009610342345369727, 0.0009656468947373265, 0.0030861481615822584])
  end

  @trixi_testset "elixir_euler_vortex_mortar_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_shockcapturing.jl"),
      l2   = [2.1203212494737245e-6, 2.8056721454951553e-5, 3.759510818832351e-5, 8.842116438932004e-5],
      linf = [5.9340649220973596e-5, 0.0007552376287922602, 0.0008152432005811283, 0.0022069702233196153])
  end

  @trixi_testset "elixir_euler_vortex_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_amr.jl"),
      l2   = [2.1325938318497577e-6, 0.003281545066125044, 0.0032806321457482615, 0.00464605975938406],
      linf = [4.667437348149228e-5, 0.03175420871507906, 0.0318039789241531, 0.04561735256198318])
  end

  @trixi_testset "elixir_euler_ec.jl with boundary_condition_slip_wall" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.06176845600430613, 0.05020624043492084, 0.05021389111189423, 0.22592682624517807],
      linf = [0.29347582879609024, 0.3108124923286465, 0.3107380389949771, 1.054035804988522],
      periodicity = false, boundary_conditions = boundary_condition_slip_wall)
  end
end

# Coverage test for all initial conditions
@testset "Compressible Euler: Tests for initial conditions" begin
  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_density_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.003489659044164644, 0.0034896590441646494, 0.0034896590441646502, 0.003489659044164646],
      linf = [0.04761180654650543, 0.04761180654650565, 0.047611806546505875, 0.04761180654650454],
      maxiters = 1,
      initial_condition = initial_condition_density_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.00021747693728234874, 0.0022010142997830533, 0.0022010142997830485, 0.010855273768135729],
      linf = [0.005451116856088789, 0.03126448432601536, 0.03126448432601536, 0.14844305553724624],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_pressure_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_density_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.003473649182284682, 0.005490887132955628, 0.005490887132955635, 0.015625074774949926],
      linf = [0.046582178207169145, 0.07332265196082899, 0.07332265196082921, 0.2107979471941368],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_density_pressure_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [7.89034964747135e-17, 8.095575651413758e-17, 1.0847287658433571e-16, 1.2897732640029767e-15],
      linf = [2.220446049250313e-16, 3.191891195797325e-16, 4.163336342344337e-16, 3.552713678800501e-15],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl one step" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0021037031798961914, 0.010667428589443025, 0.01066742858944302, 0.10738893384136498],
      linf = [0.11854059147646778, 0.7407961272348982, 0.7407961272348981, 3.92623931433345],
      maxiters=1)
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl one step with initial_condition_medium_sedov_blast_wave" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0021025532272874827, 0.010661548568022292, 0.010661548568022284, 0.10734939168392313],
      linf = [0.11848345578926645, 0.7404217490990809, 0.7404217490990809, 3.9247328712525973],
      maxiters=1, initial_condition=initial_condition_medium_sedov_blast_wave)
  end
end

end # module

ext: jl


text:
# ---
# title: 206. Reverse Linked List
# id: problem206
# author: zhwang
# date: 2022-02-17
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-linked-list/description/>
# hidden: true
# ---
#
# Reverse a singly linked list.
#
# **Example:**
#
#
#
#     Input: 1->2->3->4->5->NULL
#     Output: 5->4->3->2->1->NULL
#
#
# **Follow up:**
#
# A linked list can be reversed either iteratively or recursively. Could you
# implement both?
#
#
## @lc code=start
using LeetCode

function reverse_list(head::Union{ListNode, Nothing})::Union{ListNode, Nothing}
    new_head = pre_node = nothing
    while !isnothing(head)
        new_head = ListNode(head.val)
        pre_node, new_head.next = new_head, pre_node
        head = head.next
    end
    new_head
end
## @lc code=end
## @lc test=start
@testset "206.reverse-linked-list.jl" begin
    @test reverse_list(ListNode{Int}([1,2,3,4,5])) == ListNode{Int}([5,4,3,2,1])
    @test reverse_list(ListNode{Int}([1,2])) == ListNode{Int}([2,1])
    @test reverse_list(nothing) === nothing
end
## @lc test=end

## @lc info=start
# link: [solution to 206](https://leetcode-cn.com/problems/reverse-linked-list/solution/206-fan-zhuan-lian-biao-by-rex-rs-zon0/)
## @lc info=end


text:
using Images, MXNet #,ImageView

# start by loading a single image
IMAGES_PATH = "data/train"

single_img = load(joinpath(IMAGES_PATH, "cat.1.jpg"))
size(single_img)

preview_img = zeros(100, 0);

for i = 1:10
    seq_x_img = vcat(
        imresize(load(joinpath($IMAGES_PATH, "cat.$i.jpg")), (50, 50)),
        imresize(load(joinpath($IMAGES_PATH, "dog.$i.jpg")), (50, 50))
    )
    preview_img = hcat(preview_img, seq_x_img)
end

imshow(preview_img)

### PREPARING THE DATASET

files = readdir(IMAGES_PATH);
data_x = zeros((32, 32, 3, size(files, 1)));
data_y = zeros(size(files, 1));

for idx = 1:size(files, 1)
    file_name = joinpath(IMAGES_PATH, files[idx])

    if endswith(file_name, ".jpg")
        img = imresize(load(file_name), (32, 32))
        try
            data_x[:, :, :, idx] = permuteddimsview(Float16.(channelview(img)), (2, 3, 1))
            data_y[idx] = 1 * contains(files[idx], "dog")
        catch
            data_x[:, :, :, idx] = permuteddimsview(Float16.(channelview(RGB.(img))), (2, 3, 1))
            data_y[idx] = 1 * contains(files[idx], "dog")
        end
    end
end

total_count = size(data_y, 1);
indices = shuffle(1:total_count);

total_train_count = Int(total_count * 0.8);
total_validation_count = total_count - total_train_count

train_data_array  = mx.zeros((size(data_x, 1, 2, 3)..., total_train_count...));
train_label_array = mx.zeros(total_train_count);

validation_data_array  = mx.zeros((size(data_x, 1, 2, 3)..., total_validation_count...));
validation_label_array = mx.zeros(total_validation_count);

for idx = 1:total_train_count
    train_data_array[idx:idx] = reshape(data_x[:, :, :, indices[idx]], (size(data_x, 1, 2, 3 )..., 1...))
    train_label_array[idx:idx] = data_y[indices[idx]]
end

for idx = 1:total_validation_count
    validation_data_array[idx:idx] = reshape(data_x[:, :, :, indices[total_train_count + idx]], (size(data_x, 1, 2,3 )..., 1...))
    validation_label_array[idx:idx] = data_y[indices[total_train_count + idx]]
end

train_data_provider = mx.ArrayDataProvider(:data => train_data_array, :label => train_label_array, batch_size = 100, shuffle = true);
validation_data_provider = mx.ArrayDataProvider(:data => validation_data_array, :label => validation_label_array, batch_size = 100, shuffle = true);

arch = @mx.chain mx.Variable(:data) =>
        mx.Convolution(kernel=(3, 3), num_filter=32) =>
        mx.Activation(act_type=:relu) =>
        mx.Dropout(p = 0.25) =>
        mx.Pooling( kernel=(2, 2), pool_type=:max) =>
        mx.Flatten() =>
        mx.FullyConnected(num_hidden=256) =>
        mx.Activation(act_type=:relu) =>
        mx.FullyConnected(num_hidden=10) =>
        mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 50, initializer = mx.XavierInitializer());

### Improved neural network

arch = @mx.chain mx.Variable(:data) =>
        mx.Convolution(kernel=(3, 3), num_filter=32) =>
        mx.Activation(act_type=:relu) =>
        mx.Convolution(kernel=(3, 3), num_filter=32) =>
        mx.Activation(act_type=:relu) =>
        mx.Dropout(p = 0.25) =>
        mx.Pooling( kernel=(3, 3), pool_type=:max) =>
        mx.Flatten() =>
        mx.FullyConnected(num_hidden=256) =>
        mx.Activation(act_type=:relu) =>
        mx.FullyConnected(num_hidden=10) =>
        mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 50, initializer = mx.XavierInitializer());


text:
module TwoDimensionalTests

using TwoDimensional
using TwoDimensional: WeightedPoint, half, getaxisbounds
import TwoDimensional: distance
using Test, LinearAlgebra
import Base.MathConstants: φ

distance(a::Real, b::Real) = distance(promote(a, b)...)
distance(a::T, b::T) where {T<:Unsigned} = ifelse(a < b, b - a, a - b)
distance(a::T, b::T) where {T<:Real} = abs(a - b)

distance(A::NTuple{2,Real}, B::NTuple{2,Real}) =
    hypot(A[1] - B[1], A[2] - B[2])

distance(A::AffineTransform, B::AffineTransform) =
    max(abs(A.xx - B.xx), abs(A.xy - B.xy), abs(A.x - B.x),
        abs(A.yx - B.yx), abs(A.yy - B.yy), abs(A.y - B.y))

function randomize!(A::Matrix{Bool}, bits::Int)
    @inbounds @simd for i in eachindex(A)
        A[i] = (bits&1) == 1
        bits >>= 1
    end
    return A
end

# Naive implementation of the bounding box algorithm.
naiveboundingbox(A::AbstractMatrix{Bool}) = naiveboundingbox(identity, A)
function naiveboundingbox(f::Function, A::AbstractMatrix)
    I, J = axes(A)
    imin = jmin = typemax(Int)
    imax = jmax = typemin(Int)
    @inbounds for j in J, i in I
        if f(A[i,j])
            imin = min(imin, i)
            imax = max(imax, i)
            jmin = min(jmin, j)
            jmax = max(jmax, j)
        end
    end
    return BoundingBox(imin, imax, jmin, jmax)
end

@testset "Points and bounding-boxes" begin
    @testset "Miscellaneous" begin
        # Iterators.
        pnt = Point(2,3)
        @test Point(Tuple(pnt)...) === pnt
        @test Point(Tuple(pnt)) === pnt
        x,y = pnt
        @test (x,y) === (pnt.x, pnt.y)
        @test Tuple(pnt) === (x,y)

        wpnt = WeightedPoint(1.2,sqrt(2),3)
        @test WeightedPoint(Tuple(wpnt)...) === wpnt
        @test WeightedPoint(Tuple(wpnt)) === wpnt
        w,x,y = wpnt
        @test (w,x,y) === (wpnt.w, wpnt.x, wpnt.y)
        @test Tuple(wpnt) === (w,x,y)

        box = BoundingBox(1.2,sqrt(2),-3,11)
        @test BoundingBox(Tuple(box)...) === box
        @test BoundingBox(Tuple(box)) === box
        xmin,xmax,ymin,ymax = box
        @test (xmin,xmax,ymin,ymax) === (box.xmin,box.xmax,box.ymin,box.ymax)
        @test Tuple(box) === (xmin,xmax,ymin,ymax)
    end
    @testset "Simple points" begin
        types = (Int8, Int32, Int64, Float32, Float64)
        for T1 in types, T2 in types
            @test promote_type(Point{T1}, Point{T2}) ===
                Point{promote_type(T1,T2)}
        end
        P1 = Point(1, 1.3)
        @test eltype(P1) == Float64
        @test Point(P1) === P1
        @test Point{eltype(P1)}(P1) === P1
        @test eltype(Point(1,2)) == Int
        @test Point(1,2) === Point(y=2, x=1)
        @test Point(CartesianIndex(7,8)) === Point(7,8)
        @test CartesianIndex(Point(2,3)) === CartesianIndex(2,3)
        @test Point(TwoDimensional.WeightedPoint(2,3,7)) === Point(3.0,7.0)
        #@test convert(CartesianIndex, Point(2,3)) === CartesianIndex(2,3)
        #@test convert(CartesianIndex{2}, Point(2,3)) === CartesianIndex(2,3)
        #@test convert(Tuple, Point(2,3)) === (2,3)
        #@test convert(Tuple{Float64,Int}, Point(2,3)) === (2.0,3)
        @test Point{Float32}(P1) === Float32.(P1)
        @test promote(P1) === (P1,)
        @test Point{Int}(x=Int16(8),y=Int16(9)) === Point(8,9)
        @test Point{Float64}(CartesianIndex(8,9)) === Point(8.0,9.0)
        @test Point{Float64}((8,9)) === Point(8.0,9.0)
        P2 = Point(3,11)
        P3 = Point{Float32}(-2.1,11.8)
        @test promote(P1, P2) === (P1,Float64.(P2))
        @test promote(P1, P2, P3) === (P1,Float64.(P2),Float64.(P3))
        @test convert(Point, P1) === P1
        @test convert(Point{Float32}, P1) === Float32.(P1)
        for T in (Float32, Int16)
            @test zero(Point{T}) === Point(zero(T),zero(T))
            @test one(Point{T}) === Point(one(T),one(T))
            @test oneunit(Point{T}) === Point(oneunit(T),oneunit(T))
        end
        @test ntuple(i -> Point(3,5)[i], 2) == (3,5)
        # round
        @test round(Point(1.2,-0.7)) === Point(1.0,-1.0)
        @test round(Point(1,-7)) === Point(1,-7)
        @test round(Float32, Point(1,-7)) === Point{Float32}(1,-7)
        @test round(Float32, Point(1.0,-7.6)) === Point{Float32}(1,-8)
        @test round(Int32, Point(1,-7)) === Point{Int32}(1,-7)
        @test round(Int, Point(1.2,-0.7)) === Point(1,-1)
        @test round(Int, Point(1,-4)) === Point(1,-4)
        @test round(Point{Float64}, P1) === Point(map(round, Tuple(P1)))
        @test round(Point{Int}, Point(1.6,-0.7)) === Point(2,-1)
        @test round(Int, Point(1.2,-0.7)) === Point(1,-1)
        @test round(Point{Int}, Point(1.6,-0.7)) === Point(2,-1)
        # floor
        @test floor(Point(-1,3)) === Point(-1,3)
        @test floor(Int, Point(-1,3)) === Point(-1,3)
        @test floor(Point{Int}, Point(-1,3)) === Point(-1,3)
        @test floor(Int16, Point(-1,3)) === Point{Int16}(-1,3)
        @test floor(Point{Int16}, Point(-1,3)) === Point{Int16}(-1,3)
        @test floor(Point(-1.7,3.2)) === Point(-2.0,3.0)
        @test floor(Int, Point(-1.7,3.2)) === Point(-2,3)
        @test floor(Float32, Point(-1,3)) === Point{Float32}(-1,3)
        @test floor(Point{Float32}, Point(-1.7,3.2)) === Point{Float32}(-2,3)
        # ceil
        @test ceil(Point(-1,3)) === Point(-1,3)
        @test ceil(Int, Point(-1,3)) === Point(-1,3)
        @test ceil(Point{Int}, Point(-1,3)) === Point(-1,3)
        @test ceil(Int16, Point(-1,3)) === Point{Int16}(-1,3)
        @test ceil(Point{Int16}, Point(-1,3)) === Point{Int16}(-1,3)
        @test ceil(Point(-1.7,3.2)) === Point(-1.0,4.0)
        @test ceil(Int, Point(-1.7,3.2)) === Point(-1,4)
        @test ceil(Float32, Point(-1,3)) === Point{Float32}(-1,3)
        @test ceil(Point{Float32}, Point(-1.7,3.2)) === Point{Float32}(-1,4)
        # other methods
        @test clamp(Point(-1.1, 6.3), BoundingBox(1,4,1,5)) === Point(1.0,5.0)
        @test hypot(P1) == hypot(P1.x, P1.y)
        @test atan(P1) == atan(P1.y, P1.x)
        @test distance(P1, Point(0,0)) == hypot(P1)
        @test distance(Point(0x02,0x05), Point(0x00,0x00)) == hypot(0x02,0x05)
    end
    @testset "Weighted points" begin
        types = (Float32, Float64)
        for T1 in types, T2 in types
            @test promote_type(WeightedPoint{T1}, WeightedPoint{T2}) ===
                WeightedPoint{promote_type(T1,T2)}
        end
        P2 = WeightedPoint(x=0.1, y=-6, w=0.1)
        @test eltype(P2) == Float64
        @test WeightedPoint(P2) === P2
        @test WeightedPoint{eltype(P2)}(P2) === P2
        @test WeightedPoint{Float64}(P2) === P2
        @test WeightedPoint(P2.w, P2.x, P2.y) === P2
        @test WeightedPoint{Float32}(P2) === Float32.(P2)
        @test WeightedPoint(Point(3.1,4.2)) ===
            WeightedPoint(w=1, x=3.1, y=4.2)
    end
    @testset "Bounding boxes" begin
        types = (Int8, Int32, Int64, Float32, Float64)
        for T1 in types, T2 in types
            @test promote_type(BoundingBox{T1}, BoundingBox{T2}) ===
                BoundingBox{promote_type(T1,T2)}
        end
        B = BoundingBox(2,3,4,5)
        δ = 0.1
        @test eltype(B) === Int
        @test BoundingBox(B) === B
        @test BoundingBox{eltype(B)}(B) === B
        @test BoundingBox{Float32}(B) ===
            BoundingBox{Float32}(B.xmin, B.xmax, B.ymin, B.ymax)
        @test BoundingBox(2,3,4,5.0) === BoundingBox(2.0,3.0,4.0,5.0)
        @test convert(BoundingBox, B) === B
        @test convert(BoundingBox{eltype(B)}, B) === B
        @test convert(BoundingBox{Int16}, B) === BoundingBox{Int16}(B)
        # Construct a bounding-box from 4-tuple.
        @test BoundingBox((2,3,4,5)) === BoundingBox(2,3,4,5)
        @test BoundingBox{Float64}((2,3,4,5)) === BoundingBox(2.0,3.0,4.0,5.0)
        @test convert(BoundingBox, (2,3,4,5)) === BoundingBox(2,3,4,5)
        @test convert(BoundingBox{Float64}, (2,3,4,5)) ===
            BoundingBox(2.0,3.0,4.0,5.0)
        # Construct a bounding-box from 2 points.
        @test BoundingBox(Point(2,3),Point(4,5)) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}(Point(2,3),Point(4,5)) ===
            BoundingBox{Int16}(2,4,3,5)
        @test BoundingBox((Point(2,3),Point(4,5))) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}((Point(2,3),Point(4,5))) ===
            BoundingBox{Int16}(2,4,3,5)
        @test convert(BoundingBox,(Point(2,3),Point(4,5))) ===
            BoundingBox(2,4,3,5)
        @test convert(BoundingBox{Int16}, (Point(2,3),Point(4,5))) ===
            BoundingBox{Int16}(2,4,3,5)
        # Construct a bounding-box from 2 Cartesian indices.
        @test BoundingBox(CartesianIndex(2,3),CartesianIndex(4,5)) ===
            BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}(CartesianIndex(2,3),CartesianIndex(4,5)) ===
            BoundingBox{Int16}(2,4,3,5)
        @test BoundingBox((CartesianIndex(2,3),CartesianIndex(4,5))) ===
            BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}((CartesianIndex(2,3),CartesianIndex(4,5))) ===
            BoundingBox{Int16}(2,4,3,5)
        # Construct a bounding-box from 2 2-tuple.
        @test BoundingBox((2,3),(4,5)) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}((2,3),(4,5)) === BoundingBox{Int16}(2,4,3,5)
        # Construct a bounding-box by keywords.
        @test BoundingBox(xmin=2,ymin=3,xmax=4,ymax=5) ===
            BoundingBox(2,4,3,5)
        @test BoundingBox{Float64}(xmin=2,ymin=3,xmax=4,ymax=5) ===
            BoundingBox(2.0,4.0,3.0,5.0)
        # Construct a bounding-box from CartesianIndices instance.
        R = CartesianIndices((2:4,3:5))
        @test BoundingBox(R) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}(R) === BoundingBox{Int16}(2,4,3,5)
        @test convert(BoundingBox, R) === BoundingBox(2,4,3,5)
        @test convert(BoundingBox{Int16}, R) === BoundingBox{Int16}(2,4,3,5)
        # Construct a bounding-box from 2 ranges.
        R = (2:4, 3:5)
        @test BoundingBox(R...) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}(R...) === BoundingBox{Int16}(2,4,3,5)
        @test BoundingBox(R) === BoundingBox(2,4,3,5)
        @test BoundingBox{Int16}(R) === BoundingBox{Int16}(2,4,3,5)
        @test convert(BoundingBox,R) === BoundingBox(2,4,3,5)
        @test convert(BoundingBox{Int16},R) === BoundingBox{Int16}(2,4,3,5)

        @test first(B) === Point(B.xmin,B.ymin)
        @test last(B) === Point(B.xmax,B.ymax)
        @test isempty(typemax(BoundingBox{Int})) == false
        @test isempty(typemin(BoundingBox{Float64})) == true
        @test BoundingBox{Float32}(nothing) === typemin(BoundingBox{Float32})
        @test size(BoundingBox{Int16}(nothing)) === (0,0)
        @test ntuple(i -> BoundingBox(3:7,5:8)[i], 4) == (3,7,5,8)
        Rx, Ry = 3:7, -6:-3
        B = BoundingBox(Rx,Ry)
        @test size(B) === (length(Rx),length(Ry))
        @test size(BoundingBox{Int16}(B)) === size(B)
        @test axes(B) === (Rx,Ry)
        @test axes(BoundingBox{Int16}(B)) === (Rx,Ry)
        for k in (1,2)
            @test size(B, k) === size(B)[k]
            @test axes(B, k) === axes(B)[k]
            @test size(BoundingBox{Int16}(B), k) === size(B)[k]
            @test axes(BoundingBox{Int16}(B), k) === axes(B)[k]
        end
        @test size(B,40) === 1
        @test_throws ErrorException size(B, 0)
        B1 = BoundingBox(2,3,4,5)
        B2 = BoundingBox(1.1,3.2,4.5,5.8)
        B3 = BoundingBox{Float32}(1.1,3.2,4.5,5.8)
        @test promote(B1, B2) === (Float64.(B1),B2)
        @test promote(B1, B2, B3) === (Float64.(B1),B2,Float64.(B3))
        @test_deprecated BoundingBox(ones(5,7)) === BoundingBox(1,5, 1,7)
        @test BoundingBox(-2:6,8:11) === BoundingBox(-2,6, 8,11)
        @test BoundingBox((2:4,-1:7)) === BoundingBox(2,4, -1,7)
        @test CartesianIndices(BoundingBox(2:4,-1:7)) ===
            CartesianIndices((2:4,-1:7))
        @test getaxisbounds(9:13) === (9,13)
        @test getaxisbounds(13:9) === (13,12)
        @test getaxisbounds(13:-1:9) === (9,13)
        @test_throws ArgumentError getaxisbounds(13:-2:9)
        let A = ones(7,8)
            A[1:2,:] .= 0
            A[4,:] .= 0
            A[end,:] .= 0
            @test BoundingBox(x -> x > 0, A) === BoundingBox(3,6, 1,8)
            @test BoundingBox(A .> 0) === BoundingBox(3,6, 1,8)
            A[:,1] .= 0
            A[:3:4] .= 0
            A[:,end-1:end] .= 0
            A[2,2] = 1
            @test BoundingBox(x -> x > 0, A) === BoundingBox(2,6, 2,6)
            @test BoundingBox(A .> 0) === BoundingBox(2,6, 2,6)
        end
        @test BoundingBox{Float32}(nothing) ===
            BoundingBox{Float32}(Inf,-Inf,Inf,-Inf)
        @test typemin(BoundingBox{Float64}) === BoundingBox(Inf,-Inf,Inf,-Inf)
        @test typemax(BoundingBox{Float64}) === BoundingBox(-Inf,Inf,-Inf,Inf)
        # round
        @test round(BoundingBox(1.1,2.1,-3.6,7.7)) ===
            BoundingBox(1.0,2.0,-4.0,8.0)
        @test round(BoundingBox{Int16}, BoundingBox(1.1,2.1,-3.6,7.7)) ===
            BoundingBox{Int16}(1,2,-4,8)
        @test round(Int, BoundingBox(1,2,-3,7)) === BoundingBox(1,2,-3,7)
        @test round(Int16, BoundingBox(1.1,2.1,-3.6,7.7)) ===
            BoundingBox{Int16}(1,2,-4,8)
        @test round(Int16, BoundingBox(1,2,-4,8)) ===
            BoundingBox{Int16}(1,2,-4,8)
        @test round(Float32, BoundingBox{Int16}(1,2,-4,8)) ===
            BoundingBox{Float32}(1,2,-4,8)
        @test round(Float32, BoundingBox(1.1,2.7,-4.6,8.3)) ===
            BoundingBox{Float32}(1,3,-5,8)
        # exterior
        @test exterior(B) === B
        @test exterior(Int, B) === B
        @test exterior(Int16, B) === BoundingBox{Int16}(B)
        @test exterior(Float32, B) === BoundingBox{Float32}(B)
        @test exterior(BoundingBox{Int}, B) === B
        @test exterior(B + δ) === B + 1.0
        @test exterior(Int, B + δ) === B + 1
        @test exterior(BoundingBox{Int}, B) === BoundingBox{Int}(exterior(B))
        @test exterior(Float32, B + δ) ===
            BoundingBox{Float32}(exterior(B + δ))
        # interior
        @test interior(B) === B
        @test interior(Int, B) === B
        @test interior(Int16, B) === BoundingBox{Int16}(B)
        @test interior(Float32, B) === BoundingBox{Float32}(B)
        @test interior(BoundingBox{Int}, B) === B
        @test interior(B + δ) === Float64.(B)
        @test interior(Int, B - δ) === B - 1
        @test interior(BoundingBox{Int}, B) === BoundingBox{Int}(interior(B))
        @test interior(Float32, B + δ) === BoundingBox{Float32}(interior(B + δ))
        # other methods
        @test area(BoundingBox(2,4,5,8)) == 6
        @test area(BoundingBox(2.0,4.0,5.0,8.0)) == 6.0

        @test (Point(2,4) ∈ BoundingBox(1,2,3,4)) == true
        @test (Point(3,4) ∈ BoundingBox(1,2,3,4)) == false
        @test (Point(2,5) ∈ BoundingBox(1,2,3,4)) == false
        @test (Point(3,5) ∈ BoundingBox(1,2,3,4)) == false
        @test (WeightedPoint(1,2,4) ∈ BoundingBox(1,2,3,4)) == true
        @test (WeightedPoint(1,3,4) ∈ BoundingBox(1,2,3,4)) == false
        @test (WeightedPoint(1,2,5) ∈ BoundingBox(1,2,3,4)) == false
        @test (WeightedPoint(1,3,5) ∈ BoundingBox(1,2,3,4)) == false

        @test (BoundingBox(1,-2,3,4) ⊆ BoundingBox{Float32}(nothing)) == true
        @test (BoundingBox(1,-2,3,4) ⊆ BoundingBox(1,2,3,4)) == true
        @test (BoundingBox(1,2,3,4) ⊆ BoundingBox(1,2,3,4)) == true
        @test (BoundingBox(1,2,3,4) ⊆ BoundingBox(1,4,0,3)) == false
        @test (BoundingBox(0,2,3,4) ⊆ BoundingBox(1,2,3,5)) == false
        @test_deprecated (BoundingBox(0,2,3,4) ∈ BoundingBox(1,2,3,5)) == false

        C = Point(x = 0.5*(B.xmin + B.xmax), y = 0.5*(B.ymin + B.ymax))
        @test center(B) === C
        @test center(BoundingBox{Float64}(B)) === C
        B1 = BoundingBox(-2:6, -7:-1)
        B2 = BoundingBox(1:8, -9:3)
        @test B1 ∪ B2 === BoundingBox(-2:8, -9:3)
        @test B1 ∩ B2 === BoundingBox(1:6, -7:-1)
        A = rand(7,8)
        X, Y = 2:4, 1:3
        B = BoundingBox(X, Y)
        @test A[X,Y] == A[B]
        @test view(A,X,Y) === view(A, B)
    end
    @testset "Arithmetic" begin
        @test half(Float64) == 0.5
        @test half(Float32) == one(Float32)/2
        @test 2*Point(3,4) === Point(6,8)
        @test Point(3,4)*3 === Point(9,12)
        @test 3\Point(3,4) === Point(3/3,4/3)
        @test Point(3,4)/3 === Point(3/3,4/3)
        @test +Point(2,-9) === Point(2,-9)
        @test -Point(2,-9) === Point(-2,9)
        @test +BoundingBox(2,3,4,5) === BoundingBox(2,3,4,5)
        @test -BoundingBox(2,3,4,5) === BoundingBox(-3,-2,-5,-4)
        @test Point(2,-9) + Point(3,7) === Point(5,-2)
        @test Point(2,-9) - Point(3,7) === Point(-1,-16)
        α = 3
        δ = 0.1
        t = Point(-0.2,0.7)
        B = BoundingBox(2,3,4,5)
        @test α*B === BoundingBox(xmin = α*B.xmin, xmax = α*B.xmax,
                                  ymin = α*B.ymin, ymax = α*B.ymax)
        @test B*α === BoundingBox(xmin = α*B.xmin, xmax = α*B.xmax,
                                  ymin = α*B.ymin, ymax = α*B.ymax)
        @test α\B === BoundingBox(xmin = (1/α)*B.xmin, xmax = (1/α)*B.xmax,
                                  ymin = (1/α)*B.ymin, ymax = (1/α)*B.ymax)
        @test B/α === BoundingBox(xmin = (1/α)*B.xmin, xmax = (1/α)*B.xmax,
                                  ymin = (1/α)*B.ymin, ymax = (1/α)*B.ymax)
        @test B + δ === BoundingBox(xmin = B.xmin - δ, xmax = B.xmax + δ,
                                    ymin = B.ymin - δ, ymax = B.ymax + δ)
        @test B - δ === BoundingBox(xmin = B.xmin + δ, xmax = B.xmax - δ,
                                    ymin = B.ymin + δ, ymax = B.ymax - δ)
        @test B + t === BoundingBox(xmin = B.xmin + t.x, xmax = B.xmax + t.x,
                                    ymin = B.ymin + t.y, ymax = B.ymax + t.y)
        @test B - t === BoundingBox(xmin = B.xmin - t.x, xmax = B.xmax - t.x,
                                    ymin = B.ymin - t.y, ymax = B.ymax - t.y)
    end
end

@testset "Bounding-box algorithm" begin
    if true
        # Exhaustively test all possibilities.
        A = Array{Bool,2}(undef, (5,4))
        n = (1 << length(A)) # number of possibilities
        for bits in 0:n-1
            randomize!(A, bits)
            @test BoundingBox(A) === naiveboundingbox(A)
        end
    else
        # Semi-exhaustive testing of the bounding box algorithm.
        A = zeros(Bool, (7,8))
        fill!(A,false)[1,3] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,true)[2:end-1,2:end-1] .= false
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[1,1] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[1,end] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[end,1] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[end,end] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[3,4] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        A[2,7] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        A[5,end] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[end,4] = true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[1,:] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[:,1] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[end,:] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[:,end] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[1,2:end-1] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[2:end-1,1] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[end,2:end-1] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
        fill!(A,false)[2:end-1,end] .= true
        @test BoundingBox(A) === naiveboundingbox(A)
    end
end

@testset "AffineTransforms" begin
    tol = 1e-14
    I = AffineTransform()
    A = AffineTransform(1, 0, -3, 0.1, 1, +2)
    B = AffineTransform(-0.4,  0.1, -4.2, -0.3,  0.7,  1.1)
    C = AffineTransform( 2.3, -0.9, -6.1,  0.7, -3.1, -5.2)
    vectors = ((0.2,1.3), (-1,π), (-sqrt(2),3//4))
    scales = (2, 0.1, φ)
    angles = (-2π/11, π/7, 0.1)
    types = (BigFloat, Float64, Float32, Float16)

    @testset "construction/conversion" begin
        @test_deprecated BigFloat(A) === BigFloat.(A)
        @test_deprecated Float32(A) === Float32.(A)
        for T1 in types, T2 in types
            @test promote_type(AffineTransform{T1}, AffineTransform{T2}) ===
                AffineTransform{promote_type(T1,T2)}
        end
        @test AffineTransform(A) === A
        @test AffineTransform{eltype(A)}(A) === A
        for G in (I, A, B)
            @test eltype(G) == Float64
            for T in types
                @test typeof(AffineTransform{T}(G)) == AffineTransform{T}
                @test typeof(convert(AffineTransform{T}, G)) ==
                    AffineTransform{T}
                @test eltype(AffineTransform{T}(G)) == T
                @test eltype(T.(G)) == T
            end
        end
    end

    @testset "identity" begin
        @test det(I) == 1
        @test distance(inv(I), I) ≤ 0
        for v in vectors
            @test distance(I(v), eltype(I).(v)) ≤ 0
        end
    end

    @testset "apply" begin
        for G in (I, A, B),
            v in vectors
            @test distance(G(v...), G(v)) ≤ 0
            @test distance(G*v, G(v)) ≤ 0
            @test distance(G(Point(v)), Point(G(v))) ≤ 0
            @test distance(G*Point(v), Point(G(v))) ≤ 0
        end
    end

    @testset "composition" begin
        for G in (I, B, A)
            @test distance(G, compose(G)) ≤ 0
            for H in (A, B)
                @test distance(G*H, compose(G,H)) ≤ 0
                @test distance(G⋅H, compose(G,H)) ≤ 0
                @test distance(G∘H, compose(G,H)) ≤ 0
                for v in vectors
                    @test distance((G*H)(v), G(H(v))) ≤ tol
                end
            end
        end
        for T1 in types, T2 in types
            T = promote_type(T1, T2)
            @test eltype(T1.(A)*T2.(B)) == T
        end
        for v in vectors
            @test distance(compose(A,B,C)(v), A(B(C(v)))) ≤ tol
            @test distance((A*B*C)(v), A(B(C(v)))) ≤ tol
            @test distance(compose(A,C,B,C)(v), A(C(B(C(v))))) ≤ tol
            @test distance((A*C*B*C)(v), A(C(B(C(v))))) ≤ tol
        end
    end

    @testset "jacobian" begin
        for M in (I, B, A)
            @test jacobian(M) == abs(det(M))
        end
    end

    @testset "inverse" begin
        for M in (B, A)
            if det(M) == 0
                continue
            end
            @test distance(det(inv(M)), 1/det(M)) ≤ tol
            @test distance(M/M, M*inv(M)) ≤ tol
            @test distance(M\M, inv(M)*M) ≤ tol
            @test distance(M\M, I) ≤ tol
            @test distance(M/M, I) ≤ tol
            for v in vectors
                @test distance(M(inv(M)(v)), v) ≤ tol
                @test distance(inv(M)(M(v)), v) ≤ tol
                @test distance((M\M)(v), v) ≤ tol
                @test distance((M/M)(v), v) ≤ tol
            end
        end
        for T1 in types, T2 in types
            T = promote_type(T1, T2)
            @test eltype(T1.(A)/T2.(B)) == T
            @test eltype(T1.(A)\T2.(B)) == T
        end
    end

    @testset "scale" begin
        for M in (A, B, C),
            α in scales,
            v in vectors
            @test distance((α*M)(v), α.*M(v)) ≤ tol
            @test distance((M*α)(v), M(α.*v)) ≤ tol
        end
        for G in (A, B, C),
            α in scales,
            T in types
            @test eltype(T(α)*G) == eltype(G)
            @test eltype(G*T(α)) == eltype(G)
            H = T.(G)
            @test eltype(α*H) == eltype(H)
            @test eltype(H*α) == eltype(H)
        end
    end

    @testset "translation" begin
        for M in (B, A),
            t in vectors,
            v in vectors
            @test distance(translate(t, M)(v), t .+ M(v)) ≤ tol
            @test distance(translate(t, M)(v), (t + M)(v)) ≤ tol
            @test distance(translate(M, t)(v), M(v .+ t)) ≤ tol
            @test distance(translate(M, t)(v), (M + t)(v)) ≤ tol
        end
        for v in vectors
            @test distance(A - v, A + (-v[1], -v[2])) ≤ tol
            @test distance(A + Point(v), translate(A, v...)) ≤ tol
            @test distance(Point(v) + A, translate(v..., A)) ≤ tol
            @test distance(A - Point(v), A - v) ≤ tol
        end
        for G in (A, B, C),
            v in vectors,
            T in types
            @test eltype(T.(v) + G) == eltype(G)
            @test eltype(G + T.(v)) == eltype(G)
            H = T.(G)
            @test eltype(v + H) == eltype(H)
            @test eltype(H + v) == eltype(H)
        end
    end

    @testset "rotation" begin
        for θ in angles,
            v in vectors
            R = rotate(+θ, I)
            Q = rotate(-θ, I)
            @test distance(R*Q, I) ≤ tol
            @test distance(Q*R, I) ≤ tol
            @test distance(rotate(θ, B)(v), (R*B)(v)) ≤ tol
            @test distance(rotate(B, θ)(v), (B*R)(v)) ≤ tol
        end
        for G in (A, B, C),
            θ in angles,
            T in types
            @test eltype(rotate(T(θ), G)) == eltype(G)
            @test eltype(rotate(G, T(θ))) == eltype(G)
            H = T.(G)
            @test eltype(rotate(T.(θ), H)) == eltype(H)
            @test eltype(rotate(H, T.(θ))) == eltype(H)
        end
    end

    @testset "intercept" begin
        for M in (I, A, B)
            x, y = intercept(M)
            @test distance(M(x, y), (0,0)) ≤ tol
            P = intercept(Point, M)
            @test distance(M*P, Point(0,0)) ≤ tol
        end
    end

    @testset "show" begin
        for M in (I, A, B, C)
            @test occursin(r"AffineTransform{Float64}", string(M))
        end
    end
end

end # module


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

using StaticArrays

const _x1 = Grids._x1
const _x2 = Grids._x2
const _x3 = Grids._x3

@doc """
    knl_min_neighbor_distance!(::Val{N}, ::Val{dim}, direction,
                             min_neighbor_distance, vgeo, topology.realelems)

Computational kernel: Computes the minimum physical distance between node
neighbors within an element.

The `direction` in the reference element controls which nodes are considered
neighbors.
""" knl_min_neighbor_distance!
@kernel function knl_min_neighbor_distance!(
    ::Val{N},
    ::Val{dim},
    direction,
    min_neighbor_distance,
    vgeo,
    elems,
) where {N, dim}

    @uniform begin
        FT = eltype(min_neighbor_distance)
        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq

        if direction isa EveryDirection
            mininξ = (true, true, true)
        elseif direction isa HorizontalDirection
            mininξ = (true, dim == 2 ? false : true, false)
        elseif direction isa VerticalDirection
            mininξ = (false, dim == 2 ? true : false, dim == 2 ? false : true)
        end
    end

    e = @index(Group, Linear)
    i, j, k = @index(Local, NTuple)

    md = typemax(FT)

    ijk = i + Nq * (j - 1) + Nq * Nq * (k - 1)
    x = SVector(vgeo[ijk, _x1, e], vgeo[ijk, _x2, e], vgeo[ijk, _x3, e])

    if mininξ[1]
        @unroll for î in (i - 1, i + 1)
            if 1 ≤ î ≤ Nq
                îjk = î + Nq * (j - 1) + Nq * Nq * (k - 1)
                x̂ = SVector(
                    vgeo[îjk, _x1, e],
                    vgeo[îjk, _x2, e],
                    vgeo[îjk, _x3, e],
                )
                md = min(md, norm(x - x̂))
            end
        end
    end

    if mininξ[2]
        @unroll for ĵ in (j - 1, j + 1)
            if 1 ≤ ĵ ≤ Nq
                iĵk = i + Nq * (ĵ - 1) + Nq * Nq * (k - 1)
                x̂ = SVector(
                    vgeo[iĵk, _x1, e],
                    vgeo[iĵk, _x2, e],
                    vgeo[iĵk, _x3, e],
                )
                md = min(md, norm(x - x̂))
            end
        end
    end

    if mininξ[3]
        @unroll for k̂ in (k - 1, k + 1)
            if 1 ≤ k̂ ≤ Nqk
                ijk̂ = i + Nq * (j - 1) + Nq * Nq * (k̂ - 1)
                x̂ = SVector(
                    vgeo[ijk̂, _x1, e],
                    vgeo[ijk̂, _x2, e],
                    vgeo[ijk̂, _x3, e],
                )
                md = min(md, norm(x - x̂))
            end
        end
    end

    min_neighbor_distance[ijk, e] = md
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using SymPy
using Test

@testset "Logical" begin
    ## More logical expressions
    @vars x y
    (x ≪ 0) ∧ Le(x*y,1) ∨ (x ⩵ y) ∧ (¬(x ≫ 3))
end


text:
using Multroot
using Base.Filesystem
using Polynomials
using LinearAlgebra

import Multroot: diffmult, diffzeros

using Test

# assuming pwd = package/test
const ts = "testsuit"

function runtest(fun::Function)
    println("testing multroot with ", fun)
    @testset "$fun" begin
        p, z, job_exp, bkerr_exp, mult_exp, znorm_exp = wrap_fun(fun)
        zz, bkerr, pjcond, job = multroot(p)
        z = sort(z)
        zz = sort(z)
        job == 0 && return
        @test job == job_exp
        @test bkerr <= bkerr_exp
        @test length(zz.mult) == length(z.mult)
        lc = length(zz.mult) == length(z.mult)
        @test lc && diffmult(zz, z) <= mult_exp
        @test lc && diffzeros(zz, z) <= znorm_exp
    end
end

function wrap_fun(fun::Function, args...)
    p, z, = fun(args...)
    pp  = p isa Poly ? reverse(p.a) : p
    job_exp = 1
    bkerr_exp = Inf
    mult_exp = 0
    znorm_exp = norm(z.z, Inf) * 1e-8
    pp, z, job_exp, bkerr_exp, mult_exp, znorm_exp
end

functpath, io = mktemp(cleanup = true)
println(io, "const functions = Function[]")
for fname in filter(x -> endswith(x, "jl"), readdir(ts))
    fstring = read(abspath(ts, fname), String)
    if endswith(fname, ".m.jl")
        println(io, "fun = ")
        write(io, fstring)
        println(io)
        println(io, "push!(functions, fun)")
    else
        println(io)
        write(io, fstring)
        println(io)
    end
end
close(io)
include(functpath)

@testset "testsuite" begin
    for fun in functions
        runtest(fun)
    end
end


text:
resp_file_0 = path*"/SampleFiles/SEED/YUK7.RESP"
resp_file_1 = path*"/SampleFiles/SEED/RESP.cat"
resp_file_2 = path*"/SampleFiles/SEED/RESP.*"
rtol = eps(Float32)

printstyled("  SEED RESP\n", color=:light_green)
printstyled("    single-record file\n", color=:light_green)
S = read_seed_resp(resp_file_0, units=true)

printstyled("    multi-record file\n", color=:light_green)
S = read_seed_resp(resp_file_1, units=true)

# Channel 1 =================================================================
# Station info
@test S.id[1] == "CO.HAW.00.HHZ"
@test S.units[1] == "m/s"

R = S.resp[1]
for f in fieldnames(MultiStageResp)
  @test length(getfield(R, f)) == 11
end

# First stage
@test typeof(R.stage[1]) == PZResp64
@test isapprox(R.stage[1].a0, 5.71404E+08, rtol=rtol)
@test isapprox(R.stage[1].f0, 1.0)
@test length(R.stage[1].z) == 2
@test length(R.stage[1].p) == 5
@test isapprox(R.stage[1].z[1], 0.0+0.0im)
@test isapprox(R.stage[1].z[2], 0.0+0.0im)
@test isapprox(R.stage[1].p[1], -3.70080E-02 + 3.70080E-02im)
@test isapprox(R.stage[1].p[2], -3.70080E-02 - 3.70080E-02im)
@test isapprox(R.stage[1].p[3], -5.02650E+02 + 0.0im)
@test isapprox(R.stage[1].p[4], -1.00500E+03 + 0.0im)
@test isapprox(R.stage[1].p[5], -1.13100E+03 + 0.0im)

# Second-to-last stage
@test typeof(R.stage[10]) == CoeffResp
@test length(R.stage[10].b) == 251
@test length(R.stage[10].a) == 0
@test isapprox(R.stage[10].b[1:5], [+2.18133E-08, +1.07949E-07, +2.97668E-07, +6.73280E-07, +1.29904E-06])
@test R.fac[10] == 5
@test R.os[10] == 0
@test R.delay[10] ≈ 2.5000E-01
@test R.corr[10] ≈ 2.5000E-01
@test R.gain[10] ≈ 1.0
@test R.fg[10] ≈ 1.0

# Last stage
@test typeof(R.stage[11]) == CoeffResp
@test length(R.stage[11].b) == 251
@test length(R.stage[11].a) == 0
@test isapprox(R.stage[11].b[end-4:end], [-2.22747E-02, -1.03605E-01, +2.25295E-02, +3.17473E-01, +4.77384E-01])
@test R.fac[11] == 2
@test R.os[11] == 0
@test R.delay[11] ≈ 1.25
@test R.corr[11] ≈ 1.25
@test R.gain[11] ≈ 1.0
@test R.fg[11] ≈ 1.0
@test S.fs[1] == R.fs[end]/R.fac[end]

# Channel 5 =================================================================
@test S.id[5] == "PD.NS04..HHZ"
@test S.units[5] == "m/s"
@test S.gain[5] ≈ 1.7e8

# Total stages
R = S.resp[5]
for f in fieldnames(MultiStageResp)
  @test length(getfield(R, f)) == 4
end

# First stage
@test typeof(R.stage[1]) == PZResp64
@test isapprox(R.stage[1].a0, 1.4142, rtol=rtol)
@test isapprox(R.stage[1].f0, 1.0)
@test length(R.stage[1].z) == 2
@test length(R.stage[1].p) == 2
@test isapprox(R.stage[1].z[1], 0.0+0.0im)
@test isapprox(R.stage[1].z[2], 0.0+0.0im)
@test isapprox(R.stage[1].p[1], -4.442900E+00 + 4.442900E+00im)
@test isapprox(R.stage[1].p[2], -4.442900E+00 - 4.442900E+00im)

# Second stage
@test typeof(R.stage[2]) == CoeffResp
@test isapprox(R.fs[2], 3.000000E+03)
@test R.fac[2] == 1
@test R.os[2] == 0
@test R.delay[2] ≈ 0.0
@test R.corr[2] ≈ 0.0
@test R.gain[2] ≈ 1.0e6
@test R.fg[2] ≈ 0.0

# Third stage
@test typeof(R.stage[3]) == CoeffResp
@test length(R.stage[3].b) == 180
@test length(R.stage[3].a) == 0
@test isapprox(R.stage[3].b[1:3], [1.327638E-08, 4.137208E-08, 9.662694E-08])
@test R.fac[3] == 6
@test R.os[3] == 0
@test R.delay[3] ≈ 0.0
@test R.corr[3] ≈ 0.0
@test R.gain[3] ≈ 1.0
@test R.fs[3] ≈ 3000.0
@test R.fg[3] ≈ 0.0

# Fourth stage
@test typeof(R.stage[4]) == CoeffResp
@test length(R.stage[4].b) == 160
@test length(R.stage[4].a) == 0
@test isapprox(R.stage[4].b[1:3], [3.863808E-09, 2.261888E-09, -2.660399E-08])
@test R.fs[4] ≈ 500.0
@test R.fac[4] == 5
@test R.os[4] == 0
@test R.delay[4] ≈ 0.0
@test R.corr[4] ≈ 0.0
@test R.gain[4] ≈ 1.0
@test R.fg[4] ≈ 0.0

# Last Channel ===============================================================
n = lastindex(S.id)
@test S.id[n] == "XX.NS236..SHZ"
@test S.units[n] == "m/s"
@test S.gain[n] ≈ 3.450000e+02

# Total stages
R = S.resp[n]
for f in fieldnames(MultiStageResp)
  @test length(getfield(R, f)) == 2
end

# First stage
@test typeof(R.stage[1]) == PZResp64
@test isapprox(R.stage[1].a0, 1.0, rtol=rtol)
@test isapprox(R.stage[1].f0, 5.0)
@test length(R.stage[1].z) == 2
@test length(R.stage[1].p) == 2
@test isapprox(R.stage[1].z[1], 0.0 + 0.0im)
@test isapprox(R.stage[1].z[2], 0.0 + 0.0im)
@test isapprox(R.stage[1].p[1], -4.44 + 4.44im)
@test isapprox(R.stage[1].p[2], -4.44 - 4.44im)
@test R.fac[1] == 0
@test R.os[1] == 0
@test R.delay[1] ≈ 0.0
@test R.corr[1] ≈ 0.0
@test R.gain[1] ≈ 345.0
@test R.fg[1] ≈ 5.0

# Second stage
@test typeof(R.stage[2]) == CoeffResp
@test length(R.stage[2].b) == 1
@test length(R.stage[2].a) == 0
@test isapprox(R.stage[2].b[1], 1.0)
@test R.fac[2] == 1
@test R.os[2] == 0
@test R.delay[2] ≈ 0.0
@test R.corr[2] ≈ 0.0
@test R.gain[2] ≈ 1.0
@test R.fs[2] ≈ 40.0
@test R.fg[2] ≈ 5.0

# ===========================================================================
# Test that response start and end times are used correctly to determine
# which responses are kept

# These control values were extracted manually from RESP.cat
gain  = [1.007E+09, 8.647300E+08]
a0    = [62695.2, 86083]
f0    = [0.02, 0.02]
p1    = [-8.985E+01 + 0.0im, -5.943130E+01 + 0.0im]

#=
Target window 1:
B050F03     Station:     ANMO
B050F16     Network:     IU
B052F03     Location:    ??
B052F04     Channel:     BHZ
B052F22     Start date:  1989,241
B052F23     End date:    1991,023,22:25

Target window 2:
B050F03     Station:     ANMO
B050F16     Network:     IU
B052F03     Location:    ??
B052F04     Channel:     BHZ
B052F22     Start date:  1995,080,17:16
B052F23     End date:    1995,195
=#

for (i,t0) in enumerate([632188800000000,  # SeisIO.mktime(1990, 013, 0, 0, 0, 0)
                         801792000000000,  # SeisIO.mktime(1995, 150, 0, 0, 0, 0)
                         ])
  nx = 20000
  C = SeisChannel(id = "IU.ANMO..BHZ")
  C.t = [1 t0; nx 0]
  C.x = randn(nx)
  S = SeisData(C)
  read_meta!(S, "resp", resp_file_1)
  j = findid(S, "IU.ANMO..BHZ")
  @test S.gain[j] ≈ gain[i]
  @test S.resp[j].stage[1].a0 ≈ a0[i]
  @test S.resp[j].stage[1].f0 ≈ f0[i]
  @test S.resp[j].stage[1].p[1] ≈ p1[i]
end

# ===========================================================================
# Test that mutli-file read commands work
n += 1
printstyled("    multi-file read\n", color=:light_green)
S = read_meta("resp", resp_file_2, units=true)

# Channel 1, file 2 =========================================================
# Station info
@test S.id[n] == "AZ.DHL..BS1"
@test S.units[n] == "m/m"

R = S.resp[n]
for f in fieldnames(MultiStageResp)
  @test length(getfield(R, f)) == 9
end

# First stage
@test typeof(R.stage[1]) == PZResp64
@test R.fac[1] == 0
@test R.os[1] == 0
@test R.delay[1] ≈ 0.0
@test R.corr[1] ≈ 0.0
@test R.gain[1] ≈ 3.23
@test R.fg[1] ≈ 1.0
@test R.fs[1] ≈ 0.0

# Second stage
@test typeof(R.stage[2]) == CoeffResp
@test length(R.stage[2].b) == 0
@test length(R.stage[2].a) == 0
@test R.fac[2] == 1
@test R.os[2] == 0
@test R.delay[2] ≈ 0.0
@test R.corr[2] ≈ 0.0
@test R.gain[2] ≈ 5.263200E+05
@test R.fg[2] ≈ 1.0
@test R.fs[2] ≈ 1.28e5

# Third stage
@test typeof(R.stage[3]) == CoeffResp
@test length(R.stage[3].b) == 29
@test length(R.stage[3].a) == 0
@test R.stage[3].b[1:3] ≈ [2.441410E-04, 9.765620E-04, 2.441410E-03]
@test R.fac[3] == 8
@test R.os[3] == 0
@test R.delay[3] ≈ 8.750000E-04
@test R.corr[3] ≈ 8.750000E-04
@test R.gain[3] ≈ 1.0
@test R.fg[3] ≈ 1.0
@test R.fs[3] ≈ 1.28e5

# Strain channel
n = findid(S, "PB.CHL1.LM.LS1")
@test S.units[n] == "m/m"

R = S.resp[n]
for f in fieldnames(MultiStageResp)
  @test length(getfield(R, f)) == 4
end

# First stage
@test typeof(R.stage[1]) == PZResp64
@test R.stage[1].a0 ≈ 1.0
@test R.stage[1].f0 ≈ 0.0
@test isempty(R.stage[1].z)
@test isempty(R.stage[1].p)
@test R.fac[1] == 0
@test R.os[1] == 0
@test R.delay[1] ≈ 0.0
@test R.corr[1] ≈ 0.0
@test R.gain[1] ≈ 5.901050E-07
@test R.fg[1] ≈ 0.0
@test R.fs[1] ≈ 0.0

# Second stage
@test typeof(R.stage[2]) == PZResp64
@test R.stage[2].a0 ≈ 1.0
@test R.stage[2].f0 ≈ 0.0
@test length(R.stage[2].z) == 0
@test length(R.stage[2].p) == 4
@test R.stage[2].p[1:2] ≈ [-8.3E-02+0.0im, -8.4E-02+0.0im]
@test R.fac[2] == 0
@test R.os[2] == 0
@test R.delay[2] ≈ 0.0
@test R.corr[2] ≈ 0.0
@test R.gain[2] ≈ 1.0
@test R.fg[2] ≈ 0.0
@test R.fs[2] ≈ 0.0

# Third stage
@test typeof(R.stage[3]) == CoeffResp
@test length(R.stage[3].b) == 0
@test length(R.stage[3].a) == 0
@test R.fac[3] == 1
@test R.os[3] == 0
@test R.delay[3] ≈ 0.0
@test R.corr[3] ≈ 0.0
@test R.gain[3] ≈ 3.052E-04
@test R.fg[3] ≈ 0.0
@test R.fs[3] ≈ 10.0

# Fourth stage
@test typeof(R.stage[4]) == CoeffResp
@test length(R.stage[4].b) == 10
@test length(R.stage[4].a) == 1
@test R.stage[4].a[1] ≈ 1.0
for i = 1:10
  @test R.stage[4].b[i] ≈ 0.1
end
@test R.fac[4] == 10
@test R.os[4] == 4
@test R.delay[4] ≈ 0.5
@test R.corr[4] ≈ 0.0
@test R.gain[4] ≈ 1.0
@test R.fg[4] ≈ 0.0
@test R.fs[4] ≈ 10.0


text:
module Diffusion

export
  Problem,
  updatevars!,
  set_c!

using
  FFTW,
  Reexport

@reexport using FourierFlows

using FourierFlows: varsexpression
using LinearAlgebra: mul!, ldiv!

"""
    Problem(; parameters...)

Construct a constant diffusivity problem with steady or time-varying flow.
"""
function Problem(;
            nx = 128,
            Lx = 2π,
         kappa = 0,
            dt = 0.01,
       stepper = "RK4",
             T = Float64,
           dev = CPU()
  )

    grid = OneDGrid(dev, nx, Lx; T=T)
  params = Params(dev, kappa)
    vars = Vars(dev, grid)
     eqn = DiffusionEquation(dev, kappa, grid)

  FourierFlows.Problem(eqn, stepper, dt, grid, vars, params, dev)
end

struct Params{T} <: AbstractParams
  kappa::T
end

Params(dev, kappa::Number) = Params(kappa)
Params(dev, kappa::AbstractArray) = Params(ArrayType(dev)(kappa))

"""
    DiffusionEquation(dev, kappa, grid)

Returns the equation for constant diffusivity problem with diffusivity kappa and grid.
"""
function DiffusionEquation(dev::Device, kappa::T, grid) where T<:Number
  L = zeros(dev, T, grid.nkr)
  @. L = -kappa * grid.kr^2
  FourierFlows.Equation(L, calcN!, grid)
end

function DiffusionEquation(dev::Device, kappa::T, grid::AbstractGrid{Tg}) where {T<:AbstractArray, Tg}
  FourierFlows.Equation(0, calcN!, grid; dims=(grid.nkr,), T=cxtype(Tg))
end

struct Vars{Aphys, Atrans} <: AbstractVars
    c :: Aphys
   cx :: Aphys
   ch :: Atrans
  cxh :: Atrans
end

"""
    Vars(dev, grid)

Returns the vars for constant diffusivity problem on grid g.
"""
function Vars(::Dev, grid::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T grid.nx c cx
  @devzeros Dev Complex{T} grid.nkr ch cxh
  Vars(c, cx, ch, cxh)
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term for the 1D heat equation.
"""
function calcN!(N, sol, t, cl, v, p::Params{T}, g) where T<:Number
  @. N = 0
  nothing
end

function calcN!(N, sol, t, cl, v, p::Params{T}, g) where T<:AbstractArray
  @. v.cxh = im * g.kr * sol
  ldiv!(v.cx, g.rfftplan, v.cxh)
  @. v.cx *= p.kappa
  mul!(v.cxh, g.rfftplan, v.cx)
  @. N = im*g.kr*v.cxh
  nothing
end

"""
    updatevars!(v, g, sol)

Update the vars in v on the grid g with the solution in `sol`.
"""
function updatevars!(v, g, sol)
  ldiv!(v.c, g.rfftplan, deepcopy(sol))
  @. v.ch = sol
  nothing
end

updatevars!(prob) = updatevars!(prob.vars, prob.grid, prob.sol)

"""
    set_c!(prob, c)

Set the solution as the transform of `c`.
"""
function set_c!(prob, c)
  mul!(prob.sol, prob.grid.rfftplan, c)
  updatevars!(prob)
end

set_c!(prob, c::Function) = set_c!(prob, c.(prob.grid.x))

end # module


text:
function duffing(du,u,p,t)
	x = u[1]
	dx = u[2]
	du[1] = dx
	du[2] = -p[3]*dx - p[1]*x - p[2]*x^3 + p[4]*cos(p[5]*t)
end
parameters = [-1.0,1.0,0.3,0.2,1.2] #[alpha,beta,delta,gamma,omega]
tspan = (0.0,10.0)
dx0 = 0.0
x0 = 1.0
prob = ODEProblem(duffing,[x0,dx0],tspan,parameters)
sol = solve(prob)
sol = sol(0.0:0.001:10.0)
X = getindex.(sol.u,1)
Y = getindex.(sol.u,2)

plot(X,Y,label = false)

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
using ApproxFun
using LinearAlgebra
using Plots; gr()
using Distributions
using OrdinaryDiffEq

println("Fourier Stuff")
S = Fourier()
N = 512
x = points(S, N)
y = points(S, N)
D2x = Derivative(S,2)[1:n,1:n]
D2 = kron(D2x, D2x)
T = ApproxFun.plan_transform(S, N)
Ti = ApproxFun.plan_itransform(S, N)

println("Caches")
u = zeros(N,N)
v = similar(u)
u_n = similar(u)
v_n = similar(u)
û_n = similar(u)
v̂_n = similar(u)
û_1 = similar(u)
v̂_1 = similar(u)
u_d2 = zeros(N*N)
v_d2 = zeros(N*N)

println("Setting up problem parameters")
Mᵤ =  1e3
Mᵥ =  1.0
γᵤ =  1.0
γᵥ =  1.0
a = 0.5
c = 0.00125

println("Initial Conditions")
r01 = 0.01 .+ rand(Uniform(-0.01/100, +0.01/100), N, N)
r02 = 0.25 .+ rand(Uniform(-0.25/100, +0.25/100), N, N)
r0 = cat(r01, r02, dims=3)

cache_ops = (D2,T,Ti, u, v, u_n, v_n, û_n, v̂_n, û_1, v̂_1, u_d2, v_d2)
p = (Mᵤ, Mᵥ, γᵤ, γᵥ, a, c, cache_ops)
function NTF_nl!(dr̂,r̂,p,t)
    Mᵤ, Mᵥ, γᵤ, γᵥ, a, c, cache_ops = p
    D2,T,Ti, u, v, u_n, v_n, û_n, v̂_n, û_1, v̂_1 = cache_ops
    û = @view r̂[:,:,1]
    v̂ = @view r̂[:,:,2]
    dû = @view dr̂[:,:,1]
    dv̂ = @view dr̂[:,:,2]
    mul!(u, Ti, û)
    mul!(v, Ti, v̂)

    @. u_n = 2u*(u-1)*(2u-1)
    mul!(û_n, T, u_n)
    mul!(u_d2, D2, vec(û))
    @. û_1 = (û_n-γᵤ*reshape(u_d2,N,N)) # Figure out how to do this so you don't have to reshape?
    mul!(û_1, D2, û_1)
    @. u_n = u*v
    mul!(u_n, T, u_n)
    dû = Mᵤ*û_1 - a*u_n + c

    @. v_n = 2v*(v-1)*(2v-1)
    mul!(v̂_n, T, v_n)
    mul!(v_d2, D2, vec(v̂))
    @. û_1 = (v̂_n-γᵥ*reshape(v_d2, N, N))
    mul!(v̂_1, D2, v̂_1)
    @. v_n = u*v #Make efficient, need only once
    mul!(v_n, T, v_n)
    dv̂ = Mᵥ*v̂_1 - a*v_n + c
end

using Symbolics
#println("Jacobian sparsity")
#dr0 = copy(r0)
#@time jac_sparsity = Symbolics.jacobian_sparsity((dr,r)->NTF_nl!(dr,r,p,0.0),dr0,r0)


text:
# This file contains the exports for ClimateBase module.
# Having all exports at a single file is overall cleaner
# and leads to easier-to-identify public API.

# TODO: Actually move all exports here.

# Temporal
export monthday_indices, maxyearspan, daymonth, realtime_days, realtime_milliseconds,
temporal_sampling, timemean, timeagg, monthlyagg, yearlyagg, temporalrange, seasonalyagg,
season, DAYS_IN_ORBIT, HOURS_IN_ORBIT, seasonality, sametimespan

# Spatial
export transform_to_coord


text:
using Yeppp
using Test
import LinearAlgebra

for T in (Float32, Float64)
    x::Vector{T} = rand(1000)
    y::Vector{T} = rand(1000)
    z::Vector{T} = randn(1000)
    res::Vector{T} = rand(1000)

    @test Yeppp.add!(res, x, y) == x .+ y
    @test Yeppp.multiply!(res, x, y) == x .* y
    @test Yeppp.subtract!(res, x, y) == x .- y
    @test Yeppp.multiply(x, y) == x .* y
    @test Yeppp.subtract(x, y) == x .- y
    @test Yeppp.min!(res, x, y) == min.(x, y)
    @test Yeppp.max!(res, x, y) == max.(x, y)
    @test Yeppp.min(x, y) == min.(x, y)
    @test Yeppp.max(x, y) == max.(x, y)

    @test Yeppp.sum(x) ≈ sum(x)
    @test Yeppp.sumabs2(x) ≈ sum(abs2, x)
    @test Yeppp.dot(x, y) ≈ LinearAlgebra.dot(x, y)

end

x = rand(1000)
y = rand(1000)
z = randn(1000)
res = rand(1000)

@test Yeppp.exp!(res, x) ≈ exp.(x)
@test Yeppp.log!(res, x) ≈ log.(x)
@test Yeppp.sin!(res, x) ≈ sin.(x)
@test Yeppp.cos!(res, x) ≈ cos.(x)
@test Yeppp.tan!(res, x) ≈ tan.(x)

@test Yeppp.sum(x) ≈ sum(x)
@test Yeppp.sumabs(x) ≈ sum(abs, x)
@test Yeppp.sumabs2(x) ≈ sum(abs2, x)


text:
using SemanticReleaseDemo
using Test

@testset "SemanticReleaseDemo.jl" begin
    # Write your tests here.
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Criminisi(px, py)

Examplar-based inpainting based on confidence
and isophote maps. `(px, py)` is the patch size
as a tuple of integers.

## Notes

The Criminisi algorithm is defined for 2D images.

## References

Criminisi, A., Pérez, P., Toyama, K., 2004. Region Filling
and Object Removal by Examplar-based Image Inpainting.
"""
struct Criminisi <: InpaintAlgo
  px::Int
  py::Int
end

# implementation follows the notation in the paper
function inpaint_impl(img::AbstractArray{T,2}, mask::AbstractArray{Bool,2},
                      algo::Criminisi) where {T}
  # use all CPU cores in FFT
  set_num_threads(cpucores())

  # patch (or tile) size
  patchsize = algo.px, algo.py

  # rotation matrix: gradient -> isophote
  # (this is what makes the implementation 2D)
  R = [ cos(π/2) sin(π/2)
       -sin(π/2) cos(π/2)]

  # pad input arrays
  prepad  = @. (patchsize - 1) ÷ 2
  postpad = @. (patchsize    ) ÷ 2
  I = padarray(img, Pad(:symmetric, prepad, postpad))
  M = padarray(mask, Fill(false, prepad, postpad))

  # fix any invalid pixel value in masked region
  replace!(I, NaN => 0)

  # initialize filled region
  ϕ = .!M

  # initialize confidence map
  C = Float64.(ϕ)

  # inpainting frontier
  δΩ = findall(dilate(ϕ) .& .!ϕ)

  while !isempty(δΩ)
    # update confidence values in frontier
    for p in δΩ
      c = selectpatch(C, patchsize, p)
      b = selectpatch(ϕ, patchsize, p)
      C[p] = sum(c[b]) / prod(patchsize)
    end

    # isophote map
    ∇I = pointgradients(I, δΩ) * R
    nϕ = pointgradients(ϕ, δΩ)
    D  = vec(abs.(sum(∇I .* nϕ, dims=2)))
    D /= maximum(D)

    # select patch in frontier
    p  = δΩ[argmax(C[δΩ].*D)]
    ψₚ = selectpatch(I, patchsize, p)
    bₚ = selectpatch(ϕ, patchsize, p)

    # compute distance to all other patches
    Δ = convdist(I, ψₚ, bₚ)

    # only consider patches in filled region
    Δ[mask] .= Inf

    # select best candidate
    q = argmin(Δ)
    ψᵦ = selectpatch(I, patchsize, q)
    bᵦ = selectpatch(ϕ, patchsize, q)

    # paste patch and mark pixels as painted
    b = bᵦ .& .!bₚ
    ψₚ[b] .= ψᵦ[b]
    bₚ[b] .= true

    # update frontier
    δΩ = findall(dilate(ϕ) .& .!ϕ)
  end

  # return unpadded image
  start  = CartesianIndex(ntuple(i->1, ndims(img)))
  finish = CartesianIndex(size(img))
  I[start:finish]
end


text:
function default_algorithm(prob::DiffEqBase.AbstractRODEProblem;kwargs...)
  o = Dict{Symbol,Any}(kwargs)
  extra_kwargs = Any[]; alg=RandomEM() # Standard default

  alg_hints = get_alg_hints(o)

  # If adaptivity is not set and the tType is not a float, turn off adaptivity
  # Bad interaction with ForwardDiff
  #!(tType <: AbstractFloat) && (:adaptive ∉ keys(o)) && push!(extra_kwargs,:adaptive=>false)

  alg,extra_kwargs
end

max_stars_repo_name: ChrisRackauckas/DiffEq


text:
using SafeTestsets

@safetestset "solving 1 FDE" begin include("1FDE_test.jl") end

@safetestset "solving 2 parametric FDEs" begin include("2FDE_par_test.jl") end

@safetestset "solving 3 parametric FDEs with inner function" begin include("3FDE_par_test.jl") end

@safetestset "solving FDE of second order)" begin include("order2_test.jl") end


text:
using DemoPackageBXW
using Documenter

DocMeta.setdocmeta!(DemoPackageBXW, :DocTestSetup, :(using DemoPackageBXW); recursive=true)

makedocs(;
    modules=[DemoPackageBXW],
    authors="binxu.wang <binxu.wang@wustl.edu> and contributors",
    repo="https://github.com/Animadversio/DemoPackageBXW.jl/blob/{commit}{path}#{line}",
    sitename="DemoPackageBXW.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Animadversio.github.io/DemoPackageBXW.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Animadversio/DemoPackageBXW.jl",
    devbranch="main",
)


text:
@testset "Tokens" begin

    using DependencyTrees: Token, has_head

    @testset "Untyped Dependencies" begin
        sent = [
            ("The", 2),
            ("cat", 3),
            ("slept", 0),
            (".", 3)
        ]
        for (i, word) in enumerate(sent)
            form, head = word
            token = Token(form, head, id=i)
            @test token.form == form
            @test token.id == i
            @test token.head == head
        end
        @test !has_head(Token())
        @test !has_head(Token("hi"))
    end

    @testset "Typed Dependencies" begin
        sent = [
            ("The", "DT", 2),
            ("cat", "NN", 3),
            ("slept", "VBD", 0),
            (".", ".", 3)
        ]
        for (i, word) in enumerate(sent)
            form, deprel, head = word
            token = Token(form, head, deprel=deprel, id=i)
            @test token.form == form
            @test token.deprel == deprel
            @test token.head == head
            @test token.id == i
        end
    end

    @testset "Multi-headed tokens" begin
        t = Token("hi", (0, 2))
        @test has_head(t)
        @test has_head(t, 0)
        @test has_head(t, 2)
        @test !has_head(t, 1)
        @test DependencyTrees.headisroot(t)
    end

    @testset "Errors" begin
        t = Token("hi")
        @test_throws Exception t.x
        tx = Token("hi"; x="x")
        @test tx.x == "x"
    end
end


text:
import Base: isless

const DEFAULT_PRIORITY = 0

"""
    Priority(time_, priority)

Priority of events.

Comparison is first done by time, and after (if same time)
using priority value.

As in UNIX, lower priority numbers mean higher priority.
"""
struct Priority
    time_
    priority
end

function isless(p1::Priority, p2::Priority)
    if p1.time_ > p2.time_
        true
    elseif p1.time_ < p2.time_
        false
    else
        p1.priority > p2.priority
    end
end


text:
# Autogenerated wrapper script for spglib_jll for armv6l-linux-musleabihf
export libsymspg

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("spglib")
JLLWrappers.@declare_library_product(libsymspg, "libsymspg.so.1")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsymspg,
        "lib/libsymspg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using ParametricSurfaces
using StaticArrays

@testset "Intersection and composite surfaces" begin
    clear_entities!()
    l1   = line((0,0),(1,0));
    l2   = line((1,0),(1,1));
    l3   = line((1,1),(0,1));
    l4   = line((0,1),(0,0));
    l5   = line((1,1),(1,2))
    l6   = line((1,2),(0,2))
    l7   = line((0,2),(0,1))
    ent1 = ElementaryEntity(;boundary=[l1,l2,l3,l4])
    ent2 = ElementaryEntity(;boundary=[l5,l6,l7,flip_normal(l3)])
    ent3 = ElementaryEntity(;boundary=[l1,l2,l5,l6,l7,l4])
    Ω1 = Domain(ent1)
    Ω2 = Domain(ent2)
    Ω  = union(Ω1,Ω2)
    # intersect. When entities intersect, return them. When they don't, recurse
    # on the boundaries
    @test intersect(Ω,Ω1) == Ω1
    @test intersect(Ω,Ω2) == Ω2
    Γ12 = intersect(Ω1,Ω2)
    Γ21 = intersect(Ω2,Ω1)
    @test Γ12 == Domain(l3)
    @test Γ12 == Γ21
    # test that the orientation is different for the entity in Γ12 and Γ21
    @test normal(Γ12[1],0.5) == SVector(0,1)
    @test normal(Γ21[1],0.5) == SVector(0,-1)
    # check that skeleton includes the interior boundary l3
    Σ  = skeleton(Ω)
    @test l3 ∈ Σ
    @test l3 ∈   internal_boundary(Ω)
    @test !(l3 ∈ external_boundary(Ω))
end

@testset "curves" begin
    f   = (x) -> SVector(cos(x[1]),sin(x[1]))
    d   = HyperRectangle(0,2π)
    ent = ParametricEntity(f,d)
    s  =  SVector(rand())
    @test ent(s) == f(s)
    jac = jacobian(ent,s)
    @test jac[1] ≈ -sin(s[1]) && jac[2] ≈ cos(s[1])
    @test ent(s) ≈ normal(ent,s)
end

@testset "surfaces" begin
    f      = (x) -> SVector(x[1],x[2],sin(x[1]))
    d      = HyperRectangle((0,0),(1,1))
    ent    = ParametricEntity(f,d)
    s      = SVector{2}(rand(2))
    @test ent(s) == f(s)
    @test jacobian(ent,s) ≈ [1 0; 0 1; cos(s[1]) 0]
end

@testset "Parametric entity tests" begin
    @testset "curves" begin
        f   = (x) -> SVector(cos(x[1]),sin(x[1]))
        d   = ReferenceLine()
        ent = ParametricEntity(f,d)
        s  =  SVector(rand())
        @test ent(s) == f(s)
        jac = jacobian(ent,s)
        @test jac[1] ≈ -sin(s[1]) && jac[2] ≈ cos(s[1])
        @test ent(s) ≈ normal(ent,s)
    end
    @testset "surfaces" begin
        f      = (x) -> SVector(x[1],x[2],sin(x[1]))
        d      = ReferenceSquare()
        ent    = ParametricEntity(f,d)
        s      =  SVector{2}(rand(2))
        @test ent(s) == f(s)
        @test jacobian(ent,s) ≈ [1 0; 0 1; cos(s[1]) 0]
    end
end

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
using CUDAdrv, CUDAnative
using Base.Test

function kernel_vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    c[i] = a[i] + b[i]

    return nothing
end

dev = CuDevice(0)
ctx = CuContext(dev)

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)

d_a = CuArray(a)
d_b = CuArray(b)
d_c = similar(d_a)

len = prod(dims)
@cuda (1,len) kernel_vadd(d_a, d_b, d_c)
c = Array(d_c)
@test a+b ≈ c

destroy(ctx)


text:
#
println("Ag) Test of the QED model corrections to the level structure of atoms and ions.")
#
if   true
wa = Atomic.Computation(Atomic.Computation(), name="QED estimates for xenon-like Sn", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(50.),
                        configs=[Configuration("1s^2 2s^2 2p^6")],
                        asfSettings=AsfSettings(true, false, Basics.DFSField(), "hydrogenic", Dict{Subshell, Orbital}(), [1],    40, 1.0e-6, JAC.Subshell[], JAC.Subshell[],
                                                true, false, QedPetersburg(), LSjjSettings(false),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

elseif  false
wa = Atomic.Computation(Atomic.Computation(), name="QED estimates for hydrogen-like Sn", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(50.),
                        configs=[Configuration("1s"), Configuration("2s"), Configuration("2p")],
                        asfSettings=AsfSettings(true, false, Basics.DFSField(), "hydrogenic", Dict{Subshell, Orbital}(), [1],    0, 1.0e-6, JAC.Subshell[], JAC.Subshell[],
                                                true, false, QedPetersburg(),  LSjjSettings(false),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

end

wb = perform(wa)

# NoneQed() QedPetersburg()
#


text:
#!/usr/bin/env julia

import JuMP.optimize!

# We need this structure to prevent memory allocations while
# the shortest path algorithm runs. Other option would be
# named tuples, but I think struct makes more sense here.
mutable struct Subproblem
    inst::Instance
    depot_id::Int

    dist::Vector{Float64}
    pred::Vector{Int}

    inqueue::Vector{Bool}
    queue::Vector{Int}

    _sol_path::Vector{Int}
    _sol_cost::Int

    function Subproblem(inst::Instance, depot_id::Int)
        new(
            inst,
            depot_id,

            fill(Inf64, inst.num_trips + 2),
            fill(-1, inst.num_trips + 2),

            fill(false, inst.num_trips + 2),
            [],

            [],
            0.0
        )
    end
end

function optimize!(prob::Subproblem, depot_dual::Float64, task_duals::Vector{Float64})
    # Artificial nodes representing source and sink
    O = prob.inst.num_trips + 1
    D = prob.inst.num_trips + 2

    # Puts data structures to initial state.
    fill!(prob.dist, Inf64)
    fill!(prob.pred, -1)
    fill!(prob.inqueue, false)
    empty!(prob.queue)

    # This is only used for debugging purpose, to see if the algorithm
    # entered a cicle of negative cost.
    cnt = zeros(Int, prob.inst.num_trips + 2)

    # Set the initial state of each task, departing from the
    # depot when the path is feasible.
    for i in 1:prob.inst.num_trips
        cost = pullout_cost(prob.inst, prob.depot_id, i)
        if cost != -1
            prob.dist[i] = cost - depot_dual
            prob.pred[i] = O
            cnt[i] += 1
            prob.inqueue[i] = true
            push!(prob.queue, i)
        end
    end

    # Repeats the relaxation step until nothing changes in the path.
    while !isempty(prob.queue)
        v = pop!(prob.queue)
        prob.inqueue[v] = false
        idual = task_duals[v]

        # Relax all deadheading arcs
        for i in prob.inst.dh_forward[v]
            succ = i.succ
            cost = i.cost - idual

            if prob.dist[v] + cost < prob.dist[succ]
                prob.dist[succ] = prob.dist[v] + cost
                prob.pred[succ] = v

                if !prob.inqueue[succ]
                    push!(prob.queue, succ)
                    prob.inqueue[succ] = true
                    cnt[succ] += 1

                    if cnt[succ] >= prob.inst.num_trips + 2
                        @error "Negative cycle detected in subproblem of depot #$(prob.depot_id)."
                    end
                end
            end
        end

        # Relax the arc that conects the current task to the sink node
        cost = pullin_cost(prob.inst, prob.depot_id, v)
        if cost != -1
            cost -= idual

            if prob.dist[v] + cost < prob.dist[D]
                prob.dist[D] = prob.dist[v] + cost
                prob.pred[D] = v
            end
        end
    end

    # After solving the shortest path problem, extracts the path
    # from the pred list, and compute the original column cost.

    # Cleans the old solution if any.
    empty!(prob._sol_path)
    prob._sol_cost = 0.0

    # Only extracts the path if the cost to sink node is not too close to zero.
    if prob.dist[D] <= -0.0000001
        last = prob.pred[D]
        prob._sol_cost += pullin_cost(prob.inst, prob.depot_id, last)
        append!(prob._sol_path, last)

        while prob.pred[last] != O
            prob._sol_cost += deadhead_cost(prob.inst, prob.pred[last], last)
            append!(prob._sol_path, prob.pred[last])

            last = prob.pred[last]
        end
        prob._sol_cost += pullout_cost(prob.inst, prob.depot_id, last)
    end

    return (;
        cost=prob._sol_cost,
        path=prob._sol_path
    )
end

max_stars_repo_name: afkummer/mdvsp-colgen-study

max_issues_repo_name: afkummer/mdvsp-colgen-study

max_forks_repo_name: afkummer/mdvsp-colgen-study


text:
using Documenter
using Literate
using Photosynthesis

# define default docs pages
pages = Any[
    "Home" => "index.md",
    "API"  => "API.md"
]

# add example pages
gen_example = true;
gen_dir     = joinpath(@__DIR__, "src/generated");
rm(gen_dir, force=true, recursive=true);
mkpath(gen_dir);

if gen_example
    # array of example pages
    ex_pages = Any[];
    for _ex in ["temperature", "parasets", "photosynthesis"]
        filename    = joinpath(@__DIR__, "src/examples/$(_ex).jl");
        script      = Literate.script(filename, gen_dir);
        code        = strip(read(script, String));
        mdpost(str) = replace(str, "@__CODE__" => code);
        Literate.markdown(filename, gen_dir, postprocess=mdpost);
        push!(ex_pages, "$(uppercasefirst(_ex))" => "generated/$(_ex).md");
    end
    # add example pages to pages
    push!(pages, "Examples" => ex_pages);
end

@show pages;

# format the docs
mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
    assets = ["assets/favicon.ico"]
)

# build the docs
makedocs(
    sitename = "Photosynthesis",
    format = format,
    clean = false,
    modules = [Photosynthesis],
    pages = pages,
)

# deploy the docs to Github gh-pages
deploydocs(
    repo = "github.com/Yujie-W/Photosynthesis.jl.git",
    target = "build",
    devbranch = "main",
    push_preview = true,
)


text:
using BinaryBuilder

name = "SCALAPACK"
version = v"2.1.0"

sources = [
  ArchiveSource("http://www.netlib.org/scalapack/scalapack-2.1.0.tgz",
                "61d9216cf81d246944720cfce96255878a3f85dec13b9351f1fa0fd6768220a6"),
  DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${libdir}
cd $WORKSPACE/srcdir/scalapack-2.1.0

# the patch prevents running foreign executables, which fails on most platforms
# we instead set CDEFS manually below
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
  atomic_patch -p1 ${f}
done

CMAKE_FLAGS=(-DCMAKE_INSTALL_PREFIX=${prefix}
             -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}"
             -DCMAKE_BUILD_TYPE=Release
             -DBUILD_SHARED_LIBS=ON)

if [[ "${target}" == i686-*  ]] || [[ "${target}" == x86_64-*  ]]; then
  CMAKE_FLAGS+=(-DCMAKE_EXE_LINKER_FLAGS="-lgfortran -lquadmath")
else
  if [[ "${target}" == powerpc64le-linux-gnu ]]; then
    # special case for CMake to discover MPI_Fortran
    CMAKE_FLAGS+=(-DCMAKE_EXE_LINKER_FLAGS="-lgfortran -L/opt/${target}/${target}/sys-root/usr/lib64 -lpthread -lrt" \
                  -DCMAKE_SHARED_LINKER_FLAGS="-lgfortran -L/opt/${target}/${target}/sys-root/usr/lib64 -lpthread -lrt" \
                  -DMPI_Fortran_LINK_FLAGS="-Wl,-rpath -Wl,/workspace/destdir/lib -Wl,--enable-new-dtags -L/workspace/destdir/lib -Wl,-L/opt/${target}/${target}/sys-root/usr/lib64 -Wl,-lpthread -Wl,-lrt")
  else
    CMAKE_FLAGS+=(-DCMAKE_EXE_LINKER_FLAGS="-lgfortran")
  fi
fi

OPENBLAS=(-lopenblas)
FFLAGS=(-cpp -ffixed-line-length-none)

if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
  OPENBLAS=(-lopenblas64_)
  if [[ "${target}" == powerpc64le-linux-gnu ]]; then
    OPENBLAS+=(-lgomp)
  fi
  syms=(caxpy cbdsqr ccopy cdotc cdotu cgbmv cgbtrf cgemm cgemv cgerc cgeru cgetrf cgetrs chbmv chemm chemv cher cher2 cher2k cherk chetrd clacgv clacpy cladiv clanhs clarfg clartg claset clasr classq claswp cpbtrf cpotrf cpttrf crot csbmv cscal csscal cswap csymm csyr2k csyrk ctbtrs ctrmm ctrmv ctrsm ctrsv ctrtrs dasum daxpy dsbmv dbdsqr dcopy ddot dgbmv dgbtrf dgemm dgemv dger dgetrf dgetrs dhbmv disnan dlabad dlacpy dlae2 dlaebz dlaed4 dlaev2 dlagtf dlagts dlahqr dlamc3 dlamch dlange dlanst dlanv2 dlapy2 dlapy3 dlaqr0 dlaqr1 dlaqr3 dlaqr4 dlarfg dlarfx dlarnv dlarra dlarrb dlarrc dlarrd dlarrk dlarrv dlartg dlaruv dlascl dlaset dlasq2 dlasr dlasrt dlassq dlaswp dlasy2 dnrm2 dpbtrf dpotrf dpttrf drot dscal dstedc dsteqr dsterf dswap dsymm dsymv dsyr dsyr2 dsyr2k dsyrk dsytrd dtbtrs dtrmm dtrmv dtrsm dtrsv dtrtrs dzasum dznrm2 dzsum1 icamax icmax1 idamax ieeeck ilaenv isamax izamax izmax1 lsame lsamen sasum saxpy sbdsqr scasum scnrm2 scopy scsum1 sdot sgbmv sgbtrf sgemm sgemv sger sgetrf sgetrs shbmv sisnan slabad slacpy slae2 slaebz slaed4 slaev2 slagtf slagts slahqr slamc3 slamch slange slanst slanv2 slapy2 slapy3 slaqr0 slaqr1 slaqr3 slaqr4 slarfg slarfx slarnv slarra slarrb slarrc slarrd slarrk slarrv slartg slaruv slascl slaset slasq2 slasr slasrt slassq slaswp slasy2 snrm2 spbtrf spotrf spttrf srot ssbmv sscal sstedc ssteqr ssterf sswap ssymm ssymv ssyr ssyr2 ssyr2k ssyrk ssytrd stbtrs strmm strmv strsm strsv strtrs xerbla zaxpy zbdsqr zcopy zdotc zdotu zdscal zgbmv zgbtrf zgemm zgemv zgerc zgeru zgetrf zgetrs zhbmv zhemm zhemv zher zher2 zher2k zherk zhetrd zlacgv zlacpy zladiv zlanhs zlarfg zlartg zlaset zlasr zlassq zlaswp zpbtrf zpotrf zpttrf zrot zdbmv zscal zswap zsymm zsyr2k zsyrk ztbtrs ztrmm ztrmv ztrsm ztrsv ztrtrs)
  for sym in ${syms[@]}
  do
    FFLAGS+=("-D${sym}=${sym}_64")
    FFLAGS+=("-D${sym}_=${sym}_64_")  # due to some evil #defines in SCALAPACK
    FFLAGS+=("-D${sym^^}=${sym}_64")
  done

  CMAKE_FLAGS+=(-DCMAKE_Fortran_FLAGS=\"${FFLAGS[*]}\" \
                -DCMAKE_C_FLAGS=\"${FFLAGS[*]}\")
fi

CMAKE_FLAGS+=(-DBLAS_LIBRARIES=\"${OPENBLAS[*]}\" \
              -DLAPACK_LIBRARIES=\"${OPENBLAS[*]}\")

export CDEFS="Add_"

mkdir build && cd build
cmake .. "${CMAKE_FLAGS[@]}"

make -j${nproc} all
make install
"""

# OpenMPI and MPICH are not precompiled for Windows
# Can't get the code to build for PowerPC with libgfortran3
platforms = expand_gfortran_versions(filter!(p -> !isa(p, Windows) && arch(p) != :powerpc64le, supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libscalapack", :libscalapack),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("MPICH_jll"),
    Dependency("OpenBLAS_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using ECC, Test

tests = ["helper", "primefield", "infinity", "point", "ecc"]

for t ∈ tests
  include("$(t)_tests.jl")
end


text:
module Stheno

    using Distributions, Distances, BlockArrays, FillArrays, IterTools, Statistics, Random

    const AV{T} = AbstractVector{T}
    const AM{T} = AbstractMatrix{T}
    const AVM{T} = AbstractVecOrMat{T}

    # Various bits of utility that aren't inherently GP-related.
    include("util/covariance_matrices.jl")
    include("util/woodbury.jl")
    include("util/block_arrays.jl")
    include("util/abstract_data_set.jl")
    include("util/toeplitz.jl")
    include("util/eachindex_util.jl")
    include("util/io.jl")

    # All mean function and kernel related functionality.
    include("mean_and_kernel/mean.jl")
    include("mean_and_kernel/kernel.jl")
    include("mean_and_kernel/compose.jl")
    include("mean_and_kernel/finite.jl")
    include("mean_and_kernel/conditional.jl")
    include("mean_and_kernel/block.jl")
    # include("mean_and_kernel/transform.jl")
    include("mean_and_kernel/input_transform.jl")
    include("mean_and_kernel/degenerate.jl")
    include("mean_and_kernel/derivative.jl")
    include("mean_and_kernel/zero.jl")
    include("mean_and_kernel/printing.jl")

    # Basic Gaussian process definitions.
    include("gp/abstract_gp.jl")
    include("gp/gp.jl")
    include("gp/block_gp.jl")

    # Affine transformations of GPs.
    include("linops/indexing.jl")
    include("linops/addition.jl")
    include("linops/product.jl")
    include("linops/compose.jl")
    include("linops/project.jl")
    include("linops/conditioning.jl")
    include("linops/gradient.jl")
    # include("linops/integrate.jl")

    # The @model macro.
    include("util/model.jl")
end # module