text:
"""
    iscolumn(df::AbstractDataFrame, col::Symbol)

Return true if `col` is in prepertynames of `df`
"""

function iscolumn(df::AbstractDataFrame, col::Symbol)
    col in propertynames(df)
end

"""
    check_changes(v::AbstractVector)

Given a vector it return a boolean vector of the same length checking if the i value is equal to the i-1 value. By default the first value is false
"""

function check_changes(v::AbstractVector)
    r = v .!= lag(v, default = :placeholder)
    r[1] = false
    return r
end

"""
    count_same(count::T, change) where {T <: Number}
    count_same(v:AbstractVector)

if `change` is false it adds 1 to `count` otherwise it returns exactly 1.
On a vector it is used in conjunction with accumulate to count consecutive same value on the vector `v`
and reset to 1 when a change occurs between `v[i]` and `v[i+1]`
"""

function count_same(count::T, change) where {T <: Number}
    if !change
        count+1
    else
        1
    end
end

function count_same(changes::BitArray)
    res = accumulate(count_same,changes;init=0)
    return res
end

"""
    count_different(count::T, same) where {T <: Number}
    count_same(v:AbstractVector)

if `change` is true it adds 1 to `count` otherwise it returns exactly 1.
On a vector it is used in conjunction with accumulate to count progress a count when `v[i]` is different than `v[i-1]`
"""

function count_different(count::T, change) where {T <: Number}
    if change
        count+1
    else
        count
    end
end

function count_different(changes::BitArray)
    res = accumulate(count_different,changes;init=1)
    return res
end

"""
    `get_hierarchy(v)`
Return a vector indicating how many consecutive falses or trues
preceded an event
"""

function nextcount(count::T, rewarded) where {T <: Number}
    if rewarded
        count > 0 ? count + 1 : one(T)
    else
        count < 0 ? count - 1 : -one(T)
    end
end
signedcount(v::AbstractArray{Bool}) = accumulate(nextcount, v;init = 0)
get_hierarchy(v) = lag(signedcount(v), default = 0)

"""
    gen(mouse::String; dir = joinpath(dirname(@__DIR__), "genotypes"))

Given an animal name `mouse` look in the genotypes folder for a match.
Returns a string indicating the genotype if a match is found otherwise `missing`
"""
function gen(mouse::AbstractString; dir = joinpath(dirname(@__DIR__), "genotypes"))
    genotype = "missing"
    for file in readdir(dir)
        if endswith(file, ".csv")
            df = CSV.read(joinpath(dir, file), DataFrame)
            n = names(df)[1]
            if mouse in df[!,n]
                genotype = string(n)
            end
        end
    end
    genotype
end

"""
    set_results_dir(path::AbstractString; exp_name = nothing)

Given a `path` create a directory where to store results of raw data processing.
The  destination folder is named Results followed by the current date inside `dir`.
If `path` is a file the Results folder will be created in its parent directory.
If an experiment name is provided the Results folder will be created in a folder named`exp_name` located in the parent folder 1 level above.
"""
function set_results_dir(path::AbstractString; exp_name = nothing)
    if isfile(path)
        dir = dirname(path)
    elseif isdir(path)
        dir = path
    else
        error("couldn't resolve path destination")
    end
    results_dir = "Results_" * string(today())
    if exp_name != nothing
        parent_dir = dirname(dir)
        destination_dir = joinpath(parent_dir,exp_name,results_dir)
    else
        destination_dir = joinpath(dir, results_dir)
    end
    if !isdirpath(destination_dir)
        mkpath(destination_dir)
    end
    return destination_dir
end

"""
    exp_calendar(df::AbstractDataFrame)

Given a pokes, bouts or streaks DataFrame `df` uses Day variable to count consecutive days of experiment for each animals.
Add a new column storing the values
"""

function exp_calendar!(df::AbstractDataFrame)
    df[!,:ExpSession] = Vector{Int64}(undef,nrow(df))
    combine(groupby(df,:MouseID)) do dd
        ExpCalendar = Dict(d => n for (n,d) in enumerate(sort(union(dd.Session))))
        dd[:,:ExpSession] = [get(ExpCalendar,x,Date(2000,12,31)) for x in dd.Session]
    end
    return df
end

"""
    conditional_calendar!(df::AbstractDataFrame,condition::Symbol)

Given a pokes, bouts or streaks DataFrame `df` uses ExpDay variable to count consecutive days of the same value specified by `condition` for each animals
Add a new column storing the values
"""

function conditional_calendar!(df::AbstractDataFrame,condition::Symbol)
    gd = groupby(df,[:MouseID, :Session])
    # synthesis = combine([condition,:Session] => (p,s) ->
    #         (
    #             Flexi = length(union(p)) > 1,
    #             Condition = join(union(p),"*"),
    #             Session = s[1]
    #         )
    #     ,gd)
    synthesis = combine(groupby(df,[:MouseID, :Session])) do dd
        (Flexi = length(union(dd[:,condition])) > 1,
        Condition = join(union(dd[:,condition]),"*"),
        Session = dd[1, :Session])
    end
    rename!(synthesis,:Condition => condition)
    conditionSession = Symbol(string(condition)*"Session")

    synthesis[:,conditionSession] = Vector{Int64}(undef,nrow(synthesis))
    combine(groupby(synthesis,:MouseID)) do dd
        dd[:,conditionSession] = count_same(check_changes(dd[:,condition]))
    end

    condition_dict = Dict(session => day for (session,day) in zip(synthesis[:,:Session],synthesis[:,conditionSession]))
    df[!,conditionSession] = [get(condition_dict,x,Date(2000,12,31)) for x in df.Session]
    return df
end

"""
    add_calendars!(df::AbstractDataFrame)

Given a pokes, bouts or streaks DataFrame `df` adds variables abut Experiment Day and consecutive Protocol Day for each animals
"""

function add_calendars!(df::AbstractDataFrame)
    exp_calendar!(df)
    conditional_calendar!(df,:Protocol)
end


text:
"""
    abstract type AbstractKernelMatrix{T} <: AbstractMatrix{T}

Interface for abstract matrices represented through a kernel function `f`, target
elements `X`, and source elements `Y`. The matrix entry
`i,j` is given by `f(X[i],Y[j])`. Concrete subtypes should implement at least

    `Base.getindex(K::AbstractKernelMatrix,i::Int,j::Int)`

If a more efficient implementation of `getindex(K,I::UnitRange,I::UnitRange)`,
`getindex(K,I::UnitRange,j::Int)` and `getindex(adjoint(K),I::UnitRange,j::Int)`
is available (e.g. with SIMD vectorization), implementing such methods can
improve the speed of assembling an [`HMatrix`](@ref).
"""
abstract type AbstractKernelMatrix{T} <: AbstractMatrix{T} end

"""
    KernelMatrix{Tf,Tx,Ty,T} <:: AbstractKernelMatrix{T}

Generic kernel matrix representing a kernel function acting on two sets of
elements. If `K` is a `KernelMatrix`, then `K[i,j] = f(X[i],Y[j])` where
`f::Tf=kernel(K)`, `X::Tx=rowelements(K)` and `Y::Ty=colelements(K)`.

# Examples
```julia
X = rand(Geometry.Point2D,100)
Y = rand(Geometry.Point2D,100)
K = KernelMatrix(X,Y) do x,y
    sum(x+y)
end
```
"""
struct KernelMatrix{Tf,Tx,Ty,T} <: AbstractKernelMatrix{T}
    f::Tf
    X::Tx
    Y::Ty
end

Base.size(K::KernelMatrix)                   = length(K.X), length(K.Y)
Base.getindex(K::KernelMatrix,i::Int,j::Int) =  K.f(K.X[i],K.Y[j])

rowelements(K::KernelMatrix) = K.X
colelements(K::KernelMatrix) = K.Y
kernel(K::KernelMatrix) = K.f

function KernelMatrix(f,X,Y)
    T = Base.promote_op(f,eltype(X),eltype(Y))
    KernelMatrix{typeof(f),typeof(X),typeof(Y),T}(f,X,Y)
end

ext: jl

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
#### Create states

function create_state(
    balance_law,
    grid,
    st::AbstractStateType;
    fill_nan = false,
)
    topology = grid.topology
    Np = dofs_per_element(grid)

    h_vgeo = Array(grid.vgeo)
    FT = eltype(h_vgeo)
    DA = arraytype(grid)

    weights = view(h_vgeo, :, grid.Mid, :)
    weights = reshape(weights, size(weights, 1), 1, size(weights, 2))

    # TODO: Clean up this MPIStateArray interface...
    ns = number_states(balance_law, st)
    V = vars_state(balance_law, st, FT)
    if st isa GradientLaplacian
        ns = 3ns
        # Since the names do not match the storage we do not set them
        V = @vars()
    end
    state = MPIStateArray{FT, V}(
        topology.mpicomm,
        DA,
        Np,
        ns,
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
        weights = weights,
    )
    fill_nan && fill!(state, NaN)
    return state
end

function init_state(state, balance_law, grid, direction, ::Auxiliary)
    init_state_auxiliary!(balance_law, state, grid, direction)
    return state
end


text:
module ArrayInterfaceTracker

using ArrayInterfaceCore
using Tracker

ArrayInterfaceCore.ismutable(::Type{<:Tracker.TrackedArray}) = false
ArrayInterfaceCore.ismutable(T::Type{<:Tracker.TrackedReal}) = false
ArrayInterfaceCore.can_setindex(::Type{<:Tracker.TrackedArray}) = false
ArrayInterfaceCore.fast_scalar_indexing(::Type{<:Tracker.TrackedArray}) = false
ArrayInterfaceCore.aos_to_soa(x::AbstractArray{<:Tracker.TrackedReal,N}) where {N} = Tracker.collect(x)

end # module


text:
include("spaces.jl")
include("spec.jl")

mutable struct GymEnv
    name
    spec
    action_space
    observation_space
    reward_range
    gymenv
end

function GymEnv(id::String)
    gymenv = nothing
    try
        gymenv = gym[:make](id)
    catch e
        error("Error received during the initialization of $id\n$e")
    end

    spec = Spec(gymenv[:spec][:id],
                gymenv[:spec][:trials],
                gymenv[:spec][:reward_threshold],
                gymenv[:spec][:nondeterministic],
                gymenv[:spec][:tags],
                gymenv[:spec][:max_episode_steps],
                gymenv[:spec][:timestep_limit]
               )
    action_space = julia_space(gymenv[:action_space])
    observation_space = julia_space(gymenv[:observation_space])

    env = GymEnv(id, spec, action_space,
                 observation_space, gymenv[:reward_range], gymenv)
    return env
end

reset!(env::GymEnv) = env.gymenv[:reset]()

function render(env::GymEnv; mode="human", close=false)
    env.gymenv[:render](; mode=mode, close=close)
end

function step!(env::GymEnv, action)
    ob, reward, done, information = env.gymenv[:step](action)
    return ob, reward, done, information
end

close!(env::GymEnv) = env.gymenv[:close]()

seed!(env::GymEnv, seed=nothing) = env.gymenv[:seed](seed)


text:
module Procrustes

export Coords, procAlign

#define the fundimental type
type Coords
    rawCoords::Array
    alignCoords::Array
    centroidSize::Vector
end

#main function
"""
Implements a [generalized procrustes analysis](https://en.wikipedia.org/wiki/Procrustes_analysis#Generalized_Procrustes_analysis_.28GPA.29).

# Arguments
'data::Array': A 3D array with landmarks on the first dimension
                x,y,z values on the 2nd dimension
                samples on the 3rd dimension

'scale::Bool': Scale samples to unit size

'tol::Float': Tolerance for convergence of the iterative alignment
"""
function procAlign(data::Array,scale=true,tol=10e-5)

    centSize=centriodSize(data)
    rawCoords=copy(data)
    procAlign!(data,scale,tol)
    return(Coords(rawCoords,data,centSize));

end

function procAlign!(x::Array, scale, tol)

    centerShape!(x)

    if scale
        scaleShape!(x)
    end

    #pick random sample for frist iterative alignment
    ind=rand(1:(size(x)[3]));
    meanShapeIn=x[:,:,ind];
    meanShapeOut=zeros(size(meanShapeIn));
    procError=100;

    #iterate until meanshape converges
    while procError > tol

        for sample in 1:size(x)[3]
            x[:,:,sample]=align_two_shapes(meanShapeIn[:,:,1],x[:,:,sample])
        end

        meanShapeOut=meanShape(x)
        procError=procDist(meanShapeIn,meanShapeOut)

        meanShapeIn=meanShapeOut
    end

end

function align_two_shapes(refShape::Array, movingShape::Array)
    H=transpose(refShape)*movingShape

    U,S,V=svd(H)
    R=V*transpose(U)

    return(movingShape*R)
end

function centerShape!(x::Array)
    #for each sample
    for k in 1:size(x)[3];
        for j in 1:size(x)[2];

            colMean=mean(x[:,j,k])
            for i in 1:size(x)[1]
                x[i,j,k]=x[i,j,k]-colMean
            end
        end
    end
end

function scaleShape!(x::Array)
    for k in 1:size(x)[3]
        x[:,:,k]=x[:,:,k]/sqrt(sum(var(x[:,:,k], 1)))
    end
end

function meanShape(x::Array)
    return(mean(x,3))
end

function procDist(x::Array, y::Array)
    return(sqrt(sum((x-y).^2)))
end

"""
Calculates the centroid size of landmark configurations.

# Arguments
'data::Array': A 3D array with landmarks on the first dimension
                x,y,z values on the 2nd dimension
                samples on the 3rd dimension
"""
function centriodSize(x::Array)

    centSize=zeros(size(x)[3])

    #for each sample
    for k in 1:size(x)[3];

        colMean=mean(x[:,:,k],1)
        for i in 1:size(x)[1]
            centSize[k]+=(sum((x[i,:,k]-colMean).^2))
        end
    end

    map!(sqrt,centSize)
    return(centSize)
end

end # module


text:
v1 = SVector{3,Float64}(0.0, 0.0, 0.0)
v2 = SVector{3,Float64}(1.0, 0.0, 0.0)
v3 = SVector{3,Float64}(0.0, 1.0, 0.0)
v4 = SVector{3,Float64}(0.0, 0.0, 1.0)

@testset "triangle" begin
    tri_Sep = (v1, v2, v3)
    tri_SV  = (SVector{3,SVector{3,Float64}}(tri_Sep...),)
    tri_Tup  = ((tri_Sep...),)
    for tri_rep = (tri_Sep, tri_SV, tri_Tup)
        @test area(tri_rep...) ≈ 0.5
        @test centroid(tri_rep...) ≈ SVector{3,Float64}(1/3,1/3,0)
        @test triangleNormal(tri_rep...) ≈ SVector{3,Float64}(0,0,1)
    end
    @test triangle_area(tri_Sep, triangleNormal(tri_Sep...)) ≈ 0.5
end

@testset "tetrahedron" begin
    tet_Tup = (v1, v2, v3, v4)
    tet_SV  = (SVector{4,SVector{3,Float64}}(tet_Tup...),)
    for tet_rep = (tet_SV, tet_Tup)
        @test volume(tet_rep...) ≈ 1/6
        @test centroid(tet_rep...) ≈ SVector{3,Float64}(1,1,1)/4
    end
end


text:
import JSON #2190G
const f_en = open(joinpath(@__DIR__(), "entries.json"))
const arr_en = JSON.parse(f_en)
const words_en = [String(s[1]) for s in arr_en]
const regexes_en = [Regex(s[2]) for s in arr_en]
const subst_en = [Base.SubstitutionString(replace(s[3], "\$" => "\\")) for s in arr_en]

"""
    wikiscan(text, start = 1, regexes = regexes_en, words = words_en, subst = subst_en)

Scan for misspellings using wikipedia regexes.

# Arguments:
- `text`: The text to look in
- `start`: Line number to start at
- `regexes`: list of regexes to use
- `words`: Dictionary
- `subst`: List of substitution strings

# Credit: Jacob Wikmark, github.com/lancebeet
"""
function wikiscan(text, start = 1, regexes = regexes_en, words = words_en, subst = subst_en)
    lines = split(text, '\n')
    for (i,ln) in enumerate(lines[start:end])
        i % 100 == 0 && print("\1 ", i)
        for ind = 1:length(words)
            ma = eachmatch(regexes[ind], ln)
            out1 = "\nLine $i Regex $ind\n$ln\nMatch: $ma\n"
            out2 = ""
            for m in ma
                m = m.match
                sub = replace(m, regexes[ind] => subst[ind])
                if m != sub
                    out2 *= "$m : $sub\n"
                end
            end
            out2 == "" && continue
            print(out1)
            print(out2)
        end
    end
end


text:
module BSP

importall ImmutableArrays
import Base.search

include("bsp/Mesh.jl")
include("bsp/Tree.jl")

export Bsp
type Bsp
	tree::Tree.Node
	entities::Vector{Dict{String,String}}
	vertices::Vector{Vector3{Float32}}
	max_lights::Int
end

include("bsp/File.jl")

search(bsp::Bsp, pos::Vector3) = search(bsp.tree, pos)

end


text:
using BinDeps

@BinDeps.setup

sptk = library_dependency("libSPTK", aliases=["libSPTK", "SPTK-3"])

# TODO: needs update for windows
const version = is_windows() ? "3.8.9" : "3.10.2"

github_root = "https://github.com/r9y9/SPTK"
arch = Sys.WORD_SIZE == 64 ? "x86_64" : "i686"
major = version[1]
provides(Binaries,
         URI("$(github_root)/releases/download/v$(version)/sptk-$(major)_mingw$(Sys.WORD_SIZE)_$(arch).zip"),
         sptk, unpacked_dir = "usr/lib", os = :Windows)

provides(Sources,
         URI("$(github_root)/archive/v$(version).tar.gz"),
         sptk,
         unpacked_dir="SPTK-$(version)")

prefix = joinpath(BinDeps.depsdir(sptk), "usr")
srcdir = joinpath(BinDeps.depsdir(sptk), "src", "SPTK-$(version)")

provides(SimpleBuild,
          (@build_steps begin
              GetSources(sptk)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `./waf configure --prefix=$prefix`
                  `./waf build`
                  `./waf install`
              end
           end), sptk, os = :Unix)

@BinDeps.install Dict(:libSPTK => :libSPTK)


text:
using PyCall
using Conda

Conda.add("seaborn")
Conda.add("scikit-learn")
Conda.add("qt")
Conda.add("graphviz")
Conda.add("pydot")
Conda.add("platypus-opt", channel="conda-forge")

Conda.pip_interop(true)
Conda.pip("install", ["SALib", "mpldatacursor", "git+https://github.com/Project-Platypus/PRIM.git#egg=prim", "git+https://github.com/Project-Platypus/Rhodium.git#egg=rhodium"])


text:
export
    nutation_constants,
    nut_in_lon,
    nut_in_obl,
    nut_in_ra,
    obliquity_high,
    obliquity,
    true_obliquity

#    D, M, M1, F, omega, psiK, psiT, epsK, epsT
nutation_table = Array[
    [ 0,  0,  0,  0,  1, -171996, -1742, 92025,  89],
    [-2,  0,  0,  2,  2,  -13187,   -16,  5736, -31],
    [ 0,  0,  0,  2,  2,   -2274,    -2,   977,  -5],
    [ 0,  0,  0,  0,  2,    2062,     2,  -895,   5],
    [ 0,  1,  0,  0,  0,    1426,   -34,    54,  -1],
    [ 0,  0,  1,  0,  0,     712,     1,    -7,   0],
    [-2,  1,  0,  2,  2,    -517,    12,   224,  -6],
    [ 0,  0,  0,  2,  1,    -386,    -4,   200,   0],
    [ 0,  0,  1,  2,  2,    -301,     0,   129,  -1],
    [-2, -1,  0,  2,  2,     217,    -5,   -95,   3],
    [-2,  0,  1,  0,  0,    -158,     0,     0,   0],
    [-2,  0,  0,  2,  1,     129,     1,   -70,   0],
    [ 0,  0, -1,  2,  2,     123,     0,   -53,   0],
    [ 2,  0,  0,  0,  0,      63,     0,     0,   0],
    [ 0,  0,  1,  0,  1,      63,     1,   -33,   0],
    [ 2,  0, -1,  2,  2,     -59,     0,    26,   0],
    [ 0,  0, -1,  0,  1,     -58,    -1,    32,   0],
    [ 0,  0,  1,  2,  1,     -51,     0,    27,   0],
    [-2,  0,  2,  0,  0,      48,     0,     0,   0],
    [ 0,  0, -2,  2,  1,      46,     0,   -24,   0],
    [ 2,  0,  0,  2,  2,     -38,     0,    16,   0],
    [ 0,  0,  2,  2,  2,     -31,     0,    13,   0],
    [ 0,  0,  2,  0,  0,      29,     0,     0,   0],
    [-2,  0,  1,  2,  2,      29,     0,   -12,   0],
    [ 0,  0,  0,  2,  0,      26,     0,     0,   0],
    [-2,  0,  0,  2,  0,     -22,     0,     0,   0],
    [ 0,  0, -1,  2,  1,      21,     0,   -10,   0],
    [ 0,  2,  0,  0,  0,      17,    -1,     0,   0],
    [ 2,  0, -1,  0,  1,      16,     0,    -8,   0],
    [-2,  2,  0,  2,  2,     -16,     1,     7,   0],
    [ 0,  1,  0,  0,  1,     -15,     0,     9,   0],
    [-2,  0,  1,  0,  1,     -13,     0,     7,   0],
    [ 0, -1,  0,  0,  1,     -12,     0,     6,   0],
    [ 0,  0,  2, -2,  0,      11,     0,     0,   0],
    [ 2,  0, -1,  2,  1,     -10,     0,     5,   0],
    [ 2,  0,  1,  2,  2,      -8,     0,     3,   0],
    [ 0,  1,  0,  2,  2,       7,     0,    -3,   0],
    [-2,  1,  1,  0,  0,      -7,     0,     0,   0],
    [ 0, -1,  0,  2,  2,      -7,     0,     3,   0],
    [ 2,  0,  0,  2,  1,      -7,     0,     3,   0],
    [ 2,  0,  1,  0,  0,       6,     0,     0,   0],
    [-2,  0,  2,  2,  2,       6,     0,    -3,   0],
    [-2,  0,  1,  2,  1,       6,     0,    -3,   0],
    [ 2,  0, -2,  0,  1,      -6,     0,     3,   0],
    [ 2,  0,  0,  0,  1,      -6,     0,     3,   0],
    [ 0, -1,  1,  0,  0,       5,     0,     0,   0],
    [-2, -1,  0,  2,  1,      -5,     0,     3,   0],
    [-2,  0,  0,  0,  1,      -5,     0,     3,   0],
    [ 0,  0,  2,  2,  1,      -5,     0,     3,   0],
    [-2,  0,  2,  0,  1,       4,     0,     0,   0],
    [-2,  1,  0,  2,  1,       4,     0,     0,   0],
    [ 0,  0,  1, -2,  0,       4,     0,     0,   0],
    [-1,  0,  1,  0,  0,      -4,     0,     0,   0],
    [-2,  1,  0,  0,  0,      -4,     0,     0,   0],
    [ 1,  0,  0,  0,  0,      -4,     0,     0,   0],
    [ 0,  0,  1,  2,  0,       3,     0,     0,   0],
    [ 0,  0, -2,  2,  2,      -3,     0,     0,   0],
    [-1, -1,  1,  0,  0,      -3,     0,     0,   0],
    [ 0,  1,  1,  0,  0,      -3,     0,     0,   0],
    [ 0, -1,  1,  2,  2,      -3,     0,     0,   0],
    [ 2, -1, -1,  2,  2,      -3,     0,     0,   0],
    [ 0,  0,  3,  2,  2,      -3,     0,     0,   0],
    [ 2, -1,  0,  2,  2,      -3,     0,     0,   0]]

# Constant terms.

_kD  = [deg2rad(297.85036), deg2rad(445267.111480), deg2rad(-0.0019142), deg2rad( 1.0/189474)]
_kM  = [deg2rad(357.52772), deg2rad( 35999.050340), deg2rad(-0.0001603), deg2rad(-1.0/300000)]
_kM1 = [deg2rad(134.96298), deg2rad(477198.867398), deg2rad( 0.0086972), deg2rad( 1.0/ 56250)]
_kF  = [deg2rad( 93.27191), deg2rad(483202.017538), deg2rad(-0.0036825), deg2rad( 1.0/327270)]
_ko  = [deg2rad(125.04452), deg2rad( -1934.136261), deg2rad( 0.0020708), deg2rad( 1.0/450000)]

# Return some values needed for both nut_in_lon() and nut_in_obl()

"""
    nutation_constants(T)

Functions to calculate nutation and obliquity values.

[Meeus-1998: table 22.A]
"""
function nutation_constants(T)
    D     = mod2pi(polynomial(_kD,  T))
    M     = mod2pi(polynomial(_kM,  T))
    M1    = mod2pi(polynomial(_kM1, T))
    F     = mod2pi(polynomial(_kF,  T))
    omega = mod2pi(polynomial(_ko,  T))
    return (D, M, M1, F, omega)
end

"""
    nut_in_lon(jd)

Return the nutation in longitude.

High precision. [Meeus-1998: pg 144]

Parameters:
    jd : Julian Day in dynamical time

Returns:
    nutation in longitude, in radians
"""
function nut_in_lon(jd)
    # TODO nut_in_lon() factor the /1e5 and /1e6 adjustments into the table.
    T = jd_to_jcent(jd)
    D, M, M1, F, omega = nutation_constants(T)
    deltaPsi = 0.0
    for v in nutation_table
        (tD, tM, tM1, tF, tomega, tpsiK, tpsiT, tepsK, tepsT) = v
        arg = D*tD + M*tM + M1*tM1 + F*tF + omega*tomega
        deltaPsi = deltaPsi + (tpsiK/10000.0 + tpsiT/100000.0 * T) * sin(arg)
    end
    return deg2rad(deltaPsi/3600)
end

"""
    nut_in_obl(jd)

Return the nutation in obliquity.

High precision. [Meeus-1998: pg 144]

Parameters:
    jd : Julian Day in dynamical time

Returns:
    nutation in obliquity, in radians
"""
function nut_in_obl(jd)
    #
    # Future optimization: factor the /1e5 and /1e6 adjustments into the table.
    #
    # Could turn the loop into a generator expression. Too messy?
    #
    T = jd_to_jcent(jd)
    D, M, M1, F, omega = nutation_constants(T)
    deltaEps = 0.0;
    for v in nutation_table
    		tD, tM, tM1, tF, tomega, tpsiK, tpsiT, tepsK, tepsT = v
      arg = D*tD + M*tM + M1*tM1 + F*tF + omega*tomega
      deltaEps = deltaEps + (tepsK/10000.0 + tepsT/100000.0 * T) * cos(arg)
    end
    return deg2rad(deltaEps/3600)
end

# Constant terms

_el0 = Array([deg2rad(dms_to_d(23, 26,  21.448   )),
              deg2rad(dms_to_d( 0,  0, -46.8150  )),
              deg2rad(dms_to_d( 0,  0,  -0.00059 )),
              deg2rad(dms_to_d( 0,  0,   0.001813))
             ])

"""
    obliquity(jd)

Return the mean obliquity of the ecliptic.

Low precision, but good enough for most uses. [Meeus-1998: equation 22.2].

Accuracy is 1" over 2000 years and 10" over 4000 years.

Parameters:
    jd : Julian Day in dynamical time

Returns:
    obliquity, in radians
"""
function obliquity(jd)
    T = jd_to_jcent(jd)
    return polynomial(_el0, T)
end

currentobliquity = obliquity(Dates.datetime2julian(now()))

# Constant terms

_el1 = Array([deg2rad(dms_to_d(23, 26,    21.448)),
              deg2rad(dms_to_d( 0,  0, -4680.93)),
              deg2rad(dms_to_d( 0,  0,    -1.55)),
              deg2rad(dms_to_d( 0,  0,  1999.25)),
              deg2rad(dms_to_d( 0,  0,   -51.38)),
              deg2rad(dms_to_d( 0,  0,  -249.67)),
              deg2rad(dms_to_d( 0,  0,   -39.05)),
              deg2rad(dms_to_d( 0,  0,     7.12)),
              deg2rad(dms_to_d( 0,  0,    27.87)),
              deg2rad(dms_to_d( 0,  0,     5.79)),
              deg2rad(dms_to_d( 0,  0,     2.45))])

"""
    obliquity_high(jd)

Return the mean obliquity of the ecliptic.

High precision [Meeus-1998: equation 22.3].

Accuracy is 0.01" between 1000 and 3000, and "a few arc-seconds
after 10,000 years".

Parameters:
    jd : Julian Day in dynamical time

Returns:
    obliquity, in radians
"""
function obliquity_high(jd)
    U = jd_to_jcent(jd) / 100
    return polynomial(_el1, U)
end

function true_obliquity(jd)
    return obliquity_high(jd) + nut_in_obl(jd)
end

"""
    nut_in_ra(jd)

Return the nutation in right ascension (also called equation of the equinoxes.)

Meeus-1998: page 88.

Parameters:
    jd : Julian Day in dynamical time

Returns:
    nutation, in radians
"""
function nut_in_ra(jd)
    global days_per_second
    deltapsi = rad2deg(nut_in_lon(jd)) * 3600     # deltapsi in seconds
    epsilon  = true_obliquity(jd)                 # Epsilon kept in radians
    c = deltapsi * cos(epsilon)/15                # result in seconds...
    return (c * (pi * 2) * days_per_second)       # converted radians
end


text:
using Polyester

function oneHot(input::Array{Int64}, depth::Int64, dict::Dict)
    output_data = zeros(Float32, depth, size(input, 1))
    @batch for i in 1:length(input)
        output_data[dict[input[i]],i] = 1.0f0
    end
    return output_data
end


text:
# ---
# title: Legend tweaking
# cover: assets/legend_tweaking.png
# description: Setting legend attributes
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# To tweak the position and appearance of the legend, simply use the `legend` keyword when plotting. For example

df = (x=rand(100), y=rand(100), group=rand(["a looooooong label", "an even loooooonger label", "and one more long label"], 100))
layers = linear() + mapping(color=:group)
plt = data(df) * layers * mapping(:x, :y)
draw(plt)

#

fg = draw(plt, legend=(position=:top, titleposition=:left, framevisible=true, padding=5))

# save cover image #src
mkpath("assets") #src
save("assets/legend_tweaking.png", fg) #src


text:
using Test
using SDPT3

@testset "Semidefinite Programming example" begin
    @testset "Min offdiag" begin
        blk = ["s" 2.0]
        c = [0.0 1.0
             1.0 0.0]
        A = [1.0 0.0 0.0
             0.0 0.0 1.0]
        b = [1.0
             1.0]
        obj, X, y, Z, info, runhist = sdpt3(blk, [Matrix(A')], [c], b)
        tol = 1e-5
        @test obj[1] ≈ -2 atol=tol rtol=tol
        @test obj[2] ≈ -2 atol=tol rtol=tol
        @test X[1] ≈ [ 1.0 -1.0
                      -1.0  1.0] atol=tol rtol=tol
        @test y ≈ -ones(2) atol=tol rtol=tol
        @test Z[1] ≈ ones(2, 2) atol=tol rtol=tol
        @test info["termcode"] == 0.0
    end
    @testset "Min trace" begin
        blk = ["s" 2.0]
        c = [1.0 0.0
             0.0 1.0]
        A = [0.0 1/√2 0.0]
        b = [1.0]
        obj, X, y, Z, info, runhist = sdpt3(blk, [Matrix(A')], [c], b)
        tol = 1e-5
        @test obj[1] ≈ 2 atol=tol rtol=tol
        @test obj[2] ≈ 2 atol=tol rtol=tol
        @test X[1] ≈ ones(2, 2) atol=tol rtol=tol
        @test y ≈ [2.0] atol=tol rtol=tol
        @test Z[1] ≈ [ 1.0 -1.0
                      -1.0  1.0] atol=tol rtol=tol
        @test info["termcode"] == 0.0
    end
end


text:
function grey(decisionMat::DataFrame, weights::Array{Float64,1}, fs::Array{Function,1}, zeta::Float64=0.5)::GreyResult

    mat = convert(Matrix, decisionMat)

    nrows, ncols = size(mat)

    w = unitize(weights)

    referenceRow = apply_columns(fs, mat)

    normalizedReferenceRow = zeros(Float64, ncols)

    normalizedMat = similar(mat)
    for col in 1:ncols
        mmax = maximum(mat[:, col])
        mmin = minimum(mat[:, col])
        mrange = mmax - mmin
        for row in 1:nrows
            if fs[col] == maximum
                normalizedMat[row, col] = (mat[row, col] - mmin) / mrange
            elseif fs[col] == minimum
                normalizedMat[row, col] = (mmax - mat[row, col]) / mrange
            else
                @error fs[col]
                error("Function not defined")
            end
        end
        if fs[col] == maximum
            normalizedReferenceRow[col] = (referenceRow[col] - mmin) / mrange
        elseif fs[col] == minimum
            normalizedReferenceRow[col] = (mmax - referenceRow[col]) / mrange
        end
    end

    absoluteValueMat = similar(mat)
    for row in 1:nrows
        absoluteValueMat[row,:] = normalizedReferenceRow .- normalizedMat[row, :]  .|> abs
    end

    deltamin = absoluteValueMat |> minimum
    deltamax = absoluteValueMat |> maximum
    zetadeltamax = zeta * deltamax

    greytable = similar(absoluteValueMat)
    for row in 1:nrows
        for col in 1:ncols
            greytable[row, col] = (deltamin + zetadeltamax) / (absoluteValueMat[row, col] + zetadeltamax)
        end
    end

    weightedsums = zeros(Float64, nrows)
    for i in 1:nrows
        weightedsums[i] = w .* greytable[i, :] |> sum
    end

    orderings = sortperm(weightedsums)
    bestIndex = orderings |> last

    result = GreyResult(
        referenceRow,
        normalizedMat,
        absoluteValueMat,
        greytable,
        weightedsums,
        orderings,
        bestIndex
    )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
module PhysicalFESpacesTests

using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using Gridap.Fields
using Gridap.FESpaces
using Gridap.Polynomials
using Gridap.CellData
using Gridap.TensorValues
using Test

u(x) = x[1]+2*x[2] + 3

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
order = 1

Ω = Triangulation(model)
dΩ = Measure(Ω,order)

cell_fe = FiniteElements(ReferenceDomain(),model,lagrangian,Float64,order)
V = FESpace(model,cell_fe)
@test num_free_dofs(V) == 9
cell_fe = FiniteElements(ReferenceDomain(),model,lagrangian,Float64,order,conformity=:L2)
V = FESpace(model,cell_fe)
@test num_free_dofs(V) == 16

@test DomainStyle(get_fe_basis(V)) == ReferenceDomain()
@test DomainStyle(get_fe_dof_basis(V)) == ReferenceDomain()

uh = interpolate(u,V)
e = u - uh
@test sqrt(sum(∫(e*e)*dΩ)) < 10e-8

cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)
V = FESpace(model,cell_fe)
@test num_free_dofs(V) == 9

cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order,conformity=:L2)
V = FESpace(model,cell_fe)
@test num_free_dofs(V) == 16

@test DomainStyle(get_fe_basis(V)) == PhysicalDomain()
@test DomainStyle(get_fe_dof_basis(V)) == PhysicalDomain()

uh = interpolate(u,V)
e = u - uh
@test sqrt(sum(∫(e*e)*dΩ)) < 10e-8

end #module


text:
@testset "indentation" begin
    mds = """
        A
            B
            C
        D"""

    tokens = F.find_tokens(mds, F.MD_TOKENS, F.MD_1C_TOKENS)
    F.find_indented_blocks!(tokens, mds)
    toks = deepcopy(tokens)

    @test tokens[1].name == :LR_INDENT
    @test tokens[2].name == :LR_INDENT
    @test tokens[3].name == :LINE_RETURN

    ocp = F.OCProto(:CODE_BLOCK_IND, :LR_INDENT, (:LINE_RETURN,), false)

    blocks, tokens = F.find_ocblocks(tokens, ocp)

    @test blocks[1].name == :CODE_BLOCK_IND
    @test F.content(blocks[1]) == "B\n    C"

    blocks, tokens = F.find_all_ocblocks(toks, [ocp])

    @test blocks[1].name == :CODE_BLOCK_IND
    @test F.content(blocks[1]) == "B\n    C"

    mds = """
        @def indented_code = true
        A

            B
            C
        D"""
    steps = explore_md_steps(mds)
    toks = steps[:tokenization].tokens
    @test toks[4].name == toks[5].name == :LR_INDENT
    blk = steps[:ocblocks].blocks
    @test blk[2].name == :CODE_BLOCK_IND
    b2i = steps[:b2insert].b2insert
    @test b2i[2].name == :CODE_BLOCK_IND
    @test isapproxstr(mds |> fd2html_td, """

A

B
        C

D

        """)
end

@testset "ind+lx" begin
    s = raw"""
        \newcommand{\julia}[1]{
            ```julia
            #1
            ```
        }
        Hello
        \julia{a=5
        x=3}
        """ |> fd2html_td
    @test isapproxstr(s, """

Hello
a=5
        x=3

        """)
end


text:
@auto_hash_equals struct AccessRecordTable
    access::Vector{UInt}
    time::Vector{TimeNS}
    threadid::Vector{Int}
    tag::Vector{UInt}
    referenceid::Vector{Int}
end

"""
    MemoryAccesses.table(record = MemoryAccesses.view(); [transform])

Convert the memory access `record` to a table.

# Keyword Arguments

- `transform::Symbol = :collapse`: Following options are available:

  * `:collapse` (default) means to shift all pointers to the base pointers
    recorded via `MemoryAccesses.reference`. Roughly speaking, all pointers are
    "collapsed" to indices (of arrays of 1-byte objects).

  * `:pack` is similar to `:collapse` but pointers coming from different arrays
    can still be distinguished. The base pointer from the `i+1`-th array
    recorded via `MemoryAccesses.reference` is mapped to the access location
    next to the last pointer of the `i`-th array.

  * `:none` does not transform anything.
"""
MemoryAccesses.table

function transform_option(transform::Symbol)
    if transform === :collapse
        return collapse_pointers
    elseif transform === :pack
        return pack_pointers
    elseif transform === :none
        # TODO: this should still set `referenceid`
        return Returns(raw_accessid)
    else
        error("invalid option: transform = $transform")
    end
end

function MemoryAccesses.table(record = MemoryAccesses.view(); transform::Symbol = :collapse)
    accesses = record.accesses
    references = record.references
    accessid_factory = transform_option(transform)

    n = sum(length, accesses)
    table = AccessRecordTable(
        Vector{UInt}(undef, n),
        Vector{TimeNS}(undef, n),
        Vector{Int}(undef, n),
        Vector{UInt}(undef, n),
        Vector{Int}(undef, n),
    )

    fill_table!(accessid_factory(references), table, accesses)

    p = sortperm(table.time; alg = MergeSort)
    permute!(table.access, p)
    permute!(table.time, p)
    permute!(table.threadid, p)
    permute!(table.referenceid, p)

    return table
end

function fill_table!(accessid, table, accesses)
    j = 1
    for (tid, a) in enumerate(accesses)
        for i in eachindex(a)
            @inbounds begin
                r = a[i]
                table.access[j], table.referenceid[j] = accessid(r.ptr)
                table.time[j] = r.time
                table.threadid[j] = tid
                table.tag[j] = r.tag
            end
            j += 1
        end
    end
end

@inline raw_accessid(ptr) = (ptr, 0)

"""
    collapse_pointers(references) -> ptr -> (access, referenceid)

Compile a ptr-to-(access, metadata) mapping from an iterable of
`ReferenceRecord`-like values.
"""
function collapse_pointers(references)
    accessid = raw_accessid
    id = 1
    for refrec in references
        accessid = let id = id, fallback = accessid, accessid
            @inline function accessid(ptr)
                if refrec.first <= ptr <= refrec.last
                    return (ptr - refrec.first, id)
                else
                    return fallback(ptr)
                end
            end
        end
        id += 1
    end
    return accessid
end

function pack_pointers(references)
    offset = 0
    accessid = raw_accessid
    for (id, refrec) in enumerate(references)
        accessid =
            let id = id,
                offset = offset,
                base = refrec.first,
                bound = refrec.last + 8,
                fallback = accessid,
                accessid

                @inline function accessid(ptr)
                    if base <= ptr < bound
                        return (ptr - base + offset, id)
                    else
                        return fallback(ptr)
                    end
                end
            end
        offset += (refrec.last + 8) - refrec.first
    end
    return accessid
end

function Base.show(io::IO, ::MIME"text/plain", table::AccessRecordTable)
    print(io, MemoryAccesses)
    print(io, ".table()")
    println(io)
    print(io, "  #records: ", length(table.access))
    println(io)
    print(io, "  duration: ", table.time[end] - table.time[1], " ns")
end

const COLUMN_NAMES = fieldnames(AccessRecordTable)
const COLUMN_TYPES = Tuple{fieldtypes(AccessRecordTable)...}

Tables.istable(::Type{AccessRecordTable}) = true
Tables.columnaccess(::Type{AccessRecordTable}) = true
Tables.columns(table::AccessRecordTable) = table
Tables.schema(::AccessRecordTable) = Tables.Schema(COLUMN_NAMES, COLUMN_TYPES)


text:
include("preconditioners.jl")
include("optimization_loops.jl")
include("Losses.jl")
include("Iterations.jl")
include("equations.jl")
include("containers.jl")
using Flux.Optimise

function basic_optimization(params::optimization_parameters)
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, B, parameters = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    # initial_condition_number = measure_condition_number(parameters, prec)
    loss_function(x::AbstractArray{<:Real, 1}) = prec_loss_function(x, prec, params)
    history = standard_optimization!(parameters, loss_function, params.epochs, params.optim)
    # final_condition_number = measure_condition_number(parameters, prec)
    return history, prec, parameters # (initial_condition_number, final_condition_number), history, prec, parameters
end

function stopping_optimization(params::optimization_parameters; check_time=10, tol=1e-2)
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, B, parameters = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    pre_opt, post_opt = deepcopy(params.optim), deepcopy(params.optim)

    # θ = parameters[end:end]
    # M = assemble_matrix(parameters, prec)
    # pre_loss_function(x::AbstractArray{<:Real, 1}) = scalar_loss_function(x, M, prec, params)
    # history_pre = early_stopping_optimization!(θ, pre_loss_function, params.epochs, pre_opt; check_time=check_time, tol=tol)
    # parameters[end] = θ[1]

    loss_function(x::AbstractArray{<:Real, 1}) = prec_loss_function(x, prec, params)
    history = early_stopping_optimization!(parameters, loss_function, params.epochs, params.optim; check_time=check_time, tol=tol)

    # θ = parameters[end:end]
    # M = assemble_matrix(parameters, prec)
    # post_loss_function(x::AbstractArray{<:Real, 1}) = scalar_loss_function(x, M, prec, params)
    # history_post = early_stopping_optimization!(θ, post_loss_function, params.epochs, post_opt; check_time=check_time, tol=tol)
    # parameters[end] = θ[1]

    return history, prec, parameters
end

function stopping_optimization_continuation(params::optimization_parameters, parameters::Array{Float64, 1}; check_time=10, tol=1e-2)
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, _, _ = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)

    loss_function(x::AbstractArray{<:Real, 1}) = prec_loss_function(x, prec, params)
    history = early_stopping_optimization!(parameters, loss_function, params.epochs, params.optim; check_time=check_time, tol=tol)

    return history, prec, parameters
end

function modified_optimization(params::optimization_parameters)
    # A = params.equation(params.L, params.BCs, params.eq_params)
    # prec, B, parameters = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    # initial_condition_number = measure_s_condition_number(parameters, prec, params)
    # loss_function(x::AbstractArray{<:Real, 1}) = spectral_loss_function(x, prec, params)
    # history = standard_optimization!(parameters, loss_function, params.epochs, params.optim)
    # final_condition_number = measure_s_condition_number(parameters, prec, params)
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, B, parameters = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    loss_function(x::AbstractArray{<:Real, 1}) = prec_loss_function(x, prec, params)
    ρ(x::AbstractArray{<:Real, 1}) = spectral_radius(x, prec, params)
    history = spectral_optimization!(parameters, loss_function, ρ, params.epochs, params.optim)
    return history, prec, parameters
end

function modified_optimization_continuation(params::optimization_parameters, parameters::Array{Float64, 1})
    # A = params.equation(params.L, params.BCs, params.eq_params)
    # prec, B, parameters = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    # initial_condition_number = measure_s_condition_number(parameters, prec, params)
    # loss_function(x::AbstractArray{<:Real, 1}) = spectral_loss_function(x, prec, params)
    # history = standard_optimization!(parameters, loss_function, params.epochs, params.optim)
    # final_condition_number = measure_s_condition_number(parameters, prec, params)
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, B, _ = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    loss_function(x::AbstractArray{<:Real, 1}) = prec_loss_function(x, prec, params)
    ρ(x::AbstractArray{<:Real, 1}) = spectral_radius(x, prec, params)
    history = spectral_optimization!(parameters, loss_function, ρ, params.epochs, params.optim)
    return history, prec, parameters
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
"""
    blockmean(cmd0::String="", arg1=nothing; kwargs...)

Block average (x,y,z) data tables by L2 norm.

Full option list at [`blockmean`]($(GMTdoc)blockmean.html)

Parameters
----------

- $(GMT.opt_R)
- **I** | **inc** :: [Type => Str | Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)blockmean.html#i)
- **A** | **fields** :: [Type => Str]

    Select which fields to write to individual grids. Append comma-separated codes for available
    fields: **z** (the mean data z, but see **-S**), **s** (standard deviation), **l** (lowest value),
    **h** (highest value) and **w** (the output weight; requires **-W**). [Default is just **z**].
    ($(GMTdoc)blockmean.html#a)
- **C** | **center** :: [Type => Bool]

    Use the center of the block as the output location [Default uses the mean location]. Not used when **-A**
    ($(GMTdoc)blockmean.html#c)
- **E** | **extend** :: [Type => Str | []]

    Provide Extended report which includes s (the standard deviation about the mean), l, the lowest
    value, and h, the high value for each block. Output order becomes x,y,z,s,l,h[,w]. [Default
    outputs x,y,z[,w]. See -W for w output. If -Ep is used we assume weights are 1/(sigma squared)
    and s becomes the propagated error of the mean.
    ($(GMTdoc)blockmean.html#e)
- **G** | **outgrid** :: [Type => Str | []]

    Write one or more fields directly to grids on disk; no table data are return. If more than one
    fields are specified via **A** then grdfile must contain the format flag %s so that we can embed the
    field code in the file names. If not provided but **A** is used, return 1 or more GMTgrid type(s).
    ($(GMTdoc)blockmean.html#g)
- **S** | **npts** | **number_of_points** [Type => Str | []]

    Use S=:n to report the number of points inside each block, S=:s to report the sum of all z-values
    inside a block, S=:w to report the sum of weights [Default (or S=:m reports mean value].
    ($(GMTdoc)blockmean.html#s)
- **W** | **weights** :: [Type => Str | []]

    Unweighted input and output have 3 columns x,y,z; Weighted i/o has 4 columns x,y,z,w. Weights can
    be used in input to construct weighted mean values for each block.
    ($(GMTdoc)blockmean.html#w)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function blockmean(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("blockmean", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd = parse_these_opts("", d, [[:E :extended], [:S :npts :number_of_points]])
	return common_blocks(cmd0, arg1, d, cmd, "blockmean", kwargs...)
end

# ---------------------------------------------------------------------------------------------------
"""
    blockmedian(cmd0::String="", arg1=nothing; kwargs...)

Block average (x,y,z) data tables by L1 norm.

Full option list at [`blockmedian`]($(GMTdoc)blockmedian.html)
"""
function blockmedian(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("blockmedian", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd = parse_these_opts("", d, [[:E :extended], [:Q :quick], [:T :quantile]])
	return common_blocks(cmd0, arg1, d, cmd, "blockmedian", kwargs...)
end

# ---------------------------------------------------------------------------------------------------
"""
    blockmode(cmd0::String="", arg1=nothing; kwargs...)

Block average (x,y,z) data tables by mode estimation.

Full option list at [`blockmode`]($(GMTdoc)blockmode.html)
"""
function blockmode(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("blockmode", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd = parse_these_opts("", d, [[:E :extended], [:D :histogram_binning], [:Q :quick]])
	return common_blocks(cmd0, arg1, d, cmd, "blockmode", kwargs...)
end

# ---------------------------------------------------------------------------------------------------
function common_blocks(cmd0, arg1, d, cmd, proggy, kwargs...)

	cmd = parse_these_opts(cmd, d, [[:A :fields], [:C :center], [:G :outgrid :grid], [:I :inc], [:W :weights]])
	if (occursin("-G", cmd) && !occursin("-A", cmd))
		cmd = cmd * " -Az"					# So that we can use plain -G to mean write grid
	end
	cmd, = parse_common_opts(d, cmd, [:R :V_params :bi :di :e :f :h :i :o :r :yx])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)
    (occursin("-G", cmd)) && (cmd = replace(cmd, "-G" => ""))   # GMT API does not allow a -G from externals
	return common_grd(d, proggy * " " * cmd, arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
blockmean(arg1, cmd0::String=""; kw...) = blockmean(cmd0, arg1; kw...)
# ---------------------------------------------------------------------------------------------------
blockmedian(arg1, cmd0::String=""; kw...) = blockmedian(cmd0, arg1; kw...)
# ---------------------------------------------------------------------------------------------------
blockmode(arg1, cmd0::String=""; kw...) = blockmode(cmd0, arg1; kw...)


text:
include("varset.jl")
include("continuousvar.jl")
include("discretevar.jl")

initv(T) = typemin(T)

#=
function initv(T)
    T == Float16 && return NaN16
    T == Float32 && return NaN32
    T == Float64 && return NaN
    typemin(T)
end
=#

Base.values(v::AbstractVariable) = v.vpop

#model(v) = v.modl

function update_variable!(v, b, s, e)
    # Store best value
    v.bstv = v.vpop[b]
    # Fit model with the best selected individuals (solutions)
    v.modl = fit_model(v, s)
    # Retain the elite individuals for the next generation
    v.vpop = v.vpop[e]
    v
end

Base.rand(v::AbstractVariable) = convert(eltype(v.vpop), v.modl |> rand)


text:
module Flu

using TreeKnit
using TreeTools

#
export nt_seq_lengths
export read_flu_trees

# Constants

const flu_lineages = ["h3n2"]
const flu_segments = ["ha", "na", "pb1", ]

const nt_seq_lengths = Dict(
	"h3n2" => Dict(
		"ha" => 1698,
		"na" => 1410,
		"pb1" => 2311,
	),
)

const substitution_rates = Dict(
	"h3n2" => Dict(
		"ha" => 3.66e-3,
		"na" => 2.87e-3,
	)
)

#
"""
	read_flu_trees(files::Dict, lineage = "h3n2")

Read trees stored in `values(files)` and return them as a dictionary indexed with
  `keys(files)`.
  Short branches are removed from trees.
"""
function read_flu_trees(files::AbstractDict, lineage = "h3n2")
	if !in(lineage, flu_lineages)
		error("Unknown flu lineage $lineage. Options are $(flu_lineages)")
	end

	trees = Dict{String, Tree}()
	for (segment, f) in files
		if !in(segment, flu_segments)
			error("Unknown flu segment $segment. Options are $(flu_segments)")
		end
		trees[segment] = read_tree(f)
		L = nt_seq_lengths[lineage][segment]
		TreeTools.delete_null_branches!(trees[segment]; threshold = 1/L/2)
	end

	return trees
end

end # module

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
##############################################################################
##
## Allocations for AllocatedDogleg
##
##############################################################################

struct AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2} <: AbstractAllocatedOptimizer
    δgn::Tx1
    δgr::Tx2
    δx::Tx3
    dtd::Tx4
    ftrial::Ty1
    fpredict::Ty2
    function AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}(δgn, δgr, δx, dtd, ftrial, fpredict) where {Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}
        length(δgn) == length(δgr) || throw(DimensionMismatch("The lengths of δgn and δgr must match."))
        length(δgn) == length(δx) || throw(DimensionMismatch("The lengths of δgn and δx must match."))
        length(δgn) == length(dtd) || throw(DimensionMismatch("The lengths of δgn and dtd must match."))
        length(ftrial) == length(fpredict) || throw(DimensionMismatch("The lengths of ftrial and fpredict must match."))
        new(δgn, δgr, δx, dtd, ftrial, fpredict)
    end
end

function AllocatedDogleg(δgn::Tx1, δgr::Tx2, δx::Tx3, dtd::Tx4, ftrial::Ty1, fpredict::Ty2) where {Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}
    AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}(δgn, δgr, δx, dtd, ftrial, fpredict)
end

function AbstractAllocatedOptimizer(nls::LeastSquaresProblem, optimizer::Dogleg)
    AllocatedDogleg(_zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x),
    _zeros(nls.y), _zeros(nls.y))
end

##############################################################################
##
## Method for AllocatedDogleg
##
##############################################################################

const MIN_Δ = 1e-16 # maximum trust region radius
const MAX_Δ = 1e16 # minimum trust region radius
const MIN_STEP_QUALITY = 1e-3
const MIN_DIAGONAL = 1e-6
const MAX_DIAGONAL = 1e32
const DECREASE_THRESHOLD = 0.25
const INCREASE_THRESHOLD = 0.75

function optimize!(
    anls::LeastSquaresProblemAllocated{Tx, Ty, Tf, TJ, Tg, Toptimizer, Tsolver};
    x_tol::Number = 1e-8, f_tol::Number = 1e-8, g_tol::Number = 1e-8,
    iterations::Integer = 1_000, Δ::Number = 1.0, store_trace = false, show_trace = false, show_every = 1, lower = Vector{eltype(Tx)}(undef, 0), upper = Vector{eltype(Tx)}(undef, 0)) where {Tx, Ty, Tf, TJ, Tg, Toptimizer <: AllocatedDogleg, Tsolver}

     δgn, δgr, δx, dtd = anls.optimizer.δgn, anls.optimizer.δgr, anls.optimizer.δx, anls.optimizer.dtd
     ftrial, fpredict = anls.optimizer.ftrial, anls.optimizer.fpredict
     x, fcur, f!, J, g! = anls.x, anls.y, anls.f!, anls.J, anls.g!

     #istempty
     ((isempty(lower) || length(lower)==length(x)) && (isempty(upper) || length(upper)==length(x))) ||
             throw(ArgumentError("Bounds must either be empty or of the same length as the number of parameters."))
     ((isempty(lower) || all(x .>= lower)) && (isempty(upper) || all(x .<= upper))) ||             throw(ArgumentError("Initial guess must be within bounds."))

    # initialize
    reuse = false
    wnorm_δgn = 0.0
    wnorm_δgr = 0.0
    α = 0.0
    f_calls,  g_calls, mul_calls = 0, 0, 0
    converged, x_converged, f_converged, g_converged, converged =
        false, false, false, false, false
    f!(fcur, x)
    f_calls += 1
    ssr = sum(abs2, fcur)
    maxabs_gr = Inf

    iter = 0
    eTx, eTy = eltype(x), eltype(fcur)

    tr = OptimizationTrace()
    tracing = store_trace || show_trace
    tracing && update!(tr, iter, ssr, maxabs_gr, store_trace, show_trace, show_every)
    while !converged && iter < iterations
        iter += 1
        check_isfinite(x)
        # compute step
        if !reuse
            #update gradient
            g!(J, x)
            g_calls += 1
            colsumabs2!(dtd, J)
            clamp!(dtd, MIN_DIAGONAL, MAX_DIAGONAL)

            if iter == 1
                wnorm_x = wnorm(x, dtd)
                if wnorm_x > 0
                    Δ *= wnorm_x
                end
            end
            # compute (opposite) gradient
            mul!(δgr, J', fcur, one(eTx), zero(eTx))
            mul_calls += 1
            maxabs_gr = maximum(abs, δgr)
            wnorm_δgr = wnorm(δgr, dtd)

            # compute Cauchy point
            map!((x, y) -> x * sqrt(y), δgn, δgr, dtd)
            mul!(fpredict, J, δgn, one(eTy), zero(eTy))
            mul_calls += 1
            α = wnorm_δgr^2 / sum(abs2, fpredict)

            # compute Gauss Newton step δgn
            fill!(δgn, zero(eTx))
            δgn, ls_iter = ldiv!(δgn, J, fcur, anls.solver)
            mul_calls += ls_iter
            wnorm_δgn = wnorm(δgn, dtd)
        end
        # compute δx
        if wnorm_δgn <= Δ
            #  Case 1. The Gauss-Newton step lies inside the trust region
            copyto!(δx, δgn)
            wnorm_δx = wnorm_δgn
        elseif wnorm_δgr * α >= Δ
            # Case 2. The Cauchy point and the Gauss-Newton steps lie outside
            # the trust region.
            # rescale Cauchy step within the trust region and return
            copyto!(δx, δgr)
            rmul!(δx, Δ / wnorm_δgr)
            wnorm_δx = Δ
        else
            # Case 3. The Cauchy point is inside the trust region nd the
            # Gauss-Newton step is outside
            b_dot_a = α * wdot(δgr, δgn, dtd)
            a_squared_norm = (α * wnorm_δgr)^2
            b_minus_a_squared_norm =
                  a_squared_norm - 2 * b_dot_a + wnorm_δgn^2
            c =  b_dot_a - a_squared_norm
            d = sqrt(c^2 + b_minus_a_squared_norm * (Δ^2 - a_squared_norm))
            β = (c <= 0) ? (d - c)/ b_minus_a_squared_norm : (Δ^2 - a_squared_norm) / (d + c)
            copyto!(δx, δgn)
            rmul!(δx, β)
            axpy!(α * (1 - β), δgr, δx)
            wnorm_δx = wnorm(δx, dtd)
        end

        # apply box constraints
        if !isempty(lower)
            @simd for i in 1:length(x)
               @inbounds δx[i] = min(δx[i], x[i] - lower[i])
            end
        end
        if !isempty(upper)
            @simd for i in 1:length(x)
               @inbounds δx[i] = max(δx[i], x[i] - upper[i])
            end
        end

        #update x
        axpy!(-one(eTx), δx, x)

        f!(ftrial, x)
        f_calls += 1

        # trial ssr
        trial_ssr = sum(abs2, ftrial)

        # predicted ssr
        mul!(fpredict, J, δx, one(eTx), zero(eTx))
        mul_calls += 1
        axpy!(-one(eTy), fcur, fpredict)
        predicted_ssr = sum(abs2, fpredict)

        ρ = (ssr - trial_ssr) / (ssr - predicted_ssr)
        x_converged, f_converged, g_converged, converged =
            assess_convergence(δx, x, maxabs_gr, ssr, trial_ssr, x_tol, f_tol, g_tol)

        if ρ >= MIN_STEP_QUALITY
            # Successful iteration
            reuse = false
            copyto!(fcur, ftrial)
            ssr = trial_ssr
        else
            # unsucessful iteration
            reuse = true
            axpy!(one(eTx), δx, x)
        end

        if ρ < DECREASE_THRESHOLD
           Δ = max(MIN_Δ, Δ * 0.5)
        elseif ρ > INCREASE_THRESHOLD
           Δ = max(Δ, 3.0 * wnorm_δx)
       end
       tracing && update!(tr, iter, ssr, maxabs_gr, store_trace, show_trace, show_every)
    end
    LeastSquaresResult("Dogleg", x, ssr, iter, converged,
                        x_converged, x_tol, f_converged, f_tol, g_converged, g_tol, tr,
                        f_calls, g_calls, mul_calls)
end

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
@testset "ForwardShare" begin
    @testset "testGetDataParallel.jl" begin
        include("testGetDataParallel.jl")
    end
	@testset "testPrepareMesh2Mesh.jl" begin
        include("testPrepareMesh2Mesh.jl")
    end
    @testset "testTest.jl" begin
        include("testTest.jl")
    end
    @testset "testGetSensMat.jl" begin
        include("testGetSensMat.jl")
    end
    @testset "testClear.jl" begin
        include("testClear.jl")
    end
end


text:
# Implement the ScikitLearnBase.jl interface

import ScikitLearnBase

ScikitLearnBase.is_classifier(::GPE) = false

ScikitLearnBase.fit!(gp::GPE, X::AbstractMatrix, y::AbstractVector) = fit!(gp, permutedims(X), y)

function ScikitLearnBase.predict(gp::GPE, X::AbstractMatrix; eval_MSE::Bool=false)
    mu, Sigma = predict_y(gp, X'; full_cov=false)
    if eval_MSE
        return mu, Sigma
    else
        return mu
    end
end

# This is a default - arbitrary scoring functions can be passed to GridSearchCV
function ScikitLearnBase.score(gp::GPE, x, y)
    # I use MSE here for simplicity. Is that equivalent to the likelihood for
    # gaussian processes?
    # It would make sense to use the marginal likelihood for hyperparameter
    # tuning (like optimize.jl does), but GridSearchCV isn't designed for that.
    return -mean((ScikitLearnBase.predict(gp, x) - y) .^ 2)
end

ScikitLearnBase.clone(gp::GPE) = GPE(; m=ScikitLearnBase.clone(gp.m),                                   k=ScikitLearnBase.clone(gp.k),
                                   logNoise=gp.logNoise)

# Means and Kernels are small objects, and they are not fit to the data (except
# through hyperparameter tuning), so `deepcopy` is appropriate (if a bit lazy)
ScikitLearnBase.clone(m::Mean) = deepcopy(m)
ScikitLearnBase.clone(k::Kernel) = deepcopy(k)

################################################################################
## GaussianProcesses.jl has a get_params function, which returns a
## Vector{Float64}. We use it to implement ScikitLearnBase.get_params, which
## must return a Symbol => value dictionary. Likewise for set_params!

function add_prefix(pref, di)
    newdi = typeof(di)()
    for (param,value) in di
        newdi[Symbol(pref, param)] = value
    end
    return newdi
end

function ScikitLearnBase.get_params(gp::GPE)
    merge(add_prefix(:m_, ScikitLearnBase.get_params(gp.m)),
          add_prefix(:k_, ScikitLearnBase.get_params(gp.k)),
          Dict(:logNoise=>gp.logNoise))
end

function ScikitLearnBase.get_params(obj::Union{Mean, Kernel})
    params = get_params(obj) # GaussianProcesses' get_params
    names = get_param_names(obj)
    @assert length(params) == length(names) # sanity check
    return Dict(zip(names, params))
end

function ScikitLearnBase.set_params!(gp::GPE; params...)
    m_params = Dict()
    k_params = Dict()
    for (name, value) in params
        sname = string(name)
        if startswith(sname, "m_")
            m_params[Symbol(sname[3:end])] = value
        elseif startswith(sname, "k_")
            k_params[Symbol(sname[3:end])] = value
        else
            @assert name == :logNoise "Unknown parameter passed to set_params!: $name"
            gp.logNoise = value
        end
    end

    ScikitLearnBase.set_params!(gp.m; m_params...)
    ScikitLearnBase.set_params!(gp.k; k_params...)
    gp
end

function ScikitLearnBase.set_params!(obj::Union{Mean, Kernel}; params...)
    # In ScikitLearn, params need not contain all the parameters, so the
    # code is different from get_params
    names = get_param_names(obj)
    # Get current parameter vector
    hyp = copy(get_params(obj))
    # Update it
    for (name, value) in params
        ind = findfirst(names, name)
        @assert ind > 0 "Unknown parameter passed to set_params!: $name"
        hyp[ind] = value
    end
    set_params!(obj, hyp)
    obj
end


text:
"""
$(TYPEDEF)

Dynamic pose variable with velocity components: `x, y, theta, dx/dt, dy/dt`
"""
mutable struct DynPose2 <: IncrementalInference.InferenceVariable
  ut::Int64 # microsecond time
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol,Symbol,Symbol,Symbol,Symbol}
  DynPose2(;ut::Int64=-9999999999, labels::Vector{<:AbstractString}=String[]) = new(ut, 5, labels, (:Euclid,:Euclid,:Circular,:Euclid,:Euclid))
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 998407c2-ed46-11ea-1c81-e33bc9e340f7
using Pkg, DrWatson

# ╔═╡ c30a6ffc-ed46-11ea-1036-57f7da453ef2
begin
	#@quickactivate "StatisticalRethinkingStan"
	using Distributions
	using StatisticalRethinking
end

# ╔═╡ 115d3284-ed4c-11ea-366e-33869314ff98
md"## Distributions.jl"

# ╔═╡ c3119548-ed46-11ea-34e2-efd188f0457c
# snippet 2.1
begin
	ways = [0, 3, 8, 9, 0]
	ways = ways / sum(ways)
end

# ╔═╡ edb843be-ed46-11ea-1f85-6b7a7bff28ec
md"##### Working with distributions is a little different between Julia and R.

Instead of having `rbinom`, `dbinom` and `pbinom` we just
have the `Binomial` distribution which to work with:"

# ╔═╡ c31a57f2-ed46-11ea-1303-c3bf7b3b7c8f
d = Binomial(9, 0.5)

# ╔═╡ a0b518c8-ed51-11ea-06b1-cd86fe8d1229
md"##### Parameters of Binomial."

# ╔═╡ 5f0e19e6-ed51-11ea-3130-ebd24a5f62bc
fieldnames(Binomial)

# ╔═╡ 3be503f6-ed51-11ea-0958-5fa25106d3f5
md"##### Number of trials parameter."

# ╔═╡ c3272656-ed46-11ea-3971-a3cc353218cb
d.n

# ╔═╡ c872e37e-ed47-11ea-141c-f7a43395df4d
md"##### Singe random draw."

# ╔═╡ c22522ca-ed47-11ea-1fb8-a7e582a8969e
rand(d)

# ╔═╡ 5ff49ae0-ed4c-11ea-0994-ebf5f6c17922
md"##### 9 draws."

# ╔═╡ c3298d88-ed46-11ea-1a28-578706bf1724
rand(d, 9)                          #

# ╔═╡ c335d4ee-ed46-11ea-3c3b-637d1a3d5bc8
pdf(d, 6)                        # probability density of getting a 6

# ╔═╡ c3403d4e-ed46-11ea-3dce-7906235096c6
cdf(d, 6)                        # cumulative probability of getting 6

# ╔═╡ c354cad4-ed46-11ea-1244-9d97d407103f
md"## End of distributions.jl"

# ╔═╡ Cell order:
# ╟─115d3284-ed4c-11ea-366e-33869314ff98
# ╠═998407c2-ed46-11ea-1c81-e33bc9e340f7
# ╠═c30a6ffc-ed46-11ea-1036-57f7da453ef2
# ╠═c3119548-ed46-11ea-34e2-efd188f0457c
# ╟─edb843be-ed46-11ea-1f85-6b7a7bff28ec
# ╠═c31a57f2-ed46-11ea-1303-c3bf7b3b7c8f
# ╟─a0b518c8-ed51-11ea-06b1-cd86fe8d1229
# ╠═5f0e19e6-ed51-11ea-3130-ebd24a5f62bc
# ╟─3be503f6-ed51-11ea-0958-5fa25106d3f5
# ╠═c3272656-ed46-11ea-3971-a3cc353218cb
# ╟─c872e37e-ed47-11ea-141c-f7a43395df4d
# ╠═c22522ca-ed47-11ea-1fb8-a7e582a8969e
# ╟─5ff49ae0-ed4c-11ea-0994-ebf5f6c17922
# ╠═c3298d88-ed46-11ea-1a28-578706bf1724
# ╠═c335d4ee-ed46-11ea-3c3b-637d1a3d5bc8
# ╠═c3403d4e-ed46-11ea-3dce-7906235096c6
# ╟─c354cad4-ed46-11ea-1244-9d97d407103f


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG

const mrg63norm = 1.0842021724855052e-19
const mrg63m1 = 9223372036854769163
const mrg63m2 = 9223372036854754679
const mrg63a12 = 1754669720
const mrg63q12 = 5256471877
const mrg63r12 = 251304723
const mrg63a13n = 3182104042
const mrg63q13 = 2898513661
const mrg63r13 = 394451401
const mrg63a21 = 31387477935
const mrg63q21 = 293855150
const mrg63r21 = 143639429
const mrg63a23n = 6199136374
const mrg63q23 = 1487847900
const mrg63r23 = 985240079
const mrg63seed = 123456789

mutable struct MRG63k3a <: AbstractRNG{UInt64}
    s::NTuple{6,Int64} #The seeds for s[1:3] must be integers in [0, m1 - 1] and not all 0.
    #The seeds for s[4:6] must be integers in [0, m2 - 1] and not all 0.
end

function MRG63k3a()
    return MRG63k3a((mrg63seed, mrg63seed, mrg63seed, mrg63seed, mrg63seed, mrg63seed))
end

@inline rand(rng::MRG63k3a, ::Type{Float64}) = rand(rng, UInt64) * mrg63norm

#63-bits of randomness
@inline function rand(rng::MRG63k3a, ::Type{UInt64})

    h = rng.s[1] ÷ mrg63q13
    p13 = mrg63a13n * (rng.s[1] - h * mrg63q13) - h * mrg63r13
    h = rng.s[2] ÷ mrg63q12
    p12 = mrg63a12 * (rng.s[2] - h * mrg63q12) - h * mrg63r12
    if (p13 < 0)
        p13 += mrg63m1
    end
    if (p12 < 0)
        p12 += mrg63m1 - p13
    else
        p12 -= p13
    end
    if (p12 < 0)
        p12 += mrg63m1
    end
    h = rng.s[4] ÷ mrg63q23
    p23 = mrg63a23n * (rng.s[4] - h * mrg63q23) - h * mrg63r23
    h = rng.s[6] ÷ mrg63q21
    p21 = mrg63a21 * (rng.s[6] - h * mrg63q21) - h * mrg63r21
    if (p23 < 0)
        p23 += mrg63m2
    end
    if (p21 < 0)
        p21 += mrg63m2 - p23
    else
        p21 -= p23
    end
    if (p21 < 0)
        p21 += mrg63m2
    end
    rng.s = (rng.s[2], rng.s[3], p12, rng.s[5], rng.s[6], p21)
    if (p12 > p21)
        return convert(UInt64, p12 - p21)
    end
    return convert(UInt64, p12 - p21 + mrg63m1)
end

# mtx1 :: JumpMatrix Int64
# mtx1 = JM (0, 1, 0) (0, 0, 1) (m1 - a13n, a12, 0)
#
# mtx2 :: JumpMatrix Int64
# mtx2 = JM (0, 1, 0) (0, 0, 1) (m2 - a23n, 0, a21)


text:
module TryCatch
    export @try

    specials = (:(var"@catch"), :(var"@else"), :(var"@finally"))

    "
    Test if an expression is one of these special annotations
    "
    isspecial(ex) = begin
        if ex isa Expr && ex.head == :macrocall && ex.args[1] ∈ specials
            true
        else
            false
        end
    end

    "
    Copy an expression-like object
    "
    copyex(ex::Expr) = copy(ex)
    copyex(ex) = ex

    "
    A truncate-able one-line string representation of an expression for display purposes
    "
    strexpr(ex, truncate=nothing) = begin
        exstr = string(linenumberremove(ex))
        exstr = join(strip.(split(exstr, "\n")), " ")
        if truncate !== nothing && lastindex(exstr) > truncate-3
            exstr = join([i for i in exstr][begin:truncate])*"..."
        end
        return exstr
    end

    "
    Given a code block labelled with @catch, @else and/or @finally annotation, like `@try 3 @catch e @else 21 @finally 1`, unnest these inner macros-like
    annotations from the outer nests in order to achieve a sequential structure and not a nested one.
    "
    unnestmacrolabels(exs) = begin
        exs = [copyex(i) for i in exs]

        i=0; while (i+=1)<=length(exs)
            if isspecial(exs[i])
               for j ∈ 2:length(exs[i].args)
                    if isspecial(exs[i].args[j])
                        # append the inner nest to underneath this macro
                        exs[i].args, leftovers = exs[i].args[1:j-1], exs[i].args[j:end]
                        for k ∈ 1:length(leftovers)
                            insert!(exs, i+k, leftovers[k])
                        end

                        # advance i to continue after this insersion
                        i+=length(leftovers)-1
                    end
                end
            end
        end

        return exs
    end

    "
    Given the expression send to the @try macro as an array of 'lines-of-code', sort them into the following groups
    groups: :head, @catch, @else and @finally
    "
    split_try_labels_into_groups(exs) = begin

        groups = [:head, specials...]
        d = Dict()

        for g in groups
            d[g] = []
        end

        gᵢ = 1
        macroname = :(var"@try")
        seen = Set([:head])
        for i ∈ 1:length(exs)
            ex = exs[i]

            if ex isa Expr && ex.head == :macrocall && ex.args[1] in specials
                macronameᵢ = ex.args[1]
                if macronameᵢ ∈ groups[gᵢ:end]
                    gᵢ = findfirst(x->(macronameᵢ==x), groups)

                    # may have multiple catch blocks
                    if groups[gᵢ] == :(var"@catch")
                        #push!(d[groups[gᵢ]], [])

                        if length(ex.args) >=5
                            throw(ErrorException("syntax: $(string(groups[gᵢ])) unexpected code block `$(strexpr(ex.args[5], 23))`"))
                        end

                        # linenumber+condition => block
                        if length(ex.args) < 3
                            cond = (ex.args[2], Expr(:->, gensym(:e), true))
                            args = Vector{Any}(ex.args[2:end])
                        else
                            cond = (ex.args[2], ex.args[3])
                            args = Vector{Any}([ex.args[2], ex.args[4:end]...])
                        end

                        # Also push the linenumbernode
                        push!(d[groups[gᵢ]], cond=>args)

                    # may have single other blocks
                    else
                        if macronameᵢ == macroname
                            throw(ErrorException("syntax: unexpected $(macronameᵢ) following $(macroname)"))
                        end

                        if length(ex.args) >=4
                            throw(ErrorException("syntax: $(string(groups[gᵢ])) unexpected code block `$(strexpr(ex.args[4], 23))`"))
                        end

                        # Also push the linenumbernode
                        append!(d[groups[gᵢ]], ex.args[2:end])

                    end
                else
                    throw(ErrorException("syntax: unexpected $(macronameᵢ) following $(macroname)"))
                end

                macroname = macronameᵢ
                push!(seen, macroname)

            else
                if groups[gᵢ] == :(var"@catch")
                    push!(d[groups[gᵢ]][end][2], ex)
                else
                    push!(d[groups[gᵢ]], ex)
                end
            end
        end

        # Don't record unseen macro labels that was never seen
        for key in keys(d)
            if !(key ∈ seen)
                d[key] = nothing
            end
        end

        return d
    end

    "
    Function to generate if-elseif-else type od expressions
    "
    ifgenerator(conds=nothing, else_=nothing) = begin
        if conds === nothing || length(conds) == 0
            return else_
        end

        (x, y) = first(conds)
        orig_expr = Expr(:if, x, y)

        expr = orig_expr
        for (i,(x,y)) in enumerate(conds)
            if i==1
                continue
            end
            push!(expr.args, Expr(:elseif, x, y))
            expr = expr.args[end]
        end

        if else_ !== nothing
            push!(expr.args, else_)
        end

        return orig_expr
    end

    "
    Replace all instances of a symbol within a given expression
    "
    symbolrename(ex::Expr, find, replace) = begin
        Expr((symbolrename(i, find, replace) for i in [ex.head; ex.args])...)
    end

    symbolrename(ex::Symbol, find, replace) = ex == find ? replace : ex

    symbolrename(ex::Expr, findreplace::Dict) = begin
        Expr((symbolrename(i, findreplace) for i in [ex.head; ex.args])...)
    end

    symbolrename(ex::Symbol, findreplace::Dict) = get(findreplace, ex, ex)

    symbolrename(ex, find, replace=nothing) = ex

    "
    Remove a symbol from expression
    "
    symbolremove(ex::Expr, slist) = begin
        Expr((symbolremove(i, slist) for i ∈ [ex.head; ex.args] if !(i ∈ slist))...)
    end
    symbolremove(ex, slist) = ex

    "
    Remove any LineNumberNodes from an expression
    "
    linenumberremove(ex::Expr) = begin
        if ex.head == :macrocall
            ex = copyex(ex)
            ex.args[2] = nothing
        end
        Expr((linenumberremove(i) for i ∈ [ex.head; ex.args] if !(i isa LineNumberNode))...)
    end
    linenumberremove(ex::LineNumberNode) = Expr(:block)
    linenumberremove(ex) = ex

    "
    Convert the condition-like expression of the catch block, e.g. `catch e isa MethodError`, into a proper boolean expression
    `catch e->e isa MethodError`
    "
    function conditionhelper(ex, linenumber=nothing)
        if ex isa Expr && ex.head == :->
            larg = ex.args[1]
            if !(larg isa Symbol)
                throw(ErrorException("@catch lambda condition must have a single argument like `e->true`, got `$(strexpr(ex, 23))`"))
            end

            return ex.args[2], larg

        elseif ex isa Expr && ex.head == :(::)
            if !(length(ex.args)==2) || !(ex.args[1] isa Symbol)
                throw(ErrorException("@catch condition must be in the form `e::TypeOfError`, got `$(strexpr(ex, 23))`"))
            end

            return Expr(:block,
                        linenumber,
                        :($(ex.args[1]) isa $(ex.args[2]))), ex.args[1]

        elseif ex isa Expr
            ex_copy = ex
            while !(ex_copy isa Symbol)
                if !(ex_copy isa Expr)
                    throw(ErrorException("@catch condition must have a symbol at the leftmost location, got `$(strexpr(ex_copy, 23))`"))
                end

                i₁ = ex_copy.head ∈ (:macrocall, :call, :ref) ? 2 : 1
                ex_copy = ex_copy.args[i₁]
            end

            return Expr(:block, linenumber, ex), ex_copy

        elseif ex isa Symbol
            return true, ex
        end

    end

    """
    A @try macro to mimickes how error handling is done in Python, using @catch @except and @finally as labels.

    Example:
    ```
    @try begin
        sqrt("34")
    @catch e::MethodError
        println("Oops")
    @finally
        println("Oh, well")
    end
    ```
    """
    macro try_ end

    # Rename  @try_ to @try in order to bypass reserved keyword restriction
    var"@try" = var"@try_"
    typeof(var"@try").name.mt.name = Symbol("@try")
    macro try_(exs...)
        # Check correct number of arguments
        if length(exs) >=2 && !isspecial(exs[2])
            throw(ErrorException("syntax: @try unexpected code block `$(strexpr(exs[2], 23))`"))
        end

        # Turn :block entry into a list
        if length(exs) == 1
            if exs[1] isa Expr && exs[1].head == :block
                exs = exs[1].args
            end
        end

        exs = unnestmacrolabels(exs)
        d = split_try_labels_into_groups(exs)

        # Convert the groups into expressions
        @gensym exception successful successresult
        conditions = []
        if d[Symbol("@catch")] !== nothing
            for ((linenum, cond), block) ∈ d[Symbol("@catch")]
                # refactor into proper conditionals
                cond, sym = conditionhelper(cond, linenum)

                cond = Expr(:block, linenum, cond)
                block = Expr(:block, block...)

                # find-and-replace the error symbol in block
                cond = symbolrename(cond, sym, exception)
                block = symbolrename(block, sym, exception)

                push!(
                    conditions,
                    cond => block
                )
            end
        end

        headₓₓ = Expr(:block, d[:head]...)
        catchₓₓ = ifgenerator(conditions, :(rethrow($exception)))
        elseₓₓ = if d[Symbol("@else")]===nothing
                     nothing
                 else
                     Expr(:(=), successresult, Expr(:if, successful, Expr(:block, d[Symbol("@else")]...)))
                 end

        finallyₓₓ = if d[Symbol("@finally")]===nothing
                        nothing
                    else
                        Expr(:block, d[Symbol("@finally")]...)
                    end

        @gensym headₓ elseₓ finallyₓ catchₓ
        template = :(
            $successresult = nothing;
            $successful=true;
            try
                $headₓ
            catch $exception
                $successful=false
                $catchₓ
            finally
                try
                    $elseₓ
                finally
                    $finallyₓ
                end
            end;
            $successresult
        )

        # Only return `successresults` if @else annotation is present
        if d[Symbol("@else")] === nothing
            template.args = template.args[1:end-1]
        end

        # Strip the Trycatch file lines from this macro
        template = linenumberremove(template)

        replacements = Dict(headₓ=>headₓₓ,
                            catchₓ=>catchₓₓ,
                            elseₓ=>elseₓₓ,
                            finallyₓ=>finallyₓₓ)

        # Remove any `nothing` expression
        template = symbolremove(template, [key for (key,val) ∈ replacements if val === nothing])

        # Fill the expression with
        template = symbolrename(template, replacements)

        esc(template)
    end
end


text:
"""
    writecopyconstructor(mod, T; evaluate=true)

Constructs and compiles a simple copy constructor for struct type `T`. Each field
is copied if a corresponding method for `copy` exists,
and is passed unaltered otherwise. The code creating the copy
constructor is evaluated in module `mod`. If `evaluate` is false,
the expression defining the copy method is returned and is not evaluted.

Note that the `copy` methods for the fields must exist at the time
`writecopyconstructor` is called.

# Example
```julia
struct A
 a::Int
 b::String
 c::Vector
end

julia> DataUtils.writecopyconstructor(Main, A)

julia> x = A(1, "hi", [1, 2])
A(1, "hi", [1, 2])

julia> y = copy(x)
A(1, "hi", [1, 2])

julia> x.c[1] = 7;

julia> x
A(1, "hi", [7, 2])

julia> y
A(1, "hi", [1, 2])
```
"""
function writecopyconstructor(_module, T::DataType; evaluate=true)
    _writecopyconstructor(_module, T, T, evaluate)
end

function writecopyconstructor(_module, T::UnionAll; evaluate=true)
    _writecopyconstructor(_module, T, T.body, evaluate)
end

function _writecopyconstructor(_module, T, Tdata, evaluate::Bool)
    isstructtype(T) || throw(ArgumentError("writecopyconstructor expects a struct type"))
    field_names = fieldnames(T)
    field_types = Tdata.types
    type_name =  Tdata.name
    args = Any[Symbol(type_name)]
    for (_name, _type) in zip(field_names, field_types)
        arg = if hasmethod(copy, (_type,))
            :(copy(obj.$_name))
        else
            :(obj.$_name)
        end
        push!(args, arg)
    end
    rhs = Expr(:call, args...)
    expr = Meta.parse(string(:(Base.copy(obj::$type_name) = $rhs)))
    if evaluate
        @eval _module $expr
        return nothing
    else
        return expr
    end
end


text:
# Deflate compression of empty string.
empty_deflate = [0x03, 0x00]

# Zlib compression of empty string.
empty_zlib = [0x78, 0x9c, 0x03, 0x00, 0x00, 0x00, 0x00, 0x01]

# Gzip compression of empty string with no headers other than header crc.
empty_gzip = [0x1f, 0x8b, 0x08, 0x03, 0x00, 0x00, 0x00, 0x00,
              0x02, 0x03, 0x91, 0x1e, 0x03, 0x00, 0x00, 0x00,
              0x00, 0x00, 0x00, 0x00, 0x00, 0x00]

@testset "Empty messages" begin
    @test inflate(empty_deflate) == UInt8[]
    @test inflate_zlib(empty_zlib) == UInt8[]
    @test inflate_gzip(empty_gzip) == UInt8[]
end

@testset "Deflate corruption" begin
    d1 = UInt8[0x01, 0x00, 0x00, 0x00, 0x00] # corrupted compression mode 0 data
    d2 = UInt8[0x07]                         # invalid compression mode 3
    d3 = UInt8[0xed, 0x1c, 0xed, 0x72,
               0xdb, 0x48, 0xf2, 0x3f]       # incomplete code table
    for d in [d1, d2, d3]
        @test_throws ErrorException inflate(d)
        @test_throws ErrorException read(InflateStream(IOBuffer(d)))
    end
end

@testset "Zlib header corruption" begin
    z1 = copy(empty_zlib); z1[1] = 0x77   # unsupported compression method
    z2 = copy(empty_zlib); z2[1] = 0x88   # invalid LZ77 window size
    z3 = copy(empty_zlib); z3[2] = 0xbc   # preset dictionary not supported
    z4 = copy(empty_zlib); z4[2] = 0x9d   # header checksum error
    for z in [z1, z2, z3, z4]
        @test_throws ErrorException inflate_zlib(z)
        @test_throws ErrorException InflateZlibStream(IOBuffer(z))
    end
end

@testset "Zlib trailer corruption" begin
    z5 = copy(empty_zlib); z5[5] = 0x01   # adler checksum error
    z6 = copy(empty_zlib); z6[6] = 0x01   # adler checksum error
    z7 = copy(empty_zlib); z7[7] = 0x01   # adler checksum error
    z8 = copy(empty_zlib); z8[8] = 0x02   # adler checksum error
    for z in [z5, z6, z7, z8]
        @test_throws ErrorException inflate_zlib(z)
        @test_throws ErrorException read(InflateZlibStream(IOBuffer(z)))
        @test isempty(inflate_zlib(z, ignore_checksum = true))
        @test eof(InflateZlibStream(IOBuffer(z), ignore_checksum = true))
    end
end

@testset "Gzip header corruption" begin
    g1 = copy(empty_gzip); g1[1] = 0x1e   # not gzipped data
    g2 = copy(empty_gzip); g2[2] = 0x8a   # not gzipped data
    g3 = copy(empty_gzip); g3[3] = 0x07   # unsupported compression method
    g4 = copy(empty_gzip); g4[4] = 0x23   # reserved FLG bit set
    g5 = copy(empty_gzip); g5[11] = 0x92  # header crc error
    g6 = copy(empty_gzip); g6[12] = 0x1f  # header crc error
    for g in [g1, g2, g3, g4, g5, g6]
        @test_throws ErrorException inflate_gzip(g)
        @test_throws ErrorException InflateGzipStream(IOBuffer(g))
    end
    for g in [g5, g6]
        @test isempty(inflate_gzip(g, ignore_checksum = true))
        @test eof(InflateGzipStream(IOBuffer(g), ignore_checksum = true))
    end
end

@testset "Gzip trailer corruption" begin
    g7 = copy(empty_gzip); g7[15] = 0x01   # crc check failed
    g8 = copy(empty_gzip); g8[16] = 0x01   # crc check failed
    g9 = copy(empty_gzip); g9[17] = 0x01   # crc check failed
    g10 = copy(empty_gzip); g10[18] = 0x01   # crc check failed
    g11 = copy(empty_gzip); g11[19] = 0x01   # length check failed
    g12 = copy(empty_gzip); g12[20] = 0x01   # length check failed
    g13 = copy(empty_gzip); g13[21] = 0x01   # length check failed
    g14 = copy(empty_gzip); g14[22] = 0x01   # length check failed
    for g in [g7, g8, g9, g10, g11, g12, g13, g14]
        @test_throws ErrorException inflate_gzip(g)
        @test_throws ErrorException read(InflateGzipStream(IOBuffer(g)))
    end
    for g in [g7, g8, g9, g10]
        @test isempty(inflate_gzip(g, ignore_checksum = true))
        @test eof(InflateGzipStream(IOBuffer(g), ignore_checksum = true))
    end
end

@testset "Reading past end of file" begin
    s1 = InflateStream(IOBuffer(empty_deflate))
    s2 = InflateZlibStream(IOBuffer(empty_zlib))
    s3 = InflateGzipStream(IOBuffer(empty_gzip))
    for s in [s1, s2, s3]
        @test eof(s)
        @test_throws EOFError read(s, UInt8)
    end
end


text:
function ode_determine_initdt{tType,uType,F}(u0,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,true,F},order)
  f = prob.f
  f₀ = similar(u0./t); f₁ = similar(u0./t); u₁ = similar(u0)
  sk = abstol+abs.(u0)*reltol
  d₀ = internalnorm(u0./sk)
  f(t,u0,f₀)
  d₁ = internalnorm((f₀./sk*tType(1))/tType(1))
  T0 = typeof(d₀)
  T1 = typeof(d₁)
  if d₀ < T0(1//10^(5)) || d₁ < T1(1//10^(5))
    dt₀ = tType(1//10^(6))
  else
    dt₀ = tType((d₀/d₁)/100)
  end
  dt₀ = min(dt₀,tdir*dtmax)
  @inbounds for i in eachindex(u0)
     u₁[i] = u0[i] + tdir*dt₀*f₀[i]
  end
  f(t+tdir*dt₀,u₁,f₁)
  tmp = (f₁.-f₀)./(abstol+abs.(u0)*reltol)*tType(1)
  d₂ = internalnorm(tmp)/dt₀
  # Hairer has d₂ = sqrt(sum(abs2,tmp))/dt₀, note the lack of norm correction
  if max(d₁,d₂)<=T1(1//10^(15))
    dt₁ = max(tType(1//10^(6)),dt₀*1//10^(3))
  else
    dt₁ = tType(10.0^(-(2+log10(max(d₁,d₂)/T1(1)))/(order)))
  end
  dt = tdir*min(100dt₀,dt₁)
end

function ode_determine_initdt{uType,tType,F}(u0::uType,t,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,false,F},order)
  f = prob.f
  sk = abstol+abs.(u0)*reltol
  d₀ = internalnorm(u0./sk)
  f₀ = f(t,u0)
  d₁ = internalnorm(f₀./sk)
  T0 = typeof(d₀)
  T1 = typeof(d₁)
  if d₀ < T0(1//10^(5)) || d₁ < T1(1//10^(5))
    dt₀ = tType(1//10^(6))
  else
    dt₀ = tType((d₀/d₁)/100)
  end
  dt₀ = min(dt₀,tdir*dtmax)
  u₁ = u0 + tdir*dt₀*f₀
  f₁ = f(t+tdir*dt₀,u₁)
  d₂ = internalnorm((f₁-f₀)./(abstol+abs.(u0)*reltol))/dt₀*tType(1)
  if max(d₁,d₂) <= T1(1//10^(15))
    dt₁ = max(tType(1//10^(6)),dt₀*1//10^(3))
  else
    dt₁ = tType(10.0^(-(2+log10(max(d₁,d₂)/T1(1)))/(order)))
  end
  dt = tdir*min(100dt₀,dt₁)
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ b58c89f1-518d-4008-81a5-1557b2852711
include("structs.jl")

# ╔═╡ 04bbf20d-ebb9-4acc-b821-27ee3d5b3b0a
md"""
# Boundary Conditions Documentation

We want to create a general way to update the BC for readability and in order to minimize errors.

``
\left. \left( af(x, y) + b\frac{\partial f}{\partial n} \right) \right\vert_{(x, y) \in \partial \Omega} = g
``

#### Discretization

``
a \frac{f(\vec{r} - \vec{h})  + f(\vec{r} + \vec{h})}{2} +
b \frac{f(\vec{r} + \vec{h}) - f(\vec{r} - \vec{h})}{2h} = g
``

##### Example: Right Boundary

``
\vec{h} = (\frac{\Delta x}{2}, 0)
``

``
a \frac{f_{nx} + f_{nx - 1}}{2} + b \frac{f_{nx} - f_{nx - 1}}{\Delta x} = g
``

``
f_{nx} = \frac{2 \Delta x g + f_{nx - 1} (2b - a \Delta x)}{a \Delta x + 2b}
``

##### General Formula

``
f_{ghost} = \frac{2 \Delta h g + f_{inside} (2b - a \Delta h)}{a \Delta h + 2b}
``

#### Implementation

We keep a and b in an unmutable struct because they are set as a parameter in the beginning of the simulation, and they don't change. Therefore, they are kept with the other simulation parameters

g, on the other hand, will get updated on every step, so it seems logical to keep together with all the other system parameters that change over time. We will denote it as f_RHS, since it is then easy to infer that we are refering the right hand side of the BC equation.

The struct BoundaryCondtions holds 2 3D arrays, a and b.

Each of this holds 4 subarrays for each component, one for each border: left, right, top, bottom, in the form of a 2D array.

h[:][1][k] -> left

h[:][2][k] -> right

h[:][3][k] -> top

h[:][4][k] -> bottom

Pressure is the same for both components, but for consistency we will also use the same structure, making use of the first index k = 1 only.

"""

# ╔═╡ fc2c2996-60f6-476f-8cdf-bdde873628c4
NamedTuple{fieldnames(BoundaryCondition)}(fieldtypes(BoundaryCondition))

# ╔═╡ 550f6393-9a8a-4ce1-9f7f-dd5d0012fac1
md"""
The structure Parameters stores bc::BoundaryCondtions
"""

# ╔═╡ 47a99444-4329-4236-ad36-6c0e9491c250
NamedTuple{fieldnames(BoundaryConditions)}(fieldtypes(BoundaryConditions))

# ╔═╡ 72246b26-fa0f-49c2-9f72-9ee300d34d41
md"""
The struture Systems stores bc_RHS::BoundaryRHS
"""

# ╔═╡ c121c421-38b0-49e2-893a-f2645d9b2071
NamedTuple{fieldnames(BoundaryRHS)}(fieldtypes(BoundaryRHS))

# ╔═╡ bf125994-a82d-409e-a195-8b77d5c2cd3e

# ╔═╡ 923b99aa-52dc-4068-b888-b5cc9d8c516a
md"""
#### Algorithm

1. Enforce the BC: ``\frac{\partial P}{\partial n} = 0`` on the walls, ``P_{in}`` and ``P_{out}`` on inlet and outlet, respectively.

2. Obtain ``\hat\rho_{\alpha, \partial \Omega} = \hat\rho_{\alpha, \partial \Omega}(T, P_{\partial \Omega}) `` from their respective equations of state. ``\alpha = 1, 2`` - component. On the inlet and outlet, the initial pressure is given as ``P_{in}`` and ``P_{out}``, respectively. On the walls, we derived it in the previous step.

3. Obtain the boundary saturation.

On the inlet:

``
s_{in} = s_{in}(\psi_{in},
\hat \rho_{1, in}, \hat \rho_{2, in})
= \left(
        \frac{\hat \rho_{1, in} M_2(1 - \psi_{in}) }{\hat \rho_{2, in} M_1 \psi_{in}}
    + 1 \right)^{-1}
``

``
\left. \left( 1 \cdot s +
0 \cdot \frac{\partial s}{\partial n} \right) \right\vert_{inlet}
= s_{in}
``

On the outlet and walls:

``
\frac{\partial s}{\partial n} = 0
``

``
\left. \left( 0 \cdot s +
1 \cdot \frac{\partial s}{\partial n} \right) \right\vert_{outlet, walls}
= 0
``

4. Find ``\rho_{\alpha, \partial \Omega} = \hat\rho_{\alpha, \partial \Omega} s_\alpha``.

``
\left. \left( 1 \cdot \rho_\alpha +
0 \cdot \frac{\partial \rho_\alpha}{\partial n} \right) \right\vert_{inlet}
= \rho_{\alpha, \partial \Omega}
``

5. Update densities at boundary.

``
\rho_{k, \partial \Omega} = \hat \rho_{k, \partial \Omega} \cdot
s_{k, \partial \Omega}
``

6. Find ``P_{\partial \Omega}`` on the inlet and outlet from the densities at the boundary (using Newton-Raphson).

7. Enforce BC for velocity using Darcy's Law on inlet and outlet

``
    v_{in} = -\frac{K}{\mu} f_\alpha (s_{in})
    \frac{P_{i, 2} - P_{i, 1}^n}{\Delta y}
``

``
    v_{out} = -\frac{K}{\mu} f_\alpha (s_{out})
    \frac{P_{i, ny} - P_{i, ny - 1}^n}{\Delta y}
``

``
\left. \left( 1 \cdot v +
0 \cdot \frac{\partial v}{\partial n} \right) \right\vert_{in, out}
= v_{in, out}
``

"""

# ╔═╡ 2ef6d608-3aef-4647-a523-47fab6ff1154

# ╔═╡ f732b89b-5361-4dc5-9051-62fa8ae698cb

# ╔═╡ c427e959-cc5e-4eeb-b50c-8c06aec042e8

# ╔═╡ 6f6caa2a-f2ed-4fa0-9964-0af022384749

# ╔═╡ ce3afddf-9b25-4468-8777-b77a6898da5f

# ╔═╡ Cell order:
# ╟─04bbf20d-ebb9-4acc-b821-27ee3d5b3b0a
# ╠═b58c89f1-518d-4008-81a5-1557b2852711
# ╠═fc2c2996-60f6-476f-8cdf-bdde873628c4
# ╟─550f6393-9a8a-4ce1-9f7f-dd5d0012fac1
# ╠═47a99444-4329-4236-ad36-6c0e9491c250
# ╟─72246b26-fa0f-49c2-9f72-9ee300d34d41
# ╠═c121c421-38b0-49e2-893a-f2645d9b2071
# ╠═bf125994-a82d-409e-a195-8b77d5c2cd3e
# ╟─923b99aa-52dc-4068-b888-b5cc9d8c516a
# ╠═2ef6d608-3aef-4647-a523-47fab6ff1154
# ╠═f732b89b-5361-4dc5-9051-62fa8ae698cb
# ╠═c427e959-cc5e-4eeb-b50c-8c06aec042e8
# ╠═6f6caa2a-f2ed-4fa0-9964-0af022384749
# ╠═ce3afddf-9b25-4468-8777-b77a6898da5f

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
@testset "turing compat" begin
    k = SqExponentialKernel()
    y = randn(3)
    X = randn(3, 1)
    x = [rand(1) for _ in 1:3]
    @model function GPRegression(y, X)
        # Priors.
        α ~ LogNormal(0.0, 0.1)
        ρ ~ LogNormal(0.0, 1.0)
        σ² ~ LogNormal(0.0, 1.0)

        # Realized covariance function
        kernel = α * transform(SqExponentialKernel(), 1 / ρ)
        f = GP(kernel)

        # Sampling Distribution.
        y ~ f(X, σ²)
    end
    # Test for matrices
    m = GPRegression(y, RowVecs(X))
    @test length(sample(m, HMC(0.5, 1), 5)) == 5
    # Test for vectors of vector
    m = GPRegression(y, x)
    @test length(sample(m, HMC(0.5, 1), 5)) == 5
end