text:
################################################################################

"""
    process(self::HMM, d::Array{T, 2}, splitSw::Bool;
    args::HMMParams) where T <: Number

# Description
Process hidden Markov model object

# Examples
```julia-repl
julia>

```

See also: [`setup!`](@ref), [`HMM`](@ref), [`HMMParams`](@ref)
"""
function process!(self::HMM, d::Array{T, 2}, splitSw::Bool; args::HMMParams) where T <: Number

  # reset
  reset!(self)

  # feed frame
  for ix in axes(d, 1)
    feed!(self, ix, d, args = args)
  end

  # backtrace
  tb = backTrace(self)

  divider = fill(1, size(self.dataM, 1))
  orig = deepcopy(self.dataM)
  pp = StructArrays.StructArray{ScorePair}(undef, 0)

  # update model
  mdist = zeros(Float64, size(self.dataM))
  mcount = zeros(Float64, size(self.dataM))

  for ix in axes(d, 1)
    self.dataM[tb[ix]] .+= d[ix, :]
    divider[tb[ix]] += 1
    pair = ScorePair(args.distance(orig[tb[ix]], d[ix, :]), ix)

    mdist[tb[ix]] += pair.score
    mcount[tb[ix]] += 1

    push!(pp, pair)
  end

  scores = sort(pp.score, rev = true)
  ixs = map(x -> findall(x .== pp.score), scores)
  pp = pp[vcat(ixs...)]

  for ix in eachindex(self.dataM)
    self.dataM[ix] /= divider[ix]
  end

  sortHMM!(self)

  if !splitSw
    return tb, self.dataM
  end

  max = 0.
  toSplit = 1
  minFreq = 20 # hard coded!

  for ix in eachindex(mdist)
    if mcount[ix] > minFreq
      avdist = mdist[ix] / mcount[ix]
      if avdist > max
        max = avdist
        toSplit = ix
      end
    end
  end

  half = 1 + mcount[toSplit] / 4

  extra = fill(0, size(self.dataM[1], 1))

  count = 1
  for ix in eachindex(pp)
    if tb[pp[ix].index] != toSplit
      continue
    end
    extra += d[pp[ix].index, :]
    count += 1
    if count >= half
      break
    end
  end

  extra ./= (count - 1)

  push!(self.dataM, extra)

  push!(self.tbM, fill(0, size(self.tbM[1], 1)))

  return tb, self.dataM

end

################################################################################

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
#!/usr/bin/julia
# Julia program
#=

This program attempts to organize the X-Plane 11 (more likely do most
of the grunt work) scenery_packs.ini.

X-Plane 11 is a Flight simulation program (some call it a game)

I think "attempt" to organize the config file is the best I can hope for.
X-Plane 11 requires you to manually organize it's scenery_packs.ini when
you want to add/improve features to/in it.
The folders listed within need to be ordered in a specific order, if not
the program will not function properly.
Unfortunately X-Plane 11 has NO folder name standards.
So for example if airports (and airports are just 1 of many catagories)
are required to be ordered first, it would be nice if the word airport or
some other constant like 001 was required to be somewhere in the
folder name...it is NOT! A folder name can be anything
'Fort Lauderdale Executive Airport' or 'ft laundry dale plain takeoff area' or
'KFXE_AP' (KFXE is a valid airport ICAO code) or 'Field of Dreams'
...whatever that developer wanted to name it.

    On: 01/29/2019
    By: Bill Blasingim

01/29/2019 - Currently just laying the ground work!
02/07/2019 - Significant work on identifing  and sorting catagories
02/13/2019 - Added code to allow for testing if a filename contains
a airport ICAO code. I found a airport ICAO file and created a
SQLite airport ICAO code database to check if necessary!

=#
using Dates
using SQLite, DataFrames

fi="/home/bill/Mystuff/SQLite/AirportCodes.db3"
db=SQLite.DB(fi)
#=
sel="select count(*) from AirportCodes where ident='KFXx'"
ex=DataFrame(SQLite.Query(db, sel))
cnt=ex[1][1]
#println(">>>",cnt)
if cnt==0
    println("Not found")
end
exit()
=#

bslash=string('\\')
bslash=Char(92)
#sceneryDir="/media/bill/1450FC356E30C66C/SteamLibrary/steamapps/common/X-Plane 11/Custom Scenery/"
sceneryDir="/home/bill/Documents/Aviation/X-Plane/scenery_packs.ini/"
sceneryFln="scenery_packs11.ini"
suffix=string(now())
inFile=sceneryDir*sceneryFln    #"scenery_packs.ini"
filePieces=splitdir(inFile)

backup=sceneryDir*sceneryFln[1:end-4]*"_"*suffix*".ini"
#cp(inFile,backup)
#println("Backup created..."*backup)

println("Reading: "*inFile)
fi=open(inFile,"r")

outFl=sceneryDir*sceneryFln[1:end-4]*".new"
println("Writing: "*outFl)
fo=open(outFl, "w")

write(fo,"I\n")
write(fo,"1000 Version\n")
write(fo,"SCENERY\n\n")

cnt=recin = 0

###################################################
# Define and initialize arrays                    #
###################################################

noElements=200
apt1Array = Array{String}(undef, noElements)
apt2Array = Array{String}(undef, noElements) #Aerosoft
libArray = Array{String}(undef, noElements)
mshArray = Array{String}(undef, noElements)
ortArray = Array{String}(undef, noElements)
orbArray = Array{String}(undef, noElements)
objArray = Array{String}(undef, noElements)
othArray = Array{String}(undef, noElements)

i=0
for i=1:noElements
    apt1Array[i]=""
    apt2Array[i]="" #Aerosoft
    libArray[i]=""
    mshArray[i]=""
    ortArray[i]=""
    orbArray[i]=""
    objArray[i]=""
    othArray[i]=""
end

libIdx=0
apt1Idx=0
apt2Idx=0
mshIdx=0
ortIdx=0
orbIdx=0
objIdx=0
othIdx=0

###################################################
# Main read loop                                  #
###################################################

while !eof(fi)

    global cnt,recin,apt1Idx,apt2Idx,libIdx,mshIdx,ortIdx,orbIdx,objIdx,othIdx,aptArray
    ln=readline(fi)
    #println(">> "*ln)

    recin += 1

    if length(ln)<13
        #println("Continuing...")
        continue
    end

    if ln[1:13] != "SCENERY_PACK "
        continue
    end

    scnDir=ln[29:end-1] #scenery file
    #println("> "*scnDir[5])

    if length(scnDir)>21 && scnDir[1:21]=="zzz_hd_global_scenery"
        mshIdx+=1
        mshArray[mshIdx]=string(ln)
        continue
    end

###################################################
# Check for library                               #
###################################################
f=findfirst("librar",lowercase(ln))
if f!=nothing
    libIdx+=1
    libArray[libIdx]=string(ln)
    continue
end

###################################################
# Check for objects                               #
###################################################
f=findfirst("cars",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("roads",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("aircraft",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("texture",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("scenery",lowercase(scnDir))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("landmarks",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

f=findfirst("flags",lowercase(ln))
#show(f)
if f!=nothing
    objIdx+=1
    objArray[objIdx]=string(ln)
    continue
end

###################################################
# Check for airport                               #
###################################################

#if length(scnDir)>10 && scnDir[1:11]=="Aerosoft - "
if length(scnDir)>8 && scnDir[1:9]=="Aerosoft "
    apt2Idx+=1
    apt2Array[apt2Idx]=string(ln)
    continue
end

    if scnDir[1]=='K' && scnDir[5]=='_'
        #println("Airport File: "*ln)
        apt1Idx+=1
        apt1Array[apt1Idx]=string(ln)
        continue
    end

    f=findfirst("airport",lowercase(ln))
    if f!=nothing
        apt1Idx+=1
        apt1Array[apt1Idx]=string(ln)
        continue
    end
#=
    f=findfirst("demo",lowercase(ln))
    if f!=nothing
        apt1Idx+=1
        apt1Array[apt1Idx]=string(ln)
        continue
    end
=#
###################################################
# Check for mesh                                  #
###################################################
    f=findfirst("mesh",lowercase(ln))
    if f!=nothing
        mshIdx+=1
        mshArray[mshIdx]=string(ln)
        continue
    end

###################################################
# Check for orbx                                  #
###################################################
    f=findfirst("orbx",lowercase(ln))
    if f!=nothing
        orbIdx+=1
        orbArray[orbIdx]=string(ln)
        continue
    end

###################################################
# Check for ortho                                 #
###################################################
    if findfirst("northo",lowercase(ln))==nothing
        f=findfirst("ortho",lowercase(ln))
        if f!=nothing
            ortIdx+=1
            ortArray[ortIdx]=string(ln)
            continue
        end
    end

    #=
    The following code takes the file name and looks for 4 characters words
    because they might be an airport ICAO code. It takes into account
    that the filename may use "_" instead of a spaces.
    The ICAO code SQLite database, containing 54,000+ records, is then searched.
    =#
    apt=false
    newStr=replace(scnDir, "-" => " ")
    scnDir=newStr
    newStr=replace(scnDir, "_" => " ")
    words=split(newStr) #split string into array
    #println("\nFilename :"*scnDir)
    for word in words
        if length(word)==4
            if all(isletter,word)
                #println("looking for..."*word)
                # I could program exceptions like below but where does it end?
                #if lowercase(word)!="orbx" && lowercase(word)!="mesh"
                sel="select count(*) from AirportCodes where ident='"*word*"'"
                ex=DataFrame(SQLite.Query(db, sel))
                cnt=ex[1][1]
                if cnt>0
                    #println(word*" is an airport!!!!!!!")
                    apt1Idx+=1
                    apt1Array[apt1Idx]=string(ln)
                    apt=true
                    continue
                end
            end
        end
    end

    if apt
        continue
    end

    ############################################
    # Other...everything else                  #
    ############################################
    othIdx+=1
    othArray[othIdx]=string(ln)

end
############################################
# Finished reading file...on to sorting!   #
############################################

###################################################
# Libraries                                       #
###################################################

write(fo,"Libraries [sorted]...",'\n')
noElements=libIdx
arraySrt = Array{String}(undef, noElements)
libIdx=0

while libIdx<noElements
    global libIdx
    libIdx+=1
    arraySrt[libIdx]=libArray[libIdx]
end
sort!(arraySrt)
noElements=libIdx
libIdx=0
while libIdx<noElements
    global libIdx
    libIdx+=1
    write(fo,arraySrt[libIdx],'\n')
end
###################################################
# Airports                                        #
###################################################
#println("Airports [sorted]...")
write(fo,"\nAirports [sorted]...\n")
noElements=apt1Idx
arraySrt = Array{String}(undef, noElements)
apt1Idx=0

#I'm sure this is a bad way to copy an array

while apt1Idx<noElements
    global apt1Idx
    apt1Idx+=1
    arraySrt[apt1Idx]=apt1Array[apt1Idx]
end

sort!(arraySrt)
noElements=apt1Idx
apt1Idx=0
apt=""
glblSw=false
glbl="SCENERY_PACK Custom Scenery/Global Airports/"
while apt1Idx<noElements
    global glblSw,apt1Idx
    apt1Idx+=1
    apt=arraySrt[apt1Idx]
    if apt==glbl
        glblSw=true
    else
        write(fo,apt,'\n')
    end
end

################
# Aerosoft     #
################
noElements=apt2Idx
arraySrt = Array{String}(undef, noElements)
apt2Idx=0

#I'm sure this is a bad way to copy an array

while apt2Idx<noElements
    global apt2Idx
    apt2Idx+=1
    arraySrt[apt2Idx]=apt2Array[apt2Idx]
end

sort!(arraySrt)
noElements=apt2Idx
apt2Idx=0
apt=""
while apt2Idx<noElements
    global glblSw,apt2Idx
    apt2Idx+=1
    apt=arraySrt[apt2Idx]
    write(fo,apt,'\n')
end

if glblSw
    write(fo,glbl,'\n')
end

###################################################
# Objects                                         #
###################################################

write(fo,"\nObjects [sorted]...\n")
noElements=objIdx
arraySrt = Array{String}(undef, noElements)
objIdx=0

#I'm sure this is a bad way to copy an array

while objIdx<noElements
    global objIdx
    objIdx+=1
    arraySrt[objIdx]=objArray[objIdx]
end

sort!(arraySrt)
noElements=objIdx
objIdx=0
while objIdx<noElements
    global objIdx
    objIdx+=1
    write(fo,arraySrt[objIdx],'\n')
end

write(fo,"\nORBX [sorted]...",'\n')
noElements=orbIdx
arraySrt = Array{String}(undef, noElements)
orbIdx=0

###################################################
# ORBX                                            #
###################################################

while orbIdx<noElements
    global orbIdx
    orbIdx+=1
    arraySrt[orbIdx]=orbArray[orbIdx]
end
sort!(arraySrt)
noElements=orbIdx
orbIdx=0
while orbIdx<noElements
    global orbIdx
    orbIdx+=1
    write(fo,arraySrt[orbIdx],'\n')
end

###################################################
# Ortho                                           #
###################################################

write(fo,"\nOrtho [sorted]...",'\n')
noElements=ortIdx
arraySrt = Array{String}(undef, noElements)
ortIdx=0

while ortIdx<noElements
    global ortIdx
    ortIdx+=1
    arraySrt[ortIdx]=ortArray[ortIdx]
end
sort!(arraySrt)
noElements=ortIdx
ortIdx=0
while ortIdx<noElements
    global ortIdx
    ortIdx+=1
    write(fo,arraySrt[ortIdx],'\n')
end

###################################################
# Mesh                                            #
###################################################

write(fo,"\nMesh [sorted]...",'\n')
noElements=mshIdx
arraySrt = Array{String}(undef, noElements)
mshIdx=0

while mshIdx<noElements
    global mshIdx
    mshIdx+=1
    arraySrt[mshIdx]=mshArray[mshIdx]
end
sort!(arraySrt)
noElements=mshIdx
mshIdx=0
while mshIdx<noElements
    global mshIdx
    mshIdx+=1
    write(fo,arraySrt[mshIdx],'\n')
end

###################################################
# Other...unknown!                                #
###################################################

write(fo,"\n\n\n########## Other [sorted]...##########",'\n')
noElements=othIdx
arraySrt = Array{String}(undef, noElements)
othIdx=0

while othIdx<noElements
    global othIdx
    othIdx+=1
    arraySrt[othIdx]=othArray[othIdx]
end
sort!(arraySrt)
noElements=othIdx
othIdx=0
while othIdx<noElements
    global othIdx
    othIdx+=1
    write(fo,arraySrt[othIdx],'\n')
end

close(fi)
close(fo)
#println(recin)
#println(cnt)
#sleep(5)
println("End!")
#exit()

ext: jl


text:
# ------------- #
# Imports/Setup #
# ------------- #

const API_ROOT = "https://api.plot.ly/v2/"
const _VERSION = string(Pkg.installed()["Plotly"])
import JSON: json
import HTTP, JSON
import HTTP: post
function validate_response(res::HTTP.Response)
    if res.status > 204
        url = res.request.headers["Host"] * "/" * res.request.target
        throw(PlotlyAPIError("Request $uri failed with code $(res.status)", res))
    end
end
get_json_data(res::HTTP.Response) = JSON.parse(String(res.body))
get_headers(res::HTTP.Response) = res.headers

# --------------- #
# Tools/Utilities #
# --------------- #
struct PlotlyAPIError <: Exception
    msg
    res
end

function basic_auth(username, password)
    # ref https://github.com/plotly/plotly.py/blob/master/plotly/api/utils.py
    return string("Basic ", base64encode(string(username, ":", password)))
end

function get_req_headers()
    creds = get_credentials()
    return Dict{Any,Any}(
        "Plotly-Client-Platform" => "Julia $(_VERSION)",
        "Content-Type" => "application/json",
        "content-type" => "application/json",
        "Accept" =>  "application/json",  # TODO: for some reason I had to do this to get it to work???
        "authorization" => basic_auth(creds.username, creds.api_key),
    )
end

function get_json(;kwargs...)
    out = Dict()
    for (k, v) in kwargs
        if v !=nothing
            out[k] = v
        end
    end
    out
end

function _get_uid(uids)
    length(uids) == 0 && error("Must supply at least one uid")
    uid = join(map(string, uids), ",")
end

function api_url(endpoint; fid=nothing, route=nothing)
    extra = []
    fid !== nothing && push!(extra, fid)
    route !== nothing && push!(extra, route)
    out = string(API_ROOT, endpoint)
    if length(extra) > 0
        return out *  "/" * join(extra, "/")
    else
        return out
    end
end

function request(method, endpoint; fid=nothing, route=nothing, json=nothing, kwargs...)
    url = api_url(endpoint; fid=fid, route=route)
    query_params = Dict()
    for (k, v) in kwargs
        if v !== nothing
            query_params[string(k)] = v
        end
    end
    # build kwargs for HTTP.request
    kw = Any[:query => query_params, :headers => get_req_headers()]
    json !== nothing && push!(kw, :body => JSON.json(json))

    res = HTTP.request(method, url; kw...)
    validate_response(res)
    res
end

function request_data(method, endpoint; kwargs...)
    res = request(method, endpoint; kwargs...)
    content_type = get(Dict(get_headers(res)), "Content-Type", "application/json")
    if startswith(content_type, "application/json")
        return get_json_data(res)
    else
        return res
    end
end

# ---------------- #
# # API wrappers # #
# ---------------- #

struct ApiCall
    funname::Symbol
    method::Symbol
    endpoint::Symbol
    fid::Bool
    route::Union{Nothing,Symbol}
    required::Vector{Symbol}
    optional::Vector{Symbol}
    json::Vector{Symbol}
    pre_json::Vector{Symbol}
    required_json::Vector{Symbol}
    required_pre_json::Vector{Symbol}
    uids::Bool
    data_out::Bool
end

function ApiCall(
        funname::Symbol, method::Symbol, endpoint::Symbol, fid::Bool=false,
        route=nothing; required=Symbol[], optional=Symbol[], json=Symbol[],
        pre_json=Symbol[], required_json=Symbol[], required_pre_json=Symbol[],
        uids::Bool=false, data_out::Bool=true
    )
    ApiCall(
        funname, method, endpoint, fid, route, required, optional, json,
        pre_json, required_json, required_pre_json, uids, data_out
    )
end

function make_method(api::ApiCall)
    request_fun = api.data_out ? :request_data : :request
    request_kwargs = []

    all_kw = vcat(api.json, api.optional, api.pre_json)
    sig = Expr(:call,
        api.funname,
        Expr(:parameters, [Expr(:kw, name, nothing) for name in all_kw]...),
    )
    if api.fid
        push!(sig.args, :_fid)
        push!(request_kwargs, Expr(:kw, :fid, :_fid))
    end

    # add rest of required
    append!(sig.args, api.required)
    append!(sig.args, api.required_json)
    append!(sig.args, api.required_pre_json)
    if api.uids
        push!(sig.args, Expr(:(...), :uids))
        push!(request_kwargs, Expr(:kw, :uid, :(_get_uid(uids))))
    end

    if api.route != nothing
        push!(request_kwargs, Expr(:kw, :route, string(api.route)))
    end

    call_get_json = Expr(:call, :get_json)
    for name in vcat(api.pre_json, api.required_pre_json)
        push!(call_get_json.args, Expr(:kw, name, Expr(:call, :json, name)))
    end
    for name in vcat(api.json, api.required_json)
        push!(call_get_json.args, Expr(:kw, name, name))
    end

    if length(call_get_json.args) > 1  # we had some json args
        push!(request_kwargs, Expr(:kw, :json, call_get_json))
    elseif api.method in (:PUT, :POST, :PATCH, :delete)
        # need to add empty json argument on put these request methods when
        # no json data is needed.
        push!(request_kwargs, Expr(:kw, :json, :(Dict())))
    end

    # add rest of kwargs
    append!(request_kwargs, [Expr(:kw, name, name) for name in api.optional])
    append!(request_kwargs, [Expr(:kw, name, name) for name in api.required])

    body = Expr(:block,
        Expr(:call,
            request_fun,
            string(api.method),
            string(api.endpoint),
            request_kwargs...,
        )
    )

    if api.data_out
        raw_sig = deepcopy(sig)
        raw_sig.args[1] = Symbol(api.funname, "_", "raw")

        raw_body = deepcopy(body)
        raw_body.args[1].args[1] = :request
        return Expr(:block, Expr(:function, sig, body), Expr(:function, raw_sig, raw_body))

    else
        return Expr(:function, sig, body)
    end
end

file_writeable_metadata = [
    :parent_path, :filename, :parent, :share_key_enabled, :world_readable
]
grid_writable_metadata = vcat(file_writeable_metadata, [:])
dashboard_writeable_metadata = [
    :content, :filename, :parent, :parent_path, :world_readable
]
for _api in [
        # search
        ApiCall(:search_list, :GET, :search, false, required=[:q])

        # files
        ApiCall(:file_retrieve, :GET, :files, true)
        ApiCall(:file_content, :GET, :files, true, :content)  # failing
        ApiCall(:file_update, :PUT, :files, true, json=file_writeable_metadata)
        ApiCall(:file_partial_update, :PATCH, :files, true, json=file_writeable_metadata)
        ApiCall(:file_image, :GET, :files, true, :image)
        ApiCall(:file_copy, :GET, :files, true, :copy, optional=[:deep_copy])  # failing
        ApiCall(:file_path, :GET, :files, true, :path)
        ApiCall(:file_drop_reference, :POST, :files, true, :drop_reference, json=[:fid])
        ApiCall(:file_trash, :POST, :files, true, :trash)
        ApiCall(:file_restore, :POST, :files, true, :restore)
        ApiCall(:file_permanent_delete, :POST, :files, true, :permanent_delete, data_out=false)
        ApiCall(:file_lookup, :GET, :files, false, :lookup, required=[:path], optional=[:parent, :user, :exists])
        ApiCall(:file_star, :POST, :files, true, :star)
        ApiCall(:file_remove_star, :DELETE, :files, true, :star, data_out=false)
        ApiCall(:file_sources, :GET, :files, true, :sources)

        # grids
        ApiCall(:grid_create, :POST, :grids, false, required_json=[:data], json=file_writeable_metadata)
        # ApiCall(:grid_upload)  # failing
        ApiCall(:grid_row, :POST, :grids, true, :row, required_json=[:rows], data_out=false)
        ApiCall(:grid_get_col, :GET, :grids, true, :col, uids=true)
        ApiCall(:grid_put_col, :PUT, :grids, true, :col, uids=true, required_pre_json=[:cols])
        ApiCall(:grid_post_col, :POST, :grids, true, :col, required_pre_json=[:cols])
        ApiCall(:grid_retrieve, :GET, :grids, true)
        ApiCall(:grid_content, :GET, :grids, true, :content)
        ApiCall(:grid_destroy, :DELETE, :grids, true, data_out=false)
        ApiCall(:grid_partial_update, :PATCH, :grids, true, json=file_writeable_metadata)
        ApiCall(:grid_update, :PUT, :grids, true, json=file_writeable_metadata)
        ApiCall(:grid_drop_reference, :POST, :grids, true, :drop_reference, json=[:fid])
        ApiCall(:grid_trash, :POST, :grids, true, :trash)
        ApiCall(:grid_restore, :POST, :grids, true, :restore)
        ApiCall(:grid_permanent_delete, :POST, :grids, true, :permanent_delete, data_out=false)
        ApiCall(:grid_lookup, :GET, :grids, false, :lookup, required=[:path], optional=[:parent, :user, :exists])

        # plots
        ApiCall(:plot_list, :GET, :plots, false, optional=[:order_by, :min_quality, :max_quality])
        ApiCall(:plot_feed, :GET, :plots, false, :feed)
        ApiCall(:plot_create, :POST, :plots, false; required_json=[:figure], json=file_writeable_metadata)
        ApiCall(:plot_detail, :GET, :plots, true)
        ApiCall(:plot_content, :GET, :plots, true, :content, optional=[:inline_data, :map_data])
        ApiCall(:plot_update, :PUT, :plots, true, json=vcat(file_writeable_metadata, :figure))
        ApiCall(:plot_partial_update, :PATCH, :plots, true, json=file_writeable_metadata)

        # extras
        ApiCall(:extra_create, :POST, :extras, false, required_json=[:referencers], json=[:filename, :content])
        ApiCall(:extra_content, :POST, :extras, true, :content)
        ApiCall(:extra_partial_update, :PATCH, :extras, true, json=[:filename, :content])
        ApiCall(:extra_delete, :DELETE, :extras, true, data_out=false)
        ApiCall(:extra_detail, :GET, :extras, true)

        # folders
        ApiCall(:folder_create, :POST, :folders, false, required_json=[:path], json=[:parent])
        ApiCall(:folder_detail, :GET, :folders, true)
        ApiCall(:folder_home, :GET, :folders, false, :home, optional=[:user])
        ApiCall(:folder_shared, :GET, :folders, false, :shared)
        ApiCall(:folder_starred, :GET, :folders, false, :starred)
        ApiCall(:folder_trashed, :GET, :folders, false, :trashed)
        ApiCall(:folder_all, :GET, :folders, false, :all, optional=[:user, :filetype, :order_by])
        ApiCall(:folder_trash, :POST, :folders, true, :trash)
        ApiCall(:folder_restore, :POST, :folders, true, :restore)
        ApiCall(:folder_permanent_delete, :POST, :folders, true, :permanent_delete)

        # images
        ApiCall(:image_generate, :POST, :images, false, required_json=[:figure], json=[:width, :height, :format, :scale, :encoded])

        # comments
        ApiCall(:comment_create, :POST, :comments, false, required_json=[:fid, :comment])
        ApiCall(:comment_delete, :DELETE, :comments, true)

        # dashboards
        ApiCall(:dashboard_create, :POST, :dashboards, false, required_json=[:content])
        ApiCall(:dashboard_list, :GET, :dashboards, false)
        ApiCall(:dashboard_retrieve, :GET, :dashboards, true)
        ApiCall(:dashboard_update, :PUT, :dashboards, true, json=dashboard_writeable_metadata)
        ApiCall(:dashboard_partial_update, :PATCH, :dashboards, true, json=dashboard_writeable_metadata)
        ApiCall(:dashboard_trash, :POST, :dashboards, true, :trash)
        ApiCall(:dashboard_permanent_delete, :DELETE, :dashboards, true, :permanent_delete, data_out=false)
        ApiCall(:dashboard_schema, :GET, :dashboards, false, :schema)

        # plot-schema
        ApiCall(:plot_schema_get, :GET, Symbol("plot-schema"), required=[:sha1])
    ]
    eval(make_method(_api))
end

# --------------------- #
# Convenience functions #
# --------------------- #

function try_me(func, args...; kwargs...)
    try
        func(args...; kwargs...)
    catch
        return nothing
    end
end

try_lookup(path; kwargs...) = try_me(file_lookup, path; kwargs...)

function grid_overwrite!(cols::AbstractDict; fid::String="", path::String="")
    !isempty(fid) && !isempty(path) && error("Can't pass both fid and path")
    if !isempty(fid)
        grid_info = try_me(grid_retrieve, fid)
    elseif !isempty(path)
        grid_info = try_me(grid_lookup, path)
    else
        error("must pass one of fid or path")
    end

    grid_info == nothing && error("can't overwrite a grid that doesn't exit")
    grid_overwrite!(grid_info, cols)
end

function grid_overwrite!(grid_info::AbstractDict, cols::AbstractDict)
    col_name_uid = Dict()
    for col in grid_info["cols"]
        col_name_uid[col["name"]] = col["uid"]
    end

    uids_to_replace = []
    data_to_replace = []
    to_create = []
    out = Dict()
    for (name, data) in cols
        pop!(data, "temp", nothing)
        if haskey(col_name_uid, name)
            # prep this column for replacement
            push!(uids_to_replace, col_name_uid[name])
            push!(data_to_replace, data)
            out[name] = col_name_uid[name]
        else
            # will be adding this column
            push!(to_create, Dict("name" => name, "data" => data["data"]))
        end
    end

    if length(to_create) > 0
        # need to create some columns
        res_post = grid_post_col(grid_info["fid"], to_create)
        for res_col in res_post["cols"]
            out[res_col["name"]] = res_col["uid"]
        end
    end

    if length(data_to_replace) > 0
        res_put = grid_put_col(grid_info["fid"], data_to_replace, uids_to_replace...)
        # we've already added these uids to the output...
    end

    out
end

"""
    grid_overwrite(cols::AbstractDict; fid::String="", path::String="")

Replace the data in the grid assocaited with fid `fid` or at the path `path`
with data in `cols`. `cols` should be an AbstractDict mapping from column names
to column data. The output of this function is an AbstractDict mapping from
column names to column uids in the updated grid.

There are three possible scenarios for the data:

1. The column appears both in the grid and in `cols`. In this case the data in
   that column of the grid will be updated to match the data in `cols`
2. The column appears only in `cols`. In this case a new column will be created
   in the grid
3. The column appears only in the grid. Nothing happens...

NOTE: only one of `fid` and `path` can be passed

"""
grid_overwrite!


text:
export ProbCircuit, StructProbCircuit, StructProbLeafNode, StructProbInnerNode,
    StructProbLiteralNode, StructMulNode, StructSumNode, check_parameter_integrity

#####################
# Prob circuits that are structured,
# meaning that each conjunction is associated with a vtree node.
#####################

"Root of the plain structure probabilistic circuit node hierarchy"
abstract type StructProbCircuit <: ProbCircuit end

"A plain structured probabilistic leaf node"
abstract type StructProbLeafNode <: StructProbCircuit end

"A plain structured probabilistic inner node"
abstract type StructProbInnerNode <: StructProbCircuit end

"A plain structured probabilistic literal leaf node, representing the positive or negative literal of its variable"
mutable struct StructProbLiteralNode <: StructProbLeafNode
    literal::Lit
    vtree::Vtree
    data
    counter::UInt32
    StructProbLiteralNode(l,v) = begin
        @assert lit2var(l) ∈ v
        new(l, v, nothing, 0)
    end
end

"A plain structured probabilistic conjunction node"
mutable struct StructMulNode <: StructProbInnerNode
    prime::StructProbCircuit
    sub::StructProbCircuit
    vtree::Vtree
    data
    counter::UInt32
    StructMulNode(p,s,v) = begin
        @assert isinner(v) "Structured conjunctions must respect inner vtree node"
        @assert varsubset_left(vtree(p),v) "$p does not go left in $v"
        @assert varsubset_right(vtree(s),v) "$s does not go right in $v"
        new(p,s, v, nothing, 0)
    end
end

"A plain structured probabilistic disjunction node"
mutable struct StructSumNode <: StructProbInnerNode
    children::Vector{StructProbCircuit}
    log_probs::Vector{Float64}
    vtree::Vtree # could be leaf or inner
    data
    counter::UInt32
    StructSumNode(c, v) =
        new(c, log.(ones(Float64, length(c)) / length(c)), v, nothing, 0)
end

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension
@inline GateType(::Type{<:StructProbLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:StructMulNode}) = ⋀Gate()
@inline GateType(::Type{<:StructSumNode}) = ⋁Gate()

#####################
# methods
#####################

import LogicCircuits: children, vtree, vtree_safe, respects_vtree # make available for extension
@inline children(n::StructSumNode) = n.children
@inline children(n::StructMulNode) = [n.prime,n.sub]

"Get the vtree corresponding to the argument, or nothing if the node has no vtree"
@inline vtree(n::StructProbCircuit) = n.vtree
@inline vtree_safe(n::StructProbInnerNode) = vtree(n)
@inline vtree_safe(n::StructProbLiteralNode) = vtree(n)

# ProbCircuit has a default argument for respects: its root's vtree
respects_vtree(circuit::StructProbCircuit) =
    respects_vtree(circuit, vtree(circuit))

@inline num_parameters_node(n::StructSumNode) = num_children(n)

#####################
# constructors and compilation
#####################

multiply(arguments::Vector{<:StructProbCircuit};
        reuse=nothing, use_vtree=nothing) =
        multiply(arguments...; reuse, use_vtree)

function multiply(a1::StructProbCircuit,
                 a2::StructProbCircuit;
                 reuse=nothing, use_vtree=nothing)
    reuse isa StructMulNode && reuse.prime == a1 && reuse.sub == a2 && return reuse
    !(use_vtree isa Vtree) && (reuse isa StructProbCircuit) &&  (use_vtree = reuse.vtree)
    !(use_vtree isa Vtree) && (use_vtree = find_inode(vtree_safe(a1), vtree_safe(a2)))
    return StructMulNode(a1, a2, use_vtree)
end

function summate(arguments::Vector{<:StructProbCircuit};
                 reuse=nothing, use_vtree=nothing)
    @assert length(arguments) > 0
    reuse isa StructSumNode && reuse.children == arguments && return reuse
    !(use_vtree isa Vtree) && (reuse isa StructProbCircuit) &&  (use_vtree = reuse.vtree)
    !(use_vtree isa Vtree) && (use_vtree = mapreduce(vtree_safe, lca, arguments))
    return StructSumNode(arguments, use_vtree)
end

# claim `StructProbCircuit` as the default `ProbCircuit` implementation that has a vtree

compile(::Type{ProbCircuit}, a1::Union{Vtree, StructLogicCircuit}, args...) =
    compile(StructProbCircuit, a1, args...)

compile(n::StructProbCircuit, args...) =
    compile(typeof(n), root(vtree(n)), args...)

compile(::Type{<:StructProbCircuit}, c::StructLogicCircuit) =
    compile(StructProbCircuit, root(vtree(c)), c)

compile(::Type{<:StructLogicCircuit}, c::StructProbCircuit) =
    compile(StructLogicCircuit, root(vtree(c)), c)

compile(::Type{<:StructProbCircuit}, ::Vtree, ::Bool) =
    error("Probabilistic circuits do not have constant leafs.")

compile(::Type{<:StructProbCircuit}, vtree::Vtree, l::Lit) =
    StructProbLiteralNode(l,find_leaf(lit2var(l),vtree))

function compile(::Type{<:StructProbCircuit}, vtree::Vtree, circuit::LogicCircuit)
    f_con(n) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")
    f_lit(n) = compile(StructProbCircuit, vtree, literal(n))
    f_a(n, cns) = multiply(cns...) # note: this will use the LCA as vtree node
    f_o(n, cns) = summate(cns) # note: this will use the LCA as vtree node
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, StructProbCircuit)
end

function fully_factorized_circuit(::Type{<:ProbCircuit}, vtree::Vtree)
    ff_logic_circuit = fully_factorized_circuit(PlainStructLogicCircuit, vtree)
    compile(StructProbCircuit, vtree, ff_logic_circuit)
end


text:
using ErgodicControl
using ErgodicControlPlots

# Generate the distribution
N = 80
d = Domain([1,1], [100,100])
cov = 0.010 * eye(2)
phi = zeros(100,100,N+1)
for i = 1:N+1
    mui = (.7*(i-1)/N + .15) * ones(2)
    phi[:,:,i] = gaussian(d, mui, cov)
end
ErgodicControl.normalize!(phi, d.cell_size / (N+1))

# Now let's create the ergodic manager in R3
K = 5
em = ErgodicManagerR2T(d, phi, K)

# Set up first trajectory manager
x0 = [0.49,0.01]
h = 0.015
ci = ConstantInitializer([0.0, 0.0])
tm1 = TrajectoryManager(em, x0, h, N, ci)
tm1.xy = [0.49, 0.01]
dynamics!(tm1, SingleIntegrator(2,h))

# second tm is like the first, but different starting point
tm2 = deepcopy(tm1)
tm2.x0= [.79,.99]
tm2.xy= [.79,.99]

# tm3 = deepcopy(tm1)
# tm3.x0 = [.10,.24]
# tm3.xy = [.10,.24]
#
# tm4 = deepcopy(tm1)
# tm4.x0 = [.5,.02]
# tm4.xy = [.5,.02]
#
# tm5 = deepcopy(tm1)
# tm5.x0 = [.23,.056]
# tm5.xy = [.23,.056]
#
# tm6 = deepcopy(tm1)
# tm6.x0 = [.06,.1]
# tm6.xy = [.06,.1]
#
# tm7 = deepcopy(tm1)
# tm7.x0 = [.23,.65]
# tm7.xy = [.23,.65]
#
# tm8 = deepcopy(tm1)
# tm8.x0 = [.7,.11]
# tm8.xy = [.7,.11]
#
# tm9 = deepcopy(tm1)
# tm9.x0 = [.8,.8]
# tm9.xy = [.8,.8]
#
# tm10 = deepcopy(tm1)
# tm10.x0 = [.1,.41]
# tm10.xy = [.1,.41]

# array of trajectory managers
vtm = [tm1, tm2]

# # Generate the trajectories
# ddc = 1e-4
# tm1.barrier_cost = 1
# tm2.barrier_cost = 1
xd, ud = smc_trajectory_st(em, vtm, 0.2)

# # trajectory generation and plotting
# v = true
# xd,ud = smc_trajectory_s(em, vtm, 0.3, verbose=v)

# generating the gif

gif(em, xd, vtm)


text:
using InteractiveDynamics, DynamicalSystems, GLMakie
using OrdinaryDiffEq
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)

ds = Systems.qbh()
# Grid of initial conditions at given energy:
energy(x,y,px,py) = 0.5(px^2 + py^2) + potential(x,y)
potential(x, y) = 0.5(x^2 + y^2) + 0.55/(√2)*x*(3y^2 - x) + 0.1(x^2 + y^2)^2
function generate_ics(E, n)
    ys = range(-4, stop = 4, length = n)
    pys = range(-10, stop = 10, length = n)
    ics = Vector{Vector{Float64}}()
    for y in ys
        V = potential(0.0, y)
        V ≥ E && continue
        for py in pys
            Ky = 0.5*(py^2)
            Ky + V ≥ E && continue
            px = sqrt(2(E - V - Ky))
            ic = [0.0, y, px, py]
            push!(ics, [px, py, 0., y])
        end
    end
    return ics
end

density = 18
tpsos = 1000.0
tλ = 5000
ttr = 200.0
ics = generate_ics(120.0, density)

tinteg = tangent_integrator(ds, 1; diffeq...)
λ = Float64[]; psos = Dataset{2, Float64}[]
trs = Dataset{3, Float64}[]

@time for u in ics
    # compute Lyapunov exponent (using advanced usage)
    reinit!(tinteg, u, orthonormal(4,1))
    push!(λ, lyapunovs(tinteg, tλ, 1, 1.0)[1])
    push!(psos, poincaresos(ds, (3, 0.0), tpsos; u0 = u, idxs = [4, 2]))
    tr = trajectory(ds, ttr, u)[:, SVector(3, 4, 2)]
    push!(trs, tr)
end

# %%
trajectory_highlighter(psos, λ;
markersize = 0.1, nbins = 10, α = 0.01, hname = "λ")

# trajectory_highlighter(trs[1:(length(trs)÷10):end], λ[1:(length(trs)÷10):end];
# nbins = 10, α = 0.001, linewidth = 2.0, hname = "λ",
# transparency = true)


text:
mutable struct Double_DeepQPolicy <: QPolicy
    primaryNetwork      # Neural Network for Deep Q function
    targetNetwork       # Holds target values - call update_target to copy into this
end
Double_DeepQPolicy(primary) = Double_DeepQPolicy(primary, deepcopy(primary))

get_QValues(policy::Double_DeepQPolicy, inputs; primary=true) = primary ? policy.primaryNetwork(inputs) : policy.targetNetwork(inputs)

function get_target(policy::Double_DeepQPolicy, γ, r, done, s′)
    a_batch = argmax(policy.primaryNetwork(s′), dims=1)
    target = dropdims(r .+ γ.*(1.0 .- done).*policy.targetNetwork(s′)[a_batch], dims=1)
end

get_params(policy::Double_DeepQPolicy) = Flux.params(policy.primaryNetwork)

update_target(policy::Double_DeepQPolicy) = Flux.loadparams!(policy.targetNetwork, Flux.params(policy.primaryNetwork))

action(policy::MooreDDQN.Double_DeepQPolicy, r, s::S, A) where {S} = argmax(policy.primaryNetwork(s()))[1]


text:
facts("CLBLAS.L1") do
    context("CLBLAS axpy") do
        for T in [Float32,
                  Float64,
                  # Complex64,
                  Complex128]
            alpha = T(2.0)
            hx = rand(T, 32)
            hy = rand(T, 32)
            x = CLArray(q, hx)
            y = CLArray(q, hy)

            cl.wait(axpy!(alpha, x, y; queue=q))
            axpy!(alpha, hx, hy)

            @fact cl.to_host(y) --> roughly(hy)
        end
    end

    context("CLBLAS scal") do
        for T in [Float32,
                  Float64,
                  # Complex64,
                  Complex128]

            n = 32
            DA = T(2.0)
            hDX = rand(T, 32)
            DX = CLArray(q, hDX)
            incx = 1

            cl.wait(scal!(n, DA, DX, incx; queue=q))
            scal!(n, DA, hDX, incx)

            @fact cl.to_host(DX) --> roughly(hDX)
        end
    end
end


text:
#### Soil model

export SoilModel, SoilParamFunctions, WaterParamFunctions

"""
    AbstractSoilParameterFunctions{FT <: AbstractFloat}
"""
abstract type AbstractSoilParameterFunctions{FT <: AbstractFloat} end

"""
    WaterParamFunctions{FT, TK, TS, TR} <: AbstractSoilParameterFunctions{FT}

Necessary parameters for the soil water model. These can be floating point
 parameters or functions of space (via the auxiliary variable `aux`). Internally,
they will be converted to type FT or altered so that the functions return type FT.

This is not a complete list - the hydraulic parameters necessary for specifying
the van Genuchten or Brooks and Corey functions are stored in the hydraulics model.
# Fields
$(DocStringExtensions.FIELDS)
"""
struct WaterParamFunctions{FT, TK, TS, TR} <: AbstractSoilParameterFunctions{FT}
    "Saturated conductivity. Units of m s-1."
    Ksat::TK
    "Specific storage. Units of m s-1."
    S_s::TS
    "Residual Water Fraction - default is zero; unitless."
    θ_r::TR
end

"""
    HeatParamFunctions{FT, TK, TS, TR} <: AbstractSoilParameterFunctions{FT}

Necessary parameters for the soil water model. These can be floating point
 parameters or functions of space (via the auxiliary variable `aux`). Internally,
they will be converted to type FT or altered so that the functions return type FT.

This is not a complete list - the hydraulic parameters necessary for specifying
the van Genuchten or Brooks and Corey functions are stored in the hydraulics model.
# Fields
$(DocStringExtensions.FIELDS)
"""

function WaterParamFunctions(
    ::Type{FT};
    Ksat::Union{Function, AbstractFloat} = (aux) -> eltype(aux)(0.0),
    S_s::Union{Function, AbstractFloat} = (aux) -> eltype(aux)(1e-3),
    θ_r::Union{Function, AbstractFloat} = (aux) -> eltype(aux)(0.0),
) where {FT}
    fKsat = Ksat isa AbstractFloat ? (aux) -> FT(Ksat) : (aux) -> FT(Ksat(aux))
    fS_s = S_s isa AbstractFloat ? (aux) -> FT(S_s) : (aux) -> FT(S_s(aux))
    fθ = θ_r isa AbstractFloat ? (aux) -> FT(θ_r) : (aux) -> FT(θ_r(aux))
    args = (fKsat, fS_s, fθ)
    return WaterParamFunctions{FT, typeof.(args)...}(args...)
end

"""
    SoilParamFunctions{FT, F1, F2, F3, F4, F5, F6, F7, F8, F9, F10, F11, F12, WP}
               <: AbstractSoilParameterFunctions{FT}

Necessary parameters for the soil model. Heat parameters to be moved to their
own structure in next iteration.
# Fields
$(DocStringExtensions.FIELDS)
"""
struct SoilParamFunctions{
    FT,
    F1,
    F2,
    F3,
    F4,
    F5,
    F6,
    F7,
    F8,
    F9,
    F10,
    F11,
    F12,
    WP,
} <: AbstractSoilParameterFunctions{FT}
    "Aggregate porosity of the soil"
    porosity::F1
    "Volume fraction of gravels, relative to soil solids only; unitless."
    ν_ss_gravel::F2
    "Volume fraction of SOM, relative to soil solids only; unitless."
    ν_ss_om::F3
    "Volume fraction of quartz, relative to soil solids only; unitless."
    ν_ss_quartz::F4
    "Bulk volumetric heat capacity of dry soil. Units of J m-3 K-1."
    ρc_ds::F5
    "Particle density for soil solids. Units of kg m-3"
    ρp::F6
    "Thermal conductivity of the soil solids. Units of W m-1 K-1."
    κ_solid::F7
    "Saturated thermal conductivity for unfrozen soil. Units of W m-1 K-1."
    κ_sat_unfrozen::F8
    "Saturated thermal conductivity for frozen soil. Units of W m-1 K-1."
    κ_sat_frozen::F9
    "Adjustable scale parameter for determining Kersten number in the Balland and Arp formulation; unitless."
    a::F10
    "Adjustable scale parameter for determining Kersten number in the Balland and Arp formulation; unitless."
    b::F11
    "Parameter used in the Balland and Arp formulation for κ_dry; unitless"
    κ_dry_parameter::F12
    "Hydrology parameter functions"
    water::WP
end

function SoilParamFunctions(
    ::Type{FT};
    porosity::FT = FT(NaN),
    ν_ss_gravel::FT = FT(NaN),
    ν_ss_om::FT = FT(NaN),
    ν_ss_quartz::FT = FT(NaN),
    ρc_ds::FT = FT(NaN),
    ρp::FT = FT(NaN),
    κ_solid::FT = FT(NaN),
    κ_sat_unfrozen::FT = FT(NaN),
    κ_sat_frozen::FT = FT(NaN),
    a::FT = FT(0.24),
    b::FT = FT(18.1),
    κ_dry_parameter::FT = FT(0.053),
    water::AbstractSoilParameterFunctions{FT} = WaterParamFunctions(FT;),
) where {FT}
    args = (
        porosity,
        ν_ss_gravel,
        ν_ss_om,
        ν_ss_quartz,
        ρc_ds,
        ρp,
        κ_solid,
        κ_sat_unfrozen,
        κ_sat_frozen,
        a,
        b,
        κ_dry_parameter,
        water,
    )

    return SoilParamFunctions{FT, typeof.(args)...}(args...)
end

"""
    SoilModel{PF, W, H} <: BalanceLaw

A BalanceLaw for soil modeling.
Users may over-ride prescribed default values for each field.

# Usage

    SoilModel(
        param_functions,
        water,
        heat,
    )

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SoilModel{PF, W, H} <: BalanceLaw
    "Soil Parameter Functions"
    param_functions::PF
    "Water model"
    water::W
    "Heat model"
    heat::H
end

function vars_state(soil::SoilModel, st::Prognostic, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::Auxiliary, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::Gradient, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::GradientFlux, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    compute_gradient_argument!(land, soil, soil.heat, transform, state, aux, t)
    compute_gradient_argument!(land, soil, soil.water, transform, state, aux, t)
end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    compute_gradient_flux!(
        land,
        soil,
        soil.water,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )
    compute_gradient_flux!(
        land,
        soil,
        soil.heat,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )

end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    land_nodal_update_auxiliary_state!(land, soil, soil.water, state, aux, t)
    land_nodal_update_auxiliary_state!(land, soil, soil.heat, state, aux, t)
end

function land_init_aux!(
    land::LandModel,
    soil::SoilModel,
    aux::Vars,
    geom::LocalGeometry,
)
    soil_init_aux!(land, soil, soil.water, aux, geom)
    soil_init_aux!(land, soil, soil.heat, aux, geom)
end

"""
    abstract type AbstractSoilComponentModel <: BalanceLaw
"""
abstract type AbstractSoilComponentModel <: BalanceLaw end

## When PrescribedModels are chosen, all balance law functions and boundary
## condition functions should do nothing. Since these models are of super
## type AbstractSoilComponentModel, use AbstractSoilComponentModel in
## argument type. The more specific Water and Heat models will have
## different methods (see soil_water.jl and soil_heat.jl).

vars_state(m::AbstractSoilComponentModel, st::AbstractStateType, FT) = @vars()

function soil_init_aux!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    aux::Vars,
    geom::LocalGeometry,
) end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    state::Vars,
    aux::Vars,
    t::Real,
)

end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)

end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

end


text:
struct ABOBA{TGV, TF<:AbstractFloat, TM} <: SecondOrderNonMetropolisSampler
    ∇V!::TGV
    β::TF
    γ::TF
    M::TM
    Δt::TF
    c₀::TF
    c₁::TF
    sqrtM::TM
end

"""
    ABOBA(∇V!, β, γ, M, Δt)

Set up the ABOBA integrator for inertial Langevin.

### Fields

* ∇V!   - In place gradient of the potential
* β     - Inverse temperature
* γ     - Damping Coefficient
* M     - Mass (either scalar or vector)
* Δt    - Time step
"""
function ABOBA(∇V!::TGV, β::TF, γ::TF, M::TM, Δt::TF) where {TGV, TF<:AbstractFloat,TM}

    c₀ = exp(-Δt * γ);
    c₁ = sqrt((1 - exp(-2*γ*Δt))/β);
    sqrtM = sqrt.(M);
    return ABOBA(∇V!, β, γ, M, Δt, c₀, c₁, sqrtM)
end

mutable struct ABOBAState{Tq, Tx<:AbstractVector{Tq}} <:SecondOrderNonMetropolisSamplerState
    x::Tx
    q_mid::Tq
    p_mid::Tq
    p̂_mid::Tq
    ∇V_mid::Tq
end

function InitState!(x₀, sampler::ABOBA)
    return ABOBAState(x₀, similar(x₀[1]), similar(x₀[1]), similar(x₀[1]), similar(x₀[1]));
end

function InitState(x₀, sampler::ABOBA)
    return ABOBAState(deepcopy(x₀), similar(x₀[1]), similar(x₀[1]),similar(x₀[1]), similar(x₀[1]));
end

function UpdateState!(state::ABOBAState, sampler::ABOBA)

    @. state.q_mid = state.x[1] + 0.5 * sampler.Δt * state.x[2]/sampler.M;
    sampler.∇V!(state.∇V_mid,state.q_mid);
    @. state.p_mid = state.x[2] - 0.5 * sampler.Δt * state.∇V_mid;
    @. state.p̂_mid = sampler.c₀ * state.p_mid + sampler.c₁ * sampler.sqrtM * randn()
    @. state.x[2] = state.p̂_mid - 0.5 * sampler.Δt * state.∇V_mid;
    @. state.x[1] = state.q_mid + 0.5 * sampler.Δt * state.x[2]/sampler.M;

    state
end

max_stars_repo_name: cortner/JuBasicMD


text:
module MultiClass

using CatBoost

train_data = [["summer", 1924, 44], ["summer", 1932, 37], ["winter", 1980, 37],
              ["summer", 2012, 204]]

eval_data = [["winter", 1996, 197], ["winter", 1968, 37], ["summer", 2002, 77],
             ["summer", 1948, 59]]

cat_features = [0]

train_label = ["France", "USA", "USA", "UK"]
eval_label = ["USA", "France", "USA", "UK"]

train_dataset = Pool(; data=train_data, label=train_label, cat_features=cat_features)

eval_dataset = Pool(; data=eval_data, label=eval_label, cat_features=cat_features)

# Initialize CatBoostClassifier
model = CatBoostClassifier(; iterations=10, learning_rate=1, depth=2,
                           loss_function="MultiClass")

# Fit model
fit!(model, train_dataset)

# Get predicted classes
preds_class = predict(model, eval_dataset)

# Get predicted probabilities for each class
preds_proba = predict_proba(model, eval_dataset)

# Get predicted RawFormulaVal
preds_raw = predict(model, eval_dataset; prediction_type="RawFormulaVal")

end # module


text:
"""
    fitonecycle!(learner, nepochs[, lrmax])

Fit `learner` for `nepochs` using a one-cycle learning rate schedule.
"""
function fitonecycle!(learner, nepochs, lrmax = 0.01; div = 25, divfinal = 1e5, pct_start = 0.25)
    FluxTraining.initlearner!(learner, [TrainingPhase()])
    nsteps = length(learner.data.training)
    schedules = (
        LearningRate => onecycle(nepochs * nsteps, lrmax; div = div, divfinal = divfinal, pct_start = pct_start),
        # TODO: add momentum schedule
    )
    oldscheduler = setschedules!(learner, TrainingPhase(), schedules...)

    try
        fit!(learner, nepochs)
    catch e
        rethrow(e)
    finally
        # reset to previous Scheduler
        if !isnothing(oldscheduler)
            FluxTraining.replacecallback!(learner, oldscheduler)
        end
    end
end


text:
using BlochSim, Test, MAT

function testB2c()

    answer = matread("matlabtestdata/testB2c.mat")

    spin = Spin(1, 600, 100, 0)
    sig = mese!(spin, 1000, 50, 1)[1]

    return sig ≈ answer["sig"]

end

function testB2d()

    answer = matread("matlabtestdata/testB2d.mat")

    spin = Spin(1, 600, 100, 0)
    sig = mese!(spin, 1000, 50, 8)

    return sig ≈ vec(answer["sig"])

end

function testB3a()

    answer = matread("matlabtestdata/testB3a.mat")

    gradz = 0.3 # G/cm
    Tg = 3 # ms
    z = π/2 / (GAMMA * gradz * Tg/1000) # cm
    spin = Spin(1, 600, 100, 0, [0,0,z])
    spgr!(spin, 10, 2, π/3, [0,0,gradz], Tg; Δθinc = 0)

    return spin.M ≈ vec(answer["M"])

end

function testB3b()

    answer = matread("matlabtestdata/testB3b.mat")

    gradz = 0.3 # G/cm
    Tg = 3 # ms
    zmax = 4π / (GAMMA * gradz * Tg/1000) # cm
    z = ((1:100)/100 .- 0.5) * zmax
    spins = map(z -> Spin(1, 600, 100, 0, [0,0,z]), z)
    map(spin -> spgr!(spin, 10, 2, π/3, [0,0,gradz], Tg; Δθinc = 0), spins)
    M = zeros(ComplexF64, 3)
    for spin in spins
        M .+= spin.M
    end
    M ./= 100

    return M ≈ vec(answer["M"])

end

function testB3c()

    answer = matread("matlabtestdata/testB3c.mat")

    spin = Spin(1, 600, 100, 0)
    spgr!(spin, 10, 2, π/3)

    return spin.M ≈ vec(answer["M"])

end

function testB5a()

    answer = matread("matlabtestdata/testB5a.mat")

    gradz = 0.3 # G/cm
    Tg = 3 # ms
    zmax = 2π / (GAMMA * gradz * Tg/1000) # cm
    z = (1:100)/100 * zmax
    spins = map(z -> Spin(1, 600, 100, 0, [0,0,z]), z)
    s = map(spin -> spgr!(spin, 10, 2, π/6, [0,0,gradz], Tg, nTR = 99), spins)
    sig = sum(s) / 100

    return sig ≈ answer["sig"]

end

function testB5b()

    answer = matread("matlabtestdata/testB5b.mat")

    α = LinRange(0, π/2, 51) # rad
    gradz = 0.3 # G/cm
    Tg = 3 # ms
    zmax = 2π / (GAMMA * gradz * Tg/1000) # cm
    z = (1:100)/100 * zmax
    srf = zeros(ComplexF64, 51)
    for i = 1:51
        spins = map(z -> Spin(1, 600, 100, 0, [0,0,z]), z)
        s = map(spin -> spgr!(spin, 10, 2, α[i], [0,0,gradz], Tg, nTR = 99), spins)
        srf[i] = sum(s) / 100
    end
    spin = Spin(1, 600, 100, 0)
    sideal = map(α -> spgr!(spin, 10, 2, α), α)

    return srf ≈ vec(answer["sig1"]) && sideal ≈ vec(answer["sig2"])

end

@testset "Sequences" begin

    @testset "MESE" begin

        @test testB2c()
        @test testB2d()

    end

    @testset "SPGR" begin

        @test testB3a()
        @test testB3b()
        @test testB3c()
        @test testB5a()
        @test testB5b()

    end

end


text:
#####
##### utilities
#####

_edf_repr(value::Union{String,Char}) = value
_edf_repr(date::Date) = uppercase(Dates.format(date, dateformat"dd-u-yyyy"))
_edf_repr(date::DateTime) = Dates.format(date, dateformat"dd\.mm\.yyHH\.MM\.SS")

# XXX this is really really hacky and doesn't support use of scientific notation
# where appropriate; keep in mind if you do improve this to support scientific
# notation, that scientific is NOT allowed in EDF annotation onset/duration fields
function _edf_repr(x::Real)
    result = missing
    if isinteger(x)
        str = string(trunc(Int, x))
        if length(str) <= 8
            result = str
        end
    else
        fpart, ipart = modf(x)
        ipart_str = string('-'^signbit(x), Int(abs(ipart))) # handles `-0.0` case
        fpart_str = @sprintf "%.7f" abs(fpart)
        fpart_str = fpart_str[3:end] # remove leading `0.`
        if length(ipart_str) < 7
            result = ipart_str * '.' * fpart_str[1:(7 - length(ipart_str))]
        elseif length(ipart_str) <= 8
            result = ipart_str
        end
    end
    if !ismissing(result)
        if all(c -> c in ('0', '.', '-'), result)
            x == 0 && return result
        else
            return result
        end
    end
    error("failed to fit number into EDF's 8 ASCII character limit: $x")
end

_edf_metadata_repr(::Missing) = 'X'
_edf_metadata_repr(x) = _edf_repr(x)

function _edf_repr(metadata::T) where T<:Union{PatientID,RecordingID}
    header = T <: RecordingID ? String["Startdate"] : String[]
    return join([header; [_edf_metadata_repr(getfield(metadata, name)) for name in fieldnames(T)]], ' ')
end

function edf_write(io::IO, value, byte_limit::Integer)
    edf_value = _edf_repr(value)
    sizeof(edf_value) > byte_limit && error("EDF value exceeded byte limit (of $byte_limit bytes) while writing: $value")
    bytes_written = Base.write(io, edf_value)
    while bytes_written < byte_limit
        bytes_written += Base.write(io, UInt8(' '))
    end
    return bytes_written
end

# NOTE: The fast-path in `Base.write` that uses `unsafe_write` will include alignment
# padding bytes, which is fine for `Int16` but causes `Int24` to write an extra byte
# for each value. To get around this, we'll fall back to a naive implementation when
# the size of the element type doesn't match its aligned size. (See also `read_to!`)
function write_from(io::IO, x::AbstractArray{T}) where {T}
    if sizeof(T) == Base.aligned_sizeof(T)
        return Base.write(io, x)
    else
        n = 0
        for xi in x
            n += Base.write(io, xi)
        end
        return n
    end
end

#####
##### `write_header`
#####

function write_header(io::IO, file::File)
    length(file.signals) <= 9999 || error("EDF does not allow files with more than 9999 signals")
    expected_bytes_written = BYTES_PER_FILE_HEADER + BYTES_PER_SIGNAL_HEADER * length(file.signals)
    bytes_written = 0
    bytes_written += edf_write(io, file.header.version, 8)
    bytes_written += edf_write(io, file.header.patient, 80)
    bytes_written += edf_write(io, file.header.recording, 80)
    bytes_written += edf_write(io, file.header.start, 16)
    bytes_written += edf_write(io, expected_bytes_written, 8)
    bytes_written += edf_write(io, file.header.is_contiguous ? "EDF+C" : "EDF+D", 44)
    bytes_written += edf_write(io, file.header.record_count, 8)
    bytes_written += edf_write(io, file.header.seconds_per_record, 8)
    bytes_written += edf_write(io, length(file.signals), 4)
    signal_headers = SignalHeader.(file.signals)
    for (field_name, byte_limit) in SIGNAL_HEADER_FIELDS
        for signal_header in signal_headers
            field = getfield(signal_header, field_name)
            bytes_written += edf_write(io, field, byte_limit)
        end
    end
    bytes_written += edf_write(io, ' ', 32 * length(file.signals))
    @assert bytes_written == expected_bytes_written
    return bytes_written
end

#####
##### `write_signals`
#####

function write_signals(io::IO, file::File)
    bytes_written = 0
    for record_index in 1:file.header.record_count
        for signal in file.signals
            bytes_written += write_signal_record(io, signal, record_index)
        end
    end
    return bytes_written
end

function write_signal_record(io::IO, signal::Signal, record_index::Int)
    record_start = 1 + (record_index - 1) * signal.header.samples_per_record
    record_stop = record_index * signal.header.samples_per_record
    record_stop = min(record_stop, length(signal.samples))
    return write_from(io, view(signal.samples, record_start:record_stop))
end

function write_signal_record(io::IO, signal::AnnotationsSignal, record_index::Int)
    bytes_written = 0
    for tal in signal.records[record_index]
        bytes_written += write_tal(io, tal)
    end
    bytes_per_record = 2 * signal.samples_per_record
    while bytes_written < bytes_per_record
        bytes_written += Base.write(io, 0x00)
    end
    return bytes_written
end

function write_tal(io::IO, tal::TimestampedAnnotationList)
    bytes_written = 0
    if !signbit(tal.onset_in_seconds) # otherwise, the `-` will already be in number string
        bytes_written += Base.write(io, '+')
    end
    bytes_written += Base.write(io, _edf_repr(tal.onset_in_seconds))
    if tal.duration_in_seconds !== nothing
        bytes_written += Base.write(io, 0x15)
        bytes_written += Base.write(io, _edf_repr(tal.duration_in_seconds))
    end
    if isempty(tal.annotations)
        bytes_written += Base.write(io, 0x14)
        bytes_written += Base.write(io, 0x14)
    else
        for annotation in tal.annotations
            bytes_written += Base.write(io, 0x14)
            bytes_written += Base.write(io, annotation)
            bytes_written += Base.write(io, 0x14)
        end
    end
    bytes_written += Base.write(io, 0x00)
    return bytes_written
end

#####
##### API functions
#####

"""
    EDF.write(io::IO, file::EDF.File)
    EDF.write(path::AbstractString, file::EDF.File)

Write `file` to the given output, returning the number of bytes written.
"""
function write(io::IO, file::File)
    if !file.header.is_contiguous && !any(s -> s isa AnnotationsSignal, file.signals)
        message = """
                  `file.header.is_contiguous` is `false` but `file.signals` does not contain
                  an `AnnotationsSignal`; this is required as per the EDF+ specification for
                  noncontiguous files in order to specify the start time of each data record
                  (see section 2.2.4 for details).
                  """
        throw(ArgumentError(message))
    end
    return write_header(io, file) + write_signals(io, file)
end

write(path::AbstractString, file::File) = Base.open(io -> write(io, file), path, "w")


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # `MeshArrays.jl` defines the `MeshArray` type
#
# Each `MeshArray` contains an array of elementary arrays that (1) are connected at their edges and (2) collectively form a global grid. Overall grid specifications are contained within `gcmgrid` instances, which merely define array sizes and how e.g. grid variables are represented in memory. Importantly, it is only when e.g. grid variables are read from file that sizable memory is allocated.

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Initialize Framework
#
# 1. import `MeshArrays` and plotting tools
# 2. choose e.g. a standard `MITgcm` grid
# 3. download the grid if needed
#

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
using MeshArrays, Plots

pth="../inputs/GRID_LLC90/"
γ=GridSpec("LatLonCap",pth)

http="https://github.com/gaelforget/GRID_LLC90"
!isdir(pth) ? run(`git clone $http $pth`) : nothing;

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Read Example
#
# A `MeshArray` instance, on the chosen grid, can be obtained from a file path (argument #1). Format conversion occur inside the `read` function based on a propotype argument (#2). `read` / `write` calls then convert back and forth between `MeshArray` and `Array` formats.

# + {"slideshow": {"slide_type": "subslide"}}
D=γ.read(γ.path*"Depth.data",MeshArray(γ,Float64))
tmp1=write(D); tmp2=read(tmp1,D)
show(D)

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Subdomain Arrays
#
# The heatmap function as specialized in `../examples/Plots.jl` operates on each `inner-array` sequentially, one after the other.

# + {"slideshow": {"slide_type": "subslide"}}
p=dirname(pathof(MeshArrays))
include(joinpath(p,"../examples/Plots.jl"))
heatmap(D,title="Ocean Depth",clims=(0.,6000.))

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## A `MeshArray` Behaves Like A `Array`
#
# Here are a few examples that would be coded similarly in both cases

# +
size(D)
eltype(D)
view(D,:)

D .* 1.0
D .* D
1000*D
D*1000

D[findall(D .> 300.)] .= NaN
D[findall(D .< 1.)] .= NaN

D[1]=0.0 .+ D[1]
tmp=cos.(D);
# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Let's switch grid now

# +
pth="../inputs/GRID_CS32/"
γ=GridSpec("LatLonCap",pth)

http="https://github.com/gaelforget/GRID_CS32"
!isdir(pth) ? run(`git clone $http $pth`) : nothing

γ=GridSpec("CubeSphere",pth)
D=γ.read(γ.path*"Depth.data",MeshArray(γ,Float32))
show(D)

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## The `exchange` Function
#
# It adds neighboring points at face edges to slightly extend the computational domain as often needed e.g. to compute partial derivatives.
# -

Dexch=exchange(D,4)
show(Dexch)

# + {"slideshow": {"slide_type": "subslide"}}
P=heatmap(D.f[6],title="Ocean Depth (D, Face 6)",lims=(-4,36))
Pexch=heatmap(Dexch.f[6],title="...(Dexch, Face 6)",lims=(0,40))
plot(P,Pexch)

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # Diffusion-Based Smoothing
#
# The unit testing of `MeshArrays.jl` uses the `smooth()` function. Starting from a random noise field, the smoothing efficiency is predictable and can be set via a smoothing scale parameter [(see Weaver and Courtier, 2001)](https://doi.org/10.1002/qj.49712757518).
#
# This example also illustrates the generality of the `MeshArrays` approach, where the same code in `demo2` is readily applicable to any `PeriodicDomain`, `PeriodicChanel`, `CubeSphere`, or `LatLonCap` grid. Here the chosen grid maps onto the `6` faces of a cube with `16*16` points per face, with distances, areas, etc all set to `1.0`.

# + {"slideshow": {"slide_type": "subslide"}}
p=dirname(pathof(MeshArrays))
include(joinpath(p,"../examples/Demos.jl"))
Γ=GridOfOnes("CubeSphere",6,16)
Δ=demo2(Γ);

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# The initial noise field is `D[1]` while the smoothed one is `D[2]`. After `smooth()` has been applied via `demo2()`, the noise field is visibly smoother and more muted.

# + {"cell_style": "split"}
heatmap(Δ[1],title="initial noise",clims=(-0.5,0.5))

# + {"cell_style": "split"}
heatmap(Δ[2],title="smoothed noise",clims=(-0.5,0.5))

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# The computational cost of `smooth()` predictably  increases with the decorrelation scale. For more about how this works, please refer to **Weaver and Courtier, 2001** _Correlation modelling on the sphere using a generalized diffusion equation_ https://doi.org/10.1002/qj.49712757518
# -

Rini=Δ[1]
DXCsm=Δ[3]
DYCsm=Δ[4]
@time Rend=smooth(Rini,DXCsm,DYCsm,Γ);
@time Rend=smooth(Rini,2DXCsm,2DYCsm,Γ);

max_stars_repo_name: abelolivelanna/GlobalOceanNotebooks

max_issues_repo_name: abelolivelanna/GlobalOceanNotebooks

max_forks_repo_name: abelolivelanna/GlobalOceanNotebooks


text:
# ps3000_set_channel.jl

export ps3000_set_channel, PS3000Channel, PS3000Range

function ps3000_set_channel(handle::Int16, channel, enabled, dc, range)
  ps_status = ccall((:ps3000_set_channel, ps3000driver), Int16,
    (Int16, Int16, Int16, Int16, Int16), handle, channel, enabled, dc, range)
  if ps_status == 0
    error("ps3000 Error: one or more of the arguments are out of range")
  end
  return nothing
end

function ps3000_set_channel(handle::Int16, channel, enabled::Bool, dc::Bool, range)
  ps3000_set_channel(handle::Int16, channel, enabled?0x0001:0x0000,
                     dc?0x0001:0x0000, range)
end

baremodule PS3000Channel
  const A = 0
  const B = 1
  const C = 2
  const D = 3
end # module

baremodule PS3000Range
  const mV20 = 1  # +/-20mV
  const mV50 = 2  # +/-50mV
  const mV100 = 3 # +/- 100mV
  const mV200 = 4  # +/- 200mV
  const mV500 = 5  # +/- 500mV
  const V1 = 6 # +/- 1V
  const V2 = 7 # +/- 2V
  const V5 = 8 # +/- 5V
  const V10 = 9 # +/- 10V
  const V20 = 10 # +/- 20V
  const V50 = 11 # +/- 50V
  const V100 = 12 # +/- 100V
  const V200 = 13 # +/- 200V
  const V400 = 14 # +/- 400V
end # module


text:
module TestPursesOperators

using Test

using Purses: Purse, value

@testset "Unary operators" begin
    x = Purse(1.0)

    for op in [
        Base.complex,
        Base.float,
        # Elements
        Base.zero,
        Base.one,
        Base.oneunit,
        # Basic math operations
        Base.:+,
        Base.:-,
        Base.inv,
        Base.sqrt,
        Base.cbrt,
        Base.log,
        Base.exp,
        Base.conj,
        # Trigonomotry
        Base.cos,
        Base.sin,
        Base.tan,
        Base.cot,
        Base.sec,
        Base.csc,
        Base.cosh,
        Base.sinh,
        Base.tanh,
        Base.coth,
        Base.sech,
        Base.csch,
        Base.acos,
        Base.asin,
        Base.atan,
        Base.acot,
        Base.asec,
        Base.acsc,
        Base.acosh,
        Base.asinh,
        Base.atanh,
        Base.acoth,
        Base.asech,
        Base.acsch,
        # Iterator operations
        Base.sum,
        Base.prod,
        Base.length,
        Base.size,
        Base.iterate,
        Base.keys,
        Base.values,
        Base.pairs,
        Base.broadcastable,
        Base.maximum,
        Base.minimum,
        # Linear algebra operations
        Base.adjoint,
        Base.transpose,
    ]
        @test op(x) === op(value(x))
    end
end

@testset "Binary operators" begin
    x = Purse(1.0)

    for op in [
        Base.:+,
        Base.:*,
        Base.:-,
        Base.:/,
        Base.:\,
        Base.:^,
        Base.:(==),
        Base.isapprox,
        Base.isless,
        Base.max,
        Base.min,
    ]
        @test op(x, x) == op(x, value(x)) == op(value(x), x) == op(value(x), value(x))
    end
end

@testset "Iteration protocol" begin
    purse = Purse(1:5, sum)
    @test map(x -> x^2, purse) == map(x -> x^2, value(purse))

    @test iterate(purse) === (1, 1)
    @test iterate(purse, 1) === (2, 2)
    @test iterate(purse, 5) === nothing

    for idx in eachindex(purse)
        @test purse[idx] == value(purse)[idx]
    end

    for (x1, x2) in zip(purse, value(purse))
        @test x1 == x2
    end

    collected_purse = Purse(collect(purse), sum)
    @test_throws ErrorException purse[1] = 2
end

@testset "Function wrapping" begin
    sum_purse = Purse(sum)
    v = [1.0, 4.0, 3.5]
    @test sum_purse(v) == 8.5

    map_curry(f) = xs -> map(f, xs)
    sqrt_curry_purse = Purse(sqrt, map_curry)
    @test map_curry(sqrt_curry_purse)(v) == map(sqrt, v)
end

end # module


text:
module Photosynthesis

using CLIMAParameters
using CLIMAParameters.Planet
using DocStringExtensions
using Parameters
using PkgUtility
using WaterPhysics

# define constants here
struct EarthParameterSet <: AbstractEarthParameterSet end
const EARTH = EarthParameterSet()

T_25(FT=Float64)  = FT(T_freeze(EARTH)) + 25;
RT_25(FT=Float64) = FT(gas_constant()) * T_25(FT);

# export public types
export AbstractCalculationMode,
       AbstractFluoModelParaSet,
       AbstractPhotoModelParaSet,
       AbstractTDParameterSet,
       AirLayer,
       ArrheniusPeakTD,
       ArrheniusTD,
       C3ParaSet,
       C4ParaSet,
       FluoParaSet,
       GCO₂Mode,
       Leaf,
       PCO₂Mode,
       Q10TD

# export parasets
export C3Bernacchi,
       C3CLM,
       C4CLM,
       FluorescenceVanDerTol,
       FluorescenceVanDerTolDrought,
       JmaxTDBernacchi,
       JmaxTDCLM,
       JmaxTDLeuning,
       KcTDBernacchi,
       KcTDCLM,
       KoTDBernacchi,
       KoTDCLM,
       KpepTDBoyd,
       KpepTDCLM,
       Q10TDAngiosperm,
       Q10TDGymnosperm,
       RespirationTDBernacchi,
       RespirationTDCLM,
       VcmaxTDBernacchi,
       VcmaxTDCLM,
       VcmaxTDLeuning,
       VomaxTDBernacchi,
       VpmaxTDBoyd,
       VtoRCollatz,
       ΓStarTDBernacchi,
       ΓStarTDCLM

# export functions
export temperature_correction,
       leaf_ETR!,
       leaf_fluorescence!,
       leaf_jmax!,
       leaf_kc!,
       leaf_km!,
       leaf_ko!,
       leaf_kpep!,
       leaf_photosynthesis!,
       leaf_rd!,
       leaf_temperature_dependence!,
       leaf_vcmax!,
       leaf_vpmax!,
       leaf_Γstar!,
       light_limited_rate!,
       photo_TD_from_set,
       photo_TD_from_val,
       product_limited_rate!,
       rubisco_limited_rate!

include("types/environment.jl" )
include("types/fluorescence.jl")
include("types/leaf.jl"        )
include("types/mode.jl"        )
include("types/temperature.jl" )
include("types/photomodel.jl"  )
include("types/parasets.jl"    )

include("temperature/correction.jl")
include("temperature/dependency.jl")

include("photosynthesis/etr.jl"           )
include("photosynthesis/lightlimited.jl"  )
include("photosynthesis/productlimited.jl")
include("photosynthesis/rubiscolimited.jl")
include("photosynthesis/model.jl"         )

include("fluorescence/fluorescence.jl")

end # module

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
function tet_perm_by_num(n::Int64)
    (n == 1) && (return SVector{4,Int64}(2,4,3,1))
    (n == 2) && (return SVector{4,Int64}(4,1,3,2))
    (n == 3) && (return SVector{4,Int64}(1,4,2,3))
    (n == 4) && (return SVector{4,Int64}(1,2,3,4))
    error("n must be 1:4")
end

get_projections(p::SVector{3,SVector{3,Float64}}, ê) = SVector(dot(p[1], ê), dot(p[2], ê), dot(p[3], ê))
get_projections(p::SVector{4,SVector{3,Float64}}, ê) = SVector(dot(p[1], ê), dot(p[2], ê), dot(p[3], ê), dot(p[4], ê))

function make_obb(p::SVector{N,SVector{3,Float64}}, i_start::Int64) where {N}
    (N == 3) || (N == 4) || error("must have 3 or 4 points")
    ê_1 = normalize(p[mod1(i_start + 1, 3)] - p[i_start])
    ê_3 = triangleNormal(p[1], p[2], p[3])
    ê_2 = cross(ê_3, ê_1)
    proj_1 = get_projections(p, ê_1)
    proj_2 = get_projections(p, ê_2)
    proj_3 = get_projections(p, ê_3)
    p_min = SVector(minimum(proj_1), minimum(proj_2), minimum(proj_3))
    p_max = SVector(maximum(proj_1), maximum(proj_2), maximum(proj_3))
    c, e = minMaxToCenterExtent(p_min, p_max)
    R = hcat(ê_1, ê_2, ê_3)
    return OBB(R * c, e, R)
end

fit_tri_obb(p::SVector{3,SVector{3,Float64}}) = make_obb(p, 1)  # all bounding boxes have the same surface area
function fit_tet_obb(p::SVector{4,SVector{3,Float64}}, ϵ_tet::SVector{4,Float64})
    (0.0 < volume(p)) || error("inverted tet")
    p = sort_so_big_ϵ_last(ϵ_tet, p)
    obb_1 = make_obb(p, 1)
    obb_2 = make_obb(p, 2)
    obb_3 = make_obb(p, 3)
    area_1 = area(obb_1)
    area_2 = area(obb_2)
    area_3 = area(obb_3)
    (max(        area_2, area_3) <= area_1) && (return obb_1)
    (max(area_1,         area_3) <= area_2) && (return obb_2)
    (max(area_1, area_2        ) <= area_3) && (return obb_3)
end


text:
export to_dice_ir_string

using DirectedAcyclicGraphs: DirectedAcyclicGraphs, DAG, num_parents, Inner, Leaf, foldup_aggregate

struct IrMgr <: DiceManager end

abstract type IrNode <: DAG end
abstract type IrLeafNode <: IrNode end
abstract type IrInnerNode <: IrNode end

DirectedAcyclicGraphs.NodeType(::Type{<:IrInnerNode}) = Inner()
DirectedAcyclicGraphs.NodeType(::Type{<:IrLeafNode}) = Leaf()

# Leaf IR nodes
mutable struct IrFlip <: IrLeafNode
    # note: must be mutable to have reference equality
    prob::Float64
end

new_var(::IrMgr, prob) =
    IrFlip(prob)

struct IrConst <: IrLeafNode
    value
end

constant(::IrMgr, c::Bool) =
    IrConst(c)

struct IrId <: IrLeafNode
    id::Symbol
end

# Inner IR nodes
mutable struct IrNegate <: IrInnerNode
    x::IrNode
end

negate(::IrMgr, x) =
    IrNegate(x)

DirectedAcyclicGraphs.children(n::IrNegate) = [n.x]

mutable struct IrConjoin <: IrInnerNode
    x::IrNode
    y::IrNode
end

conjoin(::IrMgr, x, y) =
    IrConjoin(x, y)

DirectedAcyclicGraphs.children(n::IrConjoin) = [n.x, n.y]
mutable struct IrDisjoin <: IrInnerNode
    x::IrNode
    y::IrNode
end

disjoin(::IrMgr, x, y) =
    IrDisjoin(x, y)

DirectedAcyclicGraphs.children(n::IrDisjoin) = [n.x, n.y]

mutable struct IrIte <: IrInnerNode
    cond::IrNode
    then::IrNode
    elze::IrNode
end

ite(::IrMgr, x, y, z) =
    IrIte(x, y, z)

DirectedAcyclicGraphs.children(n::IrIte) = [n.cond, n.then, n.elze]

mutable struct IrLet <: IrInnerNode
    id::IrId
    x::IrNode
    y::IrNode
end

DirectedAcyclicGraphs.children(n::IrLet) = [n.id, n.x, n.y]

rundice(::IrMgr, code::IrNode) =
    rundice(code)

rundice(code::IrNode) =
    rundice(string(code))

infer(mgr::IrMgr, code::IrNode) =
    infer(code)

infer(code::IrNode) =
    infer_ocml(string(code))

###################################

function identify_reuse(root::IrNode)
    num_ids = 0
    identified = Dict{IrNode,IrId}()
    for (irnode, num_par) in num_parents(root)
        if !(irnode isa IrConst) && num_par > 1
            identified[irnode] = IrId(Symbol("x$(num_ids += 1)"))
        end
    end
    identified
end

function treeify(root::IrNode)
    ids = identify_reuse(root)
    lets = Vector{Tuple{IrId, IrNode}}()

    substitute(n_orig,n_sub) = begin
        if haskey(ids, n_orig)
            id = ids[n_orig]
            push!(lets, (id, n_sub))
            id
        else
            n_sub
        end
    end

    f_leaf(n::IrLeafNode) = substitute(n, n)
    f_inner(n::IrInnerNode, cs) = substitute(n, typeof(n)(cs...))
    root_sub = foldup_aggregate(root, f_leaf, f_inner, IrNode)

    foldr(lets; init=root_sub) do le, ir
        IrLet(le[1], le[2], ir)
    end
end

Base.show(io::IO, ::IrMgr, ir::IrNode) =
    print(io, "(", ir, ")")

Base.show(io::IO, ir::IrFlip) =
    print(io, "flip ", ir.prob)

Base.show(io::IO, ir::IrId) =
    print(io, ir.id)

Base.show(io::IO, ir::IrConst) =
    print(io, ir.value)

Base.show(io::IO, ir::IrNegate) =
    print(io, "!(", ir.x, ")")

Base.show(io::IO, ir::IrConjoin) =
    print(io, "(", ir.x, ") && (", ir.y, ")")

Base.show(io::IO, ir::IrDisjoin) =
    print(io, "(", ir.x, ") || (", ir.y, ")")

Base.show(io::IO, ir::IrIte) = begin
    print(io, "if ", ir.cond,
             " then ", ir.then,
             " else ", ir.elze)
end

Base.show(io::IO, ir::IrLet) = begin
    print(io, "let ", ir.id,
             " = ", ir.x,
             " in ", ir.y)
end


text:
@testset ExtendedTestSet "chamfer" begin
    @testset ExtendedTestSet "chamfer" begin
        x = [
            1 1 0 0
            0 1 1 0
            0 1 0 1
            0 1 0 0
        ]
        tfm = Chamfer(x)
        answer = [
            0 0 3 4
            3 0 0 3
            3 0 3 0
            3 0 3 3
        ]
        @test transform(x, tfm) == answer
    end

    @testset ExtendedTestSet "chamfer" begin
        x = [
            1 1 0 0
            0 1 0 0
            0 1 0 0
            0 1 0 0
        ]
        tfm = Chamfer(x)
        answer = [
            0 0 3 6
            3 0 3 6
            3 0 3 6
            3 0 3 6
        ]
        @test transform(x, tfm) == answer
    end

    @testset ExtendedTestSet "chamfer" begin
        x1 = [
            1 1 0 0
            0 1 0 0
            0 1 0 0
            0 1 0 0
        ]
        x = cat(x1, x1; dims=3)
        tfm = Chamfer(x)

        a1 = [
            0 0 3 6
            3 0 3 6
            3 0 3 6
            3 0 3 6
        ]
        answer = cat(a1, a1; dims=3)
        @test transform(x, tfm) == answer
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test, Statistics

datadir = first(splitext(basename(@__FILE__)))

mesh = abaqus_read_mesh(joinpath(datadir, "tet4.inp"))

upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 3*288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
push!(upper.postprocess_fields, "strain", "stress")
@info("upper postprocess: $(upper.postprocess_fields)")

lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
push!(lower.postprocess_fields, "strain", "stress")
@info("lower postprocess: $(lower.postprocess_fields)")

bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", -0.4)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)

# point-wise boundary conditions to prevent free body move
nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1]
nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1]
nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1]
nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1]

bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
# nodes in X2=0 plane
bc_sym13_elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
update!(bc_sym13_elements, "geometry", mesh.nodes)
update!(bc_sym13_elements, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements)

bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
# nodes in X1=0 plane
bc_sym23_elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
update!(bc_sym23_elements, "geometry", mesh.nodes)
update!(bc_sym23_elements, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements)

push!(bc_upper.postprocess_fields, "reaction force")
push!(bc_lower.postprocess_fields, "reaction force")
push!(bc_sym13.postprocess_fields, "reaction force")
push!(bc_sym23.postprocess_fields, "reaction force")

interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
interface.properties.dual_basis = true
interface.properties.contact_state_in_first_iteration = :AUTO
push!(interface.postprocess_fields, "contact pressure")

analysis = Analysis(Nonlinear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower,
              bc_sym13, bc_sym23, interface)

xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
# node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
upper_X = [0.5, 0.5, 0.75]
lower_X = [0.5, 0.5, 0.25]
strain_upper = upper("strain", upper_X, 0.0)
stress_upper = upper("stress", upper_X, 0.0)
strain_lower = lower("strain", lower_X, 0.0)
stress_lower = lower("stress", lower_X, 0.0)
@debug("strain at $upper_X = $strain_upper")
@debug("stress at $upper_X = $stress_upper")
@debug("strain at $lower_X = $strain_lower")
@debug("stress at $lower_X = $stress_lower")


text:
using Distances
using NearestNeighbors
import StateSpaceReconstruction: cembed

@testset "Kraskov (kNN) 3D $i" for i in 1:n_realizations
    # Some random points that we embed
    pts = rand(3, 100)
    E = cembed(pts)

    # How many neighbors to consider?
    k1 = 2
    k2 = 3
    # Which variables go into which marginals?
    target_future = [1]
    target_presentpast = [2]
    source_presentpast = [3]
    conditioned_presentpast = Int[]

    v = TEVars(target_future,
            target_presentpast,
            source_presentpast,
            conditioned_presentpast)

    @testset "On raw points, default metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, target_future,
                                                    target_presentpast,
                                                    source_presentpast,
                                                conditioned_presentpast) |> isfinite
    end

    @testset "On embeddings, default metric" begin
        @test transferentropy_kraskov(E, k1, k2, v) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, target_future,
                                                    target_presentpast,
                                                    source_presentpast,
                                                    conditioned_presentpast) |> isfinite
    end

    @testset "On points, specifying the distance metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
    end

    @testset "On embeddings, specifying the distance metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, v; metric = Chebyshev()) |> isfinite
    end
end

@testset "Kraskov (kNN) 4D $i" for i in 1:n_realizations
    # Some random points that we embed
    pts = rand(4, 200)
    E = cembed(pts)

    # How many neighbors to consider?
    k1 = 2
    k2 = 3
    # Which variables go into which marginals?
    target_future = [1]
    target_presentpast = [2, 3]
    source_presentpast = [4]
    conditioned_presentpast = Int[]

    v = TEVars(target_future,
            target_presentpast,
            source_presentpast,
            conditioned_presentpast)

    @testset "On raw points, default metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, target_future,
                                                    target_presentpast,
                                                    source_presentpast,
                                                conditioned_presentpast) |> isfinite
    end

    @testset "On embeddings, default metric" begin
        @test transferentropy_kraskov(E, k1, k2, v) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, target_future,
                                                    target_presentpast,
                                                    source_presentpast,
                                                    conditioned_presentpast) |> isfinite
    end

    @testset "On points, specifying the distance metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
    end

    @testset "On embeddings, specifying the distance metric" begin
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(pts, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, v; metric = Chebyshev()) |> isfinite
        @test transferentropy_kraskov(E, k1, k2, v; metric = Chebyshev()) |> isfinite
    end
end


text:
# --------- Control ---------
abstract type Control end

get_value(c::Control) = c.value

# --------- RangeControl ---------
mutable struct RangeControl<:Control
    value::Number
    value_min::Number
    value_max::Number
end

RangeControl(value, value_range) = RangeControl(value, value_range[1], value_range[2])

get_value_range(c::RangeControl) = (c.value_min, c.value_max)

copy(rc::RangeControl) = RangeControl(rc.value, rc.value_min, rc.value_max)

# TODO: for trimmer to work, values out of the domain are allowed.
# an optimization method with boundaries must be sought
function set_value!(c::RangeControl, val, allow_out_of_range=false, throw_error=false)
    if allow_out_of_range
        c.value = val
    else
        min, max = get_value_range(c)
        if min <= val <= max
            c.value = val
        else
            if throw_error
                throw(DomainError(val, "val must be between min and max"))
            elseif val < min
                c.value = min
            else
                c.value = max
            end
        end
    end
end

# --------- DiscreteControl ---------
mutable struct DiscreteControl<:Control
    value::Int
    value_choices::Array{Int, 1}
end

get_value_choices(c::DiscreteControl) = c.value_choices

function set_value!(c::DiscreteControl, val)
    if val in get_value_choices(c)
        c.value = val
    else
        throw(DomainError(val, "val not in discrete control options"))
    end
end

copy(dc::DiscreteControl) = DiscreteControl(dc.value, copy(dc.value_choices))


text:
using Mods

door, card = parse.(Int, readlines("day25.txt"))

function part1(door, card)
    # We need to find p such that 7^p = door (20202127)
    start = Mod{20201227}(7)
     p = 1

    while start != door
        start *= 7
        p += 1
    end
    # Niw the solution is card^p (20201227)

    solution = Mod{20201227}(card)^p
    return value(solution)
end

println("Part 1: $(part1(door, card))")


text:
function generate_imagelist(root)
    filenames = sort(TestImages.remotefiles, by=x->lowercase(splitext(x)[1]))

    N = length(filenames)
    names = [splitext(filename)[1] for filename in filenames]
    imgs = [testimage(filename) for filename in filenames]
    colors = [eltype(img) for img in imgs]
    sizes = size.(imgs)
    paths = ones(String,N)

    # Generate and save thumbnails
    mkpath(joinpath(root, "thumbnails"))
    HEIGHT = 200
    for i in 1:N

        if length(sizes[i]) == 2
            # Normal image
            width = round(Int, HEIGHT / sizes[i][1] * sizes[i][2])
            thumbnail = imresize(imgs[i], HEIGHT, width)
            path = joinpath(root, "thumbnails", names[i]*".png")
            paths[i] = path
            save(path, thumbnail)

        elseif length(sizes[i]) == 3
            # 3-dimensional image, will be converted animated gif
            width = round(Int, HEIGHT / sizes[i][1] * sizes[i][2])
            img = collect(imgs[i])
            thumbnail = Array{eltype(img),3}(undef,HEIGHT,width,sizes[i][3])
            for j in 1:sizes[i][3]
                thumbnail[:,:,j] = imresize(img[:,:,j], HEIGHT, width)
            end
            path = joinpath(root, "thumbnails", names[i]*".gif")
            paths[i] = path
            fps = ceil(Int, size(thumbnail, 3)/3) # we want the image loop within 3 seconds
            save(path, thumbnail; fps=fps)

        elseif length(sizes[i]) == 4
            # multi-channel image, will be converted animated gif
            width = round(Int, HEIGHT / sizes[i][1] * sizes[i][2])
            n_frames = sizes[i][3] * sizes[i][4]
            img = reshape(permutedims(collect(imgs[i]),(1,2,4,3)),sizes[i][1],sizes[i][2],n_frames)
            thumbnail = Array{eltype(img),3}(undef,HEIGHT,width,n_frames)
            for j in 1:n_frames
                thumbnail[:,:,j] = imresize(img[:,:,j], HEIGHT, width)
            end
            path = joinpath(root, "thumbnails", names[i]*".gif")
            paths[i] = path
            save(path, thumbnail)

        end
    end

    # Load the metadata
    path_toml = TestImages.image_path("metadata.toml")
    metadata = if isfile(path_toml)
        Dict{String,String}.(TOML.parsefile(path_toml)["images"])
    else
        # backward compat: old artifacts doesn't have "metadata.toml"
        Dict{String,String}[]
    end

    # Generate markdown, including a table of images
    script = """
    # [List of test images](@id imagelist)

    | Image | Name | Color | Size | Note |
    | :---- | :--- | :---- | :--- | :--- |
    """

    for i in 1:N
        filename = filenames[i]
        path_thumbnail = joinpath("thumbnails", basename(paths[i]))
        color = colors[i]
        size = sizes[i]
        name = extract_name(metadata, filename)
        note = extract_note(metadata, filename)
        script *= "| ![]($(path_thumbnail)) | $(name) | `$(color)` | `$(size)` | $(note) |\n"
    end

    write(joinpath(root, "imagelist.md"), script)
end

function extract_name(metadata, filename)
    path_original = "https://raw.githubusercontent.com/JuliaImages/TestImages.jl/images/images/" * filename
    link = "[`$(filename)`]($(path_original))"

    isempty(metadata) && return link

    # Extract name section in the table
    name = splitext(filename)[1]
    j = findfirst(data -> data["name"]==name, metadata)
    if !isnothing(j)
        # If full name is available, full name in bold and its link will be returned
        data = metadata[j]
        if "fullname" in keys(data)
            fullname = data["fullname"]
            return "**$(fullname)**, " * link
        end
    end
    # If full name isn't available, a filename with link to its original image will be returned
    return link
end

function extract_note(metadata, filename)
    isempty(metadata) && return ""
    # Extract note section in the table

    name = splitext(filename)[1]
    j = findfirst(data -> data["name"]==name, metadata)
    note = String[]
    if !isnothing(j)
        data = metadata[j]
        if "url" in keys(data)
            url = data["url"]
            push!(note, "[origin]($(url))")
        end
        if "author" in keys(data)
            author = data["author"]
            push!(note, "by $(author)")
        end
        if "lisence" in keys(data)
            lisence = data["lisence"]
            push!(note, "$(lisence)")
        end
    end
    return join(note, ", ")
end


text:
"""
    TypeM()

Lorem ipsum dolor sit amet.
"""
struct TypeM

end


text:
"""
    augment_batch_dim(input, n)

Repeat each sample in input batch n-times along batch dimension.
This turns arrays of size `(..., B)` into arrays of size `(..., B*n)`.

## Example
```julia-repl
julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> augment_batch_dim(A, 3)
2×6 Matrix{Int64}:
 1  1  1  2  2  2
 3  3  3  4  4  4
```
"""
function augment_batch_dim(input::AbstractArray{T,N}, n) where {T,N}
    return repeat(input; inner=(ntuple(Returns(1), N - 1)..., n))
end

"""
    reduce_augmentation(augmented_input, n)

Reduce augmented input batch by averaging the explanation for each augmented sample.
"""
function reduce_augmentation(input::AbstractArray{T,N}, n) where {T<:AbstractFloat,N}
    # Allocate output array
    in_size = size(input)
    in_size[end] % n != 0 &&
        throw(ArgumentError("Can't reduce augmented batch size of $(in_size[end]) by $n"))
    out_size = (in_size[1:(end - 1)]..., div(in_size[end], n))
    out = similar(input, eltype(input), out_size)

    axs = axes(input, N)
    colons = ntuple(Returns(:), N - 1)
    for (i, ax) in enumerate(first(axs):n:last(axs))
        view(out, colons..., i) .= sum(view(input, colons..., ax:(ax + n - 1)); dims=N) / n
    end
    return out
end

"""
    augment_indices(indices, n)

Strip batch indices and return inidices for batch augmented by n samples.

## Example
```julia-repl
julia> inds = [CartesianIndex(5,1), CartesianIndex(3,2)]
2-element Vector{CartesianIndex{2}}:
 CartesianIndex(5, 1)
 CartesianIndex(3, 2)

julia> augment_indices(inds, 3)
6-element Vector{CartesianIndex{2}}:
 CartesianIndex(5, 1)
 CartesianIndex(5, 2)
 CartesianIndex(5, 3)
 CartesianIndex(3, 4)
 CartesianIndex(3, 5)
 CartesianIndex(3, 6)
```
"""
function augment_indices(inds::Vector{CartesianIndex{N}}, n) where {N}
    indices_wo_batch = [i.I[1:(end - 1)] for i in inds]
    return map(enumerate(repeat(indices_wo_batch; inner=n))) do (i, idx)
        CartesianIndex{N}(idx..., i)
    end
end

"""
    NoiseAugmentation(analyzer, n, [std=1, rng=GLOBAL_RNG])
    NoiseAugmentation(analyzer, n, distribution, [rng=GLOBAL_RNG])

A wrapper around analyzers that augments the input with `n` samples of additive noise sampled from `distribution`.
This input augmentation is then averaged to return an `Explanation`.
"""
struct NoiseAugmentation{A<:AbstractXAIMethod,D<:Sampleable,R<:AbstractRNG} <:
       AbstractXAIMethod
    analyzer::A
    n::Int
    distribution::D
    rng::R
end
function NoiseAugmentation(analyzer, n, distr::Sampleable, rng=GLOBAL_RNG)
    return NoiseAugmentation(analyzer, n, distr::Sampleable, rng)
end
function NoiseAugmentation(analyzer, n, σ::Real=0.1f0, args...)
    return NoiseAugmentation(analyzer, n, Normal(0.0f0, Float32(σ)^2), args...)
end

function (aug::NoiseAugmentation)(input, ns::AbstractNeuronSelector)
    # Regular forward pass of model
    output = aug.analyzer.model(input)
    output_indices = ns(output)

    # Call regular analyzer on augmented batch
    augmented_input = add_noise(augment_batch_dim(input, aug.n), aug.distribution, aug.rng)
    augmented_indices = augment_indices(output_indices, aug.n)
    augmented_expl = aug.analyzer(augmented_input, AugmentationSelector(augmented_indices))

    # Average explanation
    return Explanation(
        reduce_augmentation(augmented_expl.attribution, aug.n),
        output,
        output_indices,
        augmented_expl.analyzer,
        Nothing,
    )
end

function add_noise(A::AbstractArray{T}, distr::Distribution, rng::AbstractRNG) where {T}
    return A + T.(rand(rng, distr, size(A)))
end

"""
    InterpolationAugmentation(model, [n=50])

A wrapper around analyzers that augments the input with `n` steps of linear interpolation
between the input and a reference input (typically `zero(input)`).
The gradients w.r.t. this augmented input are then averaged and multiplied with the
difference between the input and the reference input.
"""
struct InterpolationAugmentation{A<:AbstractXAIMethod} <: AbstractXAIMethod
    analyzer::A
    n::Int
end

function (aug::InterpolationAugmentation)(
    input, ns::AbstractNeuronSelector, input_ref=zero(input)
)
    size(input) != size(input_ref) &&
        throw(ArgumentError("Input reference size doesn't match input size."))

    # Regular forward pass of model
    output = aug.analyzer.model(input)
    output_indices = ns(output)

    # Call regular analyzer on augmented batch
    augmented_input = interpolate_batch(input, input_ref, aug.n)
    augmented_indices = augment_indices(output_indices, aug.n)
    augmented_expl = aug.analyzer(augmented_input, AugmentationSelector(augmented_indices))

    # Average gradients and compute explanation
    expl = (input - input_ref) .* reduce_augmentation(augmented_expl.attribution, aug.n)

    return Explanation(expl, output, output_indices, augmented_expl.analyzer, Nothing)
end

"""
    interpolate_batch(x, x0, nsamples)

Augment batch along batch dimension using linear interpolation between input `x` and a reference input `x0`.

## Example
```julia-repl
julia> x = Float16.(reshape(1:4, 2, 2))
2×2 Matrix{Float16}:
 1.0  3.0
 2.0  4.0

julia> x0 = zero(x)
2×2 Matrix{Float16}:
 0.0  0.0
 0.0  0.0

julia> interpolate_batch(x, x0, 5)
2×10 Matrix{Float16}:
 0.0  0.25  0.5  0.75  1.0  0.0  0.75  1.5  2.25  3.0
 0.0  0.5   1.0  1.5   2.0  0.0  1.0   2.0  3.0   4.0
```
"""
function interpolate_batch(
    x::AbstractArray{T,N}, x0::AbstractArray{T,N}, nsamples
) where {T,N}
    in_size = size(x)
    outs = similar(x, (in_size[1:(end - 1)]..., in_size[end] * nsamples))
    colons = ntuple(Returns(:), N - 1)
    for (i, t) in enumerate(range(zero(T), oneunit(T); length=nsamples))
        outs[colons..., i:nsamples:end] .= x0 + t * (x - x0)
    end
    return outs
end