text:
module Y
    f(x) = x + 1
    g(x) = x + 2

    export  g
end

module Z

    f(x) = x + 10
    g(x) = x + 20
    h(x) = x + 30
    k(x) = x + 40
    m(x) = x + 50
    export  h, k

end

module X

    using Reexport

    @reexport using ..Y
    # all of Y's exported symbols available here
     @reexport using ..Z
    # Z's x and y symbols available here

    f(x) = x + 100;
    a(x) = x + 200;
    g(x) = x + 300;

    #export  a

end

module A

    f(x) = "A.f($x)"
end

module B
    g(x) = "B.g($x)"
end
end
using .X

f(1)


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5]
g = LG.grid(dims)
roots = [1, 5]

rsf = PRS.RootedSpanningForest(g, roots)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_prs(rng, rsf)

path = joinpath(@__DIR__, "output", "rooted_spanning_tree", "rooted_spanning_forest.pdf")
p = plot(rsf, sample, dims, path)

ext: jl


text:
module AutoViz

using Reexport
using Parameters
using StaticArrays
using AutomotiveDrivingModels
using Printf
using Rsvg

@reexport using Colors
using Cairo

const DEFAULT_CANVAS_WIDTH = 1000
const DEFAULT_CANVAS_HEIGHT = 600
export DEFAULT_CANVAS_WIDTH, DEFAULT_CANVAS_HEIGHT

global rendermode = :basic # :fancy

function set_render_mode(m::Symbol)
    global rendermode
    rendermode = m
end

include("colorscheme.jl")
export
    COLOR_ASPHALT,
    COLOR_LANE_MARKINGS_WHITE,
    COLOR_LANE_MARKINGS_YELLOW,
    COLOR_CAR_EGO,
    COLOR_CAR_OTHER,
    MONOKAY,
    OFFICETHEME,
    LIGHTTHEME,
    set_color_theme,
    get_pastel_car_colors,
    id_to_color

# Cameras
include("cameras.jl")
export
    CameraState,
    Camera,
    update_camera!,
    StaticCamera,
    TargetFollowCamera,
    ZoomingCamera,
    ComposedCamera,
    SceneFollowCamera,
    set_camera!,
    camera_move!,
    camera_move_pix!,
    camera_rotate!,
    camera_zoom!,
    reset_camera!

include("rendermodels.jl")
export
        RenderModel,
        render_to_canvas,
        add_instruction!,
        reset_instructions!,
        set_background_color!

# Cairo drawing utilities
include("render_instructions.jl")
export
        render_paint,
        render_text,
        render_circle,
        render_arc,
        render_rect,
        render_round_rect,
        render_car,
        render_vehicle,
        render_point_trail,
        render_line,
        render_closed_line,
        render_fill_region,
        render_line_segment,
        render_dashed_line,
        render_arrow,
        render_colormesh,
        grayscale_transform,
        render_fancy_car,
        render_fancy_pedestrian

include("fancy_render.jl")

# renderable interface
include("renderable.jl")

export
    Renderable,
    render,
    render!,
    add_renderable!,
    isrenderable,
    ArrowCar,
    EntityRectangle,
    VelocityArrow,
    FancyCar,
    FancyPedestrian

# Overlays
include("overlays.jl")
export  TextOverlay,
        Overwash,
        NeighborsOverlay,
        CarFollowingStatsOverlay,
        MarkerDistOverlay,
        HistogramOverlay,
        IDOverlay,
        TextParams,
        drawtext,
        LineToCenterlineOverlay,
        LineToFrontOverlay,
        BlinkerOverlay

# Convenient implementation for roadway and vehicle rendering
include("roadways.jl")

end # module

max_stars_repo_name: tawheeler/AutoViz

max_issues_repo_name: tawheeler/AutoViz

max_forks_repo_name: tawheeler/AutoViz


text:
module BlockLinkedLists

mutable struct Node{Block}
    block::Block
    i::Int
    prev::Union{Nothing,Node{Block}}
    next::Union{Nothing,Node{Block}}
end

mutable struct BlockLinkedList{T,Block<:AbstractVector{T},Constructor}
    constructor::Constructor
    blocksize::Int
    head::Node{Block}
    tail::Node{Block}
end

_typeof(::T) where {T} = T
_typeof(::Type{T}) where {T} = Type{T}

function BlockLinkedList(constructor, blocksize::Integer)
    block = constructor(undef, blocksize)
    node = Node(block, firstindex(block), nothing, nothing)
    return BlockLinkedList{eltype(block),typeof(block),_typeof(constructor)}(
        constructor,
        blocksize,
        node,
        node,
    )
end

"""
    alloclast!(list::BlockLinkedList{T}) -> reference::AbstractArray{T,0}

Allocate one element at the end of the `list` and return the `reference` to it
as a 0-dimensional vector.
"""
function alloclast!(lst::BlockLinkedList)
    node = lst.tail
    i = node.i
    if i <= lastindex(node.block)
        block = node.block
        node.i = i + 1
    else
        block = lst.constructor(undef, lst.blocksize)
        i = firstindex(block)
        prev = node
        node = Node(block, i + 1, prev, nothing)
        prev.next = node
        lst.tail = node
    end
    return view(block, i)
end

function Base.push!(lst::BlockLinkedList{T}, x) where {T}
    x = convert(T, x)
    ref = alloclast!(lst)
    ref[] = x
    return lst
end

Base.eltype(::Type{<:BlockLinkedList{T}}) where {T} = T
Base.IteratorEltype(::Type{<:BlockLinkedList}) = Base.HasEltype()
Base.IteratorSize(::Type{<:BlockLinkedList}) = Base.SizeUnknown()

function Base.iterate(lst::BlockLinkedList)
    return iterate(lst, (lst.head, firstindex(lst.head.block), lst.head.i))
end

function Base.iterate(::BlockLinkedList, (node, j, i))
    if j == i
        node = node.next
        node === nothing && return nothing
        j = firstindex(node.block)
        i = node.i
        j == i && return nothing
    end
    return (node.block[j], (node, j + 1, i))
end

end  # module


text:
@testset "Replacement(::String, ::Int, ::Char, ::Bool)" begin
    @test Replacement("Gag", 357, 'S', true) isa Replacement
end

@testset "targets(::Replacement, ::AbstractHLAData)" begin
    fasta_path = joinpath(@__DIR__, "data", "test.fasta")
    hla_types = rand(HLAType, 5)
    hla_data = HLAData("test", fasta_path, hla_types, missing, missing)

    replacement = Replacement("test", 2, 'S', false)
    @test all(Escape.targets(replacement, hla_data) .== [0, 0, 1, 1, 1])

    replacement = Replacement("test", 2, 'G', false)
    @test all(Escape.targets(replacement, hla_data) .== [1, 1, 0, 0, 0])

    replacement = Replacement("test", 3, 'A', false)
    @test all(Escape.targets(replacement, hla_data) .=== [missing, 1, 1, 1, 1])

    replacement = Replacement("test", 3, '-', false)
    @test all(Escape.targets(replacement, hla_data) .== [1, 0, 0, 0, 0])

    replacement = Replacement("test", 10, '-', false)
    @test all(Escape.targets(replacement, hla_data) .== [1, 1, 0, 1, 1])

    replacement = Replacement("test", 10, 'S', false)
    @test all(Escape.targets(replacement, hla_data) .===
        [missing, missing, 1, missing, missing])

    replacement = Replacement("test", 13, 'S', false)
    @test all(Escape.targets(replacement, hla_data) .=== [missing, missing, 1, 1, 1])

    replacement = Replacement("test", 2, 'S', true)
    @test all(Escape.targets(replacement, hla_data) .== [1, 1, 0, 0, 0])
end


text:
function earthquake_exercise_setup()
    # explosion detector (using spiral coordinate system)
    # define the coordinate system:
    S=2000 # number of points on the spiral
    rate=25 # angular rate of spiral
    sd=0.2 # standard deviation of the sensor Gaussian noise

    x=zeros(S); y=zeros(S)
    for s=1:S
        theta=rate*2*pi*s/S;  r=s/S
        x[s]=r*cos(theta); y[s]=r*sin(theta)
    end

    # define the locations of the detection stations on the surface
    # Also define what value on each sensor would be generated by an explostion at internal location s
    N=30 # number of stations
    x_sensor=zeros(N); y_sensor=zeros(N)
    v=zeros(S,N)
    for sensor=1:N
        theta_sensor=2*pi*sensor/N
        x_sensor[sensor]=cos(theta_sensor); y_sensor[sensor]=sin(theta_sensor)
        for s=1:S
            v[s,sensor]=value(x[s],y[s],x_sensor[sensor],y_sensor[sensor]) # explosion value for some value function
        end
    end
end

max_stars_repo_path: Earthqauke - Assignment 1/.ipynb_checkpoints/earthquakeExerciseSetup-checkpoint.jl


text:
function _doi_to_url(doi)
    if startswith(doi, "http")
        return doi
    elseif startswith(doi, "doi.org")
        return string("https://", doi)
    elseif startswith(doi, r"\d{1,10}\.\d{3,20}")
        return string("https://doi.org/", doi)
    end
    return ""
end

function _format_doi_for_filename(doi)
    url = _doi_to_url(doi)
    replace(url, r"[^a-zA-Z0-9]" => "_")
end

# Search a multilayer Dict for a given key
function _find_key(dict::AbstractDict, key)
    for (ikey, val) in dict
        (ikey == key) && return val
        if (val isa AbstractDict)
            val = _find_key(val, key)
            !isnothing(val) && return val
        end
    end
    return nothing
end

# Search a multilayer Dict for a given set of keys
function _find_key(dict::AbstractDict, keys::Vector)
    for key in keys
        dict = _find_key(dict, key)
        isnothing(dict) && break
    end
    return dict
end

function _till_unchanged(f::Function, str0)
    str0 = string(str0)
    while true
        str = f(str0)::String
        (str == str0) && return str
        str0 = str
    end
end

const _CLEAR_BIBTEX_ENTRY_STRIP_SET = ['{', '}', '\n', ' ', ',', ';']
function _clear_bibtex_entry(str0::String)
    return _till_unchanged(str0) do str
        str = strip(str, _CLEAR_BIBTEX_ENTRY_STRIP_SET)
        str = replace(str, "{" => "")
        str = replace(str, "}" => "")
        str = replace(str, "\\" => "")
        str = replace(str, " ," => ",")
        str = replace(str, "`" => "")
        str = replace(str, "'" => "")
        str = replace(str, "  " => " ")
        str = replace(str, "\"" => "")
        return str
    end
end


text:
# add, merge or list split points to an existing table
function table_split(session::AccumuloSession, tablename::String, splits...)
    splitset = Set([bytes(s) for s in splits])
    addSplits(client(session), handle(session), tablename, splitset)
end
table_merge(session::AccumuloSession, tablename::String, start_split, end_split) = mergeTablets(client(session), handle(session), tablename, bytes(start_split), bytes(end_split))
table_splits(session::AccumuloSession, tablename::String, max_splits::Integer=1024) = listSplits(client(session), handle(session), tablename, Int32(max_splits))

# Initiates a major compaction on tablets within the specified range that have one or more files.  If no file selection options are specified, then all files will be compacted.  Options that configure output
# settings are only applied to this compaction and not later compactions.  If multiple concurrent user initiated compactions specify iterators or a compaction strategy, then all but one will fail to start.
#function compact()
#end

# adds, deletes, or lists constraints for a table
constraints(session::AccumuloSession, tablename::String) = listConstraints(client(session), handle(session), tablename)
function add_constraints(session::AccumuloSession, tablename::String, names::String...)
    for name in names
        addConstraint(client(session), handle(session), tablename, name)
    end
end
function remove_constraints(session::AccumuloSession, tablename::String, ids::Integer...)
    for id in ids
        removeConstraint(client(session), handle(session), tablename, Int32(id))
    end
end

# flushes a tables data that is currently in memory to disk
flush(session::AccumuloSession, tablename::String; start_row=UInt8[], end_row=UInt8[], wait::Bool=true) = flushTable(client(session), handle(session), tablename, bytes(start_row), bytes(end_row), wait)

#=
# locality groups not supported yet
# gets the locality groups for a given table
function getgroups()
end

# sets the locality groups for a given table (for binary or commas, use Java API)
function setgroups()
end
=#


text:
module StructC14N

using DataStructures

export canonicalize

######################################################################
# Private functions
######################################################################

"""
  `findabbrv(v::Vector{Symbol})`

  Find all unique abbreviations of symbols in `v`.  Return a tuple of
  two `Vector{Symbol}`: the first contains all possible abbreviations;
  the second contains the corresponding non-abbreviated symbols.
"""
function findabbrv(symLong::Vector{Symbol})
    @assert length(symLong) >= 1
    symStr = String.(symLong)

    outAbbrv = Vector{Symbol}()
    outLong  = Vector{Symbol}()

    # Max length of string representation of keywords
    maxLen = maximum(length.(symStr))

    # Identify all abbreviations
    for len in 1:maxLen
        for i in 1:length(symStr)
            s = symStr[i]
            if length(s) >= len
                s = s[1:len]
                push!(outAbbrv, Symbol(s))
                push!(outLong , symLong[i])
            end
        end
    end

    # Identify unique abbreviations
    for sym in outAbbrv
        i = findall(outAbbrv .== sym)
        count = length(i)
        if count > 1
            deleteat!(outAbbrv, i)
            deleteat!(outLong , i)
        end
    end
    @assert length(unique(outLong)) == length(symLong) "Input symbols have ambiguous abbreviations"

    i = sortperm(outAbbrv)
    return OrderedDict(Pair.(outAbbrv[i], outLong[i]))
end

function myconvert(template, vv)
    if isa(template, Type)
        tt = template
    else
        tt = typeof(template)
    end

    if isa(tt, Union)  &&  (tt != nonmissingtype(tt))
        ismissing(vv)  &&  (return missing)
    end

    if typeof(vv) <: AbstractString  &&  tt <: Number
        return convert(tt, Meta.parse(vv))
    end

    if typeof(vv) <: Number  &&  tt <: AbstractString
        return string(vv)
    end

    if length(methods(parse, (Type{tt}, typeof(vv)))) > 0
        return parse(tt, vv)
    end

    if length(methods(convert, (Type{tt}, typeof(vv)))) > 0
        return convert(tt, vv)
    end

    return tt(vv)
end

######################################################################
# Public functions
######################################################################

# Template is a NamedTuple
"""
  `canonicalize(template::NamedTuple, input::NamedTuple)`

  Canonicalize an input NamedTuple according to a NamedTuple template.
  Returns a NamedTuple.

  Default values must be specified in `template`.  To avoid specifying
  a default value the corresponding item must be a `Type` object, and
  its default value will be `missing` (unless overriden in `input`).

# Examples
```julia-repl
julia> template = (xrange=NTuple{2, Number},
                   yrange=NTuple{2, Number},
                   title="Default string");

julia> c = canonicalize(template, (xr=(1,2),))
(xrange = (1, 2), yrange = missing, title = "Default string")

julia> c = canonicalize(template, (xr=(1,2), yrange=(3, 4.), tit="Foo"))
(xrange = (1, 2), yrange = (3, 4.0), title = "Foo")
```
"""
function canonicalize(template::NamedTuple, input::NamedTuple, dconvert=Dict{Symbol, Function}())
    abbrv = findabbrv(collect(keys(template)))

    # Default values are explicitly provided in `template`.  If a slot
    # in `template` is a `Type` object the corresponding default value
    # is `missing`.
    output = OrderedDict{Symbol, Any}()
    for key in keys(template)
        val = deepcopy(getproperty(template, key))
        if isa(val, Type)
            output[key] = missing  # default value is missing
        else
            output[key] = val      # default value is copied from the template
        end
    end

    # Override with input values
    for key in keys(input)
        val = deepcopy(getproperty(input, key))
        @assert haskey(abbrv, key) "Unexpected key: $key"
        longkey = abbrv[key]
        if haskey(dconvert, longkey)
            output[longkey] = dconvert[longkey](val)
        else
            output[longkey] = myconvert(getproperty(template, longkey), val)
        end
    end
    return NamedTuple(output)
end

# Template is a structure definition
"""
  `canonicalize(template::DataType, input::NamedTuple)`

  Canonicalize an input NamedTuple according to a template in a
  structure definition.  Returns a structure instance.

  If `input` contains less values than required in the template an
  attempt will be made to create a structure with `missing` values.
  If this is forbidden by the structure definition an error is raised.

# Examples
```julia-repl
julia> struct AStruct
           xrange::NTuple{2, Number}
           yrange::NTuple{2, Number}
           title::Union{Missing, String}
       end

julia> canonicalize(AStruct, (xr=(1,2), yr=(3,4.)))
AStruct((1, 2), (3, 4.0), missing)

julia> canonicalize(AStruct, (xr=(1,2), yr=(3,4.), tit="Foo"))
AStruct((1, 2), (3, 4.0), "Foo")
```
"""
function canonicalize(template::DataType, input::NamedTuple, dconvert=Dict{Symbol, Function}())
    @assert isstructtype(template)
    abbrv = findabbrv(collect(fieldnames(template)))

    # Default value is `missing` for all fields
    output = OrderedDict{Symbol, Any}()
    for key in fieldnames(template)
        output[key] = missing
    end

    # Override with input values
    for key in keys(input)
        val = deepcopy(getproperty(input, key))
        @assert haskey(abbrv, key) "Unexpected key: $key"
        longkey = abbrv[key]
        if haskey(dconvert, longkey)
            output[longkey] = dconvert[longkey](val)
        else
            output[longkey] = myconvert(fieldtype(template, longkey), val)
        end
    end
    return template(collect(values(output))...)
end

# Template is a structure instance
"""
  `canonicalize(template, input::NamedTuple)`

  Canonicalize an input NamedTuple according to a template in a
  structure instance.  Returns a structure instance.

  All values in the template are taken as default values.

# Examples
```julia-repl
julia> struct AStruct
           xrange::NTuple{2, Number}
           yrange::NTuple{2, Number}
           title::Union{Missing, String}
       end

julia> template = AStruct((0,0), (0.,0.), missing);

julia> canonicalize(template, (xr=(1,2),))
AStruct((1, 2), (0.0, 0.0), missing)

julia> canonicalize(template, (yr=(1,2), xr=(3.3, 4.4), tit="Foo"))
AStruct((3.3, 4.4), (1, 2), "Foo")
```
"""
function canonicalize(instance, input::NamedTuple, dconvert=Dict{Symbol, Function}())
    template = typeof(instance)
    @assert isstructtype(template)
    abbrv = findabbrv(collect(fieldnames(template)))

    # Default values are copied from template
    output = OrderedDict{Symbol, Any}()
    for key in fieldnames(template)
        output[key] = deepcopy(getfield(instance, key))
    end

    # Override with input values
    for key in keys(input)
        val = deepcopy(getproperty(input, key))
        @assert haskey(abbrv, key) "Unexpected key: $key"
        longkey = abbrv[key]
        if haskey(dconvert, longkey)
            output[longkey] = dconvert[longkey](val)
        else
            output[longkey] = myconvert(fieldtype(template, longkey), val)
        end
    end
    return template(collect(values(output))...)
end

# Input is a Tuple: convert to a NamedTuple using template names and
# invoke previous methods
"""
  `canonicalize(template, input::Tuple)`

  Canonicalize an input tuple according to a template.  Returns a
  NamedTuple or a structure instance according to the type of
  `template`.  The input tuple must have same number of elements as
  the template.

# Examples
```julia-repl
julia> template = (xrange=NTuple{2,Number},
                   yrange=NTuple{2,Number},
                   title="Default string");

julia> canonicalize(template, ((1,2), (3, 4.), "Foo"))
(xrange = (1, 2), yrange = (3, 4.0), title = "Foo")

julia> struct AStruct
           xrange::NTuple{2, Number}
           yrange::NTuple{2, Number}
           title::Union{Missing, String}
       end;

julia> canonicalize(AStruct, ((1,2), (3, 4.), "Foo"))
AStruct((1, 2), (3, 4.0), "Foo")

julia> template = AStruct((0,0), (0.,0.), missing);

julia> canonicalize(template, ((1,2), (3, 4.), "Foo"))
AStruct((1, 2), (3, 4.0), "Foo")
```
"""
function canonicalize(template, input::Tuple, dconvert=Dict{Symbol, Function}())
    if isa(template, NamedTuple)
        @assert length(template) == length(input) "Input tuple must have same length as template"
        return canonicalize(template, NamedTuple{keys(template)}(input), dconvert)
    elseif isa(template, DataType)
        @assert isstructtype(template)
        @assert length(fieldnames(template)) == length(input) "Input tuple must have same number of fields as the template"
        return canonicalize(template, NamedTuple{fieldnames(template)}(input), dconvert)
    elseif isstructtype(typeof(template))
        @assert isstructtype(typeof(template))
        @assert length(fieldnames(typeof(template))) == length(input) "Input tuple must have same number of fields as the template"
        return canonicalize(template, NamedTuple{fieldnames(typeof(template))}(input), dconvert)
    end
    error("Template must be a NamedTuple, a structure definition or a structure instance")
end

# Inputs are given as keywords
"""
  `canonicalize(template, kwargs...)`

  Canonicalize the key/value pairs given as keywords according to the
  `template` structure or named tuple.
"""
function canonicalize(template, dconvert=Dict{Symbol, Function}(); kwargs...)
    return canonicalize(template, NamedTuple(kwargs), dconvert)
end

end # module


text:
using NLsolve, LinearAlgebra
"""
`Base.length(tab::ODERKTableau)`

Defines the length of a Runge-Kutta method to be the number of stages.
"""
Base.length(tab::ODERKTableau) = tab.stages

"""
`stability_region(z,tab::ODERKTableau)`

Calculates the stability function from the tableau at `z`. Stable if <1.

```math
r(z) = \\frac{\\det(I-zA+zeb^T)}{\\det(I-zA)}
```
"""
stability_region(z,tab::ODERKTableau) = det(Matrix{Float64}(I,tab.stages,tab.stages)- z*tab.A + z*ones(tab.stages)*tab.α')/det(Matrix{Float64}(I,tab.stages,tab.stages)-z*tab.A)

"""
`stability_region(tab::ODERKTableau; initial_guess=-3.0)`

Calculates the length of the stability region in the real axis.
"""
function stability_region(tab::ODERKTableau; initial_guess=-3.0)
  residual! = function (resid, x)
    resid[1] = abs(stability_region(x[1], tab)) - 1
  end
  sol = nlsolve(residual!, [initial_guess])
  sol.zero[1]
end


text:
# This file was generated, do not modify it. # hide
using MLJ
import RDatasets: dataset
using PrettyPrinting
using Random
MLJ.color_off() # hide


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays

filter_hash = fink_filter_hash(1, 16, nx=64, pc=2, wd=2)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=128)

function S1_iso_matrix3d(fhash)
    # fhash is the filter hash output by fink_filter_hash
    # The output matrix converts an S1 coeff vector to S1iso by
    #   summing over l
    # Matrix is stored in sparse CSC format using SparseArrays.
    # DPF 2021-Feb-18
    # AKS 2021-Feb-22

    # Does hash contain Omega filter?
    #Omega   = haskey(fhash, "Omega_index")
    #if Omega Ω_ind = fhash["Omega_index"] end

    # unpack fhash
    Nl      = length(fhash["2d_theta_value"])
    Nj      = length(fhash["2d_j_value"])
    Nk      = length(fhash["k_value"])
    Nf      = length(fhash["filt_value"])
    ψ_ind   = fhash["psi_index"]
    #ϕ_ind   = fhash["phi_index"]

    # number of iso coefficients
    #Niso    = Omega ? Nj+2 : Nj+1
    Niso    = (Nj+1)*Nk
    Mat     = zeros(Int32, Niso, Nf)

    # first J elements of iso
    for j = 1:Nj
        for l = 1:Nl
            for k=1:Nk
                λ = ψ_ind[j,l,k]
                Mat[j+(k-1)*Nj, λ] = 1
            end
        end
    end

    j = Nj+1
    l = 1
    for k=1:Nk
        λ = ψ_ind[j,l,k]
        Mat[j+(k-1)*Nj, λ] = 1
    end

    return sparse(Mat)
end

S1_iso_matrix3d(filt_3d)

function S2_iso_matrix3d(fhash)
    # fhash is the filter hash output by fink_filter_hash
    # The output matrix converts an S2 coeff vector to S2iso by
    #   summing over l1,l2 and fixed Δl.
    # Matrix is stored in sparse CSC format using SparseArrays.
    # DPF 2021-Feb-18
    # AKS 2021-Feb-22

    # Does hash contain Omega filter?
    #Omega   = haskey(fhash, "Omega_index")
    #if Omega Ω_ind = fhash["Omega_index"] end

    # unpack fhash
    Nl      = length(fhash["2d_theta_value"])
    Nj      = length(fhash["2d_j_value"])
    Nk      = length(fhash["k_value"])
    Nf      = length(fhash["filt_value"])
    ψ_ind   = fhash["psi_index"]
    #ϕ_ind   = fhash["phi_index"]

    # number of iso coefficients
    #Niso    = Omega ? Nj*Nj*Nl+4*Nj+4 : Nj*Nj*Nl+2*Nj+1
    Niso    = Nj*Nj*Nk*Nk*Nl+2*Nj*Nk*Nk+Nk*Nk
    Mat     = zeros(Int32, Niso, Nf*Nf)

    #should probably think about reformatting to (Nj+1) so phi
    #phi cross terms are in the right place

    # first J*J*K*K*L elements of iso
    for j1 = 1:Nj
        for j2 = 1:Nj
            for l1 = 1:Nl
                for l2 = 1:Nl
                    for k1 = 1:Nk
                        for k2 = 1:Nk
                            DeltaL = mod(l1-l2, Nl)
                            λ1     = ψ_ind[j1,l1,k1]
                            λ2     = ψ_ind[j2,l2,k2]

                            Iiso   = k1+Nk*((k2-1)+Nk*((j1-1)+Nj*((j2-1)+Nj*DeltaL)))
                            Icoeff = λ1+Nf*(λ2-1)
                            Mat[Iiso, Icoeff] = 1
                        end
                    end
                end
            end
        end
    end

    # Next J elements are λϕ, then J elements ϕλ
    for j = 1:Nj
        for l = 1:Nl
            for k1 = 1:Nk
                for k2 =1:Nk
                    λ      = ψ_ind[j,l,k1]
                    Iiso   = Nj*Nj*Nk*Nk*Nl+k1+Nk*((k2-1)+Nk*(j-1))
                    Icoeff = λ+Nf*(ψ_ind[Nj+1,1,k2]-1)  # λϕ
                    Mat[Iiso, Icoeff] = 1

                    Iiso   = Nj*Nj*Nk*Nk*Nl+Nj*Nk*Nk+k1+Nk*((k2-1)+Nk*(j-1))
                    Icoeff = (ψ_ind[Nj+1,1,k2]-1)+Nf*(λ-1)  # ϕλ
                    Mat[Iiso, Icoeff] = 1
                end
            end
        end
    end

    # Next Nk*Nk elements are ϕ(k)ϕ(k)
    j = Nj+1
    l = 1
    for k1=1:Nk
        for k2=1:Nk
            λ1 = ψ_ind[j,l,k1]
            λ2 = ψ_ind[j,l,k2]

            Iiso = Nj*Nj*Nk*Nk*Nl+2*Nj*Nk*Nk+k1+Nk*(k2-1)
            Icoeff = λ1+Nf*(λ2-1)
            Mat[Iiso, Icoeff] = 1
        end
    end

    # # If the Omega filter exists, add more terms
    # if Omega
    #     # Next J elements are λΩ, then J elements Ωλ
    #     for j = 1:Nj
    #         for l = 1:Nl
    #             λ      = ψ_ind[j,l]
    #             Iiso   = I0+j
    #             Icoeff = λ+Nf*(Ω_ind-1)  # λΩ
    #             Mat[Iiso, Icoeff] = 1
    #
    #             Iiso   = I0+Nj+j
    #             Icoeff = Ω_ind+Nf*(λ-1)  # Ωλ
    #             Mat[Iiso, Icoeff] = 1
    #         end
    #     end
    #     # Next 3 elements are ϕΩ, Ωϕ, ΩΩ
    #     Iiso   = I0+Nj+Nj
    #     Mat[Iiso+1, ϕ_ind+Nf*(Ω_ind-1)] = 1
    #     Mat[Iiso+2, Ω_ind+Nf*(ϕ_ind-1)] = 1
    #     Mat[Iiso+3, Ω_ind+Nf*(Ω_ind-1)] = 1
    # end
    #
    return sparse(Mat)
end

S2_iso_matrix3d(filt_3d)

Nj = 4
Nk = 5
Nl = 16
Nj*Nj*Nk*Nk*Nl+2*Nj*Nk*Nk+Nk*Nk

function S1_iso_matrix(fhash)
    # fhash is the filter hash output by fink_filter_hash
    # The output matrix converts an S1 coeff vector to S1iso by
    #   summing over l
    # Matrix is stored in sparse CSC format using SparseArrays.
    # DPF 2021-Feb-18

    # Does hash contain Omega filter?
    Omega   = haskey(fhash, "Omega_index")
    if Omega Ω_ind = fhash["Omega_index"] end

    # unpack fhash
    Nl      = length(fhash["theta_value"])
    Nj      = length(fhash["j_value"])
    Nf      = length(fhash["filt_value"])
    ψ_ind   = fhash["psi_index"]
    ϕ_ind   = fhash["phi_index"]

    # number of iso coefficients
    Niso    = Omega ? Nj+2 : Nj+1
    Mat     = zeros(Int32, Niso, Nf)

    # first J elements of iso
    for j = 1:Nj
        for l = 1:Nl
            λ = ψ_ind[j,l]
            Mat[j, λ] = 1
        end
    end

    # Next elements are ϕ, Ω
    I0     = Nj+1
    Mat[I0, ϕ_ind] = 1
    if Omega Mat[I0+1, Ω_ind] = 1 end

    return sparse(Mat)
end

function S2_iso_matrix(fhash)
    # fhash is the filter hash output by fink_filter_hash
    # The output matrix converts an S2 coeff vector to S2iso by
    #   summing over l1,l2 and fixed Δl.
    # Matrix is stored in sparse CSC format using SparseArrays.
    # DPF 2021-Feb-18

    # Does hash contain Omega filter?
    Omega   = haskey(fhash, "Omega_index")
    if Omega Ω_ind = fhash["Omega_index"] end

    # unpack fhash
    Nl      = length(fhash["theta_value"])
    Nj      = length(fhash["j_value"])
    Nf      = length(fhash["filt_value"])
    ψ_ind   = fhash["psi_index"]
    ϕ_ind   = fhash["phi_index"]

    # number of iso coefficients
    Niso    = Omega ? Nj*Nj*Nl+4*Nj+4 : Nj*Nj*Nl+2*Nj+1
    Mat     = zeros(Int32, Niso, Nf*Nf)

    # first J*J*L elements of iso
    for j1 = 1:Nj
        for j2 = 1:Nj
            for l1 = 1:Nl
                for l2 = 1:Nl
                    DeltaL = mod(l1-l2, Nl)
                    λ1     = ψ_ind[j1,l1]
                    λ2     = ψ_ind[j2,l2]

                    Iiso   = j1+Nj*((j2-1)+Nj*DeltaL)
                    Icoeff = λ1+Nf*(λ2-1)
                    Mat[Iiso, Icoeff] = 1
                end
            end
        end
    end

    # Next J elements are λϕ, then J elements ϕλ
    for j = 1:Nj
        for l = 1:Nl
            λ      = ψ_ind[j,l]
            Iiso   = Nj*Nj*Nl+j
            Icoeff = λ+Nf*(ϕ_ind-1)  # λϕ
            Mat[Iiso, Icoeff] = 1

            Iiso   = Nj*Nj*Nl+Nj+j
            Icoeff = ϕ_ind+Nf*(λ-1)  # ϕλ
            Mat[Iiso, Icoeff] = 1
        end
    end

    # Next 1 element is ϕϕ
    I0     = Nj*Nj*Nl+Nj+Nj+1
    Icoeff = ϕ_ind+Nf*(ϕ_ind-1)
    Mat[I0, Icoeff] = 1

    # If the Omega filter exists, add more terms
    if Omega
        # Next J elements are λΩ, then J elements Ωλ
        for j = 1:Nj
            for l = 1:Nl
                λ      = ψ_ind[j,l]
                Iiso   = I0+j
                Icoeff = λ+Nf*(Ω_ind-1)  # λΩ
                Mat[Iiso, Icoeff] = 1

                Iiso   = I0+Nj+j
                Icoeff = Ω_ind+Nf*(λ-1)  # Ωλ
                Mat[Iiso, Icoeff] = 1
            end
        end
        # Next 3 elements are ϕΩ, Ωϕ, ΩΩ
        Iiso   = I0+Nj+Nj
        Mat[Iiso+1, ϕ_ind+Nf*(Ω_ind-1)] = 1
        Mat[Iiso+2, Ω_ind+Nf*(ϕ_ind-1)] = 1
        Mat[Iiso+3, Ω_ind+Nf*(Ω_ind-1)] = 1
    end

    return sparse(Mat)
end

S1_iso_matrix(filter_hash)
S2_iso_matrix(filter_hash)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using AlphaGo

function translate_gtp_color(gtp_color)
  if lowercase(gtp_color) ∈ ["b", "black"] return go.BLACK end
  if lowercase(gtp_color) ∈ ["w", "white"] return go.WHITE end
  throw(error("invalid color $gtp_color")
end

# GTP command handler for basic play commands
mutable struct GTPHandler
  _komi
  _player
  _courtesy_pass

  function GTPHandler(player, courtesy_pass = false)
    new(6.5, player, courtesy_pass)
  end
end

function cmd_boardsize(gtp_handler::GTPHandler, n::Int)
  n == go.N ? nothing : throw(AssertionError("unsupported board size: $n"))
end

cmd_clear_board(gtp_handler::GTPHandler) =
          initialize_game!(gtp_handler._player, go.Position(komi = self._komi))

cmd_komi(gtp_handler::GTPHandler, komi::Float32)
  gtp_handler._komi = komi
  get_position(gtp_handler._player).komi = komi

function cmd_play(gtp_handler::GTPHandler, arg0::String, arg1 = nothing)
  if arg1 == nothing
    move = arg0
  else
    _accomodate_out_of_turn(gtp_handler, translate_gtp_color(arg0))
    move = arg1
  end
  return play_move!(gtp_handler._player, go.from_kgs(move))
end

function cmd_genmove(gtp_handler::GTPHandler, color = nothing)
  if color is not None:
    self._accomodate_out_of_turn(color)
  end

  if gtp_handler._courtesy_pass
    # If courtesy pass is True and the previous move was a pass, we'll
    # pass too, regardless of score or our opinion on the game.
    pos = get_position(gtp_handler._player)
    if pos.recent && position.recent[end].move == nothing
      return "pass"
    end
  end

  move = suggest_move(gtp_handler._player)
  if gtp_handler._player.should_resign()
    set_result!(gtp_handler._player, -get_position(gtp_handler._player).to_play,
                            was_resign = true)
    return "resign"
  end

  play_move(gtp_handler._player, move)
  if gtp_handler._player.root.is_done()
    set_result!(result(get_position(gtp_handler._player)), was_resign = false)
  end
  return go.to_kgs(move)
end

cmd_undo(get_position(gtp_handler) = throw(error("Not Implemented"))

cmd_final_score(gtp_handler) = gtp_handler._player.result_string

function _accomodate_out_of_turn(gtp_handler::GTPHandler, color::String)
  pos = get_position(gtp_handler._player)
  if translate_gtp_color(color) != pos.to_play
    flip_playerturn!(pos; mutate = true)
  end
end

mutable struct KGSHandler
  _player
end

function cmd_showboard(kgs_handler::KGSHandler)
  println("\n\n" * string(get_position(kgs_handler._player)) * "\n\n")
  return true
end


text:
module LocalizationMicroscopy

using CSV
using Plots

include("data/DataEntity.jl")
include("data/Localization.jl")
include("data/LocalizationGroup.jl")
include("data/InputFormat.jl")
include("loadlocalizations.jl")
include("plotlocalizations.jl")

export loadlocalizations
export InputFormat, nikonelementstext, nikonelementsbinary
export DataEntity, Localization, LocalizationGroup
export extractcoordinates
export plotlocalizations

end # module


text:
module PAGE2020

include("main_model.jl")

end


text:
export mat_back!, mat_back
############### Primitive
"""
The matrix gradient of a rotation block.
"""
@inline function rotgrad(::Type{T}, rb::RotationGate{D}) where {D,T}
    -sin(rb.theta / 2) / 2 * IMatrix{D^nqudits(rb)}() +
    im / 2 * cos(rb.theta / 2) * conj(mat(T, rb.block))
end

function mat_back!(::Type{T}, rb::RotationGate{D,RT}, adjy, collector) where {T,D,RT}
    pushfirst!(collector, projection(rb.theta, sum(adjy .* rotgrad(T, rb))))
end

function mat_back!(::Type{T}, rb::TimeEvolution, adjy, collector) where {T}
    pushfirst!(
        collector,
        projection(rb.dt, im * _sum_A_Bconj(adjy, mat(T, rb.H) * mat(T, rb))),
    )
end

#=
function mat_back!(::Type{T}, A::GeneralMatrixBlock, adjy, collector) where T
    for i in length(adjy):-1:1
        pushfirst!(collector, adjy[i])
    end
end
=#

function mat_back!(::Type{T}, rb::PhaseGate, adjy, collector) where {T}
    s = exp(-im * rb.theta)
    res = -1im * (adjy[1, 1] * s + adjy[2, 2] * s)
    pushfirst!(collector, projection(rb.theta, res))
end

function mat_back!(::Type{T}, rb::ShiftGate, adjy, collector) where {T}
    res = -1im * adjy[2, 2] * exp(-1im * rb.theta)
    pushfirst!(collector, projection(rb.theta, res))
end

######################## Composite
function mat_back!(::Type{T}, rb::AbstractBlock, adjy, collector) where {T}
    nparameters(rb) == 0 && return collector
    throw(MethodError(mat_back!, (T, rb, adjy, collector)))
end

function mat_back!(::Type{T}, rb::PutBlock{D,C,RT}, adjy, collector) where {T,D,C,RT}
    nparameters(rb) == 0 && return collector
    adjm = adjcunmat(adjy, nqudits(rb), (), (), mat(T, content(rb)), rb.locs)
    mat_back!(T, content(rb), adjm, collector)
end

function mat_back!(::Type{T}, rb::Subroutine, adjy, collector) where {T}
    nparameters(rb) == 0 && return collector
    adjm = adjcunmat(adjy, nqudits(rb), (), (), mat(T, content(rb)), rb.locs)
    mat_back!(T, content(rb), adjm, collector)
end

function mat_back!(::Type{T}, rb::CachedBlock, adjy, collector) where {T}
    mat_back!(T, content(rb), adjy, collector)
end

function mat_back!(::Type{T}, rb::Daggered, adjy, collector) where {T}
    mat_back!(T, content(rb), adjy', collector)
end

function mat_back!(::Type{T}, rb::ControlBlock, adjy, collector) where {T}
    nparameters(rb) == 0 && return collector
    adjm = adjcunmat(adjy, nqudits(rb), rb.ctrl_locs, rb.ctrl_config, mat(T, content(rb)), rb.locs)
    mat_back!(T, content(rb), adjm, collector)
end

function mat_back!(::Type{T}, rb::ChainBlock{D}, adjy, collector) where {T,D}
    np = nparameters(rb)
    np == 0 && return collector
    length(rb) == 1 && return mat_back!(T, rb[1], adjy, collector)

    # cache the tape
    mi = mat(T, rb[1])
    cache = Any[mi]
    for b in rb[2:end-1]
        mi = mat(T, b) * mi
        push!(cache, mi)
    end
    adjb = adjy * cache[end]'
    for ib = length(rb):-1:1
        b = rb[ib]
        #adjb = ib==1 ? adjy : adjy*cache[ib]'
        mat_back!(T, b, adjb, collector)
        if ib != 1
            if adjb isa DenseMatrix
                adjb = apply!(ArrayReg{D}(apply!(ArrayReg{D}(adjb), b').state'), rb[ib-1]').state'
            else
                adjb = mat(T, b)' * adjb * mat(T, rb[ib-1])
            end
        end
    end
    return collector
end

function mat_back!(::Type{T}, rb::KronBlock, adjy, collector) where {T}
    nparameters(rb) == 0 && return collector
    mat_back!(T, chain(nqudits(rb), [put(loc => rb[loc]) for loc in rb.locs]), adjy, collector)
    return collector
end

function mat_back!(::Type{T}, rb::Add, adjy, collector) where {T}
    nparameters(rb) == 0 && return collector
    for b in subblocks(rb)[end:-1:1]
        mat_back!(T, b, adjy, collector)
    end
    return collector
end

function mat_back!(::Type{T}, rb::Scale, adjy, collector) where {T}
    np = nparameters(rb)
    np == 0 && return collector
    mat_back!(T, content(rb), factor(rb) * adjy, collector)
    if niparams(rb) > 0
        pushfirst!(collector, projection(rb.alpha, _sum_A_Bconj(adjy, mat(T, content(rb)))))
    end
    return collector
end

# ∑ A .* B*, A is mutated.
_sum_A_Bconj(A::AbstractMatrix, B::AbstractMatrix) = sum(A .* conj.(B))
function _sum_A_Bconj(A::OuterProduct, B::AbstractMatrix)
    res = conj(A.left' * (B * conj(A.right)))
    return ndims(res) != 0 ? sum(diag(res)) : res
end

"""
    mat_back([::Type{T}, ]block::AbstractBlock, adjm::AbstractMatrix) -> Vector

The backward function of `mat`. Returns the gradients of parameters.
"""
mat_back(block::AbstractBlock, adjm::AbstractMatrix{T}) where {T} =
    mat_back!(T, block, adjm, parameters_eltype(block)[])
mat_back(::Type{T}, block::AbstractBlock, adjm::AbstractMatrix) where {T} =
    mat_back!(T, block, adjm, parameters_eltype(block)[])


text:
function process_bout(df::AbstractDataFrame)
    combine(groupby(df,:Patch)) do dd
        @. !ismissing(dd.RewardAvailable)
    end
end


text:
@testset "Validators" begin
    # Defining our missing datasets
    a = allowmissing(1.0:1.0:20.0)
    a[[2, 3, 7]] .= missing
    mask = map(!ismissing, a)

    # We call collect to not have a wrapper type that references the same data.
    m = collect(reshape(a, 5, 4))

    aa = AxisArray(
        deepcopy(m),
        Axis{:time}(DateTime(2017, 6, 5, 5):Hour(1):DateTime(2017, 6, 5, 9)),
        Axis{:id}(1:4)
    )

    table = DataFrame(
        :sin => allowmissing(sin.(1.0:1.0:20.0)),
        :cos => allowmissing(sin.(1.0:1.0:20.0)),
    )

    table.sin[[2, 3, 7, 12, 19]] .= missing

    @testset "Base" begin
        t = Threshold(; limit=0.1)
        @test_throws ThresholdError validate(a, t)
        @test_throws ThresholdError validate(m, t)
        @test_throws ThresholdError validate(aa, t)
        @test_throws ThresholdError validate(table, t)

        # Test showerror
        msg = try
            validate(a, t)
        catch e
            sprint(showerror, e)
        end

        @test msg == "ThresholdError: Missing data limit exceeded 0.1 (0.15)\n"

        t = Threshold(; limit=0.8)
        # Use isequal because we expect the results to contain missings
        @test isequal(validate(a, t), a)
        @test isequal(validate(m, t), m)
        @test isequal(validate(aa, t), aa)
        @test isequal(validate(table, t), table)
    end

    @testset "Weighted" begin
        # If we use an exponentially weighted context then we won't pass the limit
        # because missing earlier observations is less important than later ones.
        t = WeightedThreshold(; limit=0.8, weights=eweights(20, 0.3))
        @test isequal(validate(a, t), a)
        @test isequal(validate(table, t), table)

        @test isequal(wthreshold(m; limit=0.8, weights=eweights(5, 0.3), dims=:cols), m)
        @test isequal(wthreshold(m; limit=0.8, weights=eweights(5, 0.3), dims=:cols), aa)

        # If we reverse the weights such that earlier observations are more important
        # then our previous limit of 0.2 won't be enough to succeed.
        t = WeightedThreshold(; limit=0.1, weights=reverse!(eweights(20, 0.3)))
        @test_throws ThresholdError validate(a, t)
        @test_throws ThresholdError validate(table, t)

        t = WeightedThreshold(; limit=0.1, weights=reverse!(eweights(5, 0.3)))
        @test_throws ThresholdError validate(m, t; dims=:cols)
        @test_throws ThresholdError validate(aa, t; dims=:cols)

        @test_throws DimensionMismatch validate(a[1:10], t)
        @test_throws DimensionMismatch validate(m[1:3, :], t; dims=:cols)
    end

    @testset "functional" begin
        @test_throws ThresholdError Impute.threshold(a; limit=0.1)
        # Use isequal because we expect the results to contain missings
        @test isequal(Impute.threshold(a; limit=0.8), a)
    end
end


text:
"""
    gcmarray{T, N}

gcmarray data structure. Available constructors:

```
gcmarray{T,N}(grid::gcmgrid,f::Array{Array{T,2},N},
         fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1})

gcmarray(grid::gcmgrid,f::Array{Array{T,2},N}) where {T,N}
gcmarray(grid::gcmgrid,f::Array{Array{T,N},1}) where {T,N}

gcmarray(grid::gcmgrid,fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1})
gcmarray(<same as above>,n3::Int)
gcmarray(<same as above>,n3::Int,n4::Int)

gcmarray(grid::gcmgrid)
gcmarray(grid::gcmgrid,::Type{T})
gcmarray(grid::gcmgrid,::Type{T},n3::Int)
gcmarray(grid::gcmgrid,::Type{T},n3::Int,n4::Int)
```
"""
struct gcmarray{T, N} <: AbstractMeshArray{T, N}
   grid::gcmgrid
   f::Array{Array{T,2},N}
   fSize::Array{NTuple{2, Int}}
   fIndex::Array{Int,1}
end

function gcmarray(grid::gcmgrid,f::Array{Array{T,2},N}) where {T, N}
  gcmarray{T,N}(grid,f,grid.fSize,collect(1:grid.nFaces))
end

function gcmarray(grid::gcmgrid,f::Array{Array{T,N},1}) where {T, N}
  nFaces=grid.nFaces
  if N>2
    n3=size(f[1],3); n4=size(f[1],4);
    g=Array{Array{T,2},3}(undef,nFaces,n3,n4)
    for I in eachindex(view(g,1:nFaces,1:n3,1:n4))
      g[I]=view(f[I[1]],:,:,I[2],I[3])
    end
    n4==1 ? g=dropdims(g,dims=3) : nothing
    gcmarray{T,ndims(g)}(grid,g,grid.fSize,collect(1:nFaces))
  else
    gcmarray{T,1}(grid,f,grid.fSize,collect(1:nFaces))
  end
end

function gcmarray(grid::gcmgrid,::Type{T},
        fSize::Union{Array{NTuple{2, Int},1},NTuple{2, Int}},
        fIndex::Union{Array{Int,1},Int}) where {T}
  nFaces=length(fIndex)
  f=Array{Array{T,2},1}(undef,nFaces)
  isa(fSize,NTuple) ? fSize=[fSize] : nothing
  isa(fIndex,Int) ? fIndex=[fIndex] : nothing
  for a=1:nFaces
    f[a]=Array{T}(undef,fSize[a])
  end
  gcmarray{T,1}(grid,f,fSize,fIndex)
end

function gcmarray(grid::gcmgrid,::Type{T},
        fSize::Union{Array{NTuple{2, Int},1},NTuple{2, Int}},
        fIndex::Union{Array{Int,1},Int},n3::Int) where {T}
  nFaces=length(fIndex)
  f=Array{Array{T,2},2}(undef,nFaces,n3)
  isa(fSize,NTuple) ? fSize=[fSize] : nothing
  isa(fIndex,Int) ? fIndex=[fIndex] : nothing
  for a=1:nFaces; for i3=1:n3;
    f[a,i3]=Array{T}(undef,fSize[a]...)
  end; end;
  gcmarray{T,2}(grid,f,fSize,fIndex)
end

function gcmarray(grid::gcmgrid,::Type{T},
        fSize::Union{Array{NTuple{2, Int},1},NTuple{2, Int}},
        fIndex::Union{Array{Int,1},Int},n3::Int,n4::Int) where {T}
  nFaces=length(fIndex)
  f=Array{Array{T,2},3}(undef,nFaces,n3,n4)
  isa(fSize,NTuple) ? fSize=[fSize] : nothing
  isa(fIndex,Int) ? fIndex=[fIndex] : nothing
  for a=1:nFaces; for i4=1:n4; for i3=1:n3;
    f[a,i3,i4]=Array{T}(undef,fSize[a]...)
  end; end; end;
  gcmarray{T,3}(grid,f,fSize,fIndex)
end

# +
function gcmarray(grid::gcmgrid)
  nFaces=grid.nFaces
  fSize=grid.fSize
  fIndex=collect(1:grid.nFaces)
  T=grid.ioPrec
  gcmarray(grid,T,fSize,fIndex)
end

function gcmarray(grid::gcmgrid,::Type{T}) where {T}
  nFaces=grid.nFaces
  fSize=grid.fSize
  fIndex=collect(1:grid.nFaces)
  gcmarray(grid,T,fSize,fIndex)
end

function gcmarray(grid::gcmgrid,::Type{T},n3::Int) where {T}
  nFaces=grid.nFaces
  fSize=grid.fSize
  fIndex=collect(1:grid.nFaces)
  gcmarray(grid,T,fSize,fIndex,n3)
end

function gcmarray(grid::gcmgrid,::Type{T},n3::Int,n4::Int) where {T}
  nFaces=grid.nFaces
  fSize=grid.fSize
  fIndex=collect(1:grid.nFaces)
  gcmarray(grid,T,fSize,fIndex,n3,n4)
end

# -

# # Interface Methods

# +
Base.size(A::gcmarray) = size(A.f)
Base.size(A::gcmarray, dim::Integer) = size(A)[dim]

# +
function Base.getindex(A::gcmarray{T, N}, I::Vararg{Union{Int,Array{Int},AbstractUnitRange,Colon}, N}) where {T,N}
  J=1:length(A.fIndex)
  !isa(I[1],Colon) ? J=J[I[1]] : nothing
  nFaces=length(J)

  tmpf=A.f[I...]
  if isa(tmpf,Array{eltype(A),2})
    tmp=tmpf
  else
    n3=Int(length(tmpf)/nFaces)
    K=(A.grid,eltype(A),A.fSize[J],A.fIndex[J])
    n3>1 ? tmp=gcmarray(K...,n3) : tmp=gcmarray(K...)
    for I in eachindex(tmpf); tmp.f[I] = view(tmpf[I],:,:); end
  end

  return tmp
end

"""
    getindexetc(A::gcmarray, I::Vararg{_}) where {T,N}

Same as getindex but also returns the face size and index
"""
function getindexetc(A::gcmarray{T, N}, I::Vararg{Union{Int,Array{Int},AbstractUnitRange,Colon}, N}) where {T,N}
    f=A[I...]
    fSize=A.fSize[I[1]]
    fIndex=A.fIndex[I[1]]
    return f,fSize,fIndex
end
# -

function Base.setindex!(A::gcmarray{T, N}, v, I::Vararg{Int, N}) where {T,N}
  return (A.f[I...] = v)
end

function Base.view(A::gcmarray{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
  J=1:length(A.fIndex)
  !isa(I[1],Colon) ? J=J[I[1]] : nothing
  nFaces=length(J)

  tmpf=view(A.f,I...)
  n3=Int(length(tmpf)/nFaces) #length(tmpf)>nFaces ? n3=Int(length(tmpf)/nFaces) : n3=1

  K=(A.grid,eltype(A),A.fSize[J],A.fIndex[J])
  n3>1 ? tmp=gcmarray(K...,n3) : tmp=gcmarray(K...)
  for I in eachindex(tmpf); tmp.f[I] = view(tmpf[I],:,:); end

  return tmp
end

# ### Custom pretty-printing, similar, and broadcast

function Base.show(io::IO, z::gcmarray{T, N}) where {T,N}
    printstyled(io, " gcmarray \n",color=:normal)
    printstyled(io, "  grid type   = ",color=:normal)
    printstyled(io, "$(z.grid.class)\n",color=:blue)
    printstyled(io, "  data type   = ",color=:normal)
    printstyled(io, "$(eltype(z))\n",color=:blue)
    printstyled(io, "  tile array  = ",color=:normal)
    printstyled(io, "$(size(z))\n",color=:blue)
    printstyled(io, "  tile sizes  = ",color=:normal)
    printstyled(io, "$(size(z[1]))\n",color=:blue)
    for iFace=2:length(z.fIndex)
      printstyled(io, "                ",color=:normal)
      printstyled(io, "$(size(z[iFace]))\n",color=:blue)
    end
  return
end

function Base.similar(A::gcmarray)
    if ndims(A)==1
        B=gcmarray(A.grid,eltype(A),A.fSize,A.fIndex)
    else
        B=gcmarray(A.grid,eltype(A),A.fSize,A.fIndex,size(A,2))
    end
    return B
end

# ### Customize broadcasting

Base.BroadcastStyle(::Type{<:gcmarray}) = Broadcast.ArrayStyle{gcmarray}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{gcmarray}}, ::Type{ElType}) where ElType
  # Scan the inputs for the gcmarray:
  A = find_gcmarray(bc)
  # Create the gcmarray output:
  if ndims(A)==1
        B=gcmarray(A.grid,ElType,A.fSize,A.fIndex)
  else
        B=gcmarray(A.grid,ElType,A.fSize,A.fIndex,size(A,2))
  end
  return B
end

find_gcmarray(bc::Base.Broadcast.Broadcasted) = find_gcmarray(bc.args)
find_gcmarray(args::Tuple) = find_gcmarray(find_gcmarray(args[1]), Base.tail(args))
find_gcmarray(x) = x
find_gcmarray(a::gcmarray, rest) = a
find_gcmarray(::Any, rest) = find_gcmarray(rest)

####

import Base: copyto!

# Specialize this method if all you want to do is specialize on typeof(dest)
@inline function copyto!(dest::MeshArrays.gcmarray, bc::Broadcast.Broadcasted{Nothing})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    bc′ = Broadcast.preprocess(dest, bc)
    @simd for I in eachindex(bc′)
        #@inbounds dest[I] = bc′[I]
        @inbounds dest[I] = gcmarray_getindex_evalf(bc′,I)
    end
    return dest
end

function gcmarray_getindex_evalf(bc,I)
  @boundscheck checkbounds(bc, I)
  args = Broadcast._getindex(bc.args, I)
  return bc.f.(args...)
end

###

"""
    nFacesEtc(a::gcmarray)

Return nFaces, n3 (1 in 2D case; >1 otherwise)
"""
function nFacesEtc(a::gcmarray)
  nFaces=length(a.fIndex)
  ndims(a.f)>1 ? n3=size(a.f,2) : n3=1
  ndims(a.f)>2 ? n4=size(a.f,3) : n4=1
  return nFaces, n3, n4
end


text:
function linecount(filename)
    count = 0
    for line in eachline(filename)
        count += 1
    end
    count
end

print(linecount("src/wc.jl"))


text:
using OperatorFlux
using Test
using GaussQuadrature
using LinearAlgebra

import OperatorFlux:
    legendre_transform_forward_matrix, legendre_transform_inverse_matrix

@testset "LegendreTransform" begin
    # auxiliary data / legendre polynomials
    P₀(x) = 1
    P₁(x) = x
    P₂(x) = 0.5 * (3 * x^2 - 1)
    P₃(x) = 0.5 * (5 * x^3 - 3 * x)
    P₄(x) = 0.125 * (35 * x^4 - 30 * x^2 + 3)
    P₅(x) = 0.125 * (63 * x^5 - 70 * x^3 + 15 * x)
    polynomials = [P₀, P₁, P₂, P₃, P₄, P₅]

    # test constructors
    trafo = LegendreTransform(modes = (4, 5, 6))
    @test trafo.modes == (4, 5, 6)

    # test forward & inverse
    for mode in [(5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        trafo = LegendreTransform(modes = mode)
        @test norm(
            OperatorFlux.inverse(trafo, OperatorFlux.forward(trafo, z)) - z,
        ) ≤ 20 * eps(norm(z))
    end

    # test truncate_modes & pad_modes
    M = 32
    trafo = LegendreTransform(modes = (3, 7, 6))
    c = rand(7, M, M, M, 3)
    inds = [[1, 2, 3], [1, 2, 3, 4, 5, 6, 7], [1, 2, 3, 4, 5, 6]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = LegendreTransform(modes = (3, 6))
    c = rand(7, M, M, 3)
    inds = [[1, 2, 3], [1, 2, 3, 4, 5, 6]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = LegendreTransform(modes = (3,))
    c = rand(7, M, 3)
    inds = [[1, 2, 3]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = LegendreTransform(modes = (12, 5))
    c = rand(27, 16, 14, 1)
    inds = [collect(1:m) for m in size(c)[2:(end - 1)]]
    @test all(OperatorFlux.pad_modes(trafo, c, (32, 14))[:, inds..., :] .== c)

    # rountrip integration test
    for mode in [(5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        trafo = LegendreTransform(modes = mode)
        a = OperatorFlux.forward(trafo, z)
        b = OperatorFlux.truncate_modes(trafo, a)
        c = OperatorFlux.pad_modes(trafo, b, size(a)[2:(end - 1)])
        d = OperatorFlux.inverse(trafo, c)
        @test norm(d - z) ≤ eps(10 * norm(z))
    end

    # test utils
    for n in 6:16
        x, _ = GaussQuadrature.legendre(n, GaussQuadrature.both)

        inv_transform = legendre_transform_inverse_matrix(n)
        transform = inv(inv_transform)
        @test all(transform .== legendre_transform_forward_matrix(n))

        for (i, p) in enumerate(polynomials)
            @test (transform * p.(x))[i] - 1 ≤ eps(10 * norm(p.(x)))
            @test norm(transform * p.(x)) - 1 ≤ eps(10 * norm(p.(x)))
        end
        for (i, p) in enumerate(polynomials)
            @test norm(inv_transform * transform * p.(x) - p.(x)) ≤
                  eps(10 * norm(p.(x)))
        end
    end

    # test Base extensions
    trafo = LegendreTransform(modes = (5, 6, 12))
    @test ndims(trafo) == 3
    @test eltype(trafo) == Float32
    @test size(trafo) == (5, 6, 12)
end


text:
struct Configserror <:Exception
    msg::String
end

function deleteenvkey!(key::String)
    if haskey(ENV, key)
        delete!(ENV, key)
    end
end

function parseenvkey(key::String, value::Union{String, Number, Bool})
    haskey(ENV, key) ? ENV[key] : value
end

function parsecustomenv!(tree::Dict, key::String, value)
    haskey(ENV, value) ? (tree[key] = ENV[value]) : delete!(tree, key)
end
function parsecustomenv!(tree::Dict)
    for (key, value) in tree
        isa(value, Dict) ? parsecustomenv!(tree[key]) : parsecustomenv!(tree, key, value)
    end
end

makeimmutable(value) = value
function makeimmutable(array::Array)
    shadow = Array{Any, 1}()
    for value in array
        push!(shadow, makeimmutable(value))
    end
    tuple(shadow...)
end
function makeimmutable(dict::Dict)
    shadow = Dict{Symbol, Any}()
    for key in keys(dict)
        shadow[Symbol(key)] = makeimmutable(dict[key])
    end
    (; shadow...)
end

function override!(baseconf::Dict, newconf::Dict)
    for (key, value) in newconf
        isa(value, Dict) && isempty(value) && continue
        !isa(value, Dict) && (baseconf[key] = value; continue)
        (!haskey(baseconf, key) || !isa(baseconf[key], Dict)) && (baseconf[key] = Dict{String, Any}())
        override!(baseconf[key], value)
    end
end

function readconffile(directory::String, file::String)::Dict{String, Any}
    filepath = joinpath(directory, file)
    splitfile = splitext(file)
    ext = splitfile[2]
    if ext === ".json"
        conf = JSON.parsefile(filepath)
    elseif ext === ".yml"
        conf = YAML.load_file(filepath)
        conf = json(conf)
        conf = JSON.parse(conf)
    elseif ext === ".jl"
        conf = include(filepath)
        conf = json(conf)
        conf = JSON.parse(conf)
    end
    filename = splitext(file)[1]
    filename === "custom-environment-variables" && parsecustomenv!(conf)
    conf
end

function getconffiles(path::String)::Array{String, 1}
    !isdir(path) && (path = joinpath(pwd(), path))
    !isdir(path) && throw(Configserror("no such config directory: " * path))
    readdir(path)
end

function parseconfigs(deployment_key::String, configs_directory::String)
    configs_order = [
        "default.json",
        "default.yml",
        "default.jl",
        "custom-environment-variables.json",
        "custom-environment-variables.yml",
        "custom-environment-variables.jl"
    ]
    configs_directory = parseenvkey("CONFIGS_DIRECTORY", configs_directory)
    deployment_key = parseenvkey("DEPLOYMENT_KEY", deployment_key)

    configs_hasfiles = getconffiles(configs_directory)
    configs_files = filter((file)-> file in configs_hasfiles, configs_order)
    deployment = parseenvkey(deployment_key, false)

    if deployment != false
        deployment = lowercase(deployment)
        possiblefiles = [deployment*".json", deployment*".yml", deployment*".jl"]
        possiblefiles = intersect(possiblefiles, configs_hasfiles)
        length(possiblefiles) > 0 && insert!(configs_files, 2, possiblefiles[1])
    end
    (; files = configs_files, directory = configs_directory, )
end

function configpathtodict(path::String, value)::Dict{String, Any}
    path === "" && throw(Configserror("a path is required to set a config"))
    subpaths = split(path, ".")
    base = Dict{String, Any}()
    ref = base
    for i in eachindex(subpaths)
        subpath = subpaths[i]
        length(subpaths) === i && (ref[subpath] = value; return base)
        ref[subpath] = Dict{String, Any}()
        ref = ref[subpath]
    end
end

function parseconfigpath(tree::NamedTuple, path::String)
    paths = split(path, ".")
    key = Symbol(paths[1])
    length(paths) === 1 && haskey(tree, key) && return tree[key]
    (!haskey(tree, key) || !isa(tree[key], NamedTuple)) && throw(Configserror("no such config: " * path))
    key = Symbol(popfirst!(paths))
    parseconfigpath(tree[key], join(paths, "."))
end

function makeconfig!(state::Confstate, value::Dict)
    !state.init && initconfig()
    state.isimmutable && throw(Configserror("""config is immutable. Please set all values before calling "get" """))
    override!(state.configs, value)
end

max_stars_repo_name: citkane/Config

max_issues_repo_name: citkane/Config

max_forks_repo_name: citkane/Config


text:
using Statistics
using RDatasets
data = dataset("datasets", "iris")

# Summary Statistics
describe(data)

# Unique elements in Species
unique(data.Species)
# Standard Deviation of PetalWidth
std(data.PetalWidth)

# Variance of  PetalWidth
var(data.PetalWidth)

#Mean of PetalWidth
mean(data.PetalWidth)

#Median of PetalWidth
median(data.PetalWidth)

#Correlation between PetalWidth & PetalLength
cor(data.PetalWidth, data.PetalLength)

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
using BenchmarkTools

files = [
    "atlas"
    "axes"
    "chain_in_chain"
    "dice_nofriction"
    "dice_tiltedplane"
    "dice"
    "disconnected_bodies"
    "doublependulum_3d"
    "doublependulum_disconnection"
    "doublependulum_urdf"
    "elliptic_joint"
    "fourbar_disconnection"
    "inverted_pyramid_plane"
    "joint_force"
    "joint_torque"
    "linear_pendulum"
    "n_fourbars"
    "n_pendulum"
    "nutation"
    "pendulum_forced"
    "pendulum"
    "planar_example"
    "scissor_lift"
    "slider_crank"
    "slider_crank3d"
    "wheel"
]

controlled = [
    "doublependulum_disconnection"
    "joint_force"
    "joint_torque"
    "pendulum_forced"
    "nutation"
    "fourbar_disconnection"
    "slider_crank"
]

for file in files
    include("examples/"*file*".jl")
    if file ∈ controlled
        storage = simulate!(mech, 10., eval(Meta.parse(file*"_control!")), record = true)
    else
        storage = simulate!(mech, 10., record = true)
    end
    n = length(storage.x)
    for i=1:n
        @test !any(ConstrainedDynamics.sisnan.(storage.x[i]))
        @test !any(ConstrainedDynamics.sisnan.(storage.q[i]))
        @test !any(ConstrainedDynamics.sisnan.(storage.v[i]))
        @test !any(ConstrainedDynamics.sisnan.(storage.ω[i]))
    end
end

include("examples/dummycontroller.jl")
storage = simulate!(mech, 10., dummycontroller_controller, record = true, debug = true)
n = length(storage.x)
for i=1:n
    @test !any(ConstrainedDynamics.sisnan.(storage.x[i]))
    @test !any(ConstrainedDynamics.sisnan.(storage.q[i]))
    @test !any(ConstrainedDynamics.sisnan.(storage.v[i]))
    @test !any(ConstrainedDynamics.sisnan.(storage.ω[i]))
end


text:
loadcolorscheme(:tableau_10, [
    colorant"#4E79A7", # blue
    colorant"#F28E2B", # orange
    colorant"#E15759", # red
    colorant"#76B7B2", # turquoise
    colorant"#59A14F", # green
    colorant"#EDC948", # yellow
    colorant"#B07AA1", # purple
    colorant"#FF9DA7", # pink
    colorant"#9C755F", # brown
    colorant"#BAB0AC", # gray
], "tableau", "sequential")

# loadcolorscheme(:tableau_20, [ # original ordering
#     colorant"#4E79A7", colorant"#A0CBE8", # blue
#     colorant"#F28E2B", colorant"#FFBE7D", # orange
#     colorant"#59A14F", colorant"#8CD17D", # green
#     colorant"#B6992D", colorant"#F1CE63", # yellow
#     colorant"#499894", colorant"#86BCB6", # turquoise
#     colorant"#E15759", colorant"#FF9D9A", # red
#     colorant"#79706E", colorant"#BAB0AC", # gray
#     colorant"#D37295", colorant"#FABFD2", # pink
#     colorant"#B07AA1", colorant"#D4A6C8", # purple
#     colorant"#9D7660", colorant"#D7B5A6", # gray
# ], "tableau", "sequential")

loadcolorscheme(:tableau_20, [
    colorant"#4E79A7", colorant"#A0CBE8", # blue
    colorant"#F28E2B", colorant"#FFBE7D", # orange
    colorant"#E15759", colorant"#FF9D9A", # red
    colorant"#499894", colorant"#86BCB6", # turquoise
    colorant"#59A14F", colorant"#8CD17D", # green
    colorant"#B6992D", colorant"#F1CE63", # yellow
    colorant"#B07AA1", colorant"#D4A6C8", # purple
    colorant"#D37295", colorant"#FABFD2", # pink
    colorant"#9D7660", colorant"#D7B5A6", # brown
    colorant"#79706E", colorant"#BAB0AC", # gray
], "tableau", "sequential")

loadcolorscheme(:tableau_colorblind, [
    colorant"#1170aa",
    colorant"#fc7d0b",
    colorant"#a3acb9",
    colorant"#57606c",
    colorant"#5fa2ce",
    colorant"#c85200",
    colorant"#7b848f",
    colorant"#a3cce9",
    colorant"#ffbc79",
    colorant"#c8d0d9",
], "tableau", "sequential")

loadcolorscheme(:tableau_seattle_grays, [
    colorant"#767f8b",
    colorant"#b3b7b8",
    colorant"#5c6068",
    colorant"#d3d3d3",
    colorant"#989ca3",
], "tableau", "sequential")

loadcolorscheme(:tableau_traffic, [
    colorant"#b60a1c",
    colorant"#e39802",
    colorant"#309143",
    colorant"#e03531",
    colorant"#f0bd27",
    colorant"#51b364",
    colorant"#ff684c",
    colorant"#ffda66",
    colorant"#8ace7e",
], "tableau", "sequential")

loadcolorscheme(:tableau_miller_stone, [
    colorant"#4f6980",
    colorant"#849db1",
    colorant"#a2ceaa",
    colorant"#638b66",
    colorant"#bfbb60",
    colorant"#f47942",
    colorant"#fbb04e",
    colorant"#b66353",
    colorant"#d7ce9f",
    colorant"#b9aa97",
    colorant"#7e756d",
], "tableau", "sequential")

loadcolorscheme(:tableau_superfishel_stone, [
    colorant"#6388b4",
    colorant"#ffae34",
    colorant"#ef6f6a",
    colorant"#8cc2ca",
    colorant"#55ad89",
    colorant"#c3bc3f",
    colorant"#bb7693",
    colorant"#baa094",
    colorant"#a9b5ae",
    colorant"#767676",
], "tableau", "sequential")

loadcolorscheme(:tableau_nuriel_stone, [
    colorant"#8175aa",
    colorant"#6fb899",
    colorant"#31a1b3",
    colorant"#ccb22b",
    colorant"#a39fc9",
    colorant"#94d0c0",
    colorant"#959c9e",
    colorant"#027b8e",
    colorant"#9f8f12",
], "tableau", "sequential")

loadcolorscheme(:tableau_jewel_bright, [
    colorant"#eb1e2c",
    colorant"#fd6f30",
    colorant"#f9a729",
    colorant"#f9d23c",
    colorant"#5fbb68",
    colorant"#64cdcc",
    colorant"#91dcea",
    colorant"#a4a4d5",
    colorant"#bbc9e5",
], "tableau", "sequential")

loadcolorscheme(:tableau_summer, [
    colorant"#bfb202",
    colorant"#b9ca5d",
    colorant"#cf3e53",
    colorant"#f1788d",
    colorant"#00a2b3",
    colorant"#97cfd0",
    colorant"#f3a546",
    colorant"#f7c480",
], "tableau", "sequential")

loadcolorscheme(:tableau_winter,[
    colorant"#90728f",
    colorant"#b9a0b4",
    colorant"#9d983d",
    colorant"#cecb76",
    colorant"#e15759",
    colorant"#ff9888",
    colorant"#6b6b6b",
    colorant"#bab2ae",
    colorant"#aa8780",
    colorant"#dab6af",
], "tableau", "sequential")

loadcolorscheme(:tableau_green_orange_teal, [
    colorant"#4e9f50",
    colorant"#87d180",
    colorant"#ef8a0c",
    colorant"#fcc66d",
    colorant"#3ca8bc",
    colorant"#98d9e4",
    colorant"#94a323",
    colorant"#c3ce3d",
    colorant"#a08400",
    colorant"#f7d42a",
    colorant"#26897e",
    colorant"#8dbfa8",
], "tableau", "sequential")

loadcolorscheme(:tableau_red_blue_brown, [
    colorant"#466f9d",
    colorant"#91b3d7",
    colorant"#ed444a",
    colorant"#feb5a2",
    colorant"#9d7660",
    colorant"#d7b5a6",
    colorant"#3896c4",
    colorant"#a0d4ee",
    colorant"#ba7e45",
    colorant"#39b87f",
    colorant"#c8133b",
    colorant"#ea8783",
], "tableau", "sequential")

loadcolorscheme(:tableau_purple_pink_gray, [
    colorant"#8074a8",
    colorant"#c6c1f0",
    colorant"#c46487",
    colorant"#ffbed1",
    colorant"#9c9290",
    colorant"#c5bfbe",
    colorant"#9b93c9",
    colorant"#ddb5d5",
    colorant"#7c7270",
    colorant"#f498b6",
    colorant"#b173a0",
    colorant"#c799bc",
], "tableau", "sequential")

loadcolorscheme(:tableau_hue_circle, [
    colorant"#1ba3c6",
    colorant"#2cb5c0",
    colorant"#30bcad",
    colorant"#21B087",
    colorant"#33a65c",
    colorant"#57a337",
    colorant"#a2b627",
    colorant"#d5bb21",
    colorant"#f8b620",
    colorant"#f89217",
    colorant"#f06719",
    colorant"#e03426",
    colorant"#f64971",
    colorant"#fc719e",
    colorant"#eb73b3",
    colorant"#ce69be",
    colorant"#a26dc2",
    colorant"#7873c0",
    colorant"#4f7cba",
], "tableau", "sequential")

loadcolorscheme(:tableau_orange_blue,[
    colorant"#9e3d22",
    colorant"#d45b21",
    colorant"#f69035",
    colorant"#d9d5c9",
    colorant"#77acd3",
    colorant"#4f81af",
    colorant"#2b5c8a",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_green, [
    colorant"#a3123a",
    colorant"#e33f43",
    colorant"#f8816b",
    colorant"#ced7c3",
    colorant"#73ba67",
    colorant"#44914e",
    colorant"#24693d",
], "tableau", "diverging")

loadcolorscheme(:tableau_green_blue,[
    colorant"#24693d",
    colorant"#45934d",
    colorant"#75bc69",
    colorant"#c9dad2",
    colorant"#77a9cf",
    colorant"#4e7fab",
    colorant"#2a5783",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_blue, [
    colorant"#a90c38",
    colorant"#e03b42",
    colorant"#f87f69",
    colorant"#dfd4d1",
    colorant"#7eaed3",
    colorant"#5383af",
    colorant"#2e5a87",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_black, [
    colorant"#ae123a",
    colorant"#e33e43",
    colorant"#f8816b",
    colorant"#d9d9d9",
    colorant"#a0a7a8",
    colorant"#707c83",
    colorant"#49525e",
], "tableau", "diverging")

loadcolorscheme(:tableau_gold_purple,[
    colorant"#ad9024",
    colorant"#c1a33b",
    colorant"#d4b95e",
    colorant"#e3d8cf",
    colorant"#d4a3c3",
    colorant"#c189b0",
    colorant"#ac7299",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_green_gold, [
    colorant"#be2a3e",
    colorant"#e25f48",
    colorant"#f88f4d",
    colorant"#f4d166",
    colorant"#90b960",
    colorant"#4b9b5f",
    colorant"#22763f",
], "tableau", "diverging")

loadcolorscheme(:tableau_sunset_sunrise, [
    colorant"#33608c",
    colorant"#9768a5",
    colorant"#e7718a",
    colorant"#f6ba57",
    colorant"#ed7846",
    colorant"#d54c45",
    colorant"#b81840",
], "tableau", "diverging")

loadcolorscheme(:tableau_orange_blue_white, [
    colorant"#9e3d22",
    colorant"#e36621",
    colorant"#fcad52",
    colorant"#ffffff",
    colorant"#95c5e1",
    colorant"#5b8fbc",
    colorant"#2b5c8a",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_green_white, [
    colorant"#ae123a",
    colorant"#ee574d",
    colorant"#fdac9e",
    colorant"#ffffff",
    colorant"#91d183",
    colorant"#539e52",
    colorant"#24693d",
], "tableau", "diverging")

loadcolorscheme(:tableau_green_blue_white, [
    colorant"#24693d",
    colorant"#529c51",
    colorant"#8fd180",
    colorant"#ffffff",
    colorant"#95c1dd",
    colorant"#598ab5",
    colorant"#2a5783",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_blue_white, [
    colorant"#a90c38",
    colorant"#ec534b",
    colorant"#feaa9a",
    colorant"#ffffff",
    colorant"#9ac4e1",
    colorant"#5c8db8",
    colorant"#2e5a87",
], "tableau", "diverging")

loadcolorscheme(:tableau_red_black_white, [
    colorant"#ae123a",
    colorant"#ee574d",
    colorant"#fdac9d",
    colorant"#ffffff",
    colorant"#bdc0bf",
    colorant"#7d888d",
    colorant"#49525e",
], "tableau", "diverging")

loadcolorscheme(:tableau_orange_blue_light, [
    colorant"#ffcc9e",
    colorant"#f9d4b6",
    colorant"#f0dccd",
    colorant"#e5e5e5",
    colorant"#dae1ea",
    colorant"#cfdcef",
    colorant"#c4d8f3",
], "tableau", "diverging")

loadcolorscheme(:tableau_temperature, [
    colorant"#529985",
    colorant"#6c9e6e",
    colorant"#99b059",
    colorant"#dbcf47",
    colorant"#ebc24b",
    colorant"#e3a14f",
    colorant"#c26b51",
], "tableau", "diverging")

loadcolorscheme(:tableau_blue_green, [
    colorant"#feffd9",
    colorant"#f2fabf",
    colorant"#dff3b2",
    colorant"#c4eab1",
    colorant"#94d6b7",
    colorant"#69c5be",
    colorant"#41b7c4",
], "tableau", "sequential")

loadcolorscheme(:tableau_blue_light, [
    colorant"#e5e5e5",
    colorant"#e0e3e8",
    colorant"#dbe1ea",
    colorant"#d5dfec",
    colorant"#d0dcef",
    colorant"#cadaf1",
    colorant"#c4d8f3",
], "tableau", "sequential")

loadcolorscheme(:tableau_orange_light, [
    colorant"#e5e5e5",
    colorant"#ebe1d9",
    colorant"#f0ddcd",
    colorant"#f5d9c2",
    colorant"#f9d4b6",
    colorant"#fdd0aa",
    colorant"#ffcc9e",
], "tableau", "sequential")


text:
using STAC
using Test
using Dates

function testshow(s,substr)
    io = IOBuffer()
    show(io,s)
    @test occursin(substr,String(take!(io)))
end

@testset "STAC.jl" begin

    url = "https://raw.githubusercontent.com/sat-utils/sat-stac/master/test/catalog/catalog.json"

    cat = STAC.Catalog(url)
    @show id(cat)
    @show type(cat)
    @show stac_version(cat)
    @show stac_extensions(cat)
    @show title(cat)
    @show description(cat)
    @show keywords(cat)
    @show license(cat)
    @show providers(cat)
    @show extent(cat)
    @test summaries(cat, default = "foobar") == "foobar"

    @test length(keys(cat)) > 0
    @show keys(cat)

    subcat = cat["stac-catalog-eo"]
    @show subcat
    @show keys(subcat)

    subcat1 = subcat["landsat-8-l1"]

    @show subcat1

    item = subcat1.items["LC08_L1TP_152038_20200611_20200611_01_RT"]

    testshow(item,"box")

    bbox(item)
    geometry(item)
    links(item)
    properties(item)

    assetB4 = item.assets["B4"]
    @show href(assetB4)

    testshow(assetB4,"type")

    STAC.set_cache_max_size(10000)

    for c in eachcatalog(cat)
        @show id(c)
    end

    for item in eachitem(cat)
        @show id(item)
        @test length(bbox(item)) == 4
    end

    STAC.set_title_color(:light_red)
    @test STAC.title_color[] == :light_red

    STAC.set_catalog_color(:green)
    @test STAC.catalog_color[] == :green

    STAC.set_item_color(:cyan)
    @test STAC.item_color[] == :cyan

    STAC.empty_cache()
    @test length(STAC.CACHE) == 0
end

@testset "search" begin
    collections = "landsat-8-c2-l2"
    time_range = (DateTime(2018,01,01), DateTime(2018,01,02))
    lon_range = (2.51357303225, 6.15665815596)
    lat_range = (49.5294835476, 51.4750237087)
    cat = STAC.Catalog("https://planetarycomputer.microsoft.com/api/stac/v1")

    search_results = collect(search(cat, collections, lon_range, lat_range, time_range))
    @test length(search_results) == 2
end


text:
# ------------------------------------------------------------------------------------------
# # 多重派发
# 在这个 notebook 中我们将探索 Julia 的重要特性之一：**多重派发**。
#
# 多重派发让软件更**通用（generic）、更快**！
#
# ## 目录
#
# - [从熟悉的开始讲](#从熟悉的开始讲)
# - [指定传入变量的类型](#指定传入变量的类型)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 从熟悉的开始讲
#
# 从我们已经接触过的知识来理解多重派发。
#
# 我们可以声明一个 Julia 函数而不用给出传入变量的任何信息：
# ------------------------------------------------------------------------------------------

f(x) = x^2

# ------------------------------------------------------------------------------------------
# Julia 自会判断哪些传入变量是合理的，哪些又是不合理的：
# ------------------------------------------------------------------------------------------

f(10)

f([1, 2, 3])

# ------------------------------------------------------------------------------------------
# ## 指定传入变量的类型
#
# 但是，我们*也可以*明确指定传入变量的类型。
#
# 比如，我们写一个只接受字符串类型传入变量的函数 `foo`。
# ------------------------------------------------------------------------------------------

foo(x::String, y::String) = println("My inputs x and y are both strings!")

# ------------------------------------------------------------------------------------------
# 可以看到，要限制 `x` 和 `y` 的类型为 `String`，只要在形参名后面加上两个冒号和关键字 `String` 就行了。
#
# 现在 `foo` 只有在传入变量为 `String` 类型才正常工作。
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# 要让 `foo` 函数接受整型（`Int`）的传入变量，只要在声明 `foo` 函数时在形参后面加上 `::Int`。
# ------------------------------------------------------------------------------------------

foo(x::Int, y::Int) = println("My inputs x and y are both integers!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# 现在函数 `foo` 能接受整型变量了！但是注意，函数 `foo` 仍然能接受字符串变量！
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

# ------------------------------------------------------------------------------------------
# 这就涉及到多重派发的核心了。当我们声明
#
# ```julia
# foo(x::Int, y::Int) = println("My inputs x and y are both integers!")
# ```
# 我们并没有重写或覆盖
# ```julia
# foo(y::String, y::String)
# ```
# 而是给叫做 `foo` 的**通用函数（generic function）**添加了一个**新方法（method）**。
#
# 一个**通用函数** 是一个表示特定操作的抽象概念。
# 举个例子，通用函数 `+` 代表了加法这个概念
#
# 一个**方法** 是一个**通用函数** 接受*特定参数类型* 的一个具体实现。
# 举个例子，`+` 有不同的方法以接受浮点数、整型、矩阵等等。
#
# 我们可以通过函数 `methods` 来查看函数 `foo` 有多少种方法。
# ------------------------------------------------------------------------------------------

methods(foo)

# ------------------------------------------------------------------------------------------
# 猜猜看：加法有多少种方法？
# ------------------------------------------------------------------------------------------

methods(+)

# ------------------------------------------------------------------------------------------
# 所以，我们现在可以在调用 `foo` 的时候传入整型或字符串了。
# 当你调用 `foo` 时传入特定组合的参数时，Julia 将根据传入参数的类型派发对应的方法。*这* 就是多重派发。
#
# 多重派发让我们的代码更加通用更加快。
# 代码更加通用更灵活是因为我们是在写代码描述抽象操作例如加法和乘法，而不是描述具体的实现。
# 同时，因为 Julia 可以为相关类型调用高效的方法，代码可以运行得很快。
#
# 想查看我们调用一个通用函数时派发的是哪个方法，我们可以用宏 `@which`：
# ------------------------------------------------------------------------------------------

@which foo(3, 4)

# ------------------------------------------------------------------------------------------
# 让我们用 `@which` 来看看加号！
# ------------------------------------------------------------------------------------------

@which 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# 可以继续给通用函数 `foo` 增加方法。试试接受**抽象类型** `Number` 作为参数，它包括子类型 `Int`、`Float64` 等数字类型：
# ------------------------------------------------------------------------------------------

foo(x::Number, y::Number) = println("My inputs x and y are both numbers!")

# ------------------------------------------------------------------------------------------
# `foo` 的这个方法可以接受如浮点数作为参数：
# ------------------------------------------------------------------------------------------

foo(3.0, 4.0)

# ------------------------------------------------------------------------------------------
# 我们还可以给 `foo` 加上一个保险的、鸭子类型的方法，让它可以接受任意类型输入：
# ------------------------------------------------------------------------------------------

foo(x, y) = println("I accept inputs of any type!")

# ------------------------------------------------------------------------------------------
# 到目前为止，这个方法会在我们给 `foo` 传入非数字参数时被调用：
# ------------------------------------------------------------------------------------------

v = rand(3)
foo(v, v)

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 9.1
#
# 拓展通用函数 `foo`，给它加上一个只接受1个 `Bool` 类型参数的方法，使它打印"foo with one boolean!"
# ------------------------------------------------------------------------------------------

# 在此作答
function foo()

end

# ------------------------------------------------------------------------------------------
# 确保这个你写的方法在你执行 `foo(true)` 被派发。
# ------------------------------------------------------------------------------------------

@assert foo(true) == "foo with one boolean!"

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
using RegisterQD
using LinearAlgebra
using ImageMagick
using Distributions
using RegisterQD.StaticArrays
using RegisterQD.Interpolations
using RegisterQD.Images
using RegisterQD.CoordinateTransformations
using RegisterQD.Rotations
using RegisterQD.OffsetArrays
using RegisterQD.RegisterMismatch

using Test, TestImages
using Random

#Helper to generate test image pairs
function fixedmov(img, tfm)
    img = float(img)
    img2 = warp(img,tfm)
    inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2)))
    fixed = img[inds...]
    moving = img2[inds...]
    return fixed, moving
end

#helpers to convert Transformations to AffineMaps
to_affine(tfm::Translation) = AffineMap(Matrix{Float64}(LinearAlgebra.I, length(tfm.translation), length(tfm.translation)), tfm.translation)
to_affine(tfm::LinearMap) = AffineMap(Matrix{Float64}(LinearAlgebra.I, length(tfm.translation), length(tfm.translation)), tfm.translation)
to_affine(tfm::AffineMap) = tfm

#Helper to test that a found transform is (roughly) the inverse of the original transform
function tfmtest(tfm, tfminv)
    comp = to_affine(tfm ∘ tfminv)  #should be the identity transform
    diagtol = 0.005
    offdiagtol = 0.005
    vtol = 0.1
    @test all(x->(1-diagtol < x < 1+diagtol), diag(comp.linear))
    @test all(x->(-offdiagtol < x < offdiagtol), comp.linear.-Matrix(Diagonal(diag(comp.linear))))
    @test all(abs.(comp.translation) .< vtol)
end

# tests with standard images
# (Also, unlike the tests above these tests set up the problem so that the correct
# answer to is the inverse of an input transformation.  This seems to catch
# a different set of errors than the tests above)
@testset "QuadDIRECT tests with standard images" begin
    img = testimage("cameraman");

    #Translation (subpixel)
    tfm = Translation(@SVector([14.3, 17.6]))
    fixed, moving = fixedmov(img, tfm)
    mxshift = (100,100) #make sure this isn't too small
    tform, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, rtol=0, fvalue=0.0003)
    tfmtest(tfm, tform)

    #Rigid transform
    SD = Matrix{Float64}(LinearAlgebra.I, 2, 2)
    tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now
    fixed, moving = fixedmov(centered(img), tfm)
    mxshift = (100,100) #make sure this isn't too small
    mxrot = (0.5,)
    minwidth_rot = fill(0.002, 3)
    tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002)
    tfmtest(tfm, tform)
    #with anisotropic sampling
    SD = Matrix(Diagonal([0.5; 1.0]))
    tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.3)*SD)
    fixed, moving = fixedmov(centered(img), tfm)
    tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002)
    tfmtest(tfm, arrayscale(tform, SD))

    #Affine transform
    tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.01))
    #make it harder with nonuniform scaling
    scale = @SMatrix [1.005 0; 0 0.995]
    SD = Matrix{Float64}(LinearAlgebra.I, 2, 2)
    tfm = AffineMap(tfm.linear*scale, tfm.translation)
    mxshift = (100,100) #make sure this isn't too small
    fixed, moving = fixedmov(centered(img), tfm)
    tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=1000, rtol=0, fvalue=0.0002)
    tfmtest(tfm, tform)

    #with anisotropic sampling
    SD = Matrix(Diagonal([0.5; 1.0]))
    tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD)
    scale = @SMatrix [1.005 0; 0 0.995]
    tfm = AffineMap(tfm.linear*scale, tfm.translation)
    tfm = arrayscale(tfm, SD)
    fixed, moving = fixedmov(centered(img), tfm)
    tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=10000, rtol=0, fvalue=0.0002, ndmax = 0.25)
    tform2 = arrayscale(tform, SD)
    tfmtest(tfm, tform2)
end #tests with standard images

@testset "Quadratic interpolation (issue #7)" begin
    samplefrom(n) = rand(Poisson(n))

    Random.seed!(222)
    img = restrict(restrict(testimage("cameraman")))[2:end-1,2:end-1]
    # Convert to "photons" so we can mimic shot noise
    np = 100  # maximum number of photons per pixel
    img = round.(Int, np.*gray.(img))
    fixed  = samplefrom.(img)
    moving = samplefrom.(img)
    ff = qsmooth(fixed)

    tform, mm = qd_translate(fixed, moving, (5, 5); print_interval=typemax(Int))
    tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int))
    @test all(abs.(tformq.translation) .< abs.(tform.translation))

    tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval=typemax(Int))
    tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int))
    @test norm(tformq.linear-I) < norm(tform.linear-I)
    @test norm(tformq.translation) < norm(tform.translation)

    tform, mm = qd_affine(fixed, moving, (5, 5); print_interval=typemax(Int))
    tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int))
    @test mmq < mm

    # Test that we exactly reconstruct `qsmooth` with `presmoothed=true`
    tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int))
    @test all(iszero, tformq.translation)
    @test mmq < 1e-10
    tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int))
    @test mmq < 1e-8
    tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int))
    @test mmq < 1e-6   # on 32-bit systems this can't be 1e-8, not quite sure why
end


text:
export ActionEntropyWrapper, action_entropy_env, SurvivalRewardWrapper, AbstractRewardWrapper

##################################
# Reward Modifying Wrappers
##################################

abstract type AbstractRewardWrapper <: AbstractGymWrapper end
macro reward_fields()
    quote
        wrp_reward::Float64 # reward given by this wrapper at prev `step!`
        wrp_total_reward::Float64 # sum of rewards for the entire ep given only by this wrapper
    end |> esc
end

const reward_defaults = (0.0, 0.0)

"""
`step_reward(wrpenv::AbstractRewardWrapper, r::Real)`
Called at each `step!` to allow modifying the reward for the step
Overload this to specify custom handling of rewards at each `step!`, should
return only the reward that this wrapper applies
"""
function step_reward(wrpenv::AbstractRewardWrapper, r::Real)
    error("step_reward(wrpenv, $r) not implemented for wrpenv of type ", typeof(wrpenv))
end

function reward(wrpenv::AbstractRewardWrapper)
    wrp_reward(wrpenv) + reward(wrpenv.env)
end

function total_reward(wrpenv::AbstractRewardWrapper)
    wrp_total_reward(wrpenv) + total_reward(wrpenv.env)
end

wrp_reward(wrpenv::AbstractRewardWrapper)       = wrpenv.wrp_reward
wrp_total_reward(wrpenv::AbstractRewardWrapper) = wrpenv.wrp_total_reward

function Reinforce.step!(wrpenv::AbstractRewardWrapper, args...)
    r, s = Reinforce.step!(wrpenv.env, args...)
    wrp_reward = step_reward(wrpenv, r) # subtypes implement
    wrpenv.wrp_reward = wrp_reward
    wrpenv.wrp_total_reward += wrp_reward
    wrp_reward, s
end

function Reinforce.reset!(wrpenv::AbstractRewardWrapper)
    wrpenv.wrp_reward = wrpenv.wrp_total_reward = 0.0
    Reinforce.reset!(wrpenv.env)
end

"""
The total bonus/penalty on rewards this wrapper has applied for the current episode
"""
function wrp_total_reward(::Type{T}, env) where {T <: AbstractRewardWrapper}
    wrpenv = getwrapper(T, env)
    wrp_total_reward(wrpenv)
end

survival_reward(env) = wrp_total_reward(SurvivalRewardWrapper, env)
actentropy_reward(env)  = wrp_total_reward(ActionEntropyWrapper, env)
rage_reward(env)     = wrp_total_reward(RageQuitWrapper, env)

gym_reward(env) = total_reward(gymenv(env))

##################################
# SurvivalRewardWrapper()
##################################
"""
```
SurvivalRewardWrapper(env::AbstractGymEnv, reward_thresh::Float64,
                      reward_cap::Float64, min_steps::Int, rstep::Float64)
```
if `total_reward(env)` is less than `reward_thresh`, adds a survival reward of
`rstep` to the reward of the wrapped enviornment at each timestep

The total_reward given in an episode will be between 0.0 and `reward_cap` and
can be accessed with `survival_reward(env)`

`reward_cap` should prob be smaller than the smallest reward an agent can get,
to avoid over encouraging staying alive vs winning.
"""
mutable struct SurvivalRewardWrapper <: AbstractRewardWrapper
    env::AbstractGymEnv
    reward_thresh::Float64
    reward_cap::Float64
    min_steps::Int
    max_steps::Int
    rstep::Float64 # reward per step
    @reward_fields()
end

function SurvivalRewardWrapper(env::AbstractGymEnv, reward_thresh::Float64,
                               reward_cap::Float64, min_steps::Int, rstep::Float64)
    max_steps = floor(Int, 1.0/rstep) + min_steps
    SurvivalRewardWrapper(env, reward_thresh, reward_cap, min_steps, max_steps,
                          rstep, reward_defaults...)
end

"""
SurvivalRewardWrapper(env::AbstractGymEnv, reward_thresh::Float64,
                      reward_cap::Float64, min_steps::Int, max_steps::Int)

Convenience constructor to set `rstep` (reward per step) as `1.0/(max_steps - min_steps)`
"""
function SurvivalRewardWrapper(env::AbstractGymEnv, reward_thresh::Float64,
                               reward_cap::Float64, min_steps::Int, max_steps::Int)
    rstep = 1.0/(max_steps - min_steps)
    SurvivalRewardWrapper(env, reward_thresh, reward_cap, min_steps, max_steps,
                          rstep, reward_defaults...)
end

function step_reward(wrpenv::SurvivalRewardWrapper, r::Real)
    rstep = 0.0
    if total_reward(wrpenv.env) < wrpenv.reward_thresh &&
       wrpenv.min_steps < elapsed_steps(wrpenv) <= wrpenv.max_steps
           rstep = wrpenv.rstep
    end
    rstep
end

elapsed_steps(wrpenv::AbstractGymWrapper) = gymenv(wrpenv).pyenv[:_elapsed_steps]
elapsed_steps(env::GymEnv) = env.pyenv[:_elapsed_steps]

##################################
# ActionEntropyWrapper
##################################
"""
Gives reward for diversity of actions. Variety is the spice of (artificial)
life. The Action entropy wrapper isn't implemented like a standard
AbstractRewardWrapper in that it only changes the total_reward for an episode,
not the reward for each step
"""
mutable struct ActionEntropyWrapper{T} <: AbstractRewardWrapper
    env::AbstractGymEnv
    reward_thresh::Float64
    reward_max::Float64
    min_acceptable_ent_ratio::Float64
    ignore_actions::Vector{<:Any}
    num_actions::Int
    max_ent::Float64
    ep_actions::Vector{T}
    @reward_fields()
end

function ActionEntropyWrapper(
        env::AbstractGymEnv, reward_thresh::Float64, reward_max::Float64,
        min_acceptable_ent_ratio::Float64, ignore_actions::Vector{<:Any},
        actionT::Type=Int
    )
    num_actions = actions(env) |> length
    ActionEntropyWrapper(env, reward_thresh, reward_max, min_acceptable_ent_ratio,
                         ignore_actions, num_actions, -log(1/num_actions),
                         actionT[], reward_defaults...)
end

"""
Get reward between 0.0 and wrpenv.reward_max for mixing it up
"""
function total_reward(wrpenv::ActionEntropyWrapper)
    base_reward = total_reward(wrpenv.env)
    entr_bonus = 0.0
    base_reward < wrpenv.reward_thresh &&
        (entr_bonus = entropy_reward(wrpenv))
    monotony_penalty = 1.0
    if base_reward >= wrpenv.reward_thresh &&
      entropy_ratio(wrpenv) < wrpenv.min_acceptable_ent_ratio
        monotony_penalty = entropy_ratio(wrpenv)
    end
    monotony_penalty * (base_reward + entr_bonus)
end

# Override wrp_total_reward for action entropy since unlike other wrappers, it's
# a non-seperable function of total_reward(wrpenv.env), and so is not accumulated
# in wrpenv.wrp_total_reward at each step
wrp_total_reward(wrpenv::ActionEntropyWrapper) =
    total_reward(wrpenv) - total_reward(wrpenv.env)

"""
More diverse actions => higher entropy => entropy_ratio closer to 1.0
Less diverse actions => lower entropy  => entropy_ratio closer to 0.0
"""
function entropy_ratio(wrpenv::ActionEntropyWrapper)
    entr = proportionmap(wrpenv.ep_actions) |> values |> collect |> entropy
    entr/wrpenv.max_ent # ∈ [0,1]
end

function entropy_reward(wrpenv::ActionEntropyWrapper)
    wrpenv.reward_max * entropy_ratio(wrpenv)
end

function Reinforce.step!(wrpenv::ActionEntropyWrapper, s, a)
    r, s′ = Reinforce.step!(wrpenv.env, s, a)
    !(a in wrpenv.ignore_actions) && push!(wrpenv.ep_actions, a)
    r, s′
end

Reinforce.step!(wrpenv::ActionEntropyWrapper, a) =
    Reinforce.step!(wrpenv::ActionEntropyWrapper, state(wrpenv), a)

function Reinforce.reset!(wrpenv::ActionEntropyWrapper)
    empty!(wrpenv.ep_actions)
    Reinforce.reset!(wrpenv.env)
end

action_entropy_env(wrpenv::ActionEntropyWrapper) = wrpenv
action_entropy_env(wrpenv::AbstractGymWrapper) = action_entropy_env(wrpenv.env)
action_entropy_env(env::OpenAIGym.GymEnv) = env


text:
p = plot(rand(10));
filename = "asdf.tex";

function pgfsave(p, filename)
    filename = abspath(filename)
    savefig(p, filename)
    lines = readlines(filename)
    line = join(lines, "\n")
    @assert occursin(r", width=\{[0-9]+\.[0-9]*mm\}", line)
    @assert occursin(r", height=\{[0-9]+\.[0-9]*mm\}", line)
    new_line = replace(line, r", width=\{[0-9]+\.[0-9]*mm\}" => ", width=\\figwidth")
    new_line = replace(new_line, r", height=\{[0-9]+\.[0-9]*mm\}" => ", height=\\figheight")
    @assert !occursin(r", height=\{[0-9]+\.[0-9]*mm\}", new_line)
    @assert !occursin(r", width=\{[0-9]+\.[0-9]*mm\}", new_line)
    write(filename, new_line)
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos