text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
A = rand(0:9, 3, 4)

# %%
display(A)

# %%
print("A = "); display(A)

# %%
print("A = "); flush(stdout); display(A)

# %%
print("A = "); show(stdout, MIME"text/plain"(), A)

# %%
print("A = "); show(stdout, MIME("text/plain"), A)

# %%
?MIME

# %%

max_stars_repo_path: 0013/show MIME example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/show MIME example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/show MIME example.jl

max_forks_repo_name: genkuroki/public


text:
__precompile__(true)
module JuliaPackageTemplate

# Usings
using Reexport

# Base Module Includes
include("submodule.jl")
include("rubber_ducks.jl")

# Export Module Contents
@reexport using JuliaPackageTemplate.YourSubmodule
@reexport using JuliaPackageTemplate.RubberDucks

end # module

ext: jl


text:
using CairoMakie
using naumer_Dimensionality_2022
using LinearAlgebra: tr, svd, dot, I
using CSV
using Random: MersenneTwister

################################################
########### Plot Parameters ####################
################################################
mrng = MersenneTwister(1234)

max_steps = 1000
state     = 2*randn(mrng,3)#5*ones(3)
state_cp  = copy(state)
crlb      = I(3)
τ = 0.01
γ = 0.9
actionSpace = [randn(mrng,3) for i in 1:50]
for i in 1:50
    actionSpace[i] /= sqrt(sum(abs2,actionSpace[i]))
end

λ = 1
psdSampleCount = [50]
trajectorySampleCount = 100
timestepSampleCount = 1
σ² = 1
max_iterations = 1000
d_max = 2

T = 1

#################################################
################# Data Generation ###############
#################################################

system = LorenzSystem(10.0, 28.0, 8.0/3.0)

optimal_sampling_trace      = zeros(max_steps)
random_sampling_trace       = zeros(max_steps)
approx_sampling_trace       = zeros(max_steps)
optimal_sampling_trace_ekf  = zeros(max_steps)
approx_sampling_trace_ekf   = zeros(max_steps)

optimal_sampling_trace[1]   = tr(crlb)
random_sampling_trace[1]    = tr(crlb)
random_sampling_trace[1]    = tr(crlb)

@info("Approximating Value Function")

values, psdSamples, stateSamples = ValueFunctionApproximation_LocalAverage_precompute(  mrng,
                                                                                system,
                                                                                τ,
                                                                                γ,
                                                                                actionSpace,
                                                                                λ,
                                                                                psdSampleCount[1],
                                                                                trajectorySampleCount,
                                                                                timestepSampleCount,
                                                                                σ²,
                                                                                max_iterations,
                                                                                d_max)

@info("Solving Dynamic Programming Optimal Sampling")
NN_Policy_tuple = (system, σ², τ, values, psdSamples, stateSamples, actionSpace, d_max)

for i in 2:max_steps
    global state
    global crlb
    global NN_Policy_tuple
    action, index, state, crlb = LocalAverage_OptimalPolicy(    state,
                                                                crlb,
                                                                NN_Policy_tuple...)

    optimal_sampling_trace[i] = tr(crlb)
end

@info("Solving Optimal Sampling Dynamic Programming Problem - EKF")

NN_Policy_tuple = (system, σ², τ, values, psdSamples, stateSamples, actionSpace, d_max)
beliefState = zeros(3)
beliefCRLB  = I(3)
crlb        = I(3)
state       = copy(state_cp)

for i in 2:max_steps
    global state
    global crlb
    global beliefState
    global beliefCRLB

    # Choose Action
    action, ~, ~, ~ = LocalAverage_OptimalPolicy(   beliefState,
                                                    beliefCRLB,
                                                    NN_Policy_tuple...)

    # Update True System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Update Belief System
    observation = dot(action, state) + sqrt(σ²) * randn(mrng)
    jacobian = flowJacobian(beliefState, τ, system)
    beliefCRLB = updateCRLB(beliefCRLB, action, jacobian, σ²)
    beliefState = flow(beliefState, τ, system)
    beliefState = stateupdate_EKF(beliefState, beliefCRLB, observation, action, σ²)

    # Record cost
    optimal_sampling_trace_ekf[i] = tr(crlb)
end

@info("Solving Optimal Sampling 1D Approximation Problem")

crlb = I(3)
state = copy(state_cp)

for i in 2:max_steps
    global state
    global crlb
    global T
    global σ²
    global system

    # Decide action
    ~, action = optimalaction_1DApprox(actionSpace, system, state, crlb, σ², T)

    # Update System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Record Cost
    approx_sampling_trace[i] = tr(crlb)
end

@info("Solving Optimal Sampling 1D Approximation Problem - EKF")

crlb = I(3)
state = copy(state_cp)
beliefState = zeros(3)
beliefCRLB = I(3)

for i in 2:max_steps
    global state
    global crlb
    global beliefState
    global beliefCRLB

    # Decide action
    ~, action = optimalaction_1DApprox(actionSpace, system, beliefState, beliefCRLB, σ², T)

    # Update True System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Update Belief System
    observation = dot(action, state) + sqrt(σ²) * randn(mrng)
    jacobian = flowJacobian(beliefState, τ, system)
    beliefCRLB = updateCRLB(beliefCRLB, action, jacobian, σ²)
    beliefState = flow(beliefState, τ, system)
    beliefState = stateupdate_EKF(beliefState, beliefCRLB, observation, action, σ²)

    # Record Cost
    approx_sampling_trace_ekf[i] = tr(crlb)
end

@info("Solving Random Sampling")

crlb = I(3)
state = copy(state_cp)

for i in 2:max_steps
    global state
    global crlb
    action = actionSpace[rand(mrng, 1:length(actionSpace))]

    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)

    state = flow(state, τ, system)
    random_sampling_trace[i] = tr(crlb)
end

#################################################
################ Save Data CSV ##################
#################################################

CSV.write("LorenzSampling.csv",(random_sampling_trace       = random_sampling_trace,
                                optimal_sampling_trace      = optimal_sampling_trace,
                                optimal_sampling_trace_ekf  = optimal_sampling_trace_ekf,
                                approx_sampling_trace       = approx_sampling_trace,
                                approx_sampling_trace_ekf   = approx_sampling_trace_ekf,
                                samples = Array(1:max_steps)))


text:
using OneDimensionalNodes
using FastGaussQuadrature
using Test

include("operators.jl")
include("quadratures.jl")


text:
using FirstOrderLogic
if VERSION ≥ v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

# test all boolean operations
@testset "¬" begin
    @test ¬true  ≡ false
    @test ¬false ≡ true
end

@testset "∧" begin
    @test true  ∧ true  ≡ true
    @test true  ∧ false ≡ false
    @test false ∧ true  ≡ false
    @test false ∧ false ≡ false
end

@testset "∨" begin
    @test true  ∨ true  ≡ true
    @test true  ∨ false ≡ true
    @test false ∨ true  ≡ true
    @test false ∨ false ≡ false
end

@testset "⊻" begin
    @test true  ⊻ true  ≡ false
    @test true  ⊻ false ≡ true
    @test false ⊻ true  ≡ true
    @test false ⊻ false ≡ false
end

@testset "⇒" begin
    @test true  ⇒ true  ≡ true
    @test true  ⇒ false ≡ false
    @test false ⇒ true  ≡ true
    @test false ⇒ false ≡ true
end

@testset "⇔" begin
    @test true  ⇔ true  ≡ true
    @test true  ⇔ false ≡ false
    @test false ⇔ true  ≡ false
    @test false ⇔ false ≡ true
end

max_stars_repo_name: TotalVerb/FirstOrderLogic

max_issues_repo_name: TotalVerb/FirstOrderLogic

max_forks_repo_name: TotalVerb/FirstOrderLogic


text:
# these are in style not recommended by Julia, but they
# should not generate lint messages.
s = """
function f(Y::Array)
  X = zeros(100, 1)
  for i = 1:100
    X[i] = Y[i]
  end
  X
end

function f(Y::Array)
  function f2(Y=Y)
    Y*2
  end
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)


text:
module TestDataFrameRow
    using Base.Test
    using DataFrames, Compat

    df = DataFrame(a=@data([1,   2,   3,   1,   2,   2 ]),
                   b=@data([2.0, NA,  1.2, 2.0, NA,  NA]),
                   c=@data(["A", "B", "C", "A", "B", NA]),
                   d=PooledDataArray(
                     @data([:A,  NA,  :C,  :A,  NA,  :C])))
    df2 = DataFrame(a = @data([1, 2, 3]))

    #
    # Equality
    #
    @test_throws ArgumentError isequal(DataFrameRow(df, 1), DataFrameRow(df2, 1))
    @test !isequal(DataFrameRow(df, 1), DataFrameRow(df, 2))
    @test !isequal(DataFrameRow(df, 1), DataFrameRow(df, 3))
    @test isequal(DataFrameRow(df, 1), DataFrameRow(df, 4))
    @test isequal(DataFrameRow(df, 2), DataFrameRow(df, 5))
    @test !isequal(DataFrameRow(df, 2), DataFrameRow(df, 6))

    # hashing
    @test !isequal(hash(DataFrameRow(df, 1)), hash(DataFrameRow(df, 2)))
    @test !isequal(hash(DataFrameRow(df, 1)), hash(DataFrameRow(df, 3)))
    @test isequal(hash(DataFrameRow(df, 1)), hash(DataFrameRow(df, 4)))
    @test isequal(hash(DataFrameRow(df, 2)), hash(DataFrameRow(df, 5)))
    @test !isequal(hash(DataFrameRow(df, 2)), hash(DataFrameRow(df, 6)))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 38253cd6-f36c-11ea-397f-3d3dce9bac78
using Pkg, DrWatson

# ╔═╡ 38257d6a-f36c-11ea-3d0a-c19b131da11c
begin
	@quickactivate "SR2TuringPluto"
	using StatisticalRethinking
end

# ╔═╡ 80e5d0da-f36b-11ea-19fa-d9c7e930c6b1
md"## Clip-03-11-16t.jl"

# ╔═╡ 3825f68c-f36c-11ea-07a7-470cc146cc6e
md"### snippet 3.11"

# ╔═╡ 382ee7ec-f36c-11ea-198a-eb0a137d43d1
begin
	p_grid = range(0, step=0.001, stop=1)
	prior = ones(length(p_grid))
	likelihood = [pdf(Binomial(3, p), 3) for p in p_grid]
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
end;

# ╔═╡ 382f603a-f36c-11ea-1f05-7716e241d580
md"##### Draw 10000 samples from this posterior distribution."

# ╔═╡ 383b80ce-f36c-11ea-23a0-97b799c3315a
begin
	N = 10000
	samples = sample(p_grid, Weights(posterior), N);
end;

# ╔═╡ 383c1926-f36c-11ea-3f27-5562a255ce30
md"### snippet 3.13"

# ╔═╡ 38455cac-f36c-11ea-3634-61af92f66f5c
hpdi(samples, alpha=0.11)

# ╔═╡ 3846bb8a-f36c-11ea-086f-23705940b3d4
md"### snippet 3.14"

# ╔═╡ 38517d70-f36c-11ea-1dd8-83f037833053
mode(samples)

# ╔═╡ 3852341a-f36c-11ea-20ed-75fbf044d1d2
md"### snippet 3.15"

# ╔═╡ 38581a4c-f36c-11ea-381e-8dc2daabea5a
mean(samples)

# ╔═╡ 3863c502-f36c-11ea-1f72-ed1eb6c862fb
md"### snippet 3.16"

# ╔═╡ 386500de-f36c-11ea-307d-f3d81413fa54
median(samples)

# ╔═╡ 386c6f04-f36c-11ea-38e1-ebc6cdda48ef
md"##### Plot density."

# ╔═╡ 387396e4-f36c-11ea-3f58-ed848b00749a
begin
	density(samples, lab="density")
	vline!(hpdi(samples, alpha=0.5), line=:dash, lab="hpdi")
	vline!(quantile(samples, [0.25, 0.75]), line=:dash, lab="quantile (pi)")
end

# ╔═╡ 387a79fa-f36c-11ea-3c91-379b9c7b4ab3
md"## End of clip-03-11-16t.jl"

# ╔═╡ Cell order:
# ╟─80e5d0da-f36b-11ea-19fa-d9c7e930c6b1
# ╠═38253cd6-f36c-11ea-397f-3d3dce9bac78
# ╠═38257d6a-f36c-11ea-3d0a-c19b131da11c
# ╟─3825f68c-f36c-11ea-07a7-470cc146cc6e
# ╠═382ee7ec-f36c-11ea-198a-eb0a137d43d1
# ╟─382f603a-f36c-11ea-1f05-7716e241d580
# ╠═383b80ce-f36c-11ea-23a0-97b799c3315a
# ╟─383c1926-f36c-11ea-3f27-5562a255ce30
# ╠═38455cac-f36c-11ea-3634-61af92f66f5c
# ╟─3846bb8a-f36c-11ea-086f-23705940b3d4
# ╠═38517d70-f36c-11ea-1dd8-83f037833053
# ╟─3852341a-f36c-11ea-20ed-75fbf044d1d2
# ╠═38581a4c-f36c-11ea-381e-8dc2daabea5a
# ╟─3863c502-f36c-11ea-1f72-ed1eb6c862fb
# ╠═386500de-f36c-11ea-307d-f3d81413fa54
# ╟─386c6f04-f36c-11ea-38e1-ebc6cdda48ef
# ╠═387396e4-f36c-11ea-3f58-ed848b00749a
# ╟─387a79fa-f36c-11ea-3c91-379b9c7b4ab3


text:
using Random

function RLCore.Experiment(::Val{:JuliaRL}, ::Val{:Minimax}, ::Val{:OpenSpiel}, game;)
    env = OpenSpielEnv(string(game))
    agents = (
        Agent(policy = MinimaxPolicy(), role = 0),
        Agent(policy = MinimaxPolicy(), role = 1),
    )
    hooks = (TotalRewardPerEpisode(), TotalRewardPerEpisode())
    description = """
      # Play `$game` in OpenSpiel with Minimax
      """
    Experiment(agents, env, StopAfterEpisode(1), hooks, description)
end

function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:TabularCFR},
    ::Val{:OpenSpiel},
    game;
    n_iter = 300,
    seed = 123,
)
    env = OpenSpielEnv(
        game;
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    rng = StableRNG(seed)
    π = TabularCFRPolicy(; rng = rng)

    description = """
      # Play `$game` in OpenSpiel with TabularCFRPolicy
      """
    Experiment(π, env, StopAfterStep(300), EmptyHook(), description)
end

function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:DeepCFR},
    ::Val{:OpenSpiel},
    game;
    n_iter = 300,
    seed = 123,
)
    @assert game == "leduc_poker"
    env = OpenSpielEnv(
        "leduc_poker";
        default_state_style = RLBase.Information{Array}(),
        is_chance_agent_required = true,
    )

    #= using CxxWrap =#

    p = DeepCFR(
        Π = NeuralNetworkApproximator(
            model = Chain(Dense(30, 64, relu), Dense(64, 64, relu), Dense(64, 3)) |> gpu,
            optimizer = ADAM(0.001),
        ),
        MΠ = ReservoirTrajectory(
            3_000_000,
            :I => Any, #= CxxWrap.StdLib.StdVectorAllocated{Float64}, =#
            :t => Int,
            :σ => Vector{Float32},
            :m => Vector{Bool},
        ),
        V = Dict(
            p => NeuralNetworkApproximator(
                model = Chain(Dense(30, 64, relu), Dense(64, 64, relu), Dense(64, 3)) |> gpu,
                optimizer = ADAM(0.001),
            ) for p in get_players(env) if p != get_chance_player(env)
        ),
        MV = Dict(
            p => ReservoirTrajectory(
                3_000_000,
                :I => Any,#=CxxWrap.StdLib.StdVectorAllocated{Float64}=#
                :t => Int,
                :r̃ => Vector{Float32},
                :m => Vector{Bool},
            ) for p in get_players(env) if p != get_chance_player(env)
        ),
        K = 1500,
        n_training_steps_V = 750,
        n_training_steps_Π = 2000,
        batch_size_V = 2048,
        batch_size_Π = 2048,
        initializer = glorot_normal(CUDA.CURAND.default_rng()),
    )
    # nash_conv ≈ 0.23
    Experiment(p, env, StopAfterStep(500), EmptyHook(), "# run DeepcCFR on leduc_poker")
end


text:
function branching_sim(pomdp::POMDP, policy::Policy, b::ScenarioBelief, steps::Integer, fval)
    S = statetype(pomdp)
    O = obstype(pomdp)
    odict = Dict{O, Vector{Pair{Int, S}}}()

    if steps <= 0
        return length(b.scenarios)*fval(pomdp, b)
    end

    a = action(policy, b)

    r_sum = 0.0
    for (k, s) in b.scenarios
        if !isterminal(pomdp, s)
            rng = get_rng(b.random_source, k, b.depth)
            sp, o, r = gen(DDNOut(:sp, :o, :r), pomdp, s, a, rng)

            if haskey(odict, o)
                push!(odict[o], k=>sp)
            else
                odict[o] = [k=>sp]
            end

            r_sum += r
        end
    end

    next_r = 0.0
    for (o, scenarios) in odict
        bp = ScenarioBelief(scenarios, b.random_source, b.depth+1, o)
        if length(scenarios) == 1
            next_r += rollout(pomdp, policy, bp, steps-1, fval)
        else
            next_r += branching_sim(pomdp, policy, bp, steps-1, fval)
        end
    end

    return r_sum + discount(pomdp)*next_r
end

# once there is only one scenario left, just run a rollout
function rollout(pomdp::POMDP, policy::Policy, b0::ScenarioBelief, steps::Integer, fval)
    @assert length(b0.scenarios) == 1
    disc = 1.0
    r_total = 0.0
    scenario_mem = copy(b0.scenarios)
    (k, s) = first(b0.scenarios)
    b = ScenarioBelief(scenario_mem, b0.random_source, b0.depth, b0._obs)

    while !isterminal(pomdp, s) && steps > 0
        a = action(policy, b)

        rng = get_rng(b.random_source, k, b.depth)
        sp, o, r = gen(DDNOut(:sp, :o, :r), pomdp, s, a, rng)

        r_total += disc*r

        s = sp
        scenario_mem[1] = k=>s
        b = ScenarioBelief(scenario_mem, b.random_source, b.depth+1, o)

        disc *= discount(pomdp)
        steps -= 1
    end

    if steps == 0 && !isterminal(pomdp, s)
        r_total += disc*fval(pomdp, b)
    end

    return r_total
end


text:
# wait for ReverseDiffSource to be updated to v0.5
# import ReverseDiffSource
# import ReverseDiffSource:rdiff

using Calculus: differentiate

export AnalyticPotential

import FunctionWrappers
import FunctionWrappers: FunctionWrapper
typealias F64fun FunctionWrapper{Float64, Tuple{Float64}}

# ===========================================================================
#     macro that takes an expression, differentiates it
#     and returns anonymous functions for the derivatives
# ===========================================================================

"""
`diff2(ex::Expr)` :

takes an expression, differentiates it twice,
     and returns anonymous functions for the derivatives
"""
function diff2(ex::Expr, sym=:r)
   ex_d = differentiate(ex, sym)
   ex_dd = differentiate(ex_d, sym)
   return eval( :( (r->$ex, r->$ex_d, r->$ex_dd)) )
end

# ================== Analytical Potentials ==========================
#
# TODO: this construction should not be restricted to pair potentials
#       but need a better model

# documentation attached below
@pot type AnalyticPotential <: PairPotential
   f::F64fun
   f_d::F64fun
   f_dd::F64fun
   id::AbstractString
   cutoff::Float64
end

"""
`type AnalyticPotential <: PairPotential`

### Usage:
```julia
lj = AnalyticPotential(:(r^(-12) - 2.0*r^(-6)), "LennardJones")
println(lj)   # will output `LennardJones`
A = 4.0; r0 = 1.234
morse = AnalyticPotential("exp(-2.0*\$A*(r/\$r0-1.0)) - 2.0*exp(-\$A*(r/\$r0-1.0))",
                           "Morse(A=\$A,r0=\$r0)")
```

use kwarg `cutoff` to set a cut-off, default is `Inf`
"""
AnalyticPotential

evaluate(p::AnalyticPotential, r) = p.f(r)
evaluate_d(p::AnalyticPotential, r) = p.f_d(r)
evaluate_dd(p::AnalyticPotential, r) = p.f_dd(r)
Base.print(io::Base.IO, p::AnalyticPotential) = print(io, p.id)
cutoff(p::AnalyticPotential) = p.cutoff

# construct from string or expression
AnalyticPotential(s::AbstractString; id = s, cutoff=Inf) =
                        AnalyticPotential(parse(s), id=id, cutoff=cutoff)

function AnalyticPotential(ex::Expr; id = string(ex), cutoff=Inf)
   @assert isa(id, AbstractString)
   f, f_d, f_dd = diff2(ex)
   return AnalyticPotential(F64fun(f), F64fun(f_d), F64fun(f_dd), id, cutoff)
end

# this is a hack to make tight-binding work; but it should be reconsidered
# right now I am thinking it is actually ok as is!
evaluate(p::AnalyticPotential, r, R) = evaluate(p, r)
evaluate_d(p::AnalyticPotential, r, R) = evaluate_d(p, r)


text:
@inline function Atomix.get(ref, order)
    ptr = Atomix.pointer(ref)
    root = Atomix.gcroot(ref)
    GC.@preserve root begin
        UnsafeAtomics.load(ptr, order)
    end
end

@inline function Atomix.set!(ref, v, order)
    v = Atomix.asstorable(ref, v)
    ptr = Atomix.pointer(ref)
    root = Atomix.gcroot(ref)
    GC.@preserve root begin
        UnsafeAtomics.store!(ptr, v, order)
    end
end

@inline function Atomix.replace!(ref, expected, desired, success_ordering, failure_ordering)
    expected = Atomix.asstorable(ref, expected)
    desired = Atomix.asstorable(ref, desired)
    ptr = Atomix.pointer(ref)
    root = Atomix.gcroot(ref)
    GC.@preserve root begin
        UnsafeAtomics.cas!(ptr, expected, desired, success_ordering, failure_ordering)
    end
end

@inline function Atomix.modify!(ref, op::OP, x, ord) where {OP}
    x = Atomix.asstorable(ref, x)
    ptr = Atomix.pointer(ref)
    root = Atomix.gcroot(ref)
    GC.@preserve root begin
        UnsafeAtomics.modify!(ptr, op, x, ord)
    end
end

Atomix.asstorable(ref, v) = convert(eltype(ref), v)


text:
macro ResultFields_Start()
    ex = quote
        nbins::Int
        dbulk::Float64
        cutoff::Float64
        d::Vector{Float64} = zeros(nbins)

        # Data to compute the MDDF distribution and corresponding KB integral
        md_count::Vector{Float64} = zeros(nbins)
        md_count_random::Vector{Float64} = zeros(nbins)
        sum_md_count::Vector{Float64} = zeros(nbins)
        sum_md_count_random::Vector{Float64} = zeros(nbins)
        mddf::Vector{Float64} = zeros(nbins)
        kb::Vector{Float64} = zeros(nbins)

        # Properties of the solute and solvent selections
        autocorrelation::Bool
        solvent::SolSummary
        solute::SolSummary
    end
    esc(ex)
end

macro ResultFields_AtomsMatrix()
    ex = quote
        # Atomic contributions to the MDDFs
        solute_atom::Array{Float64,2} = zeros(nbins, solute.natomspermol)
        solvent_atom::Array{Float64,2} = zeros(nbins, solvent.natomspermol)
    end
    esc(ex)
end

macro ResultFields_AtomsVector()
    ex = quote
        # Atomic contributions to the MDDFs
        solute_atom::Array{Float64} = zeros(nbins, solute.natomspermol)
        solvent_atom::Array{Float64} = zeros(nbins, solvent.natomspermol)
    end
    esc(ex)
end

macro ResultFields_End()
    ex = quote
        # Data to compute a RDF and the KB integral from this count
        rdf_count::Vector{Float64} = zeros(nbins)
        rdf_count_random::Vector{Float64} = zeros(nbins)
        sum_rdf_count::Vector{Float64} = zeros(nbins)
        sum_rdf_count_random::Vector{Float64} = zeros(nbins)
        rdf::Vector{Float64} = zeros(nbins)
        kb_rdf::Vector{Float64} = zeros(nbins)

        # Overall densities and volumes
        density::Density = Density()
        volume::Volume = Volume(nbins)

        # Options of the calculation
        options::Options
        irefatom::Int
        lastframe_read::Int
        nframes_read::Int

        # File name(s) of the trajectories in this results
        files::Vector{String}
        weights::Vector{Float64}
    end
    esc(ex)
end

"""

$(TYPEDEF)

Structure to contain the results of the MDDF calculation.

$(TYPEDFIELDS)

"""
@with_kw_noshow struct Result
    @ResultFields_Start()
    @ResultFields_AtomsMatrix()
    @ResultFields_End()
end

# The mutable version is used for reading saved data, because some vectors
# need to be reshaped
@with_kw_noshow mutable struct MutableResult
    @ResultFields_Start()
    @ResultFields_AtomsVector()
    @ResultFields_End()
end

#
# Initialize the data structure that is returned from the computation, and checks some
# input parameters for consistency
#

function Result(trajectory::Trajectory, options::Options; irefatom = -1)

    # Check for simple input errors
    if options.stride < 1
        error("in MDDF options: stride cannot be less than 1. ")
    end
    if options.lastframe > 0 && options.lastframe < options.firstframe
        error("in MDDF options: lastframe must be greater or equal to firstframe. ")
    end
    if options.lastframe > trajectory.nframes
        error("in MDDF options: lastframe is greater than trajectory.nframes. ")
    end

    # Check for problems in dbulk and cutoff definitions
    cutoff = options.dbulk
    if (options.dbulk / options.binstep) % 1 > 1.e-5
        error("in MDDF options: dbulk must be a multiple of binstep.")
    end
    if options.usecutoff
        if options.dbulk >= options.cutoff
            error(
                " in MDDF options: The bulk volume is zero (dbulk must be smaller than cutoff). ",
            )
        end
        if (options.cutoff / options.binstep) % 1 > 1.e-5
            error("in MDDF options: cutoff must be a multiple of binstep.")
        end
        if ((options.cutoff - options.dbulk) / options.binstep) % 1 > 1.e-5
            error("in MDDF options: (cutoff-dbulk) must be a multiple of binstep. ")
        end
        cutoff = options.cutoff
    end
    nbins = setbin(cutoff, options.binstep) - 1

    if options.irefatom > trajectory.solvent.natoms
        error(
            "in MDDF options: Reference atom index",
            options.irefatom,
            " is greater than number of " *
            "                 atoms of the solvent molecule. ",
        )
    end

    # Set reference atom as the closest one to the center of coordinates of the molecule, as default
    if irefatom == -1
        if options.irefatom == -1
            nextframe!(trajectory)
            xfirst = trajectory.x_solvent[1:trajectory.solvent.natomspermol]
            cm = centerofcoordinates(xfirst)
            dmin, one, irefatom = minimumdistance(cm, xfirst)
            firstframe(trajectory)
        else
            irefatom = options.irefatom
        end
    end

    # Last frame to be considered
    if options.lastframe == -1
        lastframe_read = trajectory.nframes
    else
        lastframe_read = options.lastframe
    end

    # Actual number of frames that are read considering lastframe and stride
    nframes_read = round(Int, (lastframe_read - options.firstframe + 1) / options.stride)
    if nframes_read == 0
        error("Number of frames to read is zero. Check input parameters.")
    end

    # Return data structure built up

    return Result(
        options = options,
        nbins = nbins,
        dbulk = options.dbulk,
        cutoff = cutoff,
        irefatom = irefatom,
        lastframe_read = lastframe_read,
        nframes_read = nframes_read,
        autocorrelation = isautocorrelation(trajectory),
        solute = SolSummary(trajectory.solute),
        solvent = SolSummary(trajectory.solvent),
        files = [trajectory.filename],
        weights = [1.0],
    )

end

#
# What to show at the REPL
#

Base.show(io::IO, R::Result) = Base.show(overview(R))


text:
"""
    jd2cal(dj1, dj2)

Julian Date to Gregorian year, month, day, and fraction of a day.

### Given ###

- `dj1`, `dj2`: Julian Date (Notes 1, 2)

### Returned ###

- `iy`: Year
- `im`: Month
- `id`: Day
- `fd`: Fraction of day

### Notes ###

1. The earliest valid date is -68569.5 (-4900 March 1).  The
   largest value accepted is 1e9.

2. The Julian Date is apportioned in any convenient way between
   the arguments dj1 and dj2.  For example, JD=2450123.7 could
   be expressed in any of these ways, among others:

   | `dj1`     |     `dj2`   | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

3. In early eras the conversion is from the "proleptic Gregorian
   calendar";  no account is taken of the date(s) of adoption of
   the Gregorian calendar, nor is the AD/BC numbering convention
   observed.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992),
    Section 12.92 (p604).

"""
function jd2cal(d1, d2)
    iy = Ref{Cint}()
    imo = Ref{Cint}()
    id = Ref{Cint}()
    fd = Ref{Cdouble}()
    i = ccall((:eraJd2cal, liberfa), Cint,
              (Cdouble, Cdouble, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cdouble}),
              d1, d2, iy, imo, id, fd)
    if i == -1
        throw(ERFAException("unacceptable date (Note 1)"))
    end
    iy[], imo[], id[], fd[]
end

"""
    jdcalf(ndp, dj1, dj2)

Julian Date to Gregorian Calendar, expressed in a form convenient
for formatting messages:  rounded to a specified precision.

### Given ###

- `ndp`: Number of decimal places of days in fraction
- `dj1`, `dj2`: Dj1+dj2 = Julian Date (Note 1)

### Returned ###

- `iymdf`: Year, month, day, fraction in Gregorian calendar

### Notes ###

1. The Julian Date is apportioned in any convenient way between
   the arguments dj1 and dj2.  For example, JD=2450123.7 could
   be expressed in any of these ways, among others:

   | `dj1`     |     `dj2`   | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

2. In early eras the conversion is from the "Proleptic Gregorian
   Calendar";  no account is taken of the date(s) of adoption of
   the Gregorian Calendar, nor is the AD/BC numbering convention
   observed.

3. Refer to the function [`jd2cal`](@ref).

4. NDP should be 4 or less if internal overflows are to be
   avoided on machines which use 16-bit integers.

### Called ###

- [`jd2cal`](@ref): JD to Gregorian calendar

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992),
    Section 12.92 (p604).

"""
function jdcalf(ndp, d1, d2)
    iymdf = zeros(Cint, 4)
    i = ccall((:eraJdcalf, liberfa), Cint,
              (Cint, Cdouble, Cdouble, Ptr{Cint}),
              ndp, d1, d2, iymdf)

    if i == -1
        throw(ERFAException("date out of range"))
    elseif i == 1
        @warn "NDP not 0-9 (interpreted as 0)"
    end
    iymdf[1], iymdf[2], iymdf[3], iymdf[4]
end


text:
using ManifoldsBase
using Test

struct ErrorTestManifold <: Manifold{ℝ} end

function ManifoldsBase.check_manifold_point(::ErrorTestManifold, x)
    if any(u -> u < 0, x)
        return DomainError(x, "<0")
    end
    return nothing
end
function ManifoldsBase.check_tangent_vector(M::ErrorTestManifold, x, v)
    mpe = check_manifold_point(M, x)
    mpe === nothing || return mpe
    if any(u -> u < 0, v)
        return DomainError(v, "<0")
    end
    return nothing
end

@testset "Domain errors" begin
    M = ErrorTestManifold()
    @test isa(check_manifold_point(M, [-1, 1]), DomainError)
    @test check_manifold_point(M, [1, 1]) === nothing
    @test !is_manifold_point(M, [-1, 1])
    @test is_manifold_point(M, [1, 1])
    @test_throws DomainError is_manifold_point(M, [-1, 1], true)

    @test isa(check_tangent_vector(M, [1, 1], [-1, 1]), DomainError)
    @test check_tangent_vector(M, [1, 1], [1, 1]) === nothing
    @test !is_tangent_vector(M, [1, 1], [-1, 1])
    @test is_tangent_vector(M, [1, 1], [1, 1])
    @test_throws DomainError is_tangent_vector(M, [1, 1], [-1, 1], true)
end


text:
using Distributed
@everywhere module Dalinar

using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Plots
using JLD: save, load
using Debugger
using Random
using DataStructures
using LinearAlgebra
using PyCall
using IterTools
using Distributions
using Distributed
using GaussianProcesses
using Optim

abstract type State end
abstract type Action end
abstract type Planner end
abstract type Parameters end
abstract type Transition end

include("graph/astar.jl")
# --------- gridworld ------------------
include("env/gridworld/gridworld.jl")
include("planner/gridworld/gridworld_planner.jl")
include("agent/gridworld/gridworld_agent.jl")
include("experiment/gridworld/run_all.jl")
# --------- mountain car ---------------
include("env/mountain_car/mountain_car.jl")
include("planner/mountain_car/mountain_car_planners.jl")
include("agent/mountain_car/mountain_car_agents.jl")
include("experiment/mountain_car/run_all.jl")
# --------- wide tree ------------------
include("env/wide_tree/wide_tree.jl")
include("agent/wide_tree/wide_tree_agnostic_sysid_agent.jl")
include("experiment/wide_tree/run_wide_tree_agnostic_sysid.jl")
end


text:
io = IOBuffer("""A,B,C
1,1,10
6,1""")

@test_throws CSV.ExpectedMoreColumnsError CSV.validate(io)

io = IOBuffer("""A;B;C
1,1,10
2,0,16""")
@test_throws CSV.TooManyColumnsError CSV.validate(io)

io = IOBuffer("""A;B;C
1,1,10
2,0,16""")
@test_throws CSV.ExpectedMoreColumnsError CSV.validate(io; delim=';')

io = IOBuffer("""a b c d e
1 2  3 4 5
1 2 3  4 5
1  2 3  4 5""")
@test_throws CSV.TooManyColumnsError CSV.validate(io; delim=' ')


text:
using Jultran
using PyPlot

νs = 25:.001:3000
κs = compute_κs(:ch4, νs, 260, 1.e4)
semilogy(νs,κs)


text:
module MPI

const tag_name = "mpi"

const augment = raw"""
    # Can't use Preferences since we might be running this very early with a non-existing Manifest
    MPIPreferences_UUID = Base.UUID("3da0fdf6-3ccc-4f1b-acd9-58baa6c99267")
    const preferences = Base.get_preferences(MPIPreferences_UUID)

    # Keep logic in sync with MPIPreferences.jl
    # FIXME: When MPIPreferences is registered both `binary` and `abi` should be const
    #        and the jll packages using this tag shall depend on MPIPreferences.jl
    function augment_mpi!(platform)
        binary = get(preferences, "binary", Sys.iswindows() ? "MicrosoftMPI_jll" : "MPICH_jll")

        abi = if binary == "system"
            get(preferences, "abi")
        elseif binary == "MicrosoftMPI_jll"
            "MicrosoftMPI"
        elseif binary == "MPICH_jll"
            "MPICH"
        elseif binary == "OpenMPI_jll"
            "OpenMPI"
        elseif binary == "MPItrampoline_jll"
            "MPIwrapper"
        else
            error("Unknown binary: $binary")
        end

        if !haskey(platform, "mpi")
            platform["mpi"] = abi
        end
        return platform
    end
"""

using BinaryBuilder, Pkg
using Base.BinaryPlatforms

mpi_abis = (
    ("MPICH", PackageSpec(name="MPICH_jll"), "", !Sys.iswindows) ,
    ("OpenMPI", PackageSpec(name="OpenMPI_jll"), "", !Sys.iswindows),
    ("MicrosoftMPI", PackageSpec(name="MicrosoftMPI_jll"), "", Sys.iswindows),
    ("MPIwrapper", PackageSpec(name="MPItrampoline_jll"), "", !Sys.iswindows)
)

function augment_platforms(platforms)
    all_platforms = AbstractPlatform[]
    dependencies = []
    for (abi, pkg, compat, f) in mpi_abis
        pkg_platforms = deepcopy(filter(f, platforms))
        foreach(pkg_platforms) do p
            p[tag_name] = abi
        end
        append!(all_platforms, pkg_platforms)
        push!(dependencies, Dependency(pkg; compat, platforms=pkg_platforms))
    end
    return all_platforms, dependencies
end

end


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

using ForwardDiff

abstract type Constraint{ T <: Real } end # do I really need the `T` for the supertype?

draw_instance(::Constraint{T}, ::Predictor{In,POut,T}, ::PredictionParam{T},
    ::Vector{Instance{In,Out,T}}) where {In,POut,Out,T} = #::ConstraintInstance{PIn,Out,T}
  error("unimplemented")

apply_operator(::Constraint{T}, ::Predictor{In,POut,T}, ::PredictionParam{T},
    ::ConstraintInstance{In,COut,T}) where {In,POut,COut,T} =
  error("unimplemented")

constraint_instance_type(::Constraint{T}, ::PredictionParam{T}, ::Predictor{In,POut,T}
    ) where {In,POut,T} =
  error("unimplemented")

struct ConservationConstraint{ T } <: Constraint{T}
  input_indices :: Matrix{Int} # batched as columns so that gradient only has to be taken once
  coefficients :: Matrix{T}
  positive :: Vector{Bool} # per column
  absolute :: Vector{Bool}
  mimicked_indices :: Vector{Int} # indices to grab from a random instance; recycled
end

constraint_target(::LomaxParam, (alpha, beta)::AbstractVector) =
  (alpha - 1) / beta # mode of the rate parameter

@generated function constraint_instance_type(::ConservationConstraint{T},
    ::Predictor{In,POut,T}, ::PredictionParam{T}) where {In,POut,T} # make this a macro? (harder to know types of arguments)
  ConstraintInstance{In,0,T} # do I need to quote this even though it's just a type value?
end

function draw_instance(constraint::ConservationConstraint{T}, predictor::Predictor{In,POut,T},
    parametrization::PredictionParam{T}, instances::Vector{Instance{In,Out,T}}
    )::ConstraintInstance{In,0,T} where {In,POut,Out,T}
  grabbed_instance = sample(instances)
  input = @MVector rand(T, In)
  for i in constraint.mimicked_indices # just overwrite previously randomized
    input[i] = grabbed_instance.input[i]
  end
  # I guess I could always subclass the ConstraintInstance?
  # generic structure seems to work fine for now.
  ConstraintInstance(
    convert(SVector, input),
    SVector{0,T}() )
end

function apply_operator(constraint::ConservationConstraint{T}, predictor::Predictor{In,POut,T},
    parametrization::PredictionParam{T}, instance::ConstraintInstance{In,0,T})::T where {In,POut,T}
  n_indices, n_laws = size(constraint.input_indices)
  @assert (n_indices, n_laws) == size(constraint.coefficients)
  n_targets = n_instance_outputs(predictor, parametrization)
  sum(1:n_targets) do target_index
    # recent versions of Zygote can ALMOST differentiate their own gradients
    derivatives = ForwardDiff.gradient(instance.input) do input # input -> begin ... end
      map_prediction(params -> constraint_target(parametrization, params),
        predictor, parametrization, input, target_index)
    end
    #derivatives, = grad(central_fdm(diff_order, 1), differentiable, instance.input)
    sum(1:n_laws) do law
      operator = sum(1:n_indices) do index
        input = constraint.input_indices[index, law]
        coefficient = constraint.coefficients[index, law]
        metric = derivatives[input]
        if constraint.absolute[law]
          metric = abs(metric)
        end
        coefficient * metric
      end
      if constraint.positive[law]
        min(operator, 0) ^ 2
      else
        operator ^ 2 # penalize nonzero-ness
      end
    end
  end
end

max_stars_repo_name: marmarelis/TreatmentCurves

max_issues_repo_name: marmarelis/TreatmentCurves

max_forks_repo_name: marmarelis/TreatmentCurves


text:
module CreateTableDiaries

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:diaries) do
    [
      primary_key()
      column(:date, :date)
      column(:weather, :string)
      column(:note, :text)
    ]
  end

  add_index(:diaries, :date)
end

function down()
  drop_table(:diaries)
end

end

max_stars_repo_name: ky-s/DiaryApp

max_issues_repo_name: ky-s/DiaryApp

max_forks_repo_name: ky-s/DiaryApp


text:
using Photometry
using BenchmarkTools
using CSV
using DataFrames
using ProgressLogging
using Random

rng = Random.seed!(11256)

data = randn(rng, 512, 512) .+ 10

rows = []
using Statistics
@progress for N in [1, 10, 50, 100, 200, 400, 500, 1000, 2000]
    xs, ys = size(data) .* (rand(rng, N), rand(rng, N))
    apers = CircularAperture.(xs, ys, 10)
    time = @belapsed photometry($apers, $data)
    push!(rows, (nt=Threads.nthreads(), N=N, time=time))
end

path = joinpath(@__DIR__, "julia_num_apertures.csv")
results = CSV.File(path)

rows_to_update = @. results[:nt] == Threads.nthreads()
results[rows_to_update] = DataFrame(rows)

CSV.write(path, results)

nothing


text:
import ..Datasets: dataset, datafile, reader, Mode, Test

qqp_init() = register(DataDep(
    "GLUE-QQP",
    """
    Quora Question Pairs (QQP) task (GLUE version)
    """,
    "https://dl.fbaipublicfiles.com/glue/data/QQP-clean.zip",
    "40e7c862c04eb26ee04b67fd900e76c45c6ba8e6d8fab4f8f1f8072a1a3fbae0";
    post_fetch_method = fn -> begin
      mv(fn, "QQP-clean.zip")
      DataDeps.unpack("QQP-clean.zip")
      innerdir = "QQP"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct QQP <: Dataset end

function dataset(::Type{M}, d::QQP) where M <: Mode
    ds = reader(datafile(M, d))
    header = split(take!(ds), '\t')
    field_num = length(header)
    needed_field = M == Test ? (2,3) : (4,5,6)
    rds = get_channels(String, length(needed_field))

    task = @async begin
        for s ∈ ds
            s = split(s, '\t')
            length(s) != field_num && continue
            for (i, j) ∈ enumerate(needed_field)
                put!(rds[i], s[j])
            end
        end
    end
    for rd ∈ rds
        bind(rd, task)
    end

    rds
end

function trainfile(::QQP)
    datadep"GLUE-QQP/train.tsv"
end

function devfile(::QQP)
    datadep"GLUE-QQP/dev.tsv"
end

function testfile(::QQP)
    datadep"GLUE-QQP/test.tsv"
end

get_labels(::QQP) = ("0", "1")


text:
using CUDA
using Adapt
using JSON
using StructArrays

struct Interpolate{CuArray}
    xs::CuArray
    ys::CuArray
end
#=
function (itp::Interpolate)(x)
    i = searchsortedfirst(itp.xs, x)
    i = clamp(i, firstindex(itp.ys), lastindex(itp.ys))
    @inbounds itp.ys[i]
end
=#
function Adapt.adapt_structure(to, itp::Interpolate)
    xs = Adapt.adapt_structure(to, itp.xs)
    ys = Adapt.adapt_structure(to, itp.ys)
    Interpolate(xs, ys)
end

function kernel(itp::Interpolate)
    @cuprintln("1")
    return
end

#@cuda kernel(itp)

struct TestStruct{A}
    z::A
    x::A
    y::A
end
Adapt.Adapt.@adapt_structure TestStruct

struct struct2{A}
    b::Bool
    a::A
end
Adapt.Adapt.@adapt_structure struct2

function kernel(test::struct2)
    @cushow(test.b)
    @cushow(test.a[1].x)
    return
end

struct pt
    x::Float64
    y::Float64
end

function initpt!(a)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(a)
        a[i] = pt(1.0, 2.0)
    end
    return nothing
end

ndata = 10240
thdx  = 256
numblocks = ceil(Int, ndata/thdx)

println("Starting TEST: Memory allocation and transfer")
#test = StructArray{TestStruct}((1,2))
#display(test)
cua1 = StructArray{TestStruct}((zeros(1),ones(1), ones(1)))
#display(cua1)
#test2 = replace_storage(CuArray,test)
#display(test2)
#@cuda kernel(test2)
s1 = struct2(false,cua1)

cua2 = replace_storage(CuArray, cua1)
s2 = struct2(false,cua2)
display(s2)
@cuda kernel(s2)


text:
module YAMLCalibrationFiles

using EasyFITS: FitsHeader

using ScientificDetectors
using YAML, FITSIO
using ScientificDetectors:CalibrationFrameSampler

"""
	fill_filedict!(filedict,catdict,dir)

fill dictionary `filedict` with fitsheader of all files in `dir` according to the keywords in `catdict`
The dictionary  `filedict` key is the filepath.

"""
function fill_filedict!(filedict::Dict{String, FitsHeader},
						catdict::Dict{String, Any},
						dir::String)
	for filename in readdir(dir; join=true, sort=false)
		if !contains(filename,catdict["exclude files"])
			if isfile(filename)
				if endswith(filename,catdict["suffixes"])
					get!(filedict, filename) do
						read(FitsHeader, filename)
					end
				end
			end
		elseif isdir(filename) && filedict["include sub directory"]
			fill_filedict!(catdict,filedict,filename)
		end
	end
end

function fill_filedict!(filedict::Dict{String, FitsHeader},
						catdict::Dict{String, Any},
						dirs::Vector{String})
	for dir in dirs
		fill_filedict!(filedict,catdict,dir)
	end
end

"""
	contains(chain::Union{AbstractString,Vector{AbstractString}}, pattern::Vector{AbstractString})

Overloading of Base.contains for vectors of string. Return `true` if a `chain`
contains one of the patterns given in `pattern`
"""
function Base.contains(chain::Union{String,Vector{String}}, pattern::Vector{String})
	for str in pattern
		if contains(chain,str)
			return true
		end
	end
	return false
end

"""
	contains(chain::Vector{String}, pattern::String)

Overloading of Base.contains for vectors of string. Return `true` if one of the `chains`
contains `pattern`
"""
function Base.contains(chains::Vector{String}, pattern::AbstractString)
	for chain in chains
		if contains(chain,pattern)
			return true
		end
	end
	return false
end

"""
	endswith(chain::Union{String,Vector{String}}, pattern::Vector{String})

Overloading of Base.endswith for vectors of string. Return `true` if `chain`
ends with one of the patterns given in `pattern`
"""
function Base.endswith(chain::Union{String,Vector{String}}, pattern::Vector{String})
	for str in pattern
		if endswith(chain,str)
			return true
		end
	end
	return false
end

"""
	endswith(chain::Vector{String}, pattern::String)

Overloading of Base.endswith for vectors of string. Return `true` if `chain`
ends with `pattern`
"""
function Base.endswith(chains::Vector{String},pattern::AbstractString)
	for chain in chains
		if endswith(chain,pattern)
			return true
		end
	end
	return false
end

"""
	filterfilename(filelist,name)

Remove all files from `filelist` that do not match the strings in `name`.
"""
function filterfilename(filelist::Dict{String, FitsHeader},name::String)
	return filterfilename(filelist,[name])
end

function filterfilename(filelist::Dict{String, FitsHeader},name::Vector{String})
	pattern =Regex("("*join(name,")|(")*")")
	return filter(p->match(pattern,p.first) !== nothing,filelist)
end

"""
	newlist = filtercat(filelist,keyword,value)

Build a `newlist` dictionnary of all files where `fitsheader[keyword] == value`.
"""
function filtercat(filelist::Dict{String, FitsHeader},
					keyword::String,
					values::Union{Vector{String}, Vector{Bool}, Vector{Integer}, Vector{AbstractFloat}})
	newlist = Dict{String, FitsHeader}()
	for value in values
		merge!(newlist, filtercat(filelist,keyword,value))
	end
	return newlist
end

function filtercat(filelist::Dict{String, FitsHeader},
					keyword::String,
					value::Union{String, Bool, Integer, AbstractFloat, Nothing})
	try tmp = filter(p->p.second[keyword] == value,filelist)
		return  tmp
	catch
		return Dict{String, FitsHeader}()
	end
end

"""
	newlist = filtercat(filelist::Dict{String, FitsHeader},catdict::Dict{String, Any})

Build a `newlist` dictionnary of all files where `fitsheader[keyword] == value` for all keywords contained in `catdict`
"""
function  filterkeyword(filelist::Dict{String, FitsHeader},
						catdict::Dict{String, Any})
	filteredkeywords = "(dir)|(files)|(suffixes)|(include subdirectory)|(exclude files)|(exptime)|(hdu)|(sources)"
	keydict =  filter(p->match(Regex(filteredkeywords), p.first) === nothing,catdict)
	if length(keydict)>0
		for (keyword,value) in keydict
			filelist =  filtercat(filelist,keyword,value)
		end
	end
	return filelist
end

function default_calibdict(dir::AbstractString)
	calibdict = Dict{String, Any}()
	calibdict["dir"] = dir
	calibdict["hdu"] = 1
	calibdict["suffixes"] = [".fits", ".fits.gz","fits.Z"]
	calibdict["include subdirectory"] = true
	calibdict["exclude files"] = Vector{String}()
	return calibdict
end

function default_category_dict(calibdict::Dict{String, Any})
	filteredkeywords = "(categories)|(title)";
 	catdict  = Dict{String, Any}()
	merge!(catdict,filter(p->match(Regex(filteredkeywords), p.first) === nothing,calibdict));
	return catdict
end

"""
	ReadCalibrationFiles(yaml_file::AbstractString; part::NTuple{2} = (:,:),  dir=pwd())

Process calibration files according to the YAML configuration file `yaml_file`.

- `part` keyword can be used to consider only a part of the detector (e.g. `part=(1:100,1:100)`) default `part=(:,:)`

- `dir` is the directory containing the files. By default `dir=pwd()`. This keyword is overriden by the `dir` in the YAML config file

Return an instance of `CalibrationData` with all information statistics needed to calibrate the detector.
"""
function ReadCalibrationFiles(yaml_file::AbstractString; part::NTuple{2} = (:,:),  dir = pwd())

	calibdict = default_calibdict(dir)
	#merge!(calibdict,vararg)
	merge!(calibdict,YAML.load_file(yaml_file; dicttype=Dict{String,Any}))

	filedict = Dict{String, FitsHeader}()

	catarr =  [CalibrationCategory(cata,Meta.parse(value["sources"])) for (cata,value) in calibdict["categories"] ]
	local caldat::CalibrationData{Float64}
	local roi::DetectorAxes
	local inds::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}
	isfirst = true
	width, height = -1, -1

	for (cat,value) in calibdict["categories"]
		catdict =default_category_dict(calibdict)
		merge!(catdict, value)
		fill_filedict!(filedict,calibdict,catdict["dir"])
		haskey(catdict,"files") && fill_filedict!(filedict,calibdict,catdict["files"])
		filescat = filterkeyword(filedict,catdict)
		if !isempty(filescat)
			for (filename,fitshead) in filescat
				hdu = FITS(filename)[catdict["hdu"]] :: ImageHDU

				if isfirst
					width, height = size(hdu)
					inds = (Base.OneTo(width)[part[1]],
					Base.OneTo(height)[part[2]])
					roi = DetectorAxes(inds)
					caldat = CalibrationData{Float64}(roi,catarr)
					isfirst = false
				else
					width == size(hdu,1)|| error("incompatible sizes")
					height == size(hdu,2) || error("incompatible sizes")
				end
				if ndims(hdu)==2
					data = read(hdu, inds...)
					sampler =  CalibrationDataFrame(cat,fitshead[catdict["exptime"]],data;roi=roi)
				else
				data = read(hdu, (inds...,Base.OneTo(size(hdu,3)))...)

				if size(hdu,3)>1
					sampler = CalibrationFrameSampler(data,cat,fitshead[catdict["exptime"]])
				else
					sampler =  CalibrationDataFrame(cat,fitshead[catdict["exptime"]],view(data, inds..., 1);roi=roi)
				end
			end
				push!(caldat, sampler)
			end
		end
	end
	return caldat
end

end


text:
using Test
using LRUSets

@testset "LRUSets.jl" begin

    s = LRUSet(3)

    @test s isa LRUSet{Any}
    @test s.max_size === 3

    @test length(s) === 0
    @test issetequal(values(s), [])

    push!(s, "foo")
    @test length(s) === 1
    @test issetequal(values(s), ["foo"])

    push!(s, "foo")
    @test length(s) === 1
    @test issetequal(values(s), ["foo"])

    push!(s, "bar")
    push!(s, "baz")
    push!(s, "faz")
    @test length(s) === 3
    @test issetequal(values(s), ["bar", "baz", "faz"])

    s = LRUSet(
        3,
        Dict(
            42 => UInt64(2),
            10 => UInt64(3),
        ),
    )

    @test s isa LRUSet{Int}
    @test eltype(LRUSet{Int}) === Int
    push!(s, 1)
    @test issetequal(values(s), [1, 42, 10])

    s = LRUSet{Integer}(2)
    @test s isa LRUSet{Integer}

    push!(s, Int16(21))
    push!(s, Int128(42))
    push!(s, UInt64(84))
    push!(s, UInt8(255))
    @test issetequal(values(s), [UInt64(84), UInt8(255)])

    @test issetequal([s...], [UInt64(84), UInt8(255)])

end


text:
function wordcount(sentence)

    tmp=lowercase(sentence)
    tmp=replace(tmp, "' "  =>"")
    tmp=replace(tmp," '" =>"")
    tmp=replace(tmp,"'" =>"T")
    words=split(tmp)
    for w in words
        tst=replace(w, !isletter =>"")
        println(tst)
    end
end
function wordcounts(sentence)
    prev = ""
    prevclass = 0
    res = Dict{String,Int32}()

    for c in lowercase(sentence)

        if isspace(c) || c==','
        # separate word
            if length(prev) > 0
                tmp = get!(res, prev, 0)
                push!(res, prev => tmp + 1)
                prev=""
            end

        # elseif isnumeric(c)
        #     print(c * " is numeric")
        #     if prevclass == 0
        #         # we shall add previous word

        #     elseif prevclass == 2
        #         # push!(res, "sdf" => [1, 2])
        #         res[prev] = res[prev] + 1
        #         prev = c
        #         prevclass = 1
        #     end

        elseif isletter(c) || isnumeric(c) || (c==''' && length(prev)>0)
            print(c * " is ascii")
            if prevclass == 0
                prev = prev * c
           elseif prevclass == 1
                throw("not implemented")
           end

        end
    end
    # adding remaining word to dictionary
    if length(prev) > 0
        tmp = get!(res, prev, 0)
        push!(res, prev => tmp + 1)
        # res[prev] = res[prev] + 1
    end

    return res
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-in

max_forks_repo_name: stepinski/julia-in


text:
using Statistics
import Statistics: centralize_sumabs2!

Statistics.var(xs::CuArray; corrected::Bool=true, dims=:, mean=mean(xs, dims=dims)) = Statistics.varm(xs, mean; corrected=corrected, dims=dims)
Statistics.std(xs::CuArray; corrected::Bool=true, dims=:, mean=mean(xs, dims=dims)) = sqrt.(var(xs, corrected=corrected, dims=dims, mean=mean))

function Statistics.centralize_sumabs2!(R::CuArray{S}, A::CuArray, means::CuArray) where S
    indsAt, indsRt = Base.safe_tail(axes(A)), Base.safe_tail(axes(R))
    keep, Idefault = Broadcast.shapeindexer(indsRt)
    len = length(CartesianIndices(indsAt))

    function kernel1!(r::CuDeviceArray{S}, a::CuDeviceArray, m::CuDeviceArray, i1)
        li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if li <= len && i <= size(a, 1)
            ia = CartesianIndices(indsAt)[li]
            ir = Tuple(Broadcast.newindex(ia, keep, Idefault))
            x = abs2(a[i, Tuple(ia)...] - m[i1, ir...])
            CUDA.atomic_add!(
                pointer(r,
                        Base._to_linear_index(r,
                                              i1,
                                              ir...
                                              )
                        ),
                x
            )
        end
        return
    end

    function kernel!(r::CuDeviceArray{S}, a::CuDeviceArray, m::CuDeviceArray)
        li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if li <= len && i <= size(a, 1)
            ia = CartesianIndices(indsAt)[li]
            ir = Tuple(Broadcast.newindex(ia, keep, Idefault))
            x = abs2(a[i, Tuple(ia)...] - m[i, ir...])
            CUDA.atomic_add!(
                pointer(r,
                        Base._to_linear_index(r,
                                              i,
                                              ir...
                                              )
                        ),
                x
            )
        end
        return
    end

    max_threads = 256
    thread_x = min(max_threads, size(A, 1))
    thread_y = min(max_threads ÷ thread_x, len)
    threads = (thread_x, thread_y)
    blocks = ceil.(Int, (size(A, 1), len) ./ threads)

    if Base.reducedim1(R, A)
        i1 = first(Base.axes1(R))
        CUDA.@cuda blocks=blocks threads=threads kernel1!(R, A, means, i1)
    else
        CUDA.@cuda blocks=blocks threads=threads kernel!(R, A, means)
    end

    return R
end


text:
#_______________________________________________________________________________
# Writting results

    # requirements fulfillement
RequirSheet = getSheet( CPworkbook, "Requirements" )
if RequirSheet.ptr === Ptr{Void}( 0 )
    RequirSheet = createSheet( CPworkbook, "Requirements" )
end

rowIdx = 0
Row = createRow(RequirSheet, rowIdx)
rowIdx += 1

cellIdx = 1
for i in 1:length(MPObjectives)
    obj = ""
    for j in 1:length(MPObjectives[i].Objectives)
        obj = obj * MPObjectives[i].Objectives[j]
    end
    Cell=createCell(Row, cellIdx); setCellValue(Cell, obj)
    cellIdx += 1
end
for i in 1:NBYears
    Row = createRow(RequirSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
    cellIdx += 1
    for j in 1:length(MPObjectives)
        Cell=createCell(Row, cellIdx); setCellValue(Cell, round(Int, RequirementFulfillment[i, j]))
        cellIdx += 1
    end
end

    # Initial manpower distribution
InitDistSheet = getSheet( CPworkbook, "Initial distribution" )
if InitDistSheet.ptr === Ptr{Void}( 0 )
    InitDistSheet = createSheet( CPworkbook, "Initial distribution" )
end

rowIdx = 0
Row = createRow(InitDistSheet, rowIdx)
rowIdx += 1

cellIdx = 2
#for i in 1:length(GuyCareerPaths)
#    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
#    cellIdx += 1
#end

index = 1
for i in 1:length(InitManpower)
    Row = createRow(InitDistSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, InitManpower[i].Academiclvl)
    cellIdx += 1
    Cell=createCell(Row, cellIdx); setCellValue(Cell, InitManpower[i].Seniority)
    cellIdx += 1
    for j in 1:length(GuyCareerPaths)
        if j in InitMPPartCP[i]
            if sol.sol[index] != 0
                Cell=createCell(Row, cellIdx); setCellValue(Cell, j); cellIdx += 1
                Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[index]); cellIdx += 1
            end
            index += 1
        end
#        Cell=createCell(Row, cellIdx);

#        if j in InitMPPartCP[i]
#            setCellValue(Cell, sol.sol[index])
#            index += 1
#        else
#            setCellValue(Cell, 0)
    end
end

    # Annual Recruitment
RecruitmentSheet = getSheet( CPworkbook, "Perforemed Recruitment" )
if RecruitmentSheet.ptr === Ptr{Void}( 0 )
    RecruitmentSheet = createSheet( CPworkbook, "Perforemed Recruitment" )
end

rowIdx = 0
#Row = createRow(RecruitmentSheet, rowIdx)
rowIdx += 1

#cellIdx = 1
#for i in 1:length(GuyCareerPaths)
#    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
#    cellIdx += 1
#end
#Cell=createCell(Row, cellIdx); setCellValue(Cell, "sum")

for i in 1:NBYears
    sum = 0
    Row = createRow(RecruitmentSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
    cellIdx += 2
    for j in 1:length(GuyCareerPaths)
        if sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j] != 0
            sum += sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j]
            Cell=createCell(Row, cellIdx); setCellValue(Cell, j)
            cellIdx += 1
            Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j])
            cellIdx += 1
        end
    end
    Cell=createCell(Row, 1); setCellValue(Cell, sum)
    cellIdx += 1
end

    # rectuirment Variables

# rectuirmentVariables = getSheet( CPworkbook, "rectuirmentVariables" )
# if rectuirmentVariables.ptr === Ptr{Void}( 0 )
#     rectuirmentVariables = createSheet( CPworkbook, "rectuirmentVariables" )
# end

    # Deviation Variables

DeviationVariables = getSheet( CPworkbook, "DeviationVariables" )
if DeviationVariables.ptr === Ptr{Void}( 0 )
    DeviationVariables = createSheet( CPworkbook, "DeviationVariables" )
end

rowIdx = 0
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(GuyCareerPaths) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + i])
    cellIdx += 1
    if i%length(GuyCareerPaths) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end
        #Positive deviation
rowIdx = 0
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(MPObjectives) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + AnnualRecDivNb + i])
    cellIdx += 1
    if i%length(MPObjectives) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end
        # Negative deviation
rowIdx += 2
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(MPObjectives) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + AnnualRecDivNb + (length(MPObjectives) * NBYears) + i])
    cellIdx += 1
    if i%length(MPObjectives) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end

#try rm(joinpath( dirname( Base.source_path() ), "..", "..", "..", "GuyCPResultsAll.xlsx" )) end
write(joinpath( dirname( Base.source_path() ), "..", "..", "..", "GuyCPResultsAllPop.xlsx" ), CPworkbook)

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using MatrixDepot

const spd_small = ("HB/bcsstk01", "FIDAP/ex5", "HB/494_bus", "HB/plat362")
const spd_med = ("HB/bcsstk09", "Bates/Chem97ZtZ", "HB/plat1919",
                      "Boeing/nasa1824", "Oberwolfach/t3dl_e")
const spd_large = ("GHS_psdef/torsion1", "Cannizzo/sts4098", "GHS_psdef/obstclae")
const matrix_path = dirname(pathof(MatrixDepot))

function get_matrices(matrices)
    for matrix in matrices
        try
            # throws an error if already downloaded
            matrixdepot(matrix, :get)
        catch
        end
    end
end

get_matrices(spd_small)
get_matrices(spd_med)
get_matrices(spd_large)


text:
@info "loading packages - DiffEq ecosystem"
using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
@info "loading packages - DynamicHMC ecosystem"
using DiffEqDynamicHMC, TransformVariables, DynamicHMC
@info "loading packages - rest"
using  Test, Distributions, LinearAlgebra, Random

@info "starting tests"
Random.seed!(1)

reporter = LogProgressReport(nothing, 100, 60.0) # frequent reporting prevents CI timeout

@testset "Lotka-Volterra 1 parameter" begin

    # define model
    f1 = @ode_def LotkaVolterraTest1 begin
        dx = a*x - x*y
        dy = -3*y + x*y
    end a

    # generate data
    p = [1.5]
    u0 = [1.0, 1.0]
    tspan = (0.0, 10.0)
    prob1 = ODEProblem(f1, u0, tspan, p)
    σ = 0.01                         # noise, fixed for now
    t = collect(range(1, 10; length = 10))   # observation times
    sol = solve(prob1, Tsit5())
    randomized = VectorOfArray([(sol(t[i]) + σ * randn(2)) for i in 1:length(t)])
    data = convert(Array, randomized)

    # inference
    mcmc_kwargs = (initialization = (q = zeros(1 + 2),), reporter = reporter)
    bayesian_result = dynamichmc_inference(prob1, t, data, (Normal(1.5, 1), ),
                                           as(Vector, asℝ₊, 1); mcmc_kwargs = mcmc_kwargs)

    # check
    @test mean(p.parameters[1] for p in bayesian_result.posterior) ≈ p[1] atol = 0.1
end

@testset "Lotka-Volterra 4 parameters" begin
    # define mdel
    f1 = @ode_def LotkaVolterraTest4 begin
        dx = a*x - b*x*y
        dy = -c*y + d*x*y
    end a b c d

    # generate data
    u0 = [1.0,1.0]
    tspan = (0.0,10.0)
    p = [1.5,1.0,3.0,1.0]
    prob1 = ODEProblem(f1,u0,tspan,p)
    sol = solve(prob1,Tsit5())
    t = collect(range(1,stop=10,length=10))
    randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
    data = convert(Array,randomized)

    # inference
    priors = (a = Truncated(Normal(1.5,0.01), 0, 2),
              b = Truncated(Normal(1.0,0.01), 0, 1.5),
              c = Truncated(Normal(3.0,0.01), 0, 4),
              d = Truncated(Normal(1.0, 0.01), 0, 2))

    mcmc_kwargs = (initialization = (q = zeros(4 + 2),), reporter = reporter)
    bayesian_result = dynamichmc_inference(prob1, t, data, priors, as(Vector, asℝ₊, 4);
                                           mcmc_kwargs = mcmc_kwargs)

    # check
    @test norm(mean([p.parameters for p in bayesian_result.posterior]) .- p, Inf) ≤ 0.1
end


text:
#=

# 225: Terminal flux calculation via test functions, nD
([source code](SOURCE_URL))

After calculating solutions based on the finite volume method, it
may be interesting to obtain information about the solution besides of the graphical representation.

Here, we focus on the following data:
- integrals of the solution
- flux through parts of the boundary

Let us define the following reaction - diffusion system in a domain $\Omega$:

```math
\begin{aligned}
\partial_t u_1 - \nabla \cdot \nabla u_1 + r(u_1, u_2) &= f=1.0\\
\partial_t u_2 - \nabla \cdot \nabla u_1 - r(u_1, u_2) &= 0
\end{aligned}
```

with boundary conditions $u_2=0$ on $\Gamma_2\subset\partial\Omega$ and
$r(u_1,u_2)=u_1 + 0.1 u_2$

The source $f$ creates species $u_1$ which reacts to $u_2$, $u_2$ then leaves
the domain at boundary $\Gamma_2$.

### Stationary problem

For the stationary problem, we have the following flux balances derived from the equations and from Gauss theorem:

```math
\begin{aligned}
\int_\Omega r(u_1,u_2) d\omega &= \int_\Omega f d\omega \\
\int_\Omega -r(u_1,u_2) d\omega &= \int_{\Gamma_2} \nabla u \cdot \vec n ds \\
\end{aligned}
```

The volume integrals can be approximated based on the finite volume subdivision
$\Omega=\cup_{i\in \mathcal N} \omega_i$:

```math
\begin{aligned}
\int_\Omega r(u_1,u_2) d\omega &\approx \sum_{i\in \mathcal N} |\omega_i| r(u_{1,i},u_{2,i})\\
\int_\Omega f d\omega &\approx \sum_{i\in \mathcal N} |\omega_i| f_i
\end{aligned}
```

But what about  the boundary integral ?  Here, we use a  trick to cast
the surface  integral to the  integral to  a volume integral  with the
help of a test function:

Let $T(x)$ be the solution of the Laplace problem $-\nabla^2 T =0$ in $\Omega$ and the boundary conditions

```math
\begin{aligned}
T &=0 \quad \text{at}\; \Gamma_4\\
T &=1 \quad \text{at}\; \Gamma_2\\
\partial_n T &=0\quad \text{at}\;  \Gamma_1,\Gamma_3
\end{aligned}
```

Write ``\vec j=-\nabla u``. and assume $\nabla\cdot \vec j + r =f$.

```math
\begin{aligned}
\int_{\Gamma_2} \vec j \cdot \vec n ds&=\int_{\Gamma_2} T\vec j \cdot \vec n ds \quad \text{due to}\; T=1\; \text{on}\; \Gamma_2\\
 &=\int_{\partial\Omega}  T\vec j \cdot \vec n ds\quad \text{due to}\; T=0\; \text{on}\; \Gamma_4, \quad\vec j\cdot \vec n=0\; \text{on}\; \Gamma_1, \Gamma_3\\
&= \int_\Omega \nabla \cdot (T \vec j) d\omega \quad \text{(Gauss)}\\
&= \int_\Omega \nabla T \cdot \vec j d\omega + \int_\Omega T \nabla\cdot j d\omega\\
&=  \int_\Omega \nabla T \cdot \vec j d\omega + \int_\Omega T(f-r)dω\\
\end{aligned}
```

and we approximate

```math
\begin{aligned}
\int_\Omega \nabla T \cdot \vec j d\omega \approx \sum_{k,l}
\frac{|\omega_k\cap\omega_l|}{h_{k,l}}g(u_k, u_l) (T_k-T_l)
\end{aligned}
```

where the sum runs over pairs of neigboring control volumes.

The `integrate` method with a  test function parameter returns a value
for each species, the sign convention assumes that species leaving the
domain lead to negative values.

### Transient problem

The amount  of species created via  the source term (measured  in `F`)
integrated  over time  should be  equal to  the sum  of the  amount of
species left in  the domain at the  very end of the  evolution and the
amount of species which left the domain:

$\int_{t_0}^{t_{end}} \int_\Omega f d\omega dt= \int_\Omega (u_1+u_2)dω + \int_{t_0}^{t_{end}} \int_{\Gamma_2} \nabla u_2 \cdot \vec n ds$

Literature references:

- H. Gajewski "Analysis und Numerik von Ladungstransport in Halbleitern", WIAS Berlin, Report No.6
- Yoder, P. D., K. Gärtner, and W. Fichtner. "A generalized Ramo–Shockley theorem for classical to quantum transport at arbitrary frequencies." Journal of Applied Physics 79.4 (1996): 1951-1954.
- P. Farrell, N. Rotundo, D. H. Doan, M. Kantner, J. Fuhrmann, and T. Koprucki, "Numerical methods for drift-diffusion models", in Handbook of optoelectronic device modeling and simulation: Lasers, modulators, photodetectors, solar cells, and numerical methods, vol. 2, J. Piprek, Ed. Boca Raton: CRC Press, 2017, pp. 733–771.

=#

module Example225_TestFunctions2D

using VoronoiFVM, GridVisualize,ExtendableGrids

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse,dim=2,tend=5, dt=0.2)
    n=[101,21,5]
    X=collect(range(0.0,1,length=n[dim]))
    if dim==1
        grid=simplexgrid(X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[1]
    elseif dim==2
        grid=simplexgrid(X,X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[4]
    elseif dim==3
        grid=simplexgrid(X,X,X)
        Γ_where_T_equal_1=[2]
        Γ_where_T_equal_0=[4]
    end

    function storage(f,u,node)
        f.=u
    end

    function flux(f,u,edge)
	f[1]=u[1,1]-u[1,2]
	f[2]=u[2,1]-u[2,2]
    end

    r(u1,u2)= u1-0.1*u2

    function reaction(f,u,node)
	f[1]= r(u[1],u[2])
	f[2]=-r(u[1],u[2])
    end

    function source(f,node)
	f[1]=1.0
    end

    physics=VoronoiFVM.Physics(flux=flux,
	                       storage=storage,
	                       reaction=reaction,
	                       source=source)

    system=VoronoiFVM.System(grid,physics)

    enable_species!(system,1,[1])
    enable_species!(system,2,[1])
    boundary_dirichlet!(system,2,2,0.0);

    inival=unknowns(system,inival=0.0)

    sol=solve(inival,system)

    vis=GridVisualizer(Plotter=Plotter,layout=(1,2),resolution=(600,300),fignumber=1)
    scalarplot!(vis[1,1],grid,sol[1,:],flimits=(0,1.5),title="u_1")
    scalarplot!(vis[1,2],grid,sol[2,:],flimits=(0,1.5),title="u_2",show=true)

    """
        The `integrate` method of `VoronoiFVM`  provides a possibility to calculate
        the volume integral of a function of a solution as described above.
        It returns a `num_specie` x `num_regions` matrix of the integrals
        of the function of the unknowns over the different subdomains (here, we have only one):
    """

    """
        Amount of u_1 and u_2 in the domain aka integral over identity storage function:
    """
    U=integrate(system,storage,sol)

    """
    Amount of species created by source term per unit time:
    """
    F=integrate(system,(f,u,node)->source(f,node),sol)

    """
    Amount of  reaction per unit time:
    """
    R=integrate(system,reaction,sol)

    tf=VoronoiFVM.TestFunctionFactory(system)
    T=testfunction(tf,Γ_where_T_equal_0,Γ_where_T_equal_1)

    I=integrate(system,T,sol)

    t0=0.0

    control=fixed_timesteps!(VoronoiFVM.NewtonControl(),dt)

    tsol=solve(inival,system,[t0,tend],control=control)

    vis1=GridVisualizer(Plotter=Plotter,layout=(1,2),resolution=(600,300),fignumber=4)

    for i=1:length(tsol)
        sol=tsol[i]
        scalarplot!(vis1[1,1],grid,sol[1,:],flimits=(0,1.5),clear=true)
        scalarplot!(vis1[1,2],grid,sol[2,:],flimits=(0,1.5),show=true)
    end

    outflow_rate=Float64[]
    for i=2:length(tsol)
	ofr=integrate(system,T,tsol[i],tsol[i-1],tsol.t[i]-tsol.t[i-1])
  	push!(outflow_rate,ofr[2])
    end

    vis2=GridVisualizer(Plotter=Plotter,layout=(1,1),resolution=(600,300),fignumber=2)
    scalarplot!(vis2[1,1],[0,tend],-[I[2],I[2]],label="stationary",clear=true)
    scalarplot!(vis2[1,1],tsol.t[2:end],-outflow_rate,label="transient",show=true)

    all_outflow=0.0
    for i=1:length(tsol)-1
	all_outflow-=outflow_rate[i]*(tsol.t[i+1]-tsol.t[i])
    end

    Uend=integrate(system,storage,tsol[end])
    isapprox(F[1], R[1],rtol=1.0e-12)  ? true : return false
    isapprox(I[1], 0.0, atol=1.0e-12)  ? true : return false
    isapprox(R[2], I[2],rtol=1.0e-12)  ? true : return false
    isapprox(F[1]*(tend-t0), ( Uend[1] + Uend[2] + all_outflow), rtol=1.0e-12) ? true : return false
end

function test()
    main(dim=1, unknown_storage=:sparse )  ? true : return false
    main(dim=1, unknown_storage=:dense  )  ? true : return false
    main(dim=2, unknown_storage=:sparse )  ? true : return false
    main(dim=2, unknown_storage=:dense  )  ? true : return false
    main(dim=3, unknown_storage=:sparse )  ? true : return false
    main(dim=3, unknown_storage=:dense  )  ? true : return false
end

end


text:
export QsosedModel
struct QsosedModel <: Model
    parameters::Parameters
    bh::BlackHole
    warm::Warm
    corona::Corona
end

function QsosedModel(parameters::Parameters)
    bh = BlackHole(parameters)
    corona = Corona(bh, parameters)
    warm = Warm(corona, parameters)
    return QsosedModel(parameters, bh, warm, corona)
end
function QsosedModel(config)
    parameters = Parameters(config)
    return QsosedModel(parameters)
end


text:
"""
An AbstractYield is an object which can be used as an argument to:

- zero-coupon spot rates via [`zero`](@ref)
- forward zero rates via [`forward`](@ref)
- discount factor via [`discount`](@ref)
- accumulation factor via [`accumulation`](@ref)

It can be be constructed via:

- zero rate curve with [`Zero`](@ref)
- forward rate curve with [`Forward`](@ref)
- par rate curve with [`Par`](@ref)
- typical OIS curve with [`OIS`](@ref)
- typical constant maturity treasury (CMT) curve with [`CMT`](@ref)
"""
abstract type AbstractYield end

abstract type AbstractYieldCurve <: AbstractYield end


text:
module PharmValuation
# import mean and disallow method extension
using Statistics: mean
using Printf: @printf
using Plots

struct CashProcess{T}
    annual_cash_flow::T
    drift::T
    volatility::T
    terminal_multiplier::T
    risk_premium::T
end

struct CostProcess{T}
    investment::T
    total_cost::T
    volatility::T
    failure_rate::T
end

struct SimOptions{T,U}
    time_step::T
    term::U
    runs::U
    basis::Function
end

struct ProjectProcess{T}
    CashProcess
    CostProcess
    correlation::T
    risk_free_rate::T
    SimOptions
end

function main(;runs=100_000)
    basis(x,y) = [1, x, y, x*y, x^2, y^2, x^2 * y, x * y^2, (x*y)^2]
    cash = CashProcess(20e6,0.02,35/100,5.0,0.036)
    cost = CostProcess(10e6,100e6,0.5,0.06931)
    sim = SimOptions(0.25,20,runs,basis)
    project = ProjectProcess(cash,cost,-0.1,0.05,sim)
    project_value = lsm(project)
    @printf "The project value: \$%.0f" project_value
end

""" simulate returns the cash flow and cost matrices generated by their respective processes """
function simulate(pp::ProjectProcess)
    # Set the number of periods
    periods = Int(ceil(pp.SimOptions.term / pp.SimOptions.time_step))

    # Risk adjusted cash flow drift rate
    adj_cash_flow = pp.CashProcess.drift - pp.CashProcess.risk_premium

    # Set return matrices and correlated random errors.
    net_cash = randn(pp.SimOptions.runs, periods)
    cost = pp.correlation * net_cash .+ sqrt(1 - pp.correlation^2) * randn(pp.SimOptions.runs,periods)

    # simulate the investment costs and correlated cash flows
    @inbounds for run in 1:pp.SimOptions.runs
        @inbounds for period in 1:periods
            # cash flow simulation
            prev_cash = period != 1 ? net_cash[run,period-1] : pp.CashProcess.annual_cash_flow
            net_cash[run,period] = prev_cash * exp((adj_cash_flow-0.5 * pp.CashProcess.volatility^2)*pp.SimOptions.time_step +
                pp.CashProcess.volatility * sqrt(pp.SimOptions.time_step) * net_cash[run,period])

            # cost simulation
            prev_cost = period != 1 ? cost[run,period-1] : pp.CostProcess.total_cost
            current_cost = 0.0
            if prev_cost != 0
                current_cost = prev_cost - pp.CostProcess.investment * pp.SimOptions.time_step +
                    pp.CostProcess.volatility * sqrt(pp.CostProcess.investment*prev_cost*pp.SimOptions.time_step) * cost[run,period]
                current_cost = current_cost < 0 ? 0 : current_cost
            end
            cost[run,period] = current_cost
        end
    end
    return net_cash, cost
end

""" lsm returns the least squares monte carlo estimate for the optimal stopping problem"""
function lsm(pp::ProjectProcess)
    periods = Int(ceil(pp.SimOptions.term / pp.SimOptions.time_step))

    # local for clarity
    period_investmt = pp.CostProcess.investment * pp.SimOptions.time_step

    # Simulate the net cash flows and cost values
    cash_matrix, cost_matrix = simulate(pp)

    # Set value function array to terminal value where investment has stopped.
    value_func = zeros(size(cash_matrix))
    terminal_cash_flow = cash_matrix[:,end]
    value_func[:,end] = pp.CashProcess.terminal_multiplier * (cost_matrix[:,end] .== 0) .* terminal_cash_flow

    # Discount rates for investing phase and sales phase
    cash_discount_rate = exp(-1 * pp.risk_free_rate * pp.SimOptions.time_step)
    invest_discount_rate = exp(-1 * (pp.risk_free_rate + pp.CostProcess.failure_rate)*pp.SimOptions.time_step)

    # value function iteration
    @inbounds for period in (periods-1):-1:1
        # initialize next periods value
        next_val = invest_discount_rate * value_func[:,period+1]

        # Apply basis function to cost and cash flows to approximate the discounted
        # next period value.
        basis_vecs = pp.SimOptions.basis.(cash_matrix[:,period],cost_matrix[:,period])

        # Transform vector of vector bases into a basis matrix for regression analysis.
        basis_matrix = reduce(hcat,basis_vecs)'

        # Solve regression coefficients
        coefficients = basis_matrix \ next_val

        # Estimate continuation value of investment
        est_val = basis_matrix * coefficients

        # Update value function array
        @inbounds for run in 1:pp.SimOptions.runs
            # Investing phase
            if cost_matrix[run,period] != 0
                # Investment decision criterion
                invest_val = est_val[run] - period_investmt

                # Update value function with for the investment phase.
                value_func[run,period] = invest_val > 0 ? next_val[run] - period_investmt : 0.0
            else
                # Update value function for the post-investment sales phase.
                value_func[run, period] = cash_matrix[run,period] * pp.SimOptions.time_step +
                    cash_discount_rate * value_func[run,period+1]
            end
        end
    end
    retval = invest_discount_rate * value_func[:,1]
    return mean(retval)
end

end # module

max_stars_repo_name: danzoppo/DataScienceExamples

max_issues_repo_name: danzoppo/DataScienceExamples

max_forks_repo_name: danzoppo/DataScienceExamples


text:
using Test
using MPI
using Distributed

# This uses TCP to communicate with the workers
mgr = MPI.start_main_loop(MPI.TCP_TRANSPORT_ALL)

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

refs = []
for w in workers()
    push!(refs, @spawnat w MPI.Comm_rank(MPI.COMM_WORLD))
end
ids = falses(size)
for r in refs
    id = fetch(r)
    @test !ids[id+1]
    ids[id+1] = true
end
@test ids[1] == (length(procs()) == 1)
ids[1] = true
for id in ids
    @test id
end

s = @distributed (+) for i in 1:10
    i^2
end
@test s == 385

MPI.stop_main_loop(mgr)


text:
import Base.@deprecate, Base.@deprecate_binding
import Base.depwarn

# v0.3
# ----

@deprecate ispalindrome ispalindromic

"""
Range type of biological symbols (especially `DNANucleotide`, `RNANucleotide`,
and `AminoAcid`).
"""
immutable SymbolRange{T} <: Range{T}
    start::UInt8
    stop::UInt8

    function SymbolRange(start::T, stop::T)
        warn("SymbolRange is deprecated")
        if start > stop
            # normalize empty range
            start = convert(T, 0x01)
            stop = convert(T, 0x00)
        end
        return new(start, stop)
    end
end

function SymbolRange{T}(start::T, stop::T)
    return SymbolRange{T}(start, stop)
end

function Base.show(io::IO, r::SymbolRange)
    show(io, first(r))
    write(io, ':')
    show(io, last(r))
end

# Iterator
# --------

Base.start(r::SymbolRange) = r.start
Base.done(r::SymbolRange, i) = i > r.stop
Base.next(r::SymbolRange, i) = reinterpret(eltype(r), i), i + 0x01

Base.size(r::SymbolRange) = (r.stop - r.start + 1,)
Base.first(r::SymbolRange) = convert(eltype(r), r.start)
Base.last(r::SymbolRange) = convert(eltype(r), r.stop)
Base.in{T}(x::T, r::SymbolRange{T}) = first(r) ≤ x ≤ last(r)

# Indexing
# --------

function Base.getindex(r::SymbolRange, i::Integer)
    checkbounds(r, i)
    return convert(eltype(r), r.start + i - 1)
end

function Base.getindex{T}(r::SymbolRange{T}, ir::UnitRange)
    checkbounds(r, ir)
    if isempty(ir)
        return SymbolRange(T(0x01), T(0x00))
    else
        return SymbolRange(r[first(ir)], r[last(ir)])
    end
end

immutable FASTA <: Bio.IO.FileFormat end
immutable FASTQ <: Bio.IO.FileFormat end
immutable TwoBit <: Bio.IO.FileFormat end
export FASTA, FASTQ, TwoBit