text:
__precompile__()

module OAuth

using HTTP, MbedTLS, Base64, Random

export
oauth_timestamp,
oauth_nonce,
oauth_sign_hmac_sha1,
oauth_signing_key,
oauth_signature_base_string,
oauth_percent_encode_keys!,
oauth_serialize_url_parameters,
encodeURI!,
oauth_body_hash_file,
oauth_body_hash_data,
oauth_body_hash_encode,
oauth_header,
oauth_request_resource

"""
    oauth_timestamp()

Returns current unix timestamp as String.

# Examples
```julia-repl
julia> oauth_timestamp()
"1512235859"
```
"""
function oauth_timestamp()
    "$(round(Int, time()))"
end

"""
    oauth_nonce(length::Int)

Returns a random string of a given length.

# Examples
```julia-repl
julia> oauth_nonce(10)
"aQb2FVkrYi"
```
"""
function oauth_nonce(length::Int)
    randstring(length)
end

"""
    oauth_sign_hmac_sha1(message::String, signingkey::String)

Takes a message and signing key, converts to a SHA-1 digest, then encodes to base64.

# Examples
```jldoctest
julia> oauth_sign_hmac_sha1("foo", "bar")
"hdFVxV7ShqMAvRzxJN4I2H6RTzo="
```
"""
function oauth_sign_hmac_sha1(message::String, signingkey::String)
    base64encode(digest(MD_SHA1, message, signingkey))
end

"""
    oauth_signing_key(oauth_consumer_secret::String, oauth_token_secret::String)

Returns a signing key based on a consumer secret and token secret.

# Examples
```jldoctest
julia> oauth_signing_key("foo", "bar")
"foo&bar"
```
"""
function oauth_signing_key(oauth_consumer_secret::String, oauth_token_secret::String)
    "$(oauth_consumer_secret)&$(oauth_token_secret)"
end

"""
    oauth_signature_base_string(httpmethod::String, url::String, parameterstring::String)

Returns encoded HTTP method, url and parameters.

# Examples
```jldoctest
julia> oauth_signature_base_string("POST", "https://julialang.org", "foo&bar")
"POST&https%3A%2F%2Fjulialang.org&foo%26bar"
```
"""
function oauth_signature_base_string(httpmethod::String, url::String, parameterstring::String)
    "$(httpmethod)&$(encodeURI(url))&$(encodeURI(parameterstring))"
end

"""
    oauth_percent_encode_keys!(options::Dict)

Returns dict where keys and values are URL-encoded.

# Examples
```jldoctest
julia> oauth_percent_encode_keys!(Dict("key 1" => "value1", "key    2" => "value 2"))
Dict{String,String} with 2 entries:
  "key%20%20%20%202" => "value%202"
  "key%201"          => "value1"
```
"""
function oauth_percent_encode_keys!(options::Dict)
    #options encoded
    originalkeys = collect(keys(options))

    for key in originalkeys
        key_str = string(key)
        encoded_key = encodeURI(key_str)

        options[encoded_key] = encodeURI(options[key_str])
        if encodeURI(key_str) != key
            delete!(options, key_str)
        end
    end

    options
end

@deprecate(
    oauth_percent_encode_keys(options::Dict),
    oauth_percent_encode_keys!(options::Dict)
)

"""
    oauth_serialize_url_parameters(options::Dict)

Returns query string by concatenating dictionary keys/values.

# Examples
```jldoctest
julia> oauth_serialize_url_parameters(Dict("foo" => "bar", "foo 1" => "hello!"))
"foo=bar&foo 1=hello!"
```
"""
oauth_serialize_url_parameters(options::Dict) = join(
    ["$key=$(options[key])" for key in sort!(collect(keys(options)))],
    "&"
)

# See: https://github.com/randyzwitch/OAuth.jl/issues/3
"""
    encodeURI(s)

Convenience function for `HTTP.escapeuri`.

# Examples
```jldoctest
julia> encodeURI("hello, world!")
"hello%2C%20world%21"
```
"""
encodeURI(s) = HTTP.escapeuri(s)

"""
    encodeURI!(dict_of_parameters::Dict)

Mutates dict_of_parameters using `encodeURI` on strings.

# Examples
```jldoctest
julia> encodeURI!(Dict("iv" => 10, "s" => "value!"))
Dict{String,Any} with 2 entries:
  "iv" => 10
  "s"  => "value%21"
```
"""
function encodeURI!(dict_of_parameters::Dict)
    for (k, v) in dict_of_parameters
        if typeof(v) <: String
            dict_of_parameters[k] = encodeURI(v)
        end
    end
    return dict_of_parameters
end

@deprecate(
    encodeURI(dict_of_parameters::Dict),
    encodeURI!(dict_of_parameters::Dict)
)

"""
    oauth_body_hash_file(filename::String)

Returns `oauth_body_hash=` along with base64 encoded SHA-1 from input text file.

# Examples
```jldoctest
julia> oauth_body_hash_file(joinpath(Pkg.dir(), "OAuth/test/auth_body_hash_file.txt"))
"oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="
```
"""
function oauth_body_hash_file(filename::String)
    oauth_body_hash_data(read(open(filename), String))
end

"""
    oauth_body_hash_data(data::String)

Returns `oauth_body_hash=` along with base64 encoded SHA-1 from input.

# Examples
```jldoctest
julia> oauth_body_hash_data("Hello, World!")
"oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="
```
"""
function oauth_body_hash_data(data::String)
    "oauth_body_hash=$(oauth_body_hash_encode(data))"
end

"""
    oauth_body_hash_encode(data::String)

Convenience function for SHA-1 and base64 encoding.

# Examples
```jldoctest
julia> oauth_body_hash_encode("julialang")
"Lsztg2byou89Y8lBoH3G8v3vjbw="
```
"""
function oauth_body_hash_encode(data::String)
        base64encode(digest(MD_SHA1, data))
end

"""
    function oauth_header(httpmethod, baseurl, options, oauth_consumer_key, oauth_consumer_secret, oauth_token, oauth_token_secret; oauth_signature_method = "HMAC-SHA1", oauth_version = "1.0")

Builds OAuth header, defaulting to OAuth 1.0. Function assumes `options` has already
been run through `encodeURI!`.

"""
function oauth_header(httpmethod, baseurl, options, oauth_consumer_key, oauth_consumer_secret, oauth_token, oauth_token_secret; oauth_signature_method = "HMAC-SHA1", oauth_version = "1.0")

    #keys for parameter string
    options["oauth_consumer_key"] = oauth_consumer_key
    options["oauth_nonce"] = oauth_nonce(32)
    options["oauth_signature_method"] = oauth_signature_method
    options["oauth_timestamp"] = oauth_timestamp()
    options["oauth_token"] = oauth_token
    options["oauth_version"] = oauth_version

    #options encoded
    oauth_percent_encode_keys!(options)

    #Create ordered query string
    parameterstring = oauth_serialize_url_parameters(options)

    #Calculate signature_base_string
    signature_base_string = oauth_signature_base_string(uppercase(httpmethod), baseurl, parameterstring)

    #Calculate signing_key
    signing_key = oauth_signing_key(oauth_consumer_secret, oauth_token_secret)

    #Calculate oauth_signature
    oauth_sig = encodeURI(oauth_sign_hmac_sha1(signature_base_string, signing_key))

    return "OAuth oauth_consumer_key=\"$(options["oauth_consumer_key"])\", oauth_nonce=\"$(options["oauth_nonce"])\", oauth_signature=\"$(oauth_sig)\", oauth_signature_method=\"$(options["oauth_signature_method"])\", oauth_timestamp=\"$(options["oauth_timestamp"])\", oauth_token=\"$(options["oauth_token"])\", oauth_version=\"$(options["oauth_version"])\""

end

"""
    oauth_request_resource(endpoint::String, httpmethod::String, options::Dict, oauth_consumer_key::String, oauth_consumer_secret::String, oauth_token::String, oauth_token_secret::String)

Makes `GET` or `POST` call to OAuth API.

"""
function oauth_request_resource(endpoint::String, httpmethod::String, options::Dict, oauth_consumer_key::String, oauth_consumer_secret::String, oauth_token::String, oauth_token_secret::String)
    #Build query string
    query_str = HTTP.escapeuri(options)

    #Build oauth_header
    oauth_header_val = oauth_header(httpmethod, endpoint, options, oauth_consumer_key, oauth_consumer_secret, oauth_token, oauth_token_secret)

    #Make request
    headers = Dict{String,String}(
            "Content-Type" => "application/x-www-form-urlencoded",
            "Authorization" => oauth_header_val,
            "Accept" => "*/*"
        )

    if uppercase(httpmethod) == "POST"
        return HTTP.post(endpoint; body = query_str, headers = headers)
    elseif uppercase(httpmethod) == "GET"
        return HTTP.get("$(endpoint)?$query_str"; headers = headers)
    end
end

end


text:
# # Tilted bottom boundary layer example
#
# This example simulates a two-dimensional (x-z) tilted oceanic bottom boundary layer with a
# constant background alongslope (y-direction) velocity. It demonstrates how to simulate a domain
# tilt by
#
#   * Changing the direction of the buoyancy acceleration
#   * Changing the axis of rotation for Coriolis
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, Plots"
# ```

# ## Load `Oceananigans.jl` and define problem constants
#
# We set a named tuple that defines the necessary simulation parameters:

using Oceananigans
using Oceananigans.Units
using Printf
using CUDA

params = (f₀ = 1e-4, #1/s
          V∞ = 0.1, # m/s
          N²∞ = 1e-5, # 1/s²
          θ_rad = 0.05,
          ν = 5e-4, # m²/s
          z₀ = 0.1, # m (roughness length)
          Lx = 1000, # m
          Lz = 100, # m
          )

# Here `f₀` is the Coriolis frequency, `V∞` in the constant interior `v`-velocity, `N²∞` is the
# interior stratification, `θ_rad` is the bottom slope in radians, `ν` is the eddy viscosity, and
# `z₀` is the roughness length (needed for the drag at the bottom). `Lx` and `Lz` are the domain
# length and height.
#
#
# ## Creating the grid
#
# Since we anticipate most of the energy in the flow to be located close to the bottom, we
# can save computational resources by creating a grid that has finer spacings near the bottom
# which gets progressively coarser near the top. Such a grid can be achieved with the following
# equations:

Nx = 128
Nz = 32

refinement = 1.8 # controls spacing near surface (higher means finer spaced)
stretching = 10   # controls rate of stretching at bottom

h₁(k) = ((-k+Nz)+1) / Nz

## Linear near-surface generator
ζ₁(k) = 1 + (h₁(k) - 1) / refinement

## Bottom-intensified stretching function
Σ₁(k) = (1 - exp(-stretching * h₁(k))) / (1 - exp(-stretching))

## Generating function
z_faces(k) = -params.Lz * (ζ₁(k) * Σ₁(k) - 1)

grid = RectilinearGrid(topology=(Periodic, Flat, Bounded),
                       size=(Nx, Nz),
                       x=(0, params.Lx), z=z_faces,
                       halo=(3,3),
                       )

# We plot vertical spacing versus depth to inspect the prescribed grid stretching. Note that
# the spacing near the bottom is relatively uniform with height, which is a desired property from
# a numerical perspective.
#

using Plots

plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
     marker = :circle,
     ylabel = "Depth (m)",
     xlabel = "Vertical spacing (m)",
     legend = nothing)

# ## Setting up a buoyancy model in a tilted domain
#
# We set-up our domain in a way that the coordinates align with the tilted bottom. That means that,
# from the coordinate's point of view, gravity is tilted by `θ_rad` radians, which we can simulate
# by passing the `gravity_unit_vector` parameter to `Buoyancy()`:

ĝ = [sin(params.θ_rad), 0, cos(params.θ_rad)]
buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=ĝ)

# ## Background fields
#
# Since we are simulating a period domain that is tilted, we need to set-up the background
# stratification as a `BackgroundField` and solve for the buoyancy perturbation. This avoids mean
# horizontal buoyant gradients, which cannot exist in a horizontally-periodic domain:

b∞(x, y, z, t, p) = p.N²∞ * (x * sin(p.θ_rad) + z * cos(p.θ_rad))
B_field = BackgroundField(b∞, parameters=(; params.N²∞, params.θ_rad))

# We also set the constant interior velocity `V∞` as a background field in order to avoid inertial
# oscillations due to an unbalanced resolved flow:
#

V_bg(x, y, z, t, p) = p.V∞
V_field = BackgroundField(V_bg, parameters=(; V∞=params.V∞))

# ## Bottom drag
#
# We set-up a bottom drag that follows Monin-Obukhov theory. Note that we need to include`V∞`
# in the velocity, since we will model it as a background field.
#

const κ = 0.4 # von Karman constant
z₁ = CUDA.@allowscalar znodes(Center, grid)[1] # Closest grid center to the bottom
cᴰ = (κ / log(z₁/params.z₀))^2 # Drag coefficient

@inline drag_u(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * u
@inline drag_v(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * (v + p.V∞)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(cᴰ=cᴰ, V∞=params.V∞))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(cᴰ=cᴰ, V∞=params.V∞))

u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

# ## Create the `NonhydrostaticModel`
#
# Since the coordinate system is aligned with the tilted domain, we can use a traditional f-plane
# approximation (only keep the rotation component that's aligned with gravity) by also tilting the
# rotation from the coordinate system's perpective by `θ_rad` radians using
# `ConstantCartesianCoriolis` and passing the `rotation_axis` argument:
#

coriolis = ConstantCartesianCoriolis(f=params.f₀, rotation_axis=ĝ)

#
# We are now ready to create the model. We create a `NonhydrostaticModel` with an
# `UpwindBiasedFifthOrder` advection scheme and a `RungeKutta3` timestepper.  For simplicity, we use
# a constant-diffusivity turbulence closure:

model = NonhydrostaticModel(grid = grid, timestepper = :RungeKutta3,
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = buoyancy,
                            coriolis = coriolis,
                            tracers = :b,
                            closure = ScalarDiffusivity(ν=params.ν, κ=params.ν),
                            boundary_conditions = (u=u_bcs, v=v_bcs,),
                            background_fields = (b=B_field, v=V_field,),
                           )

# ## Create and run a simulation
#
# We are now ready to create the simulation. We begin by setting the initial time step
# conservatively, based on the smallest grid size of our domain and set-up a

using Oceananigans.Grids: min_Δz
simulation = Simulation(model,
                        Δt=0.5*min_Δz(grid)/params.V∞,
                        stop_time=2days)

# We now add callbacks to adjust the time-step and to display the simulation progress:

wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))

start_time = time_ns() # so we can print the total elapsed wall time
progress_message(sim) =
    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
            iteration(sim), prettytime(time(sim)),
            prettytime(sim.Δt), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

# ## Add outputs to the simulation
#
# We add outputs to our model using the `NetCDFOutputWriter`

u, v, w = model.velocities

b_tot = Field(model.tracers.b + model.background_fields.tracers.b)
v_tot = Field(v + model.background_fields.velocities.v)
ω_y = Field(∂z(u)-∂x(w))

fields = (; u, v_tot, w, b_tot, ω_y)

simulation.output_writers[:fields] = NetCDFOutputWriter(model, fields,
                                                        filename = joinpath(@__DIR__, "out.tilted_bbl.nc"),
                                                        schedule = TimeInterval(20minutes),
                                                        overwrite_existing = true)

# Now we just run it!

run!(simulation)

# ## Visualize the results
#
# First we load the required package to load NetCDF output files and define the coordinates for
# plotting using existing objects:
#

using NCDatasets
xᶠ, zᶠ = xnodes(ω_y), znodes(ω_y);
xᶜ, zᶜ = xnodes(v_tot), znodes(v_tot) ;

# Define keyword arguments for plotting function
#
kwargs = (xlabel = "x",
          zlabel = "z",
          fill = true,
          levels = 30,
          linewidth = 0,
          color = :balance,
          colorbar = true,
          xlim = (0, params.Lx),
          zlim = (0, params.Lz)
         )

# Read in the output_writer for the two-dimensional fields and then create an animation showing the
# vorticity.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

anim = @animate for (iter, t) in enumerate(ds["time"])
    ω_y = ds["ω_y"][:,1,:, iter]
    v_tot = ds["v_tot"][:,1,:,iter]

    ω_max = maximum(abs, ω_y)
    v_max = maximum(abs, v_tot)

    plot_ω = contour(xᶠ, zᶠ, ω_y',
                     clim = (-0.015, +0.015),
                     title = @sprintf("y-vorticity, ω_y, at t = %.1f", t); kwargs...)

    plot_v = contour(xᶜ, zᶜ, v_tot',
                     clim = (-params.V∞, +params.V∞),
                     title = @sprintf("Total along-slope velocity, at t = %.1f", t); kwargs...)

    plot(plot_ω, plot_v, layout = (2, 1), size = (800, 440))
end

close(ds)

mp4(anim, "tilted_bottom_boundary_layer.mp4", fps=15)

ext: jl


text:
include("binomial.jl")

function compute()
	const n=40
	const m=20
	ans=binomial(40,20)

	return ans
end

println(compute())


text:
# Testing reference values and precisions
# Each test block of varr and parr should be followed by an append to refVals, refPrecs arrays.
# e.g.
#   refVals=[]
#   refPrecs=[]
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   etc.....
#
#   Now for real!
#
#! format: off

    refVals=[]
    refPrecs=[]

# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/cnh/projects/cm/experiments/OceanBoxGCM/ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -7.88803234115370288659e-02,  6.54328011157667283060e-02,  2.01451525264660777012e-03,  1.24588716783015467093e-02 ],
 [       "Q",   "u[2]", -8.71605274310934263760e-02,  1.47505267452495314462e-01,  5.63158201082454307890e-03,  1.28959719954838437916e-02 ],
 [       "Q",       :η, -4.73491408441974459542e-01,  4.02693687285993751068e-01, -6.49059970676089264836e-05,  2.21689928090662985438e-01 ],
 [       "Q",       :θ,  4.24292935192304451839e-04,  9.24539353455580759089e+00,  2.49938206627401893201e+00,  2.17986626392490689952e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  1.99999999999999976717e+06,  1.15573163901915703900e+06 ],
 [   "s_aux",       :w, -2.22086406767055965930e-04,  2.00575090959195941574e-04,  2.53168866096191178104e-07,  1.66257132341925427171e-05 ],
 [   "s_aux",    :pkin, -9.00869877619915215838e-01,  0.00000000000000000000e+00, -3.33171488779369640021e-01,  2.54740287894525019308e-01 ],
 [   "s_aux",     :wz0, -2.96608015794916532427e-05,  3.66759042312538147113e-05,  3.78116102589123325624e-10,  1.07073256826409470244e-05 ],
 [ "s_gflux", "ν∇u[1]", -1.98765276856145908790e-02,  1.92976912753076192131e-02,  9.31995586073550077850e-09,  2.84016926475062683050e-04 ],
 [ "s_gflux", "ν∇u[2]", -8.60188027210821715141e-03,  2.03898608650850965418e-02,  1.44279862042391183595e-05,  7.01831343135553482361e-04 ],
 [ "s_gflux", "ν∇u[3]", -2.36999031942582461820e-05,  3.39171933774427215489e-05,  1.79310417841609871239e-08,  3.03350466717279338539e-06 ],
 [ "s_gflux", "ν∇u[4]", -3.27811105688537093861e-02,  3.33776923250933485465e-02, -8.41954909635648634245e-06,  9.93067541524441646170e-04 ],
 [ "s_gflux", "ν∇u[5]", -3.19705003898926379269e-02,  3.94072309121676778859e-03, -2.09915474589233714674e-06,  4.56596518968874199012e-04 ],
 [ "s_gflux", "ν∇u[6]", -3.84552825296212839770e-05,  3.30579003830798668533e-05, -6.85753430789401392363e-08,  3.41712910764070301712e-06 ],
 [ "s_gflux", "κ∇θ[1]", -1.90663030661463762938e-03,  1.59419367775516821958e-03, -2.07346812692515548133e-05,  1.38241222633447261199e-04 ],
 [ "s_gflux", "κ∇θ[2]", -1.17846832845764997420e-03,  3.03333905110560143534e-03,  9.95380871168155762288e-04,  8.10649417649871592909e-04 ],
 [ "s_gflux", "κ∇θ[3]", -1.91676759182478058212e-05,  6.14879473708830387359e-05, -4.12675761035662034755e-07,  1.91038619046619808554e-06 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,     8,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,     8,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

# SC ========== Test number 2 reference values and precision match template. =======
# SC ========== /home/cnh/projects/cm/experiments/OceanBoxGCM/ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -3.08415268094245687003e-01,  3.31683879748177812274e-01,  2.34355645025035056947e-03,  1.38278493625595393784e-02 ],
 [       "Q",   "u[2]", -9.30330187375882494694e-02,  1.45773331490709479041e-01,  6.38739141648009701308e-03,  1.38107523198905141754e-02 ],
 [       "Q",       :η, -4.60215148246349181615e-01,  3.97155090993862369686e-01, -1.05814653709044253051e-04,  2.14356824554893188317e-01 ],
 [       "Q",       :θ,  3.90305009874997816885e-04,  9.28488770009523101123e+00,  2.49938688339290848717e+00,  2.17986650147491722862e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  1.99999999999999976717e+06,  1.15573163901915703900e+06 ],
 [   "s_aux",       :w, -2.65432487944419423249e-04,  2.51335486222488370053e-04,  3.34985017281366756292e-07,  1.98256782258218657716e-05 ],
 [   "s_aux",    :pkin, -9.01203660410079376852e-01,  0.00000000000000000000e+00, -3.33173652152760846334e-01,  2.54737286187841360796e-01 ],
 [   "s_aux",     :wz0, -3.16776291712146345464e-05,  3.89428424525183031852e-05,  4.05688221636862481177e-10,  1.10561449663397378925e-05 ],
 [ "s_gflux", "ν∇u[1]", -6.91411158656193985950e-02,  6.71894067927724242262e-02,  1.93151282317531301526e-08,  3.67871498890457999768e-04 ],
 [ "s_gflux", "ν∇u[2]", -6.41072589623675609394e-02,  5.93158959821227672188e-03, -1.72612717858599068401e-05,  1.09370391833557892812e-03 ],
 [ "s_gflux", "ν∇u[3]", -7.40852152420143496942e-05,  7.70630509235797829829e-05,  1.94096112473072846645e-08,  3.35721658566900582726e-06 ],
 [ "s_gflux", "ν∇u[4]", -2.77618725100101104708e-02,  2.76053769492804849939e-02,  8.37000168929151316390e-06,  1.04649435866464054871e-03 ],
 [ "s_gflux", "ν∇u[5]", -3.57980327012049140412e-02,  4.20082014423407961584e-03, -2.26270579934147913334e-06,  5.16063361553070450508e-04 ],
 [ "s_gflux", "ν∇u[6]", -3.87729590851969090591e-05,  3.38326906830031177748e-05, -7.64232106453937395084e-08,  3.42855062975519724920e-06 ],
 [ "s_gflux", "κ∇θ[1]", -2.63114521554794764460e-03,  3.13139570197398927828e-03, -2.58270397697685688597e-05,  1.80489106756528977114e-04 ],
 [ "s_gflux", "κ∇θ[2]", -4.29150657020392189284e-03,  3.20181945576268351533e-03,  9.94194034856019440149e-04,  8.13274467049271117605e-04 ],
 [ "s_gflux", "κ∇θ[3]", -3.08119171379054497821e-05,  9.12508944783092335127e-05, -4.12575332587978521639e-07,  1.92043362474637987870e-06 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,     8,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,     8,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
module Optimise

export train!,
  SGD, ADAM, ADAMW, AdaMax, Momentum, Nesterov,
  RMSProp, ADAGrad, ADADelta, AMSGrad, NADAM

struct Param{T}
  x::T
  Δ::T
end

Base.convert(::Type{Param}, x::AbstractArray) = Param(x, zeros(x))

include("optimisers.jl")
include("interface.jl")
include("train.jl")

using Flux.Tracker: TrackedArray

Base.convert(::Type{Param}, x::TrackedArray) = Param(x.data, x.grad)

end


text:
#= As far as I know, Github Actions uses Intel CPUs.
So it is faster to use MKL than OpenBLAS.
It is recommended to load MKL before ANY other package.=#
import MKL
using SafeTestsets: @safetestset

@time @safetestset "wnorm" begin include("wnorm.jl") end
@time @safetestset "rbm" begin include("rbm.jl") end


text:
"""
    getreference(blocks; axis=["X","Y","Z"])

Extract the mid-surface points that split the informed block model in two.
Returns a coordinate matrix with the reference points for unfolding.

## Parameters:

* `blocks` - coordinate matrix of the regular blocks centroids.
* `axis`   - axis or list of axis that will be scanned to extract reference
  points. Defaults to ["X","Y","Z"].
"""
function getreference(blocks::AbstractMatrix; axis=["X","Y","Z"])
	typeof(axis) in (String, Symbol) && (axis = [axis])
	axis[1] isa Symbol && (axis = String.(axis))
	@assert axis ⊆ ["X","Y","Z"] "invalid axis: $axis"
	coord_ids = Dict("X"=>1, "Y"=>2, "Z"=>3)

	# get cell sizes
	cells = zeros(Float64,3)
	for c in 1:3
		coords   = sort!(unique(blocks[c,:]))
		cells[c] = minimum([coords[x]-coords[x-1] for x in 2:length(coords)])
	end

	# scan sections and save skeleton points
	ref = Float64[]
	for c in axis
		ax = coord_ids[c]
	    axis_coords = sort!(unique(blocks[ax,:]))
		sec = setdiff([1,2,3],ax)

	    for s in axis_coords
	        section = view(blocks,:,view(blocks,ax,:) .== s)

	        min_i = minimum(section[sec[1],:])
	        min_j = minimum(section[sec[2],:])
	        i = section[sec[1],:] .- min_i
	        j = section[sec[2],:] .- min_j

			i = Int.(i/cells[sec[1]]) .+ 1
			j = Int.(j/cells[sec[2]]) .+ 1

	        img = zeros(Bool,maximum(i)+1,maximum(j)+1)

	        for x in 1:length(i)
	            img[i[x],j[x]]=true
	        end

	        img = thinning(img)

	        for x in 1:length(i)
	            if img[i[x],j[x]]==true
	                io = (i[x]-1)*cells[sec[1]]+min_i
	                jo = (j[x]-1)*cells[sec[2]]+min_j
					order = sortperm([ax,sec[1],sec[2]])
					append!(ref, [s, io, jo][order])
	            end
	        end
	    end
	end

	reshape(ref, 3, Int(length(ref)/3))
end

# Average duplicate points
function remove_duplicates(coords::AbstractMatrix; tol=0.01)
	groups = [union(x.core_indices,x.boundary_indices) for x in dbscan(coords, tol)]
	outcoords = zeros(Float64,3,length(groups))
	for g in 1:length(groups)
		idx = groups[g]
		if length(idx)==1
			outcoords[:,g] .= dropdims(view(coords,:,idx),dims=2)
		else
			outcoords[:,g] .= mean!(view(outcoords,:,g),view(coords,:,idx))
		end
	end
	outcoords
end

function get_resolution(ref_pts::AbstractMatrix)
	idxs, dists = get_neighbors(ref_pts, :knn, 2)
	closest_pt = [sum(x) for x in dists]
	quantile(closest_pt,0.75)
end


text:
"""
Computes first and second moments of `` \\exp[Jx_{i}x_{j} + h_{i}x_{i} +h_{j}x_{j} -\\frac12 x^t \\cdot S \\cdot x + x^t \\cdot y] ``

Input:
* `F::PairIsing <:Factor`: Factor node ``\\psi_{ij}``
* `y::Vector{Float64}`: Array of cavity fields ``\\in \\mathbb{R}^{ 2 }``
* `S::Vector{Float64}`: Cavity coupling matrix ``\\in \\mathbb{R}^{  2 \\times 2 }``
"""
function moments(F::IsingPair, y::Vector{Float64}, S::Array{Float64,2})
    mi = tanh(F.hi + y[1] + atanh(tanh(F.J - S[1,2]) * tanh(F.hj + y[2])))
    mj = tanh(F.hj + y[2] + atanh(tanh(F.J - S[1,2]) * tanh(F.hi + y[1])))
    cij = tanh(F.J - S[1,2] + atanh(tanh(F.hj + y[2]) * tanh(F.hi + y[1]))) - mi*mj
    return [mi;mj],[1-mi^2 cij; cij 1-mj^2]
end

"""
Computes first and second moments of `` \\exp[-\\psi.E(x) -\\frac12 x^t \\cdot S \\cdot x + x^t \\cdot y] ``

Input:
* `F::EnergyFun <:Factor`: Factor node ``\\psi_{a}``
* `y::Vector{Float64}`: Array of cavity fields ``\\in \\mathbb{R}^{ |\\partial a| }``
* `S::Vector{Float64}`: Cavity coupling matrix ``\\in \\mathbb{R}^{  |\\partial a| \\times |\\partial a| }``

Output: first and second moments under tilted "a" distribution -> `av, cov`
where `av` is a vector of size ``|\\partial a|``, `cov`  is a ``|\\partial a| \\times |\\partial a|`` matrix }
"""
function moments(F::EnergyFun, y::Vector{Float64}, S::Array{Float64,2})
    n = length(F.idx)
    N = 1<<n
    av = zeros(n)
    cov = zeros(n,n)
    z = 0
    X = [ [2*((i >> k)%2)-1 for k in 0:(n-1)] for i=0:(N-1) ]
    Y = [ -F.E(X[i]) + X[i]'*(y - 0.5*S*X[i]) for i=1:N ]
    m = maximum(Y)
    @inbounds for i=1:N
        px = exp(Y[i]-m)
        @assert(px>=0)
        z += px
        av .+= X[i] * px
        cov .+= (X[i]*X[i]') * px
    end
    av ./= z
    cov ./= z
    cov .-= av*av'
    return av, cov
end

pearsonmap(x) = sqrt(1-x)/(atanh(sqrt(1-x))*x)

"""
Update the parameters of the gaussian factor "a" according to closure.

Input:

* `F::EnergyFun <:Factor`: Factor node ``\\psi_{a}``
* `y::Vector{Float64}`: Array of cavity fields of size ``|\\partial a|``
* `S::Vector{Float64}`: Cavity coupling matrix of size ``|\\partial a| \\times |\\partial a|``
* `μta::Vector{Float64}`: tilted magnetizations, size ``|\\partial a|`` -> gets updated as ref.
* `Σta::Matrix{Float64}`: Correlation matrix , size ``|\\partial a| \\times |\\partial a|`` -> gets updated as ref.
* `closure::Symbol`: closure protocol
* `η::Float64`: interpolation parameter (valid on :DC closure) -> ``η=0`` gives BP fixed points
* `epsclamp::Float64` clamp tilted moments in [-1+epsclamp, 1-epsclamp]
* `epsconv::Float64`: precision convergence
* `epsmom::Float64`: error over moments at iteration t

Output:

* `ynew::Vector{Float64}`: fields of gauss. factor a @iter t+1, size ``|\\partial a|``
* `Snew::Vector{Float64}`: coupling matrix of gauss. factor a @iter t+1, ``|\\partial a| \\times |\\partial a|``
* `epsmom::Float64`: error over moments at iteration t+1

"""
function setclosure!(F::Factor, y::Vector{Float64},
                        S::Array{Float64,2},
                        µta::Vector{Float64},
                        Σta::Array{Float64,2},
                        closure::Symbol,
                        η::Float64,
                        λ::Float64,
                        epsclamp::Float64,
                        epsmom::Float64)

  av,cov = moments(F, y, S)
  epsmom = max(epsmom , update!(μta,av,0.0), update!(Σta,cov,0.0))
    if closure == :DC
        var = sqrt.(pearsonmap.(clamp.(diag(cov), epsclamp, 1-epsclamp)))
        cov .*= var * var'
        cov = η*cov+(1-η)*Diagonal(cov)
        icov = inv(cov + λ * I)
        Snew = icov-S
        ynew = icov * av - y
    elseif closure == :DC2
        cov = η*cov + (1-η) * Diagonal(cov)
        icov = inv(cov + λ * I)
        means = diag(cov).*atanh.(clamp.(av, -1+epsclamp,1-epsclamp))
        Snew = icov - S
        ynew = icov * means - y
    elseif closure == :EP
        av, cov = moments(F,y,S)
        icov = inv(cov + λ*I)
        Snew = icov - S;
        ynew = icov * av - y
    else
        throw("Select a valid closure")
    end
    return ynew,Snew,epsmom
end

max_stars_repo_name: giovact/DensityConsistency

max_issues_repo_name: giovact/DensityConsistency

max_forks_repo_name: giovact/DensityConsistency


text:
# Try getting token
token = getRepositoryToken(mPulseAPITenant, mPulseAPIToken)
@test !isempty(token)

# Get all domains
domains = getRepositoryDomain(token)
@test length(domains) == 2

domainNames = map(d -> d["name"], domains)
@test "mPulse Demo" ∈ domainNames
@test "mPulseAPI Test" ∈ domainNames

# Get a specific domain
appKey = filter(d -> d["name"] == "mPulseAPI Test", domains)[1]["attributes"]["appKey"]
@test !isempty(appKey)

domain = getRepositoryDomain(token, appKey=appKey)
@test !isempty(domain)

# Check domain parameters
@test domain["name"] == "mPulseAPI Test"
@test domain["attributes"]["appKey"] == appKey
@test domain["resource_timing"]

# Check custom metrics
@test isa(domain["custom_metrics"], Dict)
@test length(domain["custom_metrics"]) == 3

@test haskey(domain["custom_metrics"], "Conversion")
@test haskey(domain["custom_metrics"], "OrderTotal")
@test haskey(domain["custom_metrics"], "Revenue GBP")

@test domain["custom_metrics"]["Conversion"]["index"] == 0
@test domain["custom_metrics"]["Conversion"]["fieldname"] == "custommetric0"
@test domain["custom_metrics"]["Conversion"]["dataType"]["type"] == "Percentage"

@test domain["custom_metrics"]["Revenue GBP"]["index"] == 1
@test domain["custom_metrics"]["Revenue GBP"]["fieldname"] == "custommetric1"
@test domain["custom_metrics"]["Revenue GBP"]["dataType"]["type"] == "Currency"
@test domain["custom_metrics"]["Revenue GBP"]["dataType"]["currencyCode"] == "GBP"
# Bug in mPulseAPI means that `currencySymbol` will sometimes disappear
@test !haskey(domain["custom_metrics"]["Revenue GBP"]["dataType"], "currencySymbol") || domain["custom_metrics"]["Revenue GBP"]["dataType"]["currencySymbol"] == "£"
@test domain["custom_metrics"]["Revenue GBP"]["dataType"]["decimalPlaces"] == "2"

# CustomMetric2 is inactive

@test domain["custom_metrics"]["OrderTotal"]["index"] == 3
@test domain["custom_metrics"]["OrderTotal"]["fieldname"] == "custommetric3"
@test domain["custom_metrics"]["OrderTotal"]["dataType"]["type"] == "Currency"
@test domain["custom_metrics"]["OrderTotal"]["dataType"]["currencyCode"] == "USD"
@test domain["custom_metrics"]["OrderTotal"]["dataType"]["decimalPlaces"] == "2"

# Check custom timers
@test isa(domain["custom_timers"], Dict)
@test length(domain["custom_timers"]) == 1

@test haskey(domain["custom_timers"], "ResourceTimer")

@test domain["custom_timers"]["ResourceTimer"]["index"] == 0
@test domain["custom_timers"]["ResourceTimer"]["fieldname"] == "customtimer0"
@test domain["custom_timers"]["ResourceTimer"]["mpulseapiname"] == "CustomTimer0"

# Check tenant
tenant = getRepositoryTenant(token, name=mPulseAPITenant)
@test !isempty(tenant)
@test tenant["name"] == mPulseAPITenant

if !isempty(mPulseAPIAlert)
    # Check alert
    alert = getRepositoryAlert(token, alertName=mPulseAPIAlert)
    @test !isempty(alert)
    @test alert["name"] == "mPulseAPI Test Alert"
    @test alert["tenantID"] == 236904
    @test alert["tenantID"] == tenant["id"]
end

# Now check all exceptions
@test_throws ArgumentError getRepositoryToken("", "")

@test_throws mPulseAPIAuthException getRepositoryToken(mPulseAPITenant, "some-invalid-token")

try
    getRepositoryDomain("some-invalid-token", appKey=appKey * "-" * base(16, round(Int, rand()*100000), 5))
catch ex
    if isa(ex, mPulseAPIBugException)
        warn("Expected mPulseAPIAuthException, got mPulseAPIBugException", ex)
    else
        @test isa(ex, mPulseAPIAuthException) || isa(ex, mPulseAPIBugException)
    end
end

#### Dynamic Alerting ###
# Check alert
if !isempty(DA_mPulseAPIAlert)
    DAalert = getRepositoryAlert(token, alertName=DA_mPulseAPIAlert)
    @test !isempty(alert)
    @test DAalert["name"] == "mPulseAPI Dynamic Test Alert"
    @test DAalert["id"] == 2251091
    @test DAalert["tenantID"] == 236904
    @test DAalert["tenantID"] == tenant["id"]
    @test DAalert["attributes"]["dynamic"] == true
    @test DAalert["attributes"]["statisticalModelID"] == 415
    @test DAalert["attributes"]["state"] ∈ ["AutoCleared", "Updated", "Active", "ModelNotReady"]

    # Update alert via post request
    postRepositoryAlert(token, alertID = DAalert["id"], attributes = Dict("version" => 2))

    # Check statistical model
    statModel = getRepositoryStatModel(token, statModelID = DAalert["attributes"]["statisticalModelID"])
    @test !isempty(statModel)
    @test statModel["id"] == 415
    @test statModel["parentID"] == DAalert["id"]
    @test statModel["parentID"] == 2251091
    @test statModel["tenantID"] == tenant["id"]
    @test statModel["tenantID"] == 236904
    @test statModel["attributes"]["type"] == 1
    @test statModel["attributes"]["version"] == 1.0

    # Update statistical model via post request
    postRepositoryStatModel(token, statModelID = statModel["id"], attributes = Dict("type" => 1))
end


text:
module _DataStructure

import Brainstorm._DataStructure: create, combine

using Brainstorm._DataStructure
using Base.Test

include("iterator.jl")
include("vector.jl")
include("matrix.jl")
include("intset.jl")

function test_all()
    println("")
    test_iterator_all()
    test_vector_all()
    test_matrix_all()
    test_intset_all()
end

end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# This file contains configuration-related functions and types

# Structure to keep Gaston's configuration
mutable struct GastonConfig
    # default CurveConf values
    legend::String
    plotstyle::String
    color::String
    marker::String
    linewidth::Real
    linestyle::String
    pointsize::Real
    # default AxesConf values
    title::String
    xlabel::String
    ylabel::String
    zlabel::String
    fill::String
    grid::String
    box::String
    axis::String
    xrange::String
    yrange::String
    zrange::String
    palette::String
    # default terminal type
    terminal::String
    # for terminals that support filenames
    outputfile::String
    # for printing to file
    print_color::String
    print_fontface::String
    print_fontsize::Real
    print_fontscale::Real
    print_linewidth::Real
    print_size::String     # user-configured print size
    print_size_in::String  # some terminals specify size in inches
    print_size_pix::String # ... and some in pixels
    # prefix for temp data files
    tmpprefix::String
end

function GastonConfig()
    gc = GastonConfig(
	# CurveConf
	"","lines","blue","",1,"",0.5,
	# AxesConf
	"","","","","empty","off","inside vertical right top","",
	"[*:*]","[*:*]","[*:*]","",
	# terminal
	"qt",
	# output file name
	"",
    # print parameters
    "color", "Sans", 12, 0.5, 1, "", "5in,3in", "800,600",
    # tmp file prefix
    randstring(8)
    )
    # set terminal if running inside a Jupyter notebook
    if isjupyter
        gc.terminal = "ijulia"
        gc.print_fontsize = 20
    end
    return gc
end

# Set any of Gaston's configuration variables
# This function assumes that gaston_config exists
function set(;legend         = gaston_config.legend,
             plotstyle       = gaston_config.plotstyle,
             color           = gaston_config.color,
             marker          = gaston_config.marker,
             linewidth       = gaston_config.linewidth,
             linestyle       = gaston_config.linestyle,
             pointsize       = gaston_config.pointsize,
             title           = gaston_config.title,
             xlabel          = gaston_config.xlabel,
             ylabel          = gaston_config.ylabel,
             zlabel          = gaston_config.zlabel,
             fill            = gaston_config.fill,
             grid            = gaston_config.grid,
             box             = gaston_config.box,
             axis            = gaston_config.axis,
             xrange          = gaston_config.xrange,
             yrange          = gaston_config.yrange,
             zrange          = gaston_config.zrange,
             palette         = gaston_config.palette,
             terminal        = gaston_config.terminal,
             outputfile      = gaston_config.outputfile,
             print_color     = gaston_config.print_color,
             print_fontface  = gaston_config.print_fontface,
             print_fontscale = gaston_config.print_fontscale,
             print_linewidth = gaston_config.print_linewidth,
             print_size      = gaston_config.print_size,
             print_size_in   = gaston_config.print_size_in,
             print_size_pix  = gaston_config.print_size_pix)
    # Validate paramaters
    @assert valid_label(legend) "Legend must be a string."
    @assert valid_plotstyle(plotstyle) "Plotstyle $(plotstyle) not supported."
    @assert valid_label(color) "Color must be a string."
    @assert valid_marker(marker) "Marker $(marker) not supported."
    @assert valid_number(linewidth) "Invalid linewdith."
    @assert valid_number(pointsize) "Invalid pointsize."
    @assert valid_label(title) "Title must be a string."
    @assert valid_label(xlabel) "xlabel must be a string."
    @assert valid_label(ylabel) "ylabel must be a string."
    @assert valid_label(zlabel) "zlabel must be a string."
    @assert valid_fill(fill) "Fill style $(fill) not supported."
    @assert valid_grid(grid) "Grid style $(grid) not supported."
    @assert valid_label(box) "box must be a string."
    @assert valid_axis(axis) "Axis $(axis) not supported."
    @assert valid_range(xrange) "Range $(xrange) not supported."
    @assert valid_range(yrange) "Range $(yrange) not supported."
    @assert valid_range(zrange) "Range $(yrange) not supported."
    @assert valid_label(palette) "Invalid palette."
    @assert valid_terminal(terminal) "Terminal type $(terminal) not supported."
    @assert valid_label(outputfile) "Outputfile must be a string."
    @assert valid_label(print_color) "print_color must be a string."
    @assert valid_label(print_fontface) "print_fontface must be a string."
    @assert valid_number(print_fontscale) "Invalid value of print_fontscale"
    @assert valid_number(print_linewidth) "Invalid value of print_linewidth"
    @assert valid_label(print_size) "print_size must be a string."
    @assert valid_label(print_size_in) "print_size_in must be a string."
    @assert valid_label(print_size_pix) "print_size_pix must be a string."
    @assert valid_linestyle(linestyle) string("Line style pattern accepts: ",
            "space, dash, underscore and dot")

    gaston_config.legend            = legend
    gaston_config.plotstyle         = plotstyle
    gaston_config.color             = color
    gaston_config.marker            = marker
    gaston_config.linewidth         = linewidth
    gaston_config.linestyle         = linestyle
    gaston_config.pointsize         = pointsize
    gaston_config.title             = title
    gaston_config.xlabel            = xlabel
    gaston_config.ylabel            = ylabel
    gaston_config.zlabel            = zlabel
    gaston_config.fill              = fill
    gaston_config.grid              = grid
    gaston_config.box               = box
    gaston_config.axis              = axis
    gaston_config.xrange            = xrange
    gaston_config.yrange            = yrange
    gaston_config.zrange            = zrange
    gaston_config.palette           = palette
    gaston_config.outputfile        = outputfile
    gaston_config.print_color       = print_color
    gaston_config.print_fontface    = print_fontface
    gaston_config.print_fontscale   = print_fontscale
    gaston_config.print_linewidth   = print_linewidth
    gaston_config.print_size        = print_size
    gaston_config.print_size_in     = print_size_in
    gaston_config.print_size_pix    = print_size_pix
    # don't change terminal inside jupyter
    if terminal != "ijulia" && isjupyter
        @warn("Terminal cannot be changed in a Jupyter notebook.")
        gaston_config.terminal = "ijulia"
    else
        gaston_config.terminal = terminal
    end
    return nothing
end

### Encode terminal capabilities
# supports multiple windows
const term_window = ["qt", "wxt", "x11", "aqua"]
# outputs text
const term_text = ["dumb", "null", "sixelgd", "ijulia"]
# outputs to a file
const term_file = ["svg", "gif", "png", "pdf", "eps"]
# size and units
const term_sup_size = ["qt", "wxt", "x11", "sixel", "svg", "gif", "png",
                        "pdf", "eps"]
const term_size_in = ["pdf", "eps"]
const term_size_pix = ["qt", "wxt", "x11", "sixelgd", "ijulia",
                        "svg", "gif", "png"]
# font
const term_sup_font = ["qt", "wxt", "x11", "aqua", "sixel", "svg", "gif",
                       "png", "pdf", "eps"]
# font scaling
const term_sup_fontscale = ["wxt", "sixel", "gif", "png", "pdf", "eps"]
# linewidth
const term_sup_lw = ["qt", "wxt", "x11", "aqua", "sixel", "svg", "gif",
                     "png", "pdf", "eps"]

# List of valid configuration values
const supported_terminals = ["qt", "wxt", "x11", "aqua", "dumb", "null",
                             "sixelgd", "ijulia", "svg", "gif", "png", "pdf",
                             "eps"]
const supported_2Dplotstyles = ["lines", "linespoints", "points",
        "impulses", "boxes", "errorlines", "errorbars", "dots", "steps",
        "fsteps", "fillsteps", "financebars"]
const supported_3Dplotstyles = ["lines", "linespoints", "points",
        "impulses", "pm3d", "image", "rgbimage", "dots"]
const supported_plotstyles = vcat(supported_2Dplotstyles, supported_3Dplotstyles)
const supported_axis = ["", "normal", "semilogx", "semilogy", "loglog"]
const supported_markers = ["", "+", "x", "*", "esquare", "fsquare",
        "ecircle", "fcircle", "etrianup", "ftrianup", "etriandn",
        "ftriandn", "edmd", "fdmd"]
const supported_fillstyles = ["","empty","solid","pattern"]
const supported_grids = ["", "on", "off"]

#
# Validation functions
#

valid_terminal(s) = s ∈ supported_terminals
valid_plotstyle(s) = s ∈ supported_plotstyles
valid_2Dplotstyle(s) = s ∈ supported_2Dplotstyles
valid_3Dplotstyle(s) = s ∈ supported_3Dplotstyles
valid_marker(s) = s ∈ supported_markers
valid_number(s) = isa(s, Real) && s > 0
valid_label(s) = isa(s, String)
valid_fill(s) = s ∈ supported_fillstyles
valid_grid(s) = s ∈ supported_grids
valid_axis(s) = s ∈ supported_axis

function valid_linestyle(s)
    s == "" && return true # allow empty string
    c = collect(s)
    # make sure only allowed characters are passed
    issubset(c, Set([' ', '-', '_', '.'])) || return false
    # but do not allow spaces only
    unique(c) != [' ']
end

# Validate that a given range follows gnuplot's syntax.
function valid_range(s::AbstractString)
    # floating point, starting with a dot
    f1 = "[-+]?\\.\\d+([eE][-+]?\\d+)?"
    # floating point, starting with a digit
    f2 = "[-+]?\\d+(\\.\\d*)?([eE][-+]?\\d+)?"
    # floating point
    f = "($f1|$f2)"
    # autoscale directive (i.e. `*` surrounded by
    # optional bounds lb < * < ub)
    as = "(($f\\s*<\\s*)?\\*(\\s*<\\s*$f)?)"
    # full range item: a floating point, or an
    # autoscale directive, or nothing
    it = "(\\s*($as|$f)?\\s*)"

    # empty range
    er = "\\[\\s*\\]"
    # full range: two colon-separated items
    fr = "\\[$it:$it\\]"

    # range regex
    rx = Regex("^\\s*($er|$fr)\\s*\$")

    if isempty(s) || occursin(rx, s)
        return true
    end

    return false
end

# Validate coordinates
function valid_coords(x,y;err=ErrorCoords(),fin=FinancialCoords())
    length(x) != length(y) && return false
    (err.valid && length(x) != length(err.ylow)) && return false
    (err.valid && !isempty(err.yhigh) && length(x) !=
        length(err.yhigh)) && return false
    (fin.valid && length(x) != length(fin.open)) && return false
    (fin.valid && length(x) != length(fin.low)) && return false
    (fin.valid && length(x) != length(fin.high)) && return false
    (fin.valid && length(x) != length(fin.close)) && return false
    return true
end


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

# Like Set, but using IdDict
mutable struct IdSet{T} <: AbstractSet{T}
    dict::IdDict{T,Nothing}

    IdSet{T}() where {T} = new(IdDict{T,Nothing}())
    IdSet{T}(s::IdSet{T}) where {T} = new(copy(s.dict))
end

IdSet{T}(itr) where {T} = union!(IdSet{T}(), itr)
IdSet() = IdSet{Any}()

copymutable(s::IdSet) = typeof(s)(s)
emptymutable(s::IdSet{T}, ::Type{U}=T) where {T,U} = IdSet{U}()
copy(s::IdSet) = typeof(s)(s)

isempty(s::IdSet) = isempty(s.dict)
length(s::IdSet)  = length(s.dict)
in(@nospecialize(x), s::IdSet) = haskey(s.dict, x)
push!(s::IdSet, @nospecialize(x)) = (s.dict[x] = nothing; s)
pop!(s::IdSet, @nospecialize(x)) = (pop!(s.dict, x); x)
pop!(s::IdSet, @nospecialize(x), @nospecialize(default)) = (x in s ? pop!(s, x) : default)
delete!(s::IdSet, @nospecialize(x)) = (delete!(s.dict, x); s)

sizehint!(s::IdSet, newsz) = (sizehint!(s.dict, newsz); s)
empty!(s::IdSet) = (empty!(s.dict); s)

filter!(f, d::IdSet) = unsafe_filter!(f, d)

function iterate(s::IdSet, state...)
    y = iterate(s.dict, state...)
    y === nothing && return nothing
    ((k, _), i) = y
    return (k, i)
end

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


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

using REPL.TerminalMenus, Test

mutable struct DynamicMenu <: TerminalMenus._ConfiguredMenu{TerminalMenus.Config}
    pagesize::Int
    pageoffset::Int
    selected::Int
    numopts::Int
    config::TerminalMenus.Config
end

TerminalMenus.numoptions(m::DynamicMenu) = m.numopts
function TerminalMenus.writeline(buf::IO, m::DynamicMenu, idx::Int, iscursor::Bool)
    print(buf, idx)
    iscursor && print(buf, '*')
end

matchback(str) = match(r"\e\[(\d+)A", str)

function linesplitter(str)
    strs = split(str, '\n')
    s1 = strs[1]
    m = matchback(s1)
    if m === nothing
        nback, startidx = 0, 1
    else
        nback = parse(Int, m.captures[1])
        startidx = m.offset+length(m.match)
    end
    strs[1] = s1[startidx:end]  # discard the portion that moves the terminal
    @test all(s->startswith(s, "\e[2K"), strs)
    return nback, replace.(map(s->s[5:end], strs), ('\r'=>"",))
end

io = IOBuffer()

menu = DynamicMenu(4, 0, -1, 8, TerminalMenus.Config())

cursor = 1
state = TerminalMenus.printmenu(io, menu, cursor; init=true)
str = String(take!(io))
@test count(isequal('\n'), str) == state
nback, strs = linesplitter(str)
@test nback == 0
@test strs == [" > 1*", "   2", "   3", "v  4"]

cursor = TerminalMenus.move_down!(menu, cursor)
@test cursor == 2
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
@test count(isequal('\n'), str) == state
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["   1", " > 2*", "   3", "v  4"]

menu.numopts = 3   # dynamically changing the number of options
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test count(isequal('\n'), str) - parse(Int, matchback(strs[end]).captures[1]) == state
@test strs[1:end-1] == ["   1", " > 2*", "   3"]

menu.numopts = 6
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 2
@test strs == ["   1", " > 2*", "   3", "v  4"]

cursor = TerminalMenus.move_down!(menu, cursor)
cursor = TerminalMenus.move_down!(menu, cursor)
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["^  2", "   3", " > 4*", "v  5"]

cursor = TerminalMenus.move_down!(menu, cursor)
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["^  3", "   4", " > 5*", "   6"]

cursor = TerminalMenus.move_down!(menu, cursor)
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["^  3", "   4", "   5", " > 6*"]

cursor = TerminalMenus.move_up!(menu, cursor)
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["^  3", "   4", " > 5*", "   6"]

cursor = TerminalMenus.page_up!(menu, cursor)
@test cursor == 1
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == [" > 1*", "   2", "   3", "v  4"]

cursor = TerminalMenus.move_down!(menu, cursor)
cursor = TerminalMenus.move_down!(menu, cursor)
cursor = TerminalMenus.move_down!(menu, cursor)
cursor = TerminalMenus.page_down!(menu, cursor)
@test cursor == menu.numopts
state = TerminalMenus.printmenu(io, menu, cursor; oldstate=state)
str = String(take!(io))
nback, strs = linesplitter(str)
@test nback == 3
@test strs == ["^  3", "   4", "   5", " > 6*"]

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module Filters

include("Filters/firfilters.jl")
export 	FIRFilter,
		createfirfilter

include("Filters/iirfilters.jl")
export 	IIRFilter,
		createiirfilter,
		butterworth,
		cheby1,
		cheby2,
		ellip,
		bessel

export 	applyfilter,
		freqz

function applyfilter(filt::IIRFilter, x::Float64)
	y = Cfloat[0];
	ccall((:iirfilt_rrrf_execute, "libliquid"), Void, (Ptr{Void}, Cfloat, Ptr{Cfloat}), filt.filterobjPtr,
		Cfloat(x), y);
	return y[1];
end

function applyfilter(filt::IIRFilter, x::Array{Float64})
	y = Array(Float64, length(x));
	for i in range(1,length(x))
		y[i] = applyfilter(filt, x[i]);
	end
	return y;
end

function freqz(filt::IIRFilter, fc::Float64)
	h_len = ccall((:iirfilt_rrrf_get_length, "libliquid"), Int, (Ptr{Void}, ) , filt.filterobjPtr);
	response = Array(Cfloat, h_len)
	ccall((:iirfilt_rrrf_freqresponse, "libliquid"), Void, (Ptr{Void} , Cfloat, Ptr{Cfloat}),
		filt.filterobjPtr, Cfloat(fc), response)
	return response;
end

end


text:
# Standard library
using Printf
using Dates
import Base: length, +, show, zero

abstract type AbstractLocation end

"""
    Geolocation (lon, lat)

Type to hold a lon/lat pair and internally converts
to -180 to 180 lon range. Checks if supplied coordinates
are within bounds [-180 180] and [-90 90] and throws
an error if not.

# Arguments
- 'lon<:Real': Longitude of location
- 'lat<:Real': Latitude of location

# Returns
- 'Geolocation': Geolocation object
"""
struct Geolocation <: AbstractLocation

    lon::Real
    lat::Real

    function Geolocation(lon, lat)

        # Convert to -180 to 180
        lon = mod((lon + 180), 360) - 180

        # Check for range
        if (-180 <= lon <= 180) & (-90 <= lat <= 90)
            new(lon, lat)
        else
            error("Location out of range [-180 180], [-90 90].")
        end

    end
end

"""
    GeolocationTime (location, time)

Type to hold both a location {Geolocation} and a {DateTime}
time which holds the measurement time.
"""
struct GeolocationTime <: AbstractLocation
    loc::Geolocation
    time::DateTime
end

"""
    Scene

This type represents a measurement, has a location and time, as
well as the associated solar and viewing geometries, the surface
SIF value and a NIRv value.

"""
struct Scene
    # What instrument?
    instrument::String
    # Potential sampling mode (nadir, glint, SAM, ..)
    mode::String
    # Location and time of scene
    loctime::GeolocationTime
    # Measured SIF value (noisy-fied truth, mimicks real measurement)
    SIF::Real
    # Modelled SIF truth (no noise)
    SIF_truth::Real
    # SIF single-sounding uncertainty
    SIF_ucert::Real
    # solar zenith angle
    SZA::Real
    # viewing zenith angle
    VZA::Real
    # NIR radiance
    NIR::Real
    # VIS radiance
    VIS::Real
    # NDVI
    NDVI::Real
    # NIRv (reflectance)
    NIRv::Real
    # NIRv (radiance)
    NIRv_radiance::Real
    # BOA irradiance in SIF units W/m2/nm/sr
    irradiance::Real
    # BOA PPFD
    PPFD::Real
    # Surface reflectance (ideally computed from SZA, VZA via BRDF)
    # such that L_out = reflectance * L_in  * 1000 / pi
    reflectance::Real
    # Total column optical depth (e.g. clouds, aerosols)
    OD::Real
end

"""
    Some design thoughts on the various sampling types

The main goal of this tool is to provide some OSSE-type insights into
how light response curves for some region of interest are sampled given
some either existing or fictional space-based instrument. Since we
want to be able to "play" with various configurations, we are keeping
the various data of a *sampling* somewhat redundant.

For example, scene locations are separately stored in several places. This
allows to separate the sampling times and repetitions from the location
data, and e.g. allows us to use a function on a sampling object which
adds a repetition.

"""
abstract type InstrumentSampling end

"""
    length (IS)

# Arguments
- 'IS::InstrumentSampling': InstrumentSampling object

# Returns
- 'Int': Number of scenes (not scene locations) within the sampling
"""
function length(IS::InstrumentSampling)
    return length(IS.scenes)
end

# GeoCarb-type sampling for a full day
struct GeostationaryFullDaySampling<:InstrumentSampling

end

# Think "repeated Granule"
struct GeostationaryIntensiveSampling<:InstrumentSampling
    # Information about how this object was created
    info::String

    # Which instrument(s)?
    instrument::Vector{String}

    # Scene locations obtained from files
    locations::Vector{Geolocation}

    # This is merely an array of scenes which are
    # obtained through files. Can be of different
    # dimensions than scenelocs
    scenes::Vector{Scene}
end

# OCO-2/3 type sampling pattern (this uses real data)
struct OCOSampling<:InstrumentSampling
    # Information about how this object was created
    info::String

    # Which instrument(s)?
    instrument::Vector{String}

    # Scene locations obtained from files
    locations::Vector{Geolocation}

    # This is merely an array of scenes which are
    # obtained through files. Can be of different
    # dimensions than scenelocs
    scenes::Vector{Scene}
end

# Orbiter with a grating-type nadir-looking swath instrument (e.g. TROPOMI, CO2M)
struct NadirSwathOrbiterSampling<:InstrumentSampling
    # Information about how this object was created
    info::String

    # Which instrument(s)?
    instrument::Vector{String}

    # Scene locations obtained from files
    locations::Vector{Geolocation}

    # This is merely an array of scenes which are
    # obtained through files. Can be of different
    # dimensions than scenelocs
    scenes::Vector{Scene}
end

# Some new combined sampling
# (fill this in when you have a better idea of what to do with mixed samplings)
struct CombinedSampling<:InstrumentSampling
    # Information about how this object was created
    info::String

    # Which instrument(s)?
    instruments::Vector{String}

    # Scene locations obtained from files
    locations::Vector{Geolocation}

    # This is merely an array of scenes which are
    # obtained through files. Can be of different
    # dimensions than scenelocs
    scenes::Vector{Scene}
end

############################################################
# Various getter functions to obtain data from scenes etc. #
############################################################

get_loc(S::Scene) = S.loctime.loc
get_loctime(S::Scene) = S.loctime
get_lon(S::Scene) = S.loctime.loc.lon
get_lat(S::Scene) = S.loctime.loc.lat
get_time(S::Scene) = S.loctime.time
get_mode(S::Scene) = S.mode
get_sza(S::Scene) = S.SZA
get_mu0(S::Scene) = cos(deg2rad(S.SZA))
get_vza(S::Scene) = S.VZA
get_nir(S::Scene) = S.NIR
get_vis(S::Scene) = S.VIS
get_sif(S::Scene) = S.SIF
get_sif_truth(S::Scene) = S.SIF_truth
get_ppfd(S::Scene) = S.PPFD
get_reflectance(S::Scene) = S.reflectance
get_irradiance(S::Scene) = S.irradiance
get_ndvi(S::Scene) = S.NDVI
get_nirv(S::Scene) = S.NIRv
get_nirv_radiance(S::Scene) = S.NIRv_radiance
get_instrument(S::Scene) = S.instrument
get_sif_ucert(S::Scene) = S.SIF_ucert

get_instrument(IS::OCOSampling) = IS.instrument
get_instrument(IS::GeostationaryIntensiveSampling) = IS.instrument

get_loc(IS::InstrumentSampling) = IS.locations
get_loc_lons(IS::InstrumentSampling) = (p -> p.lon).(IS.locations)
get_loc_lats(IS::InstrumentSampling) = (p -> p.lat).(IS.locations)

get_time(IS::InstrumentSampling) = get_time.(IS.scenes)

get_lon(IS::InstrumentSampling) = get_lon.(IS.scenes)
get_lat(IS::InstrumentSampling) = get_lat.(IS.scenes)
get_sza(IS::InstrumentSampling) = get_sza.(IS.scenes)
get_vza(IS::InstrumentSampling) = get_vza.(IS.scenes)
get_nir(IS::InstrumentSampling) = get_nir.(IS.scenes)
get_vis(IS::InstrumentSampling) = get_vis.(IS.scenes)
get_sif(IS::InstrumentSampling) = get_sif.(IS.scenes)
get_sif_truth(IS::InstrumentSampling) = get_sif_truth.(IS.scenes)
get_mode(IS::InstrumentSampling) = get_mode.(IS.scenes)
get_ppfd(IS::InstrumentSampling) = get_ppfd.(IS.scenes)
get_reflectance(IS::InstrumentSampling) = get_reflectance.(IS.scenes)
get_irradiance(IS::InstrumentSampling) = get_irradiance.(IS.scenes)
get_ndvi(IS::InstrumentSampling) = get_ndvi.(IS.scenes)
get_nirv(IS::InstrumentSampling) = get_nirv.(IS.scenes)
get_nirv_radiance(IS::InstrumentSampling) = get_nirv_radiance.(IS.scenes)
get_sif_ucert(IS::InstrumentSampling) = get_sif_ucert.(IS.scenes)

# Define some pretty printing for our types

function show(io::IO, S::Scene)
    print(io, [get_lon(S), get_lat(S)])
end

function show(io::IO, ::MIME"text/plain", S::Scene)
    println(io, @sprintf("%s scene at [%0.4f, %0.4f]", get_instrument(S), get_lon(S), get_lat(S)))
    println(io, Dates.format(get_time(S), "yyyy u dd, HH:MM:SS"))
    println(io, @sprintf("VZA: %0.3f deg", get_vza(S)))
    println(io, @sprintf("SZA: %0.3f deg", get_sza(S)))
    println(io, @sprintf("Reflectance@755: %0.3f", get_reflectance(S)))
    println(io, @sprintf("NDVI: %0.3f", get_ndvi(S)))
    println(io, @sprintf("NIRv: %0.3f", get_nirv(S)))
end

# Define an addition (merger) of two InstrumentSamplings
# (difficult? this should work for all instrument types?)
function +(IS1::T1, IS2::T2) where {T1<:InstrumentSampling, T2<:InstrumentSampling}

    newinfo = "$(IS1.info) + $(IS2.info)"

    # We do not want any duplicates in the merged IS object,
    # so that we don't overcount measurements
    # ------------------------------------------------------

    # Simple set union works just fine, and is fast
    newlocations = union(IS1.locations, IS2.locations)

    # Scenes are a little more tricky:
    # It could happen that the same exact location is sampled by
    # multiple instruments at the same time. So we consider a scene
    # a duplicate if the location, the time, and the VZA is the same.

    # We are therefore performing a negative search - find scenes for
    # with the time is the same (IS.scenes[].loctime.time is time-ordered),
    # and then check for those duplicate scenes whether the other
    # quantities are the same. If they are, those scenes of IS2 are
    # added to the exclusion list which are not taken up into the
    # new merged IS object.

    IS1_instrument = get_instrument(IS1)
    IS2_instrument = get_instrument(IS2)

    IS1_times = get_time.(IS1.scenes)
    IS2_times = get_time.(IS2.scenes)

    intersect_times = intersect(IS1_times, IS2_times)
    srt_IS1 = searchsortedfirst.(Ref(IS1_times), intersect_times)
    srt_IS2 = searchsortedfirst.(Ref(IS2_times), intersect_times)

    exclude_list = Int[]

    # Loop through scenes, determing if location+time
    # and viewing zenith are the same. If so, push them into
    # exclusion list
    for i in 1:length(srt_IS1)

        s1 = IS1.scenes[srt_IS1[i]]
        s2 = IS2.scenes[srt_IS2[i]]

        check::Bool = (
            (get_loctime(s1) == get_loctime(s2)) &
            (get_vza(s1) == get_vza(s2))
        )

        if check
            push!(exclude_list, srt_IS2[i])
        end
    end

    # From the exclude list, build an include list via
    # set operations.
    include_list = setdiff(
        range(1, length=length(IS2)),
        exclude_list
    )

    # Depending on what type combinations we have, create a new object
    if T1 == T2

        # If both InstrumentSamplings are of the same type
        # (i.e. from the satellite kind), we create a new
        # object of the same type.

        if IS1_instrument == IS2_instrument
            # If the two instrument descriptors are the same,
            # just use the same instrument descriptor
            newinstrument = IS1_instrument
        else
            # Otherwise, build a new array with unique entries only.
            # For example ["oco2"] and ["oco3"] will become
            # ["oco2", "oco3"].
            newinstrument = sort(unique(
                vcat(IS1_instrument, IS2_instrument)
            ))
        end

        # New scenes will be also sorted by time, even though it will
        # shuffle instruments.
        newscenes = vcat(IS1.scenes, IS2.scenes[include_list])
        newsort = sortperm(get_time.(newscenes))

        newsampling = T1(
            newinfo,
            newinstrument,
            newlocations,
            newscenes[newsort] # Pass time-sorted scenes
        )

        # Last check if scenes are truly sorted by time
        if issorted(get_time(newsampling))
            return newsampling
        else
            @error "Resulting InstrumentSampling is not time-sorted. Debug."
            return nothing
        end

    else

        @error "Not implemented yet!"
        return nothing

    end

end

abstract type TemporalSampling end

# Some regular temporal sampling, think "X days/weeks/months"
struct RegularTemporalSampling<:TemporalSampling
    # Time in [s] between two aggregation boundaries
    period::Real
end

function HourlySampling(interval::Number)
    # There are 60 * 60 = 3600 seconds in an hour
    return RegularTemporalSampling(3600 * interval)
end

function WeeklySampling(interval::Number)
    # There are 60 * 60 * 24 * 7 = 604800 seconds in a week
    return RegularTemporalSampling(604800 * interval)
end

function DailySampling(interval::Number)
    # There are 60 * 60 * 24 = 86400 seconds in a day
    return RegularTemporalSampling(86400 * interval)
end

# In case we want to investigate some more complex irregular
# temporal samplings.
# (e.g. correlations with periods of stable precipitation etc.)
struct IrregularTemporalSampling<:TemporalSampling
    time_boundaries::Array{DateTime, 1}
end

abstract type SpatialSampling end

# No further subsampling in spatial dimension
# treat whole ROI of InstrumentSampling object
# as one.
struct FullROI<:SpatialSampling
end

# Regular X by Y grid cells in lon / lat
struct RegularGridCells<:SpatialSampling
    delta_lon::Real
    delta_lat::Real
end

"""
    Aggregate

Represents a SIF aggregate
"""
struct Aggregate
    # Number of measurements in aggregate
    N::Int
    # Scene indices, in case you want to refer back to other data
    scene_idx::Vector{Int}
    # Scenes
    scenes::Vector{Scene}
    # Start time
    start_time::DateTime
    # End time
    end_time::DateTime
end

# Struct for holding surface data which is then accessed
# by sampling functions

abstract type SurfaceData end

struct VNPData <: SurfaceData
    data::Array{Int, 4} ## Lat, Lon, DOY, variable
    lon_bounds::Vector{Real}
    lat_bounds::Vector{Real}
    time_bounds::Vector{DateTime}
    data_year::Int
    # How many entries per coordinate?
    names::Vector{String}
    fill_values::Vector{Real}
end

# Struct to hold SIF data given some lon/lat boundaries
# that can then be sampled by the sampling functions

abstract type SIFData end

struct TROPOMISIFData <: SIFData
    data::Array{Float32, 3} ### Lat, Lon, Month
    lon_bounds::Vector{Real}
    lat_bounds::Vector{Real}
    time_bounds::Vector{DateTime}
    fill_values::Real
end


text:
function solve()
    s = readline()
    if s[1:2:end] == lowercase(s[1:2:end]) && s[2:2:end] == uppercase(s[2:2:end])
        "Yes"
    else
        "No"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function phi_step_modelA(phi, V, t, dt, lambda, r, g, stencil_bc, xstep, T, D)
    noise = sqrt(T * dt / xstep) * sampler(size(phi))
    d = length(size(phi))
    out = similar(phi)
    for x in stencil_bc
        out[x[1]] = (
            phi[x[1]] -
            D *
            dt *
            (
                -1.0 / (xstep^2) * (sum(phi[i] for i in x[2]) - 2 * d * phi[x[1]]) +
                g * phi[x[1]]^2 +
                r +
                lambda * V[1][x[1]]
            ) *
            phi[x[1]] + noise[x[1]]
        )
    end
    @. phi = out
end

max_stars_repo_name: phyjonas/ColloidFluctuatingField


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(abs_tol=1e-15, orderT=8, orderQ=1, max_steps=50_000)),
    "X0"   => Hyperrectangle(low=[1.3, 1.3, -0.2, -0.2], high=[1.7, 1.7, 0.2, 0.2]),
    "info" => Dict(
        "name" => "robot_arm",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])


text:
using Primes
using Base.Iterators: drop, product

proper_divisors(n) = begin
    if n == 1
        return [1]
    end
    factors_of_n = factor(n)
    nfactors = prod(x -> x[2] + 1, factors_of_n)
    res = Vector{Int}(undef, nfactors)
    res .= 1
    i = 1
    for (factor, power) in factors_of_n
        if i == 1
            res[1] = 1
            res[2:power+1] .= factor .^ (1:power)
            i = power + 1
        else
            next_few = i * power
            next_few_vals = [
                existing * more
                for (existing, more) in product(res[1:i], factor .^ (1:power))
            ]
            res[i+1:i+next_few] = vec(next_few_vals)
            i += next_few
        end
    end
    res[1:end-1]
end


text:
using TextAnalysis: AveragePerceptron

@testset "Average Perceptron Tagger" begin
    tagger = PerceptronTagger(false)

    @testset "Basic" begin
        @test typeof(tagger.classes) == Set{Any}
        @test length(tagger.classes) == 0
        @test typeof(tagger.model) == AveragePerceptron

        fit!(tagger, [[("today","NN"),("is","VBZ"),("good","JJ"),("day","NN")]])
        @test length(keys(tagger.model.weights)) == 51
        @test tagger.classes == tagger.model.classes == Set(["JJ", "VBZ", "NN"])
    end

    @testset "Average Perceptron Tagger (pretrained)" begin
        tagger = PerceptronTagger(true)

        @test typeof(tagger.classes) == Set{Any}
        @test length(tagger.classes) == 75
        @test typeof(tagger.model) == AveragePerceptron
    end

    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    @testset "Tagging over sentences and documents" begin
        tagger = PerceptronTagger(true)
        text = "This is a text"
        @test tagger(text) == predict(tagger, text)

        sd = StringDocument(text)
        @test length(predict(tagger, text)) == 4
        @test length(predict(tagger, sd)) == 4

        text2 = read(sample_file, String)
        fd = FileDocument(sample_file)
        @test length(predict(tagger, fd)) == length(predict(tagger, text2))

        td = TokenDocument(text)
        @test length(predict(tagger, td)) == 4
    end
end


text:
# ordered dict

import Base: haskey, get, get!, getkey, delete!, push!, pop!, empty!,
             setindex!, getindex, sizehint, length, isempty, start,
             next, done, keys, values, setdiff, setdiff!,
             union, union!, intersect, isequal, filter, filter!,
             hash, eltype

# This is just a simple wrapper around a HashDict, which is a modified
# implementation of the Dict implementation in Base allowing ordering
# information to be maintained.

# In particular, the HashDict stored in an OrderedDict is a
# HashDict{K,V,Ordered}

# A HashDict{K,V,Unordered} would be equivalent to Base.Dict

immutable OrderedDict{K,V} <: Associative{K,V}
    d::HashDict{K,V,Ordered}

    OrderedDict() = new(HashDict{K,V,Ordered}())
    OrderedDict(kv) = new(HashDict{K,V,Ordered}(kv))
    OrderedDict(ks,vs) = new(HashDict{K,V,Ordered}(ks,vs))
end

OrderedDict() = OrderedDict{Any,Any}()

OrderedDict{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V}) = OrderedDict{K,V}(ks,vs)
OrderedDict{K,V}(::Type{K},::Type{V}) = OrderedDict{K,V}()
OrderedDict(ks,vs) = OrderedDict{eltype(ks),eltype(vs)}(ks, vs)

OrderedDict{K,V}(kv::AbstractArray{(K,V)}) = OrderedDict{K,V}(kv)

## Functions

## Most functions are simply delegated to the wrapped HashDict

@delegate OrderedDict.d [ haskey, get, get!, getkey, delete!, pop!,
                          empty!, setindex!, getindex, sizehint,
                          length, isempty, start, next, done, keys,
                          values ]

similar{K,V}(d::OrderedDict{K,V}) = OrderedDict{K,V}()
in{T<:OrderedDict}(key, v::Base.KeyIterator{T}) = key in keys(v.dict.d.d)


text:
G(q) = 2*[-q[2:4] q[1]*I - skew(q[2:4])]

q = Quaternion(rand(Quat))
qd = Quaternion(rand(Quat))
q_ = SVector(q)
G(q_)
R = [zeros(3,1) Diagonal(I,3)]

λ = rand(3)
λhat_ = [0;λ]
λhat = Quaternion(λhat_)
ω = rand(3)
qd = q*Quaternion(0.,ω/2)
qd_ = SVector(qd)

# dGdq
2R*Rmult(qd)*N == -G(qd_)
2λhat_'Rmult(qd)*N == -λ'G(qd_)
2λhat_'Rmult(qd)*N == -λ'G(qd_)
2N*Rmult(qd)'λhat_ == -G(qd_)'λ
2*SVector(qd*inv(λhat)) ≈ -G(qd_)'λ

2λ'R*Lmult(inv(q)) == (2*Lmult(q)*λhat_)'
λ'ω ≈ 2λ'R*Lmult(inv(q))*qd_
λ'ω ≈ 2*λhat_'SVector(inv(q)*qd)
λ'ω
2*qd_'SVector(q*λhat)

λhat_'Rmult(qd)
(qd_'Rmult(λhat))*N

N*qd_

# First term: λ'Gprime(qd)
ForwardDiff.jacobian(x->G(x)*qd_,q_) == -G(qd_)
-λ'G(qd_) == -2λ'R*Lmult(inv(qd))
-λ'G(qd_) ≈ λ'R*Lmult(Quaternion(0.,ω))*Lmult(q)'

# Second term: -λdot'G(q)
λ'G(q_) == 2λ'R*Lmult(inv(q))

# Third term: λ'dGdq*qdot
ForwardDiff.jacobian(x->G(x)*qd_,q_) == -2*R*Lmult(inv(qd))
λ'ForwardDiff.jacobian(x->G(x)*qd_,q_) ≈ λ'R*Lmult(Quaternion(0.,ω))*Lmult(q)'
G(q_) == 2*R*Lmult(inv(q))
-2*SVector(qd*Quaternion(0.,λ)) ≈ 2*[0 λ'; -λ skew(λ)]*qd_
ForwardDiff.gradient(x->λ'G(x)*qd_,q_) == 2*[0 λ'; -λ skew(λ)]*qd_
ForwardDiff.gradient(x->λ'G(x)*qd_,q_)
 #≈ -2*Lmult(qd)*R'λ
-2*Lmult(qd)*R'λ ≈ -Lmult(q)*Lmult(Quaternion(0.,ω))*R'λ
-2*Lmult(qd)*R'λ ≈ (λ'R*Lmult(Quaternion(0.,ω))*Lmult(q)')'
-2*Lmult(qd)*R'λ ≈ (-2*λ'R*Lmult(qd)')'

Lmult(Quaternion(0.,ω/2))'Lmult(q)' ≈ Lmult(inv(qd))

Lmult(inv(q))*q_

-2*SVector(qd*Quaternion(0.,λ))

Lmult(q)' == Lmult(inv(q))

Lmult(q)*[0;ω/2]
Lmult(inv(q*Quaternion(0.,ω/2)))

2λ'R*Lmult(inv(q))
2λ'R*SVector(inv(q))

N = Diagonal([1,-1,-1,-1])
N*q_


text:
module Dialog

using Lynx
using Lynx: Container
using Lynx: Widget

import Gtk
import Gtk.GConstants
import Gtk.GAccessor

export ask!
export error!
export file!
export input!
export save!
export warn!

"""
        file!(title::String; filters = String[], parent = NullContainer()) -> String

Opens a `file` dialog. Suitable for `open/read` file operations.

`filters` are used to filter the types of files the user can choose.
Example: `*.png`
"""
function file!(title::String; parent::Widget = NullContainer(), filters::Vector{String} = String[])
    return Gtk.open_dialog(title, parent.widget, filters)
end
"""
        save!(title::String, default = "Untitled", path = pwd(), parent = NullContainer()) -> String

Opens a dialog for saving a `file`.

- `default` - the default file name.
- `path` - the path of the directory which the dialog opens with.
- `parent` - the window the dialog belongs.

Returns a string containing the full path the user choose
to save or a empty string if the dialog was cancelled.
"""
function save!(title::String, default::String = "Untitled", path::String = pwd(), parent::Widget = NullContainer())
    return save_file_dialog(title, default, path, parent.widget)
end
"""
        info!(message::String; parent = NullContainer())

Display a informative `message`.
"""
function info!(message::String; parent::Widget = NullContainer())
    return Gtk.info_dialog(message, parent.widget)
end
"""
        warn!(message::String; parent = NullContainer())

Warn the user with a `message`.
"""
function warn!(message::String; parent::Widget = NullContainer())
    return Gtk.warn_dialog(message, parent.widget)
end
"""
        error!(message::String; parent = NullContainer())

Display a error `message`.
"""
function error!(message::String; parent::Widget = NullContainer())
    return Gtk.error_dialog(message, parent.widget)
end
"""
        ask!(message, yes::String, no::String; parent = NullContainer())

Opens a dialog asking the user for confirmation.
The dialog contains `yes` and `no` buttons which both can be customized.
"""
function ask!(message::String, no::String = "No", yes::String = "Yes"; parent::Widget = NullContainer())
    return Gtk.ask_dialog(message, no, yes, parent.widget)
end
"""
        input!(message, default, buttons = ("Cancel", "Accept", 1); parent = NullContainer()) -> String

Prompts for text input and returns what the user typed.
"""
function input!(message::AbstractString, buttons = ("Cancel", "Accept"); default::AbstractString = "",
    parent::Widget = NullContainer())
    btns = ((buttons[1], 0), (buttons[2], 1))
    return Gtk.input_dialog(message, default, btns, parent.widget)
end

# Extended Gtk Functions

function save_file_dialog(title::AbstractString, default::AbstractString,
    path::AbstractString, parent = Gtk.GtkNullContainer(); kwargs...)
    dlg = Gtk.GtkFileChooserDialog(title, parent, GConstants.GtkFileChooserAction.SAVE,
                                (("_Cancel", GConstants.GtkResponseType.CANCEL),
                                 ("_Save",   GConstants.GtkResponseType.ACCEPT)); kwargs...)
    dlgp = Gtk.GtkFileChooser(dlg)
    GAccessor.current_name(dlgp, default)
    GAccessor.current_folder(dlgp, path)
    response = Gtk.run(dlg)
    filename = ""
    if response == Gtk.GtkResponseType.ACCEPT
        filename = Gtk.bytestring(GAccessor.filename(dlgp))
    end
    Gtk.destroy(dlg)
    return filename
end

end # Dialog


text:
using JLD2, GLMakie

function kernel_plot(kernel; colormap = :thermal, colormap2 = :balance, resolution = (1800, 1300))
    (; Kˣˣ, Kˣᶻ, Kᶻˣ, Kᶻᶻ) = kernel

    fig = Figure(resolution = resolution)
    titlestring = "Kˣˣ"
    ax1 = Axis(fig[1, 1], title = titlestring, titlesize = 30)
    titlestring = "Kˣᶻ"
    ax2 = Axis(fig[1, 3], title = titlestring, titlesize = 30)
    titlestring = "Kᶻˣ"
    ax3 = Axis(fig[2, 1], title = titlestring, titlesize = 30)
    titlestring = "Kᶻᶻ"
    ax4 = Axis(fig[2, 3], title = titlestring, titlesize = 30)

    hm1 = heatmap!(ax1, Kˣˣ, colormap = colormap)
    ax1.yreversed = true

    hm2 = heatmap!(ax2, Kˣᶻ, colormap = colormap2)
    ax2.yreversed = true

    hm3 = heatmap!(ax3, Kᶻˣ, colormap = colormap2)
    ax3.yreversed = true

    hm4 = heatmap!(ax4, Kᶻᶻ, colormap = colormap)
    ax4.yreversed = true

    Colorbar(fig[1, 2], hm1, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
        labelsize = 30, ticksize = 25, tickalign = 1,)
    Colorbar(fig[2, 2], hm2, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
        labelsize = 30, ticksize = 25, tickalign = 1,)
    Colorbar(fig[1, 4], hm3, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
        labelsize = 30, ticksize = 25, tickalign = 1,)
    Colorbar(fig[2, 4], hm4, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
        labelsize = 30, ticksize = 25, tickalign = 1,)
    display(fig)
    return fig
end

function local_diffusivity_plot(analytic, numerical, z; xlimits = nothing)
    (; a_κˣˣ, a_κˣᶻ, a_κᶻˣ, a_κᶻᶻ) = analytic
    (; κˣˣ, κˣᶻ, κᶻˣ, κᶻᶻ) = numerical

    options = (; ylabel = "z", ylabelsize = 32,
        xlabelsize = 32, xgridstyle = :dash, ygridstyle = :dash, xtickalign = 1,
        xticksize = 30, ytickalign = 1, yticksize = 30,
        xticklabelsize = 30, yticklabelsize = 30)

    fig = Figure(resolution = (1800, 1300), title = "Local Operators")
    titlestring = "Kˣˣ"
    ax1 = Axis(fig[1, 1]; options..., title = titlestring, titlesize = 30)
    titlestring = "Kˣᶻ"
    ax2 = Axis(fig[1, 2]; options..., title = titlestring, titlesize = 30)
    titlestring = "Kᶻˣ"
    ax3 = Axis(fig[2, 1]; options..., title = titlestring, titlesize = 30)
    titlestring = "Kᶻᶻ"
    ax4 = Axis(fig[2, 2]; options..., title = titlestring, titlesize = 30)

    plot_string_1 = "analytic"
    plot_string_2 = "numerical"

    ln1 = lines!(ax1, a_κˣˣ, z, color = :red)
    sc1 = scatter!(ax1, κˣˣ, z, color = :blue)
    axislegend(ax1, [ln1, sc1], [plot_string_1, plot_string_2], position = :rc, labelsize = 30)

    ln2 = lines!(ax2, a_κˣᶻ, z, color = :red)
    sc2 = scatter!(ax2, κˣᶻ, z, color = :blue)
    # axislegend(ax2, [ln2, sc2], [plot_string_1, plot_string_2], position = :rt, labelsize = 30)
    hideydecorations!(ax2)

    ln3 = lines!(ax3, a_κᶻˣ, z, color = :red)
    sc3 = scatter!(ax3, κᶻˣ, z, color = :blue)
    # axislegend(ax3, [ln3, sc3], [plot_string_1, plot_string_2], position = :rb, labelsize = 30)

    ln4 = lines!(ax4, a_κᶻᶻ, z, color = :red)
    sc4 = scatter!(ax4, κᶻᶻ, z, color = :blue)
    hideydecorations!(ax4)
    # axislegend(ax4, [ln4, sc4], [plot_string_1, plot_string_2], position = :rt, labelsize = 30)
    # colgap!(ax1, 10)
    # rowgap!(ax1, 10)

    for ax in [ax1, ax2, ax3, ax4]
        if isnothing(xlimits)
        else
            xlims!(ax, xlimits)
        end
    end

    display(fig)
    return fig
end

function file_kernel_plot(file)
    jlfile = jldopen("data/" * file, "a+")
    # jlfile = jldopen("data/nonlocal_more_velocity.jld2")
    EF¹¹ = jlfile["diffusivity"]["K11"]
    EF¹² = jlfile["diffusivity"]["K12"]
    EF²¹ = jlfile["diffusivity"]["K21"]
    EF²² = jlfile["diffusivity"]["K22"]
    x = jlfile["grid"]["x"]
    z = jlfile["grid"]["z"]
    N = length(x)
    M = length(z)
    tmpE = copy(EF¹¹)
    EF¹¹ = reshape(permutedims(reshape(EF¹¹, (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
    EF¹² = reshape(permutedims(reshape(EF¹², (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
    EF²¹ = reshape(permutedims(reshape(EF²¹, (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
    EF²² = reshape(permutedims(reshape(EF²², (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
    kernel_matrix = (; Kˣˣ = EF¹¹, Kˣᶻ = EF¹², Kᶻˣ = EF²¹, Kᶻᶻ = EF²²)
    fig = kernel_plot(kernel_matrix)
    return fig
end

function file_local_diffusivity_plot(file)
    jlfile = jldopen("data/"*file, "a+")

    κ¹¹ = jlfile["averagelocaldiffusivity"]["κ11"]
    κ¹² = jlfile["averagelocaldiffusivity"]["κ12"]
    κ²¹ = jlfile["averagelocaldiffusivity"]["κ21"]
    κ²² = jlfile["averagelocaldiffusivity"]["κ22"]
    z = jlfile["grid"]["z"][:]
    x = jlfile["grid"]["x"][:]
    ##
    # Contruct local diffusivity estimate from flow field and γ
    γ = jlfile["parameters"]["γ"]
    ω = jlfile["parameters"]["ω"]
    u¹ = jlfile["velocities"]["u¹"]
    u² = jlfile["velocities"]["u²"]
    v¹ = jlfile["velocities"]["v¹"]
    v² = jlfile["velocities"]["v²"]

    # integrating out the first dimension amounts to assuming that the dominant mode
    # is the zero'th.
    N = length(x)
    scale = (2 * (ω^2 + γ^2))
    analytic_κ¹¹ = sum((γ * u¹ .* u¹ + γ * u² .* u²) ./ scale, dims = 1)[:] ./ N
    analytic_κ¹² = sum((u¹ .* (γ * v¹ + ω * v²) + u² .* (γ * v² - ω * v¹)) ./ scale, dims = 1)[:] ./ N
    analytic_κ²¹ = sum((v¹ .* (γ * u¹ + ω * u²) + v² .* (γ * u² - ω * u¹)) ./ scale, dims = 1)[:] ./ N
    analytic_κ²² = sum((γ * v¹ .* v¹ + γ * v² .* v²) ./ scale, dims = 1)[:] ./ N
    analytic = (; a_κˣˣ = analytic_κ¹¹, a_κˣᶻ = analytic_κ¹², a_κᶻˣ = analytic_κ²¹, a_κᶻᶻ = analytic_κ²²)
    numerical = (; κˣˣ = κ¹¹, κˣᶻ = κ¹², κᶻˣ = κ²¹, κᶻᶻ = κ²²)
    fig = local_diffusivity_plot(analytic, numerical, z[:])
    return fig
end

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
# load the relevant libraries
using CSV, DataFrames

# read the data set - note that you need a special package for this
msleep = CSV.read("mammal_sleep.CSV", DataFrame)
# Look at the DataFrame in the console - what are the element types
# of the columns? Notice that some have a little questionmark after
# - why do you think that is? (hint - it's not because Julia is in doubt)

# Exercise:
# Try running some basic descriptor functions for DataFrames
# including typeof, nrow, describe, names - and try indexing into the dataframe

# You can get a column from a dataframe in three ways
msleep[:, :awake] # this copies the column
msleep[!, :awake] # this gets a link to the column
msleep.awake      # this is the same as [!, :awake]

# Exercise
# Select multiple columns and rows
# You can use e.g. `Not(1)` to exclude index one,
# and the `Between(a,b)` function to choose multiple columns
# Get a dataframe with all rows except the first and
# all sleep/awakeness related columns

# You can also select multiple columns by using a regular expression
# regular expression search strings are created by `r"mysearchstring"`
# E.g. try
select(msleep, r"sleep_")

# Exercise
# use the `groupby` function to create a dataframe grouped by
# order and "vore" called `gdf`

# Split-apply-combine
# Let's use this to find the mean sleep total
# for the mean function we need to load the standard library "Statistics"
using Statistics
combine(gdf) do sdf
    mean_sleep = mean(sdf.sleep_total)
end

# We can use the Chain package to chain operations
# much like the "pipe" function in dplyr
# Exercise - what does this syntax do?
using Chain
@chain msleep begin
    groupby([:order, :vore])
    combine(nrow, :sleep_total => mean)
end

# Exercise
# adapt the code above so the dataframe ends up being
# sorted by the sleep_total_mean column

# This transforms to long format
@chain msleep begin
    groupby([:order, :vore])
    combine(nrow)
    unstack(:order, :vore, :nrow, allowmissing = true)
end

# Exercise
# Which two orders have the lowest mean body weight?

using DataFramesMeta

# The DataFramesMeta package provides metaprogramming
# macros, that lets you use a syntax very similar to dplyrs

# dplyr       DataFramesMeta    What it does

# select      @select           Selects columns
# mutate      @transform        Calculates new columns
# filter      @subset           Selects rows by a condition
# arrange     @orderby          Sorts rows
# summarize   @combine          Creates summary statistisc
# group_by    groupby           Groups

# You can select columns using this simple] format
@select msleep :name :order

@subset msleep :order .== "Rodentia" # here is how you would subset/filter rows
                                     # note the dot in `.==`
                                     # this applies to the whole column
@rsubset msleep :order == "Rodentia" # you can use the `@rsubset` (note the added `r`)
                                     # to work on a row-by-row basis
@rsubset msleep begin                # chain multiple commands like this
     :order == "Rodentia"
     :vore == "herbi"
end

# Exercise
# Selecting only the needed columns, sort the dataframe by order
# and total sleep, retaining only those with at least 16 hours of sleep

# Exercise
# Redo the exercise from line 54 with DataFramesMeta

# use @combine to calculate the mean, minimum and maximum sleep
# across different mammal orders

max_stars_repo_path: Materials Tuesday morning/Dataframes.jl

max_stars_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_issues_repo_path: Materials Tuesday morning/Dataframes.jl

max_issues_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_forks_repo_path: Materials Tuesday morning/Dataframes.jl

max_forks_repo_name: EcoJulia/Julia-PhD-course-Copenhagen


text:
################################################################################
#
#	CONCRETE TYPE
#
#   EnergyManifold{hamiltonian}
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct EnergyManifold{H} <: AbstractEnergyManifold{H}

    # the hamiltonian
    h :: H

    # energy of the cut
    E_cut :: Float64

    # the k points
    k_points :: Vector{Vector{Float64}}

end

# export the type
export EnergyManifold

################################################################################
#
#	INTERFACING / ACCESSING ENERGY MANIFOLDS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# obtaining the hamiltonian
function hamiltonian(
            em :: EnergyManifold{H}
        ) :: H where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the internal field
    return em.h
end

# obtaining the k points of the manifold
function kpoints(
            em :: EnergyManifold{H}
        ) :: Vector{Vector{Float64}} where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the internal field
    return em.k_points
end

# obtaining the energy at which the manifold cuts
function energy(
            em :: EnergyManifold{H}
        ) :: Float64 where {L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the internal field
    return em.E_cut
end

# recalculate the energy manifold (k points)
function recalculate!(
            em       :: EnergyManifold{H},
            n_points :: Integer
            ;
            hs_mesh_resolution :: Real = 100,
            bounds_lower :: Vector{<:Real} = ones(10) .* -2 .*pi,
            bounds_upper :: Vector{<:Real} = ones(10) .*  2 .*pi,
            refold_to_first_BZ :: Bool = true,
            precision_energy :: Real = 1e-6,
            kwargs...
        ) where {L,LS,D,S<:AbstractSite{LS,D},B,UC<:AbstractUnitcell{S,B},HB,H<:AbstractHamiltonian{L,UC,HB}}

    # list of k values that contribute to the BZ
    k_values = Vector{Vector{Float64}}(undef, n_points)

    # obtain the hamiltonian
    h  = hamiltonian(em)
    uc = unitcell(h)

    # obtain the energy for the cut
    energy_cut = energy(em)

    # the current index
    current_index = 1

    # maybe build
    if hs_mesh_resolution > 0
        # build high symmetry mesh
        hs_mesh = shuffle(getPointsHighSymmetryMesh(getBrillouinZone(uc), hs_mesh_resolution))
        # check all points in the mesh
        for k in hs_mesh
            # check energy difference to energy cut
            if minimum(abs.(eigvals(matrixAtK(h,k)).-energy_cut)) < precision_energy
                k_values[current_index] = k
                current_index += 1
                if current_index > n_points
                    break
                end
            end
        end
    end

    # FIND ALL REMAINING K POINTS (RANDOM STARTING POINTS)
    @simd for i in current_index:n_points
        # build starting point
        k = rand(D)
        @simd for j in 1:D
            @inbounds k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end
        # find the k value
        @inbounds k_values[i] = findKAtEnergy(
                                    h,
                                    energy_cut,
                                    k
                                    ;
                                    precision_energy = precision_energy,
                                    kwargs...
                                )
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        ruc = getReciprocalUnitcell(uc)
        @simd for i in 1:n_points
            @inbounds shiftToFirstBZ!(ruc, k_values[i])
        end
    end

    # set the list
    em.k_points = k_values

    # return nothing
    return nothing
end

# convenience function for obtaining the manifold
function getEnergyManifold(
            h     :: H,
            e_cut :: Real,
            n     :: Int64 = 1000
            ;
            recalculate :: Bool = true,
            kwargs ...
        ) :: EnergyManifold{H} where {L,LS,D,S<:AbstractSite{LS,D},B,UC<:AbstractUnitcell{S,B},HB,H<:AbstractHamiltonian{L,UC,HB}}

    # create a new object
    em = EnergyManifold{H}(h, e_cut, Vector{Float64}[])

    # recalculate the numerical energy values
    if recalculate
        recalculate!(em, n; kwargs...)
    end

    # return the new object
    return em
end

# Convinience function if only unitcell and the bond hamiltonian is known (uses concrete Hamiltonian)
function getEnergyManifold(
        uc    :: UC,
        hb    :: HB,
        e_cut :: Real,
        n     :: Int64 = 1000
        ;
        recalculate :: Bool = true,
        kwargs ...
    ) :: EnergyManifold{Hamiltonian{L,UC,HB}} where {L,N,LS,D,S<:AbstractSite{LS,D},B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS}}

    # create a new object
    em = EnergyManifold{Hamiltonian{L,UC,HB}}(Hamiltonian(uc,hb), e_cut, Vector{Float64}[])

    # recalculate the numerical energy values
    if recalculate
        recalculate!(em, n; kwargs...)
    end

    # return the new object
    return em
end

# Convinience function if only path and unitcell but not even the bond hamiltonian is known (uses concrete Hamiltonian and simple hopping hamiltonian)
function getEnergyManifold(
            uc    :: UC,
            e_cut :: Real,
            n     :: Int64 = 1000
            ;
            recalculate :: Bool = true,
            kwargs ...
        ) :: EnergyManifold{Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}} where {D,RP<:AbstractReciprocalPoint{D}, P<:AbstractReciprocalPath{RP}, L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B}}

    # create a new object
    em = EnergyManifold{Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}}(Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}(uc, getHoppingHamiltonianSimple(uc)), e_cut, Vector{Float64}[])

    # recalculate the numerical energy values
    if recalculate
        recalculate!(em, n; kwargs...)
    end

    # return the new object
    return em
end

# export the convenience function
export getEnergyManifold


text:
# This file calculates the third application in Y.-S. Chan, A. C. Fannjiang, and G. H. Paulino,
# Integral equations with hypersingular kernels -- theory and applications to fracture mechanics,
# Int. J. Eng. Sci., 41:683--720, 2003.

using ApproxFun, SingularIntegralEquations

x = Fun(identity)
w = 1/sqrt(1-x^2)
d = domain(x)
d2 = d^2
B = [Dirichlet(d);Neumann(d)]
H2 = Hilbert(d,2)
H4 = Hilbert(d,4)
ϵ = 0.2
L = -6ϵ^2*H4[w] + H2[w]
f = -Fun(one)

uSIE = [B;L]\[[0,0],[0,0],f]
@time uSIE = [B;L]\[[0,0],[0,0],f]
println("The length of uSIE is: ",ncoefficients(uSIE))
println("The extrema of uSIE are: ",extrema(uSIE))

test0 = uSIE/(1-x^2)^2
test = Fun(x->test0(x),Ultraspherical(1,d))
temp = 0.0
for i=1:ncoefficients(test)
    temp += i*test.coefficients[i]
end
println("The normalized generalized Stress Intensity Factors are: ",3ϵ*temp)