text:
import Base.length, Base.getindex, Base.iterate, Base.view, Base.size

function view(o::GBIFRecords)
    defined = filter(i -> isassigned(o.occurrences, i), eachindex(o.occurrences))
    view(o.occurrences, defined)
end

function size(o::GBIFRecords)
    length(o.occurrences)
end

function length(o::GBIFRecords)
    length(view(o))
end

function getindex(o::GBIFRecords, i::Int64)
    view(o)[i]
end

function getindex(o::GBIFRecords, r::UnitRange{Int64})
    view(o)[r]
end

function iterate(o::GBIFRecords)
    iterate(collect(view(o)))
end

function iterate(o::GBIFRecords, t::Union{Int64,Nothing})
    iterate(collect(view(o)), t)
end


text:
using Flux.Data

@test cmudict()["CATASTROPHE"] == :[K,AH0,T,AE1,S,T,R,AH0,F,IY0].args

ext: jl


text:
# # `System` library management with `PowerSystemCaseBuilder.jl`

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# [PowerSystemCaseBuilder.jl](github.com/NREL-SIIP/PowerSystemCaseBuilder.jl) provides a utility to manage a library of `System`s.
# The package has utilities to list the available system data and to create instances of each system. By keeping track of which systems
# have been constructed locally, it makes the re-instantiation of systems efficient by utilizing the serialization features and avoiding
# the parsing process for systems that have been previously constructed.

# ### Dependencies
using PowerSystemCaseBuilder

# ### List all systems in library
show_systems()

# ### Systems can be listed by category
# the available categories can be displayed with:
show_categories()

# the available systems can be displayed with:
show_systems(SIIPExampleSystems)

# ### Create a `System`
# *The first time this is run, it will parse csv data. Subsequent executions will rely on serialized data and will execute much faster*
sys = build_system(SIIPExampleSystems, "5_bus_hydro_ed_sys")


text:
# all usage of deprecated functions are warning
# all method extensions of deprecated generic functions are errors

struct DeprecateInfo
    funcname::Any
    sig::Union{Nothing, Array{Any,1}}
    message::Compat.String
    line::Int
end
const deprecates = Dict{Symbol, Vector{DeprecateInfo}}()

function initDeprecateInfo()
    sf = Base.find_source_file("deprecated.jl")
    if sf === nothing
        return
    end
    str = read(sf, String)
    linecharc = cumsum(map(x->length(x)+1, @compat(split(str, "\n", keepempty=true))))

    itr = iterate(str)
    lineabs = 1
    while itr !== nothing
        (_, i) = itr
        problem = false
        ex = nothing
        lineabs = first(searchsorted(linecharc, i))
        try
            itr = Meta.parse(str, i)
        catch
            problem = true
        end
        if !problem
            parseDeprecate(ex, lineabs)
        else
            break
        end
    end
end

# to replace node getter's convenient (easy to remember) macro to actual macro
function argReplace!(e::Expr, pattern, replace)
    for j=1:length(e.args)
        if e.args[j] == pattern
            e.args[j] = replace
        elseif typeof(e.args[j]) == Expr
            argReplace!(e.args[j], pattern, replace)
        end
    end
end

function processOneSig(s, typeHints)
    determineType = (tex) -> begin
        if typeof(tex) == Expr
            tmp = deepcopy(tex)
            for (k,v) in typeHints
                argReplace!(tmp, k, v)
            end
            return tmp
        end

        return tex
    end

    if typeof(s) == Symbol
        return (:normal, :Any)
    elseif Meta.isexpr(s, [:(...),:kw]) && typeof(s.args[1])==Symbol
        return (s.head, :Any)
    elseif Meta.isexpr(s, :(::))
        ditype = determineType(s.args[end])
        return (:normal, ditype)
    elseif Meta.isexpr(s, [:(...), :kw]) && Meta.isexpr(s.args[1], :(::))
        ditype = determineType(s.args[1].args[end])
        return (s.head, ditype)
    else
        println(STDERR, "Lint doesn't understand " * string(s) * " as an argument")
        return (:normal, :Any)
    end
end

function parseDeprecate(ex, lineabs)
    global deprecates

    if Meta.isexpr(ex, :function) || Meta.isexpr(ex, :(=)) && Meta.isexpr(ex.args[1], :call)
        callex = ex.args[1]
        (funcname, sig) = getFuncNameAndSig(callex)
        if !isa(funcname, Symbol)
            # can't deal with complex expressions like Broadcast.func yet
            return
        end
        if in(funcname, [:depwarn, :firstcaller]) || occursin("deprecate", lowercase(string(funcname)))
            # the first two are support functions.
            # Any function declaration that has "deprecate" in the name...
            # well, the user/developer should know what they are in for.
            return
        end
        if Meta.isexpr(ex.args[2], :block)
            blockcontents = filter(x->typeof(x)==Expr && x.head != :line, ex.args[2].args)
            if isempty(blockcontents) || blockcontents[1].head != :call || blockcontents[1].args[1] != :depwarn
                # most likely it's either
                # error. The function already blows up. So why bother with Lint?
                # if the first statement isn't depwarn also no-op
                return
            end
            if !haskey(deprecates, funcname)
                deprecates[funcname] =  DeprecateInfo[]
            end
            msg = string(blockcontents[1].args[2])
            push!(deprecates[funcname], DeprecateInfo(funcname, sig, msg, lineabs))
        end
    elseif Meta.isexpr(ex, :macrocall) && ex.args[1] == Symbol("@deprecate")
        if typeof(ex.args[2]) == Symbol
            funcname = ex.args[2]
            sig = nothing
            msg = string(ex.args[2], " is deprecated. Use ", ex.args[3], " instead")
            if !haskey(deprecates, funcname)
                deprecates[funcname] =  DeprecateInfo[]
            end
            push!(deprecates[funcname], DeprecateInfo(funcname, sig, msg, lineabs))
        elseif Meta.isexpr(ex.args[2], :call)
            old = ex.args[2]
            (funcname, sig) = getFuncNameAndSig(old)
            new = ex.args[3]
            oldcall = sprint(io->Base.show_unquoted(io,old))
            newcall = sprint(io->Base.show_unquoted(io,new))

            if occursin("deprecate", string(funcname))
                return
            end
            if sig == nothing
                return
            end
            if !haskey(deprecates, funcname)
                deprecates[funcname] =  DeprecateInfo[]
            end
            msg = string(oldcall, " is deprecated. Use ", newcall, " instead")
            push!(deprecates[funcname], DeprecateInfo(funcname, sig, msg, lineabs))
        end
    end
end

function getFuncNameAndSig(callex::Expr)
    typeHints = Dict{Symbol, Any}()
    if typeof(callex.args[1])==Symbol || Meta.isexpr(callex.args[1], :(.))
        funcname = callex.args[1]
    elseif Meta.isexpr(callex.args[1], :curly)
        funcname = callex.args[1].args[1]
        for i in 2:length(callex.args[1].args)
            tconstr= callex.args[1].args[i]
            if typeof(tconstr) == Symbol && length(string(tconstr))==1
                typeHints[tconstr] = Any
            elseif Meta.isexpr(tconstr, :(<:))
                typeHints[tconstr.args[1]] = tconstr.args[2]
            end
        end
    elseif Meta.isexpr(callex.args[1], :(::))
        # these kinds of call overloads are hard to understand
        # so for now, just ignore them
        return (nothing, nothing)
    else
        return (nothing, nothing)
    end
    sig = Any[]
    for i in 2:length(callex.args)
        sube = callex.args[i]
        if Meta.isexpr(sube, :parameters)
            continue
        end
        si = processOneSig(sube, typeHints)
        push!(sig, si)
    end
    return (funcname, sig)
end

# as part of Lint's tests
function addDummyDeprecates()
    global deprecates

    mypush! = (sym, sig)->begin
        push!(deprecates[sym], DeprecateInfo(sym, sig, string(sym) * " generic deprecate message", 0))
    end

    deprecates[:testDep1] = DeprecateInfo[]
    mypush!(:testDep1, nothing)
    deprecates[:testDep2] = DeprecateInfo[]
    mypush!(:testDep2, Any[(:normal, :Integer)])
    deprecates[:testDep3] = DeprecateInfo[]
    mypush!(:testDep3, Any[(:normal, :(Array{Any, 1}))])
    deprecates[:testDep4] = DeprecateInfo[]
    mypush!(:testDep4, Any[(:normal, :Integer), (:(...), :Integer)])
    deprecates[:testDep5] = DeprecateInfo[]
    mypush!(:testDep5, Any[(:normal, :(Array{Real, 1}))])
end

# returns nothing, or DeprecateInfo
function functionIsDeprecated(callex::Expr)
    global deprecates
    funcname, sig = getFuncNameAndSig(callex)
    if Meta.isexpr(funcname, :(.))
        funcname = funcname.args[2]
    end
    if typeof(funcname) != Symbol
        return nothing
    end
    if !haskey(deprecates, funcname)
        return nothing
    end
    for di in deprecates[funcname]
        if funcMatchesDeprecateInfo(sig, di)
            return di
        end
    end
    return nothing
end

function funcMatchesDeprecateInfo(sig, di::DeprecateInfo)
    if di.sig == nothing
        return true
    end

    checkIsSubType = (s1, s2) -> begin
        if s2 == :Any
            return true
        end
        if s1 == s2
            return true
        end
        if typeof(s1) == Symbol && typeof(s2) == Symbol
            ret = false
            try
                ret = eval(:($s1 <: $s2))
            catch
            end
            return ret
        elseif typeof(s1) == Expr && typeof(s2) == Expr && s1.head == s2.head && length(s1.args)==length(s2.args)
            for i in 1:length(s1.args)
                if !checkIsSubType(s1.args[i], s2.args[i])
                    return false
                end
            end
            return true
        else
            return false
        end
    end
    for (i,si) in enumerate(sig)
        if i <= length(di.sig)
            ditype = di.sig[i][2]
        elseif !isempty(di.sig) && di.sig[end][1] == :(...)
            ditype = di.sig[end][2]
        else # too many arguments
            return false
        end

        if i == length(sig) && si[1]==:(...) && di.sig[end][1] != :(...)
            return false
        end

        if !checkIsSubType(si[2], ditype)
            return false
        end
    end
    return true
end


text:
using Test

m = page_model()
include("../src/components/CH4cycle.jl")

add_comp!(m, ch4cycle, :ch4cycle)

set_param!(m, :ch4cycle, :e_globalCH4emissions, readpagedata(m,"test/validationdata/e_globalCH4emissions.csv"))
set_param!(m, :ch4cycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))
set_param!(m,:ch4cycle,:y_year_0,2008.)

p = load_parameters(m)
p["y_year"] = Mimi.dim_keys(m.md, :time)
set_leftover_params!(m, p)

#running Model
run(m)

conc=m[:ch4cycle,  :c_CH4concentration]
conc_compare=readpagedata(m,"test/validationdata/c_ch4concentration.csv")

@test conc ≈ conc_compare rtol=1e-4


text:
using Documenter
using MPCC

makedocs(
    sitename = "MPCC.jl",
    format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [MPCC],
    pages = [
             "Home" => "index.md",
             "API" => "api.md",
             "Examples and tutorials" => "tutorial.md",
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/tmigot/MPCC.jl")#
#https://juliadocs.github.io/Documenter.jl/stable/man/hosting/ ?


text:
using PyPlot
using JLD
function plot_lines()
	X = load("../data/sens/solenoid/dJds_s3_K12.jld")
	rho = X["s"]
	dJds = X["dJds"]

	X = load("../data/obj_erg_avg/solenoid/r2_s3.jld")
    s3 = X["s3"]
    J = X["J"]
    fig, ax = subplots(1,1)
    ax.plot(s3, J, ".", ms=10.0)
    ax.xaxis.set_tick_params(labelsize=28)
    ax.yaxis.set_tick_params(labelsize=28)
    ax.set_xlabel(L"$s_2$",fontsize=28)
    ax.set_ylabel(L"$\langle J\rangle$",fontsize=28)
    ax.grid(true)

    eps = 1.5e-2

	X = load("../data/obj_erg_avg/solenoid/r2_s3_sens.jld")
	J = X["J"]
	s3 = X["s3"]

	n = size(dJds)[1]
	J_pts = reshape([J .- eps*dJds J .+ eps*dJds], n,
                    2)'
    s_pts = reshape([s3 .- eps s3 .+ eps], n,
                    2)'

    ax.plot(s_pts[:,1:n-1], J_pts[:,1:n-1], "k",lw=3.0)
    ax.plot(s_pts[:,n], J_pts[:,n], "k",lw=4.0,label=L"$d_s\langle J\rangle$ from S3")

    ax.legend(fontsize=28)

end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
module CspExamples

# -----------------------------------------------
# Add Channel constructors that allow us to use multithreading if it's available!
# These constructors don't exist yet in Julia base, so we'll define them here.
# TODO: delete these once support has been added to Julia base/channel.jl.

# Add a parameter (sticky=false) to the Channel constructors that create a Task, to allow
# that task to be scheduled via the multithreaded scheduler.
function Base.Threads.Channel(func::Function; ctype=Any, csize=0, taskref=nothing, sticky=true)
    chnl = Channel{ctype}(csize)
    task = Task(() -> func(chnl))
    # If multithreading isn't safe in the current VERSION, ignore the user's sticky request
    if VERSION >= v"1.3.0-DEV.554"
        task.sticky = sticky
        bind(chnl, task)
        if sticky
            yield(task) # immediately start it
        else
            schedule(task) # start it on (potentially) another thread
        end
    else
        bind(chnl, task)
        yield(task)
    end
    isa(taskref, Ref{Task}) && (taskref[] = task)
    return chnl
end

function Base.Threads.Channel{T}(f::Function; kwargs...) where T
    return Channel(f; ctype=T, kwargs...)
end

# -----------------------------------------------

"""
    S31_COPY(west::Channel, east::Channel)

3.1 COPY

> "Problem: Write a process X to copy characters output by process west to process, east."

Julia uses Channels to communicate messages between Tasks, which we use to model Hoare's
"processes". This function takes each character out of west and puts it in east,
continuously, until west is closed.

As an addition to the paper's example, we stop when west is closed, otherwise we would just
hang at this point. Note that we don't close east. Returning from this function is signal
enough that west is closed -- at that point the caller can close(east) if desired. Or, the
user can pass this function in the Channel constructor, and when the function returns the
Channel is automatically closed.

Note that if west is never closed, we will never return, which is absolutely fine.
"""
function S31_COPY(west::Channel{Char}, east::Channel{Char})
    for c in west
        put!(east, c)
    end
end

"""
    S32_SQUASH(west::Channel, east::Channel)

3.2 SQUASH

> "Problem: Adapt the previous program [COPY] to replace every pair of
consecutive asterisks "**" by an upward arrow "↑". Assume that the final
character input is not an asterisk."

Following Hoare's example, if we get an asterisk from west, we then take another character
as well. If the second character is also an asterisk, we write the upward arrow; if not, we
put the two characters into east as we would have in COPY.
"""
function S32_SQUASH(west::Channel{Char}, east::Channel{Char})
    for c in west
        if c == '*'
            c1 = take!(west)
            if c1 == '*'
                put!(east, '↑')  # hooray for unicode
            else
                put!(east, c)
                put!(east, c1)
            end
        else
            put!(east, c)
        end
    end
end

"""
    S32_SQUASH_EXT(west::Channel, east::Channel)

Hoare adds a remark to 3.2 SQUASH: "(2) As an exercise, adapt this
process to deal sensibly with input which ends with an odd number of
asterisks."

We do this by wrapping the `take!` in a try-catch, since the naked take! will throw if the
Channel has been closed.
"""
function S32_SQUASH_EXT(west::Channel{Char}, east::Channel{Char})
    for c in west
        if c == '*'
            # TODO: Is it possible to also have a timeout like thomas11 does in Go?
            try
                c1 = take!(west)
                if c1 == '*'
                    put!(east, '↑')
                else
                    put!(east, c)
                    put!(east, c1)
                end
                continue
            catch end
        end
        put!(east, c)
    end
end

"""
    S33_DISASSEMBLE(cardfile::Channel{String}, X::Channel{Char})

3.3 DISASSEMBLE

> "Problem: to read cards from a cardfile and output to process X the stream of characters
they contain. An extra space should be inserted at the end of each card."

NOTE: I have typed the output Channel as `Channel{>:Char}`, meaning we don't care what the
element type of Channel is, as long as it includes Chars. This gives the caller more
flexibility. (e.g. they can pass a Channel{Any} if they want to.)
"""
function S33_DISASSEMBLE(cardfile::Channel{String}, X::Channel{>:Char})
    for cardimage in cardfile
        for c in cardimage
            put!(X, c)
        end
        put!(X, ' ')
    end
end

"""
    S34_ASSEMBLE(X::Channel{Char}, lineprinter::Channel{>:String}, linelength=125)

3.4 ASSEMBLE

> "Problem: To read a stream of characters from process X and print them in lines of 125
characters on a lineprinter. The last line should be completed with spaces if necessary."

NOTE: I've made the linelength a parameter, so that this method can split the input into
lines of whatever length desired. (This makes testing easier!)
"""
function S34_ASSEMBLE(X::Channel{Char}, lineprinter::Channel{>:String}, linelength=125)
    while isopen(X)
        out = Char[]
        for c in Base.Iterators.take(X, linelength)
            push!(out, c)
        end
        out = rpad(String(out), linelength)
        put!(lineprinter, out)
    end
end

"""
    S35_Reformat(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)

3.5 Reformat

> "Problem: Read a sequence of cards of 80 characters each, and print
the characters on a lineprinter at 125 characters per line. Every card
should be followed by an extra space, and the last line should be
completed with spaces if necessary."

This one's fun! We can reuse the existing functions by creating an intermediate Channel and
Task (equivalent to a Process in Hoare's paper) to act as the output and then input.

Note that this function blocks on its call to ASSEMBLE, so it doesn't return until the
reformatting is completed. This keeps with the style of the functions seen so far, which
block internally, allowing the caller to choose whether to run it asynchronously.
"""
function S35_Reformat(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)
    # This Channel constructor creates a channel and spawns a Task (and yields to it). When
    # the Task is completed, the Channel is closed.
    disassembled = Channel(ctype=Char, sticky=false) do ch
        S33_DISASSEMBLE(cardfile, ch)
    end
    S34_ASSEMBLE(disassembled, lineprinter, linelength)
end

"""
    S35_Reformat2(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)

An alternative implementation of S35_Reformat, which makes the concurrency operations more
explicit, and keeps a style more similar to Go.
"""
function S35_Reformat2(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)
    tmp = Channel{Char}(0)
    @async begin
        S33_DISASSEMBLE(cardfile, tmp);
        # Note that we must close(tmp) to signal to ASSEMBLE that it's okay to return,
        # because we've chosen to not have DISASSEMBLE and ASSEMBLE close output channels.
        close(tmp);
    end
    S34_ASSEMBLE(tmp, lineprinter, linelength)
end

"""
    S35_Reformat_non_concurrent(cardfile::Array{String})::Array{String}

As Hoare notes in the paper:
> "This elementary problem [Reformat] is difficult to solve elegantly without coroutines."

This is indeed the case, as illustrated in this best-faith attempt. What makes it tricky is
that you have two separate, but overlapping, loops that end at different times. You cannot
represent them both as real for-loops in straightline code, but you can in concurrent
programming, because a concurrent program grants you two separate Instruction Pointers to
trace through each loop.

Here, without concurrency, we emulate the second loop (printing 125-character lines) by
with an if-statement, and manually re-initializing the loop's variables.
"""
function S35_Reformat_non_concurrent(cardfile::Array{String}, linelength=125)::Array{String}
    outs = String[]
    out = Char[]
    for s in cardfile
        for c in s
            push!(out, c)
            # This is the manual implementation of the "second loop":
            if length(out) == linelength
                # "end of the loop"
                push!(outs, rpad(String(out), linelength))
                out = Char[]  # Re-initialize
            end
        end
        push!(out, ' ')
    end
    # Here we have to duplicate the end of the "second loop" since there are two ways to
    # exit the "loop condition".
    push!(outs, rpad(String(out), linelength))
    return outs
end

"""
    S36_Conway(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)

3.6 Conway's Problem

> "Problem: Adapt the above program to replace every pair of consecutive
asterisks by an upward arrow.

Here again, we just add a call to an intermediate function, just like adding a step in a
Unix pipeline of processes.
"""
function S36_Conway(cardfile::Channel{String}, lineprinter::Channel{>:String}, linelength=125)
    ch1 = Channel(ctype=Char, sticky=false) do ch1; S33_DISASSEMBLE(cardfile, ch1); end
    ch2 = Channel(ctype=Char, sticky=false) do ch2; S32_SQUASH(ch1, ch2); end
    S34_ASSEMBLE(ch2, lineprinter, linelength)
end

end


text:
export get_cad_to_mesh_error

function get_cad_to_mesh_error()
    errors = NamedTuple{(:name, :cad_mass, :mesh_mass, :percent_error),
                        Tuple{String, Float64, Float64, Float64}}[]
    for group in gmsh.model.get_physical_groups()
        gdim, gnum = group
        name = gmsh.model.get_physical_name(gdim, gnum)
        if startswith(name, "Material: ")
            tags = gmsh.model.get_entities_for_physical_group(gdim, gnum)
            mass = 0.0
            for tag in tags
                mass +=  gmsh.model.occ.get_mass(gdim, tag)
            end
            gmsh.plugin.setNumber("MeshVolume", "Dimension", gdim)
            gmsh.plugin.setNumber("MeshVolume", "PhysicalGroup",  gnum)
            vtag = gmsh.plugin.run("MeshVolume")
            data = gmsh.view.get_list_data(vtag)
            mesh_mass = data[3][1][end]
            gmsh.view.remove(vtag)
            err_percent = 100*(mesh_mass - mass)/mass
            push!(errors, (name=name, cad_mass=mass,
                           mesh_mass=mesh_mass, percent_error=err_percent))
        end
    end
    return errors
end

function Base.show(io::IO,
        err::NamedTuple{(:name, :cad_mass, :mesh_mass, :percent_error),
                        Tuple{String, Float64, Float64, Float64}})
    print(io, err.name)
    print(io, ", Error (%): "*string(@sprintf("%.4f", err.percent_error)))
    print(io, ", CAD Mass: "*string(@sprintf("%.4f", err.cad_mass)))
    print(io, ", Mesh Mass: "*string(@sprintf("%.4f", err.mesh_mass)))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Roots, EngEconomics

# Given
MARR = 0.06
A_Init = -1000
A_Year = 1500
B_Init = -2000
B_Year = 2800

# Find the better option using Incremental ROR
ΔInit = B_Init - A_Init
ΔYear = B_Year - A_Year
ΔIRR(x) = ΔInit + presentWorthFactor(x, 1) * ΔYear
ΔIRR_ans = find_zero(ΔIRR, 0.05)

# Because ΔIRR_ans > MARR, the second alternative, B, is better
println("Choose Alt. 2")


text:
@system Plant begin
    DATE ~ hold
    # added for convenience
    DOY(DATE) => Cropbox.Dates.dayofyear(DATE) ~ track::int(u"d")
    # moved from Main
    DOYP: planting_day_of_year => 121 ~ preserve::int(u"d", parameter)

    TMAX: maximum_temperature ~ hold
    TMIN: minimum_tempearture ~ hold
    TMN: mean_temperature ~ hold
    PAR: photosynthetically_active_radiation ~ hold

    SWFAC: soil_water_stress_factor ~ hold

    # distinguish initial values of variables (parameters) from actual variables being updated (i.e. LAI0, N0)
    EMP1: LAI_coeff1 => 0.104 ~ preserve(u"m^2", parameter)
    EMP2: LAI_coeff2 => 0.64 ~ preserve(parameter)
    Fc: canopy_fraction => 0.85 ~ preserve(parameter)
    sla: specific_leaf_area => 0.028 ~ preserve(u"m^2/g")
    INTOT: duration_of_reproductive_stage => 300.0 ~ preserve(u"K*d", parameter)
    LAI0: initial_leaf_area_index => 0.013 ~ preserve(u"m^2/m^2", parameter)
    Lfmax: maximum_leaf_number => 12.0 ~ preserve(parameter)
    N0: initial_leaf_number => 2.0 ~ preserve(parameter)
    nb: LAI_coeff => 5.3 ~ preserve(parameter)
    # correct units for p1 (g/K, not g)
    p1: leaf_senescence_rate => 0.03 ~ preserve(u"g/K", parameter)
    PD: plant_density => 5.0 ~ preserve(u"m^-2", parameter)
    rm: maximum_leaf_appearance_rate => 0.100 ~ preserve(u"d^-1", parameter)
    tb: base_temperature => 10.0 ~ preserve(u"°C", parameter)
    W0: initial_plant_dry_matter => 0.3 ~ preserve(u"g/m^2", parameter)
    Wc0: initial_canopy_dry_matter => 0.045 ~ preserve(u"g/m^2", parameter)
    Wr0: initial_root_dry_matter => 0.255 ~ preserve(u"g/m^2", parameter)

    # handle empirical equation for temperature (nounit)
    PT(nounit(TMIN), nounit(TMAX)): photosynthesis_reduction_factor_for_temp => begin
        1.0 - 0.0025 * ((0.25TMIN + 0.75TMAX) - 26.0)^2
    end ~ track

    # introduce RUE parameter with proper units
    RUE: radiation_use_efficiency => 2.1 ~ preserve(u"g/MJ", parameter)
    Pg(PT, SWFAC, RUE, PAR, PD, Y1, LAI): potential_growth_rate => begin
        PT * SWFAC * RUE * PAR/PD * (1 - ℯ^(-Y1 * LAI))
    end ~ track(u"g/d")

    ROWSPC: row_spacing => 60.0 ~ preserve(u"cm", parameter)
    Y1(ROWSPC, PD): canopy_light_extinction_coeff => begin
        # no need to convert units for ROWSPC
        1.5 - 0.768 * (ROWSPC^2 * PD)^0.1
    end ~ preserve

    # introduce flags for organizing rate calculations depending on development phases
    P(DOY, DOYP): planted => (DOY >= DOYP) ~ flag
    VP(P, N, Lfmax): vegetative_phase => (P && N < Lfmax) ~ flag
    RP(N, Lfmax): reproductive_phase => (N >= Lfmax) ~ flag
    #TODO: development phase code no longer needed
    FL(VP, RP): development_phase_code => (VP ? 1 : RP ? 2 : 0) ~ track

    dN(rm, PT): leaf_number_increase => rm * PT ~ track(u"d^-1", when=VP)
    N(dN): leaf_number ~ accumulate(init=N0, when=VP)

    di(TMN, tb): daily_accumulated_temperature => begin
        tb <= TMN <= 25u"°C" ? TMN - tb : 0
    end ~ track(u"K", when=RP)
    INT(di): accumulated_temperature_during_reproductive_stage ~ accumulate(u"K*d", when=RP)

    _a(EMP2, N, nb) => exp(EMP2 * (N - nb)) ~ track
    dLAI1(SWFAC, PD, EMP1, PT, _a, dN) => begin
        SWFAC * PD * EMP1 * PT * (a/(1+a)) * dN
    end ~ track(u"m^2/m^2/d", when=VP)
    dLAI2(PD, di, p1, sla) => begin
        -PD * di * p1 * sla
    end ~ track(u"m^2/m^2/d", when=RP)
    dLAI(dLAI1, dLAI2): leaf_area_index_increase => dLAI1 + dLAI2 ~ track(u"m^2/m^2/d")
    LAI(dLAI): leaf_area_index ~ accumulate(u"m^2/m^2", init=LAI0)

    E: CH2O_conversion_efficiency => 1.0 ~ preserve(u"g/g", parameter)
    dW(E, Pg, PD) => E * Pg * PD ~ track(u"g/m^2/d")
    W(dW): plant_dry_matter ~ accumulate(u"g/m^2", init=W0, when=P)
    Wc(Fc, dW): canopy_dry_matter => Fc * dW ~ accumulate(u"g/m^2", init=Wc0, when=VP)
    Wr(Fc, dW): root_dry_matter => (1-Fc) * dW ~ accumulate(u"g/m^2", init=Wr0, when=VP)
    Wf(dW): fruit_dry_matter ~ accumulate(u"g/m^2", when=RP)

    endsim(INT, INTOT): end_of_simulation => (INT > INTOT) ~ flag
end


text:
module GraphNeuralNetworks

using Statistics: mean
using LinearAlgebra, Random
using Base: tail
using CUDA
using Flux
using Flux: glorot_uniform, leakyrelu, GRUCell, @functor, batch
using MacroTools: @forward
using NNlib, NNlibCUDA
using NNlib: scatter, gather
using ChainRulesCore
using Reexport

using SparseArrays, Graphs # not needed but if removed Documenter will complain

include("GNNGraphs/GNNGraphs.jl")
@reexport using .GNNGraphs
using .GNNGraphs: COO_T, ADJMAT_T, SPARSE_T,
                  check_num_nodes, check_num_edges

export
    # utils
    reduce_nodes,
    reduce_edges,
    softmax_nodes,
    softmax_edges,
    broadcast_nodes,
    broadcast_edges,
    softmax_edge_neighbors,

    # msgpass
    apply_edges,
    aggregate_neighbors,
    propagate,
    copy_xj,
    copy_xi,
    xi_dot_xj,
    e_mul_xj,

    # layers/basic
    GNNLayer,
    GNNChain,
    WithGraph,
    DotDecoder,

    # layers/conv
    AGNNConv,
    CGConv,
    ChebConv,
    EdgeConv,
    GATConv,
    GATv2Conv,
    GatedGraphConv,
    GCNConv,
    GINConv,
    GraphConv,
    MEGNetConv,
    NNConv,
    ResGatedGraphConv,
    SAGEConv,

    # layers/pool
    GlobalPool,
    GlobalAttentionPool,
    TopKPool,
    topk_index

include("utils.jl")
include("layers/basic.jl")
include("layers/conv.jl")
include("layers/pool.jl")
include("msgpass.jl")
include("deprecations.jl")

end


text:
using DrWatson
@quickactivate "HyperspectraWithNeXL"

using NeXLSpectrum
using Gadfly

fes=loadspectrum(joinpath(datadir(),"exp_raw","FeS_and_FeS2","FeS std.msa"))
fes2=loadspectrum(joinpath(datadir(),"exp_raw","FeS_and_FeS2","FeS2 std.msa"))

set_default_plot_size(10inch,3inch)
p = plot(fes, fes2, klms=[n"Fe K-L3", n"Fe K-M3", n"Fe L3-M5", n"Fe L2-M1", n"S K-L3", n"S K-M3"], xmax=10.0e3)
mkpath(joinpath(papersdir(),"Figures", "Figure 1"))
p |> SVG(joinpath(papersdir(),"Figures", "Figure 1","FeS vs FeS2.svg"), 10inch, 3inch)


text:
module FElements

using StaticArrays
using LinearAlgebra
using MeshCore
import MeshCore: manifdim
using ..RefShapes: RefShapePoint, RefShapeInterval, RefShapeTriangle, RefShapeTetrahedron, RefShapeSquare, RefShapeCube, manifdimv

"""
    FE{RS, SD}

Abstract type of finite element, parameterized by
- `RS`: type of reference shape of the element (triangle, square, ...), and
- `SD`: shape descriptor; refer to the package `MeshCore`.
"""
abstract type FE{RS, SD} end

"""
    FEData{SD}

Type of a finite element data.

Parameterized by
- `SD` = shape descriptor; refer to the package `MeshCore`.
"""
struct FEData{SD}
    sd::SD
    _ndofperfeat::SVector{4, Int64}
end

"""
    shapedesc(fe::FE{RS, SD}) where {RS, SD}

Topological shape description.

Refer to the MeshCore library.
"""
shapedesc(fe::FE{RS, SD}) where {RS, SD} = fe.data.sd

"""
    refshape(fe::FE{RS, SD}) where {RS, SD}

Reference shape.
"""
refshape(fe::FE{RS, SD}) where {RS, SD} = RS

"""
    nfeatofdim(fe::FE{RS, SD}, m) where {RS, SD}

Number of features of manifold dimension `m`. Note that `0 <= m <= 3`.
"""
nfeatofdim(fe::FE{RS, SD}, m) where {RS, SD} = MeshCore.nfeatofdim(fe.data.sd, m)

"""
    feathasdof(fe::FE{RS, SD}, m) where {RS, SD}

How many degrees of freedom are attached to a the feature of manifold dimension
`m`?

Note that `0 <= m <= 3`.
"""
ndofperfeat(fe::FE{RS, SD}, m) where {RS, SD} = fe.data._ndofperfeat[m+1]

"""
    ndofsperel(fe::FE{RS, SD}) where {RS, SD}

Provide the number of degrees of freedom per element.

Enumerate all features of all manifold dimensions, and for each feature multiply
by the number of degrees of freedom per feature. The assumption is that this is
a *scalar* finite element.
"""
function ndofsperel(fe::FE{RS, SD}) where {RS, SD}
    md = manifdim(fe.data.sd)
    n = 0
    for m in 0:1:md
        n = n + nfeatofdim(fe, m) * ndofperfeat(fe, m)
    end
    return n
end

"""
    manifdim(fe::FE{RS, SD}) where {RS, SD}

Get the manifold dimension of the finite element.
"""
manifdim(fe::FE{RS, SD}) where {RS, SD} = MeshCore.manifdim(fe.data.sd)

"""
    Jacobian(::Val{0}, J::T) where {T}

Evaluate the point Jacobian.

- `J` = Jacobian matrix, which isn't really defined well for a 0-manifold.
"""
function Jacobian(::Val{0}, J::T) where {T}
    return 1.0;
end

"""
    Jacobian(::Val{1}, J::T) where {T}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(::Val{1}, J::T) where {T}
    sdim,  ntan = size(J);
    @assert ntan == 1 "Expected number of tangent vectors: 1"
    return norm(J)
end

"""
    Jacobian(::Val{2}, J::T) where {T}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(::Val{2}, J::T) where {T}
    sdim,  ntan = size(J);
    @assert ntan == 2 "Expected number of tangent vectors: 2"
    if sdim == ntan
        @inbounds Jac = (J[1, 1]*J[2, 2] - J[2, 1]*J[1, 2])
        return Jac;# is det(J);% Compute the Jacobian
    else
        return norm(cross(J[:, 1], J[:, 2]));
    end
end

"""
    Jacobian(fe::T, J::FFltMat)::FFlt where {T<:AbstractFESet3Manifold}

Evaluate the volume Jacobian.

`J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(::Val{3}, J::T) where {T}
    sdim,  ntan = size(J);
    @assert (ntan == 3) && (sdim == 3) "Expected number of tangent vectors: 3"
    #Jac = det(J);# Compute the Jacobian
    # The unrolled version
    return (+J[1, 1]*(J[2, 2]*J[3, 3]-J[3, 2]*J[2, 3])
    -J[1, 2]*(J[2, 1]*J[3, 3]-J[2, 3]*J[3, 1])
    +J[1, 3]*(J[2, 1]*J[3, 2]-J[2, 2]*J[3, 1]) )
end

function _jac(locs, conn, gradNpar)
    NBFPE = length(conn)
    j = 1
    J = locs[conn[j]] * gradNpar[j]
    @inbounds for j in 2:NBFPE
        J += locs[conn[j]] * gradNpar[j]
    end
    return J
end

"""
    jacjac(fe::FE{RS, SD}, locs, nodes, gradNpar) where {RS, SD}

Compute the Jacobian matrix and the Jacobian determinant.

This is the generic version suitable for isoparametric elements.
"""
function jacjac(fe::FE{RS, SD}, locs, nodes, gradNpar) where {RS, SD}
    Jac = _jac(locs, nodes, gradNpar)
    return (Jac, Jacobian(manifdimv(refshape(fe)), Jac))
end

"""
    bfun(fe::FESUBT,  param_coords)  where {FESUBT<:FE{RS, SD}}

Evaluate the basis functions for all degrees of freedom of the scalar finite
element at the parametric coordinates. Return a vector of the values.
"""
function bfun(fe::FESUBT,  param_coords)  where {FESUBT<:FE{RS, SD}} where {RS, SD}
end

"""
    bfungradpar(fe::FESUBT,  param_coords)  where {FESUBT<:FE{RS, SD}}

Evaluate the gradients of the basis functions for all degrees of freedom of
the scalar finite element with respect to the parametric coordinates, at the
parametric coordinates given. Return a vector of the gradients.
"""
function bfungradpar(fe::FESUBT,  param_coords)  where {FESUBT<:FE{RS, SD}} where {RS, SD}
end

"""
The elements that have nodal bases can be their own geometry carriers. Elements
that do not have nodal bases, for instance an L2 quadrilateral with a single
basis function associated with the cell, need another element type (in this
case the nodal quadrilateral) to serve as geometry carriers.
"""
_geometrycarrier(fe::T) where {T<:FE{RS, SD}} where {RS, SD} = fe

# L2 ==================================================================
# Linear two-node element. Only nodal basis functions.
struct FEH1_L2_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_L2_TYPE = FEH1_L2_Type{RefShapeInterval, typeof(MeshCore.L2)}

"""
    FEH1_L2()

Construct an H1 finite element of the type L2.

L2 is two-node linear segment element.
"""
FEH1_L2() = FEH1_L2_TYPE(FEData(MeshCore.L2, SVector{4}([1, 0, 0, 0])))

function bfun(fe::FEH1_L2_TYPE,  param_coords)
    return SVector{2}([(1. - param_coords[1]); (1. + param_coords[1])] / 2.0)
end

function bfungradpar(fe::FEH1_L2_TYPE,  param_coords)
    g = reshape([-1.0; +1.0]/2.0, 2, 1)
    return [SVector{1}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T3 ==================================================================
# Linear triangular element. Only nodal basis functions.
struct FEH1_T3_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_T3_TYPE = FEH1_T3_Type{RefShapeTriangle, typeof(MeshCore.T3)}

"""
    FEH1_T3()

Construct an H1 finite element of the type T3.

T3 is 3-node linear triangle element.
"""
FEH1_T3() = FEH1_T3_TYPE(FEData(MeshCore.T3, SVector{4}([1, 0, 0, 0])))

function bfun(fe::FEH1_T3_TYPE,  param_coords)
    return SVector{3}([(1 - param_coords[1] - param_coords[2]); param_coords[1]; param_coords[2]])
end

function bfungradpar(fe::FEH1_T3_TYPE,  param_coords)
    g = [-1. -1.;  +1.  0.;  0. +1.]
    return [SVector{2}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T6 ==================================================================
# Quadratic triangular element. Only nodal basis functions.
struct FEH1_T6_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_T6_TYPE = FEH1_T6_Type{RefShapeTriangle, typeof(MeshCore.T6)}

"""
    FEH1_T6()

Construct an H1 finite element of the type T6.

T6 is 6-node quadratic triangle element.
"""
FEH1_T6() = FEH1_T6_TYPE(FEData(MeshCore.T6, SVector{4}([1, 0, 0, 0])))

function bfun(fe::FEH1_T6_TYPE,  param_coords)
    r=param_coords[1];
    s=param_coords[2];
    t = 1. - r - s;
    val = [t * (t + t - 1);
           r * (r + r - 1);
           s * (s + s - 1);
           4 * r * t;
           4 * r * s;
           4 * s * t];
    return SVector{6}(val)
end

function bfungradpar(fe::FEH1_T6_TYPE,  param_coords)
    r =param_coords[1];
    s =param_coords[2];
    t = 1. - r - s;
    val = [-3+4*r+4*s  -3+4*r+4*s;
           4*r-1  0.0;
           0.0 4*s-1;
           4-8*r-4*s  -4*r;
           4*s  4*r;
           -4*s  4-4*r-8*s];
    return [SVector{2}(val[idx, :])' for idx in 1:size(val, 1)]
end

# Q4 ==================================================================
# Linear quadrilateral element. Only nodal basis functions.
struct FEH1_Q4_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_Q4_TYPE = FEH1_Q4_Type{RefShapeSquare, typeof(MeshCore.Q4)}

"""
    FEH1_Q4()

Construct an H1 finite element of the type Q4.

Q4 is 4-node linear quadrilateral element.
"""
FEH1_Q4() = FEH1_Q4_TYPE(FEData(MeshCore.Q4, SVector{4}([1, 0, 0, 0])))

function bfun(fe::FEH1_Q4_TYPE,  param_coords)
	val = [0.25 * (1. - param_coords[1]) * (1. - param_coords[2]);
	       0.25 * (1. + param_coords[1]) * (1. - param_coords[2]);
	       0.25 * (1. + param_coords[1]) * (1. + param_coords[2]);
	       0.25 * (1. - param_coords[1]) * (1. + param_coords[2])];
    return SVector{4}(val)
end

function bfungradpar(fe::FEH1_Q4_TYPE,  param_coords)
    g =   [-(1. - param_coords[2])*0.25 -(1. - param_coords[1])*0.25;
            (1. - param_coords[2])*0.25 -(1. + param_coords[1])*0.25;
            (1. + param_coords[2])*0.25 (1. + param_coords[1])*0.25;
           -(1. + param_coords[2])*0.25 (1. - param_coords[1])*0.25];
    return [SVector{2}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T3-BUBBLE ==================================================================
# Linear triangular element with a cubic interior bubble.
struct FEH1_T3_BUBBLE_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_T3_BUBBLE_TYPE = FEH1_T3_BUBBLE_Type{RefShapeTriangle, typeof(MeshCore.T3)}

"""
    FEH1_T3_BUBBLE()

Construct an H1 finite element of the type T3 with a cubic bubble.

T3 is 3-node linear triangle element with a cubic bubble. It has the usual
nodal basis functions associated with the vertices, and cubic bubble
associated with the element itself.
"""
FEH1_T3_BUBBLE() = FEH1_T3_BUBBLE_TYPE(FEData(MeshCore.T3, SVector{4}([1, 0, 1, 0])))

function bfun(fe::FEH1_T3_BUBBLE_TYPE,  param_coords)
    xi, eta = param_coords
    return SVector{4}([(1 - xi - eta);
        xi;
        eta;
        (1 - xi - eta) * xi * eta])
end

function bfungradpar(fe::FEH1_T3_BUBBLE_TYPE,  param_coords)
    xi, eta = param_coords
    g = [-1. -1.;
    +1.  0.;
    0. +1.;
    (-xi * eta + (1 - xi - eta) * eta) (-xi * eta + (1 - xi - eta) * xi)]
    return [SVector{2}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T4 ==================================================================
# Linear tetrahedral element. Only nodal basis functions.
struct FEH1_T4_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEH1_T4_TYPE = FEH1_T4_Type{RefShapeTetrahedron, typeof(MeshCore.T4)}

"""
    FEH1_T4()

Construct an H1 finite element of the type T4.

T4 is 4-node linear tetrahedral element.
"""
FEH1_T4() = FEH1_T4_TYPE(FEData(MeshCore.T4, SVector{4}([1, 0, 0, 0])))

function bfun(fe::FEH1_T4_TYPE,  param_coords)
    return SVector{4}([(1 - param_coords[1] - param_coords[2] - param_coords[3]);
            param_coords[1];
            param_coords[2];
            param_coords[3]])
end

function bfungradpar(fe::FEH1_T4_TYPE,  param_coords)
    g = [-1.0 -1.0 -1.0;
         +1.0  0.0  0.0;
          0.0 +1.0  0.0;
          0.0  0.0 +1.0]
    return [SVector{3}(g[idx, :])' for idx in 1:size(g, 1)]
end

# Q4 ==================================================================
# Quadrilateral element. Only interior basis functions.

# To do: Provide more than the constant basis function only. Perhaps by
# additional parameters to the constructor?

struct FEL2_Q4_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEL2_Q4_TYPE = FEL2_Q4_Type{RefShapeSquare, typeof(MeshCore.Q4)}

_geometrycarrier(fe::FEL2_Q4_TYPE) = FEH1_Q4()

"""
    FEL2_Q4()

Construct an L2 finite element of the type Q4.

Q4 is 4-node linear quadrilateral element.
"""
FEL2_Q4() = FEL2_Q4_TYPE(FEData(MeshCore.Q4, SVector{4}([0, 0, 1, 0])))

function bfun(fe::FEL2_Q4_TYPE,  param_coords)
    return SVector{1}(1.0)
end

function bfungradpar(fe::FEL2_Q4_TYPE,  param_coords)
    g =   [0.0 0.0];
    return [SVector{2}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T3 ==================================================================
# Triangular element. Only interior basis functions.
struct FEL2_T3_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEL2_T3_TYPE = FEL2_T3_Type{RefShapeTriangle, typeof(MeshCore.T3)}

_geometrycarrier(fe::FEL2_T3_TYPE) = FEH1_T3()

"""
    FEL2_T3()

Construct an L2 finite element of the type T3.

T3 is 3-node linear triangle element.
"""
FEL2_T3() = FEL2_T3_TYPE(FEData(MeshCore.T3, SVector{4}([0, 0, 1, 0])))

function bfun(fe::FEL2_T3_TYPE,  param_coords)
   return SVector{1}(1.0)
end

function bfungradpar(fe::FEL2_T3_TYPE,  param_coords)
    g =   [0.0 0.0];
    return [SVector{2}(g[idx, :])' for idx in 1:size(g, 1)]
end

# T4 ==================================================================
# Tetrahedral element. Only interior basis functions.

# To do: Provide more than the constant basis function only. Perhaps by
# additional parameters to the constructor?

struct FEL2_T4_Type{RS, SD} <: FE{RS, SD}
    data::FEData{SD}
end
FEL2_T4_TYPE = FEL2_T4_Type{RefShapeTetrahedron, typeof(MeshCore.T4)}

_geometrycarrier(fe::FEL2_T4_TYPE) = FEH1_T4()

"""
    FEL2_T4()

Construct an L2 finite element of the type T4.

T4 is tetrahedral element with only internal degrees of freedom.
"""
FEL2_T4() = FEL2_T4_TYPE(FEData(MeshCore.T4, SVector{4}([0, 0, 0, 1])))

function bfun(fe::FEL2_T4_TYPE,  param_coords)
   return SVector{1}(1.0)
end

function bfungradpar(fe::FEL2_T4_TYPE,  param_coords)
    g =   [0.0 0.0 0.0];
    return [SVector{3}(g[idx, :])' for idx in 1:size(g, 1)]
end

end # module


text:
using Documenter
using skalman

const ROOT = joinpath(@__DIR__, "..")
makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename = "skalman.jl",
    modules = [skalman],
    authors = "Daniel Appelo, and contributors.",
    pages = Any[
        "Home" => "index.md",
        "Types" => "types.md",
        "Methods" => "methods.md",
        "Index" => "function-index.md",
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
 deploydocs(
    repo = "github.com/appelo/skalman.jl.git",
    devbranch = "main"
)


text:
using SMM, DataFrames, CSV, HTTP, Random, LinearAlgebra, Statistics, Optim, DataFramesMeta, GLM

function PS7()
	#question 1
	url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
	df = CSV.read(HTTP.get(url).body, DataFrames)
	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = log.(df.wage);

	function ols_gmm(beta, X, y)
		return (y.-X*beta)'*(y.-X*beta)
	end

	function ols_gmm_with_sigma(beta, X, y)
		g = vcat(y .- X*beta[1:end-1],( (size(y,1)-1)/(size(y,1)-size(X,2)) )*var(y .- X*beta[1:end-1]) .- beta[end]^2)
		return g'*g
	end

	beta_gmm = optimize(b -> ols_gmm(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000))
	println(beta_gmm.minimizer)

	beta_gmm_sigma = optimize(b -> ols_gmm_with_sigma(b, X, y), rand(size(X,2)+1), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000))
	println(beta_gmm_sigma.minimizer)

	#Question 2
	df = dropmissing(df, :occupation)
	df[df.occupation.>7,:occupation] .= 7
	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = df.occupation;

	#from PS2
	function mlogit_mle(alpha, X, y)

			K = size(X,2)
			J = length(unique(y))
			N = length(y)
			bigY = zeros(N,J)
			for j=1:J
				bigY[:,j] = y.==j
			end
			A = [reshape(alpha,K,J-1) zeros(K)]

			P = exp.(X*A) ./ sum.(eachrow(exp.(X*A)))

			return -sum( bigY.*log.(P) )
	end

	alpha_hat= optimize(a -> mlogit_mle(a, X, y), rand(6*size(X,2)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))

	print(alpha_hat.minimizer)

	function mlogit_gmm(α, X, y)
			K = size(X,2)
			J = length(unique(y))
			N = length(y)
			bigY = zeros(N,J)
			for j=1:J
				bigY[:,j] = y.==j
			end
			bigα = [reshape(α,K,J-1) zeros(K)]

			P = exp.(X*bigα) ./ sum.(eachrow(exp.(X*bigα)))

			g = bigY[:] .- P[:]

			J = g'*I*g
			return J
	end

	α_hat_gmm = optimize(a -> mlogit_gmm(a, X, y), alpha_hat.minimizer, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=10_000, show_trace=true, show_every=100))
	print(α_hat_gmm.minimizer)

	td = TwiceDifferentiable(b -> mlogit_gmm(b, X, y), alpha_hat.minimizer; autodiff = :forward)
	α_hat_gmm_random=optimize(td, rand(6*size(X,2)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=10_000, show_trace=true, show_every=100))

	print(α_hat_gmm_random.minimizer)

	#Question 3

	#Question4
	MA = SMM.parallelNormal()
	dc = SMM.history(MA.chains[1])
	dc = dc[dc[:accepted].==true, :]
	println(describe(dc))

	#Question 5
end

PS7()


text:
function compare_magdb(a, b, id, minf, maxf)
    # expects common domains
    @assert nchannels(a) == 1
    @assert nchannels(b) == 1
    @assert domain(a) == domain(b)
    dommask = map(f ->  minf <= f <= maxf, domain(a))
    a, b = SpectrumArray.((a, b))

    a_ = data(amp2db.(abs.(a[dommask])))
    b_ = data(amp2db.(abs.(b[dommask])))
    errs = abs.(a_ .- b_)
    errmask = (!).(isinf.(errs) .| isnan.(errs))

    mae = mean(errs[errmask])
    maxerr = maximum(errs[errmask])
    errs_ = SpectrumArray(errs .|> db2amp, rate(a), domain(a)[dommask], [:error])
    (id = id, mae = mae, maxerr = maxerr, errs = errs_)
end

function compare_phasedeg(a, b, id, minf, maxf)
    # TODO: merge with vituix_spinorama
    # expects common domains
    @assert nchannels(a) == 1
    @assert nchannels(b) == 1
    @assert domain(a) == domain(b)
    dommask = map(f ->  minf <= f <= maxf, domain(a))
    a, b = SpectrumArray.((a, b))

    # normalization and conversion complex -> coordinates
    c2v(c) = (nc = c / abs(c); [real(nc), imag(nc)])
    a_ = data(c2v.(a[dommask]))
    b_ = data(c2v.(b[dommask]))
    errs = acosd.(clamp.(dot.(a_, b_), -1.0, 1.0)) # clamp for numerical stability

    errmask = (!).(isinf.(errs) .| isnan.(errs))
    mae = mean(errs[errmask])
    maxerr = maximum(errs[errmask])
    errs_ = SpectrumArray(errs .|> db2amp, rate(a), domain(a)[dommask], [:error])
    (id = id, mae = mae, maxerr = maxerr, errs = errs_)
end

function test_driver_transform(driver_dir, test_dir, freqrange, comparef,
    threshold_max, threshold_mae, matchFIR;
    template_hor=Regex(raw"VituixCAD_PolarFR hor (.*)\.frd"), template_ver=Regex(raw"VituixCAD_PolarFR ver (.*)\.frd"),
    transform...)

    driver_spinorama = load_spinorama_from_FRDs(driver_dir, template_hor=template_hor, template_ver=template_ver)
    target_spinorama = load_spinorama_from_FRDs(test_dir, template_hor=template_hor, template_ver=template_ver)

    drivers = [Dict(:spinorama => driver_spinorama, :transform => Dict{Symbol,Any}(transform...))]
    test_transform(drivers, target_spinorama,
        freqrange, comparef, threshold_max, threshold_mae, matchFIR)
end

function test_transform(drivers, target_spinorama, freqrange, comparef,
    threshold_max, threshold_mae,
    matchFIR) # true -> interpolate both spinorama to the same type and size as H, false -> covert H to match the target_spinorama

    minf, maxf = toHz(first(freqrange)), toHz(last(freqrange))

    H1 = get(drivers[1][:transform], :H, nothing)
    if !isnothing(H1)
        if matchFIR
            for driver in drivers
                driver[:spinorama] = interpolate(driver[:spinorama], H1)
            end
            target_spinorama = interpolate(target_spinorama, H1)
        else
            for driver in drivers
                driver[:transform][:H] = interpolate(driver[:transform][:H], parent(target_spinorama.hor))
            end
        end
    end

    errs = []
    for ref in angles(target_spinorama.hor)
        X1 = project_vituix(drivers...; lat=0.0, lon=ref, ldist=2.0)
        X2 = parent(target_spinorama.hor[ref])
        push!(errs, comparef(X1, X2, (:hor, ref), minf, maxf))
    end

    # I found to have considerable discrepancies in my measurements of hor-180 and ver-180, horizontal measurement
    # is the default one, so remove 180° from vertical
    for ref in filter(a -> a != 180, angles(target_spinorama.ver))
        lat, lon = normlat(ref)
        X1 = project_vituix(drivers...; lat=lat, lon=lon, ldist=2.0)
        X2 = parent(target_spinorama.ver[ref])
        push!(errs, comparef(X1, X2, (:ver, ref, (lat=lat, lon=lon)), minf, maxf))
    end

    sort!(errs; by=e -> e.maxerr, rev=true)

    worst_maxerr = errs[1].maxerr
    worst_mae = errs[1].mae
    println("\tMAX: $(worst_maxerr), \tMAE: $(worst_mae)\tWORST: ", errs[1].id)
    if worst_maxerr > threshold_max || worst_mae > threshold_mae
        return false
    end
    return true
end

function test_interpolation(input_dir, freqrange, comparef, threshold_max, threshold_mae)
    target_spinorama = load_spinorama_from_FRDs(input_dir)

    minf, maxf = toHz(first(freqrange)), toHz(last(freqrange))

    cond(a) = abs(a % 10) == 0
    select_angles(am, cond) = am[:, map(cond, angles(am))]
    test_spinorama = Spinorama(select_angles(target_spinorama.hor, cond), select_angles(target_spinorama.ver, cond))

    errs = []

    for ref in angles(select_angles(target_spinorama.hor, c -> !cond(c)))
        X1 = project(test_spinorama, 0.0, ref)
        X2 = parent(target_spinorama.hor[ref])
        push!(errs, comparef(X1, X2, (:hor, ref), minf, maxf))
    end

    for ref in angles(select_angles(target_spinorama.ver, c -> !cond(c)))
        X1 = project(test_spinorama, normlat(ref)...)
        X2 = parent(target_spinorama.ver[ref])
        push!(errs, comparef(X1, X2, (:ver, ref), minf, maxf))
    end

    sort!(errs; by=e -> e.mae, rev=true)

    worst_maxerr = errs[1].maxerr
    worst_mae = errs[1].mae
    println("\tMAX: $(worst_maxerr), \tMAE: $(worst_mae)\tWORST: ", NamedTuple{(:id, :maxerr, :mae)}(errs[1]))
    if worst_maxerr > threshold_max || worst_mae > threshold_mae
        return false
    end
    return true
end


text:
"""
    convertPyColumn(pycol::PyObject)

Converts a column of a pandas array to a Julia `NullableArray`.
"""
function convertPyColumn(pycol::PyObject)::NullableArray
    nptype = pycol[:dtype][:kind]
    # list of numpy kinds can be found at
    # http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
    if nptype == "O"
        o = NullableArray([x == nothing ? Nullable() : x for x in pycol])
        return o
    elseif nptype == "M"
        # doubt there are subtle efficiency differences between these, but it's possible
        o = NullableArray{DateTime}(
            [isa(t, PyObject) ? Nullable{DateTime}() : Nullable{DateTime}(t) for t in pycol])
        # (o = NullableArray(Nullable{DateTime}[isa(t, PyObject) ? Nullable() : t for t in
        #                                        pycol]))
        return o
    else
        return NullableArray(pycol[:values])
    end
end
export convertPyColumn

function _inferColumnType(col::NullableArray; max_elements::Integer=100)::DataType
    for i in 1:max_elements
        isnull(col) ? continue : nothing
        thistype = typeof(get(col[i]))
        thistype == PyObject ? continue : nothing
        return thistype
    end
    return Any
end

function _fillNewCol!{T, U}(newcol::NullableArray{T}, oldcol::NullableArray{U})
    for i in 1:length(newcol)
        if isnull(oldcol[i])
            newcol[i] = Nullable()
        else
            newcol[i] = get(oldcol[i])
        end
    end
end

"""
    fixColumnTypes!(df)

Check to see if the dataframe `df` has any columns of type `Any` and attempt to convert
them to the proper types.  This can be called from `convertPyDF` with the option
`fixtypes`.
"""
function fixColumnTypes!(df::DataTable)
    for col in names(df)
        if eltype(eltype(df[col])) ≠ Any continue end
        dtype = _inferColumnType(df[col])
        # I can't find any way around getting into these stupid loops
        newcol = NullableArray{dtype}(length(df[col]))
        _fillNewCol!(newcol, df[col])
        df[col] = newcol
    end
end
export fixColumnTypes!

"""
    convertPyDF(pydf[, fixtypes=true])

Converts a pandas dataframe to a Julia one.

Note that it is difficult to infer the correct types of columns which contain references
to Python objects.  If `fixtypes`, this will attempt to convert any column with eltype
`Any` to the proper type.
"""
function convertPyDF(pydf::PyObject; fixtypes::Bool=true)::DataTable
    df = DataTable()
    for col in pydf[:columns]
        df[Symbol(col)] = convertPyColumn(get(pydf, PyObject, col))
    end
    if fixtypes fixColumnTypes!(df) end
    return df
end
export convertPyDF

"""
    fixPyNones(dtype, a)

Attempts to convert a `NullableArray` to have eltype `dtype` while replacing all Python
`None`s with `Nullable`.
"""
function fixPyNones{T}(::Type{T}, a::NullableArray)
    # exit silently if the array can't possibly hold Nones
    if !((eltype(a) == Any) | (eltype(a) == PyObject)) return end
    pyNone = pybuiltin("None")
    newa = NullableArray(x == pyNone ? Nullable() : convert(T, x) for x in a)
end
export fixPyNones

"""
    fixPyNones!(dtype::DataType, df::DataTable, col::Symbol)

Attempts to convert a column of the dataframe to have eltype `dtype` while replacing all
Python `None`s with `Nullable()`.
"""
function fixPyNones!{T}(::Type{T}, df::DataTable, col::Symbol)
    df[col] = fixPyNones(T, df[col])
    return df
end
export fixPyNones!

"""
    fixPyNones!(df::DataTable)

Attempts to automatically convert all columns of a dataframe to have eltype `Any` while
replacing all Python `None`s with `Nullable()`.
"""
function fixPyNones!(df::DataTable)
    for col in names(df)
        if eltype(df[col]) == PyObject
            fixPyNones!(Any, df, col)
        end
    end
end
export fixPyNones!

"""
    unpickle([dtype,] filename[, fixtypes=true])

Deserializes a python pickle file and returns the object it contains.
Additionally, if `DataTable` is given as the first argument, will
attempt to convert the object to a Julia dataframe with the flag
`fixtypes` (see `convertPyDF`).
"""
function unpickle(filename::String)::PyObject
    f = py"open($filename, 'rb')"
    pyobj = PyPickle[:load](f)
end

function unpickle(::Type{DataTable}, filename::AbstractString;
                  fixtypes::Bool=true)::DataTable
    f = py"open($filename, 'rb')"
    # TODO it may be more efficient to create this from a dictionary than to convert
    pydf = pycall(PyPickle[:load], PyObject, f)
    df = convertPyDF(pydf, fixtypes=fixtypes)
end
export unpickle

"""
    pickle(filename, object)

Converts the provided object to a PyObject and serializes it in
the python pickle format.  If the object provided is a `DataTable`,
this will first convert it to a pandas dataframe.
"""
function pickle(filename::String, object::Any)
    pyobject = PyObject(object)
    f = py"open($filename, 'wb')"
    PyPickle[:dump](pyobject, f)
end

function pickle(filename::String, df::DataTable)
    pydf = pandas(df)
    f = py"open($filename, 'wb')"
    PyPickle[:dump](pydf, f)
end
export pickle

"""
    pandas(df)

Convert a dataframe to a pandas pyobject.
"""
function pandas(df::DataTable)::PyObject
    pydf = pycall(PyPandas[:DataFrame], PyObject)
    for col in names(df)
        pycol = [isnull(x) ? nothing : get(x) for x in df[col]]
        set!(pydf, string(col), pycol)
        # convert datetime to proper numpy type
        if eltype(eltype(df[col])) == DateTime
            set!(pydf, string(col),
                 get(pydf, string(col))[:astype]("<M8[ns]"))
        end
    end
    return pydf
end
export pandas


text:
d = Normal(3,2)
qp = @qp(quantile(d,0.95),0.95)

@testset "fit moments" begin
    dfit = DistributionFits.fit(Normal, moments(d))
    @test dfit ≈ d
    @test_throws ErrorException d = fit(Normal, Moments(3.0))
end;
@testset "fit two quantiles" begin
    qpl = @qp_m(3)
    qpu = @qp_u(5)
    dn = fit(Normal, qpl, qpu)
    @test quantile.(dn, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
    dn = fit(Normal, qpu, qpl) # check resorting arguments
    @test quantile.(dn, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
end;
@testset "fit to quantilepoint and mean" begin
    dfit = fit(Normal, mean(d), qp, Val(:mean))
    @test dfit ≈ d
end;
@testset "fit to quantilepoint and mode" begin
    dfit = fit(Normal, mode(d), qp, Val(:mode))
    @test dfit ≈ d
end;
@testset "fit to quantilepoint and median" begin
    dfit = fit(Normal, median(d), qp, Val(:median))
    @test dfit ≈ d
end;


text:
using Formatting

include("NKOM.jl")
using .NKOM

zs = - collect(Float64, range(0.0, stop=1000.0, length=11))
zs = [0.0, -100, -150, -222, -300, -407, -550, -600, -666.0, -700, -750]
zs = [0.0, -10, -20, -30, -40, -50, -60, -70.0]
Nz = 7

qs = zeros(Float64, length(zs)-1)
qs[:] = [277.334, 277.349, 277.349, 277.349, 277.349, 277.346, 277.345]

Δhs = zs[1:end-1] - zs[2:end]

h_ML = 16.072753
old_q_ML = qs[1]

FLDO = NKOM.getFLDO(zs=zs, h_ML=h_ML, Nz=Nz)

if FLDO != -1
    qs[1:FLDO-1] .= old_q_ML
else
    qs .= old_q_ML
end

println("##### 1. SETTING #####")

println("zs: ", zs)
println("Δhs: ", Δhs)
println("h_ML: ", h_ML)
println(format("old_q_ML: {:f}", old_q_ML))
println("Old qs: ", qs)

println("##### 2. MIX #####")
Δq = NKOM.mixFLDO!(
    qs = qs,
    zs = zs,
    hs = Δhs,
    q_ML = old_q_ML,
    h_ML = h_ML,
    FLDO = FLDO,
)

println(format("Δq: {:f}", Δq))
println("After Mix, qs = ", qs)

println("##### 3. UNMIX #####")
new_q_ML = NKOM.unmixFLDOKeepDiff!(;
    qs = qs,
    zs = zs,
    hs = Δhs,
    h_ML = h_ML,
    FLDO = FLDO,
    Nz = Nz,
    Δq = Δq,
)

if FLDO != -1
    Δq = new_q_ML - qs[FLDO]
else
    Δq = 0.0
end

println("New qs: ", qs)

println(format("new_q_ML: {:f}", new_q_ML))
println(format("new_Δq: {:f}", Δq))

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module HierachicalClusteringsController
using DataFrames, Query
using ..HierachicalClustering
using JSON

export call

function charts(records, labels, normalized)
  results = []
  for dist_type in [:euclidean, :mahalanobis, :cosine]
    for linkage in [:ward, :average]
      result_records = HierachicalClustering.fit(records, labels, dist_type, linkage, normalized)
      results = vcat(results, Dict(:dist_type => dist_type,
                                   :linkage => linkage,
                                   :records => result_records))
    end
  end
  results
end

"""
call(message)

Computes hierachical clustering from records.

Differences for euclidean, mahalanobis and cosine.
Linkages for ward and average.

# Arguments
- `message::Dict`: the arguments dict
  - `"records"`: the array of array of rows
  - `"labels"``: the labels of columns
  - `"normalized"`: needed to normalize values for columns

# Return
the json string of result object

# Examples
```jldoctest
julia> HierachicalClusteringsController.call(Dict("labels" => ["one", "two", "three"], "normalized" => true,  "records" => [[1.0, 2.0, 0.0], [2.0, 1.0, 2.0], [3.0, 0.0, 9.0]]))
"{\"charts\":[{\"dist_type\":\"euclidean\",\"linkage\":\"ward\",\"records\":[[\"one\",\"two\"],\"three\"]},{\"dist_type\":\"euclidean\",\"linkage\":\"average\",\"records\":[[\"one\",\"two\"],\"three\"]},{\"dist_type\":\"mahalanobis\",\"linkage\":\"ward\",\"records\":[[\"one\",\"two\"],\"three\"]},{\"dist_type\":\"mahalanobis\",\"linkage\":\"average\",\"records\":[[\"one\",\"two\"],\"three\"]},{\"dist_type\":\"cosine\",\"linkage\":\"ward\",\"records\":[[\"one\",\"two\"],\"three\"]},{\"dist_type\":\"cosine\",\"linkage\":\"average\",\"records\":[[\"one\",\"two\"],\"three\"]}]}"
```
"""
function call(message)
  records = hcat(convert.(Array{Any,1}, message["records"])...)'
  labels = message["labels"]
  normalized = message["normalized"]

  Dict(:charts => charts(records, labels, normalized)) |> JSON.json
end
end


text:
export EntModel, EntTfModel, EntTpModel

struct IdWeight
    id::Int32
    weight::Float32
end

const WeightedVocabulary = Dict{Symbol,IdWeight}

mutable struct EntModel <: Model
    config::TextConfig
    tokens::WeightedVocabulary
    id2token::Dict{Int32,Symbol}
    m::Int
    n::Int
end

"""
    fit(::Type{EntModel}, model::DistModel; lower=0.0, normalize_words::Function=identity)
    fit(::Type{EntModel}, config::TextConfig, corpus, y; smooth=3, weights=:balance, lower=0.0, nclasses=0)

Fits an EntModel using the already fitted DistModel. It accepts only symbols with a final weight higher or equal than `lower`.
Parameters:
    - `corpus` is text collection
    - `y` the set of associated labels (one-to-one with corpus)
    - `smooth` is a smoothing factor for the histogram.
    - `weights` accepts a list of weights (one per class) to be applied to the histogram
    - `lower` controls the minimum weight to be accepted
    - `nclasses` specifies the number of classes
"""
function fit(::Type{EntModel}, model::DistModel; lower=0.0)
    tokens = WeightedVocabulary()
    nclasses = length(model.sizes)
    maxent = log2(nclasses)

    i = 0
    @inbounds for (token, dist) in model.tokens
        i += 1
        e = 0.0
        pop = sum(dist)

        for j in 1:nclasses
            pj = dist[j] / pop

            if pj > 0.0
                e -= pj * log2(pj)
            end
        end

        e = 1.0 - e / maxent
        if e >= lower
            tokens[token] = IdWeight(i, e)
        end
    end

    id2token = Dict(w.id => t for (t, w) in tokens)
    EntModel(model.config, tokens, id2token, model.m, model.n)
end

function fit(::Type{EntModel}, config::TextConfig, corpus, y; smooth=3, weights=:balance, lower=0.0, nclasses=0)
    dmodel = fit(DistModel, config, corpus, y, nclasses=nclasses, weights=weights, fix=false, smooth=smooth)
    fit(EntModel, dmodel, lower=lower)
end

"""
    prune(model::EntModel, lower)

Prunes the model accepting only those symbols with a weight higher than `lower`

"""
function prune(model::EntModel, lower::Float64)
    tokens = WeightedVocabulary()
    for (t, w) in model.tokens
        if w.weight >= lower
            tokens[t] = IdWeight(w.id, w.weight)
        end
    end

    id2token = Dict(w.id => t for (t, w) in tokens)
    EntModel(model.config, tokens, id2token, model.m, model.n)
end

"""
    prune_select_top(model::EntModel, k::Int)
    prune_select_top(model::EntModel, ratio::AbstractFloat)

Creates a new model preserving only the best `k` terms on `model`; the size can be indicated by the ratio of the database to be kept, i.e., ``0 < ratio < 1``.
"""
function prune_select_top(model::EntModel, k::Int)
    X = sort!(collect(model.tokens), by=x->x[2].weight, rev=true)

    tokens = WeightedVocabulary()
    for i in 1:k
        t, w = X[i]
        tokens[t] = IdWeight(w.id, w.weight)
    end

    id2token = Dict(w.id => t for (t, w) in tokens)
    EntModel(model.config, tokens, id2token, model.m, model.n)
end

prune_select_top(model::EntModel, ratio::AbstractFloat) = prune_select_top(model, floor(Int, length(model.tokens) * ratio))

abstract type EntTfModel end
abstract type EntTpModel end

"""
    vectorize(model::EntModel, data)
    vectorize(model::EntModel, scheme::Type, data)

Computes a weighted bow for the given `data`; the vector is scaled to the unit if `normalize` is true;
`data` is an string or an array of strings. The weighting scheme may be any of `EntTfModel`, `EntTpModel`, or `EntModel`;
the default scheme is EntTpModel
"""
function vectorize(model::EntModel, scheme::Type{T}, data::DataType; normalize=true) where
        T <: Union{EntTfModel,EntTpModel,EntModel} where
        DataType <: Union{AbstractString, AbstractVector{S}} where
        S <: AbstractString

    bow, maxfreq = compute_bow(tokenize(model.config, data))
    vectorize(model, scheme, bow, normalize=normalize)
end

"""
    vectorize(model::EntModel, scheme::Type{T}, bow::BOW; normalize=true)

Computes a weighted bow for the given `data`; the vector is scaled to the unit if `normalize` is true;
`data` is an bag of words. The weighting scheme may be any of `EntTfModel`, `EntTpModel`, or `EntModel`

"""
function vectorize(model::EntModel, scheme::Type{T}, bow::BOW; normalize=true) where T <: Union{EntTfModel,EntTpModel,EntModel} where Tv<:Real
    len = 0

    for v in values(bow)
        len += v
    end

    vec = SVEC()
    for (token, freq) in bow
        t = get(model.tokens, token, nothing)
        if t === nothing
            continue
        end

        w = _weight(scheme, t.weight, freq,  len)
        if w > 1e-6
            vec[t.id] = w
        end
    end

    normalize && normalize!(vec)
    vec
end

vectorize(model::EntModel, data; normalize=true) = vectorize(model, EntTpModel, data, normalize=normalize)

_weight(::Type{EntTpModel}, ent, freq, n) = ent * freq / n
_weight(::Type{EntTfModel}, ent, freq, n) = ent * freq
_weight(::Type{EntModel}, ent, freq, n) = ent

function broadcastable(model::EntModel)
    (model,)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Gdbm"
version = v"1.18.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/gdbm/gdbm-1.18.1.tar.gz", "86e613527e5dba544e73208f42b78b7c022d4fa5a6d5498bf18c8d6f745b91dc"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd gdbm-1.18.1/

if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ../patches/gdbm-1.15-win32.patch
fi

if [[ "${target}" == powerpc64le-* ]] || [[ "${target}" == *-mingw* ]]; then
    # Install `autopoint` and other tools needed by `autoreconf`
    apk add gettext-dev
    # Rebuild the configure script to convince it to build the shared library
    autoreconf -vi
fi

export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static --with-libiconv-prefix=${prefix}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libgdbm", :libgdbm),
    ExecutableProduct("gdbm_load", :gdbm_load),
    ExecutableProduct("gdbmtool", :gdbmtool),
    ExecutableProduct("gdbm_dump", :gdbm_dump)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531")),
    # We need to use our Readline library to fix this linking issue when
    # building for macOS:
    #    Undefined symbols for architecture x86_64:
    #      "_history_list", referenced from:
    #          _input_history_handler in input-rl.o
    #    ld: symbol(s) not found for architecture x86_64
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module TMMOptics

export thinfilmoptics, PlaneWave, Geometrical, Optical

"""Definition of the beam type."""
abstract type LightSource end
struct PlaneWave{T1} <: LightSource where {T1<:Float64}
    λ; λ0::T1; θ;
end

"""Material type definition with geometrical (physical) and optical thickness."""
abstract type Material end
struct Geometrical{T1, T2} <: Material where {T1<:ComplexF64, T2<:Float64}
    n::Array{T1}; d::T2;
end
struct Optical{T1, T2} <: Material where {T1<:ComplexF64, T2<:Float64}
    n::Array{T1}; f::T2;
end

"""Wrap the output into different types."""
abstract type Output end
struct Spectra{T1, T2} <: Output where {T1<:Float64, T2<:ComplexF64}
    Rp::Array{T1}; Rs::Array{T1}; Tp::Array{T1}; Ts::Array{T1};
    ρp::Array{T2}; ρs::Array{T2}; τp::Array{T2}; τs::Array{T2};
end
struct Field{T1} <: Output where {T1<:Float64}
    emfp::Array{T1}; emfs::Array{T1};
end
struct Bloch{T1} <: Output where {T1<:Number}
    κp::Array{T1}; κs::Array{T1};
end
struct Misc{T1, T2, T3} <: Output where {T1<: Float64, T2<:Number, T3<:Number}
    d::Array{T1}; ℓ::Array{T2}; nλ0::Array{T3};
end
struct AdmPhase{T1} <: Output where {T1<:ComplexF64}
    ηp::Array{T1}; ηs::Array{T1}; δ::Array{T1};
end
struct Results{T1, T2, T3, T4, T5} <:Output where {T1<:Spectra, T2<:Field, T3<:Bloch, T4<:Misc, T5<:AdmPhase}
    Spectra::T1; Field::T2; Bloch::T3; Misc::T4; AdmPhase::T5;
end

"""Data conditioning and main calculations."""
function thinfilmoptics(Beam::T1, Layers::Array{T2,N2}, emfflag::T3=0, h::T3=1, pbgflag::T3=0) where {T1<:PlaneWave, T2<:Material, N2, T3<:Int64}
    λ::Vector{Float64} = vec(Beam.λ) # Work in columns
    λLen::Int64 = length(λ)
    idxλ0::Int64 = findmin(abs.(λ .- Beam.λ0))[2][1] # Find λ closest to λ0
    θ::Vector{Float64} = vec(deg2rad.(Beam.θ)) # Radians
    θLen::Int64 = length(θ)
    nLen::Int64 = size(Layers,2)
    # Build the sequence of index of refractions and the array of thickness depending on the input
    d, nseq, nλ0 = buildArrays(Layers, idxλ0, Beam.λ0, nLen, λLen)
    # Provide the multilayer depth considering the h division
    temp1::Array{Float64,2} = (d[2:end-1] / h) * ones.(1,h) # outer product
    ℓ = cumsum([0; temp1[:]], dims=1)
    ℓ = ℓ[1:end-1]
    # Call transfer matrix method
    tmmout = transferMatrix(nseq, d, λ, θ, emfflag, h, nLen, λLen, θLen)
    # Photonic band gap for crystals without defects
    if (pbgflag == 1) & (nLen > 3)
        κp, κs = pbg(λ, θ, nseq[idxλ0,2], nseq[idxλ0,3], d[2], d[3], idxλ0)
    else
        κp = [0.]; κs = [0.]
    end
    # Return results
    Results(tmmout[1], tmmout[2], Bloch(κp, κs), Misc(d[2:end-1], ℓ, nλ0), tmmout[3])
end # thinfilmoptics(...)

"""Build the sequence of index of refractions and the array of thickness depending on the input."""
function buildArrays(layers::Array{T0,N0}, idxλ0::T1, λ0::T2, nLen::T1, λLen::T1) where {T0<:Material, N0, T1<:Int64, T2<:Float64}
    d = Array{Float64,1}(undef, nLen)
    nλ0 = Array{ComplexF64,1}(undef, nLen)
    nseq = Array{ComplexF64,2}(undef, (λLen, nLen))
    @inbounds for s in eachindex(layers)
        # Refractive index
        nseq[:, s] = layers[s].n
        # Nice to return to plot the profile steps
        nλ0[s] = layers[s].n[idxλ0]
        # convert fraction of central wavelength into physical thickness or pass the input
        if typeof(layers[s]) == Optical{Complex{Float64},Float64}
            d[s] = layers[s].f * λ0 / real(nλ0[s])
        elseif typeof(layers[s]) == Geometrical{Complex{Float64},Float64}
            d[s] = layers[s].d
        end
    end # for s in eachindex(Layers)
    return d, nseq, nλ0
end # buildArrays(..)

"""Computes the reflection and transmission coefficients, and their spectra. The electromagnetic field is calculated if emfflag = 1."""
function transferMatrix(nseq::AbstractArray{T1,N1}, d::AbstractArray{T2,N2}, λ::AbstractArray{T3,N3}, θ::AbstractArray{T4,N4}, emfflag::T5, h::T5, nLen::T5, λLen::T5, θLen::T5) where {T1<:ComplexF64, N1, T2<:Float64, N2, T3<:Number, N3, T4<:Number, N4, T5<:Int64}
    hLen::Int64 = (length(d)-2) * h
    τs = Array{ComplexF64,2}(undef, (λLen, θLen))
    τp = similar(τs); ρs = similar(τs); ρp = similar(τs)
    δ = Array{ComplexF64,3}(undef, (λLen, θLen, nLen))
    ηs = similar(δ); ηp = similar(ηs)
    if emfflag == 1
        emfs = Array{Float64,3}(undef, (λLen, θLen, hLen)); emfp = similar(emfs)
    else
        emfs = [0.]; emfp = [0.]
    end
    # Calculation of complex coefficients of reflection, transmission and emf
    @inbounds for a in eachindex(θ), l in eachindex(λ)
        # Calculation of the optical transfer matrix for all layers
        ηs[l, a, :], ηp[l, a, :], Ψs, Ψp, δ[l, a, :] = totalTransferMatrix(nseq[l, :], d, λ[l], θ[a], nLen)
        # Compute the Fresnell coefficients
        ρs[l, a] = ρ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        ρp[l, a] = ρ(ηp[l, a, 1], ηp[l, a, end], Ψp)
        τs[l, a] = τ(ηs[l, a, 1], ηs[l, a, end], Ψs)
        τp[l, a] = τ(ηp[l, a, 1], ηp[l, a, end], Ψp)
        # EMF
        if emfflag == 1
            emfs[l, a, :] = emfield(nseq[l, :], d, λ[l], θ[a], nLen, h, δ[l, a, :], ηs[l, a, :], Ψs)
            emfp[l, a, :] = emfield(nseq[l, :], d, λ[l], θ[a], nLen, h, δ[l, a, :], ηp[l, a, :], Ψp)
        end
    end # for l in eachindex(λ), a in eachindex(θ)
    return (Spectra(abs2.(ρp), abs2.(ρs), real(ηp[:, :, 1] .* ηp[:, :, end]) .* abs2.(τp), real(ηs[:, :, 1] .* ηs[:, :, end]) .* abs2.(τs), ρp, ρs, τp, τs), Field(emfp, emfs), AdmPhase(ηp, ηs, δ[:, :, 2:end-1]))
end # function transferMatrix(...)

"""Computes the total transfer matrix and admittance for the whole structure at each wavelenth and angle of incidence."""
function totalTransferMatrix(N::AbstractArray{T1,N1}, d::AbstractArray{T2,N2}, λ::T3, θ::T4, NLen::T5) where {T1<:ComplexF64, N1, T2<:Float64, N2, T3<:Number, T4<:Number, T5<:Int64}
    cosϕ = Array{ComplexF64,1}(undef, NLen)
    δ = Array{ComplexF64,1}(undef, NLen)
    # Angle of incidence inside each layer according to the cosine Snell law, to avoid cutoff of total internal reflection with complex angles
    cosϕ[1] = cos(θ)
    @inbounds for c = 2 : NLen
        cosϕ[c] = cosϑ(N[c-1], N[c], cosϕ[c-1])
    end # for c = 2 : NLen
    # Admittance of the first medium for both polarizations
    ηp::Array = ζₚ.(N, cosϕ)
    ηs::Array = ζₛ.(N, cosϕ)
    # Phase shift angles
    @. δ = 2 * π * N * d * cosϕ / λ
    # Total transfer matrix
    tempp::Array = [[complex.([1. 0.; 0. 1.])]; Φ.(δ[2:end-1], ηp[2:end-1])]
    temps::Array = [[complex.([1. 0.; 0. 1.])]; Φ.(δ[2:end-1], ηs[2:end-1])]
    return ηs, ηp, reduce(*, tempp), reduce(*, temps), δ
end # totalTransferMatrix(...)

"""Calculates the optical transfer matrix of a layer. φ: phase shift of the layer y: admittance of the layer, Τ: 2x2 optical tranfer matrix."""
Φ(φ::T1, s::T1) where {T1<:ComplexF64} = [cos(φ) (-im*sin(φ)/s); (-im*sin(φ)*s) cos(φ)]

"""Computes the reflection and transmission coefficients given the admittance and transfer matrix of the whole structure per wavelenth and angle of incidence."""
# Reflection coefficient
ρ(s1::T1, sf::T1, Ψ::AbstractArray{T1,N2}) where {T1<:ComplexF64, N2} = (s1*Ψ[1,1] - Ψ[2,1] + s1*sf*Ψ[1,2] - sf*Ψ[2,2]) / (s1*Ψ[1,1] + Ψ[2,1] + s1*sf*Ψ[1,2] + sf*Ψ[2,2])
# Transmission coefficient
τ(s1::T1, sf::T1, Ψ::AbstractArray{T1,N2}) where {T1<:ComplexF64, N2} = 2 / (s1*Ψ[1,1] + Ψ[2,1] + s1*sf*Ψ[1,2] + sf*Ψ[2,2])

"""Computes the inverse total transfer matrix for the whole structure at each wavelenth and angle of incidence and return the field."""
function emfield(N::AbstractArray{T1,N1}, d::AbstractArray{T2,N2}, λ::T3, θ::T4, NLen::T5, h::T5, δ::AbstractArray{T6,N7}, η::AbstractArray{T6,N8}, Ψ::AbstractArray{T6,N9}) where {T1<:ComplexF64, N1, T2<:Float64, N2, T3<:Number, T4<:Number, T5<:Int64, T6<:ComplexF64, N7, N8, N9}
    m0 = Array{ComplexF64,2}(undef,2,2)
    m1 = complex.([1. 0.; 0. 1.])
    g11 = Array{ComplexF64,1}(undef, (NLen-2)*h)
    g12 = similar(g11)
    # Divide the phase shift by h but keep η as is for each layer
    mδ::Array{ComplexF64} = δ/h
    temp1::Array = Ξ.(mδ, η)
    @inbounds for c = 2 : NLen-1, j = 1 : h
        k::Int64 = h * (c-2) + j
        m1 = temp1[c] * m1 # I could use cumprod outside the loop, temp1 = cumprod([[temp1[2]*m1]; temp1[3:end-1]])
        m0 = m1 * Ψ
        g11[k] = m0[1,1]
        g12[k] = m0[1,2]
    end # for c = 2 : NLen-1, j = 1 : h
    return FI(g11, g12, η[1], η[end], Ψ)
end # emfield(...)

"""Compute the electric field distribution."""
FI(g11::Array{T1,N1}, g12::Array{T1,N1}, s1::T1, sf::T1, Ψ::AbstractArray{T1,N2}) where {T1<:ComplexF64, N1, N2} = abs2.((g11 + sf*g12) / (0.25*(s1*Ψ[1,1] + Ψ[2,1] + s1*sf*Ψ[1,2] + sf*Ψ[2,2])))

"""Calculates the inverse of optical transfer matrix of a layer. φ:  phase shift of the layer, y: admittance of the layer, Τ: 2x2 optical tranfer matrix."""
Ξ(φ::T1, s::T1) where {T1<:ComplexF64} = [cos(φ) (im*sin(φ)/s); (im*sin(φ)*s) cos(φ)]

"""Computes the photonic dispersion of ordered structures (crystals only) alternating two different dielectric layers (pbgflag = 1)."""
function pbg(λ::AbstractArray{T1,N1}, θ::AbstractArray{T2,N2}, n1::T3, n2::T3, d1::T4, d2::T4, idxλ0::T5) where {T1<:Number, N1, T2<:Number, N2, T3<:ComplexF64, T4<:Float64, T5<:Int64}
    ω::Vector{Float64} = 2 * π ./ λ # Angular frequency
    L::Float64 = d1 + d2 # Unit cell
    # Angle of incidence of the second layer with Snell's law of cosine
    cosθ1::Vector{ComplexF64} = cos.(θ)
    cosθ2::Vector{ComplexF64} = cosϑ.(n1, n2, cosθ1)
    # Prefactor for Bloch wavevector
    factor_s = admFactor.(ζₛ.(n1, cosθ1), ζₛ.(n2, cosθ2))
    factor_p = admFactor.(ζₚ.(n1, cosθ1), ζₚ.(n2, cosθ2))
    # Bloch wavevectors
    cosκp = Array{ComplexF64,2}(undef, (length(ω), length(cosθ1)))
    cosκs = similar(cosκp)
    @inbounds for a in eachindex(cosθ1), b in eachindex(ω)
        cosκp[b,a] = cosκ(d1, d2, n1, n2, cosθ1[a], cosθ2[a], ω[b], factor_p[a])
        cosκs[b,a] = cosκ(d1, d2, n1, n2, cosθ1[a], cosθ2[a], ω[b], factor_s[a])
    end # for a in eachindex(cosθ1), b in eachindex(ω)
    return acos.(cosκp)./L, acos.(cosκs)./L
end # function pbg(...)

"""Prefactor for bloch wavevector."""
admFactor(s1::T1, s2::T2) where {T1<:Number, T2<:Number} = 0.5 * (s1^2 + s2^2) / s1 / s2

"""Bloch wavevector."""
cosκ(d1::T1, d2::T1, n1::T2, n2::T2, cosθ1::T3, cosθ2::T4, ω::T1, f::T5) where {T1<:Float64, T2<:ComplexF64, T3<:Number, T4<:Number, T5<:Number} = cos(d1*ω*n1*cosθ1) * cos(ω*d2*n2*cosθ2) - f * sin(ω*d1*n1*cosθ1) * sin(ω*d2*n2*cosθ2)

"""Snell's law in cosine form. Returns de cosine already."""
cosϑ(n1::T1, n2::T1, cosθ::T2) where {T1<:ComplexF64, T2<:Number} = sqrt(1 - (n1/n2)^2 * (1-cosθ^2) )

"""Admittance of the medium for p and s polarizations."""
ζₚ(n::T1, cosθ::T2) where {T1<:ComplexF64,T2<:Number} = n / cosθ
ζₛ(n::T1, cosθ::T2) where {T1<:ComplexF64,T2<:Number} = n * cosθ

end # TMMOptics


text:
using Documenter
using AutomotivePOMDPs

makedocs(
	 modules = [AutomotivePOMDPs],
	 sitename="AutomotivePOMDPs.jl"
	 )

deploydocs(
    repo = "github.com/sisl/AutomotivePOMDPs.jl.git",
)

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs