text:
function LaguerreSolver()
    return LaguerreSolver(())
end

function LaguerreSolver(arg0::jdouble)
    return LaguerreSolver((jdouble,), arg0)
end

function LaguerreSolver(arg0::jdouble, arg1::jdouble)
    return LaguerreSolver((jdouble, jdouble), arg0, arg1)
end

function LaguerreSolver(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return LaguerreSolver((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function do_solve(obj::LaguerreSolver)
    return jcall(obj, "doSolve", jdouble, ())
end

function solve_all_complex(obj::LaguerreSolver, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "solveAllComplex", Vector{JComplex}, (Vector{jdouble}, jdouble), arg0, arg1)
end

function solve_complex(obj::LaguerreSolver, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "solveComplex", JComplex, (Vector{jdouble}, jdouble), arg0, arg1)
end


text:
__precompile__()

module AutomotiveDrivingModels

# using Compat
using Reexport
# using Discretizers
using Parameters

@reexport using Records
@reexport using Distributions
@reexport using Vec
# @reexport using DataFrames

export nframes

include("interface/main.jl") # core interface
include("1d/main.jl") # 1D interface
include("2d/main.jl") # 2D interface

end # module

ext: jl


text:
module OptimMPB

using Reexport

using ForwardDiff
using Calculus
@reexport using Optim
@reexport using MathProgBase
import MathProgBase: NonlinearModel, getsense, numvar, numconstr, SolverInterface.optimize!
                     MathProgBase.SolverInterface.loadproblem!
import Optim: Optimizer, OptimizationResults, optimize, MultivariateOptimizationResults,
              initial_state, minimizer, minimum, iterations, converged, x_converged,
              f_converged, f_tol, g_tol, g_converged, iteration_limit_reached, f_calls, method

include("interface.jl")

export OptimSolver

end # module


text:
dense2arg = Dict{Symbol,Any}(
#:^ => :todo, # Matrix^Number defined in dense
)

# TODO:

# eval
# scale!
# isposdef!
# isposdef
# stride1: Not exported
# mapreduce_seq_impl: Not exported
# norm
# vecnorm1: Not exported
# vecnorm2: Not exported
# triu!
# triu
# tril!
# tril
# gradient
# diagind
# diag
# diagm
# trace
# kron
# ^
# expm
# expm!: Not exported
# rcswap!: Not exported
# logm
# sqrtm
# inv
# factorize
# \
# pinv
# nullspace
# cond
# sylvester
# lyap


text:
function Δ(n::Integer, ℓbar::Integer, τ::T, σ::T,
           t::OffsetArray{T}) where T <: AbstractFloat
    ℓ = n - ℓbar
    kn = t[n] - t[n-1]
    kℓ = t[ℓ] - t[ℓ-1]
    two_D = (t[n]+t[n-1]) - (t[ℓ]+t[ℓ-1])
    return ( τ * kn - σ * kℓ ) / two_D
end


text:
# Jonas August 1/2022

#=

This problem was asked by Google.

The power set of a set is the set of all its subsets. Write a function that, given a set, generates its power set.

For example, given the set {1, 2, 3}, it should return {{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}.

You may also use a list or array to represent a set.

=#

#=

The possible subsets are equivalent to all possible binary indicator vectors, where each bit represents presence (1) or absence (0) of a given element.

Thus 2^n possibilities entries in the power set.

What about repeated elements, and thus indistinguishable subsets?

I will assume that we are to compute the power set _of a set_ (which has distinct elements, by definition), that the array representing the set has unique elements (no repeated elements).  I can test this.

=#

#=

Another approach (more complex, not implemented):

- Let n be length of v.
- Let p be the power set of v, initially empty.
- Each entry in p will be a subset of the elements of v.
- Generate all n-bit patterns, and for each bit pattern select the components of v based where there are 1's in the pattern, and add it to p.

How to generate bit pattern?
Count from 0 through 2^n-1 in an (unsigned?) int of sufficient bits, or create a length n bit array, and "count".  If we use int representation, then counting is built-in.  If we make binary array, then we have to implement increment (with carry).

=#

"Construct array of all subsets of v"
function powerset(v)
   # TODO: assert that there are no repeated elements in v
   p = [[]]
   for item in v
       n = []   # new sets with item appended
       for set in p
           push!(n, push!(copy(set), item))
       end
       push!(p, n...)   # combine p and n
   end
   p
end

using Test
@test powerset([]) == [[]]
@test powerset([1]) == [[], [1]]
@test powerset([1,2]) == [[],[1],[2],[1,2]]
# Set needed because ordering in computed powerset is different than test case.
@test Set(powerset([1,2,3])) == Set([[],[1],[2],[3],[1,2],[1,3],[2,3],[1,2,3]])

max_stars_repo_path: Power Set - Problem 37.jl

max_issues_repo_path: Power Set - Problem 37.jl

max_forks_repo_path: Power Set - Problem 37.jl


text:
function build_incidence(nodes::Array{Junction}, links::Array{T}) where T<:Link

    linkcount = length(links)
    nodecount = length(nodes)
    num_junc = Dict{String,Int}() #maps the junction number to row in incidence matrix
    link_axis = Array{String,1}()
    for (ix, junction) in enumerate(nodes)
        num_junc[junction.name] = ix
    end

    A = zeros(Int64,nodecount,linkcount)

   #build incidence matrix
   #incidence_matrix = A
    for (ix,b) in enumerate(links)
        if typeof(b) <: ControlPipe
            A[num_junc[b.pipe.connectionpoints.from.name], ix] =  1;

            A[num_junc[b.pipe.connectionpoints.to.name], ix] = -1;
            push!(link_axis, b.pipe.name)
        else
            A[num_junc[b.connectionpoints.from.name], ix] =  1;

            A[num_junc[b.connectionpoints.to.name], ix] = -1;
            push!(link_axis, b.name)
        end

    end
    return  A, link_axis
end

function build_incidence_null(A::AbstractArray{Int64})
    null_A = nullspace(A)
    return null_A
end
# function WaterSystem(links::Array{<:Link})
# function WaterSystem(nodes, junctions, tanks, reservoirs, links, pipes, valves, pumps, demands, network, simulation)
#
#         new(nodes,
#             junctions,
#             tanks,
#             reservoirs,
#             links,
#             pipes,
#             valves,
#             pumps,
#             demands,
#             network,
#             simulation)
# end

# maps name to position in matrix
function _make_ax_ref(ax::Vector)
    ref = Dict{eltype(ax), Int}()
    for (ix, el) in enumerate(ax)
        if haskey(ref, el)
            @error("Repeated index element $el. Index sets must have unique elements.")
        end
        ref[el] = ix
    end
    return ref
end

struct Incidence
    data:: Array{Int64,2}
    axes::Tuple{Array{String,1},Array{String,1}}
    lookup::Tuple{Dict{String,Int64},Dict{String,Int64}}
end

function Incidence(nodes::Array{Junction},links::Array{<:Link})
    node_axis = [junc.name for junc in nodes]
    A, link_axis = build_incidence(nodes, links)
    axes = (link_axis, node_axis)
    look_up = (_make_ax_ref(link_axis), _make_ax_ref(node_axis))
    return Incidence(A, axes, look_up)
end

return Incidence, build_incidence, build_incidence_null


text:
const valid_union = Union{
    Bool,
    Int8,
    UInt8,
    Int16,
    UInt16,
    Int32,
    UInt32,
    Int64,
    UInt64,
    Float32,
    Float64,
    ComplexF32,
    ComplexF64,
    }
const valid_vec = [
    Bool,
    Int8,
    UInt8,
    Int16,
    UInt16,
    Int32,
    UInt32,
    Int64,
    UInt64,
    Float32,
    Float64,
    ComplexF32,
    ComplexF64,
]

const gxb_union = Union{ComplexF32, ComplexF64}
const gxb_vec = [ComplexF32, ComplexF64]

struct GBType{T} <: AbstractGBType
    p::libgb.GrB_Type
end

Base.unsafe_convert(::Type{libgb.GrB_Type}, s::AbstractGBType) = s.p

GBType{T}(name::AbstractString) where {T<:valid_union} = GBType{T}(load_global(name))

function Base.show(io::IO, ::MIME"text/plain", t::GBType{T}) where T
    print(io, "GBType(" * string(T) * "): ")
    gxbprint(io, t)
end

struct GBAllType <: AbstractGBType
    p::Ptr{libgb.GrB_Index}
end
Base.show(io::IO, ::MIME"text/plain", t::GBAllType) = print(io, "GraphBLAS type: GrB_ALL")
Base.unsafe_convert(::Type{Ptr{libgb.GrB_Index}}, s::GBAllType) = s.p
Base.length(::GBAllType) = 0 #Allow indexing with ALL

function _load_globaltypes()
    global BOOL = GBType{Bool}("GrB_BOOL")
    ptrtogbtype[BOOL.p] = BOOL
    global INT8 = GBType{Int8}("GrB_INT8")
    ptrtogbtype[INT8.p] = INT8
    global UINT8 = GBType{UInt8}("GrB_UINT8")
    ptrtogbtype[UINT8.p] = UINT8
    global INT16 = GBType{Int16}("GrB_INT16")
    ptrtogbtype[INT16.p] = INT16
    global UINT16 = GBType{UInt16}("GrB_UINT16")
    ptrtogbtype[UINT16.p] = UINT16
    global INT32 = GBType{Int32}("GrB_INT32")
    ptrtogbtype[INT32.p] = INT32
    global UINT32 = GBType{UInt32}("GrB_UINT32")
    ptrtogbtype[UINT32.p] = UINT32
    global INT64 = GBType{Int64}("GrB_INT64")
    ptrtogbtype[INT64.p] = INT64
    global UINT64 = GBType{UInt64}("GrB_UINT64")
    ptrtogbtype[UINT64.p] = UINT64
    global FP32 = GBType{Float32}("GrB_FP32")
    ptrtogbtype[FP32.p] = FP32
    global FP64 = GBType{Float64}("GrB_FP64")
    ptrtogbtype[FP64.p] = FP64
    global FC32 = GBType{ComplexF32}("GxB_FC32")
    ptrtogbtype[FC32.p] = FC32
    global FC64 = GBType{ComplexF32}("GxB_FC64")
    ptrtogbtype[FC64.p] = FC64
    global NULL = GBType{Nothing}(C_NULL)
    ptrtogbtype[NULL.p] = NULL
    global ALL = GBAllType(load_global("GrB_ALL", libgb.GrB_Index))
    ptrtogbtype[ALL.p] = ALL
end

"""
    tojuliatype(x::GBType)

Determine the Julia equivalent of a GBType.

See also: [`toGBType`](@ref)
"""
tojuliatype(::GBType{T}) where {T} = T

"""
    toGBType(x)

Determine the GBType equivalent of a Julia primitive type.

See also: [`juliatype`](@ref)
"""
function toGBType(x)
    if x == Bool
        return BOOL
    elseif x == Int8
        return INT8
    elseif x == UInt8
        return UINT8
    elseif x == Int16
        return INT16
    elseif x == UInt16
        return UINT16
    elseif x == Int32
        return INT32
    elseif x == UInt32
        return UINT32
    elseif x == Int64
        return INT64
    elseif x == UInt64
        return UINT64
    elseif x == Float32
        return FP32
    elseif x == Float64
        return FP64
    elseif x == ComplexF32
        return FC32
    elseif x == ComplexF64
        return FC64
    else
        throw(ArgumentError("Not a valid GrB data type"))
    end
end


text:
function nelson_aalen(sort_ev::SortedEvents, Θ = ones(length(sort_ev.events)))
    fs, ls, events = sort_ev.fs, sort_ev.ls, sort_ev.events
    ds = ls-fs+1
    ns = after(Θ)[fs]
    chaz = cumsum(ds./ns)
    v = cumsum(ds./ns.^2)
    return [event.time for event in events[fs]], chaz, v
end

nelson_aalen(events::Array, Θ = ones(length(events))) =
nelson_aalen(SortedEvents(events), Θ)

if isdefined(:NullableArray)
    nelson_aalen(events::NullableArray, args...) = nelson_aalen(events.values[!events.isnull], args...)
end
nelson_aalen(events::AbstractArray, args...) = nelson_aalen(convert(Array,events), args...)

function nelson_aalen(res::CoxModel)
    nelson_aalen(convert(Array,res.M.df[:,1]),
    exp.(((DataFrames.ModelMatrix(res.M).m)[:,2:end])*(res.coefmat.cols[1])))
end

function survival(res::CoxModel)
    x_chaz,y_chaz,v_chaz = nelson_aalen(res::CoxModel)
    x, y, v = x_chaz, exp.(-y_chaz), exp.(-2y_chaz) .* v_chaz
    return x, y, v
end

chaz2haz(x,y) = (x, diff(vcat(0.,y)))

function chaz2haz(x,y,bw; npoints = 1000)
    pts = linspace(0., maximum(x), npoints)
    ycont = zeros(size(pts))
    index = vcat(0,round.([Int64],x/step(pts)))+1
    yext = vcat(0.,y)
    for i in 1:length(x)
        ycont[index[i]:index[i+1]] = yext[i]
    end
    dist = Normal(0.,bw)
    smoother = pdf(dist, quantile(dist,0.001):step(pts):quantile(dist,0.999))
    append!(ycont,fill(y[end],div(length(smoother),2)+1))
    ysmooth = conv(ycont,smoother)[div(length(smoother),2):(length(pts)+div(length(smoother),2))]
    smoothaz = diff(ysmooth)
    return pts, smoothaz
end


text:
type CondEdge
    vertices::Array{Int, 1}
    condSet::Array{Int, 1}
    paths::Array{Array{Int, 1}}
    parents::Array{Int, 1}
end

function CondEdge(v1::Int, v2::Int, condSet1::Array{Int, 1},
                  condSet2::Array{Int, 1})
    condSet = sort(condSet1)
    if v1 < v2
        parents = [condSet1[end], condSet2[end]]
        paths = Array{Int, 1}[condSet1, condSet2]
        vertices = [v1, v2]
    else
        parents = [condSet2[end], condSet1[end]]
        paths = Array{Int, 1}[condSet2, condSet1]
        vertices = [v2, v1]
    end
    return CondEdge(vertices, condSet, paths, parents)
end

import Base.Multimedia.display
function display(e::CondEdge)
    condSetString = "$(e.condSet[1])"
    for ii=2:length(e.condSet)
        condSetString = join([condSetString, ",", e.condSet[ii]])
    end
    println("     $(e.vertices[1]),$(e.vertices[2]) | $condSetString")
end

function display(e::Array{CondEdge, 1})
    nEdges = length(e)
    for ii=1:nEdges
        display(e[ii])
    end
end

function ==(e1::CondEdge, e2::CondEdge)
    equVertices = (e1.vertices == e2.vertices)
    equCondSets = (e1.paths == e2.paths)
    return  equVertices & equCondSets
end

function vertexTable(edgeArr::Array{CondEdge, 1}, nVars::Int)
    vTable = Array(Array{Int, 1}, nVars)
    nEdges = length(edgeArr)
    for ii=1:nEdges
        v1, v2 = edgeArr[ii].vertices
        if isdefined(vTable, v1)
            push!(vTable[v1], ii)
        else
            vTable[v1] = [ii]
        end
        if isdefined(vTable, v2)
            push!(vTable[v2], ii)
        else
            vTable[v2] = [ii]
        end
    end
    return vTable
end


text:
# read natural images

# cut patches

# normalize them

# divide on and off

# apply NNSC

# visualize results


text:
using ProximalOperators
using Test

@testset "PointwiseMinimum" begin

T = Float64

f = PointwiseMinimum(IndPoint(T[-1.0]), IndPoint(T[1.0]))
x = T[0.1]

predicates_test(f)
@test ProximalOperators.is_set(f) == true
@test ProximalOperators.is_cone(f) == false

y, fy = prox_test(f, x)
@test all(y .== T[1.0])
@test fy == T(0)

end


text:
Base.__precompile__(true)

module Delta

using Statistics: mean, std
import Statistics.mean

include("./kit/0-kit.jl")
include("./basic/0-basic.jl")
include("./block/0-block.jl")
include("./loss/0-loss.jl")
include("./optimizer/optimizer.jl")
include("./optimizer/update.jl")
include("./optimizer/regularize.jl")
include("./normalizer/normalizers.jl")
include("./scaler/Scalers.jl")

end  # moduleDelta

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using ML101
using Test

using DataFrames

@testset "ML101.jl" begin

    include("linear_regression.jl")

end


text:
function solve(solver::QMDPSolver, pomdp::POMDP; sparse::Bool=false)
    println("solving POMDP now, please wait...")
    if sparse
        vi_solver = SparseValueIterationSolver(max_iterations=solver.max_iterations, belres=solver.tolerance, verbose=solver.verbose, include_Q=true)
        vi_policy = solve(vi_solver, pomdp)
    else
        vi_solver = ValueIterationSolver(max_iterations=solver.max_iterations, belres=solver.tolerance, verbose=solver.verbose, include_Q=true)
        vi_policy = solve(vi_solver, pomdp)
    end

    return AlphaVectorPolicy(pomdp, vi_policy.qmat)
end


text:
module LibPQ

export status,
    reset!,
    execute,
    prepare,
    async_execute,
    cancel,
    num_columns,
    num_rows,
    num_params,
    num_affected_rows

using Base: Semaphore, acquire, release
using Base.Iterators: zip, product
using Base.Threads

using Dates
using DocStringExtensions
using Decimals
using FileWatching
using Tables
using Infinity: InfExtendedTime, isposinf, ∞
using Intervals
using IterTools: imap
using LayerDicts
using Memento: Memento, getlogger, warn, info, error, debug
using OffsetArrays
using SQLStrings
using TimeZones

const Parameter = Union{String,Missing}
const LOGGER = getlogger(@__MODULE__)

function __init__()
    INTERVAL_REGEX[] = _interval_regex()
    Memento.register(LOGGER)
    return nothing
end

# Docstring template for types using DocStringExtensions
@template TYPES = """
        $(TYPEDEF)

    $(DOCSTRING)

    ## Fields:

    $(TYPEDFIELDS)
    """

include(joinpath(@__DIR__, "utils.jl"))

module libpq_c
    export Oid

    using LibPQ_jll

    include(joinpath(@__DIR__, "headers", "libpq-fe.jl"))
end

using .libpq_c

include("typemaps.jl")

const DEFAULT_CLIENT_TIME_ZONE = Ref("UTC")

"""
    const LIBPQ_TYPE_MAP::PQTypeMap

The [`PQTypeMap`](@ref) containing LibPQ-level type mappings for LibPQ.jl.
Adding type mappings to this constant will override the default type mappings for all code
using LibPQ.jl.
"""
const LIBPQ_TYPE_MAP = PQTypeMap()

"""
    const LIBPQ_CONVERSIONS::PQConversions

The [`PQConversions`](@ref) containing LibPQ-level conversion functions for LibPQ.jl.
Adding conversions to this constant will override the default conversions for all code using
LibPQ.jl.
"""
const LIBPQ_CONVERSIONS = PQConversions()

const BINARY = true
const TEXT = false

include("connections.jl")
include("results.jl")
include("statements.jl")
include("exceptions.jl")

include("parsing.jl")
include("copy.jl")
include("tables.jl")

include("asyncresults.jl")

include("deprecated.jl")

end


text:
power_law_params = Dict(
    "A" => (α=0.108),
    "B" => (α=0.112),
    "C" => (α=0.120),
    "D" => (α=0.142),
    "E" => (α=0.203),
    "F" => (α=0.253)
)

"""
    windspeed(u0, z0, stability_class)
returns the windspeed function u(z) for a given Pasquill-Gifford
stability class, `z` is assumed to be in meters and `u` is in m/s

# Arguments
`u0` reference windspeed
`z0` height of reference point
`stability_class` Pasquill stability class (A, B, C, D, E, F)

# References

"""
function windspeed(u0::Number, z0::Number, stability_class::String)
    if stability_class ∈ Set(["A","B","C","D","E","F"])
        α, = power_law_params[stability_class]
    else
        err = "$stability_class is not a valid Pasquill-Gifford stability class"
        error(err)
    end

    return z -> u0*(z/z0)^α
end

"""
    windspeed(u, zR, λ, stability_class; k=0.35)
returns the windspeed function u(z) for a given Pasquill-Gifford
stability class, `z` is assumed to be in meters and `u` is in m/s

# Arguments
`u` friction velocity
`zR` surface roughness
`λ`  Monin-Obukhov length
`stability_class` Pasquill stability class (A, B, C, D, E, F)
`k`  von Karman's constant, 0.35

# References
    Businger et al. 1971
"""
function windspeed(u::Number, zR::Number, λ::Number, stability_class::String; k=0.35)
    if stability_class ∈ Set(["A","B","C"])
        a(z) = (1-15*(z/λ))^0.25
        Ψ(a) = 2*log((1+a)/2) + log((1+a^2)/2) - 2*atan(a) + π/2
        return z -> (u/k)*(log((z+zR)/zR) - Ψ(a(z)))
    elseif stability_class == "D"
        return z -> (u/k)*log((z+zR)/zR)
    elseif stability_class ∈ Set(["E", "F"])
        return z -> (u/k)*(log((z+zR)/zR) - 4.7*(z/λ))
    else
        err = "$stability_class is not a valid Pasquill-Gifford stability class"
        error(err)
    end
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

name = "cbindgen"
version = v"0.22.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/eqrion/cbindgen/archive/refs/tags/v$(version).tar.gz",
                  "f129b453df9d84e6d098a446f928961241b2a0edc29f827addca154049dcc434"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cbindgen*/
cargo build --release
install -Dvm 755 "target/${rust_target}/release/cbindgen${exeext}" "${bindir}/cbindgen${exeext}"
"""

# Our Rust toolchain for i686 is unusable.
platforms = supported_platforms(; exclude=p -> Sys.iswindows(p) && arch(p) == "i686")

# The products that we will ensure are always built
products = [
    ExecutableProduct("cbindgen", :cbindgen),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", compilers=[:c, :rust])

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
"""
Package: Forecast

    nextodd(x)

Return the smallest odd integer greater than or equal to `x`.
"""
function nextodd(x::Real)::Integer
    cx = Integer(ceil(x))
    mod(cx,2)==0 ? cx+1 : cx
end


text:
# Arrays in Julia

# declaring an Array
julia> Array{Union{Nothing, String}}(nothing, 2)
2-element Array{Union{Nothing, String},1}:
 nothing
 nothing

# converting string to Array
julia> a = "praveen"

# collect is similar to list operation
# on strings in Python
julia> l = collect(a)
7-element Array{Char,1}:
 'p'
 'r'
 'a'
 'v'
 'e'
 'e'
 'n'

# type check
julia> typeof(l)
Array{Char,1}

# display Array
julia> l
7-element Array{Char,1}:
 'p'
 'r'
 'a'
 'v'
 'e'
 'e'
 'n'


text:
module TestMixtape

using Test
using Mixtape
import Mixtape: CompilationContext,
                allow,
                transform,
                optimize!
using CodeInfoTools
using MacroTools

EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")

for ex in readdir(EXAMPLES_DIR)
    @testset "$ex" begin
        # Just check the examples run, for now.
        include(joinpath(EXAMPLES_DIR, ex))
    end
end

end # module


text:
module Dynamics
    using DifferentialEquations
    using DiffEqSensitivity
    using ForwardDiff
    using DiffResults
    using LinearAlgebra
    import StaticArrays
    import SymEngine
    import ..SymbolicUtils
    using ..RocketlandDefns
    using ..Aerodynamics

    const mass_idx = 1
    const r_idx = 2:4 # SVector(2,3,4)
    const v_idx = 5:7 # SVector(5,6,7)
    const qbi_idx = 8:11 # SVector(8,9,10,11)
    const omb_idx = 12:14 # SVector(12,13,14)

    const thr_idx = 1:3# SVector(1,2,3)

    struct IntegratorParameters{T}
        dt :: Float64
        uk :: Matrix{T}
        up :: Matrix{T}
        sigma :: T
        pinfo :: ProbInfo
    end

    @inline function DCM(quat::AbstractArray{T}) where T
        q0 = quat[1]
        q1 = quat[2]
        q2 = quat[3]
        q3 = quat[4]
        p1 = q1 * q2
        p2 = q0 * q3
        p3 = q1 * q3
        p4 = q0 * q2
        p5 = q2 * q3
        p6 = q0 * q1

        return  StaticArrays.@SMatrix([(1-2*(q2^2 + q3^2)) (2*(p1 - p2))       (2*(p3 + p4));
                        (2*(p1 + p2))        (1-2*(q1^2 + q3^2)) (2*(p5 - p6));
                        (2*(p3 - p4))        (2*(p5 + p6))       (1-2*(q1^2 + q2^2)) ])
    end

    @inline function Omega(omegab::AbstractArray{T}) where T
        z = zero(T)
        return StaticArrays.@SMatrix([z        -omegab[1] -omegab[2] -omegab[3];
                       omegab[1] z           omegab[3] -omegab[2];
                       omegab[2] -omegab[3]  z          omegab[1];
                       omegab[3]  omegab[2] -omegab[1]  z         ])
    end

    @inline function dx_static(state::AbstractArray{T} where T, u::AbstractArray{T} where T, mult, info::ProbInfo)
        qbi = @view state[qbi_idx]
        omb = @view state[omb_idx]

        aerf, bdy_trq = Aerodynamics.aero_force(info.aero, DCM(qbi) * [1.0,0.0,0.0], (@view state[5:7]),info.sos)

        fd1 = cross(DCM(qbi) * [0.0,1.0,0.0], state[5:7])
        fd1 = fd1/sqrt(sum(fd1 .* fd1))
        fd2 = cross(fd1, state[5:7])
        #= ff = u[4] * fd1 + u[5] * fd2 =#

        thr_frc = DCM(qbi) * (u[thr_idx])
        aero_frc = aerf #= + ff =#
        acc = (thr_frc + aero_frc) ./ state[mass_idx]
        rot_vel = 0.5*Omega(omb)*qbi
        aero_trq = [0.0,0.0,0.0]#= cross(info.rFB, ff) + bdy_trq=#
        rot_acc = info.jBi*(cross(info.rTB,u[thr_idx]) + aero_trq - cross(omb,info.jB*omb))
        return StaticArrays.@SVector([-info.a*sqrt(sum(u[thr_idx] .* u[thr_idx])),
                     state[5],state[6],state[7],
                     acc[1]-info.g0, acc[2], acc[3],
                     rot_vel[1], rot_vel[2], rot_vel[3], rot_vel[4],
                     rot_acc[1], rot_acc[2], rot_acc[3]]) .* mult

    end
    @inline function dx_static(state::AbstractArray{SymEngine.Basic} where T, u::AbstractArray{SymEngine.Basic} where T, mult, info::ProbInfo)
        qbi = @view state[qbi_idx]
        omb = @view state[omb_idx]

        aerf, bdy_trq = Aerodynamics.aero_force(info.aero, DCM(qbi) * [1.0,0.0,0.0], (@view state[5:7]),info.sos)

        fd1 = cross(DCM(qbi) * [0.0,1.0,0.0], state[5:7])
        fd1 = fd1/sqrt(sum(fd1 .* fd1))
        fd2 = cross(fd1, state[5:7])
        #= ff = u[4] * fd1 + u[5] * fd2 =#

        thr_frc = DCM(qbi) * (u[thr_idx])
        aero_frc = aerf #= + ff =#
        acc = (thr_frc + aero_frc) ./ state[mass_idx]
        rot_vel = 0.5*Omega(omb)*qbi
        aero_trq = [0.0,0.0,0.0]#= cross(info.rFB, ff) + bdy_trq =#
        rot_acc = info.jBi*(cross(info.rTB,u[thr_idx]) - cross(omb,info.jB*omb))
        return [-info.a*sqrt(sum(u[thr_idx] .* u[thr_idx])),
                     state[5],state[6],state[7],
                     acc[1]-info.g0, acc[2], acc[3],
                     rot_vel[1], rot_vel[2], rot_vel[3], rot_vel[4],
                     rot_acc[1], rot_acc[2], rot_acc[3]] .* mult

    end

    @inline function dx(output, state::AbstractArray{T} where T, u::AbstractArray{T} where T, mult, info::ProbInfo)
        output[:] = dx_static(state, u, mult, info)
        0
    end

    @inline function current_control(pc, start_ctrl::AbstractArray{T} where T, ed_ctrl::AbstractArray{T} where T)
        return (1.0-pc) * start_ctrl + pc * ed_ctrl
    end

    function rk4(inp::AbstractArray{T}, dt, info; npts=10) where T
        state = collect(@view inp[1:14] )
        start_ctrl = (@view inp[15:17])
        ed_ctrl = (@view inp[18:20])
        i_sigma = dt
        idt = i_sigma/(npts)
        pcs = 1.0/(npts)
        ct = 0.0
        pca = 0.0
        for i=1:npts
            ict = current_control(pca, start_ctrl, ed_ctrl)
            mct = current_control(pca + pcs/2, start_ctrl, ed_ctrl)
            ect = current_control(pca + pcs, start_ctrl, ed_ctrl)
            k1 = dx_static(state, ict, inp[21], info)
            k2 = dx_static(state .+ k1./2, mct, inp[21], info)
            k3 = dx_static(state .+ k2./2, mct, inp[21], info)
            k4 = dx_static(state .+ k3, ect, inp[21], info)
            ct += idt
            pca += pcs
            state = state .+ idt .* (k1./6 .+ k2./3 .+ k3./3 .+ k4./6)
        end
        return state
    end

    const stateC = 1:14
    const ukC = 15:17
    const upC = 18:20
    const sigmaC = 21

    function make_dynamics_module(info::ProbInfo)
        t=SymEngine.symbols("t")
        dt=SymEngine.symbols("dt")
        lkm = (dt-t)/dt
        lkp = t/dt

        #sig: dx(J, inp, dt, pinfo, t)
        dx_fun = SymbolicUtils.make_simplified(:dx, st -> dx_static(st[stateC], st[ukC]*lkm + st[upC]*lkp, st[sigmaC], info), 21)
        dxa_fun = SymbolicUtils.make_simplified_array(:dxa, st -> dx_static(st[stateC], st[ukC]*lkm + st[upC]*lkp, st[sigmaC], info), 21)
        ddx_fun = SymbolicUtils.make_jacobian(:ddx, st -> dx_static(st[stateC], st[ukC]*lkm + st[upC]*lkp, st[sigmaC], info), 21)
        #J_fun = SymbolicUtils.make_jacobian(:jac, st -> rk4(st, 1.0, info), 21)
        genmod = quote
            module Linearizer
                using LinearAlgebra
                using ForwardDiff
                using StaticArrays
                using ..RocketlandDefns
                using ..Aerodynamics

                # all of these are called by the generated code
                function clamp_aoa(cos_aoa::T, mach::T,p::ProbInfo) where T
                    int_aoa = zero(T)
                    if mach > zero(T)
                        int_aoa = clamp(cos_aoa/(mach*p.sos),-1.0,1.0)
                    end
                    return int_aoa
                end

                function drag(cos_aoa::T,mach::T,p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    dragf = direct_drag(p.aero,int_aoa,mach)::T
                    return dragf
                end
                function drag(::Type{Val{:jac}}, cos_aoa::T,mach::T, p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    return direct_drag(Val{:jac}, p.aero, int_aoa, mach)
                end

                function lift(cos_aoa::T,mach::T,p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    return direct_lift(p.aero,int_aoa,mach)::T
                end
                function lift(::Type{Val{:jac}}, cos_aoa::T,mach::T,p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    return direct_lift(Val{:jac}, p.aero, int_aoa, mach)
                end

                function trq(cos_aoa::T,mach::T,p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    return direct_trq(p.aero,int_aoa,mach)::T
                end
                function trq(::Type{Val{:jac}}, cos_aoa::T,mach::T,p::ProbInfo) where T
                    int_aoa = clamp_aoa(cos_aoa, mach, p)
                    return direct_trq(Val{:jac}, p.aero, int_aoa, mach)
                end

                @inline function ifnz(val::T, nz::V) where {T,V}
                    if !iszero(val) && !isnan(val)
                        return nz
                    else
                        return zero(V)
                    end
                end
                @inline function ifnz(::Type{Val{:jac}}, val::T, nz::V) where {T,V}
                    return StaticArrays.SArray{Tuple{2}}(0.0,1.0) # blunt approximation
                end

                $dx_fun
                $ddx_fun
                $dxa_fun
            end
        end
        Main.eval(genmod.args[2])
    end

    struct DxIntegrator{T,V,X} <: Function
        temp_state::T
        dt::Float64
        info::ProbInfo{V}
        dx::X
    end
    @inline function (ff::DxIntegrator)(du,u,p::Vector{Float64},t::Float64)
        fill!(du, zero(eltype(du)))
        ff.temp_state .= u .+ p
        res = ff.dx(du, ff.temp_state, ff.dt, ff.info, t)
        return nothing
    end

    @inline function ddx_integrator(J,u,p,t, dt, info, ddx)
        fill!(J, zero(eltype(J)))
        for i=1:21
            J[i,i] = 1.0
        end
        return nothing
    end

    struct ParamJac{T,F,F2} <: Function
        temp_state :: Vector{Float64}
        dt::Float64
        info::ProbInfo{T}
        ddx::F
        dxa::F2
    end
    @inline function (ff::ParamJac)(pJ::Matrix{Float64}, u, p::Vector{Float64}, t::Float64)
        fill!(pJ, zero(eltype(pJ)))
        ff.temp_state .= u .+ p
        lkm = (ff.dt-t)/ff.dt
        lkp = t/ff.dt
        # dx_static(inp[1:14], inp[15:17]*lkm + inp[18:20]*lkp, inp[21], ff.info)
        jacres = Zygote.forward_jacobian(inp->ff.dxa(inp, ff.dt, ff.info, t), StaticArrays.SVector{21}(ff.temp_state))
        pJ[1:14,1:21] .= transpose(jacres[2])
        #res = ff.ddx(pJ, ff.temp_state, ff.dt, ff.info, t)
        #println("inp $(ff.temp_state) outp $pJ")
        return nothing
    end

    function (::Type{IntegratorCache})(prob::DescentProblem, info::ProbInfo, lin_mod)
        return IntegratorCache([prob.mwet; prob.rIi; prob.vIi; prob.qBIi; prob.wBi; 1.0; 0.0; 0.0; 1.0; 0.0; 0.0; 1.0], 1.0, info, lin_mod)
    end
    function (::Type{IntegratorCache})(inp::Vector{Float64}, dt::Float64, info::ProbInfo, lin_mod)
        params = [dt,info]
        dxfun = lin_mod.dx
        dxafun = lin_mod.dxa
        ddxfun = lin_mod.ddx
        dx_int = DxIntegrator(zeros(21), params[1], params[2], dxfun)
        pj_int = ParamJac(zeros(21), params[1], params[2], ddxfun, dxafun)
        fun = ODEFunction(
            dx_int,
            jac=(J,u,p,t)->fill!(J,zero(eltype(J))),
            paramjac=pj_int)
        sim_prob = ODEProblem(fun,
            zeros(21), convert.(eltype(inp), (0.0,dt)),
            inp,
            autodiff=false,
            autojacvec=false)
        sim_int = init(sim_prob, BS3(), abstol=1e-2, reltol=1e-2, save_everystep=false, save_start=false)
        sense_prob = ODEForwardSensitivityProblem(fun,
            zeros(21), convert.(eltype(inp), (0.0,dt)),
            inp,
            ForwardSensitivity(autojacvec=false);
            autodiff=false,
            autojacvec=false)
        sense_int = init(sense_prob, BS3(), abstol=1e-2, reltol=1e-2, save_everystep=false, save_start=false)
        IntegratorCache(sim_prob, sense_prob, sim_int, sense_int, params, info)
    end

    function simulate(inp::Vector{Float64}, dt::Float64, cache::IntegratorCache)
        cache.sim_int = init(cache.sim_prob, BS3(), abstol=1e-2, reltol=1e-2, save_everystep=false, save_start=false)
        fill!(cache.sim_int.u, 0.0)
        reinit!(cache.sim_int, zeros(14), t0 = convert(eltype(inp), 0.0), tf = convert(eltype(inp), dt), reset_dt=true)
        cache.sim_prob.p .= inp
        cache.params[1] = dt
        solve!(cache.sim_int)
        return @view (cache.sim_int.u .+ inp)[1:14]
    end

    function sensitivity(inp::Vector{Float64}, dt::Float64, cache::IntegratorCache)
        fill!(cache.sense_int.u, 0.0)
        reinit!(cache.sense_int, zeros(462), t0 = convert(eltype(inp), 0.0), tf = convert(eltype(inp), dt), reset_dt=true)
        cache.sense_prob.p .= inp
        cache.params[1] = dt
        solve!(cache.sense_int)
        return extract_local_sensitivities(cache.sense_int.sol, 1, Val(true))
    end

    import Zygote
    function simulate_zygote(inp::Vector{Float64}, dt::Float64, cache::IntegratorCache; npts=10)
        return rk4(inp, dt, cache.info, npts=npts)
    end
    function sensitivity_zygote(inp::Vector{Float64}, dt::Float64, cache::IntegratorCache)
        return Zygote.forward_jacobian(inp->rk4(inp, dt, cache.info), StaticArrays.SVector{21}(inp))
    end

    function predict_state(initial_state, uk, up, sigma, dt, pinfo, cache)
        return simulate([initial_state; uk; up; sigma], dt, cache)
    end
    function make_state(a::LinPoint, b::LinPoint, sig::Float64)
        return vcat(a.state, a.control, b.control, sig)
    end
    function linearize_dynamics(states::Array{LinPoint,1}, tf_guess::Float64, base_dt::Float64, cache::IntegratorCache) where {N,K}
        results = Array{LinRes,1}(undef, length(states)-1)

        for i=1:length(states)-1
            ist = make_state(states[i], states[i+1], tf_guess)
            val, mat = sensitivity(ist, base_dt, cache)
            val[1:14] .+= states[i].state
            for j=1:14
                mat[j,j] += 1
            end
            results[i] = LinRes(copy(val[1:14]), copy(mat[1:14,1:21]))
        end
        return results
    end
    export linearize_dynamics, next_step
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
# Autogenerated wrapper script for JonkerVolgenant_jll for i686-w64-mingw32
export bipartite_assignement

JLLWrappers.@generate_wrapper_header("JonkerVolgenant")
JLLWrappers.@declare_library_product(bipartite_assignement, "bipartite_assignement.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        bipartite_assignement,
        "bin\\bipartite_assignement.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
immutable VLineGeometry <: Gadfly.GeometryElement
    color::Union(ColorValue, Nothing)
    size::Union(Measure, Nothing)

    function VLineGeometry(; color=nothing,
                           size::Union(Measure, Nothing)=nothing)
        new(color === nothing ? nothing : Color.color(color), size)
    end
end

const vline = VLineGeometry

    function element_aesthetics(::VLineGeometry)
    [:xintercept]
end

# Generate a form for the vline geometry
function render(geom::VLineGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics, data::Gadfly.Data,
                scales::Dict{Symbol, ScaleElement})
    Gadfly.assert_aesthetics_defined("Geom.vline", aes, :xintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size

    return compose!(
        context(),
        Compose.line([{(x, 0h), (x, 1h)} for x in aes.xintercept]),
        stroke(color),
        linewidth(size),
        svgclass("yfixed"))
end


text:
function Base.iterate(
    layer::AbstractFeatureLayer,
    state::Integer = 0,
)::Union{Nothing,Tuple{Feature,Int64}}
    layer.ptr == C_NULL && return nothing
    state == 0 && resetreading!(layer)
    ptr = GDAL.ogr_l_getnextfeature(layer.ptr)
    return if ptr == C_NULL
        resetreading!(layer)
        nothing
    else
        (Feature(ptr), state + 1)
    end
end

Base.eltype(layer::AbstractFeatureLayer)::DataType = Feature

Base.IteratorSize(::Type{<:AbstractFeatureLayer}) = Base.SizeUnknown()

Base.length(layer::AbstractFeatureLayer)::Integer = nfeature(layer, true)

struct BlockIterator{T<:Integer}
    rows::T
    cols::T
    ni::T
    nj::T
    n::T
    xbsize::T
    ybsize::T
end

function blocks(
    ::Type{T},
    raster::AbstractRasterBand,
)::BlockIterator{T} where {T<:Integer}
    (xbsize, ybsize) = blocksize(raster)
    rows = T(height(raster))
    cols = T(width(raster))
    ni = ceil(T, rows / ybsize)
    nj = ceil(T, cols / xbsize)
    return BlockIterator{T}(rows, cols, ni, nj, ni * nj, xbsize, ybsize)
end

function blocks(raster::AbstractRasterBand)::BlockIterator{Int64}
    return blocks(Int64, raster)
end

function Base.iterate(
    obj::BlockIterator{T},
    iter::T = 0,
)::Union{Nothing,Tuple{Tuple{Tuple{T,T},Tuple{T,T}},T}} where {T}
    iter == obj.n && return nothing
    j = floor(Int, iter / obj.ni)
    i = iter % obj.ni
    nrows = if (i + 1) * obj.ybsize < obj.rows
        obj.ybsize
    else
        obj.rows - i * obj.ybsize
    end
    ncols = if (j + 1) * obj.xbsize < obj.cols
        obj.xbsize
    else
        obj.cols - j * obj.xbsize
    end
    return (((i, j), (nrows, ncols)), iter + 1)
end

struct WindowIterator{T<:Integer}
    blockiter::BlockIterator{T}
end
Base.size(i::WindowIterator) = (i.blockiter.ni, i.blockiter.nj)
Base.length(i::WindowIterator) = i.blockiter.n

function Base.IteratorSize(::Type{WindowIterator{T}}) where {T<:Integer}
    return Base.HasShape{2}()
end

function Base.IteratorEltype(::Type{WindowIterator{T}}) where {T<:Integer}
    return Base.HasEltype()
end

function Base.eltype(::WindowIterator{T})::DataType where {T<:Integer}
    return Tuple{UnitRange{T},UnitRange{T}}
end

function windows(
    ::Type{T},
    raster::AbstractRasterBand,
)::WindowIterator{T} where {T<:Integer}
    return WindowIterator{T}(blocks(T, raster))
end

windows(raster::AbstractRasterBand)::WindowIterator{Int64} =
    windows(Int64, raster)

function Base.iterate(
    obj::WindowIterator{T},
    iter::T = 0,
)::Union{Nothing,Tuple{NTuple{2,UnitRange{T}},T}} where {T<:Integer}
    handle = obj.blockiter
    next = Base.iterate(handle, iter)
    next == nothing && return nothing
    (((i, j), (nrows, ncols)), iter) = next
    return (
        ((1:ncols) .+ j * handle.xbsize, (1:nrows) .+ i * handle.ybsize),
        iter,
    )
end

mutable struct BufferIterator{R<:Real,T<:Integer}
    raster::AbstractRasterBand
    w::WindowIterator{T}
    buffer::Matrix{R}
end

function bufferwindows(
    ::Type{T},
    raster::AbstractRasterBand,
)::BufferIterator{pixeltype(raster),T} where {T<:Integer}
    return BufferIterator{pixeltype(raster),T}(
        raster,
        windows(T, raster),
        Matrix{pixeltype(raster)}(undef, blocksize(raster)...),
    )
end

function bufferwindows(
    raster::AbstractRasterBand,
)::BufferIterator{pixeltype(raster),Int64}
    return bufferwindows(Int64, raster)
end

function Base.iterate(
    obj::BufferIterator{R,T},
    iter::T = 0,
)::Union{Nothing,Tuple{Matrix{R},T}} where {R<:Real,T<:Integer}
    next = Base.iterate(obj.w, iter)
    next == nothing && return nothing
    ((cols, rows), iter) = next
    rasterio!(obj.raster, obj.buffer, rows, cols)
    return (obj.buffer[1:length(cols), 1:length(rows)], iter)
end


text:
using Turing, Distributions, MCMCChains, KernelDensity,AdvancedHMC
using StatsBase, CSV, DataFrames, Plots, StatsPlots
using ReverseDiff
using Combinatorics
Turing.setadbackend(:reversediff)

# functions
function top(x,n = 3, test = false)
    """
    Generates list of every possible order of n picks from the neural net prbabilities
    e.g. the 3 most likely species (in order) for a given data point.
    Adds up the instances of each combination, in a given data set x.
    If x is training data subset for instances of a given species, then
    return proportions.
    If x is test data, then return the actual counts.
    """
top3 = []
conf = []
tops = collect(permutations([1,2,3,4,5],n))
for i in 1:size(x)[1]
    row = Array(x[i,1:end])

    push!(topx,sortperm(row, rev = true)[1:n])
    push!(conf,maximum(row))
end

    topxs = [sum([j == k for j in topx]) for k in tops]
    if test = true

        return topxs
    else
        return topxs ./ length(top3)
end

@model top3_model(a,o,p,s,v, y, n) = begin
        """
        a - v are the top n values from the 'top' function.
        y is the top n values for the test data.
        n is the number of observations in the test data
        """
    props ~ Dirichlet(5,1)
    x = (a .* props[1] .* n) .+ (o .* props[2] .* n) .+ (p .* props[3] .* n) .+ (s .* props[4] .* n) .+ (v .* props[5] .* n)
    y .~ Poisson.(x)
end;

y5 = CSV.read("5y.csv")[1:end,3:end]# actualy values
names!(y5,[:aa,:od,:pc,:sr,:vg])
pred5 = CSV.read("5pred.csv")[1:end,3:end] # y value predictions (probabilities) from neural net
names!(pred5,[:aa,:od,:pc,:sr,:vg])

n = size(y5)[1] #number of observations

# split data into training and test sets
train_pred5 = pred5[1:n÷2,1:end]
train_y5 = y5[1:n÷2,1:end]
test_pred5 = pred5[n÷2:end,1:end]
test_y5= y5[n÷2:end,1:end]

choices = train_pred5[train_y5[:aa] .== 1,1:end]
aatop3 = top(choices,3)
choices = train_pred5[train_y5[:od] .== 1,1:end]
odtop3 = top(choices,3)
choices = train_pred5[train_y5[:pc] .== 1,1:end]
pctop3 = top(choices,3)
choices = train_pred5[train_y5[:sr] .== 1,1:end]
srtop3 = top(choices,3)
choices = train_pred5[train_y5[:vg] .== 1,1:end]
vgtop3 = top(choices,3)

#r = sample(1:size(test_y5)[1],Weights(1.0 .- Array(test_y5[1:end,1]) ), n÷4)
r = sample(1:size(test_y5)[1], n÷4)
pred_sub = test_pred5[r,1:end]
y_sub5 = test_y5[r,1:end]
y_vals1 = top(pred_sub,1,test = true)
y_vals2 = top(pred_sub,2,test = true)
y_vals3 = top(pred_sub,3, test = true)

# draw samples
top_chain1 = sample(top3_model(aatop,odtop,pctop,srtop,vgtop, y_vals1, n÷4), MH(), 1000000)
top_chain2 = sample(top3_model(aatop2,odtop2,pctop2,srtop2,vgtop2, y_vals2, n÷4), MH(), 1000000)
top_chain3 = sample(top3_model(aatop3,odtop3,pctop3,srtop3,vgtop3, y_vals3, n÷4), MH(), 1000000)

# plotting
density(top_chain1["props[1]"], label = "aa1")
density!(top_chain2["props[1]"], label = "aa2")
density!(top_chain3["props[1]"], label = "aa3")
vline!([sum(y_sub5[:aa])/(n÷4)], label = "aa")

density(top_chain1["props[2]"], label = "od1")
density!(top_chain2["props[2]"], label = "od2")
density!(top_chain3["props[2]"], label = "od3")
vline!([sum(y_sub5[:od])/(n÷4)], label = "od")

density(top_chain1["props[3]"], label = "pc1")
density!(top_chain2["props[3]"], label = "pc2")
density!(top_chain3["props[3]"], label = "pc3")
vline!([sum(y_sub5[:pc])/(n÷4)], label = "pc")

density(top_chain1["props[4]"], label = "sr1")
density!(top_chain2["props[4]"], label = "sr2")
density!(top_chain3["props[4]"], label = "sr3")
vline!([sum(y_sub5[:sr])/(n÷4)], label = "sr")

density(top_chain1["props[5]"], label = "vg1")
density(top_chain2["props[5]"], label = "vg2")
density(top_chain3["props[5]"], label = "vg3")
vline!([sum(y_sub5[:vg])/(n÷4)], label = "vg")

max_stars_repo_name: EvoArt/MachineFlow

max_issues_repo_name: EvoArt/MachineFlow

max_forks_repo_name: EvoArt/MachineFlow


text:
abstract type Scaler{N} end

include("linear.jl")


text:
# convert MSEED to SAC
using SeisIO, Glob, Dates

BASEDIR = "/media/four/"
DATADIR = joinpath(BASEDIR,"LASSO/DL/")
SACDIR = joinpath(BASEDIR,"LASSO/SAC/")

if !isdir(SACDIR)
    mkpath(SACDIR)
end

files = glob("*",DATADIR)
cd(SACDIR)

for ii = 1:length(files)
    println("Writing file $(basename(files[ii])), $ii of $(length(files)) $(now())")
    S = read_data("mseed",files[ii])
    writesac(S)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Private functions and variables of the Jacchia-Roberts 1971 Atmospheric
#   Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Roberts, C. R (1971). An analytic model for upper atmosphere densities
#       based upon Jacchia's 1970 models.
#
#   [2] Jacchia, L. G (1970). New static models of the thermosphere and
#       exosphere with empirical temperature profiles. SAO Special Report #313.
#
#   [3] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [4] Long, A. C., Cappellari Jr., J. O., Velez, C. E., Fuchs, A. J (editors)
#       (1989). Goddard Trajectory Determination System (GTDS) Mathematical
#       Theory (Revision 1). FDD/552-89/0001 and CSC/TR-89/6001.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# NamedTuple that will be use when defining constants related to the species.
JR1971_NT{T} = NamedTuple{ (:N₂, :O₂, :O, :Ar, :He, :H), NTuple{6,T} }

"""
    JR1971_CONSTANTS{T}

Structure with the constants for the Jacchia-Roberts 1971 Atmospheric Model.

"""
@with_kw struct JR1971_CONSTANTS{T}
    # General constants
    # ==========================================================================

    Rstar::T = T(8.31432)     # Rstar is the universal gas-constant [joules/(K.mol)].
    Av::T    = T(6.022045e23) # Avogadro's constant (mks) [molecules/mol].
    Ra::T    = T(6356.766)    # Mean Earth radius [km].
    g₀::T    = T(9.80665)     # Gravity at sea level [m/s²].

    # Values at the sea level
    # ==========================================================================

    M₀::T = 28.960            # Mean molecular mass at sea level [g/mol].

    # Values at lower boundary
    # ==========================================================================

    T₁::T = T(183)            # Temperature at the lower bound [K].
    z₁::T = T(90)             # Altitude of the lower bound [km].
    M₁::T = T(28.82678)       # Mean molecular mass at `z₁` [g/mol].
    ρ₁::T = T(3.46e-9)        # Density at `z₁` [g/cm³].

    # Values at the second division
    # ==========================================================================

    z₂::T = T(100)

    # Values at the inflection point
    # ==========================================================================

    zx::T = T(125)           # Altitude of the inflection point [km].

    # Molecular mass [g/mol]
    # ==========================================================================

    Mi::JR1971_NT{T} = ( N₂ = 28.0134,
                         O₂ = 31.9988,
                         O  = 15.9994,
                         Ar = 39.9480,
                         He =  4.0026,
                         H  = 1.00797 )

    # Thermal diffusion coefficient
    # ==========================================================================

    αi::JR1971_NT{T} = ( N₂ =  0,
                         O₂ =  0,
                         O  =  0,
                         Ar =  0,
                         He = -0.38,
                         H  =  0 )

    # Constituent density * M₀ / ρ_₁₀₀ / Av
    # ==========================================================================

    μi::JR1971_NT{T} = ( N₂ = 0.78110,
                         O₂ = 0.161778,
                         O  = 0.095544,
                         Ar = 0.0093432,
                         He = 0.61471e-5,
                         H  = 0 )

    # Coefficients for series expansion
    # ==========================================================================

    Aa::NTuple{7,T} = ( -435093.363387,
                          28275.5646391,
                           -765.33466108,
                             11.043387545,
                             -0.08958790995,
                              0.00038737586,
                             -0.000000697444 )

    Ca::NTuple{5,T} = ( -89284375.0,
                          3542400.0,
                           -52687.5,
                              340.5,
                               -0.8 )

    la::NTuple{5,T} = (  0.1031445e+5,
                         0.2341230e+1,
                         0.1579202e-2,
                        -0.1252487e-5,
                         0.2462708e-9 )

    α::NTuple{6,T} = ( 3144902516.672729,
                       -123774885.4832917,
                          1816141.096520398,
                           -11403.31079489267,
                               24.36498612105595,
                                0.008957502869707995 )

    β::NTuple{6,T} = ( -52864482.17910969,
                          -16632.50847336828,
                              -1.308252378125,
                               0,
                               0,
                               0)

    ζ::NTuple{7,T} = (  0.1985549e-10,
                       -0.1833490e-14,
                        0.1711735e-17,
                       -0.1021474e-20,
                        0.3727894e-24,
                       -0.7734110e-28,
                        0.7026942e-32 )

    δij::JR1971_NT{NTuple{7,T}} =
        (
          N₂ = (  0.1093155e2   ,
                  0.1186783e-2  ,    # (1/K)
                 -0.1677341e-5  ,    # (1/K)²
                  0.1420228e-8  ,    # (1/K)³
                 -0.7139785e-12 ,    # (1/K)⁴
                  0.1969715e-15 ,    # (1/K)⁵
                 -0.2296182e-19 ),   # (1/K)⁶

          O₂ = (  0.9924237e1   ,
                  0.1600311e-2  ,    # (1/K)
                 -0.2274761e-5  ,    # (1/K)²
                  0.1938454e-8  ,    # (1/K)³
                 -0.9782183e-12 ,    # (1/K)⁴
                  0.2698450e-15 ,    # (1/K)⁵
                 -0.3131808e-19 ),   # (1/K)⁶

          O =  (  0.1097083e2   ,
                  0.6118742e-4  ,    # (1/K)
                 -0.1165003e-6  ,    # (1/K)²
                  0.9239354e-10 ,    # (1/K)³
                 -0.3490739e-13 ,    # (1/K)⁴
                  0.5116298e-17 ,    # (1/K)⁵
                  0.0           ),   # (1/K)⁶

          Ar = (  0.8049405e1   ,
                  0.2382822e-2  ,    # (1/K)
                 -0.3391366e-5  ,    # (1/K)²
                  0.2909714e-8  ,    # (1/K)³
                 -0.1481702e-11 ,    # (1/K)⁴
                  0.4127600e-15 ,    # (1/K)⁵
                 -0.4837461e-19 ),   # (1/K)⁶

          He = (  0.7646886e1   ,
                 -0.4383486e-3  ,    # (1/K)
                  0.4694319e-6  ,    # (1/K)²
                 -0.2894886e-9  ,    # (1/K)³
                  0.9451989e-13 ,    # (1/K)⁴
                 -0.1270838e-16 ,    # (1/K)⁵
                  0.0           ),   # (1/K)⁶

          H =  (
                 0.0,
                 0.0,                # (1/K)
                 0.0,                # (1/K)²
                 0.0,                # (1/K)³
                 0.0,                # (1/K)⁴
                 0.0,                # (1/K)⁵
                 0.0 )               # (1/K)⁶
        )
end

"""
    _jr1971_constants

Constants for the Jacchia-Roberts 1971 Atmospheric Model.

"""
const _jr1971_constants = JR1971_CONSTANTS{Float64}()

"""
    _jr1971_id

Index of the species for the Jacchia-Roberts 1971 Atmospheric Model.

"""
const _jr1971_id = ( N₂ = 1, O₂ = 2, O  = 3, Ar = 4, He = 5, H  = 6 )

################################################################################
#                              Private Functions
################################################################################

"""
    _jr1971_M(z::R) where R

Compute the mean molecular mass at altitude `z` [km] using the empirical profile
in eq. 1 [3,4].

"""
@inline function _jr1971_M(z::R) where R
    !(90 <= z <= 100.0) && @warn "The empirical model for the mean molecular mass is valid only for 90 <= z <= 100 km."
    @unpack Aa = _jr1971_constants
    @evalpoly(z, Aa[1], Aa[2], Aa[3], Aa[4], Aa[5], Aa[6], Aa[7])
end

"""
    _jr1971_roots(p::Polynomial{R}) where R

Compute the roots of the polynomial `p` necessary to compute the density below
125 km. It returns the value `r₁`, `r₂`, `x`, and `y`.

"""
function _jr1971_roots(p::Vector{R}) where R
    # Compute the roots with a first guess.
    r = roots(p, Complex{R}[166.10; 61.32; 9.91 + 1.31im; 9.91 - 1.31im];
              polish = true)

    # Take the real roots and sort them.
    r_real = sort( r[ @. abs(imag(r)) < 1e-10 ];
                   lt = (a,b) -> real(a) < real(b) )

    # There must be 2 real roots.
    (length(r_real) != 2) && error("Internal error, there are not two real roots!")

    r₁ = real(r_real[2])
    r₂ = real(r_real[1])

    # Take the complex roots and sort them.
    r_complex = sort( r[ @. abs(imag(r)) >= 1e-10 ];
                      lt = (a,b) -> imag(a) < imag(b) )

    # There must be 2 complex roots.
    (length(r_complex) != 2) && error("Internal error, there are not two complex roots!")

    x = real(r_complex[2])
    y = imag(r_complex[2])

    r₁, r₂, x, y
end

"""
    _jr1971_T(z::R, Tx::R, T∞::R) where R<:Number

Compute the temperature [K] at height `z` [km] given the temperature `Tx` [K] at
the inflection point, and the exospheric temperature `T∞` [K] according to the
theory of the model Jacchia-Roberts 1971 [1,3,4].

The inflection point is considered to by `z = 125 km`.

"""
function _jr1971_T(z::R, Tx::R, T∞::R) where R<:Number
    @unpack T₁, z₁, zx = _jr1971_constants

    # Check the parameters
    # ==========================================================================

    (z  < z₁) && error("The altitude must not be lower than $(z₁) km.")
    (T∞ < 0 ) && error("The exospheric temperature must be positive.")

    # Compute the temperature at desire altitude
    # ==========================================================================

    if z <= zx
        @unpack Ca = _jr1971_constants

        aux = @evalpoly(z, Ca[1], Ca[2], Ca[3], Ca[4], Ca[5])
        T   = Tx + (Tx - T₁)/35^4*aux
    else
        @unpack Ra, la = _jr1971_constants

        l = @evalpoly(T∞, la[1], la[2], la[3], la[4], la[5])
        T = T∞ - ( T∞ - Tx )*exp( -( Tx - T₁ )/( T∞ - Tx )*
                                   (  z - zx )/( zx - z₁ )*
                                   (    l    )/( Ra + z  ) )
    end

    T
end


text:
@defcomp climateregional begin
    regions = Index()
    inputtemp = Parameter(index=[time])
    bregtmp = Parameter(index=[regions])
    bregstmp = Parameter(index=[regions])
    scentemp = Parameter(index=[time, regions])

    temp = Variable(index=[time,regions])
    regtmp = Variable(index=[time,regions])
    regstmp = Variable(index=[time,regions])

    function run_timestep(p, v, d, t)

        for r in d.regions
            v.regtmp[t, r] = p.inputtemp[t] * p.bregtmp[r] + p.scentemp[t, r]
        end

        for r in d.regions
            v.temp[t, r] = v.regtmp[t, r] / p.bregtmp[r]
        end

        for r in d.regions
            v.regstmp[t, r] = p.inputtemp[t] * p.bregstmp[r] + p.scentemp[t, r]
        end
    end
end


text:
@testset "Block operations" begin
    n_t = 4
    full_size = (n_t, 4, 8, 1)
    block_size = (n_t, 4, 3, 1)
    padding = (0, 0, 2, 0)
    test_blocks = Blocks(full_size, block_size, (0, 0, 0, 0), (0, 0, 0, 0), padding);

    test_block_starts, test_block_ends = PaddedBlocks.starts_ends(test_blocks);
    @test all(size(test_block_starts) .== (1, 1, 2, 1, 4))
    @test all(test_block_starts[1, 1, 2, 1, :] .== [1, 1, 4, 1])
end


text:
let
    n = 50
    p = 3
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);

    K = Symmetric(tril(Ut'*Vt) + triu(Vt'*Ut,1))
    @test ( isapprox(egrss.full(Ut,Vt),K) )
    @test ( isapprox(egrss.full(Ut,Vt,1e-3),K+Diagonal(1e-3*ones(n))) )
    @test ( isapprox(egrss.full(Ut,Vt,1e-3*ones(n)),K+Diagonal(1e-3*ones(n))) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
using Revise
using JDF

a = gf()
b = Array(a.ORIGIN)
io = iow()
metadata = compress_then_write(b, io)
close(io)

io = ior()
oo = column_loader(String, io, metadata)
close(io)

all(b.==oo)

b = Array(a.FL_DATE)
io = iow()
metadata = compress_then_write(b, io)
close(io)

io = ior()
oo = column_loader(metadata.type, io, metadata)
close(io)

all(b.==oo)

@time metadatas = savejdf(a, "c:/data/a.jdf")

@time a2 = loadjdf("c:/data/a.jdf", metadatas)

using LambdaFn
fn = @λ joinpath.(_1, readdir(_1))
files = fn("d:/data/AirOnTimeCSV")

sort!(files, by = x->filesize(x))

@time listdf =reduce(vcat, CSV.read(f) for f in files[1:50]);

using DataFrames, Blosc
@time schema = savejdf(listdf, "d:/data/plsdel.jdf");
@time loadjdf("d:/data/plsdel.jdf", schema);

	bb = [Blosc.compress(Array(b)) for b in eachcol(a)]
	bbc = bb[1])

	bbc = rand(UInt8, 100)
	io = open("io.test", "w")

	w  = Array{Float32,2}(undef,n1,nw)
	read!(fileID, w[:,iw])
	write(io, bbc)

	close(io)

	io = open("io.test", "r")

	bbc_new = read(io, Vector)

	bbc_new  = Array{UInt8,1}(undef, 100)

	read!(io, bbc_new)

	close(io)

N = 100_000_000;
K = 100;

# faster string sort
svec = rand("id".*string.(1:N÷K), N);
using Blosc
# save with no compression
io = open("c:/data/no_compression.bin", "w")
@time ncw = write.(Ref(io), codeunits.(svec));
ncwtot=sum(ncw)
@time ncw_compressed_tot = write(io, Blosc.compress(UInt16.(ncw)))
close(io)

@time res = column_loader(String, io, ncwtot)

io = open("c:/data/no_compression.bin", "r")
@time tmp = Vector{UInt8}(undef, sum(ncw))
@time read!(io, tmp)
close(io)

@time tmp_compressed = Blosc.compress(tmp)
io = open("c:/data/compression.bin", "w")
write(io, tmp_compressed)
close(io)

tot = sum(ncw)
ok(tmp, ncw) = begin
   	cncw = cumsum(ncw)
   	str = unsafe_string(pointer(tmp))
   	getindex.(Ref(str), Colon().(vcat(1, cncw[1:end-1].+1), cncw))
end

io = open("c:/data/bin.bin", "r")
@time tmp = Vector{UInt8}(undef, tot)
@time read!(io, tmp)
@time ok(tmp, ncw)
close(io)

# @time nc = ncodeunits.(svec)
# @time aa = Vector{UInt8}(undef, sum(nc))

# meh!(aa, svec, nc) = begin
# 	start = 1
# 	@time for (n, s) in zip(nc, svec)
# 		aa[1:n] = codeunits(s)
# 		start += n
# 	end
# end

@time meh!(aa, svec, nc)

@time a  = codeunits.(svec)
close(io)

@time ncw = write.(Ref(io), a);

using Blosc
io = open("c:/data/compression.bin", "w")
@time ncw = write.(Ref(io), Blosc.compress.(codeunits.(svec)));
close(io)

stringit(n, io) = begin
    tmp = Vector{UInt8}(undef, n)
    read!(io, tmp)
    tmp
end

tot = sum(ncw)

ok(tmp, ncw) = begin
   	cncw = cumsum(ncw)
   	str = unsafe_string(pointer(tmp))
   	getindex.(Ref(str), Colon().(vcat(1, cncw[1:end-1].+1), cncw))
end

io = open("c:/data/bin.bin", "r")
@time tmp = Vector{UInt8}(undef, tot)
@time read!(io, tmp)
@time ok(tmp, ncw)
close(io)

a = pointer(tmp[1:8])
@time str = unsafe_string(a)
close(io)

io = open("c:/data/bin.bin", "r")
@time svec2 = stringit.(ncw, Ref(io))
close(io)

@time svec2 = unsafe_string.(pointer.(aa))
all(svec2 .== svec)

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
#* Wrapper of GameInterface
# redirects AlphaZero.MCTS calls to MuNetwork
#TODO Uptade description
"""
    MuGameEnvWrapper

Wraper aroud GameEnv that redirects GI calls
- `game::GI.AbstractGameEnv`  GameEnv that is wrapped
- `dynamics_oracle`           g(sᵏ⁻¹,aᵏ) -> (rᵏ,sᵏ) that returns reward and next state
- `curstate`                  current (hidden) state of game
- `isrootstate::Bool`         whether or not curstate is the root state of MCTS tree
- `white_playing::Bool`       whether or not white is currently playing
- `lastwr::Float64`           last white reward obtained from dynamics_oracle
"""
mutable struct MuGameEnvWrapper{State, R, RRef<:Base.RefValue} <: GI.AbstractGameEnv
  game :: GI.AbstractGameEnv
  recurrent_inference :: R # f(g(s,a)) -> p, v, s⁺¹, r
  curstate :: State
  isrootstate :: Bool
  white_playing :: Bool # MuWrapper changes player each round
  lastwr :: Float32
  result_ref :: RRef
end

GI.game_terminated(game::MuGameEnvWrapper) = false
# GI.game_terminated(game::MuGameEnvWrapper) = GI.game_terminated(GI.init(GI.spec(game.game),game.curstate))

GI.current_state(game::MuGameEnvWrapper) = game.curstate

GI.available_actions(game::MuGameEnvWrapper) = game.isrootstate ? GI.available_actions(game.game) : GI.actions(GI.spec(game.game))

GI.white_playing(game::MuGameEnvWrapper) = game.white_playing

function GI.play!(game::MuGameEnvWrapper, action)
  # dynamics returns reward from the perspective of last player (rᵂ , sᴮ) <- g(sᵂ, aᵂ)
  # (game.lastwr, game.curstate) = game.dynamics_oracle((game.curstate, action))
  game.result_ref[] = game.recurrent_inference((game.curstate, action))
  (game.curstate, game.lastwr) = game.result_ref[][3:4]
  if GI.two_players(GI.spec(game.game))
    game.white_playing || (game.lastwr = -game.lastwr) # then it is changed if black was playing and saved as lastwr
    game.white_playing = !game.white_playing # player change
  end
  game.isrootstate = false
end

GI.white_reward(game::MuGameEnvWrapper) = game.lastwr #vulnerable when call white_reward() before play!()

function GI.clone(game::MuGameEnvWrapper)
  MuGameEnvWrapper(
    GI.clone(game.game),
    game.recurrent_inference,
    game.curstate,
    game.isrootstate,
    game.white_playing,
    game.lastwr,
    game.result_ref)
end

# cache oracle (nn) results, so it can be quickly obtained without re-running network
struct CachedOracle{O, KeyType, ValueType}
  oracle :: O
  lookuptable :: Dict{KeyType,ValueType}
end

# oracle, typeof(State), typeof(Action), typeof(Reward) - used with dynamics
function CachedOracle(oracle, ValType, S, A=Int, R=Float32)
  return CachedOracle(
    oracle,
    Dict{Tuple{S,A}, ValType}()) # (sᵏ⁻¹,aᵏ) => (rᵏ,sᵏ)
end

function (oracle::CachedOracle)(arg)
  if haskey(oracle.lookuptable, arg)
    answer = oracle.lookuptable[arg]
  else
    answer = oracle.oracle(arg)
    oracle.lookuptable[arg] = answer
  end
  return answer
end


text:
module BasicTest

using Dubins3D

# Convert degreees to radians
deg2rad(x) = pi * x / 180.

# Initial and final configurations [x, y, z, heading angle, pitch angle]
qi = [200., 500., 200., deg2rad(180.), deg2rad(-5.)]
qf = [500., 350., 100., deg2rad(0.), deg2rad(-5.)]
# Minimum turning radius
rhomin = 40.
# Pich angle constraints [min_pitch, max_pitch]
pitchmax = deg2rad.([-15., 20.])

maneuver = DubinsManeuver3D(qi, qf, rhomin, pitchmax)

# Length of the 3D Dubins path
@show maneuver.length

# Sample the manever by 500 samples
samples = compute_sampling(maneuver; numberOfSamples = 500)
# First and last samples - should be equal to qi and qf
@show samples[1]
@show samples[end]

end


text:
## io library
function io_print(line::Ref{UInt}, ns::Namespace, args...)::AetherlangObject
    print(join([string(arg) for arg in args], ' '))
    AetherlangObject()
end

function io_println(line::Ref{UInt}, ns::Namespace, args...)::AetherlangObject
    io_print(line, ns, args..., '\n')
    AetherlangObject()
end

function io_input(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject
    print(arg.dataref)
    AetherlangObject(readline())
end

function io_intinput(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject
    print(arg.dataref)
    AetherlangObject(parse(Int, readline()))
end

function io_floatinput(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject
    print(arg.dataref)
    AetherlangObject(parse(Float64, readline()))
end

function io_autoinput(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject
    print(arg.dataref)
    s::String = readline()
    try
        return AetherlangObject(parse(Int, s))
    catch
        try
            return AetherlangObject(parse(Float64, s))
        catch
            return AetherlangObject(s)
        end
    end
    AetherlangObject()
end

io_namespace_modify = Namespace(
    "print" => AetherlangObject(io_print),
    "println" => AetherlangObject(io_println),
    "input" => AetherlangObject(io_input),
    "intinput" => AetherlangObject(io_intinput),
    "floatinput" => AetherlangObject(io_floatinput),
    "autoinput" => AetherlangObject(io_autoinput)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
@inbounds revtuple(t::NTuple{N,Any}) where {N} = ntuple(i->t[N+1-i],Val(N))
@inbounds revinds(t::NTuple{N,Any},ns::Int) where {N} = ntuple(i->ns+1-t[N+1-i],Val(N))
@inbounds gettuple(t1::NTuple,t2::NTuple{N,Any}) where {N} = ntuple(i->t1[t2[i]],Val(N))

setindex(t::Tuple,val,::Val{i}) where {i} = _setindex((),t,val,Val(i))
@inline _setindex(o::Tuple,t::Tuple,val,::Val{i}) where {i} = _setindex((o...,first(t)),tail(t),val,Val(i))
@inline _setindex(o::NTuple{i,Any},t::Tuple,val,::Val{i}) where {i} = (front(o)...,val,t...)

ntuple_sans_m(m,::Val{n}) where {n} = sorted_setdiff(ntuple(identity,Val(n)),(m,))

# Thanks to @mbauman on Discourse for the following
# https://discourse.julialang.org/t/type-stable-difference-of-tuples/3933/4
@inline function sorted_setdiff(t1::Tuple, t2::Tuple)
    if t1[1] == t2[1]
        sorted_setdiff(tail(t1), tail(t2))
    else
        (t1[1], sorted_setdiff(tail(t1), t2)...)
    end
end
@noinline sorted_setdiff(t1::Tuple{}, t2::Tuple) = throw(ArgumentError("duplicate or missing index $(t2[1])"))
sorted_setdiff(t1::Tuple, ::Tuple{}) = t1
sorted_setdiff(::Tuple{}, ::Tuple{}) = ()


text:
export ca_roots, ca_roots_path

"""
    ca_roots() :: Union{Nothing, String}

The `ca_roots()` function tells the caller where, if anywhere, to find a file or
directory of PEM-encoded certificate authority roots. By default, on systems
like Windows and macOS where the built-in TLS engines know how to verify hosts
using the system's built-in certificate verification mechanism, this function
will return `nothing`. On classic UNIX systems (excluding macOS), root
certificates are typically stored in a file in `/etc`: the common places for the
current UNIX kernel will be searched and if one of these paths exists, it will
be returned; if none of these typical root certificate paths exist, then the
path to the set of root certificates that are bundled with Julia is returned.

The default value returned by `ca_roots()` may be overridden by setting the
`JULIA_SSL_CA_ROOTS_PATH` environment variable to a non-empty value, in which
case this function will always return that path (whether it exists or not).
"""
ca_roots()::Union{Nothing,String} = _ca_roots(true)

"""
    ca_roots_path() :: String

The `ca_roots_path()` function is similar to the `ca_roots()` function except
that it always returns a path to a file or directory of PEM-encoded certificate
authority roots. When called on a system like Windows or macOS, where system
root certificates are not stored in the file system, it will currently return
the path to the set of root certificates that are bundled with Julia. (In the
future, this function may instead extract the root certificates from the system
and save them to a file whose path would be returned.)

If it is possible to configure a library that uses TLS to use the system
certificates that is generally preferrable: i.e. it is better to use
`ca_roots()` which returns `nothing` to indicate that the system certs should be
used. The `ca_roots_path()` function should only be used when configuring
libraries which _require_ a path to a file or directory for root certificates.

The default value returned by `ca_roots_path()` may be overridden by setting the
`JULIA_SSL_CA_ROOTS_PATH` environment variable to a non-empty value, in which
case this function will always return that path (whether it exists or not).
"""
ca_roots_path()::String = _ca_roots(false)

const BUNDLED_CA_ROOTS = normpath(Sys.BINDIR, "..", "share", "julia", "cert.pem")

const LINUX_CA_ROOTS = [
    "/etc/ssl/cert.pem"								    # Alpine Linux
    "/etc/ssl/ca-bundle.pem"			                # OpenSUSE
    "/etc/ssl/certs/ca-certificates.crt"                # Debian/Ubuntu/Gentoo etc.
    "/etc/pki/ca-trust/extracted/pem/tls-ca-bundle.pem"	# CentOS/RHEL 7
    "/etc/pki/tls/certs/ca-bundle.crt"	                # Fedora/RHEL 6
    "/etc/pki/tls/cacert.pem"			                # OpenELEC
]

const BSD_CA_ROOTS = [
    "/etc/ssl/cert.pem"                                 # OpenBSD
    "/usr/local/share/certs/ca-root-nss.crt"            # FreeBSD
    "/usr/local/etc/ssl/cert.pem"                       # FreeBSD
]

const SYSTEM_CA_ROOTS_LOCK = ReentrantLock()
const SYSTEM_CA_ROOTS = Ref{String}()

function system_ca_roots()
    lock(SYSTEM_CA_ROOTS_LOCK) do
        isassigned(SYSTEM_CA_ROOTS) && return
        search_path = Sys.islinux() ? LINUX_CA_ROOTS :
            Sys.isbsd() && !Sys.isapple() ? BSD_CA_ROOTS : String[]
        for path in search_path
            ispath(path) || continue
            SYSTEM_CA_ROOTS[] = path
            return
        end
        # TODO: extract system certs on Windows & macOS
        SYSTEM_CA_ROOTS[] = BUNDLED_CA_ROOTS
    end
    return SYSTEM_CA_ROOTS[]
end

function _ca_roots(allow_nothing::Bool)
    path = get(ENV, "JULIA_SSL_CA_ROOTS_PATH", nothing)
    path !== nothing && !isempty(path) && return path
    allow_nothing && (Sys.iswindows() || Sys.isapple()) && return nothing
    return system_ca_roots()
end


text:
const SuperiorTemporalGyrus = BilateralStructure("superior temporal gyrus")
const TemporalPole = BilateralStructure("temporal pole")
const HeschlGyrus = BilateralStructure("Heschl's gyrus")
const SuperiorTemporalGyrus = BilateralStructure("superior temporal gyrus")
const MiddleTemporalGyrus = BilateralStructure("middle temporal gyrus")
const InferiorTemporalGyrus = BilateralStructure("inferior temporal gyrus")
const FusiformGyrus = BilateralStructure("fusiform gyrus")

const Hippocampus = BilateralStructure("hippocampus")
const Amygdala = BilateralStructure("amygdala")
const Uncus = BilateralStructure("uncus")
const DentateGyrus = BilateralStructure("dentate gyrus")

const ParahippocampalGyrus = BilateralStructure("parahippocampal gyrus")

const MedialOccipitotemporalGyrus = BilateralStructure("medial occipitotemporal")

const Subiculum = BilateralStructure("subiculum")

const EntorhinalCortex = BilateralStructure("entorhinal cortex")

const LateralParietoTemporalLine = BilateralStructure("lateral parietotemporal line")

const SuperiorTemporalSulcus = BilateralStructure("superior temporal sulcus")

const MesialTemporalLobe = BilateralStructure("mesial temporal lobe")

const TemporalLobe = BilateralStructure("temporal lobe")

#= TODO functional divisions
const PrimaryAuditoryCortex = BilateralStructure(:primary_auditory_cortex)
const SecondaryAuditoryCortex = BilateralStructure(:secondary_auditory_cortex)
const InferiorTemporalCortex = BilateralStructure(:inferior_temporal_cortex)
const PosteriorInferiorTemporalCortex = BilateralStructure(:posterior_inferior_temporal_cortex)

    functional = (
        PrimaryAuditoryCortex,
        SecondaryAuditoryCortex,
        InferiorTemporalCortex,
        PosteriorInferiorTemporalCortex
    )
=#


text:
using Documenter, PhysicalOptics

DocMeta.setdocmeta!(PhysicalOptics, :DocTestSetup, :(using PhysicalOptics); recursive=true)
makedocs(modules=[PhysicalOptics],
         sitename="PhysicalOptics.jl",
         doctest = false,
         pages = Any[
                "PhysicalOptics.jl" => "index.md",
         ]
        )

deploydocs(repo = "github.com/JuliaPhysics/PhysicalOptics.jl.git", devbranch="main")

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
@inbounds for i = ix:ix-1+n
    x[i] = a
end


text:
using ThreeBodyProblem
using LinearAlgebra
using Test

sys = earth_moon()
p = [sys.μ₁, sys.μ₂, sys.d]
rv = [0.5, 0, 0, 0, 0.5, 0]
rv_dim = [1e4, 0, 0, 0, 1, 0]
μ = sys.μ
t = 0

rvdot = zeros(6)
@test CR3BPdynamics!(rvdot, rv, sys, t) == nothing
@test rvdot == CR3BPdynamics(rv, sys, t)

rvdot = zeros(6)
@test CR3BPdynamics!(rvdot, rv_dim, p, t) == nothing
@test rvdot == CR3BPdynamics(rv_dim, p, t)

rvdot = zeros(6)
@test CR3BPdynamics!(rvdot, rv_dim, sys.μ, t) == nothing
@test rvdot == CR3BPdynamics(rv_dim, sys.μ, t)

rvdot = zeros(6)
@test CR3BPinert!(rvdot, rv, sys, t) == nothing
@test rvdot == CR3BPinert(rv, sys, t)

rvdot = zeros(6)
@test CR3BPinert!(rvdot, rv_dim, p, t) == nothing
@test rvdot == CR3BPinert(rv_dim, p, t)

wdot = zeros(42)
w = [reshape(I(6),36,1); rv]
@test CR3BPstm!(wdot, w, sys, t) == nothing
@test wdot == CR3BPstm(w, sys, t)

rvdot = zeros(6)
@test R2BPdynamics!(rvdot, rv_dim, sys.prim, t) == nothing
@test rvdot == R2BPdynamics(rv_dim, sys.prim, t)

rvdot = zeros(6)
@test R2BPdynamics!(rvdot, rv_dim, sys.μ₁, t) == nothing
@test rvdot == R2BPdynamics(rv_dim, sys.μ₁, t)

rvdot = zeros(6)
n = 1e-3
@test CWdynamics!(rvdot, rv, n, t) == nothing
@test rvdot == CWdynamics(rv, n, t)

bisys = earth_moon_sun()
rvdot = zeros(6)
@test BCPdynamics!(rvdot, rv, bisys, t) == nothing
@test rvdot == BCPdynamics(rv, bisys, t)