text:
#function __init__()
#    global chr_rngs = load_chr_rngs()
#    global rng_exon = load_rng_exon()
#end

@doc """ return gene_name, exton_number
index1: chr=>ranges
(1,5) (7,10) (12,16) (18,20)
1 7 12 18
5 10 16 20
odd is st, even is ed.
from odd, find the two numbers surround Q
from even, find the two numbers surround Q
index2: ranges=>exon
""" ->
function query_exon(chr::ASCIIString,pos::ASCIIString,direct::ASCIIString)
    global chr_rngs = load_chr_rngs()
    global rng_exon = load_rng_exon()
    pos = parse(UInt64,pos)
    query_exon(chr,pos,direct)
end
function query_exon(chr::ASCIIString, pos::UInt64, direct::ASCIIString)
    #load index
    rngs_target = chr_rngs[chr]
    sts_sorted = map(x->x[1],rngs_target)
    eds_sorted = map(x->x[2],rngs_target)

    #need to check
    function ret(idx)
        rng = map(x->string(Int64(x)), rngs_target[idx])
        map(x->convert(ASCIIString,x), rng_exon[chr,rng])
    end
    num_rng = length(rngs_target)
    st_idx = searchsortedfirst(sts_sorted, pos)
    ed_idx = searchsortedfirst(eds_sorted, pos)
    if st_idx == 1
        return ret(1)
    end
    if st_idx>num_rng || ed_idx>num_rng
        return ret(st_idx-1)
    end
    st1,st2 = sts_sorted[st_idx-1:st_idx] # bugs here
    if st2 == pos
        return ret(st_idx)
    end
    if ed_idx == 1
        return ret(1)
    end
    ed1,ed2 = eds_sorted[ed_idx-1:ed_idx]
    if ed2 == pos
        return ret(ed_idx)
    end
    #@show pos,st1,st2,ed1,ed2
    if st_idx != ed_idx
        return ret(st_idx-1)
    elseif st_idx == ed_idx
        if direct == "-"
            return ret(st_idx)
        elseif direct == "+"
            return ret(st_idx-1)
        end
    end

    #=
    @assert pos <= st2
    @assert pos <= ed2
    if st1 < pos < ed1
        return rngs_target[st_idx-1]
    elseif st2 < pos < ed2
        error("pos cannot be here")
        return
    elseif ed1<pos<st2
        if direct == "+"
            return rngs_target[st_idx]
        elseif direct == "-"
            return rngs_target[st_idx-1]
        else
            error("direct not detected!")
        end
    end
    =#
    error("unexpected")
end

function load_chr_rngs()
    chr_rngs = build_index_exon()
    #chr_rngs = load(chr_rngs_exon_fl, "chr_rngs_exon_dict")
    chr_rngs
end

function load_rng_exon()
    #load(rng_exon_dict_fl, "rng_exon_dict")
    range_exon_index(read_gtf_exon())
end


text:
"""
TLeakyIntegrator: Leaky integration for the Categorical task (section 2.3.2)
"""
struct TLeakyIntegrator
    ns::Int
    etaleak::Float64
    Ns::Array{Float64, 1}
    Ns1::Array{Float64, 1}
    Ps1::Array{Float64, 1}
end
function TLeakyIntegrator(; ns = 10, etaleak = .9)
    Ns = zeros(1)
    Ns1 = zeros(ns)
    Ps1 = zeros(ns)
    TLeakyIntegrator(ns, etaleak, Ns, Ns1, Ps1)
end
export TLeakyIntegrator
""" X[t] = etaleak * X[t-1] + etaleak * I[.] """
function update!(learnerT::TLeakyIntegrator, s1)
    learnerT.Ns[1] *= learnerT.etaleak # Discount
    learnerT.Ns[1] += learnerT.etaleak # Increase
    @. learnerT.Ns1 = learnerT.Ns1 .* learnerT.etaleak # Discount all
    learnerT.Ns1[s1] += learnerT.etaleak # Increase observed

    computePs1!(learnerT)
end
export update!

function computePs1!(learnerT::TLeakyIntegrator)
    @. learnerT.Ps1 = learnerT.Ns1 ./ learnerT.Ns[1]
end

ext: jl


text:
using MieScatter
using Test

@testset "Cross sections" begin

    # These tests are made compared to the FORTRAN77 Mie code
    # written by Karri Muinonen to the precision of the output
    # of that code.

    # Test computation of scattering, extinction and absorbsion
    # cross sections.

    # Size parameter 2, refractive index 1.33
    S, Qsca, Qext, Qback = compute_mie(2.0, 1.33, 181)
    Qabs = Qext-Qsca
    @test isapprox(Qsca, 0.71295 ; atol = 6)
    @test isapprox(Qext, 0.71295 ; atol = 6)
    @test isapprox(Qabs, 0.0 ; atol = 10)

    # Size parameter 2, refractive index 1.6
    S, Qsca, Qext, Qback = compute_mie(2.0, 1.6, 181)
    Qabs = Qext-Qsca
    @test isapprox(Qsca, 2.4834 ; atol = 5)
    @test isapprox(Qext, 2.4834 ; atol = 5)
    @test isapprox(Qabs, 0.0 ; atol = 10)

    # Size parameter 20, refractive index 1.33
    S, Qsca, Qext, Qback = compute_mie(20.0, 1.33, 181)
    Qabs = Qext-Qsca
    @test isapprox(Qsca, 0.21401E+01 ; atol = 6)
    @test isapprox(Qext, 0.21401E+01 ; atol = 6)
    @test isapprox(Qabs, 0.66613E-14 ; atol = 6)

    # Size parameter 20, refractive index 1.6
    S, Qsca, Qext, Qback = compute_mie(20.0, 1.6, 181)
    Qabs = Qext-Qsca
    @test isapprox(Qsca, 0.26115E+01 ; atol = 6)
    @test isapprox(Qext, 0.26115E+01 ; atol = 6)
    @test isapprox(Qabs, 0.39968E-14 ; atol = 6)

end # cross section testset

@testset "Return types and NaN" begin
    S, Qsca, Qext, Qback = compute_mie(20.0, 1.6, [pi/2])
    @test typeof(S) == Array{Float64,2}
    @test size(S) == (1,4)
    @test typeof(Qsca) == Float64
    @test typeof(Qext) == Float64
    @test typeof(Qback) == Float64
    @test isnan(Qext)
    @test isnan(Qback)
    @test !isnan(Qsca)
end


text:
Base.(:(==))(o1::GitObject, o2::GitObject) = isequal(Oid(o1), Oid(o2))
Base.isequal(o1::GitObject, o2::GitObject) = isequal(Oid(o1), Oid(o2))
Base.isless(o1::GitObject, o2::GitObject)  = isless(Oid(o1), Oid(o2))

Base.hash(o::GitObject) = hash(hex(o))
Base.cmp(o1::GitObject, o2::GitObject) = cmp(Oid(o1), Oid(o2))

Oid(o::GitObject) = Oid(ccall((:git_object_id, libgit2), Ptr{Oid}, (Ptr{Void},), o))

Base.hex(o::GitObject) = begin
    oid_ptr = ccall((:git_object_id, libgit2), Ptr{Uint8}, (Ptr{Void},), o)
    hex_buff = Array(Uint8, OID_HEXSZ)
    @check ccall((:git_oid_fmt, libgit2), Cint, (Ptr{Uint8}, Ptr{Uint8}), hex_buff, oid_ptr)
    return bytestring(hex_buff)
end

function raw(o::GitObject)
    repo_ptr = ccall((:git_object_owner, libgit2), Ptr{Void}, (Ptr{Void},), o)
    oid_ptr  = ccall((:git_object_id, libgit2), Ptr{Oid}, (Ptr{Void},), o)
    odb_ptr  = Ptr{Void}[0]
    @check ccall((:git_repository_odb, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), odb_ptr, repo_ptr)
    obj_ptr  = Ptr{Void}[0]
    err = ccall((:git_odb_read, libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Oid}), obj_ptr, odb_ptr[1], oid_ptr)
    ccall((:git_odb_free, libgit2), Void, (Ptr{Void},), odb_ptr[1])
    if err != GitErrorConst.GIT_OK
        throw(GitError(err))
    end
    return OdbObject(obj_ptr[1])
end

function gitobj_from_ptr(ptr::Ptr{Void})
    obj_type = ccall((:git_object_type, libgit2), Cint, (Ptr{Void},), ptr)
    T = gitobj_const_type(obj_type)
    @assert T <: GitObject
    return T(ptr)
end

function gitobj_const_type(typ::Cint)
    typ == GitConst.OBJ_BLOB   && return GitBlob
    typ == GitConst.OBJ_TREE   && return GitTree
    typ == GitConst.OBJ_COMMIT && return GitCommit
    typ == GitConst.OBJ_TAG    && return GitTag
    typ == GitConst.OBJ_ANY    && return GitAnyObject
    error("Unknown git type const: $obj_type")
end


text:
module TestGradient
using Test
using LinearAlgebra
using BlockFactorizations

using CovarianceFunctions
using CovarianceFunctions: EQ, RQ, Dot, ExponentialDot, NN, Matern, MaternP,
        Lengthscale, input_trait, GradientKernel, ValueGradientKernel, GradientKernelElement,
        DerivativeKernel, ValueDerivativeKernel, DerivativeKernelElement, Cosine,
        Woodbury, LazyMatrixProduct, Constant

const AbstractMatOrFac = Union{AbstractMatrix, Factorization}

@testset "gradient kernels" begin
    # nd test
    n = 2
    d = 5
    X = randn(d, n) / sqrt(d)
    ε = 1e4eps()
    @testset "GradientKernel" begin
        kernels = [MaternP(3), Dot()^3, Cosine(randn(d)), NN()]
        # kernels = [EQ(), RQ(1.), MaternP(2), Matern(2.7), Dot()^3, ExponentialDot(), NN(), Cosine(randn(d))]
        # kernels = [EQ(), Lengthscale(EQ(), 0.1)] # at this time, test doesn't work because fallback is incorrect
        a = randn(d*n)
        b = randn(d*n)
        for k in kernels
            G = GradientKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            @test K isa BlockFactorization{<:Real}

            MK = Matrix(K)
            @test maximum(abs, MK - MK') < ε
            MK = Symmetric(MK)
            @test all(≥(-1e-12), eigvals(MK)) # positive semidefinite
            @test size(MK) == (d*n, d*n)

            # create anonymous wrapper around kernel to trigger generic fallback
            k2 = (x, y)->k(x, y)
            G2 = GradientKernel(k2)

            # this takes centuries to precompile, maybe because of nested differentiation?
            # not really important because it is a fallback anyway
            K2 = CovarianceFunctions.gramian(G2, X) # NOTE: this is the bottleneck in test suite
            MK2 = Matrix(K2)
            @test MK ≈ MK2 # compare generic and specialized evaluation of gradient kernel

            # testing matrix multiply
            Kab = deepcopy(b)
            α, β = randn(2) # trying out two-argument mul!
            mul!(Kab, K, a, α, β) # saves some memory, Woodbury still allocates a lot
            MKab = @. α * $*(MK, a) + β * b
            @test Kab ≈ MKab
        end

        # testing matrix solve
        for k in kernels # IDEA: better conditioned NN kernel with bias
            G = GradientKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            a = randn(n*d)
            Ka = K * a
            as = K \ Ka
            @test norm(K*as-Ka) / norm(Ka) < 1e-6
        end

        # testing non-lazy data-sparse representation
        x, y = X[:, 1], X[:, 2]
        a = randn(length(x))
        G = GradientKernelElement(EQ(), x, y)
        W = Woodbury(G)
        @test W*a ≈ G*a

        G = GradientKernelElement(Dot()^3, x, y)
        W = Woodbury(G)
        @test W*a ≈ G*a

        G = GradientKernelElement(Cosine(randn(d)), x, y)
        W = LazyMatrixProduct(G)
        @test W*a ≈ G*a

        # testing constant kernel
        c = Constant(1)
        g = GradientKernel(c)
        @test g(x, y) ≈ zeros(d, d)
    end

    # TODO:
    @testset "ValueGradientKernel" begin
        kernels = [MaternP(3), Dot()^3, Cosine(randn(d))]
        a = randn((d+1)*n)
        b = randn((d+1)*n)
        for k in kernels
            G = ValueGradientKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            @test K isa BlockFactorization{<:Real}

            MK = Matrix(K)
            @test maximum(abs, MK - MK') < ε
            MK = Symmetric(MK)
            if !all(≥(-1e-10), eigvals(MK))
                println("eigen-snafoo")
                display(eigvals(MK))
            end
            @test all(≥(-1e-10), eigvals(MK)) # positive semidefinite
            @test size(MK) == ((d+1)*n, (d+1)*n)

            # create anonymous wrapper around kernel to trigger generic fallback
            k2 = (x, y)->k(x, y)
            G2 = ValueGradientKernel(k2)
            K2 = CovarianceFunctions.gramian(G2, X) # NOTE: compilation here is super slow (not important for real applicaiton since it is fallback)
            MK2 = Matrix(K2)
            @test MK ≈ MK2 # compare generic and specialized evaluation of gradient kernel

            # instead, compare to already tested GradientKernel
            # K12 = K.A[1, 2]
            # @test K12 isa DerivativeKernelElement
            # @test K12.value_value ≈ k(X[:, 1], X[:, 2])
            # @test Matrix(K12.gradient_gradient) ≈ Matrix(GradientKernel(k)(X[:, 1], X[:, 2]))

            # testing matrix multiply
            Kab = deepcopy(b)
            α, β = randn(2) # trying out two-argument mul
            mul!(Kab, K, a, α, β) # saves some memory, Woodbury still allocates a lot
            MKab = @. α * $*(MK, a) + β * b
            @test Kab ≈ MKab
        end

        # testing matrix solve
        # println("div")
        for k in [EQ(), RQ(1.)] # TODO: better conditioned NN kernel with bias
            G = ValueGradientKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            a = randn(n*(d+1))
            Ka = K*a
            as = K\Ka
            @test norm(K*as-Ka) / norm(Ka) < 1e-6
        end
    end
    # matrix multiply
    # y = zeros(n*d)
    # mul!(y, K, x)

    # 1d special case
    @testset "ValueDerivativeKernel" begin
        k = CovarianceFunctions.EQ()
        G = ValueDerivativeKernel(k)
        n = 3
        x = randn(n)
        K = CovarianceFunctions.gramian(G, x)
        MK = Matrix(K)
        @test issymmetric(MK)
        @test all(≥(-1e-10), eigvals(MK)) # positive semidefinite
        @test size(K) == (2n, 2n)
        @test size(MK) == (2n, 2n)
    end
end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)
module Baz
import Foo, Qux
this = "Baz"
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using POMDPs, Crux, Flux, POMDPGym, Random, Distributions, POMDPPolicies, Plots
Crux.set_function("isfailure", POMDPGym.isfailure)

# Define the disturbance distribution based on a normal distribution
px =MvNormal2([0f0], [0.5f0])

# Construct the MDP
mdp = AdditiveAdversarialMDP(InvertedPendulumMDP(λcost=0, include_time_in_state=true), px)
S = state_space(mdp)
N = 100000

# construct the model
QSA() = ContinuousNetwork(Chain(Dense(4, 64, relu), Dense(64, 64, relu), Dense(64, 1)))
Pf() = ContinuousNetwork(Chain(Dense(4, 64, relu), Dense(64, 64, relu), Dense(64, 1, sigmoid)))
A() = ContinuousNetwork(Chain(Dense(3, 64, relu), Dense(64, 64, relu), Dense(64, 1, tanh), x -> 2f0 * x), 1)

function G()
    base = Chain(Dense(3, 64, relu), Dense(64, 64, relu))
    mu = ContinuousNetwork(Chain(base..., Dense(64, 1)))
    logΣ = ContinuousNetwork(Chain(base..., Dense(64, 1)), 1)
    GaussianPolicy2(mu, logΣ, true)
end

Protag() = ActorCritic(A(), DoubleNetwork(QSA(), QSA()))
Antag() = ActorCritic(G(), Pf())

AdvPol(p = Protag()) = AdversarialPolicy(p, Antag())

𝒮_td3 = TD3(;π=ActorCritic(A(), DoubleNetwork(QSA(), QSA())), S=S, N=50000, buffer_size=Int(1e5), buffer_init=1000)
π_td3 = solve(𝒮_td3, mdp)

# solve with IS
𝒮_isarl = ISARL_Continuous(π=AdvPol(),
                           S=S,
                           N=100000,
                           px=px,
                           buffer_size=Int(1e5),
                           buffer_init=1000,
                           desired_AP_ratio=1)
π_isarl = solve(𝒮_isarl, mdp)

# check the number of failures in the buffers (compared to successs)
sum(𝒮_isarl.buffer[:fail])
sum(𝒮_isarl.buffer[:done])

# Plot the distribution of disturbances
histogram(𝒮_isarl.buffer[:x][:])

# Plot the value function in 2D
v(θ, t) = value(antagonist(𝒮_isarl.π), [θ, 0, t, 0])[1]
heatmap(deg2rad(-25):0.01:deg2rad(25), 0:0.1:5, v)

# Show the distribution of data in the replay buffer
scatter(𝒮_isarl.buffer[:s][1, :], 𝒮_isarl.buffer[:s][3, :], marker_z = 𝒮_isarl.buffer[:done][1,:], xlabel="θ", ylabel="t", alpha=0.5)
vline!([deg2rad(20), deg2rad(-20)])

plot(-3:0.1:3, [exp.(logpdf(antagonist(𝒮_isarl.π).A, [0, 0, 0], [x]))[1] for x=-3:0.1:3], label="Px")

antagonist(𝒮_isarl.π).A.logΣ([0, 0, 0])
, [1])[1]

# Solve with DQN
𝒮_dqn = DQN(π=QS(as), S=S, N=N)
π_dqn = solve(𝒮_dqn, mdp)

pfail_dqn = Crux.failure(Sampler(mdp, π_dqn, S=S, max_steps=100), Neps=Int(1e3), threshold=100)
println("DQN Failure rate: ", pfail_dqn)

# solve with RARL
𝒮_rarl = RARL(π=AdvPol(), S=S, N=N)
π_rarl = solve(𝒮_rarl, mdp)

pfail_rarl = Crux.failure(Sampler(mdp, protagonist(π_rarl), S=S, max_steps=100), Neps=Int(1e5), threshold=100)
println("RARL Failure rate: ", pfail_rarl)

pfail_isarl = Crux.failure(Sampler(mdp, protagonist(π_isarl), S=S, max_steps=100), Neps=Int(1e5), threshold=100)
println("IS Failure rate: ", pfail_isarl)

pol = AdvPol()

𝒮_isarl.buffer

pol = AdversarialPolicy(π_dqn, Pf(xs), ϵGreedyPolicy(Crux.LinearDecaySchedule(1., 0.1, floor(Int, N/2)), xs))
𝒮_isarl = ISARL_Discrete(π=pol, S=S, N=N, xlogprobs=xlogprobs)
π_isarl = solve(𝒮_isarl, mdp)


text:
import Base.diag

function diag(a::Number)
	return a
end

function diagmmult(A,l)
    (n,m)=size(A)
    Anew=zeros(A)
    for i in 1:m
        Anew[:,i]=A[:,i]*l[i]
    end

    return Anew
end

function bdot(a,b)
    (n,m)=size(a)
    out=zeros(ComplexF64,m)
    for i in 1:m
        out[i]=dot(a[:,i],b[:,i])
    end

    return out
end

function zbicgstabBlock(A,b,x0,maxit,eps)
	#matrix multiply:

	(n,m)=size(x0)

        r=b-A*x0
	s=copy(r)
	y0=copy(r)

	x=copy(x0)

	#inner product:
	#delta=(y0'*r)[1]

	#delta=diag(y0'*r)
	delta=bdot(y0,r)

	#phi=(y0'*(A*s))[1]/delta

	for i in 1:maxit
        #println("    BCG $i")
		#matrix multiply:
		As=A*s
		#inner product
		#phi=(y0'*(As))[1]/delta
		#phi=diag(y0'*As) ./delta
		phi=bdot(y0,As) ./delta

		#omega=1/phi
		omega=1 ./phi
		w=r-As*spdiagm( 0 => omega)
		#w=r-diagmmult(As,omega)

		#matrix multiply:
		Aw=A*w
		#inner product:
		#chi=((Aw)'*w)[1]/norm(Aw)^2
		#chi=diag(Aw'*w) ./diag(Aw'*Aw)
        chi=bdot(Aw,w) ./bdot(Aw,Aw)

		#r=w-chi*Aw
		r=w-Aw*spdiagm( 0 => chi)
        #r=w-diagmmult(Aw,chi)

		rs=sqrt.(bdot(r,r))
		#println("   $i     $(maximum(abs.(rs)))   $(norm(Aw))")
		#if(norm(r)<eps)
		if(maximum(abs.(rs))<eps)
			#println("r finish: $(norm(r))")
			#println("     bicgstab its $i")
			return x
		end

		#x=x+s*omega+w*chi
		x=x+s*spdiagm( 0 =>  omega)+w*spdiagm( 0 =>  chi)
		#x=x+diagmmult(s,omega)+diagmmult(w,chi)

		deltaold=delta
		#inner product:
		#delta=(y0'*r)[1]
		delta=bdot(y0,r)#diag(y0'*r)
		#psi=-1.0*omega*delta/(deltaold*chi)
		psi=-1.0*omega.*delta ./(deltaold.*chi)
		s=r-(s-(As)*spdiagm( 0 =>  chi))*spdiagm( 0 =>  psi)
		#s=r-(s-diagmmult(As,chi.*psi))
		#print("     $(norm(s))\n")
		ss=sqrt.(bdot(s,s))
		if(maximum(abs.(ss))<eps)
		#if(norm(s)<eps)
			#println("s finish: $(norm(s))")
			#println("     bicgstab its $i")
			return x
		end

	end

	#println("          no converge: $(norm(r))")
	return x
end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
using Test
using SineFit

function wp_roughly_equal(a, b)
    ok = true
    for f in [:amplitude, :frequency, :phaseshft, :vertoffst]
        ok &= round(getfield(a, f), digits=3) == round(getfield(b, f), digits=3)
    end
    ok
end

function test_basic_estimation()
    step = 0.01
    xs = collect(range(0, 10, step=step))
    ys = 3 .* Base.sin.(2pi * 1 .* xs .+ 2.2) .- 5
    wp = SineFit.estimate_wave_parameters(xs, ys, sample_rate=round(Int32, 1/step))

    println("test_basic_estimation: ", wp)
    wp
end

@test wp_roughly_equal(test_basic_estimation(), SineFit.WaveFitParams(3.0, 6.276908398780805, 2.230908739161565, -5.0))

function test_basic_fit()
    step = 1.
    xs = collect(range(0, 1000, step=step))
    ys = .3 .* Base.sin.(2pi * .1 .* xs .+ 1.1) .- 5
    wp = SineFit.calculate_wave_shape(xs, ys, sample_rate=round(Int32, 1/step))

    println("test_basic_fit: ", wp)
    wp
end

@test wp_roughly_equal(test_basic_fit(), SineFit.WaveFitParams(0.299999999999996, 0.6283185307179554, 1.1, -5.0))


text:
struct TeeLogger{T<:NTuple{<:Any, AbstractLogger}} <: AbstractLogger
    loggers::T
end

"""
    TeeLogger(loggers...)

Send the same log message to all the loggers.

To include the current logger do:
`TeeLogger(current_logger(), loggers...)`
to include the global logger, do:
`TeeLogger(global_logger(), loggers...)`
"""
function TeeLogger(loggers::Vararg{AbstractLogger})
    return TeeLogger(loggers)
end

function handle_message(demux::TeeLogger, args...; kwargs...)
    for logger in demux.loggers
        if comp_handle_message_check(logger, args...; kwargs...)
            handle_message(logger, args...; kwargs...)
        end
    end
end

function shouldlog(demux::TeeLogger, args...)
    any(comp_shouldlog(logger, args...) for logger in demux.loggers)
end

function min_enabled_level(demux::TeeLogger)
    minimum(min_enabled_level(logger) for logger in demux.loggers)
end

function catch_exceptions(demux::TeeLogger)
    any(catch_exceptions(logger) for logger in demux.loggers)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 864b72c0-3e11-11eb-2cd6-9951cdcb78c3
begin#imports
	using Flux:gradient,params
	using Flux.Optimise:Descent,train!,@epochs
	using Statistics:mean
end

# ╔═╡ 588bf94e-3e17-11eb-216c-090eb107222e
using Flux:Chain,Dense,softmax

# ╔═╡ c30bbc58-3e11-11eb-1e98-cda1cf1337e3
begin #one-variable gradient
	f(x)=x^2
	df(x)=gradient(f,x)
end

# ╔═╡ a055c570-3e12-11eb-07d2-955dda5ae765
f(3),df(3)

# ╔═╡ e3427b4e-3e12-11eb-3779-7fa565a2bc17
begin #multivariable
	g(x,y,z)=x.*42 .+ y^2 .+ sin(z)
	dg(x,y,z)=gradient(g,x,y,z)
end

# ╔═╡ f2632f96-3e13-11eb-373c-d5bc583b99e9
dg(4,0,0)

# ╔═╡ bcf3b9d8-3e14-11eb-0f98-7f9dde938646
begin
	x=rand(42,1)
	W=rand(1,42)
	U=rand(1)
	grads = gradient(()->(W*x*U)[1], params([W,U]))
end

# ╔═╡ 40d88bd4-3e15-11eb-2628-71a41b821c59
grads[W],grads[U]

# ╔═╡ 354afdd6-3e17-11eb-2b2a-8fc10a36d820
model=Chain(Dense(42,9),
	Dense(9,1))

# ╔═╡ 6932f126-3e17-11eb-3b87-e5bdbd2a3f54
model(x)

# ╔═╡ ae7c1ad2-3e17-11eb-074a-03df157b9566
params(model)

# ╔═╡ 33b26a94-3e18-11eb-1125-af41dc7537e7
begin
X,y=rand(42,1000),rand(42,1000)
loss(x, y) = mean((model(x).-y))
@epochs 5 train!(loss, params(model), [(X,y)], Descent(1e-3))
end

# ╔═╡ 47eeeb8a-3e19-11eb-2158-5bf7ceb2ed90
loss(X,y)

# ╔═╡ Cell order:
# ╠═864b72c0-3e11-11eb-2cd6-9951cdcb78c3
# ╠═c30bbc58-3e11-11eb-1e98-cda1cf1337e3
# ╠═a055c570-3e12-11eb-07d2-955dda5ae765
# ╠═e3427b4e-3e12-11eb-3779-7fa565a2bc17
# ╠═f2632f96-3e13-11eb-373c-d5bc583b99e9
# ╠═bcf3b9d8-3e14-11eb-0f98-7f9dde938646
# ╠═40d88bd4-3e15-11eb-2628-71a41b821c59
# ╠═588bf94e-3e17-11eb-216c-090eb107222e
# ╠═354afdd6-3e17-11eb-2b2a-8fc10a36d820
# ╠═6932f126-3e17-11eb-3b87-e5bdbd2a3f54
# ╠═ae7c1ad2-3e17-11eb-074a-03df157b9566
# ╠═33b26a94-3e18-11eb-1125-af41dc7537e7
# ╠═47eeeb8a-3e19-11eb-2158-5bf7ceb2ed90

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
"""
    Code for a Gaussian mask shape for use with CCF
Author: Eric Ford
Created: August 2020
"""

import SpecialFunctions: erf

"""   GaussianCCFMask
A truncated Gaussian mask with two parameters, it's stdandard deviation and where to truncate it, both as a velocity in m/s.
Mask weights are stored separately in a line list.
"""
struct GaussianCCFMask <: AbstractCCFMaskShape
    σ_sqrt2::Float64
    half_width_truncation::Float64
    normalization::Float64
    cdf_normalization::Float64

    """ GaussianCCFMask( σ ; half_truncation_width_in_σ=2 ) """
    function GaussianCCFMask(σ::Real, w::Real=2 )
        @assert 0 < σ <= 300000   # 300km/s is arbitrary choice for an upper limit
        @assert 0 < w <= 4
        norm = 1.0/(sqrt(2π)*σ*erf(w/(sqrt(2.0))))
        cdf_norm = 0.5/(erf(w/(sqrt(2.0))))
        new(σ*sqrt(2.0),σ*w,norm,cdf_norm)
    end

end

""" `GaussianCCFMask( inst ; scale_factor )` """
function GaussianCCFMask(inst::AbstractInstrument; σ_scale_factor::Real = 1, truncation_scale_factor::Real = default_gaussian_ccf_truncation_scale_factor )
    #=
    default_σ = 5000
    σ_sqrt2 = scale_factor * default_σ * sqrt(2.0)
    w = scale_factor * default_ccf_mask_v_width(inst)
    norm = 1.0/(σ_sqrt2*sqrt(π)*erf(w/σ_sqrt2))
    GaussianCCFMask(σ_sqrt2,w/2) # ,norm)
    =#
    σ = σ_scale_factor * default_ccf_mask_v_width(inst)
    w = truncation_scale_factor
    GaussianCCFMask(σ,w/2)
end

λ_min(m::GaussianCCFMask,λ::Real) = λ/calc_doppler_factor(m.half_width_truncation)
λ_max(m::GaussianCCFMask,λ::Real) = λ*calc_doppler_factor(m.half_width_truncation)

function integrate(m::GaussianCCFMask, v_lo::Real,v_hi::Real)
    #quadgk(m, v_lo, v_hi)[1]
    m.cdf_normalization*(erf(v_hi/m.σ_sqrt2)-erf(v_lo/m.σ_sqrt2))
end

""" Functor for returning PSF for Δv <= half_width.  """
function (m::GaussianCCFMask)(Δv::Real)
    if abs2(Δv) > abs2(m.half_width_truncation)
        return zero(Δv)
    else
        return m.normalization*exp(-(Δv/m.σ_sqrt2)^2)
    end
end

function mask_with_increased_fwhm(m::GaussianCCFMask, Δfwhm::Real )
    σ = m.σ_sqrt2/sqrt(2)
    w = m.half_width_truncation/σ
    #fwhm_orig = sqrt(8*log(2)) * σ
    fwhm_orig = 2*sqrt(log(2)) * m.σ_sqrt2
    fwhm_new = sqrt(fwhm_orig^2 + Δfwhm^2)
    σ_new = 0.5*fwhm_new/sqrt(2*log(2))
    return GaussianCCFMask(σ_new, w)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
"""
Collection of data models used for content parsing
"""
module Models

using Base: @kwdef
using Dates

export AbstractStonksRecord,
  AssetInfo,
  AssetPrice,
  ExchangeRate,
  IncomeStatement,
  BalanceSheet,
  CashflowStatement,
  Earnings

"""
Abstract type inteded to be subtyped by an data model.

All subtypes are constructed from keyword args and optional values default to `missing`.
"""
abstract type AbstractStonksRecord end

"""
Stores general information about a quoted symbol / ticker.

### Constructors
```julia
AssetInfo(;
  symbol::String,
  currency::String,
  name::Union{String,Missing}=missing,
  type::Union{String,Missing}=missing,
  exchange::Union{String,Missing}=missing,
  country::Union{String,Missing}=missing,
  industry::Union{String,Missing}=missing,
  sector::Union{String,Missing}=missing,
  timezone::Union{String,Missing}=missing,
  employees::Union{Int,Missing}=missing,
)
```
"""
@kwdef struct AssetInfo <: AbstractStonksRecord
  symbol::String
  currency::String
  name::Union{String,Missing} = missing
  type::Union{String,Missing} = missing
  exchange::Union{String,Missing} = missing
  country::Union{String,Missing} = missing
  industry::Union{String,Missing} = missing
  sector::Union{String,Missing} = missing
  timezone::Union{String,Missing} = missing
  employees::Union{Int,Missing} = missing
end
"""
Stores a time series datapoint with price information. Lowest frequency is daily.

### Constructors
```julia
AssetPrice(
  symbol::String,
  date::Date,
  close::Float64,
  open::Union{Float64,Missing} = missing,
  high::Union{Float64,Missing} = missing,
  low::Union{Float64,Missing} = missing,
  close_adjusted::Union{Float64,Missing} = missing,
  volume::Union{Integer,Missing} = missing,
)
end
```
"""
@kwdef struct AssetPrice <: AbstractStonksRecord
  symbol::String
  date::Date
  close::Float64
  open::Union{Float64,Missing} = missing
  high::Union{Float64,Missing} = missing
  low::Union{Float64,Missing} = missing
  close_adjusted::Union{Float64,Missing} = missing
  volume::Union{Integer,Missing} = missing
end

"""
Stores an exchange rate datapoint. Lowest frequency is daily.

### Constructors
```julia
ExchangeRate(;
  base::String,
  target::String,
  date::Date,
  rate::Float64,
)
```
"""
@kwdef struct ExchangeRate <: AbstractStonksRecord
  base::String
  target::String
  date::Date
  rate::Float64
end

"""
Stores a datapoint containing income statement information.
Follows normalized fields mapped to GAAP and IFRS taxonomies of the SEC.

### Constructors
```julia
IncomeStatement(;
  symbol::String,
  frequency::String,
  date::Date,
  currency::Union{String,Missing} = missing,
  total_revenue::Int64,
  cost_of_revenue::Int64,
  gross_profit::Int64,
  operating_income::Int64,
  selling_general_and_administrative::Union{Int64,Missing} = missing,
  research_and_development::Union{Int64,Missing} = missing,
  depreciation::Union{Int64,Missing} = missing, # Depreciation, Total
  depreciation_and_amortization::Union{Int64,Missing} = missing,
  income_before_tax::Union{Int64,Missing} = missing,
  income_tax_expense::Union{Int64,Missing} = missing,
  interest_expense::Union{Int64,Missing} = missing,
  interest_and_debt_expense::Union{Int64,Missing} = missing,
  ebit::Union{Int64,Missing} = missing,
  ebitda::Union{Int64,Missing} = missing,
  net_income::Int64,
  net_income_common_shares::Union{Int64,Missing} = missing,
)
```
"""
@kwdef struct IncomeStatement <: AbstractStonksRecord
  symbol::String
  frequency::String
  date::Date
  currency::Union{String,Missing} = missing
  total_revenue::Int64
  cost_of_revenue::Int64
  gross_profit::Int64
  operating_income::Int64
  selling_general_and_administrative::Union{Int64,Missing} = missing
  research_and_development::Union{Int64,Missing} = missing
  depreciation::Union{Int64,Missing} = missing # Depreciation, Total
  depreciation_and_amortization::Union{Int64,Missing} = missing
  income_before_tax::Union{Int64,Missing} = missing
  income_tax_expense::Union{Int64,Missing} = missing
  interest_expense::Union{Int64,Missing} = missing
  interest_and_debt_expense::Union{Int64,Missing} = missing
  #comprehensive_income_net_of_tax::Union{Int64,Missing} = missing
  ebit::Union{Int64,Missing} = missing
  ebitda::Union{Int64,Missing} = missing
  net_income::Int64
  net_income_common_shares::Union{Int64,Missing} = missing
end

"""
Stores a datapoint containing balance sheet information.
Follows normalized fields mapped to GAAP and IFRS taxonomies of the SEC.

### Constructors
```julia
BalanceSheet(;
  symbol::String,
  frequency::String,
  date::Date,
  currency::Union{String,Missing} = missing,
  total_assets::Int64,
  total_liabilities::Int64,
  total_shareholder_equity::Int64,
  cash_and_equivalents::Union{Int64,Missing} = missing,
  current_net_receivables::Union{Int64,Missing} = missing,
  inventory::Union{Int64,Missing} = missing,
  short_term_investments::Union{Int64,Missing} = missing,
  other_current_assets::Union{Int64,Missing} = missing,
  total_current_assets::Union{Int64,Missing} = missing,
  property_plant_equipment::Union{Int64,Missing} = missing,
  goodwill::Union{Int64,Missing} = missing,
  long_term_investments::Union{Int64,Missing} = missing,
  intangible_assets::Union{Int64,Missing} = missing,
  total_noncurrent_assets::Union{Int64,Missing} = missing,
  current_accounts_payable::Union{Int64,Missing} = missing,
  deferred_revenue::Union{Int64,Missing} = missing,
  short_term_debt::Union{Int64,Missing} = missing,
  other_current_liabilities::Union{Int64,Missing} = missing,
  total_current_liabilities::Union{Int64,Missing} = missing,
  current_debt::Union{Int64,Missing} = missing,
  current_long_term_debt::Union{Int64,Missing} = missing,
  long_term_debt::Union{Int64,Missing} = missing,
  long_term_debt_noncurrent::Union{Int64,Missing} = missing,
  capital_lease_obligations::Union{Int64,Missing} = missing,
  other_noncurrent_liabilities::Union{Int64,Missing} = missing,
  total_noncurrent_liabilities::Union{Int64,Missing} = missing,
  treasury_stock::Union{Int64,Missing} = missing,
  retained_earnings::Union{Int64,Missing} = missing,
  common_stock::Union{Int64,Missing} = missing,
  common_stock_shares_outstanding::Union{Int64,Missing} = missing,
)
```
"""
@kwdef struct BalanceSheet <: AbstractStonksRecord
  symbol::String
  frequency::String
  date::Date
  currency::Union{String,Missing} = missing
  total_assets::Int64
  total_liabilities::Int64
  total_shareholder_equity::Int64
  cash_and_equivalents::Union{Int64,Missing} = missing
  current_net_receivables::Union{Int64,Missing} = missing
  inventory::Union{Int64,Missing} = missing
  short_term_investments::Union{Int64,Missing} = missing
  other_current_assets::Union{Int64,Missing} = missing
  total_current_assets::Union{Int64,Missing} = missing
  property_plant_equipment::Union{Int64,Missing} = missing
  goodwill::Union{Int64,Missing} = missing
  long_term_investments::Union{Int64,Missing} = missing
  intangible_assets::Union{Int64,Missing} = missing
  total_noncurrent_assets::Union{Int64,Missing} = missing
  current_accounts_payable::Union{Int64,Missing} = missing
  deferred_revenue::Union{Int64,Missing} = missing
  short_term_debt::Union{Int64,Missing} = missing
  other_current_liabilities::Union{Int64,Missing} = missing
  total_current_liabilities::Union{Int64,Missing} = missing
  current_debt::Union{Int64,Missing} = missing
  current_long_term_debt::Union{Int64,Missing} = missing
  long_term_debt::Union{Int64,Missing} = missing
  long_term_debt_noncurrent::Union{Int64,Missing} = missing
  capital_lease_obligations::Union{Int64,Missing} = missing
  other_noncurrent_liabilities::Union{Int64,Missing} = missing
  total_noncurrent_liabilities::Union{Int64,Missing} = missing
  treasury_stock::Union{Int64,Missing} = missing
  retained_earnings::Union{Int64,Missing} = missing
  common_stock::Union{Int64,Missing} = missing
  common_stock_shares_outstanding::Union{Int64,Missing} = missing
end

"""
Stores a datapoint containing cashflow statement information.

### Constructors
```julia
CashflowStatement(;
  symbol::String,
  frequency::String,
  date::String,
  currency::String,
  operating_cashflow::Union{Int64,Missing} = missing,
  cashflow_investment::Union{Int64,Missing} = missing,
  cashflow_financing::Union{Int64,Missing} = missing,
  change_operating_liabilities::Union{Int64,Missing} = missing,
  change_receivables::Union{Int64,Missing} = missing,
  change_inventory::Union{Int64,Missing} = missing,
  change_cash_and_equivalents::Union{Int64,Missing} = missing,
  depreciation_and_amortization::Union{Int64,Missing} = missing,
  capital_expenditures::Union{Int64,Missing} = missing,
  dividend_payout::Union{Int64,Missing} = missing,
  stock_repurchase::Union{Int64,Missing} = missing,
  net_income::Union{Int64,Missing} = missing,
)
```
"""
@kwdef struct CashflowStatement <: AbstractStonksRecord
  symbol::String
  frequency::String
  date::Date
  currency::String
  operating_cashflow::Union{Int64,Missing} = missing
  cashflow_investment::Union{Int64,Missing} = missing
  cashflow_financing::Union{Int64,Missing} = missing
  change_operating_liabilities::Union{Int64,Missing} = missing
  change_receivables::Union{Int64,Missing} = missing
  change_inventory::Union{Int64,Missing} = missing
  change_cash_and_equivalents::Union{Int64,Missing} = missing
  depreciation_and_amortization::Union{Int64,Missing} = missing
  capital_expenditures::Union{Int64,Missing} = missing
  dividend_payout::Union{Int64,Missing} = missing
  stock_repurchase::Union{Int64,Missing} = missing
  net_income::Union{Int64,Missing} = missing
end

"""
Stores a datapoint containing earnings (per share) information.
### Constructors
```julia
Earnings(;
  symbol::String,
  frequency::String,
  date::Date,
  currency::Union{String,Missing} = missing,
  actual::Float16,
  estimate::Union{Float16,Missing} = missing,
)
```
"""
@kwdef struct Earnings <: AbstractStonksRecord
  symbol::String
  frequency::String
  date::Date
  currency::Union{String,Missing} = missing
  actual::Float16
  estimate::Union{Float16,Missing} = missing
end

# @kwdef struct EconomicIndicator <: AbstractStonksRecord
#   symbol::String
#   name::String
#   date::Date
#   value::Number
# end

end


text:
module BarotropicQG

export
  Problem,
  set_zeta!,
  updatevars!,

  energy,
  enstrophy,
  meanenergy,
  meanenstrophy,
  dissipation,
  work,
  drag

using Reexport

@reexport using FourierFlows

using FFTW: rfft
using LinearAlgebra: mul!, ldiv!
using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2

abstract type BarotropicQGVars <: AbstractVars end
abstract type BarotropicQGForcedVars <: BarotropicQGVars end

const realscalarvars = [:U]
const physicalvars = [:q, :zeta, :psi, :u, :v]
const transformvars = [ Symbol(var, :h) for var in physicalvars ]
const forcedvars = [:Fqh]
const stochforcedvars = [:prevsol]

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a BarotropicQG turbulence problem.
"""

function Problem(;
    # Numerical parameters
            nx = 256,
            Lx = 2π,
            ny = nx,
            Ly = Lx,
            dt = 0.01,
    # Physical parameters
            f0 = 1.0,
          beta = 0.0,
           eta = nothing,
    # Drag and/or hyper-/hypo-viscosity
            nu = 0.0,
           nnu = 1,
            mu = 0.0,
   # Timestepper and eqn options
       stepper = "RK4",
        calcFU = nothingfunction,
        calcFq = nothingfunction,
    stochastic = false,
             T = Float64)

  # the grid
  gr  = TwoDGrid(nx, Lx, ny, Ly; T=T)
  x, y = gridpoints(gr)

  # topographic PV
  if eta==nothing
    eta = 0*x
    etah = rfft(eta)
  end

  if typeof(eta)!=Array{T, 2} #this is true if eta was passes in Problem as a function
    pr = Params(gr, f0, beta, eta, mu, nu, nnu, calcFU, calcFq)
  else
    pr = Params{T}(f0, beta, eta, rfft(eta), mu, nu, nnu, calcFU, calcFq)
  end

  if calcFq == nothingfunction && calcFU == nothingfunction
      vs = Vars(gr)
  else
      if stochastic ==  false
          vs = ForcedVars(gr)
      elseif stochastic == true
          vs = StochasticForcedVars(gr)
      end
  end
  eq = Equation(pr, gr)
  FourierFlows.Problem(eq, stepper, dt, gr, vs, pr)
end

InitialValueProblem(; kwargs...) = Problem(; kwargs...)
ForcedProblem(; kwargs...) = Problem(; kwargs...)

# ----------
# Parameters
# ----------

"""
    Params(g::TwoDGrid, f0, beta, FU, eta, mu, nu, nnu, calcFU, calcFq)

Returns the params for an unforced two-dimensional barotropic QG problem.
"""
struct Params{T} <: AbstractParams
  f0::T                      # Constant planetary vorticity
  beta::T                    # Planetary vorticity y-gradient
  eta::Array{T,2}            # Topographic PV
  etah::Array{Complex{T},2}  # FFT of Topographic PV
  mu::T                      # Linear drag
  nu::T                      # Viscosity coefficient
  nnu::Int                   # Hyperviscous order (nnu=1 is plain old viscosity)
  calcFU::Function    # Function that calculates the forcing F(t) on
                      # domain-averaged zonal flow U(t)
  calcFq!::Function   # Function that calculates the forcing on QGPV q
end

"""
    Params(g::TwoDGrid, f0, beta, eta::Function, mu, nu, nnu, calcFU, calcFq)

Constructor for Params that accepts a generating function for the topographic PV.
"""
function Params(g::AbstractGrid{T}, f0, beta, eta::Function, mu, nu, nnu, calcFU, calcFq)  where T
  x, y = gridpoints(g)
  etagrid = eta(x, y)
  etah = rfft(etagrid)
  Params{T}(f0, beta, etagrid, etah, mu, nu, nnu, calcFU, calcFq)
end

# ---------
# Equations
# ---------

"""
    Equation(p, g)

Returns the equation for two-dimensional barotropic QG problem with params p and grid g.
"""
function Equation(p::Params, g::AbstractGrid{T}) where T
  L = @. -p.mu - p.nu*g.Krsq^p.nnu + im*p.beta*g.kr*g.invKrsq
  L[1, 1] = 0
  FourierFlows.Equation(L, calcN!, g)
end

# ----
# Vars
# ----

varspecs = cat(
  getfieldspecs(realscalarvars, :(Array{T,0})),
  getfieldspecs(physicalvars, :(Array{T,2})),
  getfieldspecs(transformvars, :(Array{Complex{T},2})),
  dims=1)

forcedvarspecs = cat(varspecs, getfieldspecs(forcedvars, :(Array{Complex{T},2})), dims=1)
stochforcedvarspecs = cat(forcedvarspecs, getfieldspecs(stochforcedvars, :(Array{Complex{T},2})), dims=1)

# Construct Vars types
eval(structvarsexpr(:Vars, varspecs; parent=:BarotropicQGVars))
eval(structvarsexpr(:ForcedVars, forcedvarspecs; parent=:BarotropicQGForcedVars))
eval(structvarsexpr(:StochasticForcedVars, stochforcedvarspecs; parent=:BarotropicQGForcedVars))

"""
    Vars(g)

Returns the vars for unforced two-dimensional barotropic QG problem with grid g.
"""
function Vars(g::AbstractGrid{T}) where T
  U = Array{T,0}(undef, ); U[] = 0
  @createarrays T (g.nx, g.ny) q u v psi zeta
  @createarrays Complex{T} (g.nkr, g.nl) qh uh vh psih zetah
  Vars(U, q, zeta, psi, u, v, qh, zetah, psih, uh, vh)
end

"""
    ForcedVars(g)

Returns the vars for forced two-dimensional barotropic QG problem with grid g.
"""
function ForcedVars(g::AbstractGrid{T}) where T
  v = Vars(g)
  Fqh = zeros(Complex{T}, (g.nkr, g.nl))
  ForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., Fqh)
end

"""
    StochasticForcedVars(g)

Returns the vars for stochastically forced two-dimensional barotropic QG problem with grid g.
"""
function StochasticForcedVars(g::AbstractGrid{T}) where T
  v = ForcedVars(g)
  prevsol = zeros(Complex{T}, (g.nkr, g.nl))
  StochasticForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., prevsol)
end

# -------
# Solvers
# -------

function calcN_advection!(N, sol, t, cl, v, p, g)
  # Note that U = sol[1, 1]. For all other elements ζ = sol
  v.U[] = sol[1, 1].re
  @. v.zetah = sol
  v.zetah[1, 1] = 0

  @. v.uh =  im * g.l  * g.invKrsq * v.zetah
  @. v.vh = -im * g.kr * g.invKrsq * v.zetah

  ldiv!(v.zeta, g.rfftplan, v.zetah)
  ldiv!(v.u, g.rfftplan, v.uh)
  v.psih .= v.vh # FFTW's irfft destroys its input; v.vh is needed for N[1, 1]
  ldiv!(v.v, g.rfftplan, v.psih)

  @. v.q = v.zeta + p.eta
  @. v.u = (v.U[] + v.u)*v.q # (U+u)*q
  @. v.v = v.v*v.q # v*q

  mul!(v.uh, g.rfftplan, v.u) # \hat{(u+U)*q}
  # Nonlinear advection term for q (part 1)
  @. N = -im*g.kr*v.uh # -∂[(U+u)q]/∂x
  mul!(v.uh, g.rfftplan, v.v) # \hat{v*q}
  @. N += - im*g.l*v.uh # -∂[vq]/∂y
end

function calcN!(N, sol, t, cl, v, p, g)
  calcN_advection!(N, sol, t, cl, v, p, g)
  addforcing!(N, sol, t, cl, v, p, g)
  if p.calcFU != nothingfunction
    # 'Nonlinear' term for U with topographic correlation.
    # Note: < v*eta > = sum( conj(vh)*eta ) / (nx^2*ny^2) if fft is used
    # while < v*eta > = 2*sum( conj(vh)*eta ) / (nx^2*ny^2) if rfft is used
    N[1, 1] = p.calcFU(t) + 2*sum(conj(v.vh).*p.etah).re / (g.nx^2.0*g.ny^2.0)
  end
  nothing
end

addforcing!(N, sol, t, cl, v::Vars, p, g) = nothing

function addforcing!(N, sol, t, cl, v::ForcedVars, p, g)
  p.calcFq!(v.Fqh, sol, t, cl, v, p, g)
  @. N += v.Fqh
  nothing
end

function addforcing!(N, sol, t, cl, v::StochasticForcedVars, p, g)
  if t == cl.t # not a substep
    @. v.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    p.calcFq!(v.Fqh, sol, t, cl, v, p, g)
  end
  @. N += v.Fqh
  nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(v, s, g)

Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(sol, v, p, g)
  v.U[] = sol[1, 1].re
  @. v.zetah = sol
  v.zetah[1, 1] = 0.0

  @. v.psih = -v.zetah * g.invKrsq
  @. v.uh = -im * g.l  * v.psih
  @. v.vh =  im * g.kr * v.psih

  ldiv!(v.zeta, g.rfftplan, deepcopy(v.zetah))
  ldiv!(v.psi, g.rfftplan, deepcopy(v.psih))
  ldiv!(v.u, g.rfftplan, deepcopy(v.uh))
  ldiv!(v.v, g.rfftplan, deepcopy(v.vh))

  @. v.q = v.zeta + p.eta
  nothing
end

updatevars!(prob) = updatevars!(prob.sol, prob.vars, prob.params, prob.grid)

"""
    set_zeta!(prob, zeta)
    set_zeta!(s, v, g, zeta)

Set the solution sol as the transform of zeta and update variables v
on the grid g.
"""
function set_zeta!(sol, v::Vars, p, g, zeta)
  mul!(v.zetah, g.rfftplan, zeta)
  v.zetah[1, 1] = 0.0
  @. sol = v.zetah

  updatevars!(sol, v, p, g)
  nothing
end

function set_zeta!(sol, v::BarotropicQGForcedVars, p, g, zeta)
  v.U[] = deepcopy(sol[1, 1])
  mul!(v.zetah, g.rfftplan, zeta)
  v.zetah[1, 1] = 0.0
  @. sol = v.zetah
  sol[1, 1] = v.U[]

  updatevars!(sol, v, p, g)
  nothing
end

set_zeta!(prob, zeta) = set_zeta!(prob.sol, prob.vars, prob.params, prob.grid, zeta)

"""
    set_U!(prob, U)
    set_U!(sol, v, g, U)

Set the (kx, ky)=(0, 0) part of solution sol as the domain-average zonal flow U.
"""
function set_U!(sol, v, p, g, U::Float64)
  sol[1, 1] = U
  updatevars!(sol, v, p, g)
  nothing
end

set_U!(prob, U::Float64) = set_U!(prob.sol, prob.vars, prob.params, prob.grid, U)

"""
Calculate the domain-averaged kinetic energy.
"""
function energy(prob)
  sol, g = prob.sol, prob.grid
  0.5*(parsevalsum2(g.kr.*g.invKrsq.*sol, g)
        + parsevalsum2(g.l.*g.invKrsq.*sol, g))/(g.Lx*g.Ly)
end

"""
Returns the domain-averaged enstrophy.
"""
function enstrophy(prob)
  sol, v, g = prob.sol, prob.vars, prob.grid
  @. v.uh = sol
  v.uh[1, 1] = 0
  0.5*parsevalsum2(v.uh, g)/(g.Lx*g.Ly)
end

"""
Returns the energy of the domain-averaged U.
"""
meanenergy(prob) = real(0.5*prob.sol[1, 1].^2)

"""
Returns the enstrophy of the domain-averaged U.
"""
meanenstrophy(prob) = real(prob.params.beta*prob.sol[1, 1])

"""
    dissipation(prob)
    dissipation(s, v, p, g)

Returns the domain-averaged dissipation rate. nnu must be >= 1.
"""
@inline function dissipation(sol, v, p, g)
  @. v.uh = g.Krsq^(p.nnu-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.nu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline dissipation(prob) = dissipation(prob.sol, prob.vars, prob.params, prob.grid)

"""
    work(prob)
    work(s, v, p, g)

Returns the domain-averaged rate of work of energy by the forcing Fqh.
"""
@inline function work(sol, v::ForcedVars, g)
  @. v.uh = g.invKrsq * sol * conj(v.Fqh)
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline function work(sol, v::StochasticForcedVars, g)
  @. v.uh = g.invKrsq * (v.prevsol + sol)/2.0 * conj(v.Fqh) # Stratonovich
  # @. v.uh = g.invKrsq * v.prevsol * conj(v.Fqh)             # Ito
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline work(prob) = work(prob.sol, prob.vars, prob.grid)

"""
    drag(prob)
    drag(s, v, p, g)

Returns the extraction of domain-averaged energy by drag mu.
"""
@inline function drag(sol, v, p, g)
  @. v.uh = g.Krsq^(-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.mu/(g.Lx*g.Ly)*FourierFlows.parsevalsum(v.uh, g)
end

@inline drag(prob) = drag(prob.sol, prob.vars, prob.params, prob.grid)

end # module


text:
#
# spg2.jl --
#
# Implements Spectral Projected Gradient Method (Version 2: "continuous
# projected gradient direction") to find the local minimizers of a given
# function with convex constraints, described in:
#
# [1] E. G. Birgin, J. M. Martinez, and M. Raydan, "Nonmonotone spectral
#     projected gradient methods on convex sets", SIAM Journal on Optimization
#     10, pp. 1196-1211 (2000).
#
# [2] E. G. Birgin, J. M. Martinez, and M. Raydan, "SPG: software for
#     convex-constrained optimization", ACM Transactions on Mathematical
#     Software (TOMS) 27, pp. 340-349 (2001).
#
# ----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2014-2019: Éric Thiébaut.
#
# ----------------------------------------------------------------------------

module SPG

export
    spg2

using Printf
using ..OptimPack:
    DenseVariable,
    DenseVariableSpace,
    combine!,
    dot,
    norm2,
    norminf,
    wrap

mutable struct Info{T,N}
    x::Array{T,N}      # current iterate
    xbest::Array{T,N}  # best solution found so far
    f::Float64         # function value at `x`
    fbest::Float64     # function value at `xbest`
    pgtwon::Float64    # Euclidean norm of the projected gradient at `x`
    pginfn::Float64    # infinite norm of the projected gradient at `x`
    fcnt::Int          # number of function and gradient evaluations
    pcnt::Int          # number of projections
    iter::Int          # number of iterations
    status::Symbol     # status of the algorithm
end

Info(x::DenseArray{T,N}, xbest::DenseArray{T,N}) where {T,N} =
    Info{T,N}(x, xbest, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, :SEARCHING)

"""
# Spectral Projected Gradient Method

The `spg2` method implements the Spectral Projected Gradient Method (Version 2:
"continuous projected gradient direction") to find the local minimizers of a
given function with convex constraints, described in the references below.  A
typical use is:

```julia
spg2(fg!, prj!, x0, m=10) -> res
```

The user must supply the functions `fg!` and `prj!` to evaluate the objective
function and its gradient and to project an arbitrary point onto the feasible
region.  These functions must be defined as in the following pseudo-code:

```julia
function fg!(x::T, g::T) where {T}
   g[:] = gradient_at(x)
   return function_value_at(x)
end

function prj!(dst::T, src::T) where {T}
    dst[:] = projection_of(src)
    return dst
end
```

For instance, to constrain the variables to be all nonnegative, the following
projector should be specified:

```julia
function prj!(dst::AbstractArray{T,N},
              src::AbstractArray{T,N}) where {T,N}
    @assert axes(dst) == axes(src)
    @inbounds @simd for i in eachindex(dst, src)
        dst[i] = max(src[i], zero(T))
    end
    return dst
end
```

Argument `x0` is the initial solution and argument `m` is the number of
previous function values to be considered in the nonmonotone line search.  If
`m ≤ 1`, then a monotone line search with Armijo-like stopping criterion will
be used.

The returned value `res` is an instance of `SPG.Info` storing information about
the final iterate.

The following keywords are available:

* `eps1` specifies the stopping criterion `‖pg‖_∞ ≤ eps1` with `pg` the
  projected gradient.  By default, `eps1 = 1e-6`.

* `eps2` specifies the stopping criterion `‖pg‖_2 ≤ eps2` with `pg` the
  projected gradient.  By default, `eps2 = 1e-6`.

* `eta` specifies a scaling parameter for the gradient.  The projected gradient
  is computed as:

  ```
  (x - prj(x - eta*g))/eta
  ```

  (with `g` the gradient at `x`) instead of `x - prj(x - g)` which corresponds
  to the default behavior (same as if `eta=1`) and is usually used in
  methodological publications although it does not scale correctly (for
  instance, if you make a change of variables or simply multiply the function
  by some factor).

* `ftol` specify the parameter of the nonmonotone Armijo-like stopping
  criterion.  By default: `ftol = 1e-4`.

* `lmin` and `lmax` specify safeguard bounds for the steplength.  By default,
  `lmin = 1e-30` and `lmax = 1e+30`.

* `amin` and `amax` specify safeguardbounds for the qiadratic interpolation.
   By default, `amin = 0.1` and `amax = 0.9`.

* `maxit` specifies the maximum number of iterations.

* `maxfc` specifies the maximum number of function evaluations.

* `verb` indicates whether to print some information at each iteration.

* `printer` specifies a subroutine to print some information at each iteration.
  This subroutine will be called as `printer(io, ws)` with `io` the output
  stream and `ws` an instance of `SPG.Info` with information about the current
  iterate.

* `io` specifes the output stream for iteration information.  It is `stdout` by
  default.

The `SPG.Info` type has the following members:

* `x` is the current iterate.
* `f` is the function value.
* `xbest` is the best solution so far.
* `fbest` is the function value at `xbest`.
* `pginfn` is the infinite norm of the projected gradient.
* `pgtwon` is the Eucliddean norm of the projected gradient.
* `iter` is the number of iterations.
* `fcnt` is the number of function (and gradient) evaluations.
* `pcnt` is the number of projections.
* `status` indicates the type of termination:

  | Status                  | Reason                                          |
  |:------------------------|:------------------------------------------------|
  | `:SEARCHING`            | Work in progress                                |
  | `:INFNORM_CONVERGENCE`  | Convergence of projected gradient infinite-norm |
  | `:TWONORM_CONVERGENCE`  | Convergence of projected gradient 2-norm        |
  | `:TOO_MANY_ITERATIONS`  | Too many iterations                             |
  | `:TOO_MANY_EVALUATIONS` | Too many function evaluations                   |

## References

* E. G. Birgin, J. M. Martinez, and M. Raydan, "Nonmonotone spectral projected
  gradient methods on convex sets", SIAM Journal on Optimization 10,
  pp. 1196-1211 (2000).

* E. G. Birgin, J. M. Martinez, and M. Raydan, "SPG: software for
  convex-constrained optimization", ACM Transactions on Mathematical Software
  (TOMS) 27, pp. 340-349 (2001).
"""
function spg2(fg!::Function,
              prj!::Function,
              x0::AbstractArray{T,N},
              m::Integer = 10;
              maxit::Integer = typemax(Int),
              maxfc::Integer = typemax(Int),
              eps1::Real = 1e-6,
              eps2::Real = 1e-6,
              eta::Real = 1.0,
              ftol::Real = 1e-4,
              lmin::Real = 1e-30,
              lmax::Real = 1e+30,
              amin::Real = 0.1,
              amax::Real = 0.9,
              printer = nothing,
              io::IO = stdout,
              verb::Bool = false) where {T<:AbstractFloat,N}
    x = copyto!(Array{T}(undef, size(x0)), x0)
    _spg!(fg!, prj!, x, Int(m), Int(maxit), Int(maxfc),
          Float64(eps1), Float64(eps2), Float64(eta), Float64(ftol),
          Float64(lmin), Float64(lmax), Float64(amin), Float64(amax),
          printer, io, verb)
end

function _spg!(fg!::Function,
               prj!::Function,
               x::Array{T,N},
               m::Int,
               maxit::Int,
               maxfc::Int,
               eps1::Float64,
               eps2::Float64,
               eta::Float64,
               ftol::Float64,
               lmin::Float64,
               lmax::Float64,
               amin::Float64,
               amax::Float64,
               printer,
               io::IO,
               verb::Bool) where {T<:AbstractFloat,N}

    @assert m ≥ 1
    @assert maxit ≥ 0
    @assert maxfc ≥ 0
    @assert eps1 ≥ 0
    @assert eps2 ≥ 0
    @assert eta > 0
    @assert ftol > 0
    @assert 0 < lmin < lmax
    @assert 0 < amin < amax

    # Allocate workspace arrays.
    dims = size(x)
    g = Array{T}(undef, dims)
    d = Array{T}(undef, dims)
    x0 = Array{T}(undef, dims)
    g0 = Array{T}(undef, dims)
    xbest = Array{T}(undef, dims)
    lastfv = Array{Float64}(undef, m)

    # Wrap OptimPack variables over Julia workspace arrays (making a few
    # aliases to save memory).
    space = DenseVariableSpace(T, dims)
    Vx = wrap(space, x)
    Vg = wrap(space, g)
    Vd = wrap(space, d)
    Vs = Vx0 = wrap(space, x0)
    Vy = Vg0 = wrap(space, g0)

    # Initialization.
    ws = Info(x, xbest)

    # Project initial guess.
    prj!(x, x)
    ws.pcnt += 1

    # Evaluate function and gradient.
    f = Float64(fg!(x, g))
    ws.fcnt += 1

    # Initialize best solution and best function value.
    ws.fbest = f
    copyto!(xbest, x)

    # List of `m` previous function values.
    fill!(lastfv, Inf)
    # FIXME: should be: fill!(lastfv, f) # or -Inf
    lastfv[1] = f

    # Main loop.
    while true

        # Compute continuous projected gradient (and its norms) using the
        # storage of `d` as: `pg = (x - prj(x - eta*g))/eta`.
        combine!(Vd, 1, Vx, -eta, Vg)
        prj!(d, d)
        combine!(Vd, 1/eta, Vx, -1/eta, Vd)
        ws.pcnt += 1
        ws.pgtwon = norm2(Vd)
        ws.pginfn = norminf(Vd)

        # Print iteration information.
        if printer !== nothing
            printer(io, ws)
        end
        if verb
            @printf("ITER = %-5d  EVAL = %-5d  PROJ = %-5d  F(%s) =%24.17e  ||PG||oo = %17.10e\n",
                    ws.iter, ws.fcnt, ws.pcnt, (f ≤ ws.fbest ? "+" : "-"), f, ws.pginfn)
        end

        # Test stopping criteria.
        if ws.pginfn ≤ eps1
            # Gradient infinite-norm stopping criterion satisfied, stop.
            ws.status = :INFNORM_CONVERGENCE
            return ws
        end
        if ws.pgtwon ≤ eps2
            # Gradient 2-norm stopping criterion satisfied, stop.
            ws.status = :TWONORM_CONVERGENCE
            return ws
        end
        if ws.iter ≥ maxit
            # Maximum number of iterations exceeded, stop.
            ws.status = :TOO_MANY_ITERATIONS
            return ws
        end
        if ws.fcnt ≥ maxfc
            # Maximum number of function evaluations exceeded, stop.
            ws.status = :TOO_MANY_EVALUATIONS
            return ws
        end

        # Store function value for the nonmonotone line search and
        # find maximum function value since m last calls.
        if m > 1
            lastfv[(ws.iter%m) + 1] = f
            fmax = maximum(lastfv)
        else
            fmax = f
        end

        # Compute spectral steplength.
        if ws.iter == 0
            # Initial steplength. (FIXME: check type stability)
            lambda = clamp(1/ws.pginfn, lmin, lmax)
        else
            combine!(Vs, 1, Vx, -1, Vx0)
            combine!(Vy, 1, Vg, -1, Vg0)
            sty = dot(Vs, Vy)
            if sty > 0
                # Safeguarded Barzilai & Borwein spectral steplength.
                sts = dot(Vs, Vs)
                lambda = min(lmax, max(lmin, sts/sty))
            else
                lambda = lmax
            end
        end

        # Save current point.
        copyto!(x0, x)
        copyto!(g0, g)
        f0 = f

        # Compute the spectral projected gradient direction and <G,D>
        combine!(Vx, 1, Vx, -lambda, Vg)
        prj!(x, x)
        ws.pcnt += 1
        combine!(Vd, 1, Vx, -1, Vx0) # d = x - x0
        delta = dot(Vg0, Vd)

        # Nonmonotone line search.
        stp = 1.0 # Step length for first trial.
        while true
            # Evaluate function and gradient at trial point.
            f = Float64(fg!(x, g))
            ws.fcnt += 1

            # Compare the new function value against the best function
            # value and, if smaller, update the best function value and the
            # corresponding best point.
            if f < ws.fbest
                ws.fbest = f
                copyto!(xbest, x)
            end

            # Test stopping criteria.
            if f ≤ fmax + stp*ftol*delta
                # Nonmonotone Armijo-like stopping criterion satisfied, stop.
                break
            end
            if ws.fcnt ≥ maxfc
                # Maximum number of function evaluations exceeded, stop.
                ws.status = :TOO_MANY_EVALUATIONS
                return ws
            end

            # Safeguarded quadratic interpolation.
            q = -delta*(stp*stp)
            r = (f - f0 - stp*delta)*2
            if r > 0 && amin*r ≤ q ≤ amax*stp*r
                stp = q/r
            else
                stp /= 2
            end

            # Compute trial point.
            combine!(Vx, 1, Vx0, stp, Vd) # x = x0 + stp*d
        end

        if ws.status != :SEARCHING
            # The number of function evaluations was exceeded inside the line
            # search.
            return ws
        end

        # Proceed with next iteration.
        ws.iter += 1
    end
end

end # module


text:
mutable struct Edge
    rev::Int
    from::Int
    to::Int
    cap::Int
end

function push!(G::Vector{Vector{Edge}}, from, to, cap)
    fromrev = length(G[from]) + 1
    torev = length(G[to]) + 1
    Base.push!(G[from], Edge(torev, from, to, cap))
    Base.push!(G[to], Edge(fromrev, to, from, 0))
end

function run_flow!(G::Vector{Vector{Edge}}, e::Edge, f)
    e.cap -= f
    G[e.to][e.rev].cap += f
    return G
end

function fodfs(G, v, t, f, seens)
    v == t && return f

    seens[v] = true
    for e ∈ G[v]
        if seens[e.to] || e.cap == 0
            # Skip
        else
            flow = fodfs(G, e.to, t, min(f, e.cap), seens)
            if flow == 0
                # Skip
            else
                run_flow!(G, e, flow)
                return flow
            end
        end
    end

    return 0
end

function FordFulkerson(G, s, t)
    INF = 10^9
    res = 0

    while true
        flow = fodfs(G, s, t, INF, falses(length(G)))
        if flow == 0
            return res
        else
            res += flow
        end
    end

    return 0
end

function main()
    r, c = parse.(Int, split(readline()))
    ss = [readline() for _ = 1:r]
    N = r * c
    V = N
    G = [Vector{Edge}() for _ = 1:N+2] # 2 super nodes
    for i = 1:r, j = 1:c
        left = (i - 1) * c + j
        s = ss[i][j]
        if s == '*'
            V -= 1
        else
            if (i + j) % 2 == 0
                for (dx, dy) ∈ zip([0, 1, 0, -1], [1, 0, -1, 0])
                    i′ = i + dy
                    j′ = j + dx

                    if 1 ≤ i′ ≤ r && 1 ≤ j′ ≤ c && ss[i′][j′] == '.'
                        right = (i′ - 1) * c + j′
                        push!(G, left, right, 1)
                    end
                end
                push!(G, N+1, left, 1)
            else
                push!(G, left, N+2, 1)
            end
        end
    end

    println(V - FordFulkerson(G, N+1, N+2))
end

main()


text:
import FireflyModule.Firefly

function attractiveness(source_node::Firefly, destination_node::Firefly, λ::Float)::Float
    """Firefly’s attractiveness is proportional to the light
    intensity seen by adjacent fireflies, we can now define the
    attractiveness β of a firefly depending on distance between
    two fireflies.
    """
    β = source_node.cost
    r, _ = hamming_distance(source_node, destination_node)
    return β * exp(-λ ^ (r^2))
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
using CollisionDetection

"""
Returns a mesh on the same vertexbuffer as the input mesh. The submesh
will be a mesh of dimension k containing all the k-cells that are in mesh
and that fulfill the predicate pred.

`pred` is a function with signature `pred(cell) -> Bool` returning true if
the simplex is to be added to the submesh under construction.
"""
function submesh(pred, mesh::AbstractMesh)

    # kcells = similar(mesh.faces)
    kcells = cells(mesh)

    sub2sup = zeros(Int, 0)
    # sup2sub = zeros(Int, length(mesh))

    j = 0
    # for i in 1:length(kcells)
    for (i,kcell) = enumerate(cells(mesh))

        # kcell = mesh.faces[i]
        # kcell = kcells[i]
        if pred(kcell)
            # kcells[j] = kcell
            push!(sub2sup, i)
            # sup2sub[i] = j
            j += 1
        end

    end

    kcells = kcells[1:j]

    # sub = Mesh(mesh.vertices, kcells)
    return SubMesh(mesh, sub2sup)
end

"""
    submesh(selection, mesh)

Create a submesh from `mesh` comprising those elements that overlap with elements from `selection`. It is assumed that `selection` and `mesh` have the same dimension.
"""
function submesh(sm::Mesh, bm::Mesh)
    overlaps = overlap_gpredicate(sm)
    in_smallmesh = c -> overlaps(simplex(vertices(bm,c)))
    submesh(in_smallmesh, bm)
end

mutable struct SubMesh{U,D1,T} <: AbstractMesh{U,D1,T}
    supermesh::AbstractMesh{U,D1,T}
    sub2sup::Vector{Int}
    sup2sub::Vector{Int}
    cells::Vector{SVector{D1,Int}}
end

function SubMesh(supermesh, sub2sup)

    sup2sub = zeros(Int, length(supermesh))
    @inbounds for (i,j) in enumerate(sub2sup)
        sup2sub[j] = i
    end

    Cells = cells(supermesh)[sub2sup]
    SubMesh(supermesh, sub2sup, sup2sub, Cells)
end

vertextype(m::SubMesh) = vertextype(m.supermesh)
celltype(m::SubMesh) = celltype(m.supermesh)
coordtype(m::SubMesh) = coordtype(m.supermesh)
dimension(m::SubMesh) = dimension(m.supermesh)
universedimension(m::SubMesh) = universedimension(m.supermesh)

vertexarray(m::SubMesh) = vertexarray(m.supermesh)
function cellarray(m::SubMesh)
    [cells(m.supermesh)[i][j] for i in m.sub2sup, j in 1:dimension(m)+1]
end

vertices(m::SubMesh) = vertices(m.supermesh)
vertices(m::SubMesh, cell) = vertices(m.supermesh, cell)
numvertices(m::SubMesh) = numvertices(m.supermesh)

cells(m::SubMesh) = m.cells #m.supermesh.faces[m.sub2sup]
numcells(m::SubMesh) = length(m.sub2sup)

issubmesh(sub, sup) = (sub == sup)
issubmesh(sub::SubMesh, sup) = (sub.supermesh == sup)
issubmesh(sub::SubMesh, sup::SubMesh) = (sub == sup)

extend(sm, cell) = sm.sub2sup[cell]
restrict(sm, cell) = sm.sup2sub[cell]

chart(mesh::SubMesh, cell) = simplex(vertices(mesh,cell))


text:
function Base.showerror(error::HODException)

end


text:
abstract type PullType end

struct PullAction{pt <: PullType, ARGS <: Tuple} <: Function
    f::Function
    args::ARGS
end

PullAction(f, args, pt = StandardPull) = PullAction{pt, typeof(args)}(f, args)

pull_args(pa::PullAction) = pull_args(pa.args)
pull_args(args) = map(args) do arg
    pull!(arg)
end

valid_args(args) = all(args) do arg
    !(typeof(arg) <: Signal) ? true : valid(arg)
end
valid(pa::PullAction) = valid_args(pa.args)

(pa::PullAction)() = pa.f(pull_args(pa)...)


text:
"""
Type for storing NamedD
"""
mutable struct NamedD{T<:Union{Srcs,Recs}}
    grid::StepRangeLen{Float64}
    d::NamedArrays.NamedArray{
        Array{Float64,2},
        1,
        Array{Array{Float64,2},1},
        Tuple{OrderedCollections.OrderedDict{Symbol,Int64}},
    }
    sr::T
    #	NamedD{T}(grid,d,sr)=any([(size(dd)≠(length(grid),sr.n)) for dd in d]) ? error("NamedD construction") : new(grid,d,sr)
end

"""
Simplify getting properties
"""
function Base.getindex(obj::NamedD, sym::Symbol)
    if (sym ∈ [:nr, :ns, :n])
        return obj.sr.n
    else
        return Base.getindex(obj.d, sym)
    end
end

"""
Initialize with zeros, given grid and names.
"""
function Base.zero(
    ::Type{NamedD},
    grid::StepRangeLen,
    s::Union{Srcs,Recs},
    fields::Vector{Symbol},
)
    for f in fields
        @assert eval(f) <: Fields
    end
    return NamedD(
        grid,
        NamedArray([zeros(length(grid), s.n) for i in fields], (fields,)),
        s,
    )
end
NamedD(grid, s::Union{Srcs,Recs}, fields::Vector{Symbol}) = zero(NamedD, grid, s, fields)

VNamedD = Union{Array{NamedD{Recs},1},Array{NamedD{Srcs},1}}

function Array{NamedD{Recs},1}(
    grid::StepRangeLen,
    ss::SSrcs,
    sr::Vector{Recs},
    fields::Vector{Symbol},
)
    return [NamedD(grid, sr[i], fields) for i = 1:ss.n]
end
function Array{NamedD{Srcs},1}(
    grid::StepRangeLen,
    ss::SSrcs,
    sr::Vector{Srcs},
    fields::Vector{Symbol},
)
    return [NamedD(grid, sr[i], fields) for i = 1:ss.n]
end

function Array{NamedD{Srcs},1}(
    grid::StepRangeLen,
    ss::SSrcs,
    sr::Srcs,
    fields::Vector{Symbol},
)
    return [NamedD(grid, sr, fields) for i = 1:ss.n]
end
function Array{NamedD{Recs},1}(
    grid::StepRangeLen,
    ss::SSrcs,
    sr::Recs,
    fields::Vector{Symbol},
)
    return [NamedD(grid, sr, fields) for i = 1:ss.n]
end

"""
Print information about `NamedD`.
"""
function Base.show(io::Base.IO, src::NamedD)
    s = isa(src.sr, Srcs) ? "source" : "receiver"
    println(io, "> data with fields ", names(src.d)[1], " at ", src[:n], " ", s)
end

#function Base.isapprox(src1::NamedD, src2::NamedD)
#	vec=([(src1.nss==src2.nss), (src1.fields==src2.fields), (src1.ns==src2.ns),
#       		(isequal(src1.grid, src2.grid)),
#		])
#	return all(vec)
#end

"""
```
rmul!(srcwav[1], b)
```
Scale `srcwav` for first supersource by a scalar `b` overwriting in-place.
"""
function LinearAlgebra.rmul!(dat::NamedD, x::Number)
    for dd in dat.d
        rmul!(dd, x)
    end
end

"""
```
rmul!(srcwav, b)
```
Scale `srcwav` for all supersources by a scalar `b` overwriting in-place.
"""
function LinearAlgebra.rmul!(dat::VNamedD, x::Number)
    for d in dat
        rmul!(d, x)
    end
end

"""
```julia
update!(srcwav[1], [:p, :vx], w)
```
Populate first supersource of `srcwav` by a wavelet vector `w` for `:p` and `:vx` fields.
"""
function update!(d::NamedD, fields::Vector{Symbol}, w::AbstractVector)
    @assert length(w) == length(d.grid)
    @inbounds for f in fields
        @assert f ∈ names(d.d)[1]
        @inbounds for i = 1:d[:n]
            for it = 1:length(d.grid)
                d.d[f][it, i] = w[it]
            end
        end
    end
    return d
end

"""
```julia
update!(srcwav[1], [:p, :vx], w)
```
Populate first supersource of `srcwav` by a wavelets `w` for `:p` and `:vx` fields.
"""
function update!(d::NamedD, fields::Vector{Symbol}, w::AbstractMatrix)
    @assert size(w, 1) == length(d.grid)
    @assert size(w, 2) == d[:n]
    @inbounds for f in fields
        @assert f ∈ names(d.d)[1]
        @inbounds for i = 1:d[:n]
            for it = 1:length(d.grid)
                d.d[f][it, i] = w[it, i]
            end
        end
    end
    return d
end

"""
```julia
update!(srcwav, [:p, :vx], w)
```
Populate `srcwav` by a wavelet vector `w` for `:p` and `:vx` fields.
"""
function update!(dat::VNamedD, fields::Vector{Symbol}, w::AbstractArray)
    for d in dat
        update!(d, fields, w)
    end
end

"""
```julia
isequal(srcwav[1], srcwav[2])
```
Assert if wavelets due to first two supersources are equal.
"""
function Base.isequal(dat1::NamedD, dat2::NamedD)
    return (dat1.grid == dat2.grid) && (dat1.d == dat2.d)
end

"""
```julia
isequal(srcwav1, srcwav2)
```
Assert if `srcwav1` and `srcwav2` are equal.
"""
function Base.isequal(dat1::VNamedD, dat2::VNamedD)
    return all([isequal(dat1[i], dat2[i]) for i = 1:length(dat1)])
end

"""
Return if two `NamedD`'s have same dimensions and bounds.
"""
function issimilar(dat1::NamedD, dat2::NamedD)
    return isequal(dat1.grid, dat2.grid) &&
           isequal(names(dat1.d), names(dat2.d)) &&
           (size(dat1.d) == size(dat2.d))
end

"""
```julia
issimilar(srcwav1, srcwav2)
```
Assert if `srcwav1` and `srcwav2` have same dimensions and fields.
"""
function issimilar(dat1::VNamedD, dat2::VNamedD)
    @assert length(dat1) == length(dat2)
    return all([issimilar(dat1[i], dat2[i]) for i = 1:length(dat1)])
end

function Base.length(data::NamedD)
    return length(data.grid) * data[:n] * length(names(data.d)[1])
end

"""
Return a vec of data object sorted in the order
time, channels
"""
function Base.vec(data::NamedD)
    v = Vector{Float64}()
    for dd in data.d
        append!(v, dd)
    end
    return v
end

"""
```julia
v=vec(srcwav)
Reshape `srcwav` as a one-dimensional column vector.
```
"""
function Base.vec(data::VNamedD)
    v = Vector{Float64}()
    for d in data
        append!(v, vec(d))
    end
    return v
end

function Base.copyto!(d::NamedD, v::AbstractVector{Float64}, i0 = 1)
    #	@assert length(d)==length(v)
    for dd in d.d
        for i = 1:d[:n]
            for it = 1:length(d.grid)
                dd[it, i] = v[i0]
                i0 += 1
            end
        end
    end
    return i0
end

"""
```julia
v=vec(srcwav)
randn!(v)
copyto!(srcwav, v)
```
Copy `v` to `srcwav`. The sizes of the two objects much match.
No memory allocations.
"""
function Base.copyto!(d::VNamedD, v::AbstractVector{Float64})
    i0 = 1
    for i = 1:length(d)
        dd = d[i]
        i0 = copyto!(dd, v, i0)
    end
    return d
end

function Base.copyto!(v::AbstractVector{Float64}, d::NamedD, i0 = 1)
    for dd in d.d
        for i = 1:d[:n]
            for it = 1:length(d.grid)
                v[i0] = dd[it, i]
                i0 += 1
            end
        end
    end
    return i0
end

"""
```julia
copyto!(v, srcwav)
```
Same as `vec(srcwav)`, but in-place operation.
"""
function Base.copyto!(v::AbstractVector{Float64}, d::VNamedD)
    i0 = 1
    for i = 1:length(d)
        dd = d[i]
        i0 = copyto!(v, dd, i0)
    end
    return d
end

function Random.randn!(data::NamedD)
    for dd in data.d
        Random.randn!(dd)
    end
    return data
end

"""
Fill `srcwav` with Random values.
"""
function Random.randn!(data::VNamedD)
    for d in data
        Random.randn!(d)
    end
end

function Base.copyto!(dataout::NamedD, data::NamedD)
    for i in names(dataout.d)[1]
        d1 = dataout.d[i]
        d2 = data.d[i]
        copyto!(d1, d2)
    end
end

"""
```julia
copyto!(srcwav1, srcwav2)
```
Copy `srcwav2` to `srcwav1`, which has same size.
"""
function Base.copyto!(dataout::VNamedD, data::VNamedD)
    for i = 1:length(dataout)
        d1 = dataout[i]
        d2 = data[i]
        copyto!(d1, d2)
    end
end

function Base.iszero(data::NamedD)
    return maximum(broadcast(maximum, abs, data.d)) == 0.0 ? true : false
end

"""
```julia
iszero(srcwav)
```
Returns bool if `srcwav` has all zeros.
"""
function Base.iszero(data::VNamedD)
    return all([iszero(d) for d in data])
end

function LinearAlgebra.dot(data1::NamedD, data2::NamedD)
    dotd = 0.0
    for iff in names(data1.d)[1]
        dd1 = data1.d[iff]
        dd2 = data2.d[iff]
        dotd += LinearAlgebra.dot(dd1, dd2)
    end
    return dotd
end

"""
```julia
dot(srcwav1, srcwav2)
```
Returns dot product. The sizes must match.
"""
function LinearAlgebra.dot(data1::VNamedD, data2::VNamedD)
    dotd = 0.0
    for i = 1:length(data1)
        dd1 = data1[i]
        dd2 = data2[i]
        dotd += LinearAlgebra.dot(dd1, dd2)
    end
    return dotd
end

function Base.fill!(data::NamedD, k::Float64)
    for dd in data.d
        fill!(dd, k)
    end
end

"""
```julia
fill!(srcwav, b)
```
In-place method to fill `srcwav` with scalar `b`.
"""
function Base.fill!(data::VNamedD, k::Float64)
    for d in data
        fill!(d, k)
    end
end

function Base.reverse!(data::NamedD)
    for dd in data.d
        for i = 1:data[:n]
            dv = view(dd, :, i)
            reverse!(dv)
        end
    end
end

"""
```julia
reverse!(srcwav)
```
Perform in-place time-reversal operation for each wavelet in `srcwav`.
"""
function Base.reverse!(data::VNamedD)
    for d in data
        reverse!(d)
    end
end

include("misfits.jl")
include("statistics.jl")

"""
```julia
srcwav_new=interp(srcwav, grid_new)
```
Interpolates `srcwav` onto a new grid.
"""
function interp(data::VNamedD, grid::StepRangeLen, Battrib = :B1)
    nss = length(data)
    dataout = [NamedD(grid, data[iss].sr, names(data[iss].d, 1)) for iss = 1:nss]
    interp_spray!(data, dataout, :interp, Battrib)
    return dataout
end

function interp(data::NamedD, grid::StepRangeLen, Battrib = :B1)
    dataout = NamedD(grid, data.sr, names(data.d, 1))
    interp_spray!(data, dataout, :interp, Battrib)
    return dataout
end

function taper!(data::VNamedD, perc = 0.0; bperc = perc, eperc = perc)
    for d in data
        taper!(d, perc, bperc = bperc, eperc = eperc)
    end
    return data
end

function taper!(data::NamedD, perc = 0.0; bperc = perc, eperc = perc)
    for dd in data
        Utils.taper!(dd, bperc = bperc, eperc = eperc)
    end
    return data
end

"""
```julia
interp_spray!(srcwav_new, srcwav, attrib, Battrib)
```
Method to interpolate `srcwav` onto to a new grid.
Can reduce allocations =========

"""
function interp_spray!(
    data::NamedD,
    dataout::NamedD,
    attrib = :interp,
    Battrib = :B1;
    pa = nothing,
)
    @assert length(data.d) == length(dataout.d)
    xin = data.grid
    xout = dataout.grid
    if (pa === nothing)
        pa = Interpolation.Kernel([xin], [xout], :B1)
    end
    for i = 1:length(data.d)
        dd = data.d[i]
        ddo = dataout.d[i]
        @assert size(dd, 2) == size(ddo, 2)
        for ir = 1:size(dd, 2)
            din = view(dd, :, ir)
            dout = view(ddo, :, ir)
            Interpolation.interp_spray!(din, dout, pa, attrib)
        end
    end
    return dataout
end

function interp_spray!(
    data::VNamedD,
    dataout::VNamedD,
    attrib = :interp,
    Battrib = :B1;
    pa = nothing,
)
    @assert length(data) == length(dataout)
    for i = 1:length(data)
        d = data[i]
        dout = dataout[i]
        interp_spray!(d, dout, attrib, Battrib, pa = pa)
    end
end


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

mutable struct LoadBalancerBackendAddressPropertiesFormat <: SwaggerModel
    virtualNetwork::Any # spec type: Union{ Nothing, SubResource } # spec name: virtualNetwork
    ipAddress::Any # spec type: Union{ Nothing, String } # spec name: ipAddress
    networkInterfaceIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: networkInterfaceIPConfiguration

    function LoadBalancerBackendAddressPropertiesFormat(;virtualNetwork=nothing, ipAddress=nothing, networkInterfaceIPConfiguration=nothing)
        o = new()
        validate_property(LoadBalancerBackendAddressPropertiesFormat, Symbol("virtualNetwork"), virtualNetwork)
        setfield!(o, Symbol("virtualNetwork"), virtualNetwork)
        validate_property(LoadBalancerBackendAddressPropertiesFormat, Symbol("ipAddress"), ipAddress)
        setfield!(o, Symbol("ipAddress"), ipAddress)
        validate_property(LoadBalancerBackendAddressPropertiesFormat, Symbol("networkInterfaceIPConfiguration"), networkInterfaceIPConfiguration)
        setfield!(o, Symbol("networkInterfaceIPConfiguration"), networkInterfaceIPConfiguration)
        o
    end
end # type LoadBalancerBackendAddressPropertiesFormat

const _property_map_LoadBalancerBackendAddressPropertiesFormat = Dict{Symbol,Symbol}(Symbol("virtualNetwork")=>Symbol("virtualNetwork"), Symbol("ipAddress")=>Symbol("ipAddress"), Symbol("networkInterfaceIPConfiguration")=>Symbol("networkInterfaceIPConfiguration"))
const _property_types_LoadBalancerBackendAddressPropertiesFormat = Dict{Symbol,String}(Symbol("virtualNetwork")=>"SubResource", Symbol("ipAddress")=>"String", Symbol("networkInterfaceIPConfiguration")=>"SubResource")
Base.propertynames(::Type{ LoadBalancerBackendAddressPropertiesFormat }) = collect(keys(_property_map_LoadBalancerBackendAddressPropertiesFormat))
Swagger.property_type(::Type{ LoadBalancerBackendAddressPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerBackendAddressPropertiesFormat[name]))}
Swagger.field_name(::Type{ LoadBalancerBackendAddressPropertiesFormat }, property_name::Symbol) =  _property_map_LoadBalancerBackendAddressPropertiesFormat[property_name]

function check_required(o::LoadBalancerBackendAddressPropertiesFormat)
    true
end

function validate_property(::Type{ LoadBalancerBackendAddressPropertiesFormat }, name::Symbol, val)
end


text:
module LAHMC

using LinearAlgebra

export LogLikelihood, GaussianLogLikelihood
include("loglikelihood.jl")

export Sampler, LAHMCSampler, HMCSampler
export set_x0!, sample!
include("sampler.jl")

export mcmc!, mcmc

function mcmc!(s, N, S::Sampler{T}; burnin = 0, thin = 1) where T
    if length(s) > 0
        set_x0!(S, s[end])
    end
    for k in 1:burnin
        sample!(S)
    end
    for k in 1:N
        sample!(S)
        if k % thin == 0
            push!(s, S.x)
        end
    end
    return s
end

function mcmc(x0, N, S0::Sampler{T}; kwargs...) where T
    S = copy(S0)
    set_x0!(S, x0)
    s = Array{T, 1}[]
    mcmc!(s, N, S; kwargs...)
    return s, S
end

end # module


text:
##################################################
# Orbit-averagering of the diffusion coefficients
# In (E,L) space
##################################################
using PolynomialRoots

"""
    radius_s_bounds(E,L)

Computes the bounds of the orbit with action parameters `(E,L)`.
Bounds are in the variable `s`, defined as `s^2 = 1 + x^2`, where `x = r/b`.

Returns `"Unbounded orbit"` for unbounded orbits.
Returns `"Not a possible orbit"` for forbidden parameters.

# Remarks
- Uses `PolynomialRoots`.

# Arguments
- `E::Float64`: Energy parameter.
- `L ::Float64`: Angular momentum parameter.
"""
function radius_s_bounds(E::Float64, L::Float64)
    if (E >= 0.0)
        return "Unbounded orbit"
    elseif (L > Lc(_tE(E)))
        return "Not a possible orbit"
    else
        tE = _tE(E)
        tL = _tL(L)
        if (L != 0.0)
            rts = sort(real(roots([1,-(tE-tL^2/2),-1,tE])))
            return Float64(rts[2]), Float64(rts[3])
        else
            return 1.0, 1/tE
        end
    end
end

"""
    sma_eff(sp,sa)

Computes the effective semi-major axis of an orbit with `s`-bounds `sp` and `sa`.
Bounds are in the variable `s`, defined as `s^2 = 1 + x^2`, where `x = r/b`.

# Arguments
- `sp::Float64`: Pericenter in the radius variable `s`.
- `sa::Float64`: Apocenter in the radius variable `s`.
"""
function sma_eff(sp::Float64, sa::Float64)
    return (sa+sp)/2
end

"""
    ecc_eff(sp,sa)

Computes the effective eccentricity of an orbit with `s`-bounds `sp` and `sa`.
Bounds are in the variable `s`, defined as `s^2 = 1 + x^2`, where `x = r/b`.

# Arguments
- `sp::Float64`: Pericenter in the radius variable `s`.
- `sa::Float64`: Apocenter in the radius variable `s`.
"""
function ecc_eff(sp::Float64, sa::Float64)
    return (sa-sp)/(sa+sp)
end

"""
    tr_eff_anomaly(u,a,ecc)

Computes the non-dimensional radius of a star with effective anomaly `u`, semi-major axis `a` and eccentricity `ecc`.

# Arguments
- `u::Float64`: Effective anomaly parameter.
- `a::Float64`: Effective semi-major axis.
- `ecc::Float64`: Effective eccentricity.
"""
function tr_eff_anomaly(u::Float64, a::Float64, ecc::Float64)
    s = a*(1+ecc*u*(1.5-u^2/2))

    return sqrt(abs(s^2-1)) # use absolute value for when s->1
end

"""
    r_eff_anomaly(u,a,ecc)

Computes the radius of a star with effective anomaly `u`, semi-major axis `a` and eccentricity `ecc`.

# Arguments
- `u::Float64`: Effective anomaly parameter.
- `a::Float64`: Effective semi-major axis.
- `ecc::Float64`: Effective eccentricity.
"""
function r_eff_anomaly(u::Float64, a::Float64, ecc::Float64)
    return _b*tr_eff_anomaly(u,a,ecc)
end

"""
    dth1du(u,sp,sa)

Non-dimensional Jacobian of the transformation from angles `θ` to effective anomaly `u`.

# Arguments
- `u::Float64`: Effective anomaly.
- `sp::Float64`: Pericenter in the variable `s`.
- `sa::Float64`: Apocenter in the variable `s`.
"""
function dth1du(u::Float64, sp::Float64, sa::Float64)
    num = sa*(-2+u)*(1+u)^2-sp*(2-3*u+u^3)
    denumSq = ((sa^2*sp*(-2+u)*(1+u)^2-sp*(6-3*u+u^3)
               +sa*(-6-3*u+u^3-sp^2*(2-3*u+u^3)))/
               (sa*sp*(sa+sp)*(sa*(-2+u)*(1+u)^2-sp*(2-3*u+u^3))))
    return -3/(4*sqrt(2))*1/sqrt(4-u^2)*num/sqrt(denumSq)
end

"""
    halfperiod(E,L,[nbAvr=nbAvr_default])

Computes half the radial period of the orbit with parameters `(E,L)`.

# Remarks
- Well-defined everywhere by using an effective anomaly.

# Arguments
- `E::Float64`: Energy parameter.
- `L::Float64`: Angular momentum parameter.
- `nbAvr::Int64`: Orbit average sampling number.
"""
function halfperiod(E::Float64, L::Float64, nbAvr::Int64=nbAvr_default)
    sp, sa = radius_s_bounds(E,L)
    halfperiod = 0.0

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        jac_loc = dth1du(uloc,sp,sa)
        halfperiod += jac_loc
    end

    halfperiod *= (1/_Omega0)*(2/nbAvr)

    return halfperiod
end

"""
    avr_action_coefficients(E,L,m_field,[nbK=nbK_default,nbAvr=nbAvr_default,m_test=m_field])

Orbit-averaged diffusion coefficients in energy/angular momentum space `(E,L)`.
Returns a 5-uple `(avrDE, avrDL, avrDEE, avrDEL, avrDLL)`.

# Arguments
- `E::Float64`: Energy parameter.
- `L::Float64`: Angular momentum parameter.
- `m_field::Float64`: Mass of field star.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `nbAvr::Int64`: Sampling number of orbit average.
- `m_test::Float64`: Mass of test star.
"""
function avr_orbit_coefficients(E::Float64, L::Float64, m_field::Float64,
                                nbK::Int64=nbK_default, nbAvr::Int64=nbAvr_default,
                                m_test::Float64=m_field)
    sp, sa = radius_s_bounds(E,L)
    a = sma_eff(sp,sa)
    ecc = ecc_eff(sp,sa)

    avrDE = 0.0
    avrDL = 0.0
    avrDEE = 0.0
    avrDEL = 0.0
    avrDLL = 0.0
    halfperiod = 0.0

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        rloc = r_eff_anomaly(uloc,a,ecc)
        jac_loc = dth1du(uloc,sp,sa)
        dEloc, dE2loc, dLloc, dEdLloc, dL2loc = localOrbitChangeReg(rloc,E,L,m_field,nbK,m_test)

        avrDE += jac_loc*dEloc
        avrDL += jac_loc*dLloc
        avrDEE += jac_loc*dE2loc
        avrDEL += jac_loc*dEdLloc
        avrDLL += jac_loc*dL2loc
        halfperiod += jac_loc
    end
    avrDE /= halfperiod
    avrDL /= halfperiod
    avrDEE /= halfperiod
    avrDEL /= halfperiod
    avrDLL /= halfperiod

    return avrDE, avrDL, avrDEE, avrDEL, avrDLL
end

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
# Functions that apply the entropy and mutual information formulae

function applyentropyformula(probabilities, base)
	return -sum([p == 0 ? 0 : p .* log(p) for p in probabilities]) / log(base)
end

function applymutualinformationformula(entropyX, entropyY, entropyXY)
	return entropyX + entropyY - entropyXY
end

function applyconditionalmutualinformationformula(entropyZ, entropyXZ, entropyYZ, entropyXYZ)
	return entropyXZ + entropyYZ - entropyXYZ - entropyZ
end

function applytotalcorrelationformula(entropyX, entropyY, entropyZ, entropyXYZ)
	return entropyX + entropyY + entropyZ - entropyXYZ
end

max_stars_repo_name: Tchanders/ShrinkageEstimators

max_issues_repo_name: Tchanders/ShrinkageEstimators

max_forks_repo_name: Tchanders/ShrinkageEstimators


text:
@testset "models.jl" begin
    @testset "PoissonExpModel" begin
        @testset "Zero firing rate, infinite time constant" begin
            x = ones(10)
            model = PoissonExpModel(0, Inf, 10)
            update!(x, model, 0.01)

            @test all(x .== 1)
        end

        @testset "Zero firing rate, unit time constant" begin
            x = ones(10)
            model = PoissonExpModel(0, 1, 10)
            update!(x, model, 0.01)

            @test all(x .== exp(-0.01))
        end

        @testset "Unit firing rate, infinite time constant" begin
            x = ones(10)
            model = PoissonExpModel(1, Inf, 10)
            update!(x, model, 1.)

            @test all(x .>= 1)
        end
    end#PoissonExpModel

    @testset "BlockPoissonExpModel" begin
        @testset "" begin
            x = zeros(32)
            model = BlockPoissonExpModel(1, Inf, 8, 4, 1)

            @test_throws MethodError update!(x, model, 0.01)

            update!(x, model, 0.01, 0.)
            @test all(x[5:32] .== 0)

            update!(x, model, 0.01, 5.1)
            @test all(x[5:20] .== 0)
            @test all(x[25:32] .== 0)
        end
    end#BlockPoissonExpModel

    @testset "OUModel" begin
        @testset "Unit time constant, zero variance" begin
            w = ones(10)
            model = OUModel(1, 0)
            update!(w, model, 0.01)

            @test all(w .== 0.99)
        end

        @testset "unit time constant, unit variance" begin
            w = zeros(10)
            model = OUModel(1, 1)
            update!(w, model, 0.01)

            @test all(w .!= 0)
        end
    end#OUModel

    @testset "State" begin
        w = ones(10)
        x = ones(10)
        state = State(w, x)

        Imodel = PoissonExpModel(1, 0.1, 10)
        Smodel = OUModel(1, 1)

        wcopy = copy(w)
        update!(state, Imodel, 0.01)
        @test all(state.w .== wcopy)

        xcopy = copy(state.x)
        update!(state, Smodel, 0.01)
        @test all(state.x .== xcopy)
    end#State

    @testset "ExpGainModel" begin
        @testset "zero firing rate" begin
            w = rand(10)
            x = rand(10)

            wcopy = copy(w)
            xcopy = copy(x)

            state = State(w, x)
            model = ExpGainModel(0, 1)
            obs = update!(state, model, 0.01)

            @test all(state.w .== wcopy)
            @test all(state.x .== xcopy)
            @test all(obs.x   .== xcopy)
            @test all(obs.y   .== 0)
        end

        @testset "high firing rate" begin
            w = rand(10)
            x = rand(10)

            wcopy = copy(w)
            xcopy = copy(x)

            state = State(w, x)
            model = ExpGainModel(100, 1)
            obs = update!(state, model, 1.)

            @test all(state.w .== wcopy)
            @test all(state.x .== xcopy)
            @test all(obs.x   .== xcopy)
            @test all(obs.y   .> 0)
        end
    end#ExpGainModel

    @testset "AdaptiveExpGainModel" begin
        @testset "Infinite time constant, α = 1" begin
            w = rand(10)
            x = rand(10)
            η = [0.]

            wcopy = copy(w)
            xcopy = copy(x)
            ηcopy = copy(η)

            state = State(w, x, η)
            model = AdaptiveExpGainModel(1, 1, 1, Inf)
            obs = update!(state, model, 0.01)

            @test all(state.w .== wcopy)
            @test all(state.x .== xcopy)
            @test all(state.η .== obs.y)
            @test all(obs.x   .== xcopy)
        end

        @testset "Unit time constant, α = 0" begin
            w = rand(10)
            x = rand(10)
            η = [1.]

            wcopy = copy(w)
            xcopy = copy(x)

            state = State(w, x, η)
            model = AdaptiveExpGainModel(1, 1, 0, 1)
            obs = update!(state, model, 1.)

            @test all(state.w .== wcopy)
            @test all(state.x .== xcopy)
            @test all(state.η .== exp(-1))
            @test all(obs.x   .== xcopy)
        end
    end#AdaptiveExpGainModel

    if RUN_BENCHMARKS
    @testset "Benchmark PoissonExpModel" begin
        w = rand(1024)
        x = rand(1024)
        state = State(w, x)
        model = PoissonExpModel(1, 0.1, 10)

        println("")
        println("Benchmarking one update step for PoissonExpModel")
        display(@benchmark update!($state, $model, 0.01))
        println("")
        println("")
    end#Benchmark PoissonExpModel

    @testset "Benchmark BlockPoissonExpModel" begin
        w = rand(1024)
        x = rand(1024)
        state = State(w, x)
        model = BlockPoissonExpModel(1, 0.1, 128, 8, 1.)

        println("")
        println("Benchmarking one update step for BlockPoissonExpModel")
        display(@benchmark update!($state, $model, 0.01, 5.2))
        println("")
        println("")
    end#Benchmark BlockPoissonExpModel

    @testset "Benchmark OUModel" begin
        w = rand(1024)
        x = rand(1024)
        state = State(w, x)
        model = OUModel(1, 0.1)

        println("")
        println("Benchmarking one update step for OUModel")
        display(@benchmark update!($state, $model, 0.01))
        println("")
        println("")
    end#Benchmark OUModel

    @testset "Benchmark ExpGainModel" begin
        w = rand(1024)
        x = rand(1024)
        state = State(w, x)
        model = ExpGainModel(1, 0.1)

        println("")
        println("Benchmarking one update step for ExpGainModel")
        display(@benchmark update!($state, $model, 0.01))
        println("")
        println("")
    end#Benchmark ExpGainModel

    @testset "Benchmark AdaptiveExpGainModel" begin
        w = randn(1024)
        x = randn(1024)
        state = State(w, x)
        model = AdaptiveExpGainModel(1, 0.01, -1, 10.)

        println("")
        println("Benchmarking one update step for AdaptiveExpGainModel")
        display(@benchmark update!($state, $model, 0.01))
        println("")
        println("")
    end#Benchmark AdaptiveExpGainModel
    end#if
end#models.jl


text:
module Goxel # VoxelPlots

include("Goxel/mesh.jl")

end # module VoxelPlots.Goxel


text:
using Oceananigans.Fields: FunctionField, location
using Oceananigans.TurbulenceClosures: implicit_step!
using Oceananigans.Architectures: device_event
using Oceananigans.Utils: @apply_regionally, apply_regionally!

mutable struct QuasiAdamsBashforth2TimeStepper{FT, GT, IT} <: AbstractTimeStepper
                  χ :: FT
        previous_Δt :: FT
                 Gⁿ :: GT
                 G⁻ :: GT
    implicit_solver :: IT
end

"""
    QuasiAdamsBashforth2TimeStepper(grid, tracers,
                                    χ = 0.1;
                                    implicit_solver = nothing,
                                    Gⁿ = TendencyFields(grid, tracers),
                                    G⁻ = TendencyFields(grid, tracers))

Return a 2nd-order quasi Adams-Bashforth (AB2) time stepper (`QuasiAdamsBashforth2TimeStepper`)
on `grid`, with `tracers`, and AB2 parameter `χ`. The tendency fields `Gⁿ` and `G⁻` can be
specified via  optional `kwargs`.

The 2nd-order quasi Adams-Bashforth timestepper steps forward the state `Uⁿ` by `Δt` via

```julia
Uⁿ⁺¹ = Uⁿ + Δt * [(3/2 + χ) * Gⁿ - (1/2 + χ) * Gⁿ⁻¹]
```

where `Uⁿ` is the state at the ``n``-th timestep, `Gⁿ` is the tendency
at the ``n``-th timestep, and `Gⁿ⁻¹` is the tendency at the previous
timestep (`G⁻`).

!!! note "First timestep"
    For the first timestep, since there are no saved tendencies from the previous timestep,
    the `QuasiAdamsBashforth2TimeStepper` performs an Euler timestep:

    ```julia
    Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ
    ```
"""
function QuasiAdamsBashforth2TimeStepper(grid, tracers,
                                         χ = 0.1;
                                         implicit_solver::IT = nothing,
                                         Gⁿ = TendencyFields(grid, tracers),
                                         G⁻ = TendencyFields(grid, tracers)) where IT

    FT = eltype(grid)
    GT = typeof(Gⁿ)

    return QuasiAdamsBashforth2TimeStepper{FT, GT, IT}(χ, Inf, Gⁿ, G⁻, implicit_solver)
end

function reset!(timestepper::QuasiAdamsBashforth2TimeStepper)
    timestepper.previous_Δt = Inf
    return nothing
end

#####
##### Time steppping
#####

"""
    time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and
pressure-correction substep. Setting `euler=true` will take a forward Euler time step.
"""
function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false)
    Δt == 0 && @warn "Δt == 0 may cause model blowup!"

    # Shenanigans for properly starting the AB2 loop with an Euler step
    euler = euler || (Δt != model.timestepper.previous_Δt)

    χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ)

    if euler
        @debug "Taking a forward Euler step."
        # Ensure zeroing out all previous tendency fields to avoid errors in
        # case G⁻ includes NaNs. See https://github.com/CliMA/Oceananigans.jl/issues/2259
        for field in model.timestepper.G⁻
            !isnothing(field) && @apply_regionally fill!(field, 0)
        end
    end

    model.timestepper.previous_Δt = Δt

    # Be paranoid and update state at iteration 0
    model.clock.iteration == 0 && update_state!(model)

    @apply_regionally calculate_tendencies!(model)

    ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities.
    calculate_pressure_correction!(model, Δt)

    @apply_regionally correct_velocties_and_store_tendecies!(model, Δt)

    tick!(model.clock, Δt)
    update_state!(model)
    update_particle_properties!(model, Δt)

    return nothing
end

function correct_velocties_and_store_tendecies!(model, Δt)
    pressure_correct_velocities!(model, Δt)
    store_tendencies!(model)
end

#####
##### Time stepping in each step
#####

""" Generic implementation. """
function ab2_step!(model, Δt, χ)

    workgroup, worksize = work_layout(model.grid, :xyz)
    arch = model.architecture
    barrier = device_event(arch)
    step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize)
    model_fields = prognostic_fields(model)
    events = []

    for (i, field) in enumerate(model_fields)

        field_event = step_field_kernel!(field, Δt, χ,
                                         model.timestepper.Gⁿ[i],
                                         model.timestepper.G⁻[i],
                                         dependencies = device_event(arch))

        push!(events, field_event)

        # TODO: function tracer_index(model, field_index) = field_index - 3, etc...
        tracer_index = Val(i - 3) # assumption

        implicit_step!(field,
                       model.timestepper.implicit_solver,
                       model.closure,
                       model.diffusivity_fields,
                       tracer_index,
                       model.clock,
                       Δt,
                       dependencies = field_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step velocity fields via the 2nd-order quasi Adams-Bashforth method

    `U^{n+1} = U^n + Δt ((3/2 + χ) * G^{n} - (1/2 + χ) G^{n-1})`

"""
@kernel function ab2_step_field!(u, Δt, χ, Gⁿ, G⁻)
    i, j, k = @index(Global, NTuple)

    T = eltype(u)
    one_point_five = convert(T, 1.5)
    oh_point_five = convert(T, 0.5)

    @inbounds u[i, j, k] += Δt * ((one_point_five + χ) * Gⁿ[i, j, k] - (oh_point_five + χ) * G⁻[i, j, k])
end

@kernel ab2_step_field!(::FunctionField, args...) = nothing