text:
# # Solving the heat equation in soil

# This tutorial shows how to use CliMA code to solve the heat
# equation in soil.
# For background on the heat equation in general,
# and how to solve it using CliMA code, please see the
# [`heat_equation.jl`](../../Heat/heat_equation.md)
# tutorial.

# The version of the heat equation we are solving here assumes no
# sources or sinks and no diffusion of liquid water. It takes the form

# ``
# \frac{∂ ρe_{int}}{∂ t} =  ∇ ⋅ κ(θ_l, θ_i; ν, ...) ∇T
# ``

# Here

# ``t`` is the time (s),

# ``z`` is the location in the vertical (m),

# ``ρe_{int}`` is the volumetric internal energy of the soil (J/m^3),

# ``T`` is the temperature of the soil (K),

# ``κ`` is the thermal conductivity (W/m/K),

# ``ϑ_l`` is the augmented volumetric liquid water fraction,

# ``θ_i`` is the volumetric ice fraction, and

# ``ν, ...`` denotes parameters relating to soil type, such as porosity.

# We will solve this equation in an effectively 1-d domain with ``z ∈ [-1,0]``,
# and with the following boundary and initial conditions:

# ``T(t=0, z) = 277.15^\circ K``

# ``T(t, z = 0) = 288.15^\circ K ``

# `` -κ ∇T(t, z = -1) = 0 ẑ``

# The temperature ``T`` and
# volumetric internal energy ``ρe_{int}`` are related as

# ``
# ρe_{int} = ρc_s (θ_l, θ_i; ν, ...) (T - T_0) - θ_i ρ_i LH_{f0}
# ``

# where

# ``ρc_s`` is the volumetric heat capacity of the soil (J/m^3/K),

# ``T_0`` is the freezing temperature of water,

# ``ρ_i`` is the density of ice (kg/m^3), and

# ``LH_{f0}`` is the latent heat of fusion at ``T_0``.

# In this tutorial, we will use a [`PrescribedWaterModel`](@ref
# ClimateMachine.Land.PrescribedWaterModel). This option allows
# the user to specify a function for the spatial and temporal
# behavior of `θ_i` and `θ_l`; it does not solve Richard's equation
# for the evolution of moisture. Please see the tutorials
# in the `Soil/Coupled/` folder or the `Soil/Water/`
# folder for information on solving
# Richard's equation, either coupled or uncoupled from the heat equation, respectively.

# # Import necessary modules

# External (non - CliMA) modules
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Dierckx
using Plots
using DelimitedFiles

# CliMA Parameters
using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, ρ_cloud_ice, cp_l, cp_i, T_0, LH_f0

# ClimateMachine modules
using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Land.SoilHeatParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state
import ClimateMachine.DGMethods: calculate_dt
using ClimateMachine.ArtifactWrappers

# # Preliminary set-up

# Get the parameter set, which holds constants used across CliMA models.
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet();
# Initialize and pick a floating point precision.
ClimateMachine.init()
FT = Float32;
# Load a function that will create an interpolation of the
# simulation output, to be used in plotting:
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "tutorials",
    "Land",
    "Soil",
    "interpolation_helper.jl",
));
# # Determine soil parameters

# Below are the soil component fractions for various soil
# texture classes,  from Cosby et al. (1984) [1, 5].
# Note that these fractions are volumetric fractions, relative
# to other soil solids, i.e. not including pore space. These are denoted `ν_ss_i`; the CliMA
# Land Documentation uses the symbol `ν_i` to denote the volumetric fraction
# of a soil component `i` relative to the soil, including pore space.

ν_ss_silt_array =
    FT.(
        [5.0, 12.0, 32.0, 70.0, 39.0, 15.0, 56.0, 34.0, 6.0, 47.0, 20.0] ./
        100.0,
    )
ν_ss_quartz_array =
    FT.(
        [92.0, 82.0, 58.0, 17.0, 43.0, 58.0, 10.0, 32.0, 52.0, 6.0, 22.0] ./
        100.0,
    )
ν_ss_clay_array =
    FT.(
        [3.0, 6.0, 10.0, 13.0, 18.0, 27.0, 34.0, 34.0, 42.0, 47.0, 58.0] ./
        100.0,
    )
porosity_array =
    FT.([
        0.395,
        0.410,
        0.435,
        0.485,
        0.451,
        0.420,
        0.477,
        0.476,
        0.426,
        0.492,
        0.482,
    ]);

# The soil types that correspond to array elements above are, in order,
# sand, loamy sand, sandy loam, silty loam, loam, sandy clay loam,
# silty clay loam, clay loam, sandy clay, silty clay, and clay.

# Here we choose the soil type to be sandy.
# The soil column is uniform in space and time.
soil_type_index = 1
ν_ss_minerals =
    ν_ss_clay_array[soil_type_index] + ν_ss_silt_array[soil_type_index]
ν_ss_quartz = ν_ss_quartz_array[soil_type_index]
porosity = porosity_array[soil_type_index];

# This tutorial additionally compares the output of a ClimateMachine simulation with that
# of Supplemental Program 2, Chapter 5, of Bonan (2019) [1].
#  We found this useful as it
# allows us compare results from our code against a published version.

# The simulation code of [1] employs a formalism for the thermal
# conductivity `κ` based on Johansen, 1975, [2]. It assumes
# no organic matter, and only requires the volumetric
# fraction of soil solids for quartz and other minerals.
# ClimateMachine employs the formalism of Balland and Arp (2005)[3],
# which requires the
# fraction of soil solids for quartz, gravel,
# organic matter, and other minerals. Dai (2019)[4] found
# model [3] to better match
# measured soil properties across a range of soil types.

# To compare the output of the two simulations, we set the organic
# matter content and gravel content to zero in the CliMA model.
# The remaining soil components (quartz and other minerals) match between
# the two. We also run the simulation for relatively wet soil
# (water content at 80% of porosity). Under these conditions,
# the two formulations for `κ`, though taking different functional forms,
# are relatively consistent.
# The differences between models are important
# for soil with organic material and for soil that is relatively dry.
ν_ss_om = FT(0.0)
ν_ss_gravel = FT(0.0);

# We next calculate a few intermediate quantities needed for the
# determination of the thermal conductivity [3]. These include
# the conductivity of the solid material, the conductivity
# of saturated soil, and the conductivity of frozen saturated soil.

κ_quartz = FT(7.7) # W/m/K
κ_minerals = FT(2.5) # W/m/K
κ_om = FT(0.25) # W/m/K
κ_liq = FT(0.57) # W/m/K
κ_ice = FT(2.29); # W/m/K

# The particle density of soil solids in moisture-free soil
# is taken as a constant, across soil types, as in [1].
# This is a good estimate for organic material free soil. The user is referred to
# [3] for a more general expression.
ρp = FT(2700) # kg/m^3
κ_solid = k_solid(ν_ss_om, ν_ss_quartz, κ_quartz, κ_minerals, κ_om)
κ_sat_frozen = ksat_frozen(κ_solid, porosity, κ_ice)
κ_sat_unfrozen = ksat_unfrozen(κ_solid, porosity, κ_liq);

# The thermal conductivity of dry soil is also required, but this is
# calculated internally using the expression of [3].

# The volumetric specific heat of dry soil is chosen so as to match Bonan's simulation.
# The user could instead compute this using a volumetric fraction weighted average
# across soil components.
ρc_ds = FT((1 - porosity) * 1.926e06) # J/m^3/K

# Finally, we store the soil-specific parameters and functions
# in a place where they will be accessible to the model
# during integration.
soil_param_functions = SoilParamFunctions{FT}(
    porosity = porosity,
    ν_ss_gravel = ν_ss_gravel,
    ν_ss_om = ν_ss_om,
    ν_ss_quartz = ν_ss_quartz,
    ρc_ds = ρc_ds,
    ρp = ρp,
    κ_solid = κ_solid,
    κ_sat_unfrozen = κ_sat_unfrozen,
    κ_sat_frozen = κ_sat_frozen,
);

# # Initial and Boundary conditions
# We will be using a [`PrescribedWaterModel`](@ref
# ClimateMachine.Land.PrescribedWaterModel), where the user supplies the augmented
# liquid fraction and ice fraction as functions of space and time. Since we are not
# implementing phase changes, it makes sense to either have entirely liquid or
# frozen water. This tutorial shows liquid water.

# Because the two models for thermal conductivity agree well for wetter soil, we'll
# choose that here. However, the user could also explore how they differ by choosing
# drier soil.

# Please note that if the user uses a mix of liquid and frozen water, that they must
# ensure that the total water content does not exceed porosity.
prescribed_augmented_liquid_fraction = FT(porosity * 0.8)
prescribed_volumetric_ice_fraction = FT(0.0);

# Choose boundary and initial conditions for heat that will not lead to freezing of water:
heat_surface_state = (aux, t) -> eltype(aux)(288.15)
heat_bottom_flux = (aux, t) -> eltype(aux)(0.0)
T_init = (aux) -> eltype(aux)(275.15);

# We also need to define a function `init_soil!`, which
# initializes all of the prognostic variables (here, we
# only have `ρe_int`, the volumetric internal energy).
# The initialization is based on user-specified
# initial conditions. Note that the user provides initial
# conditions for heat based on the temperature - `init_soil!` also
# converts between `T` and `ρe_int`.

function init_soil!(land, state, aux, coordinates, time)
    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, time)
    θ_l = volumetric_liquid_fraction(ϑ_l, land.soil.param_functions.porosity)
    ρc_ds = land.soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, land.param_set)

    state.soil.heat.ρe_int = volumetric_internal_energy(
        θ_i,
        ρc_s,
        land.soil.heat.initialT(aux),
        land.param_set,
    )
end;

# # Create the model structure
soil_water_model = PrescribedWaterModel(
    (aux, t) -> prescribed_augmented_liquid_fraction,
    (aux, t) -> prescribed_volumetric_ice_fraction,
);

# The boundary value problem in this case, with two spatial derivatives,
# requires a boundary condition at the top of the domain and the bottom.
# This gives four possible combinations of boundary conditions - top Dirichlet
# and bottom either Neumann or Dirichlet, or top Neumann and bottom either
# Dirichlet or Neumann. The user should set the unused fields to `nothing`
# to indicate that they do not want to supply a boundary condition of that type.
# For example, below we indicate that we are applying (and supplying!) a Dirichlet
# condition at the top of the domain, and a Neumann condition at the bottom.

soil_heat_model = SoilHeatModel(
    FT;
    initialT = T_init,
    dirichlet_bc = Dirichlet(
        surface_state = heat_surface_state,
        bottom_state = nothing,
    ),
    neumann_bc = Neumann(
        surface_flux = nothing,
        bottom_flux = heat_bottom_flux,
    ),
);

# The full soil model requires a heat model and a water model, as well as the
# soil parameter functions:
m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model);

# The equations being solved in this tutorial have no sources or sinks:
sources = ();

# Finally, we create the `LandModel`. In more complex land models, this would
# include the canopy, carbon state of the soil, etc.
m = LandModel(
    param_set,
    m_soil;
    source = sources,
    init_state_prognostic = init_soil!,
);

# # Specify the numerical details
# These include the resolution, domain boundaries, integration time,
# Courant number, and ODE solver.

N_poly = 1
nelem_vert = 100

zmax = FT(0)
zmin = FT(-1)

driver_config = ClimateMachine.SingleStackConfiguration(
    "LandModel",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m;
    zmin = zmin,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

# In this tutorial, we determine a timestep based on a Courant number (
# also called a Fourier number in the context of the heat equation).
# In short, we can use the parameters of the model (`κ` and `ρc_s`),
# along with with the size of
# elements of the grid used for discretizing the PDE, to estimate
# a natural timescale for heat transfer across a grid cell.
# Because we are using an explicit ODE solver, the timestep should
# be a fraction of this in order to resolve the dynamics.

# This allows us to automate, to a certain extent, choosing a value for
# the timestep, even as we switch between soil types.

function calculate_dt(dg, model::LandModel, Q, Courant_number, t, direction)
    Δt = one(eltype(Q))
    CFL = DGMethods.courant(diffusive_courant, dg, model, Q, Δt, t, direction)
    return Courant_number / CFL
end

function diffusive_courant(
    m::LandModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    soil = m.soil
    ϑ_l, θ_i = get_water_content(soil.water, aux, state, t)
    θ_l = volumetric_liquid_fraction(ϑ_l, soil.param_functions.porosity)
    κ_dry = k_dry(m.param_set, soil.param_functions)
    S_r = relative_saturation(θ_l, θ_i, soil.param_functions.porosity)
    kersten = kersten_number(θ_i, S_r, soil.param_functions)
    κ_sat = saturated_thermal_conductivity(
        θ_l,
        θ_i,
        soil.param_functions.κ_sat_unfrozen,
        soil.param_functions.κ_sat_frozen,
    )
    κ = thermal_conductivity(κ_dry, kersten, κ_sat)
    ρc_ds = soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, m.param_set)
    return Δt * κ / (Δx * Δx * ρc_ds)
end

t0 = FT(0)
timeend = FT(60 * 60 * 3)
Courant_number = FT(0.5) # much bigger than this leads to domain errors

solver_config = ClimateMachine.SolverConfiguration(
    t0,
    timeend,
    driver_config;
    Courant_number = Courant_number,
    CFL_direction = VerticalDirection(),
);

# # Run the integration
ClimateMachine.invoke!(solver_config);

# # Plot results and comparison data from [1]
# We can pull out the `z` and `T` values on the DG grid, but these will be
# multi-valued at boundaries. We'll create an additional cartesian grid
# , create an interpolation of the
# DG output, and evaluate on this second grid.
# `T` and `z` are stored in aux:
aux = solver_config.dg.state_auxiliary;

# Specify interpolation grid:
zres = FT(0.02)
boundaries = [
    FT(0) FT(0) zmin
    FT(1) FT(1) zmax
]
resolution = (FT(2), FT(2), zres)
thegrid = solver_config.dg.grid
intrp_brck = create_interpolation_grid(boundaries, resolution, thegrid);

# Smooth output, and look at T vs z:
iaux = interpolate_variables([(aux)], intrp_brck)
iaux = iaux[1]
z_ind = varsindex(vars_state(m, Auxiliary(), FT), :z)
iz = Array(iaux[:, z_ind, :][:])
z = Array(aux[:, z_ind, :][:])
T_ind = varsindex(vars_state(m, Auxiliary(), FT), :soil, :heat, :T)
iT = Array(iaux[:, T_ind, :][:])

plot(
    iT,
    iz,
    label = "ClimateMachine",
    ylabel = "z (m)",
    xlabel = "T (K)",
    title = "Heat transfer in sand",
)
plot!(T_init.(z), z, label = "Initial condition")
filename = "bonan_heat_data.csv"
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
bonan_dataset = ArtifactWrapper(
    joinpath(clima_dir, "tutorials", "Land", "Soil", "Artifacts.toml"),
    "bonan_soil_heat",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/99vm8q8tlyoulext6c35lnd3355tx6bu.csv",
        filename = filename,
    ),],
)
bonan_dataset_path = get_data_folder(bonan_dataset)
data = joinpath(bonan_dataset_path, filename)
ds_bonan = readdlm(data, ',')
bonan_T = reverse(ds_bonan[:, 2])
bonan_z = reverse(ds_bonan[:, 1])
bonan_T_continuous = Spline1D(bonan_z, bonan_T)
bonan_at_clima_z = bonan_T_continuous.(z)
plot!(bonan_at_clima_z, z, label = "Bonan simulation")
plot!(legend = :bottomleft)
savefig("thermal_conductivity_comparison.png")
# ![](thermal_conductivity_comparison.png)

# The plot shows that the temperature at the top of the
# soil is gradually increasing. This is because the surface
# temperature is held fixed at a value larger than
# the initial temperature. If we ran this for longer,
# we would see that the bottom of the domain would also
# increase in temperature because there is no heat
# leaving the bottom (due to zero heat flux specified in
# the boundary condition).

# # References
# [1] Bonan, G. Climate Change and Terrestrial Ecosystem Modeling (2019),
# Cambridge University Press

# [2] Johansen, O. 1975. Thermal conductivity of soils. Ph.D. thesis,
# Trondheim, Norway. Cold Regions Research and Engineering Laboratory
# Draft Translation 637, 1977, ADA 044002.

# [3] Balland, V., and P. A. Arp (2005), Modeling soil thermal
# conductivities over a wide range of conditions, J. Env. Eng. Sci., 4, 549–558.

# [4] Dai, Y., N. W. amd Hua Yuan, S. Zhang, W. Shangguan, S. Liu, X. Lu,
# and Y. Xin (2019a), Evaluation of soil thermal conductivity schemes for
# use in land surface modeling, J. Adv. Model. Earth Sys., 11, 3454–3473.

# [5] Cosby, B. J., Hornberger, G. M., Clapp, R. B., and Ginn, T. R. (1984).
# A statistical exploration of the relationships of soil moisture
# characteristics to the physical properties of soils. Water Resources
# Research, 20, 682–690.

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
function solve()
    a = readline()
    a == uppercase(a) ? "A" : "a"
end

println(solve())

ext: jl

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Plots

function PlotReflectivity()

    ε_0, e, ħ = 8.85418782e-12, 1.60217662e-19, 1.05457e-34

    σ = 3.6e7 # Ω⋅m

    # 減衰小さい場合
    ε(ω, ω_p) = 1 - ω_p^2 / ω^2
    # 減衰大きい場合
    ε_damped(ω) = 1 - (σ) / (im * ω * ε_0)

    ε_actual(ω, ω_p, σ) = 1 - ω_p^2 / (ω^2 + im * ω * ε_0 * ω_p^2 / σ)

    Reflectivity(ε) = abs2((√complex(ε) - 1) / (√complex(ε) + 1))

    x = range(0.01, stop = 2, step = 0.002)

    ħωp_eV = 15.2
    omega_p = ħωp_eV * e / ħ

    println(omega_p)

    ω_range_max = omega_p * 2
    x2 = range(omega_p / 100, stop = ω_range_max, length = 500)

    function reflectivity_nodecay(ω::Float64, ω_p::Float64)
        if ω < ω_p
            1.0
        else
            (2ω^2 - ω_p^2 - 2ω * √(ω^2 - ω_p^2)) / (2ω^2 - ω_p^2 + 2ω * √(ω^2 - ω_p^2))
        end
    end

    reflectivity_largedecay(ω) = (ω^2 * ε_0^2 / σ^2) * (√complex(1 + σ / (im * ω * ε_0)) - 1)^2 * (√complex(1 - σ / (im * ω * ε_0)) - 1)^2

    Plots.display(plot(x, Reflectivity.(ε.(x, 1.0)), label = "Reflectivity, σ = ∞ (no decay)", ylims = (0, 1.5)))
    # Plots.display(plot!(x, reflectivity_nodecay.(x, 1.0), label = "Reflectivity, σ = ∞ (no decay); after simplified"))
    Plots.display(plot!(x2 ./ omega_p, Reflectivity.(ε_damped.(x2)), xlabel = "ω/ω_p", label = "Reflectivity, σ → 0 (large decay)"))
    # Plots.display(plot!(x2 ./ omega_p, reflectivity_largedecay.(x2), label = "Reflectivity, σ → 0 (large decay); after simplified"))

    Plots.display(plot!(x2 ./ omega_p, Reflectivity.(ε_actual.(x2, omega_p, σ)), label = "Reflectivity, theoretical (without approximation)"))
    gui()
    readline()
    savefig("Reflectivity-of-Free-Electron.svg")
    savefig("Reflectivity-of-Free-Electron.png")
end

PlotReflectivity()


text:
module MapTypes

export ConformalMap, DerivativeMap, InverseMap

abstract type ConformalMap end

struct DerivativeMap{M<:ConformalMap}
  m :: M
end
function Base.show(io::IO, dm::DerivativeMap)
  println(io, "d/dζ of $(dm.m)")
end

struct InverseMap{M<:ConformalMap}
  m :: M
end
function Base.show(io::IO, minv::InverseMap)
  println(io, "Inverse of $(minv.m)")
end

end


text:
mutable struct YesPleaseGame <: Game

    # List of Players in Turn Order. UserList[1] is first player
    UserList::Vector{String}

    # Deck Remaining
    Deck::Vector{Int}

    # Cards Player i have taken
    CardsTaken::Vector{Vector{Int}}
    # Tokens Player i have remaining
    TokensLeft::Vector{Int}

    # Num Tokens on Top Card
    TokenOnCard::Int

    # CurrentPlayer
    CurrPlayer::Int
end

numPlayers(game::YesPleaseGame) = length(game.UserList)
cardsLeft(game::YesPleaseGame) = length(game.Deck) - 1 # Doesn't include top card
topCard(game::YesPleaseGame) = game.Deck[end]
getCurrHandle(game::YesPleaseGame) = game.UserList[game.CurrPlayer]
getCurrTokens(game::YesPleaseGame) = game.TokensLeft[game.CurrPlayer]
gameEnded(game::YesPleaseGame) = (length(game.Deck) == 0)

function generateDeck(BaseDeck, Len)
    TempCards = deepcopy(BaseDeck)
    Random.shuffle!(TempCards)

    while length(TempCards) > 24
        pop!(TempCards)
    end
    TempCards
end

function initialiseYP(PlayingUsers::Vector{String};
                      BaseDeck = Vector(3:35), Len = 24)
    TempList = deepcopy(PlayingUsers)
    Random.shuffle!(TempList)

    N = length(TempList)

    YesPleaseGame(TempList,
                  generateDeck(BaseDeck, Len),
                  [Vector{Int}() for i in 1:N],
                  [10 for i in 1:N],
                  0,
                  1)
end

function moveNextPlayer(game::YesPleaseGame)
    game.CurrPlayer += 1
    if game.CurrPlayer > numPlayers(game)
        game.CurrPlayer -= numPlayers(game)
    end
    nothing
end

function addCardtoPlayer(game::YesPleaseGame)
    push!(game.CardsTaken[game.CurrPlayer], topCard(game))
    sort!(game.CardsTaken[game.CurrPlayer])
    pop!(game.Deck)
    nothing
end

# Returns a String depicting the status
function rejectCard(game::YesPleaseGame, handle::String)
    gameEnded(game)                 && return "game-ended"
    (handle != getCurrHandle(game)) && return "wrong-player"
    (0      >= getCurrTokens(game)) && return "insufficient-tokens"

    game.TokenOnCard += 1
    game.TokensLeft[game.CurrPlayer] -= 1
    moveNextPlayer(game)

    return "success"
end

function acceptCard(game::YesPleaseGame, handle::String)
    gameEnded(game)                 && return "game-ended"
    (handle != getCurrHandle(game)) && return "wrong-player"

    game.TokensLeft[game.CurrPlayer] += game.TokenOnCard
    game.TokenOnCard = 0

    addCardtoPlayer(game)
    return "success"
end

function getPoints(game::YesPleaseGame, idx::Int)
    A = game.CardsTaken[idx]
    score = 0
    if !isempty(A)
        score += A[1]
    end
    score += sum([ (A[i] == A[i-1] + 1) ? 0 : A[i] for i in 2:length(A) ])
    score -= game.TokensLeft[idx]
    score
end

function gameState(game::YesPleaseGame)
    msg = Dict{String, Any}()
    msg["player-list"] = game.UserList
    msg["num-cards-left"] = cardsLeft(game)
    msg["game-ended"] = gameEnded(game)
    msg["token-on-card"] = game.TokenOnCard
    msg["tokens-left"] = game.TokensLeft
    msg["curr-player"] = game.CurrPlayer
    msg["cards-taken"] = game.CardsTaken
    msg["current-points"] = [getPoints(game,i) for i in 1:numPlayers(game)]

    if !gameEnded(game)
        msg["top-card"] = topCard(game)
    end

    msg
end


text:
module Transformers

using Flux
using Requires
using Requires: @init

export Transformer, TransformerDecoder
export Stack, @nntopo_str, @nntopo

export dataset, datafile, get_batch, get_vocab

export todevice, enable_gpu
export Gpt
export Bert

const Abstract3DTensor{T} = AbstractArray{T, 3}
const Container{T} = Union{NTuple{N, T}, Vector{T}} where N

const ϵ = Ref(1e-8)

"""
    set_ϵ(x)

set the ϵ value
"""
set_ϵ(x) = (ϵ[] = x; x)

"""
    epsilon(T)

get the ϵ value in type T
"""
epsilon(::Type{T}) where T = convert(T, ϵ[])

using CuArrays
const has_cuda = Ref(false)
const use_cuda = Ref(false)

"""
  enable_gpu(t=true)

enable gpu for todevice, disable with `enable_gpu(false)`.
"""
enable_gpu(t::Bool=true) = !has_cuda[] && t ? error("CuArrays not functional") : (use_cuda[] = t)

"""
  todevice(x)

move data to device, only when gpu is enable with `enable_gpu`, basically equal `Flux.gpu` except `AbstractArray{Int}` become `CuArray{Int}`. Otherwise equal `Flux.cpu`
"""
todevice(args...) = use_cuda[] ? togpudevice(args...) : tocpudevice(args...)

tocpudevice(x) = cpu(x)
tocpudevice(x, xs...) = (x, map(cpu, xs)...)
togpudevice(x...) = error("CuArrays not functional")

#implement batchmul, batchtril for flux
include("./fix/batchedmul.jl")
include("./fix/batched_tril.jl")

#scatter/gather with atomic ops
include("./fix/atomic.jl")
include("./fix/scatter.jl")

include("./basic/Basic.jl")
include("./stacks/Stacks.jl")
include("./datasets/Datasets.jl")

include("./pretrain/Pretrain.jl")

include("./gpt/GenerativePreTrain.jl")
include("./bert/BidirectionalEncoder.jl")

using .Basic
using .Stacks
using .Datasets
using .Pretrain
using .GenerativePreTrain
using .BidirectionalEncoder

function __init__()
  precompiling = ccall(:jl_generating_output, Cint, ()) != 0

  # we don't want to include the CUDA module when precompiling,
  # or we could end up replacing it at run time (triggering a warning)
  precompiling && return

  if !CuArrays.functional()
    # nothing to do here, and either CuArrays or one of its dependencies will have warned
  else
    has_cuda[] = true
    include(joinpath(@__DIR__, "cuda/cuda.jl"))
  end
end

end # module


text:
export hessian_sparse
"""
Returns the hessian in a sparse format.
"""
function hessian_sparse(nlp,x)
    n = length(x)
    H=hess(nlp,x)
    tempH = (H+tril(H,-1)')
    H = tempH
    return H
end


text:
@testset "io" begin

    path = joinpath(@__DIR__, "test_temp.dist")
    try
        open(path, "w") do io
            write(io, "17\n0\n4\n2\n")
        end

        global ps = FBD.read_probabilities(path)

        @test ps == AliasTable([17, 0, 4, 2])
    finally
        rm(path)
    end

    path = joinpath(@__DIR__, "test_temp.graph")
    try
        FBD.write_graph(path, [[5,3,2],[17,33,9],[0,77,-5,12,12,0]])

        open(path) do io
            @test read(io, String) == "5 3 2\n17 33 9\n0 77 -5 12 12 0\n"
        end

    finally
        rm(path)
    end

end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
load("Toivo.jl")
load("Debug.jl")

module PatternDispatch
using Toivo, Debug
export @pattern, @qpat

include(find_in_path("PatternDispatch/src/Immutable.jl"))
include(find_in_path("PatternDispatch/src/Graph.jl"))
include(find_in_path("PatternDispatch/src/Recode.jl"))
include(find_in_path("PatternDispatch/src/CodeMatch.jl"))
include(find_in_path("PatternDispatch/src/Macros.jl"))
using Graph, Recode, CodeMatch
using Macros

macro qpat(ex)
    recode(ex)
end

end # module


text:
"""
A struct that fully describes the guiding-centre motion.

The fields are:\\
`vacuum` - Determines if the path was computed assuming a vacuum or not\\
`drift` - If true, then drift motion was included when computing path\\
`energy` - The energy array of the path\\
`pitch` - The pitch array of the path\\
`r` - The r array of the path\\
`z` - The z array of the path\\
`phi` - The phi array of the path\\
`dt` - The incremental times (array) of the path
"""
struct OrbitPath{T}
    vacuum::Bool
    drift::Bool
    energy::Vector{T}
    pitch::Vector{T}
    r::Vector{T}
    z::Vector{T}
    phi::Vector{T}
    dt::Vector{T}
#    dl::Vector{T}
end

function OrbitPath(T::DataType=Float64; vacuum=false,drift=false)
    OrbitPath(vacuum,drift,T[],T[],T[],T[],T[],T[])
end

Base.length(op::OrbitPath) = length(op.r)

function Base.show(io::IO, path::OrbitPath)
    print(io, typeof(path),"\n")
    print(io, " length = $(length(path))\n")
    print(io, " vacuum = $(path.vacuum)\n")
    print(io, " drift = $(path.vacuum)")
end

"""
A struct that describes a tokamak orbit.

The fields are:\\
`coordinate` - The coordinate of the orbit. Please see coordinates.jl.\\
`class` - The class of the orbit (co-passing, trapped, potato etc).\\
`tau_p` - The poloidal transit time of the orbit.\\
`tau_t` - The toroidal transit time of the orbit.\\
`path` - The OrbitPath of the orbit.
"""
struct Orbit{T,S<:AbstractOrbitCoordinate{Float64}}
    coordinate::S
    class::Symbol
    tau_p::T
    tau_t::T
    path::OrbitPath{T}
    gcvalid::Bool
end

function Orbit(c::AbstractOrbitCoordinate{T},class=:incomplete) where {T}
    return Orbit(c, class, zero(T), zero(T), OrbitPath(T), false)
end

Base.length(o::Orbit) = length(o.path.r)

"""
A struct that fully describes the integration status of the guiding-centre path/motion.

The fields are:\\
`errcode` - The error code of the integration\\
`ri` - The initial guiding-centre vector (r0,phi0,z0,p_para0,mu0)\\
`vi` - The initial guiding-centre vector (dr,dphi,dz,dp_para,dmu)\\
`initial_dir` - The initial direction of the guiding-centre motion\\
`nr` - A process number (please see integrate() and callbacks.jl)\\
`rm` - The rm coordinate\\
`zm` - The zm coordinate\\
`pm` - The pm coordinate\\
`tm` - The time at the rm,zm coordinate\\
`tau_p` - The poloidal transit time\\
`tau_t` - The toroidal transit time\\
`poloidal_complete` - True if the guiding-centre particle has completed one orbit poloidally\\
`hits_boundary` - True if the guiding-centre particle hit the boundaries set by M (AbstractEquilibrum)\\
`class` - The type of orbit (:co-passing, :trapped, :potato etc)
"""
mutable struct GCStatus{T<:Number}
    errcode::Int
    ri::SArray{Tuple{5},T,1,5}
    vi::SArray{Tuple{5},T,1,5}
    initial_dir::Int
    nr::Int
    rm::T
    zm::T
    pm::T
    tm::T
    tau_p::T
    tau_t::T
    poloidal_complete::Bool
    hits_boundary::Bool
    class::Symbol
end

function GCStatus(T=Float64)
    z = zero(T)
    return GCStatus(1,SVector{5}(z,z,z,z,z),SVector{5}(z,z,z,z,z),0,0,z,z,z,z,z,z,false,false,:incomplete)
end

function reset!(stat::GCStatus{T}, one_transit, r_callback) where T
    stat.errcode = 1
    stat.nr = 0
    stat.rm = zero(T)
    stat.zm = zero(T)
    stat.tm = zero(T)
    stat.tau_p = zero(T)
    stat.tau_t = zero(T)
    stat.poloidal_complete = false
    stat.hits_boundary = false
    stat.class = :incomplete
    if one_transit && !r_callback
        stat.nr = 2
        stat.rm = stat.ri[1]
        stat.zm = stat.ri[3]
    end
end

"""
    gc_velocity(M, gcp, r, z, p_para, mu, vacuum, drift)

Calculate the velocity of the guiding-centre particle gcp, given the inputs. The equations
are from https://doi.org/10.1063/1.2773702, with a correction for a factor c.
"""
@inline function gc_velocity(M::AbstractEquilibrium, gcp::GCParticle, r, z, p_para, mu, vacuum::Bool, drift::Bool)
    # Phys. Plasmas 14, 092107 (2007); https://doi.org/10.1063/1.2773702
    # NOTE: Equations 13 and 15 are incorrect. Remove c factor to correct
    q = gcp.q*e0 # gcp.q is given in units of integer number of elementary charges. e0 is the elementary charge.
    m = gcp.m
    mc2 = m*c0*c0

    B, E = fields(M,r,z) # Please check Equilibrium.jl/equil.jl/fields()
    Babs = norm(B)
    bhat = B/Babs

    if drift
        gB = gradB(M,r,z) # Please check Equilibrium.jl/equil.jl/gradB()
        if !vacuum
            cB = curlB(M,r,z) #Jfield(M,r,z)*4pi*10^-7
        else
            cB = SVector{3}(0.0,0.0,0.0)
        end
        bhatXgradB = cross(bhat,gB)
        cbhat = (cB + bhatXgradB)/Babs

        Bstar = B + (p_para/q)*cbhat
    else
        Bstar = B
    end
    Bstar_para = dot(Bstar,bhat)

    gamma = sqrt(1 + (2*mu*Babs)/mc2 + (p_para/(m*c0))^2) # The relativistic gamma factor
    if drift
        gradg = (mu/(gamma*mc2))*gB
        Estar = E - mc2*gradg/q
        Xdot = ((p_para/(gamma*m))*Bstar + cross(Estar,bhat))/Bstar_para
    else
        Estar = E
        Xdot = ((p_para/(gamma*m))*Bstar)/Bstar_para
    end

    p_para_dot = dot(q*Estar,Bstar/Bstar_para)

    Xdot = cylindrical_cocos(cocos(M), Xdot[1], Xdot[2]/r, Xdot[3])
    return SVector{5}(Xdot[1], Xdot[2], Xdot[3], p_para_dot, zero(p_para_dot))
end

"""
    make_gc_ode(M, gcp, stat, vacuum, drift)

Make and return an ordinary differential equation function for the guiding-centre motion.
"""
function make_gc_ode(M::AbstractEquilibrium, gcp::GCParticle, stat::GCStatus,vacuum::Bool,drift::Bool)
    ir::Int, iphi::Int, iz::Int = cylindrical_cocos_indices(cocos(M))
    ode = function f(y,p::Bool,t)
        stat
        v_gc = gc_velocity(M, gcp, y[ir], y[iz], y[4], y[5],vacuum,drift)
    end
    return ode
end

"""
    integrate(M, gcp, phi0, ..., maxphi, debug)

Integrate the guiding-centre particle motion given the axisymmetric equilibrium M
and lots of input.

Inputs:\\
`M` - The AbstractEquilibrium struct containing all info about the tokamak equilibrium. Please see Equilibrium.jl/equil.jl.\\
`gcp` - The guiding-centre particle whose motion is to be integrated. Please see GuidingCenterParticles.jl/particles.jl.\\
`phi0` - The initial toroidal angle of the guiding-centre particle\\
`dt` - The (initial) incremental time step of the integration (s)\\
`tmin` - The starting time for the integration (s)\\
`tmax` - The stopping time for the integration (s)\\
`integrator` - The type of ODE integrator to be used. For example Tsit5(), BS3(), Vern9() etc. See DifferentialEquations.jl for full list\\
`wall` - Boundary object defined in Equilibrium.jl. The (R,z) coordinates of the wall of the tokamak (in meters)\\
`interp_dt` - Interpolate the resulting integrated path onto a path evenly spaced in time, with time step size inter_dt. If 0.0, interpolate onto max_length time points.\\
`classify_orbit` - If true, let the algorithm classify the resulting orbit\\
`one_transit` - If true, stop the integration after particles completes one poloidal transit\\
`store_path` - If true, then the orbit path will be returned. Otherwise, an empty path object is returned with the status object.\\
`max_length` - The solve() function returns path arrays that are arbitrarily long. max_length will be the length of an interpolated path that is returned instead.\\
`maxiter` - If the first adaptive and non-adaptive integration attempt fail, the algorithm will re-try with progresively smaller time steps this number of times.\\
`toa` - Stands for 'try only adaptive'. If true, non-adaptive integration will not be attempted (safe but possibly incomplete).\\
`maxiters` - The solve() function has a default maximum number of iterations of 1e5. Use maxiters to increase this number.\\
`autodiff` - Deprecated.\\
`r_callback` - If true, then the integration will be terminated when the ratio between the r direction speed and the total speed has gone to zero 20 times (prevents infinite loop)\\
`verbose` - If true, lots of output messages will be printed during execution.\\
`vacuum` - If true, then vacuum will be assumed\\
`drift` - If true, then drift effects will be included\\
`limit_phi` - If true, then the integration will be terminated when phi direction reaches the value maxphi\\
`maxphi` - Please see 'limit_phi'\\
`debug` - If true, then de-bugging mode is activated. Function will terminate and return after first adaptive integration (if failed).
"""
function integrate(M::AbstractEquilibrium, gcp::GCParticle, phi0,
                   dt, tmin, tmax, integrator, wall::Union{Nothing,Wall}, interp_dt, classify_orbit::Bool,
                   one_transit::Bool, store_path::Bool, max_length::Int, maxiter::Int, toa::Bool, maxiters::Int, autodiff::Bool,
                   r_callback::Bool,verbose::Bool, vacuum::Bool, drift::Bool, limit_phi::Bool, maxphi, debug::Bool)

    stat = GCStatus(typeof(gcp.r))
    rlims, zlims = limits(M)
    if !((rlims[1] < gcp.r < rlims[end]) && (zlims[1] < gcp.z < zlims[end]))
        verbose && @warn "Starting point outside boundary: " r0 = (gcp.r, gcp.z)
        return OrbitPath(;vacuum=vacuum,drift=drift), stat
    end

    # Initial Conditions
    mc2 = gcp.m*c0^2
    p0 = sqrt(((1e3*e0*gcp.energy + mc2)^2 - mc2^2)/(c0*c0)) # The initial particle momentum
    # move abs(pitch) away from one to avoid issues when calculating jacobian;
    # don't try to be clever by using clamp it doesn't work with the autodiff
    if abs(gcp.pitch) == 1.0
        pitch0 = sign(gcp.pitch)*prevfloat(abs(gcp.pitch))
    else
        pitch0 = gcp.pitch
    end
    p_para0 = p0*pitch0*B0Ip_sign(M) # The initial parallel momentum
    p_perp0 = p0*sqrt(1.0 - pitch0^2) # The initial perpendicular momentum
    mu_0 = (p_perp0^2)/(2*gcp.m*norm(Bfield(M,gcp.r,gcp.z))) # The initial (and presumably constant) magnetic moment

    r0 = SVector{5}(cylindrical_cocos(cocos(M), gcp.r, one(gcp.r)*phi0, gcp.z)...,
                    one(gcp.r)*p_para0, one(gcp.r)*mu_0) # The initial guiding-centre element vector

    stat.ri = r0 # The initial status vector is the initial vector

    ir, iphi, iz = cylindrical_cocos_indices(cocos(M))
    gc_ode = make_gc_ode(M,gcp,stat,vacuum,drift) # Make the ordinary differential equation
    stat.vi = gc_ode(r0,false,tmin) # Obtain the initial velocity

    r_callback = r_callback && !iszero(stat.vi[ir])

    stat.initial_dir = abs(stat.vi[ir]) > abs(stat.vi[iz]) ? ir : iz # Is it mostly in the R- or z-direction? Element number 1 or 3?

    tspan = (one(gcp.r)*tmin,one(gcp.r)*tmax)
    ode_prob = ODEProblem(gc_ode,r0,tspan,one_transit) # Define the ODE problem

    # The callbacks. Please see callbacks.jl for specification.
    if wall != nothing
        wall_cb = wall_callback(wall)
    end
    if one_transit
        if r_callback
            cb = transit_callback
        else
            stat.nr = 2
            stat.rm = r0[ir]
            stat.zm = r0[iz]
            cb = CallbackSet(pol_cb,oob_cb, brr_cb)
        end
        if wall != nothing
            cb = CallbackSet(cb.continuous_callbacks..., wall_cb, cb.discrete_callbacks...)
        end
    else
        cb = oob_cb
        if wall != nothing
            cb = CallbackSet(wall_cb,oob_cb)
        end
    end
    if limit_phi
        phi_cb = phi_callback(maxphi)
        cb = CallbackSet(cb.continuous_callbacks..., phi_cb, cb.discrete_callbacks...)
    end

    if verbose
        println(" ")
        println("---- Initiating orbit integration ---- ")
        println("- gcp: $(gcp)")
        println("- dt: $(round(dt/cyclotron_period(M,gcp),sigdigits=4)) times a cyclotron period")
        println("- (tmin,tmax)= ($(tmin),$(round(tmax/dt))*dt)")
        println("- Integrator: $(integrator)")
        println("- interp_dt: $(interp_dt)")
        println("- classify_orbit: $(classify_orbit)")
        println("- one_transit: $(one_transit)")
        println("- Maximum orbit length (will otherwise be interpolated onto this length): $(max_length)")
        println("- The first integration will use $(maxiters) numerical iterations. ")
        r_callback && println("- r_callback will be used.")
        vacuum && println("- Vacuum is assumed.")
        drift && println("- Particle drifts are included.")
        limit_phi && println("- The integration will be toroidally limited to $(round(maxphi/(2*pi))) turns.")
        println("--------------------------------")
    end

    # Always try adaptive first
    dts = dt
    success = false
    retcode = :TotalFailure
    try
        sol = solve(ode_prob, integrator, dt=dts, reltol=1e-8, abstol=1e-12, verbose=false, force_dtmin=true,
                    callback=cb,adaptive=true,maxiters=maxiters)
        if sol.t[end] >= tspan[2] && one_transit
            reset!(stat,one_transit,r_callback)
            ode_prob = remake(ode_prob; tspan=(tspan[1],100*tspan[2]))
            sol = solve(ode_prob, integrator, dt=dts, reltol=1e-8, abstol=1e-12, verbose=false, force_dtmin=true,
                        callback=cb,adaptive=true,maxiters=Int(maxiters*10))
        end
        success = (sol.retcode == :Success || sol.retcode == :Terminated) &&
                  ((stat.poloidal_complete || !one_transit) || limit_phi) || stat.hits_boundary
        retcode = sol.retcode

    catch err
        verbose && (println("Adaptive"); println(err))
        if isa(err,InterruptException)
            rethrow(err)
        end
    end

    if !success && debug
        verbose && println("Adaptive failed. You wanted to fix bugs. Here is some extra information that might help.")
        verbose && println("Success: $(success)")
        verbose && println("sol.retcode: $(retcode)")
        verbose && println("stat.poloidal_complete: $(stat.poloidal_complete)")
        verbose && println("!one_transit: $(!one_transit)")
        verbose && println("Length(sol): $(length(sol))")
        verbose && println(gcp)
        verbose && println("The integration will terminate now.")

        ir, iphi, iz = cylindrical_cocos_indices(cocos(M))
        r = getindex.(sol.u,ir)
        phi = getindex.(sol.u,iphi)
        z = getindex.(sol.u,iz)
        ppara = getindex.(sol.u,4)
        mu = getindex.(sol.u,5)
        pitch = get_pitch.(M, gcp, ppara, mu, r, z)
        n = length(sol)
        dt = eltype(r)[(sol.t[min(i+1,n)] - sol.t[i]) for i=1:n]
        energy = get_kinetic_energy.(M, gcp, ppara, mu, r, z)

        path = OrbitPath(vacuum,drift,energy,pitch,r,z,phi,dt)
        return path, stat
    end

    if !success && !toa #Try non-adaptive
        verbose && println(" ")
        verbose && println("Adaptive did not find guiding-centre orbit. Trying non-adaptive... ")
        try
            reset!(stat,one_transit,r_callback)
            sol = solve(ode_prob, integrator, dt=dts, reltol=1e-8, abstol=1e-12, verbose=false, force_dtmin=true,
                        callback=cb,adaptive=false)
            if sol.t[end] >= tspan[2] && one_transit
                reset!(stat,one_transit,r_callback)
                ode_prob = remake(ode_prob; tspan=(tspan[1],100*tspan[2]))
                sol = solve(ode_prob, integrator, dt=dts, reltol=1e-8, abstol=1e-12, verbose=false, force_dtmin=true,
                            callback=cb,adaptive=false)
            end
            success = (sol.retcode == :Success || sol.retcode == :Terminated) &&
                      (stat.poloidal_complete || !one_transit || limit_phi) || stat.hits_boundary
            retcode = sol.retcode
        catch err
            verbose && (println("Non-Adaptive: dt = $(dts)"); println(err))
            if isa(err,InterruptException)
                rethrow(err)
            end
        end
    end

    for i=1:maxiter #Try progressivly smaller time step
        (success || toa) && break
        dts = dts/10
        verbose && println("Re-trying non-adaptive with smaller timestep... ")
        try
            reset!(stat,one_transit,r_callback)
            sol = solve(ode_prob, integrator, dt=dts, reltol=1e-8, abstol=1e-12, verbose=false,
                        callback=cb, force_dtmin=true, adaptive=false)
            success = (sol.retcode == :Success || sol.retcode == :Terminated) &&
                      (stat.poloidal_complete || !one_transit || limit_phi) || stat.hits_boundary
            retcode = sol.retcode
        catch err
            verbose && (println("Non-Adaptive: dt = $(dts)"); println(err))
            if isa(err,InterruptException)
                rethrow(err)
            end
        end
    end

    if !success
        verbose && @warn "Unable to find Guiding Center Orbit" gcp retcode
        stat.class = :incomplete
        return OrbitPath(;vacuum=vacuum,drift=drift), stat
    end
    stat.errcode=0

    if one_transit && stat.class != :lost && !stat.poloidal_complete && !toa #Try one more time
        verbose && println("Trying one last time with even smaller timestep... ")
        try
            reset!(stat,one_transit,r_callback)
            sol = solve(ode_prob, integrator, dt=dts/10, reltol=1e-8, abstol=1e-12, verbose=false,
                        callback=cb, force_dtmin=true, adaptive=false)
            success = (sol.retcode == :Success || sol.retcode == :Terminated) &&
                      (stat.poloidal_complete || !one_transit) || stat.hits_boundary
            retcode = sol.retcode
        catch err
            verbose && println(err)
            if isa(err,InterruptException)
                rethrow(err)
            end
        end
    end

    if one_transit && stat.class != :lost && !stat.poloidal_complete
        verbose && println(" ")
        verbose && println("maxiters: $(maxiters)")
        verbose && println("Solution array lengths: $(length(sol.t))")
        verbose && @warn "Orbit did not complete one transit in allotted time" gcp tmax retcode
    end

    if interp_dt > 0.0
        n = floor(Int,abs(sol.t[end]/interp_dt))
        if n > max_length
            @warn "Requested interp_dt conflicts with max_length"
        end
        n = min(max_length,n)
        sol = sol(range(tmin,sol.t[end],length=n)) # Interpolate solution onto n points evently spaced in time
    else
        n = min(max_length, length(sol))
        sol = sol(range(tmin,sol.t[end],length=n)) # Interpolate solution onto n points evently spaced in time
    end
    n = length(sol)

    #Needed for classification
    ir, iphi, iz = cylindrical_cocos_indices(cocos(M))
    r = getindex.(sol.u,ir)
    phi = getindex.(sol.u,iphi)
    z = getindex.(sol.u,iz)
    ppara = getindex.(sol.u,4)
    mu = getindex.(sol.u,5)
    pitch = get_pitch.(M, gcp, ppara, mu, r, z)

    if !r_callback # Check rm, if not automatically checked in r_callback
        ind = argmax(r)
        if r[ind] > stat.rm
            stat.rm = r[ind]
            stat.zm = z[ind]
            stat.pm = pitch[ind]
        end
    end

    # Classify the orbit
    if stat.class == :unknown && classify_orbit && one_transit
        stat.class = classify(r,z,pitch,magnetic_axis(M))
    end

    if !store_path
        return OrbitPath(;vacuum=vacuum,drift=drift), stat
    end

    # Get everything else in path
    dt = eltype(r)[(sol.t[min(i+1,n)] - sol.t[i]) for i=1:n] # Get the incremental time steps
    energy = get_kinetic_energy.(M, gcp, ppara, mu, r, z) # Get the energy of the guiding-centre particle

    # P_rz
#    prz = zeros(n)
#    dl = zeros(n)
#    @inbounds for i=1:n
#        v_gc = gc_ode([r[i],phi[i],z[i]],false,0.0)
#        v = norm(v_gc[1:2:3])
##        prz[i] = 1/(T_p*v)
#        dl[i] = v*dt[i]
#    end

    path = OrbitPath(vacuum,drift,energy,pitch,r,z,phi,dt) # Create an OrbitPath object

    if stat.class == :lost # Take care of lost orbits
        stat.tau_p = zero(stat.tau_p)
        stat.tau_t = zero(stat.tau_t)
        return path, stat
    end

    return path, stat # Return the orbit path and the integration status object
end

"""
    integrate(M, gcp; phi0=0.0, ..., debug = false)

All default values for the arguments are set in this function. They are passed into the actual integrate() function
written above.
"""
function integrate(M::AbstractEquilibrium, gcp::GCParticle; phi0=0.0, dt=cyclotron_period(M,gcp)*1e-1,
                   tmin=0.0,tmax=1e6*dt, integrator=Tsit5(), wall=nothing, interp_dt = 0.0,
                   classify_orbit=true, one_transit=false, store_path=true, max_length=500,
                   maxiter=3, toa=false, maxiters=Int(1e6), autodiff=true, r_callback=false, verbose=false,
                   vacuum=false, drift=true, limit_phi=false, maxphi=10*2*pi, debug=false)

    path, stat = integrate(M, gcp, phi0, dt, tmin, tmax, integrator, wall, interp_dt,
                           classify_orbit, one_transit, store_path, max_length, maxiter,
                           toa, maxiters, autodiff, r_callback, verbose, vacuum, drift, limit_phi, maxphi, debug)
    return path, stat
end

function integrate(M::AbstractEquilibrium, c::EPRCoordinate; kwargs...)
    gcp = GCParticle(c)
    return integrate(M, gcp; kwargs...)
end

function get_orbit(M::AbstractEquilibrium, gcp::GCParticle; kwargs...)
    path, stat = integrate(M, gcp; one_transit=true, r_callback=true, kwargs...)
    gcvalid = gcde_check(M, gcp, path)

    if stat.class == :incomplete || stat.class == :lost
        return Orbit(EPRCoordinate(typeof(gcp.r);amu=(gcp.m/mass_u),q=gcp.q),stat.class,stat.tau_p,stat.tau_t,path,gcvalid)
    end
    hc = HamiltonianCoordinate(M,gcp)
    KE = get_kinetic_energy(M,hc,stat.rm,stat.zm)
    tm = (stat.tau_p - stat.tm)/stat.tau_p
    c = EPRCoordinate(KE,stat.pm,stat.rm,stat.zm;t=tm,amu=(gcp.m/mass_u),q=gcp.q)

    return Orbit(c,stat.class,stat.tau_p,stat.tau_t,path,gcvalid)
end

function get_orbit(M::AbstractEquilibrium, c::EPRCoordinate; hard = false, rtol=1e-4, kwargs...)
    gcp = GCParticle(c)
    path, stat = integrate(M, gcp; one_transit=true, r_callback=false, kwargs...)
    gcvalid = gcde_check(M, gcp, path)

    rmax = stat.rm
    if stat.class != :incomplete && stat.class != :lost
        if rmax > c.r && (hard || !isapprox(rmax,c.r,rtol=rtol))
            stat.class = :invalid
        end
    else
        stat.tau_p=zero(stat.tau_p)
        stat.tau_t=zero(stat.tau_t)
    end
    return Orbit(c,stat.class,stat.tau_p,stat.tau_t,path,gcvalid)
end

function Base.show(io::IO, orbit::Orbit)
    classes = Dict(:trapped=>"Trapped ",:co_passing=>"Co-passing ",:ctr_passing=>"Counter-passing ",
                   :stagnation=>"Stagnation ",:potato=>"Potato ",:incomplete=>"Incomplete ",
                   :Invalid=>"Invalid ",:meta=>"Meta ",:lost=>"Lost ")
    class_str = orbit.class in keys(classes) ? classes[orbit.class] : string(orbit.class)

    println(io, class_str*"Orbit:")
    @printf(io, " τ_p= %.3f μs\n", orbit.tau_p*1e6)
    @printf(io, " τ_t= %.3f μs\n", orbit.tau_t*1e6)
    print(io,   " gc-valid: $(orbit.gcvalid)")
end


text:
using ShapeDistributions
using Documenter

makedocs(;
    modules=[ShapeDistributions],
    authors="Yusri Dwi Heryanto",
    repo="https://github.com/yusri-dh/ShapeDistributions.jl/blob/{commit}{path}#L{line}",
    sitename="ShapeDistributions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://yusri-dh.github.io/ShapeDistributions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/yusri-dh/ShapeDistributions.jl",
)


text:
module ADIOS2

using ADIOS2_jll
using MPI

### Helpers

const Maybe{T} = Union{Nothing,T}
maybe(::Nothing, other) = other
maybe(x, other) = x

function free(ptr::Ptr)
    @static Sys.iswindows() ? Libc.free(ptr) :
            ccall((:free, libadios2_c), Cvoid, (Ptr{Cvoid},), ptr)
end

include("types.jl")
include("adios.jl")
include("io.jl")
include("variable.jl")
include("attribute.jl")
include("engine.jl")

include("highlevel.jl")

end


text:
function updateGraph!(jg::JGraph, ag::AbstractGraph)
end

function updateGraph!(
    jg::JGraph,
    g,
    adjacency_matrix;
    get_vertex_property=(args...)->nothing,
    get_edge_property=(args...)->nothing,
    edge_property_limits=Dict{Symbol, Any}(),
    vertex_property_limits=Dict{Symbol, Any}()
)
end

function updateGraph!(jg::JGraph, rg::ReferenceGraph)
end

function animate_inneighbors(v::GraphVertex)
end

function animate_outneighbors(v::GraphVertex)
end

function animate_neighbors(v::GraphVertex)
end

function animate_path(jg::JGraph, path::Vector{GraphEdge})
end

function animate_bfs(jg::JGraph, vertex::GraphVertex)
end

function animate_dfs(jg::JGraph, vertex::GraphVertex)
end


text:
# https://discourse.julialang.org/t/yaml-jl-custom-tag-constructor-example-include/3900

module yaml

export load_file_relative

using YAML

const _constructors = Dict{AbstractString, Function}()

"""Loads YAML files relative to the current directory"""
function load_file_relative(filename::AbstractString;
    constructors::Dict{AbstractString, Function}=_constructors
)
    filename = abspath(filename)
    cd(dirname(filename)) do
        YAML.load_file(filename, constructors)
    end
end

"""Processes an !include constructor"""
function construct_include(constructor::YAML.Constructor, node::YAML.ScalarNode)
    filename = node.value
    if !isfile(filename)
        throw(YAML.ConstructorError(nothing, nothing, "file $(abspath(filename)) does not exist", node.start_mark))
    end
    if !(last(splitext(filename)) in (".yaml", ".yml", ".json"))
        return readstring(filename)
    end
    # pass forward custom constructors
    constructors = Dict{AbstractString, Function}(
        tag => f for (tag, f) in constructor.yaml_constructors
        if tag !== nothing && ismatch(r"^!\w+$", tag)
    )
    load_file_relative(filename; constructors=constructors)
end
construct_include(::YAML.Constructor, node::YAML.Node) = throw(
    YAML.ConstructorError(nothing, nothing, "expected a scalar node, but found $(typeof(node))", node.start_mark)
)
_constructors["!include"] = construct_include

end


text:
export zoom
function zoom(h :: LineModel,
              t₀ :: Float64,
              t₁ :: Float64;
              c₁ :: Float64=0.01,
              tol :: Float64=1e-7,
              ϵ :: Float64=1e-10,
              maxiter :: Int=50,
              verbose :: Bool=false)

        if obj(h,t₀)<obj(h,t₁)
          tl=t₀
          th=t₁
        else
          tl=t₁
          th=t₀
        end

        h₀=obj(h,t₀)
        dh₀=grad(h,t₀)

        hl=obj(h,tl)
        dhl=grad(h,tl)
        hh=obj(h,th)
        dhh=grad(h,th)

        i=0

        ti= (tl+th)/2
        hi=obj(h,ti)
        dhi=grad(h,ti)

        verbose && println("on a les paramètre de zoom")

        verbose && @printf(" iter        tₗ        tₕ         tᵢ        hₗ        hₕ         hᵢ         dhᵢ\n")
        verbose && @printf(" %7.2e %7.2e  %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e\n", i,tl,th,ti,hl,hh,hi,dhi)

        while ((i<maxiter) & (abs(dhi)>tol)) || i==0
          # if i==0
          #   verbose && println("on est dans le while de zoom")
          # end
          hi=obj(h,ti)
          if (hi>h₀+c₁*(ti-t₀)*dh₀) || (hl<=hi)
            th=ti
            hh=hi
            dhh=hi
          else
            dhi=grad(h,ti)
            if (abs(dhi)<ϵ)
              topt=ti
              iter=i
              return (topt,iter)
            elseif dhi*(th-tl)>=0
              th=tl
              hh=hl
              dhh=dhl
            end
            tl=ti
            hl=hi
            dhl=dhi
          end

          ti= (tl+th)/2

          i=i+1
          verbose && @printf(" %7.2e %7.2e  %7.2e  %7.2e %7.2e %7.2e %7.2e %7.2e\n", i,tl,th,ti,hl,hh,hi,dhi)

        end
        topt=ti
        iter=i

        return (topt,iter)

end


text:
function set_velocity_tracer_fields(arch, grid, fieldname, value, answer)
    model = IncompressibleModel(architecture=arch, float_type=eltype(grid), grid=grid)
    kwarg = Dict(fieldname=>value)
    set!(model; kwarg...)

    if fieldname ∈ propertynames(model.velocities)
        ϕ = getproperty(model.velocities, fieldname)
    else
        ϕ = getproperty(model.tracers, fieldname)
    end

    return interior(ϕ) ≈ answer
end

function initial_conditions_correctly_set(arch, FT)
    model = IncompressibleModel(grid=RegularCartesianGrid(FT, size=(16, 16, 8), extent=(1, 2, 3)),
                                architecture=arch, float_type=FT)

    # Set initial condition to some basic function we can easily check for.
    # We offset the functions by an integer so that we don't end up comparing
    # zero values with other zero values. I was too lazy to pick clever functions.
    u₀(x, y, z) = 1 + x + y + z
    v₀(x, y, z) = 2 + sin(x * y * z)
    w₀(x, y, z) = 3 + y * z
    T₀(x, y, z) = 4 + tanh(x + y - z)
    S₀(x, y, z) = 5

    set!(model, u=u₀, v=v₀, w=w₀, T=T₀, S=S₀)

    Nx, Ny, Nz = model.grid.Nx, model.grid.Ny, model.grid.Nz
    xC, yC, zC = nodes(model.tracers.T)
    xF, yF, zF = nodes((Face, Face, Face), model.grid)
    u, v, w = model.velocities.u.data, model.velocities.v.data, model.velocities.w.data
    T, S = model.tracers.T.data, model.tracers.S.data

    all_values_match = true
    for i in 1:Nx, j in 1:Ny, k in 1:Nz
        values_match = ( u[i, j, k] ≈ 1 + xF[i] + yC[j] + zC[k]       &&
                         v[i, j, k] ≈ 2 + sin(xC[i] * yF[j] * zC[k])  &&
                         w[i, j, k] ≈ 3 + yC[j] * zF[k]               &&
                         T[i, j, k] ≈ 4 + tanh(xC[i] + yC[j] - zC[k]) &&
                         S[i, j, k] ≈ 5)
        all_values_match = all_values_match & values_match
    end

    return all_values_match
end

@testset "Models" begin
    @info "Testing models..."

    @testset "Model constructor errors" begin
        grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
        @test_throws TypeError IncompressibleModel(architecture=CPU, grid=grid)
        @test_throws TypeError IncompressibleModel(architecture=GPU, grid=grid)
        @test_throws TypeError IncompressibleModel(grid=grid, boundary_conditions=1)
        @test_throws TypeError IncompressibleModel(grid=grid, forcing=2)
        @test_throws TypeError IncompressibleModel(grid=grid, background_fields=3)
    end

    topos = ((Periodic, Periodic, Periodic),
             (Periodic, Periodic,  Bounded),
             (Periodic,  Bounded,  Bounded),
             (Bounded,   Bounded,  Bounded))

    for topo in topos
        @testset "$topo model construction" begin
            @info "  Testing $topo model construction..."
            for arch in archs, FT in float_types
		        arch isa GPU && topo == (Bounded, Bounded, Bounded) && continue

                grid = RegularCartesianGrid(FT, topology=topo, size=(16, 16, 2), extent=(1, 2, 3))
                model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT)

                # Just testing that the model was constructed with no errors/crashes.
                @test model isa IncompressibleModel

                # Test that the grid didn't get mangled
                @test grid === model.grid
            end
        end
    end

    @testset "Adjustment of halos in IncompressibleModel constructor" begin
        @info "  Testing adjustment of halos in IncompressibleModel constructor..."

        default_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 2, 3))
        funny_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 2, 3), halo=(1, 3, 4))

        # Model ensures that halos are at least of size 1
        model = IncompressibleModel(grid=default_grid)
        @test model.grid.Hx == 1 && model.grid.Hy == 1 && model.grid.Hz == 1

        model = IncompressibleModel(grid=funny_grid)
        @test model.grid.Hx == 1 && model.grid.Hy == 3 && model.grid.Hz == 4

        # Model ensures that halos are at least of size 2
        for scheme in (CenteredFourthOrder(), UpwindBiasedThirdOrder())
            model = IncompressibleModel(advection=scheme, grid=default_grid)
            @test model.grid.Hx == 2 && model.grid.Hy == 2 && model.grid.Hz == 2

            model = IncompressibleModel(advection=scheme, grid=funny_grid)
            @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4
        end

        # Model ensures that halos are at least of size 3
        for scheme in (WENO5(), UpwindBiasedFifthOrder())
            model = IncompressibleModel(advection=scheme, grid=default_grid)
            @test model.grid.Hx == 3 && model.grid.Hy == 3 && model.grid.Hz == 3

            model = IncompressibleModel(advection=scheme, grid=funny_grid)
            @test model.grid.Hx == 3 && model.grid.Hy == 3 && model.grid.Hz == 4
        end

        # Model ensures that halos are at least of size 2 with AnisotropicBiharmonicDiffusivity
        model = IncompressibleModel(closure=AnisotropicBiharmonicDiffusivity(), grid=default_grid)
        @test model.grid.Hx == 2 && model.grid.Hy == 2 && model.grid.Hz == 2

        model = IncompressibleModel(closure=AnisotropicBiharmonicDiffusivity(), grid=funny_grid)
        @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4
    end

    @testset "Model construction with single tracer and nothing tracer" begin
        @info "  Testing model construction with single tracer and nothing tracer..."
        for arch in archs
            grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 2, 3))

            model = IncompressibleModel(grid=grid, architecture=arch, tracers=:c, buoyancy=nothing)
            @test model isa IncompressibleModel

            model = IncompressibleModel(grid=grid, architecture=arch, tracers=nothing, buoyancy=nothing)
            @test model isa IncompressibleModel
        end
    end

    @testset "Non-dimensional model" begin
        @info "  Testing non-dimensional model construction..."
        for arch in archs, FT in float_types
            grid = RegularCartesianGrid(FT, size=(16, 16, 2), extent=(3, 2, 1))
            model = NonDimensionalModel(architecture=arch, float_type=FT, grid=grid, Re=1, Pr=1, Ro=Inf)

            # Just testing that a NonDimensionalModel was constructed with no errors/crashes.
            @test model isa IncompressibleModel
        end
    end

    @testset "Setting model fields" begin
        @info "  Testing setting model fields..."
        for arch in archs, FT in float_types
            N = (4, 4, 4)
            L = (2π, 3π, 5π)

            grid = RegularCartesianGrid(FT, size=N, extent=L)
            x, y, z = nodes((Face, Cell, Cell), grid, reshape=true)

            u₀(x, y, z) = x * y^2 * z^3
            u_answer = @. x * y^2 * z^3

            T₀ = rand(size(grid)...)
            T_answer = deepcopy(T₀)

            @test set_velocity_tracer_fields(arch, grid, :u, u₀, u_answer)
            @test set_velocity_tracer_fields(arch, grid, :T, T₀, T_answer)
            @test initial_conditions_correctly_set(arch, FT)
        end
    end
end


text:
@testset "io.jl" begin
    imc, head = load(joinpath(@__DIR__, "../example_image.fits"), IntensityMap)
    ime = ehtim.image.load_image(joinpath(@__DIR__, "../example_image.fits"))
    data = load_data()
    @test size(ime.imarr("I")) == size(imc)
    @test flux(imc) ≈ ime.total_flux()
    @test imc.fovx ≈ ime.fovx()
    save("test.fits", imc)
    rm("test.fits")
    save("test.fits", imc, last(data))
    rm("test.fits")
end


text:
module CodecBase

export
    # base 16
    Base16Encoder,
    Base16EncoderStream,
    Base16Decoder,
    Base16DecoderStream,

    # base 32
    Base32Encoder,
    Base32EncoderStream,
    Base32Decoder,
    Base32DecoderStream,

    # base 64
    Base64Encoder,
    Base64EncoderStream,
    Base64Decoder,
    Base64DecoderStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Codec,
    Memory,
    Error

macro unreachable()
    :(@assert false "unreachable")
end

include("error.jl")
include("state.jl")
include("buffer.jl")
include("codetable.jl")
include("base16/codetable.jl")
include("base16/encoder.jl")
include("base16/decoder.jl")
include("base32/codetable.jl")
include("base32/encoder.jl")
include("base32/decoder.jl")
include("base64/codetable.jl")
include("base64/encoder.jl")
include("base64/decoder.jl")

end # module


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.bipartite" begin

@testset "maxcardmatch" begin
    relation1 = Dict(
        1 => Set([1]), 2 => Set([1, 2, 3]),
        3 => Set([1, 2]), 4 => Set([2, 3, 4])
    )
    @test maxcard(Set(1:4), Set(1:4), relation1) == 4
    relation2 = Dict(
        1 => Set([2]), 2 => Set([2, 3]),
        3 => Set(1:5), 4 => Set([3]), 5 => Set([4])
    )
    @test maxcard(Set(1:5), Set(1:5), relation2) == 4
    relation3 = Dict(
        1 => Set([1, 2]), 2 => Set([3]),
        3 => Set([4, 5]), 4 => Set([3, 4]), 5 => Set([2, 5])
    )
    @test maxcard(Set(1:5), Set(1:5), relation3) == 5
    relation4 = Dict(
        1 => Set(1:6), 2 => Set([4]), 3 => Set([4]), 4 => Set([4])
    )
    @test maxcard(Set(1:4), Set(1:6), relation4) == 2
end

end # graph.bipartite


text:
"""
    reference(d::MangalDataset)

Returns the `MangalReference` associated to a `MangalDataset`. This is a
convenience function that returns the `reference` field of the dataset.
"""
function reference(d::MangalDataset)
    return d.reference
end