text:
##### Beginning of file

module Run # Begin submodule RemoveLFS.Run

__precompile__(true)

import ArgParse
import Dates
import HTTP
import Pkg
import TimeZones

import ..Types
import ..Utils
import ..Common

import ..delayederror
import ..process_delayed_error_list

function run_removelfs_snapshots!!(
        ;
        git_lfs_repos::AbstractDict,
        dst_provider,
        git_user_name,
        git_user_email,
        include_branches,
        exclude_branches,
        task::String = "all",
        is_dry_run::Bool = false,
        do_not_push_to_these_destinations::Vector{String} =
            String[],
        do_not_try_url_list::Vector{String} =
            String[],
        try_but_allow_failures_url_list::Vector{String} =
            String[],
        time_zone::Dates.TimeZone =
            Dates.TimeZone("America/New_York"),
        )::Nothing
    @info("Running RemoveLFS.Run.run_removelfs_snapshots!!")

    all_src_dest_pairs::Vector{Types.SrcDestPair} = Types.SrcDestPair[
        Types.SrcDestPair(
            ;
            source_url = x["source_url"],
            destination_repo_name = x["destination_repo_name"],
            ) for x in values(git_lfs_repos)
        ]
    sort!(all_src_dest_pairs)
    unique!(all_src_dest_pairs)
    sort!(all_src_dest_pairs)

    @debug(
        string(
            "All source URLs in the list ",
            "($(length(all_src_dest_pairs))):",
            )
        )
    for i = 1:length(all_src_dest_pairs)
        @debug(
            "$(i). $(all_src_dest_pairs[i].source_url)"
            )
    end

    if task == "all"
        task_src_dest_pairs = all_src_dest_pairs
    elseif Types._is_interval(task)
        task_interval::Types.AbstractInterval = Types._construct_interval(
            task
            )
        @info(string("Using interval: "),task_interval,)
        task_src_dest_pairs = Common._pairs_that_fall_in_interval(
            all_src_dest_pairs,
            task_interval,
            )
    else
        @warn("not a valid task: ", task,)
        delayederror("\"$(task)\" is not a valid task")
    end

    unique!(task_src_dest_pairs)
    sort!(task_src_dest_pairs)
    @debug(
        string(
            "Source URLs in the task interval ",
            "($(length(task_src_dest_pairs))):",
            )
        )
    for i = 1:length(task_src_dest_pairs)
        @debug(
            "$(i). $(task_src_dest_pairs[i].source_url)"
            )
    end

    n = length(task_src_dest_pairs)
    for i = 1:n
        pair = task_src_dest_pairs[i]
        source_url = pair.source_url
        destination_repo_name = pair.destination_repo_name
        @debug(
            string(
                "Processing \"$(source_url)\" ",
                "(URL $(i) of $(n))",
                )
            )
        if false
        else
            if destination_repo_name in do_not_push_to_these_destinations
            else
                Common._snapshot_repo!!(
                    pair;
                    dst_provider = dst_provider,
                    include_branches = include_branches,
                    exclude_branches = exclude_branches,
                    git_user_name = git_user_name,
                    git_user_email = git_user_email,
                    time_zone = time_zone,
                    is_dry_run = is_dry_run,
                    )
            end
        end
    end

    @info(
        string(
            "SUCCESS: run_removelfs_snapshots completed ",
            "successfully :) Good-bye!",
            )
        )

    process_delayed_error_list()
    return nothing
end

end # End submodule RemoveLFS.Run

##### End of file


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module Odoo

using XMLRPC

export OdooServer, AuthenticatedOdooSession

immutable OdooServer
    url::AbstractString
end

const common_endpt = "/xmlrpc/2/common"
const object_endpt = "/xmlrpc/2/object"

immutable AuthenticatedOdooSession
    server::OdooServer
    db::AbstractString
    uid::Int
    pw::AbstractString
end

common(o::OdooServer) = XMLRPCProxy(o.url*common_endpt)
object(o::OdooServer) = XMLRPCProxy(o.url*object_endpt)
common(o::AuthenticatedOdooSession) = XMLRPCProxy(o.server.url*common_endpt)
object(o::AuthenticatedOdooSession) = XMLRPCProxy(o.server.url*object_endpt)

"""
Create a `AuthenticatedOdooSession` using the Odoo demo service.
"""
function demo_server()
    info = XMLRPCProxy("https://demo.odoo.com/start")["start"]()
    AuthenticatedOdooSession(info["host"], info["database"], info["user"], info["password"])
end

function version(o::Union{OdooServer, AuthenticatedOdooSession})
    common(o)["version"]()
end

function AuthenticatedOdooSession(o::OdooServer,
                                  db::AbstractString,
                                  username::AbstractString,
                                  pw::AbstractString)
    uid = common(o)["authenticate"](db, username, pw, [])
    AuthenticatedOdooSession(o, db, uid, pw)
end

function AuthenticatedOdooSession(o::AbstractString,
                                  db::AbstractString,
                                  username::AbstractString,
                                  pw::AbstractString)
    AuthenticatedOdooSession(OdooServer(o), db, username, pw)
end

function execute_kw(a::AuthenticatedOdooSession, args...)
    object(a)["execute_kw"](a.db, a.uid, a.pw, args...)
end

"""
List Odoo record fields, optionally constrain the
listed attributes.
"""
function Base.fieldnames(a::AuthenticatedOdooSession, d::AbstractString,
                         args::AbstractString...)
    execute_kw(a, d, "fields_get", [], isempty(args) ? Dict() :
                                       Dict("attributes" => [args...]))
end

"""
Search Odoo record fields.

Example:

```
search(odoo, "product.template", ["name", "=", "FooBarProduct"])
```
"""
function Base.search(a::AuthenticatedOdooSession, d::AbstractString,
                         args...)
    execute_kw(a, d, "search", Any[Any[args...]])
end

"""
Read Odoo record fields.

"""
function Base.read(a::AuthenticatedOdooSession, d::AbstractString,
                         args::Array)
    execute_kw(a, d, "read", Any[args])
end

Base.read(a::AuthenticatedOdooSession, d::AbstractString, args::Integer...) =
    read(a,d, [args...])

end # module

ext: jl


text:
# {{{ some steps before
# import Pkg
# Pkg.add("PyCall")
# }}}

# basic trick: use the PyCall module
using PyCall
pycdo = pyimport("cdo")
cdo = pycdo.Cdo()

print(cdo.topo())

numpy_data = cdo.topo(;returnArray="topo")
#print(numpy_data)

xarray_data = cdo.topo(;returnXArray="topo")
#print(xarray_data)

# chains
masked_data = cdo.setrtomiss(-100000,0;input="-expr,'logTopo=log(abs(topo)+0.1)' -topo,r10x10", returnArray="logTopo")
print(masked_data)

max_stars_repo_name: unartig/cdo-bindings

max_issues_repo_name: unartig/cdo-bindings

max_forks_repo_name: unartig/cdo-bindings


text:
## Define a particle consisting of all cluster variables for a Categorical mixture
mutable struct CategoricalCluster
  n::Int                          # count of cluster members
  counts::Matrix{Int}     # count of occurrence of levels
  nlevels::Array{Float64}           # The denominator in logprob calculation
  CategoricalCluster(dataFile) =
  new(0,
      zeros(Int, maximum(Int, dataFile), size(dataFile, 2)),
      # [0.5 * maximum(dataFile[:, d]) for d = 1:size(dataFile, 2)]
      0.5 * maximum(dataFile, dims = 1))
end

function copy_particle1(particle::CategoricalCluster, dataFile)
    new_particle = CategoricalCluster(dataFile)
    new_particle.n = particle.n
    for i in eachindex(particle.counts)
      new_particle.counts[i] = particle.counts[i]
    end
    for i in eachindex(particle.nlevels)
      new_particle.nlevels[i] = particle.nlevels[i]
    end
    return new_particle
end

function copy_particle(particle::CategoricalCluster, dataFile)
    return deepcopy(particle)
end

function calc_logprob(obs, cl::CategoricalCluster, featureFlag::Array)
  out = - sum(log.(cl.nlevels[featureFlag] .+ cl.n))
  @inbounds for q = 1:length(obs)
    if featureFlag[q]
      if cl.n == 0
        out += log(0.5)
      else
        out += log(0.5 + cl.counts[obs[q], q])
      end
    end
  end
  return out
end

function cluster_add!(cl::CategoricalCluster, obs, featureFlag::Array)
  @inbounds cl.n  += Int(1)
  @simd for q = 1:length(obs)
    if featureFlag[q]
      @inbounds cl.counts[obs[q], q] += 1
    end
  end
  return
end

function calc_logmarginal(cl::CategoricalCluster)
  # This returns the log of the marginal likelihood of the cluster
  # Used for feature selection
  # β_0 = 0.5
  lm = zeros(Float64, length(cl.nlevels))
  for q in 1:length(cl.nlevels)
    lm[q] += lgamma(cl.nlevels[q] * 2) - lgamma.(cl.nlevels[q] * 2 .+ cl.n)
    for r in 1:Int(2 * cl.nlevels[q])
      lm[q] += lgamma(cl.counts[r, q] + 0.5)

    end
  end
  return lm
end
using SpecialFunctions

"""
`coerce_categorical(data)`

Converts a matrix of categorical data into the format expected by particleMDI
## Input
- `data` an `n × d` matrix containing discrete values of which relatively few are unique.

## Output
- An `n × d` matrix where the discrete values are mapped to the range `1, 2, …, no. of unique values`. Ordering may not be preserved as ordering is irrelevant to this data type.
"""
function coerce_categorical(data)
  max_i = size(data, 1)
  max_j = size(data, 2)
  out = Matrix{Int}(undef, max_i, max_j)
  for j in 1:max_j
    u = unique(data[:, j])
    for i in 1:max_i
      out[i, j] = Int(findall(x -> x == data[i, j], u)[1])
    end
  end
  return out
end


text:
function tracer!(   i::Integer,
                    u::AbstractMatrix,
                    v::AbstractMatrix,
                    Prog::PrognosticVars,
                    Diag::DiagnosticVars,
                    S::ModelSetup)

    @unpack tracer_advection, tracer_consumption = S.parameters
    @unpack tracer_relaxation = S.parameters
    @unpack nadvstep_half,nadvstep = S.grid

    # mid point (in time) velocity for the advective time step
    if tracer_advection && ((i+nadvstep_half) % nadvstep) == 0
        @unpack um,vm = Diag.SemiLagrange

        copyto!(um,u)
        copyto!(vm,v)
    end

    if tracer_advection && (i % nadvstep) == 0

        @unpack sst = Prog
        @unpack ssti,sst_ref,dsst_comp = Diag.SemiLagrange

        # convert to type T for mixed precision
        sstrhs = convert(Diag.PrognosticVarsRHS.sst,sst)

        departure!(u,v,Diag,S)
        adv_sst!(sstrhs,Diag,S)

        if tracer_relaxation || tracer_consumption
            tracer_relax_consumption!(ssti,sst_ref,dsst_comp,S)
        end

        ghost_points_sst!(ssti,S)

        # copy back to sst, and convert to type Tprog if necessary
        copyto!(sst,ssti)
    end
    return nothing
end

"""Computes the departure point for semi-Lagrangian advection following Diamantakis, 2014.
u,v are assumed to be the time averaged velocities over the previous advection time step.
(Presumably need to be changed to 2nd order extrapolation in case the tracer is not passive)

Uses fixed-point iteration once to find the departure point."""
function departure!(u::AbstractMatrix,
                    v::AbstractMatrix,
                    Diag::DiagnosticVars,
                    S::ModelSetup)

    @unpack u_T,v_T,um,vm,um_T,vm_T = Diag.SemiLagrange
    @unpack uinterp,vinterp,xd,yd = Diag.SemiLagrange
    @unpack dtadvu,dtadvv,half_dtadvu,half_dtadvv = S.grid
    @unpack ep = S.grid

    # u,v is t + dtadv, um,vm are averaged over (t,t_dtadv)
    # interpolate u,um,v,vm onto the T-grid
    Ix!(u_T,u)
    Ix!(um_T,um)
    Iy!(v_T,v)
    Iy!(vm_T,vm)

    # simple mid-point in time advection
    #backtraj!(xd,dtadvu,um_T)
    #backtraj!(yd,dtadvv,vm_T)

    # fixed-point iteration for departure point
    # initial guess for departure point - mid point
    backtraj!(xd,half_dtadvu,u_T,ep)
    backtraj!(yd,half_dtadvv,v_T,ep)

    # interpolate um,vm onto mid-point
    interp_uv!(uinterp,um_T,xd,yd,S)
    interp_uv!(vinterp,vm_T,xd,yd,S)

    # update departure point
    backtraj!(xd,dtadvu,uinterp,ep)
    backtraj!(yd,dtadvv,vinterp,ep)
end

""" Solves the trajectory equation for a given time step dt and the velocity uv (this can be u or v).
One function for three cases

(i) u is interpolated from u-grid with halo onto T-grid
(ii) v is interpolated from v-grid with halo onto T-grid
(iii) u or v already on the T-grid: All matrices have same size.

Uses relative grid nodes in the departure points rd, such that actually solved is rd = 0 - dt*uv.
The indices i,j of rd determine the arrival grid point."""
function backtraj!( rd::Array{T,2},
                    dt::T,
                    uv::Array{T,2},
                    ep::Int) where {T<:AbstractFloat}
    m,n = size(rd)

    if (m,n) == size(uv)            # update departue point case, matrices have same size
        ishift = 0
        jshift = 0
    elseif (m+4,n+2) == size(uv)    # v-vel mid-point case, v has halo
        ishift = 2
        jshift = 1
    elseif (m+2+ep,n+4) == size(uv) # u-vel mid-point case, u has halo
        ishift = 1+ep
        jshift = 2
    else
        throw(BoundsError())
    end

    # relative grid means rd = 0 - dt*uv. The arrival information is stored in
    # the indices i,j of rd: rd[2,3] => -0.5,-0.5 means for the arrival grid node (2,3)
    # the departure is (2-0.5,3-0.5) = (1.5,2.5)
    for j ∈ 1:n
        for i ∈ 1:m
            rd[i,j] = -dt*uv[i+ishift,j+jshift]
        end
    end
end

""" Interpolates the matrix uv into the matrix uvi, where xx, yy specify the coordinates as indices (including fraction.
Interpolation only onto the inner entries of uvi. (They will be copied back later via the ghostpoint function).
Two cases
    (i) u velocities: from u-grid with halo to T-grid
    (ii) v velocities: from v-grid with halo to T-grid."""
function interp_uv!(uvi::AbstractMatrix,
                    uv::AbstractMatrix,
                    xx::AbstractMatrix,
                    yy::AbstractMatrix,
                    S::ModelSetup)

    @unpack ep = S.grid

    m,n = size(uvi)
    muv,nuv = size(uv)
    @boundscheck (m,n) == size(xx) || throw(BoundsError())
    @boundscheck (m,n) == size(yy) || throw(BoundsError())

    if (m+2+ep,n+4) == (muv,nuv)      # u case
        # halo sizes
        xh1 = 1+ep
        xh2 = 1
        yh1 = 2
        yh2 = 2
    elseif (m+4,n+2) == (muv,nuv)    # v case
        # halo sizes 1: left/bottom, 2: right/top
        xh1 = 2
        xh2 = 2
        yh1 = 1
        yh2 = 1
    else
        throw(BoundsError())
    end

    @unpack bc = S.parameters
    clip_or_wrap = if bc == "periodic" wrap else clip end

    for j ∈ 1:n
        for i ∈ 1:m
            xi = Int(floor(xx[i,j]))   # departure point indices lower left corner within grid cell
            yi = Int(floor(yy[i,j]))

            k,x0 = clip_or_wrap(xi,i,xx[i,j],muv,xh1,xh2)
            l,y0 = clip(yi,j,yy[i,j],nuv,yh1,yh2)

            uvi[i,j] = bilin(uv[k,l],uv[k+1,l],uv[k,l+1],uv[k+1,l+1],x0,y0)
        end
    end
end

"""Advection of sst/tracer based on the departure points xd,yd via bilinear interpolation.
Departure points are clipped/wrapped to remain within the domain. Boundary conditions either
periodic (wrap around behaviour) or no-flux (no gradient via clipping). Once the respective
4 surrounding grid points are found do bilinear interpolation on the unit square."""
function adv_sst!(  sst::AbstractMatrix,
                    Diag::DiagnosticVars,
                    S::ModelSetup)

    @unpack ssti,xd,yd = Diag.SemiLagrange
    @unpack halosstx,halossty = S.grid

    m,n = size(ssti)
    @boundscheck (m,n) == size(sst) || throw(BoundsError())
    @boundscheck (m-2*halosstx,n-2*halossty) == size(xd) || throw(BoundsError())
    @boundscheck (m-2*halosstx,n-2*halossty) == size(yd) || throw(BoundsError())

    @unpack bc = S.parameters
    clip_or_wrap = if bc == "periodic" wrap else clip end

    for j ∈ 1:n-2*halossty
        for i ∈ 1:m-2*halosstx

            xi = Int(floor(Float64(xd[i,j])))   # departure point lower left corner
            yi = Int(floor(Float64(yd[i,j])))   # coordinates

            k,x0 = clip_or_wrap(xi,i,xd[i,j],m,halosstx,halosstx)
            l,y0 = clip(yi,j,yd[i,j],n,halossty,halossty)

            ssti[i+halosstx,j+halossty] = bilin(sst[k,l],sst[k+1,l],sst[k,l+1],sst[k+1,l+1],x0,y0)
        end
    end
end

"""Clips the relative lower-left corner index xyi (for both x or y indices) to remain within the domain.
ij is the index of the underlying matrix. xy is the actual coordinate, mn (m or n) the size of the domain,
and h1,h2 are the halo sizes (left/south and right/north)."""
function clip(xyi::Int,ij::Int,xy::T,mn::Int,h1::Int,h2::Int) where {T<:AbstractFloat}

    xyis = xyi+ij+h1                # shifted index (i.e. revert the relative index to an absolute)

    if xyis < 1+h1                  # beyond left/southern boundary
        return 1+h1,zero(T)         # coordinate is then 0
    elseif xyis > mn-1              # beyond right/northern boundary
        return mn-1,one(T)          # coordinate is then 1
    else                            # normal case.
        return xyis,xy-T(xyi)       # relative coordinate ∈ [0,1]
    end
end

"""Clips the relative lower-left corner index xyi (for both x or y indices) to remain within the domain.
ij is the index of the underlying matrix. xy is the actual coordinate, mn (m or n) the size of the domain,
and h is the halo size."""
function wrap(xyi::Int,ij::Int,xy::T,mn::Int,h1::Int,h2::Int) where {T<:AbstractFloat}

    xyis = xyi+ij+h1                # shifted index (i.e. revert the relative index to an absolute)
    xy0 = xy-T(xyi)                 # relative coordinate ∈ [0,1]

    if xyis < 1+h1                  # beyond left/southern boundary
        return xyis+mn-h1-h2,xy0    # coordinate is wrapped around
    elseif xyis > mn-1              # beyond right/northern boundary
        return xyis-mn+h1+h2,xy0    # coordinate is wrapped around
    else                            # normal case.
        return xyis,xy0             # just pass
    end
end

"""Bilinear interpolation on (x,y) in the unit square [0,1]x[0,1].
The values at the corners are f00 = f(0,0), f01 = f(0,1), etc."""
function bilin(f00::T,f10::T,f01::T,f11::T,x::T,y::T) where {T<:AbstractFloat}
    # old version (more computations and larger rounding error)
    # oone = one(T)
    # return f00*(oone-x)*(oone-y) + f10*x*(oone-y) + f01*(oone-x)*y + f11*x*y

    # new version
    return x*(f10 - y*(f10-f11)) + (1-x)*(f00 - y*(f00 - f01))
end

"""Tracer relaxation."""
function tracer_relax_consumption!( sst::AbstractMatrix,
                                    sst_ref::AbstractMatrix,
                                    dsst_comp::AbstractMatrix,
                                    S::ModelSetup)
    @unpack τSST,jSST = S.constants
    @unpack halosstx,halossty = S.grid
    @unpack compensated = S.parameters

    τj = -(τSST+jSST)

    m,n = size(sst)
    @boundscheck size(sst) == size(sst_ref) || throw(BoundsError())

    if compensated
        for j ∈ 1+halossty:n-halossty
            for i ∈ 1+halosstx:m-halosstx
                dsst = (τj*sst[i,j] + τSST*sst_ref[i,j]) - dsst_comp[i,j]
                sst_new = sst[i,j] + dsst
                dsst_comp[i,j] = (sst_new - sst[i,j]) - dsst
                sst[i,j] = sst_new
            end
        end
    else
        for j ∈ 1+halossty:n-halossty
            for i ∈ 1+halosstx:m-halosstx
                sst[i,j] += τj*sst[i,j] + τSST*sst_ref[i,j]
            end
        end
    end
end


text:
using YaoBase, YaoArrayRegister
export Measure, AllLocs, ComputationalBasis

"""
    Measure{N, K, OT} <: PrimitiveBlock{N}
    Measure(n::Int; operator=ComputationalBasis(), locs=1:n, collapseto=nothing, remove=false)

Measure operator.
"""
mutable struct Measure{N, K, OT} <: PrimitiveBlock{N}
    operator::OT
    locations::Union{NTuple{K, Int}, AllLocs}
    collapseto::Union{Int, Nothing}
    remove::Bool
    results::Vector{Int}
    function Measure{N, K, OT}(operator, locations, collapseto, remove) where {N, K, OT}
        locations isa AllLocs || @assert_locs_safe N locations
        if collapseto !== nothing && remove == true
            error("invalid keyword combination, expect collapseto or remove, got (collapseto=$collapseto, remove=true)")
        end
        new{N, K, OT}(operator, locations, collapseto, remove)
    end
end

@interface nqubits_measured(::Measure{N, K}) where {N, K} = K

"""
    Measure(n::Int; operator=ComputationalBasis(), locs=AllLocs(), collapseto=nothing, remove=false)

Create a `Measure` block with number of qubits `n`.

# Example

You can create a `Measure` block on given basis (default is the computational basis).

```jldoctest
julia> Measure(4)
Measure(4)
```

Or you could specify which qubits you are going to measure

```jldoctest
julia> Measure(4; locs=1:3)
Measure(4;locs=(1, 2, 3))
```

by default this will collapse the current register to measure results.

```julia
julia> r = rand_state(3)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> state(r)
8×1 Array{Complex{Float64},2}:
  0.19864933724343375 - 0.4740335956912438im
  -0.2057912765333517 - 0.2262668486124923im
 -0.41680007712245676 - 0.13759187422609387im
 -0.24336704548326407 + 0.27343538360398184im
 -0.09308092255704317 - 0.005308959093704435im
  0.24555464152683212 + 0.02737969837364506im
  -0.3828287267256825 + 0.02401578941643196im
 0.048647936794205926 + 0.31047610497928607im

julia> r |> Measure(3)
Measure(3)

julia> state(r)
8×1 Array{Complex{Float64},2}:
                 0.0 + 0.0im
                 0.0 + 0.0im
 -0.9495962023170939 - 0.31347576069762273im
                 0.0 + 0.0im
                 0.0 + 0.0im
                 0.0 + 0.0im
                 0.0 + 0.0im
                 0.0 + 0.0im
```

But you can also specify the target bit configuration you want to collapse to with keyword `collapseto`.

```jldoctest
julia> m = Measure(4; collapseto=0b101)
Measure(4;collapseto=5)

julia> m.collapseto
5
```
"""
function Measure(n::Int; operator::OT=ComputationalBasis(), locs=AllLocs(), collapseto=nothing, remove=false) where OT
    if locs isa AllLocs
        Measure{n, n, OT}(operator, locs, collapseto, remove)
    else
        Measure{n, length(locs), OT}(operator, tuple(locs...), collapseto, remove)
    end
end

Measure(;locs=AllLocs(), operator=ComputationalBasis(), collapseto=nothing, remove=false) where K = @λ(n->Measure(n; locs=locs, operator=operator, collapseto=collapseto, remove=remove))
mat(x::Measure) = error("use BlockMap to get its matrix.")

function apply!(r::AbstractRegister, m::Measure{N}) where {N}
    _check_size(r, m)
    if m.collapseto !== nothing
        m.results = measure_collapseto!(m.operator, r, m.locations; config=m.collapseto)
    elseif m.remove
        m.results = measure_remove!(m.operator, r, m.locations)
    else
        m.results = measure!(m.operator, r, m.locations)
    end
    return m
end


text:
denSigns = [" ", "░", "▒", "▓", "█"]

type DensityCanvas <: Canvas
  grid::Array{(@compat UInt),2}
  colors::Array{(@compat UInt8),2}
  pixelWidth::Int
  pixelHeight::Int
  plotOriginX::(@compat AbstractFloat)
  plotOriginY::(@compat AbstractFloat)
  plotWidth::(@compat AbstractFloat)
  plotHeight::(@compat AbstractFloat)
  maxDensity::(@compat AbstractFloat)
end

function DensityCanvas(charWidth::Int, charHeight::Int;
                       plotOriginX::(@compat AbstractFloat) = 0.,
                       plotOriginY::(@compat AbstractFloat) = 0.,
                       plotWidth::(@compat AbstractFloat) = 1.,
                       plotHeight::(@compat AbstractFloat) = 1.)
  charWidth = max(charWidth, 5)
  charHeight = max(charHeight, 5)
  pixelWidth = charWidth
  pixelHeight = charHeight * 2
  plotWidth > 0 || throw(ArgumentError("Width has to be positive"))
  plotHeight > 0 || throw(ArgumentError("Height has to be positive"))
  grid = fill(0, charWidth, charHeight)
  colors = fill(0x00, charWidth, charHeight)
  DensityCanvas(grid, colors,
                pixelWidth, pixelHeight,
                plotOriginX, plotOriginY,
                plotWidth, plotHeight,
                1)
end

function printRow(io::IO, c::DensityCanvas, row::Int)
  nunrows = nrows(c)
  0 < row <= nunrows || throw(ArgumentError("Argument row out of bounds: $row"))
  y = row
  denSignCount = length(denSigns)
  valScale = (denSignCount - 1) / c.maxDensity
  for x in 1:ncols(c)
    #if c.grid[x,y] == 0
    #  printColor(c.colors[x,y], io, " ")
    #else
      denIndex = safeRound(c.grid[x,y] * valScale) + 1
      printColor(c.colors[x,y], io, denSigns[denIndex])
    #end
  end
end

function setPixel!(c::DensityCanvas, pixelX::Int, pixelY::Int, color::Symbol)
  0 <= pixelX <= c.pixelWidth || return nothing
  0 <= pixelY <= c.pixelHeight || return nothing
  pixelX = pixelX < c.pixelWidth ? pixelX: pixelX - 1
  pixelY = pixelY < c.pixelHeight ? pixelY: pixelY - 1
  cw, ch = size(c.grid)
  tmp = pixelX / c.pixelWidth * cw
  charX = safeFloor(tmp) + 1
  charY = safeFloor(pixelY / c.pixelHeight * ch) + 1
  c.grid[charX,charY] += 1
  c.maxDensity = max(c.maxDensity, c.grid[charX,charY])
  c.colors[charX,charY] = c.colors[charX,charY] | colorEncode[color]
  c
end


text:
using Crux
using Test
using Flux
using POMDPs
using POMDPGym
using Random
using Statistics
using Distributions
using CUDA

## Deep copy of gpu policies
π1 = ContinuousNetwork(Chain(Dense(5, 5, relu))) |> gpu
π1⁻ = deepcopy(π1)

polyak_average!(π1, π1⁻, 0.005f0)

weightbefore = cpu(π1⁻.network[1].weight)

π1.network[1].weight .= 1

@test cpu(π1⁻.network[1].weight)  == weightbefore

## Copyto and polyak average
c1 = Chain(Dense(5, 5, relu))
c2 = Chain(Dense(5, 5, relu))
c3 = Chain(Dense(5, 5, relu)) |> gpu

@test c1[1].weight != c2[1].weight
copyto!(c1, c2)

@test c1[1].weight == c2[1].weight

polyak_average!(c3, c2, 1f0)
@test c3[1].weight isa CuArray
@test cpu(c3[1].weight) == c2[1].weight

## Continuous Network
p = ContinuousNetwork(Chain(Dense(2, 32, relu), Dense(32, 4)))

@test Crux.device(p) == cpu
@test p.output_dim == 4

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test p_gpu.output_dim == 4

@test length(Flux.params(p)) == 4
@test length(layers(p)) == 2

s0 = rand(2)
s = rand(2,100)

@test size(value(p, s0)) == (4,)
@test size(value(p, s)) == (4,100)
@test size(value(p, s[1:1,:], s[2:2,:])) == (4, 100)
@test size(value(p, s0[1], s[2])) == (4,)

@test size(action(p, s0)) == (4,)
@test size(action(p, s)) == (4,100)
@test all(action(p, s) .≈ action(p_gpu, s))

@test value(p, s0) == action(p, s0)
@test value(p, s) == action(p, s)

@test action_space(p) == ContinuousSpace(4)

p2 = ContinuousNetwork(ConstantLayer(rand(10)), 10)
@test size(p2(rand(1))) == (10,)

## Discrete Network
p = DiscreteNetwork(Chain(Dense(2, 32, relu), Dense(32, 4)), [1,2, 3, 4])

@test Crux.device(p) == cpu
@test p.outputs == [1,2,3,4]

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test p_gpu.outputs == [1,2,3,4]

@test length(Flux.params(p)) == 4
@test length(layers(p)) == 2

s0 = rand(2)
s = rand(2,100)
a0 = rand([true, false], 4)
a = rand([true, false], 4, 100)

@test size(value(p, s0)) == (4,)
@test size(value(p, s)) == (4,100)

@test size(value(p, s0, a0)) == (1,)
@test size(value(p, s, a)) == (1, 100)

@test size(action(p, s0)) == (1,)
@test size(action(p, s)) == (1,100)
@test all(action(p, s) .≈ action(p_gpu, s))

a, logprob = exploration(p, s0)
@test size(a) == (1,)
@test size(logprob) == (1,)
@test all(logpdf(p, s0, Flux.onehotbatch(a, [1,2,3,4])) .≈ logprob)

a, logprob = exploration(p, s)
@test size(a) == (1,100)
@test size(logprob) == (1,100)
@test all(logpdf(p, s, Flux.onehotbatch(a[:], [1,2,3,4])) .≈ logprob)

@test size(entropy(p, s0)) == (1,)
@test size(entropy(p, s)) == (1, 100)

@test action_space(p).vals == DiscreteSpace(4).vals

# Onhotbatch - This used to error
Flux.onehotbatch(p_gpu, ones(Int, 1, 10) |> gpu)

# Testing the logit conversion function
p = DiscreteNetwork(Chain(Dense(2, 32, relu), Dense(32, 4, sigmoid)), [1,2, 3, 4], (s) -> s ./ sum(s,dims=1))

vals = value(p, s)
lgts = Crux.logits(p, s)
@test sum(lgts) ≈ 100

## Double Network
Q1 = ContinuousNetwork(Chain(Dense(2,32, relu), Dense(32, 1)))
Q2 = ContinuousNetwork(Chain(Dense(2,32, relu), Dense(32, 2)))
p = DoubleNetwork(Q1, Q2)

@test action_space(p) == action_space(p.N1)

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test Crux.device(p_gpu.N1) == gpu
@test Crux.device(p_gpu.N2) == gpu

@test length(Flux.params(p)) == 8
@test length(layers(p)) == 4

@test all(value(p, s) .≈ value(p_gpu, s))
@test value(p, s) == (value(p.N1, s), value(p.N2, s))
@test value(p, s[1:1,:], s[2:2,:]) == (value(p.N1, s[1:1,:], s[2:2,:]), value(p.N2, s[1:1,:], s[2:2,:]))
@test action(p, s) == (action(p.N1, s), action(p.N2, s))

Q1 = DiscreteNetwork(Chain(Dense(2,32, relu), Dense(32, 1)), [1,2])
Q2 = DiscreteNetwork(Chain(Dense(2,32, relu), Dense(32, 2)), [1,2])
p = DoubleNetwork(Q1, Q2)

Random.seed!(1)
e1 = exploration(p, s)
Random.seed!(1)
e2 = (exploration(p.N1, s), exploration(p.N2, s))
@test e1 == e2

a = rand([1,], 100)
@test logpdf(p, s, Flux.onehotbatch(Q1, a)) == (logpdf(p.N1, s, Flux.onehotbatch(Q1, a)), logpdf(p.N2, s, Flux.onehotbatch(Q2, a)))
@test entropy(p, s) == (entropy(p.N1, s), entropy(p.N2, s))

@test action_space(p) == action_space(p.N1)

## Actor Critic
Anet = ContinuousNetwork(Chain(Dense(2,32, relu), Dense(32, 4)))
C = ContinuousNetwork(Chain(Dense(6,32, relu), Dense(32, 1)))
p = ActorCritic(Anet, C)

@test action_space(p) == action_space(p.A)

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test Crux.device(p_gpu.A) == gpu
@test Crux.device(p_gpu.C) == gpu

@test length(Flux.params(p)) == 8
@test length(layers(p)) == 4

sfull = rand(6, 100)
a = rand(4, 100)
@test all(value(p, sfull) .≈ value(p_gpu, sfull))
@test value(p, sfull) == value(p.C, sfull)
@test value(p, s, a) == value(p.C, s, a)
@test action(p, s) == action(p.A, s)
@test all(action(p, s) .≈ action(p_gpu, s))

Anet = DiscreteNetwork(Chain(Dense(2,32, relu), Dense(32, 1)), [1,])
p = ActorCritic(Anet, C)

Random.seed!(1)
e1 = exploration(p, s)
Random.seed!(1)
e2 = exploration(p.A, s)
@test e1 == e2

a = rand([1,], 100)
@test logpdf(p, s, a) == logpdf(p.A, s, a)
@test entropy(p, s) == entropy(p.A, s)

@test action_space(p) == action_space(p.A)

## Gaussian Policy
μnet = ContinuousNetwork(Chain(Dense(2,32, relu), Dense(32, 4)))
logΣnet = ones(Float32, 4)
p = GaussianPolicy(μnet, logΣnet)
@test Crux.device(p) == cpu
@test all(value(p.logΣ, rand(2)) .== 1)

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test Crux.device(p_gpu.μ) == gpu
@test Crux.device(p_gpu.logΣ) == gpu
@test all(value(p.logΣ, rand(2)) .== 1)

@test length(Flux.params(p)) == 5
@test length(layers(p)) == 3

s0 = rand(2)
s = rand(2,100)
a0 = rand([1,2,3,4],1)
a = rand([1,2,3,4], 100)

@test size(action(p, s0)) == (4,)
@test size(action(p, s)) == (4,100)
@test all(action(p, s) .≈ action(p_gpu, s))

a, logprob = exploration(p, s0)
@test size(a) == (4,)
@test size(logprob) == (1,)
@test all(logpdf(p, s0, a) .≈ logprob)

a, logprob = exploration(p, s)
@test size(a) == (4,100)
@test size(logprob) == (1,100)
@test all(logpdf(p, s, a) .≈ logprob)

@test size(entropy(p, s0)) == ()
@test size(entropy(p, s)) == ()

@test action_space(p) == ContinuousSpace(4)

# Check empirically if the mean and std are correct
p =  GaussianPolicy(ContinuousNetwork((s)-> zeros(Float32, 6), 6), -0.5*ones(Float32, 6))
a = [exploration(p, [1])[1][1] for i=1:100000]

@test abs(std(a) .- exp(-0.5)) < 1e-2
@test abs(mean(a)) < 1e-2

## Squashed Gaussian Policy
trunk = Chain(Dense(2,32, relu))
μnet = ContinuousNetwork(Chain(trunk..., Dense(32, 4)))
logΣnet = ContinuousNetwork(Chain(trunk..., Dense(32, 4)))
p = SquashedGaussianPolicy(μnet, logΣnet)
@test Crux.device(p) == cpu

p_gpu = p |> gpu
@test Crux.device(p_gpu) == gpu
@test Crux.device(p_gpu.μ) == gpu
@test Crux.device(p_gpu.logΣ) == gpu

@test length(Flux.params(p)) == 6
@test length(layers(p))  == 3

s0 = rand(2)
s = rand(2,100)

@test size(action(p, s0)) == (4,)
@test size(action(p, s)) == (4,100)
@test all(action(p, s) .≈ action(p_gpu, s))

a, logprob = exploration(p, s0)
@test size(a) == (4,)
@test size(logprob) == (1,)
@test all(logpdf(p, s0, a) .≈ logprob)

a, logprob = exploration(p, s)
@test all(a .>= -1)
@test all(a .<= 1)
@test size(a) == (4,100)
@test size(logprob) == (1,100)
@test all(logpdf(p, s, a) .≈ logprob)

@test size(entropy(p, s0)) == (1,)
@test size(entropy(p, s)) == (1, 100)

@test action_space(p) == ContinuousSpace(4)

## Distribution Network
p = DistributionPolicy(Distributions.product_distribution([Uniform(-1,1) for i=1:6]))

@test Crux.device(p) == cpu

@test length(Flux.params(p)) == 0
@test length(layers(p)) == 0

s0 = rand(2)
s = rand(2,100)

a, logprob = exploration(p, s0)
@test size(a) == (6,)
@test size(logprob) == (1,1)
@test all(logpdf(p, s0, a) .≈ logprob)

a, logprob = exploration(p, s)
@test size(a) == (6,100)
@test size(logprob) == (1,100)
@test all(logpdf(p, s, a) .≈ logprob)

@test size(entropy(p, s0)) == (1,)
@test size(entropy(p, s)) == (1, 100)

@test action_space(p) == ContinuousSpace(6)

## ϵGreedyPolicy
p_on = DiscreteNetwork(Dense(2,2), [2,3])
p = ϵGreedyPolicy(0f0, [2,3])
a, logprob = exploration(p, s0; π_on=p_on, i=1)

@test a[1] in [2,3]
@test logprob[1] < 0

## GaussianNoiseExplorationPolicy
p_on = ContinuousNetwork((s)->[1.0], 1)
p = GaussianNoiseExplorationPolicy(10f0, ϵ_min=-0.1, ϵ_max=0.1)
a, logprob = exploration(p, s0; π_on=p_on, i=1)
@test any(a[1] .≈ [0.9, 1.1]) && isnan(logprob)
@test p.σ(1) == 10f0
@test p.a_min == [-Inf32]
@test p.a_max == [Inf32]
@test p.ϵ_min == -0.1f0
@test p.ϵ_max == 0.1f0

p = GaussianNoiseExplorationPolicy((i)->10f0, a_min=[-0.1], a_max=[0.1])
a, logprob = exploration(p, s0; π_on=p_on, i=1)
@test any(a[1] .≈ [-0.1, .1]) && isnan(logprob)
@test p.σ(1) == 10f0
@test p.ϵ_min == -Inf32
@test p.ϵ_max == Inf32
@test p.a_min == [-0.1f0]
@test p.a_max == [0.1f0]

## First exlore policy
p = FirstExplorePolicy(10, ContinuousNetwork((x)->[0,], 1))
a, logprob = exploration(p, s; π_on=ContinuousNetwork((x)->[1,], 1), i=2)
@test a==[0]

a, logprob = exploration(p, s; π_on=ContinuousNetwork((x)->[1,], 1), i=20)
@test a==[1]


text:
module Esolangs
export Boolfuck, Brainfuck, PortableMinskyMachineNotation

include("Boolfuck/src/Boolfuck.jl")
include("Brainfuck/src/Brainfuck.jl")
include("PortableMinskyMachineNotation/src/PortableMinskyMachineNotation.jl")

end # module

max_stars_repo_name: stellartux/Esolangs

max_issues_repo_name: stellartux/Esolangs

max_forks_repo_name: stellartux/Esolangs


text:
#Estimate the probability of a sample being in the interval [x1,x2] for a
#single-variable normal distribution
function normprob(x1::Real, x2::Real, mu::Real, sigma2::Real)
    sigma = sqrt(sigma2)
    0.5*(erf((x2-mu)/(sqrt(2)*sigma)) - erf((x1-mu)/(sqrt(2)*sigma)))
end

function adjustparam(x::Vector, dx_idx::Integer, dx::Real)
    dxv = zeros(length(x))
    dxv[abs(dx_idx)] = dx*sign(dx_idx)
    x + dxv
end

function isrepeat{D}(pt::Vector, model::GaussianProcess{D}, tol)
    for p = 1:length(model)
        if all(abs(getx(model,p) - pt) .<= tol*ones(D))
            return p
        end
    end
    return -1
end

#A common print helper
function _printvec(io::IO, x)
    print(io,"(")
    for i=1:length(x)
        @printf(io, "%6.3f", x[i])
        if i != length(x)
            print(io, ", ")
        end
    end
    print(io,")")
end


text:
export  initialPopPlot,
        initialPopAgePlot

function initialPopPlot( mpSim::MPsim; showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="" )

    savePlot |= filename != ""

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    report = initPopReport( mpSim )

    if isempty( report )
        return
    end  # if isempty( report )

    generateInitialPopPlot( report, showPlot, savePlot, filename )

end  # initialPopPlot( mpSim, showPlot, savePlot, filename )

function initialPopAgePlot( mpSim::MPsim, ageRes::Real, ageType::Symbol;
    showPlot::Bool=true, savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 )

    savePlot |= filename != ""

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    report = initPopAgeReport( mpSim, ageRes, ageType )

    if isempty( report )
        return
    end  # if isempty( report )

    generateInitialPopAgePlot( report, timeFactor, showPlot, savePlot,
        filename )

end  # initialPopAgePlot( mpSim, ageRes, ageType; showPlot, savePlot, filename,
     #   timeFactor )

include( "private/initpop.jl" )


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

using Test
using FiniteDifferences
using StaticArrays
using LinearAlgebra
using Suppressor
using Random
using Unitful
using Plots

using OpticSim
# curve/surface imports
using OpticSim: findspan, makemesh, knotstoinsert, coefficients, inside, quadraticroots, tobeziersegments, evalcsg, makiemesh
# interval imports
using OpticSim: α, halfspaceintersection, positivehalfspace, lower, upper, EmptyInterval, rayorigininterval, intervalcomplement, intervalintersection, intervalunion, RayOrigin, Infinity, Intersection
include("TestData/TestData.jl")
# bounding box imports
using OpticSim: doesintersect
# RBT imports
using OpticSim: rotmat, rotmatd
# lens imports
using OpticSim: reflectedray, snell, refractedray, trace, intersection, pathlength, power
#Bounding volume hierarchy imports
# using OpticSim: partition!

const COMP_TOLERANCE = 25 * eps(Float64)
const RTOLERANCE = 1e-10
const ATOLERANCE = 10 * eps(Float64)
const SEED = 12312487
const ALL_TESTS = isempty(ARGS) || "all" in ARGS || "All" in ARGS || "ALL" in ARGS

"""Optional testset: if the name of the testset is passed as an arg then it will be executed. If none are specified, or "all" is given as an arg then all optional testsets will be executed """
macro otestset(name, expr)
    quote
        if ALL_TESTS || $name in ARGS
            @testset $name begin
                $expr
            end
        end
    end
end

"""Evalute all functions not requiring arguments in a given module and test they don't throw anything"""
macro test_all_no_arg_functions(m)
    quote
        for n in names($m, all = true)
            # get all the valid functions
            if occursin("#", string(n)) || string(n) == string(nameof($m))
                continue
            end
            f = Core.eval($m, n)
            # get the methods of this function
            for meth in methods(f)
                # if this method has no args then try and evaluate it
                if meth.nargs == 1
                    # suppress STDOUT to prevent loads of stuff spamming the console
                    @test (@suppress_out begin
                        f()
                    end;
                    true)
                end
            end
        end
    end
end

include("testsets/JuliaLang.jl")
# include("testsets/BVH.jl")
include("testsets/TestData.jl")
# include("testsets/Examples.jl")  # slow
include("testsets/General.jl")
include("testsets/SurfaceDefs.jl")
include("testsets/Intersection.jl")
include("testsets/Lenses.jl")
include("testsets/Emitters.jl")  # TODO
include("testsets/Comparison.jl")
include("testsets/Visualization.jl")

include("Benchmarks/Benchmarks.jl")
include("testsets/Allocations.jl")

include("testsets/GlassCat.jl")


text:
import GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase
import Test

const MAX_STEPS = 3000
const NUM_RESETS = 3

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SingleRoomUndirectedModule.SingleRoomUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SingleRoomDirectedModule.SingleRoomDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.GridRoomsUndirectedModule.GridRoomsUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.GridRoomsDirectedModule.GridRoomsDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SequentialRoomsUndirectedModule.SequentialRoomsUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SequentialRoomsDirectedModule.SequentialRoomsDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.MazeUndirectedModule.MazeUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.MazeDirectedModule.MazeDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.GoToTargetUndirectedModule.GoToTargetUndirected} = (env.env.terminal_reward, env.env.terminal_penalty)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.GoToTargetDirectedModule.GoToTargetDirected} = (env.env.env.terminal_reward, env.env.env.terminal_penalty)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.DoorKeyUndirectedModule.DoorKeyUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.DoorKeyDirectedModule.DoorKeyDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.CollectGemsUndirectedModule.CollectGemsUndirected} = (env.env.gem_reward * env.env.num_gem_init,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.CollectGemsDirectedModule.CollectGemsDirected} = (env.env.env.gem_reward * env.env.env.num_gem_init,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.CollectGemsMultiAgentUndirectedModule.CollectGemsMultiAgentUndirected} = (env.env.num_gem_init * env.env.gem_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.DynamicObstaclesUndirectedModule.DynamicObstaclesUndirected} = (env.env.terminal_reward, env.env.terminal_penalty)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.DynamicObstaclesDirectedModule.DynamicObstaclesDirected} = (env.env.env.terminal_reward, env.env.env.terminal_penalty)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SokobanUndirectedModule.SokobanUndirected} = (typeof(env.env.reward)(length(env.env.box_positions)),)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.SokobanDirectedModule.SokobanDirected} = (env.env.env.terminal_reward,)

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.CatcherModule.Catcher} = env.env.terminal_penalty : env.env.gem_reward : MAX_STEPS * env.env.gem_reward

get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.TransportUndirectedModule.TransportUndirected} = (env.env.terminal_reward,)
get_terminal_returns(env::GW.RLBaseEnv{E}) where {E <: GW.TransportDirectedModule.TransportDirected} = (env.env.env.terminal_reward,)

function is_valid_terminal_return(env::GW.RLBaseEnv{E}, terminal_return) where {E <: GW.SnakeModule.Snake}
    terminal_reward = env.env.terminal_reward
    terminal_penalty = env.env.terminal_penalty
    food_reward = env.env.food_reward
    _, height, width = size(env.env.tile_map)

    terminal_returns_win = terminal_reward : food_reward : terminal_reward + height * width * food_reward
    terminal_returns_lose = terminal_penalty : food_reward : terminal_penalty + height * width * food_reward

    if terminal_return in terminal_returns_win
        return true
    elseif terminal_return in terminal_returns_lose
        return true
    else
        return false
    end
end

Test.@testset "GridWorlds.jl" begin
    for Env in GW.ENVS
        Test.@testset "$(Env)" begin
            R = Float32
            env = GW.RLBaseEnv(Env(R = R))
            for _ in 1:NUM_RESETS
                RLBase.reset!(env)
                Test.@test RLBase.reward(env) == zero(R)
                Test.@test RLBase.is_terminated(env) == false

                total_reward = zero(R)
                for i in 1:MAX_STEPS
                    state = RLBase.state(env)
                    action = rand(RLBase.action_space(env))
                    env(action)
                    total_reward += RLBase.reward(env)

                    if RLBase.is_terminated(env)
                        if Env == GW.SnakeModule.Snake
                            is_valid_terminal_return(env, total_reward)
                        else
                            Test.@test total_reward in get_terminal_returns(env)
                        end

                        break
                    end

                    if i == MAX_STEPS
                        @info "$Env not terminated after MAX_STEPS = $MAX_STEPS"
                    end
                end
            end
        end
    end
end


text:
@testset "context" begin
    notebook = PlutoStaticHTML._load_notebook(NOTEBOOK_PATH)
    html = notebook2html!(notebook; append_build_context=true)
    @test contains(html, r"Built with Julia 1.*")
    @test contains(html, "CairoMakie")
end


text:
export state, register_spk!

import Base.convert

abstract type AbstractSPKType end

# Type 1: Modified Difference arrays
#FIXME
struct SPKMDA <: AbstractSPKType
    coeffs::Matrix{Float64}
    epochs::Vector{Float64}
end

function convert(::Type{SPKMDA}, data::Vector{Float64})
    N = Int(data[end])
    diff_line_coeffs = hcat([data[i:(i + 70)] for i in 1:71:(71N)]...)
    epochs = [data[i] for i in (71(N + 1)):((71(N + 1)) + N)]
    return SPKMDA(diff_line_coeffs, epochs)
end

function (r::SPKMDA)(x::Real)
    return pos, vel
end

function state(entry::SPKMDA, time::Real)
    # Find the coeff line that is in range and evaluate
end

# Type 2: Chebyshev (position only)
struct ChebyshevRecord
    midpoint::Float64
    radius::Float64
    x::ChebyshevT
    y::ChebyshevT
    z::ChebyshevT
    ẋ::ChebyshevT
    ẏ::ChebyshevT
    ż::ChebyshevT
end

function (r::ChebyshevRecord)(x::Real)
    s = (x - r.midpoint) / r.radius
    pos = [r.x(s), r.y(s), r.z(s)]
    vel = [r.ẋ(s), r.ẏ(s), r.ż(s)]
    return pos, vel
end
struct SPKChebyshev <: AbstractSPKType
    records::Vector{ChebyshevRecord}
    INIT::Float64
    INTLEN::Float64
    RSIZE::Float64
    N::Float64
end

function convert(::Type{SPKChebyshev}, data::Vector{Float64})
    init = Int(data[end - 3])
    intlen = Int(data[end - 2])
    rsize = Int(data[end - 1])
    n = Int(data[end])
    degree = (rsize - 2) ÷ 3 - 1
    records = ChebyshevRecord[]
    @inbounds for record in 1:n
        offset = (record - 1) * rsize + 1
        mid = data[offset]
        radius = data[offset + 1]
        # Grab coefficients
        x_coefs = data[(offset + 2):(offset + 2 + degree)]
        y_coefs = data[(offset + 3 + degree):(offset + 3 + 2degree)]
        z_coefs = data[(offset + 4 + 2degree):(offset + 4 + 3degree)]
        # Build polynomials
        x = ChebyshevT(x_coefs)
        y = ChebyshevT(y_coefs)
        z = ChebyshevT(z_coefs)
        ẋ = derivative(x)
        ẏ = derivative(y)
        ż = derivative(x)
        # And create
        push!(records, ChebyshevRecord(mid, radius, x, y, z, ẋ, ẏ, ż))
    end
    # Sort
    sort!(records; by=x -> x.midpoint)
    return SPKChebyshev(records, init, intlen, rsize, n)
end

function state(entry::SPKChebyshev, time::Real)
    # Verify the target time is even in this range
    ending = entry.INIT + entry.INTLEN * entry.N
    @assert entry.INIT <= time <= ending "Target time not in range"
    # Binary search through segments to find the polynomial we care about
    found = false
    search_range = 1:Int(entry.N)
    current = -1
    while !found
        current = (search_range[end] - search_range[1]) ÷ 2 + search_range[1]
        # Check if time is covered by this segment
        current_record = entry.records[current]
        rstart = current_record.midpoint - current_record.radius
        rend = current_record.midpoint + current_record.radius
        if rstart <= time <= rend
            found = true
            break
        elseif length(search_range) == 1
            break
        elseif time < rstart
            search_range = search_range[1]:(current - 1)
        elseif time > rend
            search_range = (current + 1):search_range[end]
        end
    end
    @assert found "Segments don't cover requested time"
    # Evaluate time on segment
    return entry.records[current](time)
end

# There are more types, do we need them?
const SPK_TYPES = Dict(1 => SPKMDA, 2 => SPKChebyshev)
const SPK_REGISTRY = Dict{NTuple{3,Int},AbstractSPKType}()

"""
    register_spk!(daf)

Registers the ephemeris data from a DAF file into the global registry
"""
function register_spk!(daf::DAFFile)
    # Reinterpret DAF into an SPK file
    @assert daf.id == "DAF/SPK " "Input DAF is not a valid SPK file"
    # Each array in the DAF file represents a segment for a given descriptor
    # We are going to represent it as (target,center,reference) => callable
    for segment in daf.arrays
        target = segment.ints[1]
        center = segment.ints[2]
        reference = segment.ints[3]
        @assert segment.ints[4] ∈ keys(SPK_TYPES) "Unknown SPK representation"
        type = SPK_TYPES[segment.ints[4]]
        callable = convert(type, segment.data)
        SPK_REGISTRY[(target, center, reference)] = callable
    end
end

function register_spk!(filename::String)
    return register_spk!(read_daf(filename))
end

"""
    state(time,target,center)

Evaluates the state vectors (position and velocity) at a given time for a target and center system.

# Arguments
- `time::Real`: Time in TDB seconds since J2000
- `target::Int`: NAIF code for target
- `center::Int`: NAIF code for center

# Optional Arguments
- `reference::Int`: NAIF code for reference system
"""
function state(time::Real, target::Int, center::Int, reference::Int=1)
    # Check to see if a record is registered with this combination
    key = (target, center, reference)
    @assert key ∈ keys(SPK_REGISTRY) "No SPK entry found for this target/center/reference combination"
    # Evaluate state for the one that matches
    return state(SPK_REGISTRY[key], time)
end

function state(time::Real, target::String, center::String, reference::String="J2000")
    return state(time, naif_body_codes[target], naif_body_codes[center],
                 naif_reference_codes[reference])
end


text:
#simple DFT function
function DFT(x)
    N = length(x)

    # We want two vectors here for real space (n) and frequency space (k)
    n = 0:N-1
    k = n'
    transform_matrix = exp.(-2im*pi*n*k/N)
    return transform_matrix*x

end

# Implementing the Cooley-Tukey Algorithm
function cooley_tukey(x)
    N = length(x)

    if (N > 2)
        x_odd = cooley_tukey(x[1:2:N])
        x_even = cooley_tukey(x[2:2:N])
    else
        x_odd = x[1]
        x_even = x[2]
    end
    n = 0:N-1
    half = div(N,2)
    factor = exp.(-2im*pi*n/N)
    return vcat(x_odd + x_even .* factor[1:half],
                x_odd - x_even .* factor[1:half])

end

function bitreverse(a::Array)
    # First, we need to find the necessary number of bits
    digits = convert(Int,ceil(log2(length(a))))

    indices = [i for i = 0:length(a)-1]

    bit_indices = []
    for i = 1:length(indices)
        push!(bit_indices, bits(indices[i]))
    end

    # Now stripping the unnecessary numbers
    for i = 1:length(bit_indices)
        bit_indices[i] = bit_indices[i][end-digits:end]
    end

    # Flipping the bits
    for i =1:length(bit_indices)
        bit_indices[i] = reverse(bit_indices[i])
    end

    #replacing indices
    for i = 1:length(indices)
        indices[i] = 0
        for j = 1:digits
            indices[i] += 2^(j-1) * parse(string(bit_indices[i][end-j]))
        end
       indices[i] += 1
    end

    b = [float(i) for i = 1:length(a)]
    for i = 1:length(indices)
        b[i] = a[indices[i]]
    end

    return b
end

function iterative_cooley_tukey(x)
    N = length(x)
    logN = convert(Int,ceil(log2(length(x))))
    bnum = div(N,2)
    stride = 0;

    x = bitreverse(x)

    z = [Complex(x[i]) for i = 1:length(x)]
    for i = 1:logN
       stride = div(N, bnum)
       for j = 0:bnum-1
           start_index = j*stride + 1
           y = butterfly(z[start_index:start_index + stride - 1])
           for k = 1:length(y)
               z[start_index+k-1] = y[k]
           end
       end
       bnum = div(bnum,2)
    end

    return z
end

function butterfly(x)
    N = length(x)
    half = div(N,2)
    n = [i for i = 0:N-1]
    half = div(N,2)
    factor = exp.(-2im*pi*n/N)

    y = [0 + 0.0im for i = 1:length(x)]

    for i = 1:half
        y[i] = x[i] + x[half+i]*factor[i]
        y[half+i] = x[i] - x[half+i]*factor[i]
    end

    return y
end

function approx(x, y)
    val = true
    for i = 1:length(x)
        if (abs(x[i]) - abs(y[i]) > 1e-5)
            val = false
        end
    end
    println(val)
end

function main()
    x = rand(128)
    y = cooley_tukey(x)
    z = iterative_cooley_tukey(x)
    w = fft(x)
    approx(y, w)
    approx(z, w)
end

main()

max_stars_repo_name: sanstorik/algorithm-archive

max_issues_repo_name: sanstorik/algorithm-archive

max_forks_repo_name: sanstorik/algorithm-archive


text:
"""
    Geometric(p)

A *Geometric distribution* characterizes the number of failures before the first success in a sequence of independent Bernoulli trials with success rate `p`.

```math
P(X = k) = p (1 - p)^k, \\quad \\text{for } k = 0, 1, 2, \\ldots.
```

```julia
Geometric()    # Geometric distribution with success rate 0.5
Geometric(p)   # Geometric distribution with success rate p

params(d)      # Get the parameters, i.e. (p,)
succprob(d)    # Get the success rate, i.e. p
failprob(d)    # Get the failure rate, i.e. 1 - p
```

External links

*  [Geometric distribution on Wikipedia](http://en.wikipedia.org/wiki/Geometric_distribution)

"""
struct Geometric{T<:Real} <: DiscreteUnivariateDistribution
    p::T

    function Geometric{T}(p::T) where {T <: Real}
        new{T}(p)
    end
end

function Geometric(p::Real; check_args::Bool=true)
    @check_args Geometric (p, zero(p) < p < one(p))
    return Geometric{typeof(p)}(p)
end

Geometric() = Geometric{Float64}(0.5)

@distr_support Geometric 0 Inf

### Conversions
convert(::Type{Geometric{T}}, p::Real) where {T<:Real} = Geometric(T(p))
Base.convert(::Type{Geometric{T}}, d::Geometric) where {T<:Real} = Geometric{T}(T(d.p))
Base.convert(::Type{Geometric{T}}, d::Geometric{T}) where {T<:Real} = d

### Parameters

succprob(d::Geometric) = d.p
failprob(d::Geometric) = 1 - d.p
params(d::Geometric) = (d.p,)
partype(::Geometric{T}) where {T<:Real} = T

### Statistics

mean(d::Geometric) = failprob(d) / succprob(d)

median(d::Geometric) = -fld(logtwo, log1p(-d.p)) - 1

mode(d::Geometric{T}) where {T<:Real} = zero(T)

var(d::Geometric) = (1 - d.p) / abs2(d.p)

skewness(d::Geometric) = (2 - d.p) / sqrt(1 - d.p)

kurtosis(d::Geometric) = 6 + abs2(d.p) / (1 - d.p)

entropy(d::Geometric) = (-xlogx(succprob(d)) - xlogx(failprob(d))) / d.p

### Evaluations

function logpdf(d::Geometric, x::Real)
    insupport(d, x) ? log(d.p) + log1p(-d.p) * x : log(zero(d.p))
end

function cdf(d::Geometric, x::Int)
    p = succprob(d)
    n = max(x + 1, 0)
    p < 1/2 ? -expm1(log1p(-p)*n) : 1 - (1 - p)^n
end

ccdf(d::Geometric, x::Real) = ccdf_int(d, x)
function ccdf(d::Geometric, x::Int)
    p = succprob(d)
    n = max(x + 1, 0)
    p < 1/2 ? exp(log1p(-p)*n) : (1 - p)^n
end

logcdf(d::Geometric, x::Real) = logcdf_int(d, x)
logcdf(d::Geometric, x::Int) = log1mexp(log1p(-d.p) * max(x + 1, 0))

logccdf(d::Geometric, x::Real) = logccdf_int(d, x)
logccdf(d::Geometric, x::Int) =  log1p(-d.p) * max(x + 1, 0)

quantile(d::Geometric, p::Real) = invlogccdf(d, log1p(-p))

cquantile(d::Geometric, p::Real) = invlogccdf(d, log(p))

invlogcdf(d::Geometric, lp::Real) = invlogccdf(d, log1mexp(lp))

function invlogccdf(d::Geometric{T}, lp::Real) where T<:Real
    if (lp > zero(d.p)) || isnan(lp)
        return T(NaN)
    elseif isinf(lp)
        return T(Inf)
    elseif lp == zero(d.p)
        return zero(T)
    end
    max(ceil(lp/log1p(-d.p)) - 1, zero(T))
end

function mgf(d::Geometric, t::Real)
    p = succprob(d)
    p / (expm1(-t) + p)
end

function cf(d::Geometric, t::Real)
    p = succprob(d)
    # replace with expm1 when complex version available
    p / (exp(-t*im) - 1 + p)
end

### Sampling

rand(rng::AbstractRNG, d::Geometric) = floor(Int,-randexp(rng) / log1p(-d.p))

### Model Fitting

struct GeometricStats <: SufficientStats
    sx::Float64
    tw::Float64

    GeometricStats(sx::Real, tw::Real) = new(sx, tw)
end

suffstats(::Type{<:Geometric}, x::AbstractArray{T}) where {T<:Integer} = GeometricStats(sum(x), length(x))

function suffstats(::Type{<:Geometric}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    n = length(x)
    if length(w) != n
        throw(DimensionMismatch("Inconsistent argument dimensions."))
    end
    sx = 0.
    tw = 0.
    for i = 1:n
        wi = w[i]
        sx += wi * x[i]
        tw += wi
    end
    GeometricStats(sx, tw)
end

fit_mle(::Type{<:Geometric}, ss::GeometricStats) = Geometric(1 / (ss.sx / ss.tw + 1))


text:
# Create example data
data = randn(100)
# Draw a polar histogram with 19 bins
polarhistogram(data, nbins = 19, fullcircle=true)


text:
module DiffusionDefinition

    using Random, Trajectories
    using LinearAlgebra, StaticArrays
    using MacroTools

    const ℝ{N} = SVector{N,Float64}

    include("types.jl")
    include("state_space_restrictions.jl")
    include("diffusion_process.jl")
    include(joinpath("..", "examples", "example_list.jl"))
    include("buffer.jl")
    include("sampling.jl")

    export @diffusion_process
    export @load_diffusion

    export UnboundedStateSpace, LowerBoundedStateSpace, UpperBoundedStateSpace,
        BoundedStateSpace

    export AbstractBuffer
    export Wiener

    using Reexport
    @reexport using Trajectories
end # module


text:
struct JSON <: Format end

formatextension(::JSON) = "json"

dict(x::XNode) = Dict(
    :tag => tag(x),
    :attributes => attributes(x),
    :children => dict.(children(x))
)

# TODO: embed mime types
function dict(x::XLeaf{<:AbstractString})
    return Dict(:mimes => Dict("text/plain" => x[]))

end

function dict(x::XLeaf{Nothing})
    Dict(:mimes => Dict("text/plain" => ""))
end
function dict(x::XLeaf)
    mimes = IJulia.display_dict(x[])
    d = Dict(:mimes => mimes)
end

function render!(io, x::XTree, format::JSON, val)
    JSON3.write(io, dict(x))
end

dict(x::XLeaf{PreRendered}) = Dict(:mimes => x[].reprs)


text:
using Printf, DataStructures

function funcfreqs(expr::Expr)
    cnt = counter(Symbol)
    expr.head == :call &&
        push!(cnt, expr.args[1])
    for e in expr.args
        e isa Expr && merge!(cnt, funcfreqs(e))
    end
    return cnt
end

function parseall(str::AbstractString)
    exs = Any[]
    pos = start(str)
    while !done(str, pos)
        ex, pos = parse(str, pos) # returns next starting point as well as expr
        ex.head == :toplevel ? append!(exs, ex.args) : push!(exs, ex)
    end
    if isempty(exs)
        throw(ParseError("end of input"))
    elseif length(exs) == 1
        return exs[1]
    else
        return Expr(:block, exs...)
    end
end

freqs = readstring("src/Function_frequency.jl") |> parseall |> funcfreqs

for (v, f) in freqs
    @printf("%10s → %i\n", v, f)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ## Test note for RunMyNotes.jl
# This is a plain text `.jl` file formatted so that
# [Literate.jl](https://github.com/fredrikekre/Literate.jl) will understand it.

[i^2 for i=1:10]  # column
[i^2 for i=1:10]' # row

#-

using Plots
## notice that Plots is mentioned in /test/REQUIRE
plot(rand(20))

# Commented text becomes text (markdown) cells, `##` stays a comment,
# and `#-` divides two code cells to keep their output visible.


text:
mutable struct SimulationSamples
  input::InputParameters
  #output::SimResult
  sampleddataDF::DataFrame
end

function run1simulationsample(IP::InputParameters, minclonesize, maxclonesize)

  sresult = tumourgrow_birthdeathsample(IP.b, IP.d, IP.Nmax, IP.μ; numclones = IP.numclones, s = IP.selection, tevent = IP.tevent, clonalmutations = IP.clonalmutations);

  return sresult, IP
end

function getsamples(;nsample = 1, nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^3, samplesize = round(Int64, 0.1 * Nmax), s = repeat([1.0], inner = nclones), tevent = collect(1.0:0.5:100.0)[1:nclones], cellularity = 1.0, fixedmu = false)

  nclones == length(s) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(s)), these must be the same size ")

  nclones == length(tevent) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(tevent)), these must be the same size ")

  IP = InputParameters(nclones,
  Nmax,
  detectionlimit,
  ploidy,
  read_depth,
  clonalmutations,
  s,
  μ,
  b,
  d,
  tevent,
  ρ,
  cellularity,
  fixedmu)

  simresult, IP = run1simulationsample(IP, 0.0, 1.0)

  samples = sample(simresult.cells, samplesize * nsample, replace = false)

  samplesvec = []
  for i in 0:(nsample - 1)
      push!(samplesvec, samples[(i * samplesize + 1):(samplesize * (i+ 1))])
  end

  mutations = unique(CancerSeqSim.cellsconvert(samples)[1]);
  DF = DataFrame(mutations = 1:maximum(mutations))

  for i in 1:nsample
      DF[Symbol("sample$i")] = counts(CancerSeqSim.cellsconvert(samplesvec[i])[1], 1:maximum(mutations))
  end

  DFs = sampledhistmultiple(Array(DF[:sample1]), Array(DF[:mutations]), samplesize, 1)
  for i in 2:nsample
      DFtemp = sampledhistmultiple(Array(DF[Symbol("sample$i")]), Array(DF[:mutations]), samplesize, i, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
      DFs = join(DFs, DFtemp, on = :mutations)
  end

  return SimulationSamples(IP, DFs), simresult, DF

end

function sampledhistmultiple(AF, mutations, samplesize, nsample; detectionlimit = 0.1, ploidy = 2.0, read_depth = 100.0, cellularity = 1.0)

    AF = AF./ploidy

    AF = AF .* cellularity

    ind = AF.>(detectionlimit * samplesize)
    AF = AF[ind]
    mutations = mutations[ind]

    samp_percent = read_depth/samplesize

    depth = rand(Binomial(samplesize,samp_percent), length(AF))

    samp_alleles = map((n, p) -> rand(Binomial(n, p)), depth, AF/samplesize)

    VAF = samp_alleles./depth

    DF = DataFrame(mutations = mutations)
    DF[Symbol("VAF$(nsample)")] = VAF

    return DF
end


text:
module JSON2

using Dates, Parsers

# for reading/writing javascript functions as value literals
struct Function
    str::String
end

include("write.jl")
include("read.jl")
include("strings.jl")
include("pretty.jl")

## JSON2.@format
function getformats(nm; kwargs...)
    for (k, v) in pairs(kwargs)
        k == nm && return v
    end
    return NamedTuple()
end

macro format(T, expr)
    if occursin("keywordargs", string(expr)) || occursin("noargs", string(expr))
        esc(:(JSON2.@format($T, $expr, $(Expr(:block)))))
    else
        esc(:(JSON2.@format($T, "", $expr)))
    end
end

"""
    JSON2.@format T [noargs|keywordargs] begin
        _field_ => (; options...)
        _field2_ => (; options...)
    end

Specify a custom JSON formatting for a struct `T`. Options include:
    `name`:
    `jsontype`:
    `omitempty`:
    `exclude`:
    `default`:

By default, the JSON input is expected to contain each field of a type and be in the same order as the type was defined. For example, the struct:
```julia
struct T
    a::Int
    b::Int
    c::Union{Nothing, Int}
end
```
Would have JSON like:
```
{"a": 0, "b": 1, "c": null}
{"a": 0, "b": 1, "c": 2}
{"a": 0, "b": 1, "c": null, "d": 3} // extra fields are ignored
{"a": 0} // will work if T(a) constructor is defined
{"a": 0, "b": 1} // will work if T(a, b) constructor is defined
```
That is, each field _must_ be present in the JSON input and match in position to the original struct definition. Extra arguments after the struct's own fieldtypes are ignored.

Again, the default case is for JSON input that will have consistently ordered, always-present fields; for cases where the input JSON is _not_ well-ordered or if there is a possibility of a `null`, non-`Union{T, Nothing}` field, here's how to approach more advanced custom formatting:
    - If the input will always be consistenly-ordered, but fields may be missing (not `null`, but the field key isn't even available in the input), defaults can be provided like:
    ```
    JSON2.@format T begin
        c => (default=0,)
    end
    ```
    This says that, when reading from a JSON input, if field `c` isn't present, to set it's value to 0.

    - If the JSON input is not consistenly-ordered, there are two other options for allowing direct type parsing:

      ```
      T(; a=0, b=0, c=0, kwargs...) = T(a, b, c)
      JSON2.@format T keywordargs begin
        ...
      end
      ```
      Here we've defined a "keywordargs" constructor for `T` that essentially takes a default for each field as keyword arguments, then constructs `T`.
      During parsing, the JSON input will be parsed for any valid field key-values and the keyword constructor will be called
      with whatever arguments are parsed in whatever order. Note that we also included a catchall `kwargs...` in our constructor which can be used to "throw away" or ignore any extra fields in the JSON input.

      ```
      mutable struct T
          a::Int
          b::Int
          c::Union{Nothing, Int}
      end
      T() = T(0, 0, 0)
      JSON2.@format T noargs begin
        ...
      end
      ```
      In this case, we've made `T` a _mutable_ struct and defined a "noargs" constructor `T() = ...`; we then specified in `JSON2.@format T noargs` the `noargs` option.
      During parsing, an instance of `T` will first constructed using the "noargs" constructor, then fields will be set as they're parsed from the JSON input (hence why `mutable struct` is required).

"""
macro format(T, typetype, expr)
    args = filter(x->typeof(x) != LineNumberNode, expr.args)
    foreach(x->x.args[2] = QuoteNode(x.args[2]), args)
    anydefaults = any(z->:default in z, map(x->map(y->y.args[1], x.args[3].args), args))
    wr = quote
        @generated function JSON2.write(io::IO, obj::$T)
            fieldformats = [JSON2.getformats(nm; $(args...)) for nm in fieldnames($T)]
            # @show fieldformats
            inds = Tuple(i for i = 1:length(fieldformats) if !get(fieldformats[i], :exclude, false))
            names = Tuple(string(get(fieldformats[i], :name, fieldname($T, i))) for i in inds)
            omitempties = Tuple(get(fieldformats[i], :omitempty, false) for i in inds)
            converts = Tuple(JSON2.getconvert(get(fieldformats[i], :jsontype, fieldtype($T, i))) for i in inds)
            N = length(inds)
            ex = JSON2.generate_write_body(N, inds, names, omitempties, converts)
            # @show ex
            return ex
        end
    end
    if occursin("noargs", string(typetype))
        q = esc(quote
            @generated function JSON2.read(io::IO, T::Type{$T})
                N = fieldcount($T)
                fieldformats = Dict($(args...))
                names = (; ((Symbol(get(get(fieldformats, nm, NamedTuple()), :name, nm)), nm) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false))...)
                jsontypes = (; ((get(get(fieldformats, nm, NamedTuple()), :name, nm), get(get(fieldformats, nm, NamedTuple()), :jsontype, fieldtype($T, i))) for (i, nm) in enumerate(fieldnames($T)) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false))...)
                defaults = (; ((get(get(fieldformats, nm, NamedTuple()), :name, nm), get(fieldformats, nm, NamedTuple())[:default]) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false) && haskey(get(fieldformats, nm, NamedTuple()), :default))...)
                return JSON2.generate_read_body_noargs(N, names, jsontypes, defaults)
            end
            $wr
        end)
    elseif occursin("keywordargs", string(typetype))
        q = esc(quote
            @generated function JSON2.read(io::IO, T::Type{$T})
                N = fieldcount($T)
                fieldformats = Dict($(args...))
                names = (; ((Symbol(get(get(fieldformats, nm, NamedTuple()), :name, nm)), nm) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false) && haskey(get(fieldformats, nm, NamedTuple()), :name))...)
                types = (; ((get(get(fieldformats, nm, NamedTuple()), :name, nm), get(get(fieldformats, nm, NamedTuple()), :jsontype, fieldtype($T, i))) for (i, nm) in enumerate(fieldnames($T)) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false))...)
                defaults = (; ((get(get(fieldformats, nm, NamedTuple()), :name, nm), get(fieldformats, nm, NamedTuple())[:default]) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false) && haskey(get(fieldformats, nm, NamedTuple()), :default))...)
                q = quote
                    JSON2.@expect '{'
                    JSON2.wh!(io)
                    keys = Symbol[]
                    vals = Any[]
                    JSON2.peekbyte(io) == JSON2.CLOSE_CURLY_BRACE && (JSON2.readbyte(io); @goto done)
                    typemap = $types
                    while true
                        key = JSON2.read(io, Symbol)
                        push!(keys, key)
                        JSON2.wh!(io)
                        JSON2.@expect ':'
                        JSON2.wh!(io)
                        push!(vals, JSON2.read(io, get(typemap, key, Any))) # recursively reads value
                        JSON2.wh!(io)
                        JSON2.@expectoneof ',' '}'
                        b == JSON2.CLOSE_CURLY_BRACE && @goto done
                        JSON2.wh!(io)
                    end
                    @label done
                    return T(; $defaults..., NamedTuple{Tuple(keys)}(Tuple(vals))...)
                end
                # @show q
                return q
            end
            $wr
        end)
    elseif anydefaults
        q = esc(quote
            @generated function JSON2.read(io::IO, T::Type{$T})
                N = fieldcount($T)
                fieldformats = Dict($(args...))
                names = Tuple(Symbol(get(get(fieldformats, nm, NamedTuple()), :name, nm)) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false))
                types = Tuple(fieldtype($T, i) for i = 1:fieldcount($T) if !get(get(fieldformats, fieldname($T, i), NamedTuple()), :exclude, false))
                jsontypes = Tuple(get(get(fieldformats, fieldname($T, i), NamedTuple()), :jsontype, fieldtype($T, i)) for i = 1:fieldcount($T) if !get(get(fieldformats, fieldname($T, i), NamedTuple()), :exclude, false))
                defaults = (; ((get(get(fieldformats, nm, NamedTuple()), :name, nm), get(fieldformats, nm, NamedTuple())[:default]) for nm in fieldnames($T) if !get(get(fieldformats, nm, NamedTuple()), :exclude, false) && haskey(get(fieldformats, nm, NamedTuple()), :default))...)
                return JSON2.generate_missing_read_body(names, types, jsontypes, defaults)
            end
            $wr
        end)
        # @show q
    else
        q = esc(quote
            @generated function JSON2.read(io::IO, T::Type{$T})
                fieldformats = Dict($(args...))
                types = Tuple(fieldtype($T, i) for i = 1:fieldcount($T) if !get(get(fieldformats, fieldname($T, i), NamedTuple()), :exclude, false))
                jsontypes = Tuple(get(get(fieldformats, fieldname($T, i), NamedTuple()), :jsontype, fieldtype($T, i)) for i = 1:fieldcount($T) if !get(get(fieldformats, fieldname($T, i), NamedTuple()), :exclude, false))
                N = length(types)
                return JSON2.generate_default_read_body(N, types, jsontypes, $T <: NamedTuple)
            end
            $wr
        end)
    end
    # @show q
    return q
end

end # module


text:
function print_header(o::Optimizer)
    return error("print_header not implemented for optimizer $(typeof(o))")
end

function update_fg∇f!(state::OptimizerState, pb)
    return error("update_fg∇f! not implemented for optimizer state $(typeof(state)) and problem $(typeof(pb)).")
end

function update_iterate!(state, pb, optimizer)
    return error("update_iterate! not implemented for optimizer state $(typeof(state)), problem $(typeof(pb)) and optimizer $(typeof(optimizer)).")
end

#
## Print and logs
#

display_logs_header_pre(o) = nothing
display_logs_header_post(o) = nothing

function display_logs_header(o::Optimizer)
    display_logs_header_pre(o::Optimizer)
    print("it.   F(x)                    f(x)       g(x)       step         tgt ∇f+g   nml ∇f+g     Manifold                          ")
    display_logs_header_post(o::Optimizer)
    println()
    return
end

display_logs_pre(optimizer, state, pb) = "0"
display_logs_post(optimizer, state, pb) = ""

function display_logs(state, pb, optimizer, iteration, time, optimstate_extensions, tracing)
    x = get_repr(state.x)
    x_old = get_repr(state.x_old)

    normstep = norm(embed(state.x) - embed(state.x_old))
    minsubgradient_tan, minsubgradient_norm = firstorder_optimality_tangnorm(pb, get_repr(state.x), state.M, state.∇f_x)

    if tracing
        F_x = state.f_x + state.g_x

        linestyle = display_logs_pre(optimizer, state, pb)
        print("\033[$(linestyle)m")

        F_x = state.f_x + state.g_x
        @printf "%4i  %.16e  %-.3e  %-.3e  %-.3e    %-.3e  % .3e " iteration F_x state.f_x state.g_x normstep minsubgradient_tan minsubgradient_norm

        ## underlined manifold if state changed
        manstyle = (state.M == state.M_old) ? linestyle : "4"
        if length(string(state.M)) <= 35
            print(repeat(" ", 35-length(string(state.M))), "\033[$(manstyle)m",state.M, "\033[0m\033[$(linestyle)m")
        else
            print("\033[$(manstyle)m",state.M, "\033[0m\033[$(linestyle)m")
        end
        print("  ")

        display_logs_post(optimizer, state, pb)

        print("\033[0m")
        println()
    end

    return build_optimstate(optimizer, state, iteration, time, normstep, minsubgradient_tan, minsubgradient_norm, optimstate_extensions)
end

function optimize!(
    pb::CompositeProblem,
    optimizer::O,
    initial_x;
    manifold = wholespace_manifold(pb.regularizer, initial_x),
    state::S = initial_state(optimizer, initial_x, pb.regularizer, manifold=manifold),
    optimstate_extensions = [],
    optparams = OptimizerParams()
) where {O<:Optimizer,S<:OptimizerState}

    if getfield(StructuredSolvers, :timeit_debug_enabled)()
        reset_timer!()
    end

    ## Collecting parameters
    @unpack iterations_limit, time_limit = optparams
    @unpack show_trace, trace_length = optparams
    tracing = show_trace

    iteration = 0
    converged = false
    stopped = false
    time_count = 0.0

    update_fg∇f!(state, pb)

    show_trace && print_header(optimizer)
    show_trace && display_logs_header(optimizer)
    optimizationstate = display_logs(
        state,
        pb,
        optimizer,
        iteration,
        time_count,
        optimstate_extensions,
        tracing,
    )
    tr = Vector{OptimizationState}([optimizationstate])

    while !converged && !stopped && iteration < iterations_limit
        iteration += 1

        state.M_old = state.M
        state.x_old = deepcopy(state.x)

        _time = time()
        @timeit_debug "update_iterate" update_iterate!(state, pb, optimizer)

        @timeit_debug "oracles calls" update_fg∇f!(state, pb)
        time_count += time() - _time

        ## Decide if tracing this iteration
        tracing = show_trace && (mod(iteration, ceil(iterations_limit / trace_length)) == 0 || iteration==iterations_limit)

        optimizationstate = display_logs(
            state,
            pb,
            optimizer,
            iteration,
            time_count,
            optimstate_extensions,
            tracing,
        )
        push!(tr, optimizationstate)

        converged = false
        for cvchecker in optparams.cvcheckers
            converged = converged || has_converged(cvchecker, pb, optimizer, optimizationstate)
        end

        # stopped_by_empymanifold = manifold_dimension(state.M) == 0
        stopped_by_time_limit = time_count > time_limit
        stopped = stopped_by_time_limit
    end

    if show_trace && !tracing
        optimizationstate = display_logs(
            state,
            pb,
            optimizer,
            iteration,
            time_count,
            optimstate_extensions,
            true,
        )
        push!(tr, optimizationstate)
    end

    ## Optimality of last iterate
    x = state.x
    M = state.M

    if show_trace
        println("Optimality status of last iterate:")
        res_tan, res_norm = firstorder_optimality_tangnorm(pb, get_repr(state.x), M, state.∇f_x)

        println("- ||Π_tangent(∇f+ḡ)||      : ", res_tan)
        println("- ||Π_normal(∇f+ḡ)||       : ", res_norm)
    end

    if getfield(StructuredSolvers, :timeit_debug_enabled)()
        printstyled("\n\n")
        printstyled("**********************************\n", color=:red)
        print_timer()
        printstyled("**********************************", color=:red)
        printstyled("\n\n")
    end

    return tr
end


text:
using LinearAlgebra: I, inv
using Test: @test, @testset

using BlockArrays: mortar, Block, BlockMatrix
using LinearMaps: LinearMap

using DeltaRCWA

@testset "DeltaRCWA.jl" begin
    @testset "S-matrix  star properties" begin include("smat_star_test.jl") end
    # @testset "Redheffer star properties" begin include("rred_star_test.jl") end
end


text:
using RMInsurance

using Distributions
using DataFrames
include("SSTLife_Input.jl")

#################################################################
println("Start SSTLife.jl")

Random.seed!(seed) ## fix random seed for repeatable results

## Setting up capital market, investments, liabilities ----------
cap_mkt = SSTCapMkt(deepcopy(spot), deepcopy(stock_increase))

assets = Array{Asset}(undef, 0)
for 𝑟𝑜𝑤 ∈ 1:nrow(invest)
  if invest[𝑟𝑜𝑤, :kind] == "stock"
    push!(assets, StockIndex(invest[𝑟𝑜𝑤, :nominal],
                             invest[𝑟𝑜𝑤, :index]))
  else ## zero bond
    push!(assets, ZeroBond(invest[𝑟𝑜𝑤, :nominal],
                           invest[𝑟𝑜𝑤, :maturity],
                           invest[𝑟𝑜𝑤, :index]))
  end
end

liabs =   Liabilities(B_PX, qx, index_mort)

## setting up risk factors
rf = RiskFactor(σ , corr,
                [x0_spot; x0_stock; x0_mort],
                [h_spot; h_stock; h_mort],
                [add_spot; add_stock; add_mort])

## scenarios ----------------------------------------------------
Γ = gammamatrix(assets, liabs, rf, cap_mkt)
δ = delta(assets, liabs, rf, cap_mkt)
for 𝑖 = 1:stress.n
  stress.Δrtk[𝑖] =  Δrtk(vec(stress.Δx[𝑖, :]), δ, Γ)
end

## distribution of the rtk --------------------------------------
## only market risk: _mkt
## onlu insurance risk: _ins
## both market and insurance risk: __mkt_ins
## market and insurance risk, as well as stress: _mkt_ins_stress
index_mkt = [index_spot; index_stock]
index_ins = collect(index_mort)
index_mkt_ins = [index_mkt; index_ins]
r_Δrtk_mkt = rΔrtk(n_scen, assets, liabs, rf, cap_mkt, index_mkt)
r_Δrtk_ins = rΔrtk(n_scen, assets, liabs, rf, cap_mkt, index_ins)
r_Δrtk_mkt_ins =
  rΔrtk(n_scen, assets, liabs, rf, cap_mkt, index_mkt_ins)
r_Δrtk_mkt_ins_stress = aggrstress(stress, r_Δrtk_mkt_ins)

## corresponding expectes shortfall
c_mkt = es(-r_Δrtk_mkt, α)
c_ins = es(-r_Δrtk_ins, α)
c_mkt_ins = es(-r_Δrtk_mkt_ins, α)
c_mkt_ins_stress = es(-r_Δrtk_mkt_ins_stress,α)

## calculation of market value margin ---------------------------
c_fac =
  [value(𝑡, liabs, rf.x0, rf, cap_mkt) /
     value(0, liabs, rf.x0, rf, cap_mkt) for 𝑡 ∈ 1:T]
c_fut_ins_stress =
  c_fac * (c_mkt_ins_stress - c_mkt_ins + c_ins)
mvm_start =
  sum([c_fut_ins_stress[𝑡] * coc_rate / (1 + cap_mkt.spot[𝑡])^𝑡
       for 𝑡 ∈ 1:T])
mvm = mvm_start * (1 + cap_mkt.spot[1])

## target capital and sst ratio ---------------------------------
rtk_start = rtk(0, assets, liabs, rf.x0, rf, cap_mkt)
tc = ( c_mkt_ins_stress - cap_mkt.spot[1] * rtk_start + mvm) /
  (1 + cap_mkt.spot[1])
sst_ratio = rtk_start/tc

## main results ------------------------------------------------
println("rtk_start     :  $(round(rtk_start, digits = 2))")
println("target capital:  $(round(tc, digits = 2))")
println("sst ratio     :  $(round(sst_ratio, digits = 2))")

println("End SSTLife.jl")


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability,
      TaylorIntegration, Plots

include(@relpath "bouncingBallNonlinear_model.jl")
include(@relpath "bouncingBallNonlinear_specifications.jl")

S = bouncingBallNonlinear_model()
X0, options = bouncingBallNonlinear_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# reachability algorithm
algorithm_continuous = TMJets(:orderT => 5, :orderQ => 2, :abs_tol => 1e-10)
algorithm_hybrid = LazyDiscretePost(:check_invariant_intersection => true)

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options, algorithm_continuous, algorithm_hybrid)

# project flowpipe to x₁ and x₂
solution.options[:plot_vars] = [1, 2]
solution_proj = project(solution)

# plot projection
plot(solution_proj, xlab="x", ylab="v")
savefig("bouncingBallNonlinear")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
module test_test_logs

using Test

@test_logs (:info, "foo") @info "foo"

function f()
    @info "bar"
    @info "baz"
    :ok
end
@test (@test_logs (:info, "bar") (:info, "baz") f()) === :ok

end # module test_test_logs


text:
module RateDistortionDecisionMaking

#using #dependencies

using DataFrames, Colors, Gadfly,
      Patchwork,  #there seems to be an issue with Patchwork and firefox right now
      Distances.kl_divergence,
      Distributions,
      Reactive, Interact

export boltzmanndist, BAiterations, setuputilityarrays,
       compute_marginals, marginalizeo, threevarBAiterations,
       mutualinformation, conditional_mutualinformation, expectedutility, entropybits,
       log_bits, kl_divergence_bits,
       RDobjective, ThreeVArRDobjective, analyzeBAsolution, analyze_three_var_BAsolution,
       BAtheme, BAcontinuouscolorscale, BAmatrixvisscale, BAprobabilityvisscale, BAdiscretecolorscale,
       visualizeBAmarginal, visualizeBAconditional, visualizeBAsolution, visualizeMatrix,
       visualizeBA_double_conditional, visualize_three_var_BAsolution,
       plotperformancemeasures, plot_three_var_BA_convergence, plot_three_var_performancemeasures,
       boltzmannresult2DataFrame, BAmarginal2DataFrame, BAconditional2DataFrame,
       BAresult2DataFrame, performancemeasures2DataFrame,
       visualize_three_var_BAsolution, plot_three_var_performancemeasures

#include helper functions for mutual information, expected utiliy
include("InformationTheoryFunctions.jl")

#include BlahutArimoto iterations
include("BlahutArimoto.jl")

#include three-variable iteration schemes
include("ThreeVariableBlahutArimoto.jl")

#include helper functions to aggregate results (as vectors/matrices) into DataFrames
include("ConversionFunctions.jl")

#include helper functions for visualization
include("VisualizationFunctions.jl")

end


text:
#=

Cryptarithmetic puzzle in Julia

Prolog benchmark problem GNU Prolog (crypta.pl)
"""
Name           : crypta.pl
Title          : crypt-arithmetic
Original Source: P. Van Hentenryck's book
Adapted by     : Daniel Diaz - INRIA France
Date           : September 1992

Solve the operation:

B A I J J A J I I A H F C F E B B J E A
+ D H F G A B C D I D B I F F A G F E J E
-----------------------------------------
= G J E G A C D D H F A F J B F I H E E F
"""

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

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function crypta_problem(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
    # "all_optimal_solutions"=>all_solutions,
    "logging"=>[],

    "traverse_strategy"=>:BFS,
    # "traverse_strategy"=>:DFS,
    # "traverse_strategy"=>:DBFS,

    # "branch_split"=>:Smallest,
    # "branch_split"=>:Biggest,
    "branch_split"=>:InHalf,

    # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
    "branch_strategy" => :IMPS, # default
    # "branch_strategy" => :ABS, # Activity Based Search
    # "activity.decay" => 0.999, # default 0.999
    # "activity.max_probes" => 10, # default, 10
    # "activity.max_confidence_deviation" => 20, # default 20

    # "simplify"=>false,
    # "simplify"=>true, # default

    "time_limit"=>6,

    # "backtrack" => false, # default true
    # "backtrack_sorting" => false, # default true

    # "lp_optimizer" => cbc_optimizer,
    # "lp_optimizer" => glpk_optimizer,
    # "lp_optimizer" => ipopt_optimizer,
    ))
    n = 10
    @variable(model, 0 <= x[1:n] <= 9, Int)
    a,b,c,d,e,f,g,h,i,j = x
    @variable(model, sr1, Bin)
    @variable(model, sr2, Bin)

    @constraint(model, x in CS.AllDifferent())

    @constraint(model, b >= 1)
    @constraint(model, d >= 1)
    @constraint(model, g >= 1)

    @constraint(model, a+10*e+100*j+1000*b+10000*b+100000*e+1000000*f+
                      e+10*j+100*e+1000*f+10000*g+100000*a+1000000*f
                      == f+10*e+100*e+1000*h+10000*i+100000*f+1000000*b+10000000*sr1)

    @constraint(model, c+10*f+100*h+1000*a+10000*i+100000*i+1000000*j+
                        f+10*i+100*b+1000*d+10000*i+100000*d+1000000*c+sr1
                        == j+10*f+100*a+1000*f+10000*h+100000*d+1000000*d+10000000*sr2)

    @constraint(model, a+10*j+100*j+1000*i+10000*a+100000*b+
                        b+10*a+100*g+1000*f+10000*h+100000*d+sr2
                        == c+10*a+100*g+1000*e+10000*j+100000*g)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time crypta_problem()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using DICOM

testdir = joinpath(@__DIR__,"testdata")
if !isdir(testdir)
    mkdir(testdir)
end

# TEST SET 1: Simple Reading/Writing

fileMR = joinpath(testdir, "MR_Implicit_Little")
fileCT = joinpath(testdir, "CT_Explicit_Little")
fileMG = joinpath(testdir, "MG_Explicit_Little.zip")

# Don't download files if they already exist
if !isfile(fileMR) && !isfile(fileCT) && !isfile(fileMG)
    download("http://www.barre.nom.fr/medical/samples/files/MR-MONO2-16-head.gz", fileMR*".gz")
    download("http://www.barre.nom.fr/medical/samples/files/CT-MONO2-16-brain.gz", fileCT*".gz")
    download("http://www.dclunie.com/images/pixelspacingtestimages.zip", fileMG)

    run(`gunzip -f $(fileMR*".gz")`)
    run(`gunzip -f $(fileCT*".gz")`)
    run(`unzip -o $fileMG -d $testdir`)
end

# Load dicom data
dcmMR_partial = dcm_parse(fileMR, maxGrp=0x0008)
dcmMR = dcm_parse(fileMR)
dcmCT = dcm_parse(fileCT)
(dcmMG, vrMG) = dcm_parse(joinpath(testdir, "DISCIMG/IMAGES/MGIMAGEA"), true)

@testset "Loading DICOM data" begin
    @test dcmMR_partial[(0x0008,0x0060)] == "MR"
    @test haskey(dcmMR_partial, (0x7FE0,0x0010)) == false

    @test dcmMR[(0x0008,0x0060)] == "MR"
    @test dcmCT[(0x0008,0x0060)] == "CT"
    @test dcmMG[(0x0008,0x0060)] == "MG"

    @test length(dcmMR[(0x7FE0,0x0010)]) == 65536
    @test length(dcmCT[(0x7FE0,0x0010)]) == 262144
    @test length(dcmMG[(0x7FE0,0x0010)]) == 262144

    # Test lookup-by-fieldname
    @test dcmMR[(0x0008,0x0060)] == lookup(dcmMR, "Modality")
    @test dcmMR[(0x7FE0,0x0010)] == lookup(dcmMR, "Pixel Data")
end

# Define two output files for each dcm - data will be saved, reloaded, then saved again
outMR1 = joinpath(testdir,"outMR1.dcm")
outMR2 = joinpath(testdir,"outMR2.dcm")

outCT1 = joinpath(testdir,"outCT1.dcm")
outCT2 = joinpath(testdir,"outCT2.dcm")

outMG1 = joinpath(testdir,"outMG1.dcm")
outMG1b = joinpath(testdir,"outMG1b.dcm")
outMG2 = joinpath(testdir,"outMG2.dcm")

@testset "Writing DICOM data" begin
    # No specific test, just test if file is saved without errors
    outIO = open(outMR1, "w+"); dcm_write(outIO,dcmMR); close(outIO)
    outIO = open(outCT1, "w+"); dcm_write(outIO,dcmCT); close(outIO)
    outIO = open(outMG1, "w+"); dcm_write(outIO,dcmMG,vrMG); close(outIO)
    dcm_write(outMG1b,dcmMG,vrMG)
end

# Load saved DICOM data
dcmMR1 = dcm_parse(outMR1)
dcmCT1 = dcm_parse(outCT1)
(dcmMG1, vrMG1) = dcm_parse(outMG1, true)

@testset "Consistence of Reading/Writing data" begin
    # Confirm that re-loading/saving leads to same file
    outIO = open(outMR2, "w+"); dcm_write(outIO,dcmMR1); close(outIO)
    outIO = open(outCT2, "w+"); dcm_write(outIO,dcmCT1); close(outIO)
    outIO = open(outMG2, "w+"); dcm_write(outIO,dcmMG1, vrMG1); close(outIO)

    @test read(outMR1)==read(outMR2)
    @test read(outCT1)==read(outCT2)
    @test read(outMG1)==read(outMG2)

    # Repeat first tests on reloaded data
    @test dcmMR1[(0x0008,0x0060)] == "MR"
    @test dcmCT1[(0x0008,0x0060)] == "CT"
    @test dcmMG1[(0x0008,0x0060)] == "MG"

    @test length(dcmMR1[(0x7FE0,0x0010)]) == 65536
    @test length(dcmCT1[(0x7FE0,0x0010)]) == 262144
    @test length(dcmMG1[(0x7FE0,0x0010)]) == 262144

    # Test lookup-by-fieldname; cross-compare dcmMR with dcmMR1
    @test dcmMR1[(0x0008,0x0060)] == lookup(dcmMR, "Modality")
    @test dcmMR1[(0x7FE0,0x0010)] == lookup(dcmMR, "Pixel Data")
end

# TEST SET 2: Reading uncommon datasets

# 1. Loading DICOM file with missing header
fileOT = joinpath(testdir, "OT_Implicit_Little_Headless")
if !isfile(fileOT)
    download("http://www.barre.nom.fr/medical/samples/files/OT-MONO2-8-a7.gz", fileOT*".gz")
    run(`gunzip -f $(fileOT*".gz")`)
end

dcmOT = dcm_parse(fileOT, header=false)

# 2. Loading DICOM file with missing header and retired DICOM elements
fileCT = joinpath(testdir, "CT-Implicit_Little_Headless_Retired")
if !isfile(fileCT)
    download("http://www.barre.nom.fr/medical/samples/files/CT-MONO2-12-lomb-an2.gz", fileCT*".gz")
    run(`gunzip -f $(fileCT*".gz")`)
end

# 2a. With user-supplied VRs
dVR_CTa = Dict(
    (0x0008,0x0010) => "SH",
    (0x0008,0x0040) => "US",
    (0x0008,0x0041) => "LO",
    (0x0018,0x1170) => "DS",
    (0x0020,0x0030) => "DS",
    (0x0020,0x0035) => "DS",
    (0x0020,0x0050) => "DS",
    (0x0020,0x0070) => "LO",
    (0x0028,0x0005) => "US",
    (0x0028,0x0040) => "CS",
    (0x0028,0x0200) => "US")
dcmCTa = dcm_parse(fileCT, header=false, dVR=dVR_CTa);

# 2b. With a master VR which skips elements
# Here we skip any element where lookup_vr() fails
# And we also force (0x0018,0x1170) to be read as float instead of integer
dVR_CTb = Dict( (0x0000,0x0000) => "",  (0x0018,0x1170) => "DS")
dcmCTb = dcm_parse(fileCT, header=false, dVR=dVR_CTb);

# 3. Loading DICOM file containing multiple frames
fileMR_multiframe = joinpath(testdir, "MR-Explicit_Little_MultiFrame")
if !isfile(fileMR_multiframe)
    download("http://www.barre.nom.fr/medical/samples/files/MR-MONO2-8-16x-heart.gz", fileMR_multiframe*".gz")
    run(`gunzip -f $(fileMR_multiframe*".gz")`)
end
dcmMR_multiframe = dcm_parse(fileMR_multiframe)

# 4. Load DICOM with unspecified_length()
fileMR_UnspecifiedLength = joinpath(testdir, "MR-UnspecifiedLength")
if !isfile(fileMR_UnspecifiedLength)
    download("https://drive.google.com/uc?export=download&id=1lm0750H-1O22O7Bqy0yfq0FK_vDrqC7-", fileMR_UnspecifiedLength)
end
dcmMR_UnspecifiedLength = dcm_parse(fileMR_UnspecifiedLength)

@testset "Loading uncommon DICOM data" begin
    @test dcmOT[(0x0008,0x0060)] == "OT"

    @test dcmCTa[(0x0008,0x0060)] == "CT"
    @test dcmCTb[(0x0008,0x0060)] == "CT"
    @test haskey(dcmCTa, (0x0028,0x0040)) # dcmCTa should contain retired element
    @test !haskey(dcmCTb, (0x0028,0x0040)) # dcmCTb skips retired elements

    @test dcmMR_multiframe[(0x0008,0x0060)] == "MR"

    @test size(dcmMR_UnspecifiedLength[tag"Pixel Data"]) === (256, 256, 27)
end

@testset "tag" begin
    @test tag"Modality"                 === (0x0008, 0x0060) ===
        DICOM.fieldname_dict["Modality"]
    @test tag"Shutter Overlay Group"    === (0x0018, 0x1623) ===
        DICOM.fieldname_dict["Shutter Overlay Group"]
    @test tag"Histogram Last Bin Value" === (0x0060, 0x3006)
        DICOM.fieldname_dict["Histogram Last Bin Value"]

    # test that compile time error is thrown if tag does not exist
    @test macroexpand(Main,:(tag"Modality")) === (0x0008, 0x0060)
    @test_throws LoadError macroexpand(Main,:(tag"nonsense"))
end


text:
using LinearAlgebra

@testset "core_nonmutating" begin

    # Basic Linear Algebra: cross3, dot3, norm3, norm3!
    a = [1. 2. 9.; 2. 5. 6]
    b = [4. 5. 7.; 3. 2. 8.]
    xab = [-31. 29. -3.; 28. 2. -11.]
    dab = [77., 64.]
    na = [sqrt(86), sqrt(65)]
    @test xab ≈ DER.cross3(a,b)
    @test dab ≈ DER.dot3(a,b)
    @test na ≈ DER.norm3(a)

    # Edges and Tangents
    ea = [1. 3. -3.]
    ta = ea./sqrt(19)
    @test ea ≈ DER.edges(a)
    @test ta ≈ DER.tangents(a)

    # Parallel transport
    θ = [1.0, 0.01, 0.0]
    ptic_θ = [(1 - cos(1.0))/sin(1.0)^2, (1 - cos(0.01))/sin(0.01)^2, 0.5]

    s2 = sin.(θ).^2
    c = cos.(θ)
    @test ptic_θ ≈ DER.ptransport_inner_coefficient.(s2, c)

    t1 = [1. 0. 0.]
    t2 = [1. 1. 1.]/sqrt(3)
    v1 = [0. 1. 1.]
    v2 = [-2. 1. 1.]/sqrt(3)
    @test DER.ptransport(t1, t1, t2) ≈ t2
    @test DER.ptransport(DER.cross3(t1, t2), t1, t2) ≈ DER.cross3(t1, t2)
    @test DER.ptransport([v1; v1], [t1; t1], [t2; t2]) ≈ [v2; v2]

    # Rotations of orthogonal unit vectors
    v1 = [0. 1. 0.; 0. 0. 1.]
    v2 = [0. 1. 1.; 0. -1. 1.]/sqrt(2)
    θ = π/4
    t = [1. 0. 0.]
    @test DER.rotate_orthogonal_unit(v1, [t; t], [θ; θ]) ≈ v2

    # Basic Rod update
    x0 = [0. 0. 0.; 1. 0. 0.; 2. 0. 0.]
    d0 = [0. 0. 1.; 0. 0. 1.]
    Δx = [0. 0. 0.; 0. 0. 1.; 0. 0. 0.]
    Δθ = [0., π/2]
    x, d = DER.rod_update(x0, d0, Δx, Δθ)
    @test x ≈ x0 + Δx
    @test d ≈ [-1/sqrt(2) 0. 1/sqrt(2); 0. -1. 0.]

    # Full kinematics
    l, κ, τ = DER.full_kinematics(x, d)
    @test l ≈ [sqrt(2), sqrt(2)]
    @test κ ≈ [-1. 1.]
    @test τ ≈ [π/2]
end

@testset "core_mutating" begin

    # Basic Linear Algebra: cross3, dot3, norm3, norm3!
    a = [1. 2. 9.; 2. 5. 6]
    b = [4. 5. 7.; 3. 2. 8.]
    xab = [-31. 29. -3.; 28. 2. -11.]
    dab = [77., 64.]
    na = [sqrt(86), sqrt(65)]
    ta = a./na

    c = similar(a)
    d = zeros(2)
    cache = similar(a)

    DER.cross3!(c,a,b)
    @test xab ≈ c

    DER.dot3!(d,a,b)
    @test dab ≈ d

    DER.norm3!(d,a)
    @test na ≈ d

    c = copy(a)
    DER.norm3!(c)
    @test ta ≈ c

    # Edges and tangents
    eb = [-1.0 -3.0 1.0]
    c = similar(eb)
    DER.edges!(c,b)
    @test eb ≈ c

    tb = eb./sqrt(11)
    c = similar(tb)
    DER.tangents!(c,b)
    @test tb ≈ c

    # Parallel transport inner coeff: TBD!
    θ = [1.0, 0.01, 0.0]
    ptic_θ = [(1 - cos(1.0))/sin(1.0)^2, (1 - cos(0.01))/sin(0.01)^2, 0.5]

    s2 = sin.(θ).^2
    c = cos.(θ)
    out = similar(θ)
    DER.ptransport_inner_coefficient!(out, s2, c)
    @test ptic_θ ≈ out

    # Parallel transport
    t1 = [1. 0. 0.]
    t2 = [1. 1. 1.]/sqrt(3)
    v1 = [0. 1. 1.]
    v2 = [-2. 1. 1.]/sqrt(3)
    cache = zeros(1, 9)

    val = similar(t1)
    DER.ptransport!(val, t1, t1, t2, cache)
    @test val ≈ t2

    DER.ptransport!(val, DER.cross3(t1, t2), t1, t2, cache)
    @test val ≈ DER.cross3(t1, t2)

    cache = zeros(2, 9)
    val = [v1; v1]
    DER.ptransport!(val, [v1; v1], [t1; t1], [t2; t2], cache)
    @test val ≈ [v2; v2]

    # Rotations of orthogonal unit vectors
    v1 = [0. 1. 0.; 0. 0. 1.]
    v2 = [0. 1. 1.; 0. -1. 1.]/sqrt(2)
    θ = π/4
    t = [1. 0. 0.]
    cache = zeros(2, 5)
    DER.rotate_orthogonal_unit!(v1, [t; t], [θ; θ], cache)
    @test  v1 ≈ v2

    # Basic Rod update
    x1 = [0. 0. 0.; 1. 0. 0.; 2. 0. 0.]
    d1 = [0. 0. 1.; 0. 0. 1.]
    Δx = [0. 0. 0.; 0. 0. 1.; 0. 0. 0.]
    Δθ = [0., π/2]
    x2 = x1 + Δx
    d2 = [-1/sqrt(2) 0. 1/sqrt(2); 0. -1. 0.]

    xout = similar(x1)
    dout = similar(d1)
    cache = DER.allocate_cache(DER.rod_update!, Float64, 1)
    DER.rod_update!(xout, dout, x1, d1, Δx, Δθ, cache)
    @test xout ≈ x2
    @test dout ≈ d2

    # Full Kinematics
    l = zeros(2)
    κ = zeros(1,2)
    τ = zeros(1)
    caches = DER.allocate_cache(DER.full_kinematics!, Float64, 1)
    DER.full_kinematics!(l, κ, τ, x2, d2, caches)
    @test l ≈ [sqrt(2), sqrt(2)]
    @test κ ≈ [-1. 1.]
    @test τ ≈ [π/2]

end

# Next: Elastic and Gravitational energy