text:
module Overseer
    using Parameters

    abstract type ComponentData end

    abstract type AbstractComponent{T<:ComponentData} end
    abstract type AbstractGroup end

    abstract type System end

    """
        Abstract type for all entity, component and system ledgers. In order to use the interface,
        make sure that every subtype of AbstractLedger has an overload for the function `ledger` that
        points towards the fields needed for functionality (see ledger.jl for more info).
    """
    abstract type AbstractLedger end

    include("utils.jl")
    include("indices.jl")
    include("entity.jl")
    include("component.jl")
    include("group.jl")
    include("system.jl")
    include("ledger.jl")

    export AbstractLedger, Ledger, System, Stage, Component, SharedComponent, ComponentData, Entity
    export @component, @shared_component, @component_with_kw, @shared_component_with_kw
    export @entities_in

    export update, schedule_delete!, delete_scheduled!, empty_entities!, stage, components, entities, stages
    export prepare, singleton, valid_entities, groups, group, create_group!, regroup!, remove_group!

    # Components
    export swap_order!
end # module


text:
using Weave

docs = [
    "intro.jmd", "intro.jmd",   # weave two times to remove Ipopt's message
]

for d in docs
    weave("$(@__DIR__)/$d")
end

ext: jl


text:
function writevtk(outputname::String,nodes::AbstractArray{<:AbstractVec},cells::Vector{<:TriangleCell})
    info("Writing to $outputname")
    open(outputname,"w") do output
        NN = length(nodes)
        NC = length(cells)
        print(output,"# vtk DataFile Version 3.0\nThis file was generated from a neutral file\nASCII\n")
        print(output,"DATASET UNSTRUCTURED_GRID\n")
        print(output,"POINTS ",NN," double\n")
        for node in nodes
            println(output,xpos(node)," ",ypos(node)," ",zpos(node))
        end
        println(output,"CELLS ",NC," ",4*NC)
        for cell in cells
            #Juila starts indexing at 1, so we subtract 1 from node index.
            println(output,"3 ",cell[1]-1, " ", cell[2]-1, " ",cell[3]-1)
        end
        println(output,"CELL_TYPES ",NC)
        for cell in cells
            println(output,vtkcelltype(cell))
        end
    end
end

@inline writevtk(oname::String,m::AbstractMesh) = writevtk(oname,m.nodes,m.cells)

function append_points_scalar_data(vtkfile::String,dataname::String,data::AbstractArray,header::Bool)
    open(vtkfile,"a") do out
        header && println(out,"POINT_DATA ",length(data))
        println(out,"SCALARS ",dataname," double 1")
        println(out,"LOOKUP_TABLE default")
        for el in data
            println(out,el)
        end
    end
end

function append_points_vector_data(vtkfile::String,dataname::String,data::AbstractArray,header::Bool)
    open(vtkfile,"a") do out
        header && println(out,"POINT_DATA ",length(data))
        println(out,"VECTORS ",dataname," double")
        for el in data
            println(out,join(pos(el)," "))
        end
    end
end

function append_cells_scalar_data(vtkfile::String,dataname::String,data::AbstractArray,header::Bool)
    open(vtkfile,"a") do out
        header && println(out,"CELL_DATA ",length(data))
        println(out,"SCALARS ",dataname," double 1")
        println(out,"LOOKUP_TABLE default")
        for el in data
            println(out,el)
        end
    end
end

function append_cells_vector_data(vtkfile::String,dataname::String,data::AbstractArray,header::Bool)
    open(vtkfile,"a") do out
        header && println(out,"CELL_DATA ",length(data))
        println(out,"VECTORS ",dataname," double")
        for el in data
            println(out,join(pos(el)," "))
        end
    end
end


text:
function createStartingPositionsmeinsneu(x0,x1,y0,y1,N,nx,ny,alpha,verschx,verschy)
X = [];
Y = [];
ALPHA = [];
eps = 0.0000001;
unten = 0;
oben = 0;

xmid = (x0+x1)/2-verschx;
ymid = (y0+y1)/2+verschy;

    if (mod(N,2) == 0)
        unten = N/2;
        oben = (N/2)-1;
    else
        unten = (N-1)/2;
        oben = (N-1)/2;
    end
Deltas = -unten:1:oben;
nRaysPerAngle = length(Deltas);
omega = [cos(alpha),sin(alpha)];
omegaPerp = [-sin(alpha),cos(alpha)];
mycolor = rand(1,3);
x = 0;
y = 0;
r = 0;
s = 0;
t = 0;
xStart = [];
yStart = [];

    for j = 1:nRaysPerAngle

        delta = Deltas[j];

        x = xmid+delta*omegaPerp[1];
        y = ymid+delta*omegaPerp[2];

            r = (x0-x)/cos(alpha);

            if y+r*sin(alpha) < y0-eps || y+r*sin(alpha) > y1+eps
                s = (y0-y)/sin(alpha);
                if x+s*cos(alpha) < x0-eps || x+s*cos(alpha) > x1+eps
                    t = (x1-x)/cos(alpha);
                    if y+t*sin(alpha) < y0-eps || y+t*sin(alpha) > y1+eps
                        xStart = [];
                        yStart = [];
                    else
                        xStart = x1;
                        yStart = y+t*sin(alpha);
                    end
                else
                    xStart = x+s*cos(alpha);
                    yStart = y0;
                end
            else
                xStart = x0;
                yStart = y+r*sin(alpha);
            end

        X = [X;xStart];
        Y = [Y;yStart];
    end

return X,Y,Deltas;

end


text:
using SymPy

x,y = @vars x y
u = x*(1-x)*y*(1-y)
ux = diff(u, x)
uy = diff(u, y)
κ = 1+x^2+y^2
f = diff(ux*κ, x) + diff(uy*κ, y)

s = replace(replace(sympy.julia_code(f), ".*"=>"*"), ".^"=>"^")
print(s)


text:
module CircuitscapeUI

using Blink
using Tachyons
using WebIO
using InteractBase
using InteractBulma
using Circuitscape
using MbedTLS
using Knockout
using Distributed
using CSSUtil
using JSExpr
using AssetRegistry
using MacroTools

include("ui.jl")

export run_ui

end # module


text:
const g = @SMatrix [1  0  0  0;
                    0 -1  0  0;
                    0  0 -1  0;
                    0  0  0 -1] # Minkowski space metric

function dynEq(v,type,τ)
    x = v[1:4]
    u = v[5:8]

    q = type.constants.q
    mₑ = type.constants.mₑ

    dx = u
    du = q/mₑ * Fμν(x, type) * g * u

    return vcat(dx,du)
end

function Mμν(v)
    x⁰, x¹, x², x³, u⁰, u¹, u², u³ = v
    Lx = x²*u³ - x³*u²
    Ly = x³*u¹ - x¹*u³
    Lz = x¹*u² - x²*u¹
    # Nx = x¹*u⁰ - x⁰*u¹
    # Ny = x²*u⁰ - x⁰*u²
    # Nz = x³*u⁰ - x⁰*u³
    L² = Lx^2 + Ly^2 + Lz^2
    # N² = Nx^2 + Ny^2 + Nz^2
    #the physical quantities are m* [1,2,3] , m^2*[4], m/c* [5,6,7] , m^2/c^2*[8],
    return (Lx = Lx, Ly = Ly, Lz = Lz, L² = L², #Nx = Nx, Ny = Ny, Nz = Nz, N² = N²,
            x⁰ = x⁰, x¹ = x¹, x² = x², x³ = x³, u⁰ = u⁰, u¹ = u¹, u² = u², u³ = u³)
end

function σA(A,A²)
    .√abs.((A²-A.^2))
end

function Mμντ(X)
    [Mμν(xτ) for xτ in X.u]
end

function UniformAnnulus(r₁,r₂)
    θ = rand(Uniform(0,2*π))
    u = rand(Uniform(r₁,r₂))
    v = rand(Uniform(0,r₁+r₂))
    ρ = v<u ? u : r₁+r₂ - u
    [SVector{3}(ρ*cos(θ), ρ*sin(θ),0), ρ]
end

R = [UniformAnnulus(0,maxR) for i in 1:N]
RVec = VectorOfArray(R)
idx = sortperm(R, lt=(x,y)->x[2]<y[2])
R₀ = RVec[1,:][idx]
ρ₀ = RVec[2,:][idx]

function xu_cat(X,U,T)
    x0 = c*T
    x = vcat([x0],X)
    u0 = √(c*c + U⋅U)
    u = vcat([u0],U)
    SVector{8}(vcat(x,u))
end

function problem_set(R₀, p, U₀ = [0,0,0])
    V₀ = [xu_cat(r₀,U₀,τi) for r₀ in R₀]
    prob = ODEProblem(dynEq, V₀[1], (τi, τf), p)

    function prob_func(prob,i,repeat)
        ODEProblem(dynEq, V₀[i], (τi, τf), p) #return
    end

    function output_func(sol, i)
        (Mμντ(sol), false) # put sol as argument, if you want the actual trajectories
    end

    EnsembleProblem(prob; prob_func, output_func)
end

# new
function Moment_Avg(M, J; moment_func = identity, f = 1, l = N)
    N_samples = 1 + l - f
    func(i,t) = moment_func(getproperty.((M[i][t], ), J)...)
    [ThreadsX.sum(func(i,t) for i in f:l)/N_samples for t in 1:time_samples+1]
end

# function Radial_Avg(M, args_funcs, bins = 100) # J = (:Lz, (:x¹, :x²))
# data_points = zeros(bins, 1 + sum(length.(args_funcs)) - length(args_funcs))

# particle_idx = 1
# current_M = M[particle_idx][end] # current_particle data at the final time
# current_ρ = ρ₀[particle_idx] # ⟺ hypot(m[begin].x¹, m[begin].x²), but faster since in memory already

#     for i in 1:bins
#         bin_averages = zeros(1 + 2length(args_funcs))
#         tally = 0
#         while current_ρ < i * maxR/bins
#             bin_averages[1]     += current_ρ
#             bin_averages[2:end] += [func(getproperty.((current_M,), args[begin])...)
#                                                         for args in args_funcs
#                                                             for func in args[2:end]]
#             tally += 1
#             try
#                 Midx += 1
#                 current_M = M[particle_idx][end]
#                 current_ρ = ρ₀[particle_idx]
#             catch err
#                 # isa(err, BoundsError)
#                 break
#             end
#         end
#         data_points[i,:] = bin_averages/tally
#     end
#     return data_points
# end

# args_funcs = [((:Lz,), identity, x -> x^2), ((:x¹, :x²), hypot, (x,y) -> x^2 + y^2)]

function Radial_Avg(M, bins = 100)
    LzStats = zeros(bins, 3)
    ρStats  = zeros(bins, 3)
    Midx = 1
    m = M[Midx][end]
    ρ = ρ₀[Midx]
    for i in 1:bins
        ρ_avg, Lz_, Lz²_, ρ_, ρ²_ = 0., 0., 0., 0., 0.
        tally = 0
        while ρ < i * maxR/bins
            ρ_avg += ρ
            Lz_ += m.Lz
            Lz²_ += m.Lz^2
            ρ_ += hypot(m.x¹, m.x²)
            ρ²_ += m.x¹^2 + m.x²^2
            tally += 1
            try Midx < N
                Midx += 1
                m = M[Midx][end]
                ρ = ρ₀[Midx]
            catch err
                break
            end
        end
        rt = ρ_avg/tally
        LzStats[i,:] .= rt, Lz_/tally, Lz²_/tally
        ρStats[i,:] .= rt, ρ_/tally,  ρ²_/tally
    end
    return LzStats, ρStats
end

common_settings = (fontfamily = "computer modern",
                   dpi = 360,
                   label = false,
                   fillalpha=.3,
                   color_palette = :Set1_4,
                   linewidth = 2)

for type in lasers
    pol = polarization(type)

    folder_setup("$(pol)_p$(type.p)_m$(type.m)") # change folder for each simulation
    folder_setup("particles") # separate particles from fields

    file_ext = "$(pol)_p$(type.p)_m$(type.m).$fmt"

    write(io, "$(pol), p = $(type.p), m = $(type.m)\n")

    plot_data = [plot([],[], label = false) for i in 1:8]

    eprob = problem_set(R₀, type)
    M = solve(eprob, Vern9(), EnsembleThreads(), abstol = 1e-15, reltol = 1e-15,
            saveat = (τf-τi)/time_samples, trajectories = N)

    roots = rootsLUT[(p = type.p, m = abs(type.m))]
    idxs = zeros(Int, length(roots))                   # initialize index list

    root_idx = 2                                       # first root is always zero, start from root 2
    for (i,v) in enumerate(ρ₀)                         # ρ₀ is ordered!
        try
            if v > w₀*roots[root_idx]
                idxs[root_idx] = i-1                   # take the last one which respects the condition
                root_idx += 1
            end
        catch err
            # isa(err, BoundsError)
            break
        end
    end

    for j in 2:(length(roots))
        z,  Lz,  L² = Moment_Avg.((M,), (:x³, :Lz, :L²); f = 1+idxs[j-1], l = idxs[j])
        z², Lz², L⁴ = Moment_Avg.((M,), (:x³, :Lz, :L²); moment_func = x -> x^2, f = 1+idxs[j-1], l = idxs[j])
        ρ  = Moment_Avg(M, (:x¹, :x²); moment_func = hypot, f = 1+idxs[j-1], l = idxs[j])
        # ρ² = Moment_Avg(M, (:x¹, :x²); moment_func = (x, y) -> x^2 + y^2, f = 1+idxs[j-1], l = idxs[j])

        σz  = σA(z, z²)/w₀
        # σρ  = σA(ρ, ρ²)/w₀ # not worth computing because it's too wide and breaks plots
        σLz = σA(Lz, Lz²)/w₀/π₀
        σL² = σA(L², L⁴)/w₀/π₀/w₀/π₀

        write(io,
            "from $(roots[j-1]/w₀) to $(roots[j]/w₀)\n"*
            "z ± σz $(z[end]/w₀) ± $(σz[end])\n"*
            "Lz ± σLz $(Lz[end]/w₀/π₀) ± $(σLz[end])\n"*
            "L² ± σL² $(L²[end]/w₀/π₀/w₀/π₀) ± $(σL²[end])\n"
             )

        plot_data[1] = plot!(plot_data[1], Δτ, z/λ;
                                common_settings...,
                                ribbon = σz,
                                seriescolor = j - 1, # these cannot be moved outside due to variable scope
                                xlabel = L"\tau/T_0",
                                ylabel = L"\overline{z}/\lambda")

        # the average radius for particles homogeneously distributed in the annulus between root[i] and root[i+1]
        mean_ρ = 2/3*(roots[j]^3 - roots[j-1]^3)/(roots[j]^2 - roots[j-1]^2)

        plot_data[2] = plot!(plot_data[2], Δτ, ρ/w₀ .- mean_ρ;
                                common_settings...,
                                # ribbon = σρ, too big relative to ρ changes, makes plot unreadable
                                seriescolor = j - 1,
                                xlabel = L"\tau/T_0",
                                ylabel = L"\Delta\overline{\rho}/w_0")

        plot_data[3] = plot!(plot_data[3], Δτ, Lz/w₀/π₀;
                                common_settings...,
                                ribbon = σLz,
                                seriescolor = j - 1,
                                xlabel = L"\tau/T_0",
                                ylabel = L"\overline{L_z}/(w_0 \pi_0)")

        plot_data[4] = plot!(plot_data[4], Δτ, L²/w₀/π₀/w₀/π₀;
                                common_settings...,
                                ribbon = (L²/w₀/π₀/w₀/π₀, σL²),
                                seriescolor = j - 1,
                                xlabel = L"\tau/T_0",
                                ylabel = L"\overline{L^2}/(w_0^2 \pi_0^2)")

    end
    write(io,"\n")

    add_info!(plot_data[1], type)
    savefig(plot_data[1], "z_t_$file_ext")
    add_info!(plot_data[2], type)
    savefig(plot_data[2], "rho_t_$file_ext")
    add_info!(plot_data[3], type)
    savefig(plot_data[3], "Lz_t_$file_ext")
    add_info!(plot_data[4], type)
    savefig(plot_data[4], "L2_t_$file_ext")

    X = [M[i][begin].x¹ for i in 1:N]/w₀ # pulling out the initial x coordinates
    Y = [M[i][begin].x² for i in 1:N]/w₀ # pulling out the initial y coordinates
    Lz = [M[i][end].Lz for i in 1:N] # pulling out the final Lz for colouring
    Lzmax = max(abs.(extrema(Lz))...) # finding maximum
    # colours = Lz/Lzmax # doing a normalization

    write(io, "Maximum Lz transfered = $Lzmax\n\n")

    plot_data[7] = scatter(X,Y;
        marker_z = Lz,
        markershape = circ,
        markerstrokealpha = 0,
        markersize = 1.2,
        seriescolor = cgrad(:PuOr_9, rev = true),
        clims = (-Lzmax, Lzmax),
        fontfamily = "computer modern",
        dpi = 360,
        label = false,
        xlabel = L"x/w_0", ylabel = L"y/w_0",
        aspect_ratio = 1,
        xticks = -4:1:4,
        right_margin = 0.9*Plots.cm,)
    add_info!(plot_data[7], type, true)
    savefig(plot_data[7], "XYLz_$file_ext")

    @time begin
    particle_range = 1:50:N
    Xt  = [[M[i][t].x¹/w₀ for t in 1:time_samples] for i in particle_range]
    Yt  = [[M[i][t].x²/w₀ for t in 1:time_samples] for i in particle_range]
    Zt  = [[M[i][t].x³/λ  for t in 1:time_samples] for i in particle_range]
    Lzt = [[M[i][t].Lz    for t in 1:time_samples] for i in particle_range]
    plot_data[8] = plot(Xt, Yt, Zt;
        line_z = VectorOfArray(Lzt),
        right_margin = 0.9*Plots.cm,
        seriescolor = cgrad(:PuOr_9, rev = true),
        xlabel = L"x/w_0",
        ylabel = L"y/w_0",
        zlabel = L"z/\lambda",
        common_settings...)
    # add_info!(plot_data[8], type) broken for 3D plots
    savefig(plot_data[8], "XYZ_$file_ext")
    end

    LzS, ρS = Radial_Avg(M)

    plot_data[5] = plot(LzS[:,1]/w₀, LzS[:,2]/π₀/w₀,
                        label = false,
                        ribbon = σA(LzS[:,2], LzS[:,3])/π₀/w₀)
    plot_data[5] = add_rect!(plot_data[5], LzS[:,1]/w₀, LzS[:,2]/π₀/w₀, roots;
                        common_settings...,
                        ribbon = σA(LzS[:,2], LzS[:,3])/π₀/w₀,
                        seriescolor = 4,
                        xlabel = L"\rho_0/w_0",
                        ylabel = L"\overline{L_z}(\rho_0)/(w_0 \pi_0)",
                        )
    add_info!(plot_data[5], type)
    savefig(plot_data[5], "Lz_rho_$file_ext")

    plot_data[6] = plot(ρS[:,1]/w₀, ρS[:,2]/w₀ .- ρS[:,1]/w₀;
                        label = false,
                        ribbon = σA(ρS[:,2], ρS[:,2])/w₀)
    plot_data[6] = add_rect!(plot_data[6],
                        ρS[:,1]/w₀, ρS[:,2]/w₀ .- ρS[:,1]/w₀, roots;
                        common_settings...,
                        ribbon = σA(ρS[:,2], ρS[:,3])/w₀,
                        seriescolor = 4,
                        xlabel = L"\rho_0/w_0",
                        ylabel = L"\Delta\overline{\rho}_{\mathrm{final}}(\rho_0)/w_0")
    add_info!(plot_data[6], type)
    savefig(plot_data[6], "rho_rho_$file_ext")

    serialize("plots.jls", plot_data)
    serialize("radial_ρ_data.jls", LzS)
    serialize("radial_Lz_data.jls", LzS)
    # serialize("plot_data.jls", M) save at your own peril, too large to be read efficiently
    cd("../..") # return to output_* base folder
end

max_stars_repo_name: toma-vlad/electrodynamic-simulations

max_issues_repo_name: toma-vlad/electrodynamic-simulations

max_forks_repo_name: toma-vlad/electrodynamic-simulations


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FlowLogFormatParameters <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    version::Any # spec type: Union{ Nothing, Int32 } # spec name: version

    function FlowLogFormatParameters(;type=nothing, version=nothing)
        o = new()
        validate_property(FlowLogFormatParameters, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(FlowLogFormatParameters, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        o
    end
end # type FlowLogFormatParameters

const _property_map_FlowLogFormatParameters = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("version")=>Symbol("version"))
const _property_types_FlowLogFormatParameters = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("version")=>"Int32")
Base.propertynames(::Type{ FlowLogFormatParameters }) = collect(keys(_property_map_FlowLogFormatParameters))
Swagger.property_type(::Type{ FlowLogFormatParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FlowLogFormatParameters[name]))}
Swagger.field_name(::Type{ FlowLogFormatParameters }, property_name::Symbol) =  _property_map_FlowLogFormatParameters[property_name]

const _allowed_FlowLogFormatParameters_type = ["JSON"]

function check_required(o::FlowLogFormatParameters)
    true
end

function validate_property(::Type{ FlowLogFormatParameters }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "FlowLogFormatParameters", :enum, val, _allowed_FlowLogFormatParameters_type)
    end
end


text:
"""
    abstract type LagrangianRefFE{D} <: ReferenceFE{D}

Abstract type representing a Lagrangian reference FE. Lagrangian in the sense that
`get_dof_basis` returns an instance of `LagrangianDofBasis`.
The interface for this type is defined with the methods of `ReferenceFE`
plus the following ones

- [`get_face_own_nodes(reffe::LagrangianRefFE,conf::Conformity)`](@ref)
- [`get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)`](@ref)
- [`get_face_nodes(reffe::LagrangianRefFE)`](@ref)

"""
abstract type LagrangianRefFE{D} <: ReferenceFE{D} end

"""
    get_face_own_nodes(reffe::LagrangianRefFE,conf::Conformity)
"""
function get_face_own_nodes(reffe::LagrangianRefFE,conf::Conformity)
  @abstractmethod
end

"""
    get_face_own_nodes(reffe::LagrangianRefFE)
"""
function get_face_own_nodes(reffe::LagrangianRefFE)
  conf = get_default_conformity(reffe)
  get_face_own_nodes(reffe,conf)
end

function get_face_own_nodes(reffe::LagrangianRefFE,conf::L2Conformity)
  _get_face_own_nodes_l2(reffe)
end

function _get_face_own_nodes_l2(reffe::LagrangianRefFE)
  p = get_polytope(reffe)
  r = [Int[] for i in 1:num_faces(p)]
  r[end] = collect(1:num_nodes(reffe))
  r
end

"""
    get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)
"""
function get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)
  face_own_nodes = get_face_own_nodes(reffe,conf)
  _trivial_face_own_dofs_permutations(face_own_nodes)
end

"""
    get_face_own_nodes_permutations(reffe::LagrangianRefFE)
"""
function get_face_own_nodes_permutations(reffe::LagrangianRefFE)
  conf = get_default_conformity(reffe)
  get_face_own_nodes_permutations(reffe,conf)
end

"""
    get_face_nodes(reffe::LagrangianRefFE)
"""
function get_face_nodes(reffe::LagrangianRefFE)
  @abstractmethod
end

# Tester

"""
    test_lagrangian_reference_fe(reffe::LagrangianRefFE)
"""
function test_lagrangian_reference_fe(reffe::LagrangianRefFE)
  conf = get_default_conformity(reffe)
  @test isa(conf,Conformity)
  test_lagrangian_reference_fe(reffe,conf)
end

function test_lagrangian_reference_fe(reffe::LagrangianRefFE,conf::Conformity)
  test_reference_fe(reffe,conf)
  D = num_dims(reffe)
  node_coordinates = get_node_coordinates(reffe)
  @test isa(node_coordinates,Vector{<:Point{D}})
  @test length(node_coordinates) == num_nodes(reffe)
  node_and_comp_to_dof = get_node_and_comp_to_dof(reffe)
  @test isa(get_dof_to_node(reffe),Vector{Int})
  @test isa(get_dof_to_comp(reffe),Vector{Int})
  @test isa(node_and_comp_to_dof,Vector)
  dof_to_node = get_dof_to_node(reffe)
  @test isa(dof_to_node,Vector{Int})
  @test isa(get_face_own_nodes(reffe,conf),Vector{Vector{Int}})
  @test isa(get_face_own_nodes_permutations(reffe,conf),Vector{Vector{Vector{Int}}})
  @test isa(get_face_nodes(reffe),Vector{Vector{Int}})
end

# Default API

"""
    get_node_coordinates(reffe::LagrangianRefFE)
"""
function get_node_coordinates(reffe::LagrangianRefFE)
  dofs = get_dof_basis(reffe)
  dofs.nodes
end

"""
    num_nodes(reffe::LagrangianRefFE)
"""
num_nodes(reffe::LagrangianRefFE) = length(get_node_coordinates(reffe))

"""
    get_node_and_comp_to_dof(reffe::LagrangianRefFE)
"""
function get_node_and_comp_to_dof(reffe::LagrangianRefFE)
  dofs = get_dof_basis(reffe)
  dofs.node_and_comp_to_dof
end

"""
    get_dof_to_node(reffe::LagrangianRefFE)
"""
function get_dof_to_node(reffe::LagrangianRefFE)
  dofs = get_dof_basis(reffe)
  dofs.dof_to_node
end

"""
    get_dof_to_comp(reffe::LagrangianRefFE)
"""
function get_dof_to_comp(reffe::LagrangianRefFE)
  dofs = get_dof_basis(reffe)
  dofs.dof_to_comp
end

"""
    get_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)
"""
function get_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity)
  n = num_faces(reffe)
  get_face_own_nodes_permutations(reffe,conf)[n]
end

"""
    get_own_nodes_permutations(reffe::LagrangianRefFE)
"""
function get_own_nodes_permutations(reffe::LagrangianRefFE)
  conf = get_default_conformity(reffe)
  get_own_nodes_permutations(reffe,conf)
end

"""
    get_vertex_node(reffe::LagrangianRefFE,conf::Conformity) -> Vector{Int}
"""
function get_vertex_node(reffe::LagrangianRefFE,conf::Conformity)
  d = 0
  p = get_polytope(reffe)
  range = get_dimranges(p)[d+1]
  vertex_to_nodes = get_face_own_nodes(reffe,conf)[range]
  map(first, vertex_to_nodes)
end

"""
    get_vertex_node(reffe::LagrangianRefFE) -> Vector{Int}
"""
function get_vertex_node(reffe::LagrangianRefFE)
  conf = get_default_conformity(reffe)
  get_vertex_node(reffe,conf)
end

"""
    get_face_own_nodes(reffe::LagrangianRefFE,conf::Conformity,d::Integer)
"""
function get_face_own_nodes(reffe::LagrangianRefFE,conf::Conformity,d::Integer)
  p = get_polytope(reffe)
  range = get_dimrange(p,d)
  get_face_own_nodes(reffe,conf)[range]
end

"""
    get_face_own_nodes(reffe::LagrangianRefFE,d::Integer)
"""
function get_face_own_nodes(reffe::LagrangianRefFE,d::Integer)
  conf = get_default_conformity(reffe)
  get_face_own_nodes(reffe,conf,d)
end

"""
    get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity,d::Integer)
"""
function get_face_own_nodes_permutations(reffe::LagrangianRefFE,conf::Conformity,d::Integer)
  p = get_polytope(reffe)
  range = get_dimrange(p,d)
  get_face_own_nodes_permutations(reffe,conf)[range]
end

"""
    get_face_own_nodes_permutations(reffe::LagrangianRefFE,d::Integer)
"""
function get_face_own_nodes_permutations(reffe::LagrangianRefFE,d::Integer)
  conf = get_default_conformity(reffe)
  get_face_own_nodes_permutations(reffe,conf,d)
end

"""
    get_face_nodes(reffe::LagrangianRefFE,d::Integer)
"""
function get_face_nodes(reffe::LagrangianRefFE,d::Integer)
  p = get_polytope(reffe)
  range = get_dimrange(p,d)
  get_face_nodes(reffe)[range]
end

function get_face_own_dofs(reffe::LagrangianRefFE,conf::Conformity)
  dofs = get_dof_basis(reffe)
  face_own_nodes = get_face_own_nodes(reffe,conf)
  face_own_dofs = _generate_face_own_dofs(face_own_nodes, dofs.node_and_comp_to_dof)
  face_own_dofs
end

function _generate_face_own_dofs(face_own_nodes, node_and_comp_to_dof)
  faces = 1:length(face_own_nodes)
  T = eltype(node_and_comp_to_dof)
  comps = 1:num_components(T)
  face_own_dofs = [Int[] for i in faces]
  for face in faces
    nodes = face_own_nodes[face]
    # Node major
    for comp in comps
      for node in nodes
        comp_to_dofs = node_and_comp_to_dof[node]
        dof = comp_to_dofs[comp]
        push!(face_own_dofs[face],dof)
      end
    end
  end

  face_own_dofs

end

function get_face_own_dofs_permutations(reffe::LagrangianRefFE,conf::Conformity)
  dofs = get_dof_basis(reffe)
  face_own_nodes_permutations = get_face_own_nodes_permutations(reffe,conf)
  face_own_nodes = get_face_own_nodes(reffe,conf)
  face_own_dofs = get_face_own_dofs(reffe,conf)
  face_own_dofs_permutations = _generate_face_own_dofs_permutations(
    face_own_nodes_permutations, dofs.node_and_comp_to_dof, face_own_nodes, face_own_dofs)
  face_own_dofs_permutations
end

function  _generate_face_own_dofs_permutations(
  face_own_nodes_permutations, node_and_comp_to_dof, face_own_nodes, face_own_dofs)

  T = eltype(node_and_comp_to_dof)
  ncomps = num_components(T)

  face_own_dofs_permutations = Vector{Vector{Int}}[]
  for  (face, pindex_to_inode_to_pinode) in enumerate(face_own_nodes_permutations)
    idof_to_dof = face_own_dofs[face]
    inode_to_node = face_own_nodes[face]
    pindex_to_idof_to_pidof = Vector{Int}[]
    for inode_to_pinode in pindex_to_inode_to_pinode
      ninodes = length(inode_to_pinode)
      nidofs = ncomps*ninodes
      idof_to_pidof = fill(INVALID_PERM,nidofs)
      for (inode,ipnode) in enumerate(inode_to_pinode)
        if ipnode == INVALID_PERM
          continue
        end
        node = inode_to_node[inode]
        pnode = inode_to_node[ipnode]
        comp_to_pdof = node_and_comp_to_dof[pnode]
        comp_to_dof = node_and_comp_to_dof[node]
        for comp in 1:ncomps
          dof = comp_to_dof[comp]
          pdof = comp_to_pdof[comp]
          idof = findfirst(i->i==dof,idof_to_dof)
          ipdof = findfirst(i->i==pdof,idof_to_dof)
          idof_to_pidof[idof] = ipdof
        end
      end
      push!(pindex_to_idof_to_pidof,idof_to_pidof)
    end
    push!(face_own_dofs_permutations,pindex_to_idof_to_pidof)
  end
  face_own_dofs_permutations
end

function get_face_own_dofs(reffe::LagrangianRefFE,conf::L2Conformity)
  _get_face_own_dofs_l2(reffe)
end

function get_face_own_dofs_permutations(reffe::LagrangianRefFE,conf::L2Conformity)
  face_own_dofs = get_face_own_dofs(reffe,conf)
  _trivial_face_own_dofs_permutations(face_own_dofs)
end

"""
    (==)(a::LagrangianRefFE{D}, b::LagrangianRefFE{D}) where D
"""
function (==)(a::LagrangianRefFE{D}, b::LagrangianRefFE{D}) where D
  t = true
  pa = get_polytope(a)
  pb = get_polytope(b)
  t = t && (pa == pb)
  xa = get_node_coordinates(a)
  xb = get_node_coordinates(b)
  t = t && (xa ≈ xb)
  expsa = get_exponents(get_prebasis(a))
  expsb = get_exponents(get_prebasis(b))
  t = t && (expsa == expsb)
  facedofsa = get_face_dofs(a)
  facedofsb = get_face_dofs(b)
  t = t && (facedofsa == facedofsb)
  ia = get_node_and_comp_to_dof(a)
  ib = get_node_and_comp_to_dof(b)
  t = t && (ia == ib)
  t
end

function (==)(a::LagrangianRefFE, b::LagrangianRefFE)
  false
end

"""
    get_order(reffe::LagrangianRefFE)
"""
function get_order(reffe::LagrangianRefFE)
  get_order(get_prebasis(reffe))
end

"""
    get_orders(reffe::LagrangianRefFE)
"""
function get_orders(reffe::LagrangianRefFE)
  get_orders(get_prebasis(reffe))
end

# Generic implementation
"""
  struct GenericLagrangianRefFE{C,D} <: LagrangianRefFE{D}
    reffe::GenericRefFE{C,D}
    face_nodes::Vector{Vector{Int}}
  end
"""
struct GenericLagrangianRefFE{C,D} <: LagrangianRefFE{D}
  reffe::GenericRefFE{C,D}
  face_nodes::Vector{Vector{Int}}
end

# LagrangianRefFE

get_face_nodes(reffe::GenericLagrangianRefFE) = reffe.face_nodes

# Reffe

num_dofs(reffe::GenericLagrangianRefFE) = num_dofs(reffe.reffe)

get_polytope(reffe::GenericLagrangianRefFE) = get_polytope(reffe.reffe)

get_prebasis(reffe::GenericLagrangianRefFE) = get_prebasis(reffe.reffe)

get_dof_basis(reffe::GenericLagrangianRefFE) = get_dof_basis(reffe.reffe)

get_default_conformity(reffe::GenericLagrangianRefFE) = get_default_conformity(reffe.reffe)

get_face_dofs(reffe::GenericLagrangianRefFE) = get_face_dofs(reffe.reffe)

get_shapefuns(reffe::GenericLagrangianRefFE) = get_shapefuns(reffe.reffe)


text:
# file for homotopy functions
# TODO: add mechanism for selecting which homotopy function to use
# TODO: better idea: move the homotopy function to NLSolver, because it is
#                    physics agnostic

import PDESolver.evalHomotopy

"""
  This function calls the appropriate homotopy function for the Euler module.
"""
function evalHomotopy(mesh::AbstractMesh, sbp::AbstractOperator,
                      eqn::EulerData{Tsol, Tres},
                      opts::Dict, res::Abstract3DArray, t = 0.0) where {Tsol, Tres}

  fill!(res, 0)
  if opts["homotopy_function"] == "ViscousBO"
    sensor = ShockSensorBO{Tsol, Tres}(mesh, sbp, opts)
    capture = eqn.shock_capturing
    res_orig = eqn.res
    eqn.res = res
    applyShockCapturing(mesh, sbp, eqn, opts, sensor, capture)
    eqn.res = res_orig
  elseif opts["homotopy_function"] == "FirstOrderDissipation"

    calcHomotopyDiss(mesh, sbp, eqn, opts, res)
  else
    error("unrecognized homotop_function: $(opts["homotopy_function"])")
  end

  return nothing
end

"""
  Calculate a first order accurate dissipation to use as a homotopy function

  Inputs:
    mesh: a DG mesh
    sbp: an SBP operator
    eqn: an EulerData object
    opts: options dictionary

  Inputs/Outputs:
    res: 3D array to store the homotopy function in

  Note eqn.res is *not* modified by this function.

  Aliasing restrictions: none
"""
function calcHomotopyDiss(mesh::AbstractDGMesh{Tmsh}, sbp,
        eqn::EulerData{Tsol, Tres}, opts,
        res::Abstract3DArray{Tres}) where {Tsol, Tres, Tmsh}

#  println("\nentered calcHomotopyDiss")

  # some checks for when parallelism is enabled
  @assert opts["parallel_data"] == PARALLEL_DATA_ELEMENT
  for i=1:mesh.npeers
    @assert eqn.shared_data[i].recv_waited
  end

  fill!(res, 0.0)

  #----------------------------------------------------------------------------
  # volume dissipation

  t1 = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  t2 = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
  fill!(res, 0.0)
  nrm2 = zeros(Tmsh, mesh.dim)

  # compute the D operator in each direction
  D = zeros(mesh.numNodesPerElement, mesh.numNodesPerElement, mesh.dim)
  for d=1:mesh.dim
    D[:, :, d] = inv(diagm(sbp.w))*sbp.Q[:, :, d]
  end

  for el=1:mesh.numEl
    q_el = sview(eqn.q, :, :, el)
    res_el = sview(res, :, :, el)
    for d1=1:mesh.dim
      fill!(t1, 0.0)
      fill!(t2, 0.0)
      differentiateElement!(sbp, d1, q_el, t1)

      for j=1:mesh.numNodesPerElement
        q_j = sview(eqn.q, :, j, el)
        for k=1:mesh.dim
          nrm2[k] = mesh.dxidx[d1, k, j, el]
        end

        lambda_max = getLambdaMax(eqn.params, q_j, nrm2)
        for k=1:mesh.numDofPerNode
          t2[k, j] = lambda_max*t1[k, j]
        end

      end  # end loop j
     weakDifferentiateElement!(sbp, d1, t2, res_el, SummationByParts.Add(), true)
    end  # end loop dim
  end  # end loop el

  #----------------------------------------------------------------------------
  # interface dissipation

  # interpolate to face rather than using eqn.q_face, in case it hasn't been
  # updated since eqn.q was updated

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)
#  nrm2 = eqn.params.nrm2
  flux = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    qL = sview(eqn.q, :, :, iface_i.elementL)
    qR = sview(eqn.q, :, :, iface_i.elementR)
    resL = sview(res, :, :, iface_i.elementL)
    resR = sview(res, :, :, iface_i.elementR)
#    fill!(q_faceL, 0.0)
#    fill!(q_faceR, 0.0)

    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL, q_faceR)

    # calculate numerical flux at each face node
    for j=1:mesh.numNodesPerFace
      qL_j = sview(q_faceL, :, j)
      qR_j = sview(q_faceR, :, j)

      # get the face normal
      nrm2 = sview(mesh.nrm_face, :, j, i)

      lambda_max = getLambdaMaxSimple(eqn.params, qL_j, qR_j, nrm2)

      for k=1:mesh.numDofPerNode
        flux[k, j] = 0.5*lambda_max*(qL_j[k] - qR_j[k])
      end
    end  # end loop j

    # integrate over the face
    interiorFaceIntegrate!(mesh.sbpface, iface_i, flux, resL, resR)
  end  # end loop i

# the boundary term makes the predictor-corrector algorithm converge slower
  #----------------------------------------------------------------------------
  # boundary dissipation
  # use q_faceL, nrm2, flux  from interface dissipation
  if opts["homotopy_addBoundaryIntegrals"]
    qg = zeros(Tsol, mesh.numDofPerNode)  # boundary state
    for i=1:mesh.numBoundaryFaces
      bndry_i = mesh.bndryfaces[i]
      qL = sview(eqn.q, :, :, bndry_i.element)
      resL = sview(res, :, :, bndry_i.element)
      fill!(q_faceL, 0.0)

      boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, qL, q_faceL)

  #    q_faceL = sview(eqn.q_bndry, :, :, i)
      for j=1:mesh.numNodesPerFace
        q_j = sview(q_faceL, :, j)
  #      dxidx_j = sview(mesh.dxidx_bndry, :, :, j, i)

        # calculate boundary state
        coords = sview(mesh.coords_bndry, :, j, i)
        calcFreeStream(eqn.params, coords, qg)
#        calcInvChannelIC(eqn.params, coords, qg)

        # calculate face normal
        nrm2 = sview(mesh.nrm_bndry, :, j, i)

        # calculate lambda_max
        lambda_max = getLambdaMaxSimple(eqn.params, q_j, qg, nrm2)

        # calculate dissipation
        for k=1:mesh.numDofPerNode
          flux[k, j] = 0.5*lambda_max*(q_j[k] - qg[k])
        end

      end  # end loop j

      # integrate over the face
      boundaryFaceIntegrate!(mesh.sbpface, bndry_i.face, flux, resL)
    end  # end loop i
  end  # end if addBoundaryIntegrals

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

  # shared face integrals
  # use q_faceL, q_faceR, flux from above

  workarr = zeros(q_faceR)
  for peer=1:mesh.npeers
    # get data for this peer
    bndries_local = mesh.bndries_local[peer]
    interfaces_peer = mesh.shared_interfaces[peer]

    qR_peer = eqn.shared_data[peer].q_recv
    nrm_peer = mesh.nrm_sharedface[peer]
    start_elnum = mesh.shared_element_offsets[peer]

    for i=1:length(bndries_local)
      bndry_i = bndries_local[i]
      iface_i = interfaces_peer[i]
      qL_i = sview(eqn.q, :, :, bndry_i.element)
      qR_i = sview(qR_peer, :, :, iface_i.elementR - start_elnum + 1)
      resL = sview(res, :, :, bndry_i.element)

      # interpolate to face
      interiorFaceInterpolate!(mesh.sbpface, iface_i, qL_i, qR_i, q_faceL, q_faceR)
      # compute flux at every face node
      for j=1:mesh.numNodesPerFace
        qL_j = sview(q_faceL, :, j)
        qR_j = sview(q_faceR, :, j)
        nrm2 = sview(nrm_peer, :, j, i)

        # get max wave speed
        lambda_max = getLambdaMaxSimple(eqn.params, qL_j, qR_j, nrm2)

        # calculate flux
        for k=1:mesh.numDofPerNode
          flux[k, j] = 0.5*lambda_max*(qL_j[k] - qR_j[k])
        end
      end  # end loop j

      # integrate over the face
      boundaryFaceIntegrate!(mesh.sbpface, bndry_i.face, flux, resL)
    end  # end loop i
  end  # end loop peer

  # negate for consistency with the physics module
  for i=1:length(res)
    res[i] = -res[i]
  end

#  println("homotopy residual norm = ", norm(vec(res)))

  return nothing
end


text:
using DifferentialEquations, Test

prob = DiscreteProblem(zeros(2),(0.0,1.0))
sol = solve(prob)
@test typeof(sol.alg) <: FunctionMap

max_stars_repo_name: ChrisRackauckas/DiffEq


text:
using MonteCarloX
using StatsBase

function test_rng_mutable(;verbose = false)
    pass = true
    seed = 1000

    if verbose; println("\nstatic mutable random numbers"); end
    rng = MutableRandomNumbers(MersenneTwister(seed),100)
    rng_MT = MersenneTwister(seed)

    if verbose; println("uniform floats:"); end
    for (r1, r2) in zip(rand(rng, 10), rand(rng_MT,10))
        if verbose; println(r1, " vs ", r2); end
        pass &= r1==r2
    end

    if verbose; println("\nexponentially distributed floats:"); end
    for (r1, r2) in zip(randexp(rng, 10), randexp(rng_MT,10))
        if verbose; println(r1, " vs ", r2); end
        pass &= r1==r2
    end

    if verbose; println("\ndynamic mutable random numbers"); end
    rng = MutableRandomNumbers(MersenneTwister(seed),10, mode=:dynamic)
    rng_MT = MersenneTwister(seed)
    for (i, r1, r2) in zip(1:20, rand(rng, 20), rand(rng_MT,20))
        status = "old"
        if i > 10; status = "new"; end
        if verbose; println(status, " ", r1, " vs ", r2); end
        pass &= r1==r2
    end

    if verbose; println("\ntest outcome:"); end
    return pass
end

function test_histogram_set_get(;verbose = false)
    pass = true

    # bin uniform values into histogram (each bin has 10 elements)
    list_vals = [i for i = 1:100]
    hist = fit(Histogram, list_vals, 1:10:101)

    # reset the size of the histogram bins (accessing each bin via random
    # elements that are included in the bin)
    target = [i for i = 1:10]
    hist[1] = target[1]
    hist[12] = target[2]
    hist[23] = target[3]
    hist[34] = target[4]
    hist[45] = target[5]
    hist[56] = target[6]
    hist[67] = target[7]
    hist[78] = target[8]
    hist[89] = target[9]
    hist[100] = target[10]

    # check that the histgram entries correspond to the target values
    for i = 1:100
        if verbose
            println("... $(hist[i]) == $(target[1 + floor(Int, (i - 1) / 10)])")
        end
        pass &= hist[i] == target[1 + floor(Int, (i - 1) / 10)]
    end

    return pass
end


text:
#-------------------------------------------------------------------------------
# core differential equations

function dτdP(P, ::Any, param::Tuple)
    #unpack parameters
    𝒜, idx, g, 𝓂, fT, fμ = param
    #temperature from given profile [K]
    T = fT(P)
    #mean molar mass [mole/kg]
    μ = fμ(T, P)
    #sum of all cross-sections [cm^2/molecule]
    σ = Σ(𝒜, idx, T, P)
    #dτ/dlnP, scaled by the angle m = 1/cos(θ)
    𝓂*dτdP(σ, g, μ) #no Planck emission
end

function dIdP(P, I, param::Tuple)
    #unpack parameters
    𝒜, idx, g, 𝓂, fT, fμ = param
    #temperature from given profile [K]
    T = fT(P)
    #mean molar mass [mole/kg]
    μ = fμ(T, P)
    #sum of all cross-sections [cm^2/molecule]
    σ = Σ(𝒜, idx, T, P)
    #pull out wavenumber
    ν = @inbounds 𝒜.ν[idx]
    #dI/dlnP, scaled by the angle m = 1/cos(θ)
    𝓂*schwarzschild(I, ν, σ, g, μ, T)
end

function dJdP(P, I, param::Tuple)
    #unpack parameters
    𝒜, idx, g, 𝓂, fT, fμ = param
    #temperature from given profile [K]
    T = fT(P)
    #mean molar mass [mole/kg]
    μ = fμ(T, P)
    #sum of all cross-sections [cm^2/molecule]
    σ = Σ(𝒜, idx, T, P)
    #dI/dlnP, scaled by the angle m = 1/cos(θ), without emission
    𝓂*absorption(I, σ, g, μ)
end

#-------------------------------------------------------------------------------
# wrappers for log pressure coordinates

function dτdι(ι, τ, param::Tuple)
    P = ι2P(ι)
    dιfac(P)*dτdP(P, τ, param)
end

function dτdω(ω, τ, param::Tuple)
    P = ω2P(ω)
    dωfac(P)*dτdP(P, τ, param)
end

function dIdω(ω, I, param::Tuple)
    P = ω2P(ω)
    dωfac(P)*dIdP(P, I, param)
end

function dIdι(ι, I, param::Tuple)
    P = ι2P(ι)
    dιfac(P)*dIdP(P, I, param)
end

function dJdι(ι, I, param::Tuple)
    P = ι2P(ι)
    dιfac(P)*dJdP(P, I, param)
end

#-------------------------------------------------------------------------------
# functions for optical depth paths

function 𝓇depth(dτdx::Q,
                x₁::Real,
                x₂::Real,
                𝒜::R,
                idx::Int,
                g::Real,
                𝓂::Real, # 1/cos(θ)
                fT::S,
                fμ::U,
                tol::Float64
                ) where {Q,R<:AbstractAbsorber,S,U}
    #pack parameters, copying the functions to avoid race conditions
    param = (𝒜, idx, g, 𝓂, copyprofile(fT), copyprofile(fμ))
    #integrate with the ODE solver (appears to be faster than quadrature)
    radau(dτdx, zero(x₁), x₁, x₂, param, atol=tol, rtol=tol)
end

function 𝓇depth!(τ,
                 x,
                 dτdx::Q,
                 x₁::Real,
                 x₂::Real,
                 𝒜::R,
                 idx::Int,
                 g::Real,
                 𝓂::Real, # 1/cos(θ)
                 fT::S,
                 fμ::U,
                 tol::Float64
                 ) where {Q,R<:AbstractAbsorber,S,U}
    #pack parameters, copying the functions to avoid race conditions
    param = (𝒜, idx, g, 𝓂, copyprofile(fT), copyprofile(fμ))
    #integrate with the ODE solver (appears to be faster than quadrature)
    radau!(τ, x, dτdx, zero(eltype(τ)), x₁, x₂, param, atol=tol, rtol=tol)
end

#-------------------------------------------------------------------------------
# functions for streams and fluxes up/down the atmosphere, no storage

function 𝓇stream(dIdx::Q, #version of schwarzschild equation
                 I₀::Real, #initial irradiance
                 x₁::Real, #initial pressure coordinate
                 x₂::Real, #final pressure coordinate
                 𝒜::R,
                 idx::Int,
                 g::Real, #gravity [m/s^2]
                 𝓂::Real, #1/cos(θ), where θ is the stream angle
                 fT::S, #temperature profile fT(P)
                 fμ::U, #mean molar mass μ(T,P)
                 tol::Real #integrator error tolerance
                 ) where {Q,R<:AbstractAbsorber,S,U}
    #pack parameters
    param = (𝒜, idx, g, 𝓂, fT, fμ)
    #integrate the Schwarzschild equation in log pressure coords and return
    radau(dIdx, I₀, x₁, x₂, param, atol=tol, rtol=tol)
end

function 𝓇streams(dIdx::Q, #version of schwarzschild equation
                  I₀::Real, #initial irradiance
                  x₁::Real, #initial pressure coordinate
                  x₂::Real, #final pressure coordinate
                  𝒜::R,
                  idx::Int, #index of wavenumber
                  g::Real, #gravity [m/s^2]
                  fT::S, #temperature profile fT(P)
                  fμ::U, #mean molar mass μ(T,P)
                  nstream::Int,
                  tol::Real #integrator error tolerance
                  ) where {Q,R<:AbstractAbsorber,S,U}
    #setup gaussian quadrature nodes
    𝓂, 𝒲 = streamnodes(nstream)
    #copy profiles to avoid interpolator race conditions
    fT, fμ = copyprofiles(fT, fμ)
    #solve schwarzschild w multiple streams, integrating over hemisphere
    M = zero(I₀)
    for i ∈ 1:nstream
        I = 𝓇stream(dIdx, I₀, x₁, x₂, 𝒜, idx, g, 𝓂[i], fT, fμ, tol)
        # integral over hemisphere: ∫∫ I cos(θ) sin(θ) dθ dϕ, where θ∈[0,π/2], ϕ∈[0,2π]
        M += 𝒲[i]*I #𝒲 = 2π*w*cos(θ)*sin(θ), precomputed
    end
    return M
end

#-------------------------------------------------------------------------------
# functions for streams and fluxes up/down the atmosphere with in-place storage

function 𝓇stream!(I, #output/solution vector
                  x, #output/solution coordinates
                  dIdx::Q, #version of schwarzschild equation
                  I₀::Real, #initial irradiance
                  𝒜::R,
                  idx::Int,
                  g::Real, #gravity [m/s^2]
                  𝓂::Real, #1/cos(θ), where θ is the stream angle
                  fT::S, #temperature profile fT(P)
                  fμ::U, #mean molar mass μ(T,P)
                  tol::Real #integrator error tolerance
                  )::Nothing where {Q,R<:AbstractAbsorber,S,U}
    #pack parameters
    param = (𝒜, idx, g, 𝓂, fT, fμ)
    #integrate the Schwarzschild equation in log pressure coords, in-place
    radau!(I, x, dIdx, I₀, x[1], x[end], param, atol=tol, rtol=tol)
    return nothing
end

function 𝓇streams!(M, #output/solution vector
                   x, #output/solution coordinates
                   dIdx::Q, #version of schwarzschild equation
                   I₀::R, #initial irradiance
                   𝒜::S,
                   idx::Int,
                   g::Real, #gravity [m/s^2]
                   fT::U, #temperature profile fT(P)
                   fμ::V, #mean molar mass μ(T,P)
                   nstream::Int,
                   tol::Real #integrator error tolerance
                   )::Nothing where {Q,R<:Real,S<:AbstractAbsorber,U,V}
    @assert length(M) == length(x)
    #setup gaussian quadrature nodes
    𝓂, 𝒲 = streamnodes(nstream)
    #copy profiles to avoid interpolator race conditions
    fT, fμ = copyprofiles(fT, fμ)
    #solve schwarzschild w multiple streams, integrating over hemisphere
    for i ∈ 1:nstream-1
        𝓇stream!(M, x, dIdx, I₀, 𝒜, idx, g, 𝓂[i], fT, fμ, tol)
        M .*= 𝒲[i]/𝒲[i+1] #fancy way to do accumulating dot product
    end
    𝓇stream!(M, x, dIdx, I₀, 𝒜, idx, g, 𝓂[end], fT, fμ, tol)
    M .*= 𝒲[end]
    return nothing
end

#-------------------------------------------------------------------------------
# core function for whole atmosphere upward and downward monochromatic fluxes

function 𝓇monoflux!(M⁺, #downward monochromatic fluxes [W/m^2/cm^-1]
                    M⁻, #upward monochromatic fluxes [W/m^2/cm^-1]
                    P, #pressure coordinates of output
                    ω, #transformed pressure coords
                    ι, #transformed pressure coords
                    𝒜::Q,
                    idx::Int,
                    g::Real, #gravity [m/s^2]
                    fT::R, #temperature profile fT(P)
                    fμ::S, #mean molar mass μ(T,P)
                    fS::U, #incoming stellar radiation fS(ν) [W/m^2]
                    fa::V, #surface albedo fa(ν)
                    θₛ::Real, #stellar radiation angle, corresponds to cos(θ) = 2/3
                    nstream::Int, #number of streams to integrate in both directions
                    tol::Real) where {Q<:AbstractAbsorber,R,S,U,V}
    #setup
    @assert length(M⁻) == length(M⁺) == length(P) == length(ω) == length(ι)
    #surface pressure, assuming ascending pressures
    Pₛ = P[end]
    #surface temperature
    Tₛ = fT(Pₛ)
    #wavenumber
    ν = 𝒜.ν[idx]
    #angle factor for incoming stellar radiation
    𝓂 = 1/cos(θₛ)
    #downward stellar irradiance at ν
    Iₜ⁻ = fS(ν)
    #copy profiles to avoid possible interpolator race conditions
    fT, fμ, fS, fa = copyprofiles(fT, fμ, fS, fa)
    #wipe any previous values
    M⁺[:] .= zero(eltype(M⁺))
    M⁻[:] .= zero(eltype(M⁻))

    #===================================
    downgoing flux throughout atmosphere
    ===================================#
    #cosine of the stellar zenith angle
    c = cos(θₛ)
    #atmospheric contribution to downward flux
    𝓇streams!(M⁻, ι, dIdι, zero(Iₜ⁻), 𝒜, idx, g, fT, fμ, nstream, tol)
    #divide by c before adding the stellar irradiance in-place
    M⁻ ./= c
    #add downward stellar irradiance, **assuming no emission**
    𝓇stream!(M⁻, ι, dJdι, Iₜ⁻, 𝒜, idx, g, 𝓂, fT, fμ, tol)
    #multiply everything by c to get the true flux
    M⁻ .*= c

    #=================================
    upgoing flux throughout atmosphere
    =================================#
    #some of the downward stellar flux is reflected
    Iₛ⁺ = M⁻[end]*fa(ν)/π #Lambertian
    #and the surface emits some radiation
    Iₛ⁺ += planck(ν, Tₛ)
    #upward radiation streams
    𝓇streams!(M⁺, ω, dIdω, Iₛ⁺, 𝒜, idx, g, fT, fμ, nstream, tol)
    #reverse the upward flux to match the coordinate ordering of P and ι
    reverse!(M⁺)

    return nothing
end


text:
abstract type AbstractRewardNormalizer end

"""
    RewardNormalizer()

Use as the `f` reward mapping function in a `AbstractEnv` to create an environment that returns rewards with mean 0 and variance 1.
It can also directly be used in an algorithm to normalize sampled rewards without affecting the stored rewards. To do so, use
`rewardnormalizer(rewards)` to update the estimates and obtain normalized inputs, normalize without update with `rewardnormalizer(rewards, update = false)`.
A running estimate of the mean and standard deviation is kept in memory. `RewardNormalizer` gives the same weight to each sample to compute its estimate.
For a version that uses an exponential estimate, a more suited approach for non-stationary estimates, see `ExpRewardNormalizer`.
#example
```
normalized_env = AbstractEnv(env, RewardNormalizer())
```
"""
mutable struct RewardNormalizer{T} <: AbstractRewardNormalizer
    mean::T
    moment2::T
    std::T
    step_count::Int
end

RewardNormalizer() = RewardNormalizer(0.0f0, 1.0f0, 1.0f0, 0)

function (rn::RewardNormalizer)(rewards; update = true)
    if update
        N = length(rewards)
        rn.step_count += N
        tmp_mean = rn.mean
        rn.mean =
            (rn.step_count - N) / rn.step_count * rn.mean + sum(rewards) / rn.step_count
        rn.moment2 += sum((rewards .- tmp_mean) .* (rewards .- rn.mean))
        rn.std = max(sqrt(rn.moment2 / (max(1, rn.step_count - 1))), eps(rn.std))
    end
    return (rewards .- rn.std) ./ rn.std
end

"""
    ExpRewardNormalizer()

Use as the `f` reward mapping function in a `AbstractEnv` to create an environment that returns rewards with mean 0 and variance 1.
It can also directly be used in an algorithm to normalize sampled rewards without affecting the stored rewards. To do so, use
`exprewardnormalizer(rewards)` to update the estimates and obtain normalized inputs, normalize without update with `rewardnormalizer(rewards, update = false)`.
A running estimate of the mean and standard deviation is kept in memory. `RewardNormalizer` gives the same weight to each sample to compute its estimate.
For a version that does not use an exponential estimate, see `RewardNormalizer`.
#example
```
normalized_env = AbstractEnv(env, RewardNormalizer())
```
"""
mutable struct ExpRewardNormalizer{T} <: AbstractRewardNormalizer
    mean::T
    var::T
    std::T
    factor::T
    first::Bool
end

ExpRewardNormalizer(factor = 0.2f0) = ExpRewardNormalizer(0.0f0, 0.0f0, 0.0f0, factor, true)

function (rn::ExpRewardNormalizer)(rewards; update = true)
    if update
        if rn.first
            rn.first = false
            rn.mean = mean(rewards)
        else
            rn.mean = (1 - rn.factor) * rn.mean + rn.factor * sum(rewards)
        end
        rn.var = (1 - rn.factor) * (rn.var + rn.factor * sum(rewards .^ 2))
        rn.std = sqrt(rn.var)
    end
    return (rewards .- rn.std) ./ rn.std
end


text:
module CustomConfigHexTest

using Hashids
using Test

@testset "Custom Config Hex" begin

config = Hashids.configure(
    salt="this is my salt",
    min_length=30,
    alphabet="xzal86grmb4jhysfoqp3we7291kuct5iv0nd"
)

@testset "encodehex and decodehex back" begin

@test decodehex(config, encodehex(config, "deadbeef")) == "deadbeef"
@test decodehex(config, encodehex(config, "abcdef123456")) == "abcdef123456"
@test decodehex(config, encodehex(config, "ABCDDD6666DDEEEEEEEEE")) == lowercase("ABCDDD6666DDEEEEEEEEE")
@test decodehex(config, encodehex(config, "507f1f77bcf86cd799439011")) == "507f1f77bcf86cd799439011"
@test decodehex(config, encodehex(config, "f00000fddddddeeeee4444444ababab")) == "f00000fddddddeeeee4444444ababab"
@test decodehex(config, encodehex(config, "abcdef123456abcdef123456abcdef123456")) == "abcdef123456abcdef123456abcdef123456"
@test decodehex(config, encodehex(config, "f000000000000000000000000000000000000000000000000000f")) == "f000000000000000000000000000000000000000000000000000f"
@test decodehex(config, encodehex(config, "fffffffffffffffffffffffffffffffffffffffffffffffffffff")) == "fffffffffffffffffffffffffffffffffffffffffffffffffffff"

end  # of "encodehex and decodehex back"

@testset "encodehex" begin

@test encodehex(config, "deadbeef") == "0dbq3jwa8p4b3gk6gb8bv21goerm96"
@test encodehex(config, "abcdef123456") == "190obdnk4j02pajjdande7aqj628mr"
@test encodehex(config, "ABCDDD6666DDEEEEEEEEE") == "a1nvl5d9m3yo8pj1fqag8p9pqw4dyl"
@test encodehex(config, "507f1f77bcf86cd799439011") == "1nvlml93k3066oas3l9lr1wn1k67dy"
@test encodehex(config, "f00000fddddddeeeee4444444ababab") == "mgyband33ye3c6jj16yq1jayh6krqjbo"
@test encodehex(config, "abcdef123456abcdef123456abcdef123456") == "9mnwgllqg1q2tdo63yya35a9ukgl6bbn6qn8"
@test encodehex(config, "f000000000000000000000000000000000000000000000000000f") == "edjrkn9m6o69s0ewnq5lqanqsmk6loayorlohwd963r53e63xmml29"
@test encodehex(config, "fffffffffffffffffffffffffffffffffffffffffffffffffffff") == "grekpy53r2pjxwyjkl9aw0k3t5la1b8d5r1ex9bgeqmy93eata0eq0"

end  # of "encodehex"

end  # of "Custom Config Hex"

end  # module


text:
using AbstractDifferentiation
using Test

@testset "AbstractDifferentiation.jl" begin
    # Write your tests here.
end


text:
using LinearAlgebra

function gauss(a, b, N)
    λ, Q = eigen(SymTridiagonal(zeros(N), [n / sqrt(4n^2 - 1) for n = 1:N-1]))
    @. (λ + 1) * (b - a) / 2 + a, [2Q[1, i]^2 for i = 1:N] * (b - a) / 2
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Day03
export part1, part2

function parse_move(s)
    (s[1], parse(Int, s[2:end]))
end

function parse_wire(s)
    parse_move.(split(s, ","))
end

function follow_wire!(path, f, wire)
    pos = (0,0)
    dist = 0

    for (dir, length) in wire
        (x, y) = (0, 0)
        dir == 'L' && (x=-1)
        dir == 'R' && (x=1)
        dir == 'U' && (y=1)
        dir == 'D' && (y=-1)

        for k = 1:length
            dist += 1
            pos = pos .+ (x, y)
            path[pos] = f(get(path, pos, 0), dist)
        end
    end
end

function closest_intersection(wire1, wire2)
    path = Dict{Tuple{Int64,Int64},Int64}()
    wire = parse_wire(wire1)
    follow_wire!(path, (x,_)->x|1, wire)

    wire = parse_wire(wire2)
    follow_wire!(path, (x,_)->x|2, wire)

    dist = typemax(Int)
    for (pos, cnt) in path
        if cnt == 3
            d = abs(pos[1]) + abs(pos[2])
            dist = min(dist, d)
        end
    end
    dist
end

function part1()
    open("input03.dat") do f
        wire1 = readline(f)
        wire2 = readline(f)
        closest_intersection(wire1, wire2)
    end
end

function fastest_intersection(wire1, wire2)
    update(x, y) = x == 0 ? y : x

    path1 = Dict{Tuple{Int, Int}, Int}()
    wire = parse_wire(wire1)
    follow_wire!(path1, update, wire)

    path2 = Dict{Tuple{Int, Int}, Int}()
    wire = parse_wire(wire2)
    follow_wire!(path2, update, wire)

    dist = typemax(Int)
    for (pos, dist1) in path1
        dist2 = get(path2, pos, 0)
        if dist2 > 0
            dist = min(dist, dist1+dist2)
        end
    end
    dist
end

function part2()
    open("input03.dat") do f
        wire1 = readline(f)
        wire2 = readline(f)
        fastest_intersection(wire1, wire2)
    end
end

end


text:
export connect_signals

"""
```julia
GSignalQuery(
  signal_id::Cuint,
  signal_name::Ptr{Int8},
  itype::Gtk.GType,
  signal_flags::Gtk.GEnum,
  return_type::Gtk.GType,
  n_params::Cuint,
  param_types::Ptr{Gtk.GType})
```
For internal use.

A structure holding information for a specific signal. It is filled in by the
`g_signal_query` C function.

See `struct GSignalQuery` in the GObject reference manual for more information
* `signal_id` - The signal id of the signal being queried, or 0 if the signal
to be queried was unknown.
* `signal_name` - The signal name.
* `itype` - The interface/instance type that this signal can be emitted for.
* `signal_flags` - The signal flags as passed in to the `g_signal_new` C
function
* `return_type` - The return type for user callbacks.
* `n_params` - The number of parameters that user callbacks take.
* `param_types` - The individual parameter types for user callbacks, note that
the effective callback signature is:
```julia
callback(
    data1::Ptr{GObject},
    [param_name_1::ParamType1 ... param_name_n::ParamTypeN],
    data2::Ptr{GObject})::return_type
```
see also: SignalInfo, query_signal
"""
immutable GSignalQuery
  signal_id::Cuint
  signal_name::Ptr{Int8}
  itype::Gtk.GType
  signal_flags::Gtk.GEnum
  return_type::Gtk.GType
  n_params::Cuint
  param_types::Ptr{Gtk.GType}
end

"""
```julia
SignalInfo(
    itype::Type,
    return_type::Type,
    parameter_types::Array{Type, 1})
```
For internal use.

A data type with information about a particular signal.
* `itype` - The interface type that the signal is emitted for.
* `return_type` - The return type for the user callbacks.
* `parameter_types` - The parameter types for the user callbacks.
see also: GSignalQuery, query_signal
"""
immutable SignalInfo
  itype::Type
  return_type::Type
  parameter_types::Array{Type, 1}
end

# Not all types were already covered so add in an auxiliary case
const more_types = Dict{Symbol, Type}(
  :GdkEvent => Ptr{Gtk.GdkEvent},
  :GdkEventButton => Ptr{Gtk.GdkEventButton});

"""
```julia
gtype_to_jtype(t::Gtk.GType)::Type
```
For internal use.

Maps from a `GType` object to the corresponding Julia object type.
* `t` - The type to convert.
"""
function gtype_to_jtype(t::Gtk.GType)
  for (i, id) in enumerate(Gtk.GLib.fundamental_ids)
    if id == t
      return Gtk.GLib.fundamental_types[i][2]
    end
  end
  typename = Gtk.GLib.g_type_name(t)
  if typename in keys(Gtk.GLib.gtype_wrappers)
    return Ptr{Gtk.GLib.gtype_abstracts[typename]}
  end
  if typename in keys(more_types)
    return more_types[typename]
  end
  return Ptr{Void}
end

"""
```julia
query_signal(obj::GObject, signal_name::String)::SignalInfo
```
For internal use.

Looks up details about a particular signal on an object and provides it as
`SignalInfo`.
* `obj` - The GObject to get the details of the signal for.
* `signal_name` - The name of the signal to get details for
"""
function query_signal(obj::GObject, signal_name::String)
  obj_class = Gtk.GLib.G_OBJECT_CLASS_TYPE(obj)
  signal_id = ccall(
    (:g_signal_lookup, Gtk.GLib.libgobject),
    Cuint,
    (Ptr{Int8}, Gtk.GType),
    signal_name,
    obj_class)
  result = Ref{GSignalQuery}()
  ccall(
    (:g_signal_query, Gtk.GLib.libgobject),
    Void,
    (Cuint, Ptr{GSignalQuery}),
    signal_id,
    result)
  SignalInfo(
    gtype_to_jtype(result[].itype),
    gtype_to_jtype(result[].return_type),
    [gtype_to_jtype(gtype) for gtype in unsafe_wrap(Array,
      result[].param_types,
      result[].n_params)])
end

"""
```julia
SignalConnectionData
```
For internal use.

"""
type SignalConnectionData
  handlers::Dict{String, Function}
  data
  warn_pipe::IO
  passthrough::Function
end

"""
```julia
PassthroughData{T, O, P}(
  ret_type::Type{T},
  object_type::Type{O},
  func::Ptr{Void},
  data::P)
```
For internal use.
"""
type PassthroughData{T, O, P}
  ret_type::Type{T}
  object_type::Type{O}
  func::Ptr{Void}
  data::P
end

"""
```
connect_signals_c_function(
    builder,
    object_ptr,
    signal_name_ptr,
    handler_name_ptr,
    connect_object_ptr,
    flags,
    userdata_ptr)
```
For internal use.

A compilable function for configuring signal connections. Passed in its
compiled form to the `gtk_builder_connect_signals_full` C function.
The cfunction version runs but the julia version doesn't
"""
function connect_signals_c_function(
    builder,
    object_ptr,
    signal_name_ptr,
    handler_name_ptr,
    connect_object_ptr,
    flags,
    userdata_ptr)
  userdata = unsafe_pointer_to_objref(userdata_ptr)
  wpipe = userdata.warn_pipe

  handler_name = unsafe_string(handler_name_ptr)
  if !(handler_name in keys(userdata.handlers))
    warn(wpipe, "Signal handler, $handler_name, could not be found")
    return nothing
  end
  handler = userdata.handlers[handler_name]
  passthrough = userdata.passthrough

  object = GObject(object_ptr)
  signal_name = unsafe_string(signal_name_ptr)
  signal_info = query_signal(object, signal_name)

  passed_data = (connect_object_ptr == C_NULL)? userdata.data : GObject(connect_object_ptr)
  ptypes = tuple(Ref{typeof(object)}, signal_info.parameter_types..., Ref{typeof(passed_data)})
  cptr = cfunction(
      handler,
      signal_info.return_type,
      ptypes)
  data = PassthroughData(signal_info.return_type, typeof(object), cptr, passed_data)
  try
    cptr = cfunction(
        handler,
        signal_info.return_type,
        ptypes)
    signal_connect(
        passthrough,
        object,
        signal_name,
        signal_info.return_type,
        (signal_info.parameter_types...),
        false,
        data)
  catch err
    warn(wpipe, "Signal connection failed; signal, $signal_name; handler, $handler_name")
    warn(wpipe, err)
  end

  nothing
end

"""
```julia
connect_signals(
    built::GtkBuilderLeaf,
    handlers::Dict{String, Function},
    userdata,
    passthrough::Function;
    wpipe=Base.STDERR)
```
Connects signals specified within a `GtkBuilder`. Internally calls
C function `gtk_builder_connect_signals_full`. Contains lots of
magic for connecting to functions.
* `built` - A `GtkBuilder` that the signals are connected for.
* `handlers` - A mapping between signal names and callbacks.
* `userdata` - The user provided data to be passed to the callbacks.
* `passthrough` - A generic function that covers requirements for passthrough
functionality.
"""
function connect_signals(
    built::GtkBuilderLeaf,
    handlers::Dict{String, Function},
    userdata,
    passthrough::Function;
    wpipe=Base.STDERR)
  connector = cfunction(
      connect_signals_c_function,
      Void,
      (
          Ptr{GObject},
          Ptr{GObject},
          Ptr{UInt8},
          Ptr{UInt8},
          Ptr{GObject},
          Int,
          Ptr{Void}))
  ccall(
      (:gtk_builder_connect_signals_full, Gtk.libgtk),
      Void,
      (Ptr{GObject}, Ptr{Void}, Ptr{Void}),
      built,
      connector,
      pointer_from_objref(SignalConnectionData(handlers, userdata, wpipe, passthrough)))
end


text:
module AircraftGeometry

## Package imports
#==========================================================================================#

using Base.Math
using Base.Iterators
using DelimitedFiles
using StaticArrays
using CoordinateTransformations
using Rotations
using LinearAlgebra
using SplitApplyCombine
using Interpolations

# Math tools
import ..MathTools: uniform_spacing, linear_spacing, sine_spacing, cosine_spacing, cosine_interp, splitat, adj3, slope, columns, forward_sum, forward_division, forward_difference, weighted_vector, vectarray, extend_yz

# Panel geometry
import ..PanelGeometry: Panel2D, Panel3D, panel_area, normal_vector, transform, make_panels, wetted_area

import ..AeroMDAO: properties

## Types
#==========================================================================================#

abstract type AbstractAircraft end

abstract type AbstractWing <: AbstractAircraft end

## Foil geometry
#==========================================================================================#

include("Foils/foil.jl")
include("Foils/class_shape_transformation.jl")
include("Foils/naca_airfoils.jl")

## Fuselage geometry
#==========================================================================================#

include("fuselage.jl")

## Wing geometry
#==========================================================================================#

aspect_ratio(span, area) = span^2 / area
mean_geometric_chord(span, area) = area / span
taper_ratio(root_chord, tip_chord) = tip_chord / root_chord
area(span, chord) = span * chord
mean_aerodynamic_chord(root_chord, taper_ratio) = (2/3) * root_chord * (1 + taper_ratio + taper_ratio^2)/(1 + taper_ratio)
y_mac(y, b, λ) = y + b / 2 * (1 + 2λ) / 3(1 + λ)
quarter_chord(chord) = 0.25 * chord

include("Wings/halfwing.jl")
include("Wings/wing.jl")
include("Wings/mesh_tools.jl")
include("Wings/mesh_wing.jl")
include("Wings/controls.jl")
include("Wings/sections.jl")

"""
    span(wing :: AbstractWing)

Compute the planform span of an `AbstractWing`.
"""
span(wing :: AbstractWing) = span(wing)

"""
    projected_area(wing :: AbstractWing)

Compute the projected area of an `AbstractWing`` by summing the trapezoidal areas.
"""
projected_area(wing :: AbstractWing) = projected_area(wing)

"""
    mean_aerodynamic_chord(wing :: AbstractWing)

Compute the mean aerodynamic chord of an `AbstractWing`.
"""
mean_aerodynamic_chord(wing :: AbstractWing) = mean_aerodynamic_chord(wing)

"""
    mean_aerodynamic_chord(wing :: AbstractWing)

Compute the coordinates of the mean aerodynamic center of an `AbstractWing`.
"""
mean_aerodynamic_center(wing :: AbstractWing) = mean_aerodynamic_center(wing)

"""
    aspect_ratio(wing :: AbstractWing)

Compute the aspect ratio of an `AbstractWing`.
"""
aspect_ratio(wing) = aspect_ratio(span(wing), projected_area(wing))

"""
    taper_ratio(wing :: AbstractWing)

Compute the taper ratio of an `AbstractWing`, defined as the tip chord length divided by the root chord length.
"""
taper_ratio(wing :: AbstractWing) = taper_ratio(wing)

"""
    properties(wing :: AbstractWing)

Compute the generic properties of interest (span, area, etc.) of an `AbstractWing`.
"""
properties(wing :: AbstractWing) = [ aspect_ratio(wing), span(wing), projected_area(wing), mean_aerodynamic_chord(wing), mean_aerodynamic_center(wing) ]

"""
    mesh_chords(wing :: AbstractWing, n_s :: Vector{Integer}, n_c :: Integer; flip = false)

Mesh the span and chord distributions of an `AbstractWing` with ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions.
"""
mesh_chords(wing :: AbstractWing, span_num, chord_num; spacings = symmetric_spacing(wing)) = mesh_chords(wing, span_num, chord_num; spacings = spacings)

"""
    mesh_wing(wing :: AbstractWing, n_s :: Vector{Integer}, n_c :: Integer; flip = false)

Mesh the span and airfoil coordinate distributions of an `AbstractWing` with ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions.
"""
mesh_wing(wing :: AbstractWing, span_num, chord_num; spacings = symmetric_spacing(wing)) = mesh_wing(wing, span_num, chord_num; spacings = spacings)

"""
    mesh_cambers(wing :: AbstractWing, n_s :: Integer, n_c :: Integer; spacings = symmetric_spacing(wing))

Mesh the camber distribution of a `Wing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions with an `AbstractSpacing`` distribution.
"""
mesh_cambers(wing :: AbstractWing, span_num, chord_num; spacings = symmetric_spacing(wing)) = mesh_cambers(wing, span_num, chord_num; spacings = spacings)

"""
    coordinates(wing :: AbstractWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the planform coordinates of a `HalfWing` given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels..
"""
coordinates(wing :: AbstractWing, span_num, chord_num; span_spacing = Cosine(), chord_spacing = Cosine()) = chop_wing(coordinates(wing), span_num, chord_num; span_spacing = span_spacing, chord_spacing = chord_spacing)

"""
    chord_coordinates(wing :: AbstractWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the chord coordinates of an `AbstractWing` given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels, with an option to flip the signs of the ``y``-coordinates.
"""
chord_coordinates(wing :: AbstractWing, span_num, chord_num; span_spacing = Cosine(), chord_spacing = Cosine()) = chord_coordinates(wing, span_num, chord_num; span_spacing = span_spacing, chord_spacing = chord_spacing)

"""
    camber_coordinates(wing :: AbstractWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the camber coordinates of an `AbstractWing` given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels, with an option to flip the signs of the ``y``-coordinates.
"""
camber_coordinates(wing :: AbstractWing, span_num, chord_num; span_spacing = Cosine(), chord_spacing = Cosine()) = camber_coordinates(wing, span_num, chord_num; span_spacing = span_spacing, chord_spacing = chord_spacing)

"""
    surface_coordinates(wing :: AbstractWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the surface coordinates of an `AbstractWing` given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels, with an option to flip the signs of the ``y``-coordinates.
"""
surface_coordinates(wing :: AbstractWing, span_num, chord_num; span_spacing = Cosine(), chord_spacing = Cosine()) = surface_coordinates(wing, span_num, chord_num; span_spacing = span_spacing, chord_spacing = chord_spacing)

function Base.show(io :: IO, wing :: AbstractWing)
    println(io, supertype(typeof(wing)),  " with ", length(spans(wing)), " spanwise section(s).")
    println(io, "Aspect Ratio: ", aspect_ratio(wing))
    println(io, "Span (m): ", span(wing))
    println(io, "Projected Area (m): ", projected_area(wing))
    println(io, "Mean Aerodynamic Chord (m): ", mean_aerodynamic_chord(wing))
    println(io, "Mean Aerodynamic Center (m): ", mean_aerodynamic_center(wing))

    nothing
end

end


text:
const USE_GPU = true
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3);
else
    @init_parallel_stencil(Threads, Float64, 3);
end

@parallel function copy3D!(T2, T, Ci)
    @all(T2) = @all(T) + @all(Ci);
    return
end

function memcopy3D()
# Numerics
nx, ny, nz = 512, 512, 512;                              # Number of gridpoints in dimensions x, y and z
nt  = 100;                                               # Number of time steps

# Array initializations
T   = @zeros(nx, ny, nz);
T2  = @zeros(nx, ny, nz);
Ci  = @zeros(nx, ny, nz);

# Initial conditions
Ci .= 1/2.0;
T  .= 1.7;
T2 .= T;

# Time loop
for it = 1:nt
    if (it == 11) global t0=time(); end  # Start measuring time.
    @parallel copy3D!(T2, T, Ci);
    T, T2 = T2, T;
end
time_s=time()-t0

# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number);      # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it  = time_s/(nt-10);                                  # Execution time per iteration [s]
T_eff = A_eff/t_it;                                      # Effective memory throughput [GB/s]
println("time_s=$time_s T_eff=$T_eff");
end

memcopy3D()


text:
include("clamped.jl")

# ====== Converting to MAT files =======

using MAT

function write_mat(; N, damped=true)

    if damped
        sys = clamped(N=N, a=1e-6, b=1e-6, damped=true, homogeneize=false)
        name = "CB21d_$N.mat"

    else
        sys = clamped(N=N, damped=false, homogeneize=false)
        name = "CB21_$N.mat"
    end

    dat = Dict("A" => sys.A, "b" => sys.c)
    matwrite(name, dat; compress = true)
end

#=
for N in [100, 500, 1000]
    for damped in [true, false]
        write_mat(N=N, damped=damped)
    end
end
=#

# ====== Converting to SX files =======

using SpaceExParser

# WARNING requires manual modification of the xml
function write_sxmodel_clamped_constant(; N, damped=true)
    if damped
        sys = clamped(N=N, a=1e-6, b=1e-6, damped=true, homogeneize=true)
        name = "CB21Cd_$N.xml"
    else
        sys = clamped(N=N, damped=false, homogeneize=true)
        name = "CB21C_$N.xml"
    end
    writesxmodel(name, sys)
end

function write_sxmodel_clamped_timevarying(; N, U = Interval(9900, 10100), damped=true)

    X = Universe(N)

    if damped
        sys = clamped(N=N, a=1e-6, b=1e-6, damped=true, homogeneize=false)
        name = "CB21Fd_$N.xml"

    else
        sys = clamped(N=N, damped=false, homogeneize=false)
        name = "CB21F_$N.xml"
    end

    A = sys.A
    B = hcat(sys.c)
    S = @system(x' = A*x + B*u, x ∈ X, u ∈ U)
    return writesxmodel(name, S)
end

#=
for func in [write_sxmodel_clamped_constant, write_sxmodel_clamped_timevarying]
for N in [100, 500, 1000]
    for damped in [true, false]
        func(N=N, damped=damped)
    end
end
end
=#

max_stars_repo_name: JuliaReach/ClampedFreeBar

max_issues_repo_name: JuliaReach/ClampedFreeBar

max_forks_repo_name: JuliaReach/ClampedFreeBar


text:
# ---
# title: 1044. Longest Duplicate Substring
# id: problem1044
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Binary Search, Suffix Array
# link: <https://leetcode.com/problems/longest-duplicate-substring/description/>
# hidden: true
# ---
#
# Given a string `s`, consider all _duplicated substrings_ : (contiguous)
# substrings of s that occur 2 or more times. The occurrences may overlap.
#
# Return **any** duplicated substring that has the longest possible length. If
# `s` does not have a duplicated substring, the answer is `""`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "banana"
#     Output: "ana"
#
#
# **Example 2:**
#
#
#
#     Input: s = "abcd"
#     Output: ""
#
#
#
#
# **Constraints:**
#
#   * `2 <= s.length <= 3 * 104`
#   * `s` consists of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
(G::NGrid)(A::Vector{Float64}, x::Array{Float64, 2}, y = zeros(size(x, 1))) = jl_ainterp(G, A, x, y)

(G::NGrid)(A::Array{Float64, 2}, x::Array{Float64, 2}, y = zeros(size(x, 1), size(A, 2))) = jl_ainterp(G, A, x, y)


text:
include("day4_solution.jl")
using Test

# byr valid:   2002
# byr invalid: 2003
@test validate_birth_year("2002")
@test !validate_birth_year("2003")

# hgt valid:   60in
# hgt valid:   190cm
# hgt invalid: 190in
# hgt invalid: 190
@test validate_height("60in")
@test validate_height("190cm")
@test !validate_height("190in")
@test !validate_height("190")

# hcl valid:   #123abc
# hcl invalid: #123abz
# hcl invalid: 123abc
@test validate_hair_color("#123abc")
@test !validate_hair_color("#123abz")
@test !validate_hair_color("123abc")

# ecl valid:   brn
# ecl invalid: wat
@test validate_eye_color("brn")
@test !validate_eye_color("wat")

# pid valid:   000000001
# pid invalid: 0123456789
@test validate_pid("000000001")
@test !validate_pid("0123456789")

bad_pwd_1 = "eyr:1972 cid:100
hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926"

bad_pwd_2 = "iyr:2019
hcl:#602927 eyr:1967 hgt:170cm
ecl:grn pid:012533040 byr:1946"

bad_pwd_3 = "hcl:dab227 iyr:2012
ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277"

bad_pwd_4 = "hgt:59cm ecl:zzz
eyr:2038 hcl:74454a iyr:2023
pid:3556412378 byr:2007"

@test !is_valid_passport(bad_pwd_1)
@test !is_valid_passport(bad_pwd_2)
@test !is_valid_passport(bad_pwd_3)
@test !is_valid_passport(bad_pwd_4)

good_pwd_1 = "pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
hcl:#623a2f"

good_pwd_2 = "eyr:2029 ecl:blu cid:129 byr:1989
iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm"

good_pwd_3 = "hcl:#888785
hgt:164cm byr:2001 iyr:2015 cid:88
pid:545766238 ecl:hzl
eyr:2022"

good_pwd_4 = "iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719"

@test is_valid_passport(good_pwd_1)
@test is_valid_passport(good_pwd_2)
@test is_valid_passport(good_pwd_3)
@test is_valid_passport(good_pwd_4)


text:
# init.jl

println("Adding packages\n")
@time Pkg.add("Makie")
@time Pkg.add("Cairo")
@time Pkg.add("Documenter")
@time Pkg.add("ExcelReaders")

println("Done!\n")

max_stars_repo_name: anthony-wang/julia-scripts

max_issues_repo_name: anthony-wang/julia-scripts

max_forks_repo_name: anthony-wang/julia-scripts


text:
export predict
using TupleVectors
using SampleChains

predict(args...; kwargs...) = predict(Random.GLOBAL_RNG, args...; kwargs...)

# TODO: Fix this hack
predict(d::AbstractMeasure, x) = x
predict(d::Dists.Distribution, x) = x
predict(d::AbstractModel, args...; kwargs...) = predict(Random.GLOBAL_RNG, d, args...; kwargs...)

@inline function predict(rng::AbstractRNG, m::AbstractModel, nt::NamedTuple{N}) where {N}
    pred = predictive(Model(m), N...)
    rand(rng, pred(merge(argvals(m), nt)))
end

predict(rng::AbstractRNG, m::AbstractModel; kwargs...) = predict(rng, m, (;kwargs...))

@inline function predict(rng::AbstractRNG, d::AbstractModel, nt::LazyMerge)
    predict(rng, d, convert(NamedTuple, nt))
end

function predict(rng::AbstractRNG, d::AbstractModel, post::AbstractVector{<:NamedTuple{N}}) where {N}
    m = Model(d)
    pred = predictive(m, N...)
    args = argvals(d)
    y1 = rand(rng, pred(merge(args,post[1])))
    n = length(post)
    v = TupleVectors.chainvec(y1, n)
    @inbounds for j in 2:n
        newargs = merge(args,post[j])
        v[j] = rand(rng, pred(newargs))
    end

    v
end

function predict(rng::AbstractRNG, d::ConditionalModel, post::MultiChain)
    [predict(rng, d, c) for c in getchains(post)]
end


text:
function Base.show(io::IO, ist::T) where T <: Union{CleanEnergyReserve, InertiaReserve, ThermalCleanEnergy}
    print(io, string(nameof(typeof(ist))), "(")
    is_first = true
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        getter_name = Symbol("get_$name")
        if field_type <: IS.TimeSeriesContainer ||
           field_type <: InfrastructureSystemsInternal
            continue
        elseif hasproperty(PowerSystems, getter_name)
            getter_func = getproperty(PowerSystems, getter_name)
            val = getter_func(ist)
        else
            val = getproperty(ist, name)
        end
        if is_first
            is_first = false
        else
            print(io, ", ")
        end
        print(io, val)
    end
    print(io, ")")
    return
end

function Base.show(io::IO, ::MIME"text/plain", ist::T) where T <: Union{CleanEnergyReserve, InertiaReserve, ThermalCleanEnergy}
    default_units = false
    if !PSY.has_units_setting(ist)
        print(io, "\n")
        @warn(
            "SystemUnitSetting not defined, using NATURAL_UNITS for displaying device specification."
        )
        PSY.set_units_setting!(
            ist,
            PSY.SystemUnitsSettings(100.0, PSY.UNIT_SYSTEM_MAPPING["NATURAL_UNITS"]),
        )
        default_units = true
    end
    try
        print(io, summary(ist), ":")
        for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
            obj = getproperty(ist, name)
            getter_name = Symbol("get_$name")
            if (obj isa InfrastructureSystemsInternal) && !default_units
                print(io, "\n   ")
                show(io, MIME"text/plain"(), obj.units_info)
                continue
            elseif obj isa IS.TimeSeriesContainer ||
                   obj isa InfrastructureSystemsType ||
                   obj isa Vector{<:InfrastructureSystemsComponent}
                val = summary(getproperty(ist, name))
            elseif hasproperty(PowerSystems, getter_name)
                getter_func = getproperty(PowerSystems, getter_name)
                val = getter_func(ist)
            else
                val = getproperty(ist, name)
            end
            # Not allowed to print `nothing`
            if isnothing(val)
                val = "nothing"
            end
            print(io, "\n   ", name, ": ", val)
        end
    finally
        if default_units
            PSY.set_units_setting!(ist, nothing)
        end
    end
    return
end

max_stars_repo_name: sourabhdalvi/PowerSystemExtensions

max_issues_repo_name: sourabhdalvi/PowerSystemExtensions

max_forks_repo_name: sourabhdalvi/PowerSystemExtensions


text:
using Test

println("Testing...")
@test 1 == 1


text:
"A basis of Chebyshev polynomials of the second kind on the interval `[-1,1]`."
struct ChebyshevU{T} <: OPS{T}
    n			::	Int
end

ChebyshevU(n::Int) = ChebyshevU{Float64}(n)

similar(b::ChebyshevU, ::Type{T}, n::Int) where {T} = ChebyshevU{T}(n)

show(io::IO, d::ChebyshevU{Float64}) = print(io, "ChebyshevU($(length(d)))")
show(io::IO, d::ChebyshevU{T}) where T = print(io, "ChebyshevU{$(T)}($(length(d)))")

function unsafe_eval_element(b::ChebyshevU, idx::PolynomialDegree, x::Real)
    # Don't use the formula when |x|=1, because it will generate NaN's
    d = degree(idx)
    abs(x) < 1 ? sin((d+1)*acos(x))/sqrt(1-x^2) : recurrence_eval(b, idx, x)
end

first_moment(b::ChebyshevU{T}) where {T} = convert(T, pi)/2

interpolation_grid(b::ChebyshevU{T}) where {T} = ChebyshevUNodes{T}(b.n)

iscompatible(dict::ChebyshevU, grid::ChebyshevUNodes) = length(dict) == length(grid)
issymmetric(::ChebyshevU) = true
measure(dict::ChebyshevU{T}) where {T} = ChebyshevUWeight{T}()
hasmeasure(::ChebyshevU) = true

function innerproduct_native(b1::ChebyshevU, i::PolynomialDegree, b2::ChebyshevU, j::PolynomialDegree, m::ChebyshevUWeight;
			T = coefficienttype(b1), options...)
	if i == j
		convert(T, pi)/2
	else
		zero(T)
	end
end

# Parameters alpha and beta of the corresponding Jacobi polynomial
jacobi_α(b::ChebyshevU{T}) where {T} = one(T)/2
jacobi_β(b::ChebyshevU{T}) where {T} = one(T)/2

# See DLMF, Table 18.9.1
# http://dlmf.nist.gov/18.9#i
rec_An(b::ChebyshevU{T}, n::Int) where {T} = convert(T, 2)

rec_Bn(b::ChebyshevU{T}, n::Int) where {T} = zero(T)

rec_Cn(b::ChebyshevU{T}, n::Int) where {T} = one(T)

support(b::ChebyshevU{T}) where {T} = ChebyshevInterval{T}()

"A Chebyshev polynomial of the second kind"
struct ChebyshevUPolynomial{T} <: OrthogonalPolynomial{T}
    degree  ::  Int
end

ChebyshevUPolynomial{T}(p::ChebyshevUPolynomial) where {T} = ChebyshevUPolynomial{T}(p.degree)

show(io::IO, p::ChebyshevUPolynomial) = print(io, "U_$(degree(p))(x) (Chebyshev polynomial of the second kind)")

convert(::Type{TypedFunction{T,T}}, p::ChebyshevUPolynomial) where {T} = ChebyshevUPolynomial{T}(p.degree)

support(::ChebyshevUPolynomial{T}) where {T} = ChebyshevInterval{T}()

(p::ChebyshevUPolynomial{T})(x) where {T} = eval_element(ChebyshevU{T}(degree(p)+1), degree(p)+1, x)

basisfunction(dict::ChebyshevU, idx) = basisfunction(dict, native_index(dict, idx))
basisfunction(dict::ChebyshevU{T}, idx::PolynomialDegree) where {T} = ChebyshevUPolynomial{T}(degree(idx))

dictionary(p::ChebyshevUPolynomial{T}) where {T} = ChebyshevU{T}(degree(p)+1)
index(p::ChebyshevUPolynomial) = degree(p)+1