text:
using Base.Test
using RPOMDPs, RPOMDPModels
using RobustValueIteration
using SimpleProbabilitySets, RPOMDPToolbox, Distances
const RPBVI = RobustValueIteration

@testset "Robust and Info-Reward Baby POMDP" begin
    # constructor
    p = BabyPOMDP()
    ip = BabyIPOMDP()
    rp = BabyRPOMDP()
    rip = BabyRIPOMDP()
    @test (rip.r_feed, rip.r_hungry, rip.discount) == (-5., -10., 0.9)

    # s,a,z indexes
    # s = (true (hungry), false (full))
    # a = (true (feed), false (no action))
    # z = (:quiet, :crying, :yelling)
    @test state_index(p, :hungry) == 1
    @test action_index(ip, :nothing) == 2
    @test observation_index(rp, :yelling) == 3

    # transition function
    # P(hungry | full, no action) = 0.1
    @test pdf(transition(p, :full, :nothing), :hungry) == 0.1
    @test pdf(transition(ip, :full, :nothing), :hungry) == 0.1
    @test transition(rp, :full, :nothing, :hungry)[1] ≈ 0.075 atol = 1e-6
    @test transition(rip, :full, :nothing, :hungry)[1] ≈ 0.075 atol = 1e-6

    # observation function
    # P_lowercdf(yelling | hungry, feed) = 0.3
    # P_uppercdf(yelling | hungry, feed) = 0.05
    # note: p-box (not upper and lower probabilities)
    @test pdf(observation(p, :feed, :hungry), :quiet) == 0.3
    @test pdf(observation(ip, :feed, :hungry), :yelling) == 0.1
    @test observation(rp, :feed, :hungry)[1][1] ≈ 0.275 atol = 1e-6
    @test observation(rip, :feed, :hungry)[2][1] ≈ 0.325 atol = 1e-6

    # dynamics
    @test dynamics(p)[1,1,1,1] == 0.0
    @test dynamics(rip)[1][2,2,2,2] ≈ 0.065625 atol = 1e-6

    # reward
    # R(hungry, feed) = -15
    @test reward(p, :hungry, :feed) == -15
    @test reward(rp, :hungry, :feed) == -15
    @test reward(ip, [0.2, 0.8], :feed) ≈ 0.6 atol = 1e-6
    @test reward(rip, [0.2, 0.8], :feed) ≈ 0.6 atol = 1e-6

    # generate samples
    p = BabyPOMDP()
    rp = BabyRPOMDP()
    b = [0.8, 0.2]
    s = :hungry
    a = :nothing
    rng = MersenneTwister(20348)
    @test generate_sor(p, b, s, a, rng)[2] == :crying
    @test generate_sor(rp, b, s, a, rng)[2] == :crying
end

@testset "Cyber Assessment" begin
    p = CyberPOMDP()
    ip = CyberIPOMDP()
    rp = CyberRPOMDP()
    rip = CyberRIPOMDP()

    @test state_index(p, [3,3,3]) == 27
    @test action_index(ip, [2,3]) == 4
    @test observation_index(rp, [1,3]) == 3
    @test transition(ip, [1,1,1], [3,2]).probs |> sum ≈ 1.0
    @test transition(ip, [2,1,3], [3,2]).probs |> sum ≈ 1.0
    @test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,2]).probs)
    @test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,3]).probs)
    @test transition(rp, [2,2,2], [1,3], [3,2,3])[2] ≈ 0.053816
    @test RPOMDPModels.o([1,2], [2,2,2], [1,1], [0.8, 0.9]) ≈ 0.005
    @test observation(ip, [1,2], [1,1,1]).probs |> sum == 1.0
    @test observation(rp, [1,2], [2,1,3])[1] |> sum < 1.0
    @test dynamics(ip)[:,:,20,4] |> sum == 1
    @test all(dynamics(rip)[2] .>= dynamics(rip)[1])
    @test all(dynamics(rip)[2] .>= dynamics(ip))
    @test all(dynamics(rip)[1] .<= dynamics(ip))
    @test reward(rp, [1,3,2], [1,3]) ≈ 0.49
    @test reward(rip, fill(1/27,27), [3,3]) ≈ 0.0 atol = 1e-9
    @test reward(rip, vcat(1.0, zeros(26)), [1,3]) == 1.0

    b = vcat([0.8, 0.2], zeros(25))
    s = [1,2,3]
    a = [2,3]
    rng = MersenneTwister(20348)
    @test generate_sor(p, b, s, a, rng)[2] == [2,3]
    @test generate_sor(rp, b, s, a, rng)[2] == [2,2]

    nb = 5
    srand(47329342)
    bs = [psample(zeros(27), ones(27)) for i = 1:nb]

    maxiter = 5
    solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
    Vold = fill(RPBVI.AlphaVec(zeros(n_states(rip)), ordered_actions(rip)[1]),
        length(solver.beliefpoints))
    Volda = [Vold[i].alpha for i = 1:length(Vold)]
    e11 = zeros(27); e11[1] = 1.0;
    b = e11
    a = [1,2]
    ai = action_index(ip, a)
    umin, pmin = RPBVI.minutil(rip, b, a, Volda)
    @test pmin[:,:,1] |> sum ≈ 1.0 atol = 1e-9
    @test pmin[:,:,2] |> sum ≈ 1.0 atol = 1e-9
    @test pmin[:,:,13] |> sum ≈ 1.0 atol = 1e-9
    @test pmin[:,:,27] |> sum ≈ 1.0 atol = 1e-9
end

@testset "Cyber Assessment (Off Nominal)" begin
    ip = CyberTestIPOMDP()
    @test state_index(ip, [3,3,3]) == 27
    @test action_index(ip, [2,3]) == 4
    @test observation_index(ip, [1,3]) == 3
    @test transition(ip, [1,1,1], [3,2]).probs |> sum ≈ 1.0
    @test transition(ip, [2,1,3], [3,2]).probs |> sum ≈ 1.0
    @test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,2]).probs)
    @test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,3]).probs)
    @test transition(ip, [2,2,2], [1,3], [3,2,3]) ≈ 0.0495
    @test observation(ip, [1,2], [1,1,1]).probs |> sum == 1.0
    @test dynamics(ip)[:,:,20,4] |> sum == 1
    @test reward(ip, fill(1/27,27), [3,3]) ≈ 0.0 atol = 1e-9
    @test reward(ip, vcat(1.0, zeros(26)), [1,3]) == 1.0

    b = vcat([0.8, 0.2], zeros(25))
    s = [1,2,3]
    a = [2,3]
    rng = MersenneTwister(20348)
    @test generate_sor(ip, b, s, a, rng)[2] == [2,3]

    nb = 5
    srand(47329342)
    bs = [psample(zeros(27), ones(27)) for i = 1:nb]
    maxiter = 5
    solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
    polinit = RPBVI.create_policy(solver, ip)
    polip = RPBVI.solve(solver, ip, verbose = false)
    @test value(polip, b) > value(polinit, b)
end

@testset "Rock Diagnosis" begin
    ip = RockIPOMDP()
    rip = RockRIPOMDP()
    @test states(ip) == [10, 11, 20, 21]
    @test actions(ip) == [:left, :right, :check]
    @test observations(ip) == [:good, :bad, :none]
    @test states(rip) == [10, 11, 20, 21]
    @test actions(rip) == [:left, :right, :check]
    @test observations(rip) == [:good, :bad, :none]
    @test state_index.(ip, states(ip)) == [1,2,3,4]
    @test action_index.(ip, actions(ip)) == [1,2,3]
    @test observation_index.(ip, observations(ip)) == [1,2,3]

    s = 11
    a = :right
    z = :none
    sp = 21
    @test pdf(transition(ip, s, a), sp) == 1.0
    @test transition(ip, s, a, sp) == 1.0
    @test pdf(observation(ip, a, sp), z) == 1.0
    @test observation(ip, a, sp, z) == 1.0
    @test transition(rip, s, a)[1][4] == 0.999999
    @test transition(rip, s, a, sp) == (0.999999, 1.0)
    @test all(dynamics(rip)[2] .> dynamics(rip)[1])
    @test all(dynamics(rip)[2] .>= dynamics(ip))
    @test all(dynamics(ip) .>= dynamics(rip)[1])
    @test reward(ip, [0.5, 0.5, 0.0, 0.0], :left) ≈ 0.3333333 atol = 1e-5
    @test reward(rip, [0.5, 0.5, 0.0, 0.0], :left) ≈ 0.3333333 atol = 1e-5
    @test reward(ip, [1.0, 0.0, 0.0, 0.0], a) == 1.0
    @test reward(ip, [0.0, 1.0, 0.0, 0.0], a) == 1.0
    @test reward(ip, [0.0, 0.0, 1.0, 0.0], a) == 1.0
    @test reward(ip, [0.0, 0.0, 0.0, 1.0], a) == 1.0

    b = [0.7, 0.3, 0.0, 0.0]
    s = 11
    a = :check
    rng = MersenneTwister(0)
    @test generate_sor(ip, b, s, a, rng) == (11, :bad, 0.6)
    @test generate_sor(rip, b, s, a, rng) == (11, :good, 0.6)

    nb = 10
    srand(47329342)
    bs1 = [vcat(psample(zeros(2), ones(2)),zeros(2)) for i = 1:nb]
    bs2 = [vcat(zeros(2), psample(zeros(2), ones(2))) for i = 1:nb]
    bs = vcat(bs1, bs2)
    brand = rand() > 0.5 ? vcat(zeros(2), psample(zeros(2), ones(2))) : vcat(psample(zeros(2), ones(2)), zeros(2))
    @test minimum(norm(brand - bs[i]) for i=1:(2nb)) < 0.3

    maxiter = 100
    uncsize = 0.2
    solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
    Vold = fill(RPBVI.AlphaVec(zeros(n_states(rip)), ordered_actions(rip)[1]),
        length(solver.beliefpoints))
    Volda = [Vold[i].alpha for i = 1:length(Vold)]
    b = [0.620804, 0.379196, 0.0, 0.0]
    a = :left
    ai = action_index(ip, a)
    umin, pmin = RPBVI.minutil(rip, b, a, Volda)
    @test pmin[:,:,1] |> sum == 1.0
    @test pmin[:,:,2] |> sum == 1.0
    @test pmin[:,:,3] |> sum == 1.0
    @test pmin[:,:,4] |> sum == 1.0

    polinit = RPBVI.create_policy(solver, rip)
    polip = RPBVI.solve(solver, ip)
    polrip = RPBVI.solve(solver, rip)
    @test value(polip, [0.5, 0.5, 0.0, 0.0]) ≈ 18.5815 atol = 1e-3
    @test value(polip, [1.0, 0.0, 0.0, 0.0]) ≈ 19.8816 atol = 1e-3
    @test value(polinit, [0.5, 0.5, 0.0, 0.0]) ≈ 0.0 atol = 1e-3
    @test value(polinit, [1.0, 0.0, 0.0, 0.0]) ≈ 0.0 atol = 1e-3
    @show value(polrip, [0.5, 0.5, 0.0, 0.0])
    @show value(polrip, [1.0, 0.0, 0.0, 0.0])
    @test value(polrip, [0.5, 0.5, 0.0, 0.0]) ≈ 18.5051 atol = 1e-3
    @test value(polrip, [1.0, 0.0, 0.0, 0.0]) ≈ 19.8657 atol = 1e-3
end

######################
#
# Developmental Testing
#
#####################

using Base.Test
using RPOMDPs, RPOMDPModels
using RobustValueIteration
using SimpleProbabilitySets, RPOMDPToolbox, Distances
const RPBVI = RobustValueIteration

ip = CyberTestIPOMDP()

@test state_index(ip, [3,3,3]) == 27
@test action_index(ip, [2,3]) == 4
@test observation_index(ip, [1,3]) == 3
@test transition(ip, [1,1,1], [3,2]).probs |> sum ≈ 1.0
@test transition(ip, [2,1,3], [3,2]).probs |> sum ≈ 1.0
@test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,2]).probs)
@test all(transition(ip, [1,1,2], [2,1]).probs .== transition(ip, [1,1,2], [1,3]).probs)
@test transition(ip, [2,2,2], [1,3], [3,2,3]) ≈ 0.0495
@test observation(ip, [1,2], [1,1,1]).probs |> sum == 1.0
@test dynamics(ip)[:,:,20,4] |> sum == 1
@test reward(ip, fill(1/27,27), [3,3]) ≈ 0.0 atol = 1e-9
@test reward(ip, vcat(1.0, zeros(26)), [1,3]) == 1.0

b = vcat([0.8, 0.2], zeros(25))
s = [1,2,3]
a = [2,3]
rng = MersenneTwister(20348)
@test generate_sor(ip, b, s, a, rng)[2] == [2,3]

nb = 5
srand(47329342)
bs = [psample(zeros(27), ones(27)) for i = 1:nb]
maxiter = 5
solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
polinit = RPBVI.create_policy(solver, ip)
polip = RPBVI.solve(solver, ip, verbose = true)
@test value(polip, b) > value(polinit, b)


text:
const TEST_UUID = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
const TEST_DEP = PackageSpec(; name="Test", uuid=TEST_UUID)

"""
    Tests(; file="$(contractuser(default_file("test", "runtests.jl")))", project=false)

Sets up testing for packages.

## Keyword Arguments
- `file::AbstractString`: Template file for `runtests.jl`.
- `project::Bool`: Whether or not to create a new project for tests (`test/Project.toml`).
  See [here](https://julialang.github.io/Pkg.jl/v1/creating-packages/#Test-specific-dependencies-in-Julia-1.2-and-above-1)
  for more details.

!!! note
    Managing test dependencies with `test/Project.toml` is only supported
    in Julia 1.2 and later.
"""
@plugin struct Tests <: FilePlugin
    file::String = default_file("test", "runtests.jl")
    project::Bool = false
end

source(p::Tests) = p.file
destination(::Tests) = joinpath("test", "runtests.jl")
view(::Tests, ::Template, pkg::AbstractString) = Dict("PKG" => pkg)

function validate(p::Tests, t::Template)
    invoke(validate, Tuple{FilePlugin, Template}, p, t)
    p.project && t.julia < v"1.2" && @warn string(
        "Tests: The project option is set to create a project (supported in Julia 1.2 and later) ",
        "but a Julia version older than 1.2 ($(t.julia)) is supported by the template",
    )
end

function hook(p::Tests, t::Template, pkg_dir::AbstractString)
    # Do the normal FilePlugin behaviour to create the test script.
    invoke(hook, Tuple{FilePlugin, Template, AbstractString}, p, t, pkg_dir)

    # Then set up the test depdendency in the chosen way.
    f = p.project ? make_test_project : add_test_dependency
    f(pkg_dir)
end

# Create a new test project.
function make_test_project(pkg_dir::AbstractString)
    with_project(() -> Pkg.add(TEST_DEP), joinpath(pkg_dir, "test"))
end

# Add Test as a test-only dependency.
function add_test_dependency(pkg_dir::AbstractString)
    # Add the dependency manually since there's no programmatic way to add to [extras].
    path = joinpath(pkg_dir, "Project.toml")
    toml = TOML.parsefile(path)
    get!(toml, "extras", Dict())["Test"] = TEST_UUID
    get!(toml, "targets", Dict())["test"] = ["Test"]
    open(io -> TOML.print(io, toml), path, "w")

    # Generate the manifest by updating the project.
    # This also ensures that keys in Project.toml are sorted properly.
    touch(joinpath(pkg_dir, "Manifest.toml"))  # File must exist to be modified by Pkg.
    with_project(Pkg.update, pkg_dir)
end

ext: jl


text:
push!(LOAD_PATH, "../src/")

using Documenter, AnimatedDistributions

makedocs(
    sitename = "AnimatedDistributions",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [AnimatedDistributions],
    pages = [
       "index.md",
       "animations.md",
       "animate_beta.md",
       "animate_normal.md",
       "animate_cauchy.md",
       "animate_tstudent.md",
       "animate_exponential.md",
       "animate_frechet.md",
       "animate_gamma.md",
       "animate_invgamma.md",
       "animate_weibull.md",
       "animate_semicircle.md",
       "animate_rayleigh.md",
       "animate_pareto.md",
       "animate_triweight.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/gongcastro/AnimatedDistributions.jl.git",
    devbranch = "main"
)


text:
############################
## To be refined....
##
############################
function mdtest(mf::MomentFunction,
                div::Divergence,
                θ₀::Vector,
                lb::Vector, ub::Vector;
                solver=IpoptSolver())

    model = MathProgSolverInterface.model(solver)

    m = mf.nmom
    n = mf.nobs
    k = mf.npar

    u₀ = [ones(n), θ₀]

    gele = int((n+k)*(m+1)-k)
    hele = int(n*k + n + (k+1)*k/2)

    # lb = [-2., -2]
    # ub = [2, 2.]
    g_L = [zeros(m), n];
    g_U = [zeros(m), n];

    u_L = [zeros(n),  lb];
    u_U = [ones(n)*n, ub];
    # l = [1,1,1,1]
    # u = [5,5,5,5]
    # lb = [25, 40]
    # ub = [Inf, 40]

    MdProb = MinimumDivergenceProblem(mf, div, n, m, k, gele, hele, Array(Float64, n), Array(Float64, m+1))

    MathProgSolverInterface.loadnonlinearproblem!(model, n+k, m+1, u_L, u_U, g_L, g_U, :Min, MdProb)

    MathProgSolverInterface.setwarmstart!(model, u₀)

    MathProgSolverInterface.optimize!(model)
    stat = MathProgSolverInterface.status(model)

    # @test stat == :Optimal
    uᵒ = MathProgSolverInterface.getsolution(model)
    Qᵒ = MathProgSolverInterface.getobjval(model)
    (model, uᵒ, Qᵒ)

end


text:
module MoonGlitchBackgroundTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Moon Glitch Background" => Ahorn.EntityPlacement(
        Maple.MoonGlitchBackgroundTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.MoonGlitchBackgroundTrigger)
    return Dict{String, Any}(
        "duration" => Maple.moon_glitch_background_trigger_durations
    )
end

end


text:
"
    use an idmrg2 step to truncate/expand the bond dimension
"
@with_kw struct VumpsSvdCut <: Algorithm
    tol_gauge = Defaults.tolgauge
    tol_galerkin = Defaults.tol
    tol_eigenval = Defaults.tol
    trscheme = notrunc()
end

function changebonds_1(state::InfiniteMPS, H::Hamiltonian,alg::VumpsSvdCut,envs=environments(state,H)) #would be more efficient if we also repeated envs
    #the unitcell==1 case is unique, because there you have a sef-consistency condition

    #expand the one site to two sites
    nstate = InfiniteMPS(repeat(state.AL,2))
    nH = repeat(H,2)

    nstate,nenvs = changebonds(nstate,nH,alg)

    A1 = nstate.AL[1]; A2 = nstate.AL[2]
    D1 = space(A1,1); D2 = space(A2,1)

    #collapse back to 1 site
    if D2 != D1
        (nstate,nenvs) = changebonds(nstate,nH,SvdCut(trscheme = truncspace(infimum(D1,D2))),nenvs)
    end

    collapsed = InfiniteMPS([nstate.AL[1]],nstate.CR[1], tol = alg.tol_gauge);

    return collapsed, envs
end

function changebonds_n(state::InfiniteMPS, H::Hamiltonian,alg::VumpsSvdCut,envs=environments(state,H))
    meps=0.0
    for loc in 1:length(state)
        @plansor AC2[-1 -2;-3 -4] := state.AC[loc][-1 -2;1]*state.AR[loc+1][1 -4;-3]

        (vals,vecs,_) = eigsolve(@closure(x->ac2_prime(x,loc,state,envs)),AC2, 1, :SR, tol = alg.tol_eigenval; ishermitian=false )
        nAC2 = vecs[1]

        (vals,vecs,_)  = eigsolve(@closure(x->c_prime(x,loc+1,state,envs)),state.CR[loc+1], 1, :SR, tol = alg.tol_eigenval; ishermitian=false )
        nC2 = vecs[1]

        #svd ac2, get new AL1 and S,V ---> AC
        (AL1,S,V,eps) = tsvd(nAC2, trunc=alg.trscheme, alg=TensorKit.SVD())
        @plansor AC[-1 -2;-3] := S[-1;1]*V[1;-3 -2]
        meps = max(eps,meps)

        #find AL2 from AC and C as in vumps paper
        QAC,_ = leftorth(AC, alg = QRpos())
        QC ,_ = leftorth(nC2 , alg = QRpos())
        dom_map = isometry(domain(QC),domain(QAC))

        @plansor AL2[-1 -2;-3] := QAC[-1 -2;1]*conj(dom_map[2;1])*conj(QC[-3;2])

        #make a new state using the updated A's
        copied = copy(state.AL)
        copied[loc]   = AL1;
        copied[loc+1] = AL2;
        state = InfiniteMPS(copied,tol = alg.tol_gauge)
    end

    return state, envs
end

function changebonds(state::InfiniteMPS,H,alg::VumpsSvdCut,envs=environments(state,H))
    if (length(state) == 1)
        return changebonds_1(state,H,alg,envs)
    else
        return changebonds_n(state,H,alg,envs);
    end
end


text:
export DistProduct, ProdGeneric

import Distributions
import Base: prod, show

"""
    DistProduct

If inference backend cannot return an analytical solution for a product of two distributions it may fallback to the `DistProduct` structure
`DistProduct` is useful to propagate the exact forms of two messages until it hits some approximation method for form-constraint.
However `DistProduct` cannot be used to compute statistics such as mean or variance.
It has to be approximated before using in actual inference procedure.

Backend exploits form constraints specification which usually help to deal with intractable distributions products.
User may use EM form constraint with a specific optimisation algorithm or it may approximate intractable product with Gaussian Distribution
using for example Laplace approximation

See also: [`prod`](@ref)
"""
struct DistProduct{ L, R }
    left  :: L
    right :: R
end

Base.show(io::IO, product::DistProduct) = print(io, "DistProduct(", getleft(product), ",", getright(product), ")")

getleft(product::DistProduct)  = product.left
getright(product::DistProduct) = product.right

function Distributions.support(product::DistProduct)
    lsupport = Distributions.support(getleft(product))
    rsupport = Distributions.support(getright(product))
    if lsupport != rsupport
        error("Product $product has different support for left and right entries.")
    end
    return lsupport
end

Distributions.mean(product::DistProduct)      = error("mean() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.median(product::DistProduct)    = error("median() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.mode(product::DistProduct)      = error("mode() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.shape(product::DistProduct)     = error("shape() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.scale(product::DistProduct)     = error("scale() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.rate(product::DistProduct)      = error("rate() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.var(product::DistProduct)       = error("var() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.std(product::DistProduct)       = error("std() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.cov(product::DistProduct)       = error("cov() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.invcov(product::DistProduct)    = error("invcov() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.logdetcov(product::DistProduct) = error("logdetcov() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.entropy(product::DistProduct)   = error("entropy() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Distributions.params(product::DistProduct)    = error("params() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")

Distributions.pdf(product::DistProduct, x)    = Distributions.pdf(product.left, x) * Distributions.pdf(product.right, x)
Distributions.logpdf(product::DistProduct, x) = Distributions.logpdf(product.left, x) + Distributions.logpdf(product.right, x)

Base.precision(product::DistProduct) = error("precision() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Base.length(product::DistProduct)    = error("length() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Base.ndims(product::DistProduct)     = error("ndims() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
Base.size(product::DistProduct)      = error("size() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")

probvec(product::DistProduct)         = error("probvec() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")
weightedmean(product::DistProduct)    = error("weightedmean() is not defined for $(product). DistProduct structure has to be approximated and cannot be used in inference procedure.")

variate_form(::P)                         where { P <: DistProduct } = variate_form(P)
variate_form(::Type{ DistProduct{L, R} }) where { L, R }             = _check_dist_product_variate_form(variate_form(L), variate_form(R))

_check_dist_product_variate_form(::Type{ F }, ::Type{ F })   where { F <: VariateForm }                     = F
_check_dist_product_variate_form(::Type{ F1 }, ::Type{ F2 }) where { F1 <: VariateForm, F2 <: VariateForm } = error("DistProduct has different variate forms for left ($F1) and right ($F2) entries.")

value_support(::P)                         where { P <: DistProduct } = value_support(P)
value_support(::Type{ DistProduct{L, R} }) where { L, R }             = _check_dist_product_value_support(value_support(L), value_support(R))

_check_dist_product_value_support(::Type{ S }, ::Type{ S })   where { S <: ValueSupport }                      = S
_check_dist_product_value_support(::Type{ S1 }, ::Type{ S2 }) where { S1 <: ValueSupport, S2 <: ValueSupport } = error("DistProduct has different value supports for left ($S1) and right ($S2) entries.")

"""
    ProdGeneric{C}

`ProdGeneric` is one of the strategies for `prod` function. This strategy mimics `prod_constraint` constraint but does not fail in case of no analytical rule is available.

See also: [`prod`](@ref), [`ProdAnalytical`](@ref), [`ProdPreserveType`](@ref)
"""
struct ProdGeneric{C}
    prod_constraint :: C
end

get_constraint(prod_generic::ProdGeneric) = prod_generic.prod_constraint

ProdGeneric() = ProdGeneric(ProdAnalytical())

prod(::ProdGeneric, ::Missing, right)     = right
prod(::ProdGeneric, left, ::Missing)      = left
prod(::ProdGeneric, ::Missing, ::Missing) = missing

prod(generic::ProdGeneric, left::L, right::R) where { L, R } = prod(generic, prod_analytical_rule(L, R), left, right)

prod(generic::ProdGeneric, ::ProdAnalyticalRuleAvailable, left, right) = prod(get_constraint(generic), left, right)
prod(generic::ProdGeneric, ::ProdAnalyticalRuleUnknown, left, right)   = DistProduct(left, right)

# In case of ProdPointMass we want to propagate a single `DistProduct` as much as possible and do not create a big tree of product which will reduce performance significantly
# In this methods the general rule is the folowing: If we see that one of the arguments of `DistProduct` has the same function form
# as second argument of `prod` function it is better to try to `prod` them together with `NoConstraint` strategy.
prod(generic::ProdGeneric, left::DistProduct{L, R}, right::T) where { L, R, T } = prod(generic, prod_analytical_rule(L, T), prod_analytical_rule(R, T), left, right)
prod(generic::ProdGeneric, left::T, right::DistProduct{L, R}) where { L, R, T } = prod(generic, prod_analytical_rule(T, L), prod_analytical_rule(T, R), left, right)

prod(generic::ProdGeneric, ::ProdAnalyticalRuleUnknown,   ::ProdAnalyticalRuleUnknown,   left::DistProduct, right) = DistProduct(left, right)
prod(generic::ProdGeneric, ::ProdAnalyticalRuleAvailable, ::ProdAnalyticalRuleUnknown,   left::DistProduct, right) = DistProduct(prod(get_constraint(generic), getleft(left), right), getright(left))
prod(generic::ProdGeneric, ::ProdAnalyticalRuleUnknown,   ::ProdAnalyticalRuleAvailable, left::DistProduct, right) = DistProduct(getleft(left), prod(get_constraint(generic), getright(left), right))

prod(generic::ProdGeneric, ::ProdAnalyticalRuleUnknown,   ::ProdAnalyticalRuleUnknown,   left, right::DistProduct) = DistProduct(left, right)
prod(generic::ProdGeneric, ::ProdAnalyticalRuleAvailable, ::ProdAnalyticalRuleUnknown,   left, right::DistProduct) = DistProduct(prod(get_constraint(generic), left, getleft(right)), getright(right))
prod(generic::ProdGeneric, ::ProdAnalyticalRuleUnknown,   ::ProdAnalyticalRuleAvailable, left, right::DistProduct) = DistProduct(getleft(right), prod(get_constraint(generic), left, getright(right)))

function prod(generic::ProdGeneric, left::DistProduct{L1, R1}, right::DistProduct{L2, R2}) where {L1, R1, L2, R2}
    return prod(
        generic,
        prod_analytical_rule(L1, L2), prod_analytical_rule(L1, R2),
        prod_analytical_rule(R1, L2), prod_analytical_rule(R1, R2),
        left, right
    )
end

prod(::ProdGeneric, _, _, _, _, left::DistProduct, right::DistProduct) = DistProduct(left, right)

function prod(generic::ProdGeneric, ::ProdAnalyticalRuleAvailable, _, _, ::ProdAnalyticalRuleAvailable, left::DistProduct, right::DistProduct)
    return prod(generic, prod(get_constraint(generic), getleft(left), getleft(right)), prod(get_constraint(generic), getright(left), getright(right)))
end


text:
#!/usr/bin/env julia
# coding: utf-8

#rabin karp is a string search algorithm
#it is useful for string pattern matching
#particularly in plagiarism detection
#this algorithm leverages hash function instead of direct substring matching
#more details can be found in the link below
# http://www-igm.univ-mlv.fr/~lecroq/string/node5.html#SECTION0050
#in theory. a wise choice of hash function
#can make its execution hypersonic compared to naïve search
#more details of hash function can be found in the link below
# https://github.com/je-suis-tm/search-and-sort/blob/master/hash%20search.jl

pattern="Canary Wharf"

rawtext="""The main landlord in London's Canary Wharf financial district is considering scrapping plans for a new office skyscraper and building a 60-story apartment tower instead.

Canary Wharf Group is consulting with local authorities about the potential for building around 700 homes at its 1 Park Place site that would be purpose-built for tenants, according to an application on the Tower Hamlets council website. The firm had won initial permission in 2015 for a one million square foot office building, the equivalent of two Gherkin skyscrapers.

The reason for the potential switch to apartments wasn't disclosed in the document. A spokesperson for Canary Wharf declined to comment.

The pandemic has hammered London's commercial landlords, who have grappled with plunging rent collections and the ongoing uncertainty about whether white collar workers will return to the office en masse. Canary Wharf is facing a fall by as much as 15% in advertised office rents through 2021, the steepest throughout London, according to a forecast by broker Carter Jonas.

Read more: Canary Wharf Plans Bond as Offices Grapple With Empty Space

Focusing on rental homes, however, looks an increasingly safe bet as a housing shortage across the capital fuels demand for accommodation. DWS, the investment arm of Deutsche Bank AG, last month made its first acquisition in the U.K. build-to-rent sector with a development in South London, citing a lack of supply of good quality and affordable rental housing in the area.

Even before the pandemic, financial firms were reconsidering their real estate needs. HSBC Holdings Plc, which has its headquarters in Canary Wharf, said it expected to slash its global property footprint by 40%. So far, initial steps to ease lockdown restrictions haven't led to a rush back to the office: occupancy levels across the U.K. are around 30% this week, according to building software company Metrikus.

Still, a flurry of building approvals and recent leasings in the City of London, Canary's long-time rival district, suggest that there is still high demand for new offices in the capital. British Land Co. on Thursday said it had leased almost 30% of its planned 1 Broadgate redevelopment to Jones Lang LaSalle Inc. -- four years before the project's planned completion.

The Docklands area around Canary Wharf has staged a residential building boom in the past five years, with the construction of luxury apartment towers targeting young finance workers. Around 9% of new homes under construction across London at the end of last year were in the area's E14 postcode, according to data company Molior.

The news about Park Place was earlier reported by Estates Gazette."""

#simple hash function according to wikipedia
# https://en.wikipedia.org/wiki/Rabin%E2%80%93Karp_algorithm
#my preference is to choose a small prime number
#a smaller prime number will not exceed float64 limit in julia
function hash_function(word,prime_num=3)

    hash_value=sum([Float64(word[i])*(prime_num^i) for i in 1:length(word)])

    return hash_value

end

#naïve search iterates the letter one by one
function naive_search(pattern,rawtext)

    len_pattern=length(pattern)
    len_rawtext=length(rawtext)
    output=[]

    #this part is stupid
    #as julia allows us to do
    #rawtext[i:i+len_pattern-1]==pattern
    for i in 1:(len_rawtext-len_pattern+1)
        ignore=false
        j=1
        while !ignore
            if rawtext[i+j-1]!=pattern[j]
                ignore=true
            end
            j+=1
            if j>len_pattern
                if !ignore
                    push!(output,i)
                end
                ignore=true

            end
        end
    end
    return output
end

#instead of checking rawtext[i:i+len_pattern-1]==pattern
#rabin karp algorithm leverages hash function
function rabin_karp(pattern,rawtext,prime_num=3)

    output=[]
    len_pattern=length(pattern)
    len_rawtext=length(rawtext)
    hash_pattern=hash_function(pattern,prime_num)
    hash_rawtext=hash_function(rawtext[1:len_pattern],
    prime_num)

    #hash
    for i in 1:(len_rawtext-len_pattern+1)

        #rolling hash
        if i>1
            hash_rawtext-=hash_function(rawtext[i-1:i-1],prime_num)
            hash_rawtext/=prime_num
            hash_rawtext+=hash_function(rawtext[i+len_pattern-1:i+len_pattern-1],
                                        prime_num)*(prime_num^(len_pattern-1))
        end

        if hash_rawtext==hash_pattern
            push!(output,i)
        end
    end

    return output

end

print(rabin_karp(pattern,rawtext)==naive_search(pattern,rawtext))

@time rabin_karp(pattern,rawtext)

#ironically, naïve search seems to be a lot faster than rabin karp...
@time naive_search(pattern,rawtext)

max_stars_repo_path: rabin karp search.jl

max_stars_repo_name: tattooday/search-and-sort

max_issues_repo_path: rabin karp search.jl

max_issues_repo_name: tattooday/search-and-sort

max_forks_repo_path: rabin karp search.jl

max_forks_repo_name: tattooday/search-and-sort


text:
using Distributed
using LinearAlgebra #, Statistics
using Printf, WriteVTK, AbaqusReader
using PyCall, PyPlot, JLD, ProgressMeter

using AD4SM

@show Elements.setp(4)
;

mean(x) = sum(x)/length(x)

# this function is the constraint equation for the internal holes
function cnst_func(u,ucntr,r0,cntr0,Rsq)
    r    = r0+u
    cntr = cntr0+ucntr
    sum((r-cntr).^2)-Rsq
end
# this function find the centre of the constrained holes in the reference configuration
function find_centre(nodes; dTol = 1e-12, maxiter=21)

  c0    = mean(nodes)
  bdone = false
  iter  = 0

  while !bdone
    D2c0   = adiff.D2(c0)
    global deltas = [(D2c0[1]-p[1]).^2+(D2c0[2]-p[2]).^2 for p in nodes]
    r      = norm(deltas)
    grad   = adiff.grad(r)
    if (norm(grad)<dTol) || (iter>maxiter)
      bdone = true
    else
      c0  -= adiff.hess(r)\grad
    end
    iter +=1
  end
  Rsquares = adiff.val.(deltas)
  return (c0, Rsquares)
end

# some plotting functions
patch = pyimport("matplotlib.patches")
coll  = pyimport("matplotlib.collections")

function plot_model(elems, nodes;
                    u = zeros(length(nodes[1]), length(nodes)),
                    Φ = [],
                    linewidth = 0.25,
                    facecolor = :c,
                    edgecolor = :b,
                    alpha     = 1,
                    cmap      = :hsv,
                    clim      = [],
                    dTol      = 1e-6,
                    cfig      = figure(),
                    ax        = cfig.add_subplot(1,1,1))

  nodes     = [node + u[:,ii] for (ii,node) in enumerate(nodes)]
  patchcoll = coll.PatchCollection([patch.Polygon(nodes[elem.nodes])
                                    for elem ∈ elems], cmap=cmap)
  if !isempty(Φ)
    patchcoll.set_array(Φ)
    cfig.colorbar.(patchcoll, ax=ax)

    if isempty(clim)
      clim = patchcoll.get_clim()
      if abs(clim[2]-clim[1]) < dTol
        clim  = sum(clim)/2*[0.9, 1.1]
      end
    end
    patchcoll.set_clim(clim)
  else
    patchcoll.set_color(facecolor)
  end

  patchcoll.set_edgecolor(edgecolor)
  patchcoll.set_alpha(alpha)
  patchcoll.set_linewidth(linewidth)

  ax.set_aspect("equal")
  ax.add_collection(patchcoll)
  ax.autoscale()
  return (cfig, ax, patchcoll)
end
function get_I1(elems, u0)

    nElems = length(elems)
    F      = Elements.getinfo(elems,u0,info=:F)
    J      = [det(F)                       for F  in F]
    C      = [transpose(F)*F               for F  in F]
    L3     = [(C[1]C[4]-C[2]C[3])^-1       for C  in C]
    Ic     = [Materials.getInvariants(C[ii], L3[ii]) for ii in 1:nElems]

    I1     = [item[1] for item in Ic]
end
;

sMeshFile = "Pattern2D03FinerMesh02j.inp"
ϵn0       = 0.24
bfside    = true
mat       = Materials.NeoHooke(10)
bprogress = false
bwarmup   = true
blheqs    = false
bsheqs    = true
sFileName = splitext(sMeshFile)[1]
;

mymodel   = AbaqusReader.abaqus_read_mesh(sMeshFile)
@printf("\n\t loaded %s \n", sMeshFile)

(xmin, xmax) = extrema([item[2][1] for item ∈ mymodel["nodes"]])
(ymin, ymax) = extrema([item[2][2] for item ∈ mymodel["nodes"]])

nodes    = [mymodel["nodes"][ii]    for ii ∈ 1:length(mymodel["nodes"])]
el_nodes = [mymodel["elements"][ii] for ii ∈ 1:length(mymodel["elements"])]

elems = []
for item ∈ el_nodes
if length(item)==3
  push!(elems, Elements.Tria(item, nodes[item], mat=mat))
elseif length(item)==4
  push!(elems, Elements.Quad(item, nodes[item], mat=mat))
end
end

points   = hcat(nodes...)
cells    = [if length(nodes)==3
            MeshCell(VTKCellTypes.VTK_TRIANGLE, nodes)
          else length(nodes)==4
            MeshCell(VTKCellTypes.VTK_QUAD, nodes)
          end  for nodes in el_nodes ]

@show nNodes, nElems  = length(nodes), length(el_nodes)
@show Δx, Δy          = (xmax-xmin), (ymax-ymin)

nid_bndl = mymodel["node_sets"]["ID_L"]
nid_bndr = mymodel["node_sets"]["ID_R"]
nid_bndt = mymodel["node_sets"]["ID_T"]
nid_bndb = mymodel["node_sets"]["ID_B"]

nid_sh01 = mymodel["node_sets"]["ID_SH01"]
nid_sh02 = mymodel["node_sets"]["ID_SH02"]
nid_sh03 = mymodel["node_sets"]["ID_SH03"]
nid_sh04 = mymodel["node_sets"]["ID_SH04"]
nid_sh05 = mymodel["node_sets"]["ID_SH05"]
nid_sh06 = mymodel["node_sets"]["ID_SH06"]
nid_sh07 = mymodel["node_sets"]["ID_SH07"]
nid_sh08 = mymodel["node_sets"]["ID_SH08"]
nid_sh09 = mymodel["node_sets"]["ID_SH09"]

nid_shx  = [nid_sh01, nid_sh02, nid_sh03,
          nid_sh04, nid_sh05, nid_sh06,
          nid_sh07, nid_sh08, nid_sh09]
nSHs     = length(nid_shx)

nid_lh01 = mymodel["node_sets"]["ID_LH01"]
nid_lh02 = mymodel["node_sets"]["ID_LH02"]
nid_lh03 = mymodel["node_sets"]["ID_LH03"]
nid_lh04 = mymodel["node_sets"]["ID_LH04"]
nid_lh05 = mymodel["node_sets"]["ID_LH05"]
nid_lh06 = mymodel["node_sets"]["ID_LH06"]
nid_lh07 = mymodel["node_sets"]["ID_LH07"]
nid_lh08 = mymodel["node_sets"]["ID_LH08"]
nid_lh09 = mymodel["node_sets"]["ID_LH09"]
nid_lh10 = mymodel["node_sets"]["ID_LH10"]
nid_lh11 = mymodel["node_sets"]["ID_LH11"]
nid_lh12 = mymodel["node_sets"]["ID_LH12"]
nid_lh13 = mymodel["node_sets"]["ID_LH13"]
nid_lh14 = mymodel["node_sets"]["ID_LH14"]
nid_lh15 = mymodel["node_sets"]["ID_LH15"]
nid_lh16 = mymodel["node_sets"]["ID_LH16"]

nid_lhx  = [nid_lh01, nid_lh02, nid_lh03,
          nid_lh04, nid_lh05, nid_lh06,
          nid_lh07, nid_lh08, nid_lh09, nid_lh10,
          nid_lh11, nid_lh12, nid_lh13,
          nid_lh14, nid_lh15, nid_lh16]
nLHs     = length(nid_lhx)

cntrs_s = [find_centre(nodes[idx]) for idx in nid_shx]
cntrs_l = [find_centre(nodes[idx]) for idx in nid_lhx]
nodes = vcat(nodes,getindex.(cntrs_s,1),getindex.(cntrs_l,1))
;

plot_model(elems, nodes, facecolor=:c, edgecolor=:c)
PyPlot.title("undeformed model")
;

idxs  = LinearIndices(((2, nNodes+nSHs+nLHs)))
eqns  = Array{Elements.ConstEq}(undef, 0)
if bsheqs
    eqns  = vcat(eqns,
                 [ [ begin
                      r0      = nodes[nid]
                      cntr0   = nodes[nNodes+ii]
                      Rsq     = cntrs_s[ii][2][jj]
                      id_dofs = vcat(idxs[:,nid][:],idxs[:,nNodes+ii][:])
                      Elements.ConstEq(x->cnst_func(x[1:2],x[3:4],r0,cntr0,Rsq),
                                       id_dofs,adiff.D2)
                    end  for (jj,nid) in enumerate(nid_ii)]
                  for (ii,nid_ii) in enumerate(nid_shx)]...)
end
if blheqs
    eqns = vcat(eqns,
                [ [ begin
                     r0      = nodes[nid]
                     cntr0   = nodes[nNodes+nSHs+ii]
                     Rsq     = cntrs_l[ii][2][jj]
                     id_dofs = vcat(idxs[:,nid][:],idxs[:,nNodes+nSHs+ii][:])
                     Elements.ConstEq(x->cnst_func(x[1:2],x[3:4],r0,cntr0,Rsq),
                                      id_dofs,adiff.D2)
                   end  for (jj,nid) in enumerate(nid_ii)]
                 for (ii,nid_ii) in enumerate(nid_lhx)]... )
end
;

u                  = fill(NaN, 2, nNodes+nSHs+nLHs)
u[:,nid_bndb]     .= 0
u[1,nid_bndt]     .= 0
u[2,nid_bndt]     .= -ϵn0*Δy

if !bsheqs u[:,nNodes+1:nNodes+nSHs] .= NaN; end
if !blheqs u[:,nNodes+nSHs+1:end]    .= NaN; end

if !bfside
    u[1,nid_bndl] .= 0
    u[1,nid_bndr] .= 0
end

ifree     = isnan.(u)
icnst     = .!ifree
;

u0         = 1e-4Δx*randn(2, nNodes+nSHs+nLHs)
u0[icnst] .= 0
unew       = copy(u0)
@time Elements.solvestep!(elems, u0, unew, ifree, eqns=eqns, λ = zeros(length(eqns)),
                          bprogress=false, becho=true)
u[ifree] .= unew[ifree]
;

N       = 24
LF_c    = vcat(range(0.0, 0.9, length=2N÷8),
               range(0.9, 1.0, length=6N÷8))

allus_c = Elements.solve(elems, u, LF=LF_c, ifree=ifree, eqns=eqns,
                          bprogress=false, becho=true, bechoi=true)
;

rf_tot_c = [sum(item[2][2,nid_bndt])   for item in allus_c]
Δu_tot_c = [mean(item[1][2,nid_bndt])  for item in allus_c]
;

I1 = get_I1(elems, allus_c[end][1])

cfig = figure()
ax1   = cfig.add_subplot(2,1,1)

plot_model(elems, nodes, alpha=0.05,
  facecolor=:c, edgecolor=:c, cfig=cfig, ax=ax1)
plot_model(elems, nodes, u = allus_c[end][1],
  edgecolor=:c, Φ = get_I1(elems, allus_c[end][1]), cfig=cfig, ax=ax1)
title("deformed model - compression")

ax2   = cfig.add_subplot(2,1,2)
ax2.plot(Δu_tot_c, rf_tot_c)
xlabel("Δu_tot_")
ylabel("rf_tot_c")
title("force-displacement, compressive branch")

cfig.tight_layout(pad=2.0, w_pad=2.0, h_pad=2.0)
;

u                  = fill(NaN, 2, nNodes+nSHs+nLHs)
u[:,nid_bndb]     .= 0
u[1,nid_bndt]     .= 0
u[2,nid_bndt]     .= ϵn0*Δy

if !bsheqs u[:,nNodes+1:nNodes+nSHs] .= NaN; end
if !blheqs u[:,nNodes+nSHs+1:end]    .= NaN; end

if !bfside
    u[1,nid_bndl]     .= 0
    u[1,nid_bndr]     .= 0
end

ifree     = isnan.(u)
icnst     = .!ifree
;

u0         = 1e-4Δx*randn(2, nNodes+nSHs+nLHs)
u0[icnst] .= 0
unew       = copy(u0)
@time Elements.solvestep!(elems, u0, unew, ifree, eqns=eqns, λ = zeros(length(eqns)),
                          bprogress=false, becho=true)
u[ifree] .= unew[ifree]
;

N       = 10
LF_t    = vcat(range(0.0, 0.9, length=2N÷3),
               range(0.9, 1.0, length=N÷3+1))
println("\n\t starting the tensile branch \n"); flush(stdout)
allus_t = Elements.solve(elems, u, LF=LF_t, ifree=ifree, eqns=eqns,
                          bprogress=false, becho=true, bechoi=true)
;

rf_tot_t = [sum(item[2][2,nid_bndt])   for item in allus_t]
Δu_tot_t = [mean(item[1][2,nid_bndt])  for item in allus_t]
;

JLD.save(sFileName*".jld",
            "nodes", nodes, "elems", elems,
            "points", points, "cells", cells, "mat", mat,
            "nid_bndl", nid_bndl, "nid_bndr", nid_bndr,
            "nid_bndt", nid_bndt, "nid_bndb", nid_bndb,
            "allus_t", allus_t, "allus_c", allus_c,
            "LF_t", LF_t, "LF_c", LF_c,
            "rf_tot_t", rf_tot_t, "Δu_tot_t", Δu_tot_t,
            "rf_tot_c", rf_tot_c, "Δu_tot_c", Δu_tot_c,
            "cntrs_s", cntrs_s, "cntrs_l", cntrs_l,
            "nid_shx", nid_shx, "nid_lhx", nid_lhx,
            "Δx", Δx, "Δy", Δy, "nNodes", nNodes, "nElems", nElems,
            "sMeshFile", sMeshFile, "sFileName", sFileName)
@printf("results written to %s\n", sFileName); flush(stdout)
;

I1    = get_I1(elems, allus_t[end][1])

cfig  = figure()
ax1   = cfig.add_subplot(2,1,1)

plot_model(elems, nodes, alpha=0.05,
  facecolor=:c, edgecolor=:c, cfig=cfig, ax=ax1)
plot_model(elems, nodes, u = allus_t[end][1],
  edgecolor=:c, Φ = get_I1(elems, allus_t[end][1]), cfig=cfig, ax=ax1)
title("deformed model - tensile")

ax2   = cfig.add_subplot(2,1,2)
ax2.plot(Δu_tot_t, rf_tot_t)
xlabel("Δu_tot")
ylabel("rf_tot")
title("force-displacement, tensile branch")

cfig.tight_layout(pad=2.0, w_pad=2.0, h_pad=2.0)
;

PyPlot.figure()
PyPlot.plot(Δu_tot_c, rf_tot_c, Δu_tot_t, rf_tot_t)
PyPlot.xlabel("Δu_tot_")
PyPlot.ylabel("rf_tot_c")
PyPlot.title("force-displacement, total")
;


text:
# Autogenerated wrapper script for Enzyme_jll for i686-w64-mingw32-cxx03-julia_version+1.7.0
export libEnzyme

JLLWrappers.@generate_wrapper_header("Enzyme")
JLLWrappers.@declare_library_product(libEnzyme, "libEnzyme-12.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libEnzyme,
        "bin\\libEnzyme-12.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MyPackage1
using Documenter

DocMeta.setdocmeta!(MyPackage1, :DocTestSetup, :(using MyPackage1); recursive=true)

makedocs(;
    modules=[MyPackage1],
    authors="Sandesh Katakam",
    repo="https://github.com/sandeshkatakam/MyPackage1.jl/blob/{commit}{path}#{line}",
    sitename="MyPackage1.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sandeshkatakam.github.io/MyPackage1.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sandeshkatakam/MyPackage1.jl",
    devbranch="master",
)


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Financial

Financial data struct with inner constructor:
```julia
function Financial(;
    om_cost_escalation_pct::Real = 0.025,
    elec_cost_escalation_pct::Real = 0.019,
    boiler_fuel_cost_escalation_pct::Real = 0.034,
    chp_fuel_cost_escalation_pct::Real = 0.034,
    generator_fuel_cost_escalation_pct::Real = 0.027,
    offtaker_tax_pct::Real = 0.26,
    offtaker_discount_pct::Real = 0.0564,
    third_party_ownership::Bool = false,
    owner_tax_pct::Real = 0.26,
    owner_discount_pct::Real = 0.0564,
    analysis_years::Int = 25,
    value_of_lost_load_per_kwh::Union{Array{R,1}, R} where R<:Real = 1.00,
    microgrid_upgrade_cost_pct::Real = off_grid_flag ? 0.0 : 0.3,
    macrs_five_year::Array{Float64,1} = [0.2, 0.32, 0.192, 0.1152, 0.1152, 0.0576],  # IRS pub 946
    macrs_seven_year::Array{Float64,1} = [0.1429, 0.2449, 0.1749, 0.1249, 0.0893, 0.0892, 0.0893, 0.0446],
    offgrid_other_capital_costs::Real = 0.0, # only applicable when off_grid_flag is true. Straight-line depreciation is applied to this capex cost, reducing taxable income.
    offgrid_other_annual_costs::Real = 0.0 # only applicable when off_grid_flag is true. Considered tax deductible for owner. Costs are per year.
)
```

!!! note
    When `third_party_ownership` is `false` the offtaker's discount and tax percentages are used throughout the model:
    ```julia
        if !third_party_ownership
            owner_tax_pct = offtaker_tax_pct
            owner_discount_pct = offtaker_discount_pct
        end
    ```
"""
struct Financial
    om_cost_escalation_pct::Real
    elec_cost_escalation_pct::Real
    boiler_fuel_cost_escalation_pct::Real
    chp_fuel_cost_escalation_pct::Real
    generator_fuel_cost_escalation_pct::Real
    offtaker_tax_pct::Real
    offtaker_discount_pct::Real
    third_party_ownership::Bool
    owner_tax_pct::Real
    owner_discount_pct::Real
    analysis_years::Int
    value_of_lost_load_per_kwh::Union{Array{R,1}, R} where R<:Real
    microgrid_upgrade_cost_pct::Real
    macrs_five_year::Array{Float64,1}
    macrs_seven_year::Array{Float64,1}
    offgrid_other_capital_costs::Real
    offgrid_other_annual_costs::Real

    function Financial(;
        off_grid_flag::Bool = false,
        om_cost_escalation_pct::Real = 0.025,
        elec_cost_escalation_pct::Real = 0.019,
        boiler_fuel_cost_escalation_pct::Real = 0.034,
        chp_fuel_cost_escalation_pct::Real = 0.034,
        generator_fuel_cost_escalation_pct::Real = 0.027,
        offtaker_tax_pct::Real = 0.26,
        offtaker_discount_pct::Real = 0.0564,
        third_party_ownership::Bool = false,
        owner_tax_pct::Real = 0.26,
        owner_discount_pct::Real = 0.0564,
        analysis_years::Int = 25,
        value_of_lost_load_per_kwh::Union{Array{R,1}, R} where R<:Real = 1.00,
        microgrid_upgrade_cost_pct::Real = off_grid_flag ? 0.0 : 0.3,
        macrs_five_year::Array{Float64,1} = [0.2, 0.32, 0.192, 0.1152, 0.1152, 0.0576],  # IRS pub 946
        macrs_seven_year::Array{Float64,1} = [0.1429, 0.2449, 0.1749, 0.1249, 0.0893, 0.0892, 0.0893, 0.0446],
        offgrid_other_capital_costs::Real = 0.0, # only applicable when off_grid_flag is true. Straight-line depreciation is applied to this capex cost, reducing taxable income.
        offgrid_other_annual_costs::Real = 0.0 # only applicable when off_grid_flag is true. Considered tax deductible for owner.
    )

        if off_grid_flag && !(microgrid_upgrade_cost_pct == 0.0)
            @warn "microgrid_upgrade_cost_pct is not applied when off_grid_flag is true. Setting microgrid_upgrade_cost_pct to 0.0."
            microgrid_upgrade_cost_pct = 0.0
        end

        if !off_grid_flag && (offgrid_other_capital_costs != 0.0 || offgrid_other_annual_costs != 0.0)
            @warn "offgrid_other_capital_costs and offgrid_other_annual_costs are only applied when off_grid_flag is true. Setting these inputs to 0.0 for this grid-connected analysis."
            offgrid_other_capital_costs = 0.0
            offgrid_other_annual_costs = 0.0
        end

        if !third_party_ownership
            owner_tax_pct = offtaker_tax_pct
            owner_discount_pct = offtaker_discount_pct
        end

        return new(
            om_cost_escalation_pct,
            elec_cost_escalation_pct,
            boiler_fuel_cost_escalation_pct,
            chp_fuel_cost_escalation_pct,
            generator_fuel_cost_escalation_pct,
            offtaker_tax_pct,
            offtaker_discount_pct,
            third_party_ownership,
            owner_tax_pct,
            owner_discount_pct,
            analysis_years,
            value_of_lost_load_per_kwh,
            microgrid_upgrade_cost_pct,
            macrs_five_year,
            macrs_seven_year,
            offgrid_other_capital_costs,
            offgrid_other_annual_costs
        )
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
export AbstractMeasData, MeasData
export meastime, samplingrate, measdata

abstract type AbstractMeasData end

"""
`MeasData`

Structure to store data acquired from a DAQ device. It also stores metadata related to
the DAQ device, data acquisition process and daq channels.
"""
struct MeasData{T} <: AbstractMeasData
    "Device that generated the data"
    devname::String
    "Type of device"
    devtype::String
    "Time of data aquisition"
    time::DateTime
    "Sampling rate (Hz)"
    rate::Float64
    "Data acquired"
    data::T
    "Index of each channel"
    chans::OrderedDict{String,Int}
end

#DateTime(Dates.UTInstant(Millisecond(d.t)))
"Convert a DateTime object to ms"
time2ms(t::DateTime) = t.instant.periods.value
"Convert a time in ms to DateTimeObject"
ms2time(ms::Int64) = DateTime(Dates.UTInstant{Millisecond}(Millisecond(ms)))

"Device name that acquired the data"
devname(d::MeasData) = d.devname

"Device type that acquired the data"
devtype(d::MeasData) = d.devtype

"When did the data acquisition take place?"
meastime(d::AbstractMeasData) = d.time

"What was the sampling rate of the data acquisition?"
samplingrate(d::MeasData) = d.rate

"Access to the data acquired"
measdata(d::MeasData) = d.data
daqchannels(d::MeasData) = collect(keys(d.chans))

import Base.getindex

"Access the data in channel name `ch`"
getindex(d::MeasData{T},ch::AbstractString) where {T<:AbstractMatrix}=view(d.data,d.chans[ch],:)
"Access the data in channel index `i`"
getindex(d::MeasData{T}, i::Integer) where {T<:AbstractMatrix} = view(d.data, i, :)
"Access the data in channel index `i` at time index `k`"
getindex(d::MeasData{T}, i::Integer,k::Integer) where {T<:AbstractMatrix}= d.data[i,k]
"Access the data in channel name `ch` at time index `k`"
getindex(d::MeasData{T}, ch::AbstractString,k::Integer) where {T<:AbstractMatrix}= d.data[d.chans[ch],k]

getindex(d::MeasData) = d.data


text:
#=
bisect_left:
- Julia version: 1.0
- Author: qz
- Date: 2018-08-26
=#

using Test

using BisectPy: bisect_left

@test bisect_left([1, 2, 3, 4, 5], 3.5) == 4

@test bisect_left([1, 2, 3, 4, 5], 2) == 2

@test bisect_left([1, 2, 3, 3, 3, 5], 3) == 3


text:
## ExoplanetsSysSim/src/transit_observations.jl
## (c) 2015 Eric B. Ford

#using Distributions
#include("constants.jl")

#  Starting Section of Observables that are actually used
immutable TransitPlanetObs
  # ephem::ephemeris_type     # For now hardcode P and t0, see transit_observation_unused.jl to reinstate
  period::Float64             # days
  t0::Float64                 # days
  depth::Float64              # fractional
  duration::Float64           # days; Full-width, half-max-duration until further notice
  # ingress_duration::Float64   # days;  QUERY:  Will we want to use the ingress/egress duration for anything?
end
TransitPlanetObs() = TransitPlanetObs(0.0,0.0,0.0,0.0)

immutable StarObs
  radius::Float64      # in Rsol
  mass::Float64        # in Msol
end

semimajor_axis(P::Float64, M::Float64) = (grav_const/(4pi^2)*M*P*P)^(1/3)

function semimajor_axis(ps::PlanetarySystemAbstract, id::Integer)
  M = mass(ps.star) + ps.planet[id].mass   # TODO SCI DETAIL: Replace with Jacobi mass?  Probably not important unless start including TTVs at some point
  @assert(M>0.0)
  @assert(ps.orbit[id].P>0.0)
  return semimajor_axis(ps.orbit[id].P,M)
end

function calc_transit_depth(t::KeplerTarget, s::Integer, p::Integer)  # WARNING: IMPORTANT: Assumes non-grazing transit & no limb darkening
  depth = (t.sys[s].planet[p].radius/t.sys[s].star.radius)^2 # TODO SCI DETAIL: Include limb darkening?
  depth *=  flux(t.sys[s].star)/flux(t)                      # Flux ratio accounts for dillution
end

function calc_transit_duration_central_circ(ps::PlanetarySystemAbstract, pl::Integer)
  duration = rsol_in_au*ps.star.radius * ps.orbit[pl].P /(pi*semimajor_axis(ps,pl) )
end
calc_transit_duration_central_circ(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central_circ(t.sys[s],p)

function calc_transit_duration_central(ps::PlanetarySystemAbstract, pl::Integer)
  ecc = ps.orbit[pl].ecc
  vel_fac = sqrt((1+ecc)*(1-ecc))/(1+ecc*sin(ps.orbit[pl].omega))
  duration = calc_transit_duration_central_circ(ps,pl) * vel_fac
end
calc_transit_duration_central(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_duration_central(t.sys[s],p)

function calc_transit_duration(ps::PlanetarySystemAbstract, pl::Integer)
  a = semimajor_axis(ps,pl)
  @assert a>=zero(a)
  ecc = ps.orbit[pl].ecc
  @assert zero(ecc)<=ecc<=one(ecc)
  b = (a*abs(cos(ps.orbit[pl].incl))/(ps.star.radius*rsol_in_au)) * (1+ecc)*(1-ecc)/(1+ecc*sin(ps.orbit[pl].omega))
  @assert !isnan(b)
  @assert zero(b)<=b
  vel_fac = sqrt((1+ecc)*(1-ecc))/(1+ecc*sin(ps.orbit[pl].omega))
  duration = 0.0<=b<1.0 ? calc_transit_duration_central_circ(ps,pl) * vel_fac * sqrt((1-b)*(1+b)) : 0.0
end
calc_transit_duration(t::KeplerTarget, s::Integer, p::Integer ) = calc_transit_duration(t.sys[s],p)

function calc_expected_num_transits(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
 period = t.sys[s].orbit[p].P
 exp_num_transits = t.duty_cycle * t.data_span/period
 #=
 if exp_num_transits <=6
 end
 # TODO SCI DETAIL: Calculate more accurat number of transits, perhaps using star and specific window function or perhaps specific times of data gaps more given module
 =#
 return exp_num_transits
end

include("transit_detection_model.jl")
include("transit_prob_geometric.jl")

type KeplerTargetObs                        # QUERY:  Do we want to make this type depend on whether the catalog is based on simulated or real data?
  obs::Vector{TransitPlanetObs}
  sigma::Vector{TransitPlanetObs}           # Simplistic approach to uncertainties for now.  QUERY: Should estimated uncertainties be part of Observations type?
  # snr::Vector{Float64}                      # Dimensionless SNR of detection for each planet QUERY: Should we store this here?
  # phys_id::Vector{tuple(Int32,Int32)}       # So we can lookup the system's properties for 0.3.*
  phys_id::Vector{Tuple{Int32,Int32}}     # So we can lookup the system's properties for 0.4

  prob_detect::SystemDetectionProbsAbstract  # QUERY: Specialize type of prob_detect depending on whether for simulated or real data?

  has_sc::Vector{Bool}                      # TODO OPT: Make Immutable Vector or BitArray to reduce memory use?  QUERY: Should this go in KeplerTargetObs?

  star::StarObs                             # TODO SCI DETAIL: Add more members to StarObs, so can used observed rather than actual star properties
end
KeplerTargetObs(n::Integer) = KeplerTargetObs( fill(TransitPlanetObs(),n), fill(TransitPlanetObs(),n), fill(tuple(0,0),n),  ObservedSystemDetectionProbsEmpty(),  fill(false,num_quarters), StarObs(0.0,0.0) )
num_planets(t::KeplerTargetObs) = length(t.obs)

function calc_target_obs_sky_ave(t::KeplerTarget, sim_param::SimParam)
  const max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  const transit_noise_model = get_function(sim_param,"transit_noise_model")
  const min_detect_prob_to_be_included = 0.0  # get_real(sim_param,"min_detect_prob_to_be_included")
  const num_observer_samples = 1 # get_int(sim_param,"num_viewing_geometry_samples")

  np = num_planets(t)
  obs = Array{TransitPlanetObs}(np)
  sigma = Array{TransitPlanetObs}(np)
  #id = Array{Tuple{Int32,Int32}}(np)
  id = Array{Tuple{Int32,Int32}}(np)
  ns = length(t.sys)
  sdp_sys = Array{SystemDetectionProbsAbstract}(ns)
  i = 1
  for (s,sys) in enumerate(t.sys)
    pdet = zeros(num_planets(sys))
    for (p,planet) in enumerate(sys.planet)
      if get(sim_param,"verbose",false)
         println("# s=",s, " p=",p," num_sys= ",length(t.sys), " num_pl= ",num_planets(sys) )
      end
        ntr = calc_expected_num_transits(t, s, p, sim_param)
        period = sys.orbit[p].P
        # t0 = rand(Uniform(0.0,period))   # WARNING: Not being calculated from orbit
        depth = calc_transit_depth(t,s,p)
        duration_central = calc_transit_duration_central(t,s,p)
	snr_central = calc_snr_if_transit(t, depth, duration_central, sim_param, num_transit=ntr)
	pdet_ave = calc_ave_prob_detect_if_transit(t, snr_central, sim_param, num_transit=ntr)
	add_to_catalog = pdet_ave > min_detect_prob_to_be_included  # Include all planets with sufficient detection probability
	if add_to_catalog
	   const hard_max_num_b_tries = 100
	   max_num_b_tries = min_detect_prob_to_be_included == 0. ? hard_max_num_b_tries : min(hard_max_num_b_tries,convert(Int64,1/min_detect_prob_to_be_included))
           pdet_this_b = 0.0
           for j in 1:max_num_b_tries
              b = rand()  # WARNING: Making an approximation: Using a uniform distribution for b (truncated to ensure detection probability >0) when generating measurement uncertainties, rather than accounting for increased detection probability for longer duration transits
              transit_duration_factor = sqrt((1+b)*(1-b))
	      duration = duration_central * transit_duration_factor
	      snr = snr_central * transit_duration_factor
              pdet_this_b = calc_prob_detect_if_transit(t, snr, sim_param, num_transit=ntr)
              if pdet_this_b > 0.0
	         pdet[p] = pdet_ave
                 obs[i], sigma[i] = transit_noise_model(t, s, p, depth, duration, snr, ntr)   # WARNING: noise properties don't have correct dependance on b
                 id[i] = tuple(convert(Int32,s),convert(Int32,p))
      	         i += 1
                 break
              end
           end
	else
	   # Do anything for planets that are extremely unlikely to be detected even if they were to transit?
	end
    end
    resize!(obs,i-1)
    resize!(sigma,i-1)
    resize!(id,i-1)
    sdp_sys[s] = calc_simulated_system_detection_probs(sys, pdet, max_tranets_in_sys=max_tranets_in_sys, min_detect_prob_to_be_included=min_detect_prob_to_be_included, num_samples=num_observer_samples)
  end
  # TODO SCI DETAIL: Combine sdp_sys to allow for target to have multiple planetary systems
  s1 = findfirst(x->num_planets(x)>0,sdp_sys)  # WARNING IMPORTANT: For now just take first system with planets
  if s1 == 0
     s1 = 1
  end
  sdp_target = sdp_sys[s1]

  has_no_sc = fill(false,num_quarters)
  star_obs = StarObs( t.sys[1].star.radius, t.sys[1].star.mass )  # TODO SCI DETAIL: Could improve.  WARNING: ASSUMES STAR IS KNOWN PERFECTLY
  return KeplerTargetObs(obs, sigma, id, sdp_target, has_no_sc, star_obs )
end

function calc_target_obs_single_obs(t::KeplerTarget, sim_param::SimParam)
  #const max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  const transit_noise_model = get_function(sim_param,"transit_noise_model")
  const min_detect_prob_to_be_included = 0.0  # get_real(sim_param,"min_detect_prob_to_be_included")

  np = num_planets(t)
  obs = Array{TransitPlanetObs}(np)
  sigma = Array{TransitPlanetObs}(np)
  #id = Array{Tuple{Int32,Int32}}(np)
  id = Array{Tuple{Int32,Int32}}(np)
  ns = length(t.sys)
  #sdp_sys = Array{SystemDetectionProbsAbstract}(ns)
  sdp_sys = Array{ObservedSystemDetectionProbs}(ns)
  i = 1
  for (s,sys) in enumerate(t.sys)
    pdet = zeros(num_planets(sys))
    for (p,planet) in enumerate(sys.planet)
      if get(sim_param,"verbose",false)
         println("# s=",s, " p=",p," num_sys= ",length(t.sys), " num_pl= ",num_planets(sys) )
      end
        duration = calc_transit_duration(t,s,p)
	if duration <= 0.
	   continue
	end
        ntr = calc_expected_num_transits(t, s, p, sim_param)
        period = sys.orbit[p].P
        # t0 = rand(Uniform(0.0,period))   # WARNING: Not being calculated from orbit
        depth = calc_transit_depth(t,s,p)
	snr = calc_snr_if_transit(t, depth, duration, sim_param, num_transit=ntr)
	pdet[p] = calc_prob_detect_if_transit(t, snr, sim_param, num_transit=ntr)
	if pdet[p] > min_detect_prob_to_be_included
           obs[i], sigma[i] = transit_noise_model(t, s, p, depth, duration, snr, ntr)
           id[i] = tuple(convert(Int32,s),convert(Int32,p))
      	   i += 1
	end
    end
    resize!(obs,i-1)
    resize!(sigma,i-1)
    resize!(id,i-1)
    sdp_sys[s] = ObservedSystemDetectionProbs(pdet)
  end
  # TODO SCI DETAIL: Combine sdp_sys to allow for target to have multiple planetary systems
  s1 = findfirst(x->num_planets(x)>0,sdp_sys)  # WARNING: For now just take first system with planets, assumes not two stars wht planets in one target
  if s1 == 0
     s1 = 1
  end
  sdp_target = sdp_sys[s1]

  has_no_sc = fill(false,num_quarters)
  star_obs = StarObs( t.sys[1].star.radius, t.sys[1].star.mass )  # TODO SCI DETAIL: Could improve.  WARNING: ASSUMES STAR IS KNOWN PERFECTLY
  return KeplerTargetObs(obs, sigma, id, sdp_target, has_no_sc, star_obs )
end

function test_transit_observations(sim_param::SimParam; verbose::Bool=false)  # TODO TEST: Add more tests
  #transit_param = TransitParameter( EphemerisLinear(10.0, 0.0), TransitShape(0.01, 3.0/24.0, 0.5) )
  generate_kepler_target = get_function(sim_param,"generate_kepler_target")
  const max_it = 100000
  local obs
  for i in 1:max_it
    target = generate_kepler_target(sim_param)::KeplerTarget
    while num_planets(target) == 0
      target = generate_kepler_target(sim_param)::KeplerTarget
    end
    calc_transit_prob_single_planet_one_obs(target,1,1)
    calc_transit_prob_single_planet_obs_ave(target,1,1)
    obs = calc_target_obs_single_obs(target,sim_param)
    obs = calc_target_obs_sky_ave(target,sim_param)
    if verbose && (num_planets(obs) > 0)
      println("# i= ",string(i)," np= ",num_planets(obs), " obs= ", obs )
      break
    end
  end
  return obs
end

randtn() = rand(TruncatedNormal(0.0,1.0,-0.999,0.999))

function transit_noise_model_no_noise(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64)
  period = t.sys[s].orbit[p].P
  t0 = rand(Uniform(0.0,period))    # WARNING: Not being calculated from orbit
  sigma_period = 0.0
  sigma_t0 = 0.0
  sigma_depth =  0.0
  sigma_duration =  0.0
  sigma = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration )
  obs = TransitPlanetObs( period, t0, depth,duration)
  return obs, sigma
end

function transit_noise_model_fixed_noise(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Real)
  period = t.sys[s].orbit[p].P
  t0 = rand(Uniform(0.0,period))   # WARNING: Not being calculated from orbit

  sigma_period = 1e-6
  sigma_t0 = 1e-4
  sigma_depth =  0.01
  sigma_duration =  0.01

  sigma = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration)
  #obs = TransitPlanetObs( period, t0, depth, duration)
  obs = TransitPlanetObs( period*(1.0+sigma.period*randtn()), t0*(1.0+sigma.period*randtn()), depth*(1.0+sigma.depth*randtn()),duration*(1.0+sigma.duration*randtn()))
  return obs, sigma
end

function transit_noise_model_diagonal(t::KeplerTarget, s::Integer, p::Integer, depth::Float64, duration::Float64, snr::Float64, num_tr::Real)
  period = t.sys[s].orbit[p].P
  t0 = rand(Uniform(0.0,period))    # WARNING: Not being calculated from orbit

	# Use variable names from Price & Rogers
	one_minus_e2 = (1-t.sys[s].orbit[p].ecc)*(1+t.sys[s].orbit[p].ecc)
	a_semimajor_axis = semimajor_axis(t.sys[s],p)
	b = a_semimajor_axis *cos(t.sys[s].orbit[p].incl)/t.sys[s].star.radius
        b *= one_minus_e2/(1+t.sys[s].orbit[p].ecc*sin(t.sys[s].orbit[p].omega))
	tau0 = t.sys[s].star.radius*period/(a_semimajor_axis*2pi)
	tau0 *= sqrt(one_minus_e2)/(1+t.sys[s].orbit[p].ecc*sin(t.sys[s].orbit[p].omega))
	r = t.sys[s].planet[p].radius/t.sys[s].star.radius
	sqrt_one_minus_b2 = (0.0<=b<1.0) ? sqrt((1-b)*(1+b)) : 0.0
 	T = 2*tau0*sqrt_one_minus_b2
	tau = 2*tau0*r/sqrt_one_minus_b2
	Ttot = period
	I = LC_integration_time      # WARNING: Assumes LC only
	Lambda_eff = LC_rate * num_tr # calc_expected_num_transits(t, s, p, sim_param)
	delta = depth
	sigma = t.cdpp[1,1]

	# Price & Rogers Eqn A8 & Table 1 # WARNING: Someone should check Eqns
	tau3 = tau^3
	I3 = I^3
	#a1 = (10*tau3+2*I^3-5*tau*I^2)/tau3
	a2 = (5*tau3+I3-5*tau*tau*I)/tau3
	a3 = (9*I^5*Ttot-40*tau3+I*I*Ttot+120*tau^4*I*(3*Ttot-2*tau))/tau^6
	a4 = (a3*tau^5+I^4*(54*tau-35*Ttot)-12*tau*I3*(4*tau+Ttot)+360*tau^4*(tau-Ttot))/tau^5
	a5 = a2*(24T*T*(I-3*tau)-24*T*Ttot*(I-3*tau)+tau3*a4)/tau3
	#a6 = (3*tau*tau+T*(I-3*tau))/(tau*tau)
	#a7 = (-60*tau^4+12*a2*tau3*T-9*I^4+8*tau*I3+40*tau3*I)/(tau^4)
	#a8 = (2T-Ttot)/tau
	#a9 = (-3*tau*tau*I*(-10*T*T+10*T*Ttot+I*(2*I+5*Ttot))-I^4*Ttot+8*tau*I3*Ttot)/(tau^5)
	#a10 = ((a9+60)*tau*tau+10*(-9*T*T+9*T*Tot+I*(3*I+Tot))-75*tau*Ttot)/(tau*tau)
	a11 = (I*Ttot-3tau*(Ttot-2*tau))/(tau*tau)
	#a12 = (-360*tau^5-24*a2*tau3*T*(I-3tau)+9*I^5-35tau*I^4-12tau*tau*I3-40tau^3*I*I+360*tau^4*I)/(tau^5)
	a13 = (-3*I3*(8*T*T-8*T*Ttot+3*I*Ttot)+120*tau*tau*T*I*(T-Ttot)+8tau*I3*Ttot)/tau^5
	a14 = (a13*tau*tau+40*(-3*T*T+3*T*Ttot+I*Ttot)-60*tau*Ttot)/tau^5
	a15 = (2I-6tau)/tau
	b1  = (6*I*I-3*I*Ttot+tau*Ttot)/(I*I)
	b2  = (tau*T+3*I*(I-T))/(I*I)
	b3 = (tau3-12T*I*I+8I3+20tau*I*I-8*tau*tau*I)/I3
	b4 = (6*T*T-6*T*Ttot+I*(5Ttot-4I))/(I*I)
	b5 = (10I-3tau)/I
	b6 = (12b4*I3+4tau*(-6*T*T+6T*Ttot+I*(13*Ttot-30I)))/I3
	b7 = (b6*I^5+4*tau*tau*I*I*(12I-11Ttot)+tau3*I*(11Ttot-6I)-tau^4*Ttot)/I^5
	b8 = (3T*T-3*T*Ttot+I*Ttot)/(I*I)
	b9 = (8*b8*I^4+20*tau*I*I*Ttot-8*tau*tau*I*Ttot+tau3*Ttot)/I^4
	#b10 =  (-tau^4+24*T*I*I*(tau-3I)+60*I^4+52*tau*I3-44*tau*tau*I*I+11*tau3*I)/I^4
	#b11 =  (-15b4*i3+10b8*tau*I*I+15*tau*tau*(2I-Ttot))/I3
	#b12 =  (b11*I^5+2tau3*I*(4Ttot-3I)-tau^4*Ttot)/I^5
	#b13 =  (Ttot-2T)/I
	#b14 =  (6I-2tau)/I

        Q = snr/sqrt(num_tr)
        sigma_t0 = tau>=I ?  sqrt(0.5*tau*T/(1-I/(3tau)))/Q : sqrt(0.5*I  *T/(1-tau/(3I)))/Q
        sigma_period = sigma_t0/sqrt(num_tr)
        sigma_duration = tau>=I ? sigma*sqrt(abs(6*tau*a14/(delta*delta*a5)) /Lambda_eff )  : sigma*sqrt(abs(6*I*b9/(delta*delta*b7)) / Lambda_eff)
        sigma_depth = tau>=I ? sigma*sqrt(abs(-24*a11*a2/(tau*a5)) / Lambda_eff)  : sigma*sqrt(abs(24*b1/(I*b7)) / Lambda_eff)

  	sigma = TransitPlanetObs( sigma_period, sigma_t0, sigma_depth, sigma_duration )

        local obs
        if true     # Assume uncertainties uncorrelated (Diagonal)
  	    obs = TransitPlanetObs( period*(1.0+sigma.period*randtn()), t0*(1.0+sigma.period*randtn()), depth*(1.0+sigma.depth*randtn()),duration*(1.0+sigma.duration*randtn()))
        else        # TODO SCI DETAIL:  Account for correlated uncertaintties in transit parameters
            cov = zeros(4,4)
        if tau>=I
	#cov[0,0] = -3tau/(delta*delta*a15)
	cov[1,1] = 24*tau*a10/(delta*delta*a5)
	cov[1,2] = cov[2,1] = 36*a8*tau*a1/(delta*delta*a5)
	cov[1,3] = cov[3,1] = -12*a11*a1/(delta*a5)
	cov[1,4] = cov[4,1] = -12*a6*a1/(delta*a5)
	cov[2,2] = 6*tau*a14/(delta*delta*a5)
	cov[2,3] = cov[3,2] = 72*a8*a2/(delta*a5)
	cov[2,4] = cov[4,2] = 6*a7/(delta*a5)
	cov[3,3] = -24*a11*a2/(tau*a5)
	cov[3,4] = cov[4,3] = -24*a6*a2/(tau*a5)
	cov[4,4] = a12/(tau*a5)
        else
	#cov[0,0] = -3I/(delta*delta*b14)
	cov[1,1] = -24*I*I*b12/(delta*delta*b7)
	cov[1,2] = cov[2,1] = 36*I*b13*b5/(delta*delta*b7)
	cov[1,3] = cov[3,1] = 12*b5*b1/(delta*b7)
	cov[1,4] = cov[4,1] = 12*b5*b2/(delta*b7)
	cov[2,2] = 6*I*b9/(delta*delta*b7)
	cov[2,3] = cov[3,2] = 72*b13/(delta*b7)
	cov[2,4] = cov[4,2] = 6*b3/(delta*b7)
	cov[3,3] = 24*b1/(I*b7)
	cov[3,4] = cov[4,3] = 24*b2/(I*b7)
	cov[4,4] = b10/(I*b7)
	end
	cov .*= sigma*sigma/Lambda_eff
	obs_dist = MvNormal(zeros(4),cov)

	local obs_duration, obs_depth, sigma_duration, sigma_depth
        isvalid = false
	while !isvalid
	  obs_vec = rand(obs_dist)
          obs_duration = duration + obs_vec[2]
          obs_depth = depth + obs_vec[3]
          if (obs_duration>0.0) && (obs_depth>0.0)
             isvalid = true
	  end
	end
     	    obs = TransitPlanetObs( period*(1.0+sigma.period*randtn()), t0*(1.0+sigma.period*randtn()), obs_depth,obs_duration)
        end
  	return obs, sigma
end


text:
module BenchParallelHistogram

using BangBang.NoBang: SingletonDict
using BenchmarkTools
using Transducers
using Transducers: reduce_commutative

suite = BenchmarkGroup()

mergehist!(args...) = merge!(+, args...)
initializer = OnInit(() -> Dict{Float64, Int}())

hist_seq(xf, xs; kwargs...) =
    foldl(mergehist!, xf, xs; init=initializer, kwargs...)

hist_assoc(xf, xs; kwargs...) =
    reduce(mergehist!, xf, xs; init=initializer, kwargs...)

hist_comm(xf, xs; kwargs...) =
    reduce_commutative(mergehist!, xf, xs; init=initializer, kwargs...)

xf = Map() do x
    for _ in 1:10
        x = sin(x)
    end
    return SingletonDict(round(x, digits=4), 1)
end

xs = randn(2^15)

let
    ys_seq = hist_seq(xf, xs)
    ys_assoc = hist_assoc(xf, xs)
    ys_comm = hist_comm(xf, xs)
    @assert ys_seq == ys_assoc
    @assert ys_seq == ys_comm
end

suite["seq"] = @benchmarkable hist_seq($xf, $xs)

suite_assoc = suite["assoc"] = BenchmarkGroup()
for basesize in length(xs) .÷ (2 .^ (0:2) .* Threads.nthreads())
    suite_assoc["basesize=$basesize"] =
        @benchmarkable hist_assoc($xf, $xs; basesize = $basesize)
end

suite_comm = suite["comm"] = BenchmarkGroup()
for basesize in length(xs) .÷ (2 .^ (0:2) .* Threads.nthreads())
    suite_comm["basesize=$basesize"] =
        @benchmarkable hist_comm($xf, $xs; basesize = $basesize)
end

end  # module
BenchParallelHistogram.suite


text:
module CuArrays

using CUDAapi, CUDAdrv, CUDAnative

using GPUArrays

export CuArray, CuVector, CuMatrix, CuVecOrMat, cu
export CUBLAS, CUSPARSE, CUSOLVER, CUFFT, CURAND, CUDNN, CUTENSOR

import LinearAlgebra

using Adapt

using Libdl

using Requires

## source code includes

include("memory.jl")
include("array.jl")
include("subarray.jl")
include("permuteddimsarray.jl")
include("utils.jl")
include("indexing.jl")
include("broadcast.jl")
include("matmul.jl")
include("mapreduce.jl")
include("accumulate.jl")
include("linalg.jl")

include("gpuarray_interface.jl")

include("blas/CUBLAS.jl")
include("sparse/CUSPARSE.jl")
include("solver/CUSOLVER.jl")
include("fft/CUFFT.jl")
include("rand/CURAND.jl")
include("dnn/CUDNN.jl")
include("tensor/CUTENSOR.jl")

include("nnlib.jl")

include("deprecated.jl")

## initialization

const __initialized__ = Ref(false)
functional() = __initialized__[]

export has_cudnn, has_cutensor
has_cudnn() = Libdl.dlopen_e(CUDNN.libcudnn[]) !== C_NULL
has_cutensor() = Libdl.dlopen_e(CUTENSOR.libcutensor[]) !== C_NULL

function __init__()
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false")) || precompiling
    verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))

    # if any dependent GPU package failed, expect it to have logged an error and bail out
    if !CUDAdrv.functional() || !CUDAnative.functional()
        verbose && @warn "CuArrays.jl did not initialize because CUDAdrv.jl or CUDAnative.jl failed to"
        return
    end

    try
        # discover libraries
        toolkit = find_toolkit()
        for name in ("cublas", "cusparse", "cusolver", "cufft", "curand", "cudnn", "cutensor")
            mod = getfield(CuArrays, Symbol(uppercase(name)))
            lib = Symbol("lib$name")
            handle = getfield(mod, lib)

            # on Windows, the library name is version dependent
            if Sys.iswindows()
                cuda = CUDAnative.version()
                suffix = cuda >= v"10.1" ? "$(cuda.major)" : "$(cuda.major)$(cuda.minor)"
                handle[] = "$(name)$(Sys.WORD_SIZE)_$(suffix)"
            end

            # check if we can't find the library
            if Libdl.dlopen_e(handle[]) == C_NULL
                path = find_cuda_library(name, toolkit)
                if path !== nothing
                    handle[] = path
                end
            end
        end

        # library dependencies
        CUBLAS.version()
        CUSPARSE.version()
        CUSOLVER.version()
        CUFFT.version()
        CURAND.version()
        # CUDNN and CUTENSOR are optional

        # library compatibility
        if has_cutensor()
            cutensor = CUTENSOR.version()
            if cutensor < v"1"
                silent || @warn("CuArrays.jl only supports CUTENSOR 1.0 or higher")
            end

            cuda = CUDAnative.version()
            cutensor_cuda = CUDNN.cuda_version()
            if cutensor_cuda.major != cuda.major || cutensor_cuda.minor != cuda.minor
                silent || @warn("You are using CUTENSOR $cutensor for CUDA $cutensor_cuda with CUDA toolkit $cuda; these might be incompatible.")
            end
        end
        if has_cudnn()
            cudnn = CUDNN.version()
            if cudnn < v"7.6"
                silent || @warn("CuArrays.jl only supports CUDNN v7.6 or higher")
            end

            cuda = CUDAnative.version()
            cudnn_cuda = CUDNN.cuda_version()
            if cudnn_cuda.major != cuda.major || cudnn_cuda.minor != cuda.minor
                silent || @warn("You are using CUDNN $cudnn for CUDA $cudnn_cuda with CUDA toolkit $cuda; these might be incompatible.")
            end
        end

        # package integrations
        @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")

        __init_memory__()

        __initialized__[] = true
    catch ex
        # don't actually fail to keep the package loadable
        if !silent
            if verbose
                @error "CuArrays.jl failed to initialize" exception=(ex, catch_backtrace())
            else
                @info "CuArrays.jl failed to initialize and will be unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)"
            end
        end
    end
end

end # module


text:
using Plots
using Statistics

function main()

    x = range(1, 5, length=100)
    l = size(x)
    noise = rand(-2:0.5:2, l)

    y1 = exp.(x) .+ noise
    y2 = exp.(x)
    yerr = std(noise)

    # normal_plot(x, y1, y2, yerr)
    # log_plots(x, y1, y2)
    # subplots(x, y1, y2, yerr)

end

function normal_plot(x, y1, y2, yerr)

    scatter(x, y1, yerr=yerr, markersize=2, markershape=:o, alpha=0.5, label="data")
    plot!(x, y2, color="red", legend=:topleft, label="f(x)", dpi=200)

    annotate!(3, 50, Plots.text("Hello world!", 5, :red))

    xaxis!("x")
    yaxis!("y")
    title!("simple plot")

end

function log_plots(x, y1, y2)

    scatter(x, y1, markersize=1, markershape=:o, label="data")
    plot!(x, y2, xaxis="x", yaxis=("y", :log), legend=:topleft, label="f(x)", dpi=200)

    title!("log plot")

end

function subplots(x, y1, y2, yerr)

    p1 = scatter(x, y1,yerr=yerr, markershape=:o, markersize=1, alpha=0.5, label="data")
    p1 = plot!(x, y2, yaxis="y", legend=:topleft, label="f(x)", title="α data")

    p2 = scatter(x, y1, label="label")
    p2 = plot!(x, y2, xaxis="x", yaxis=("y", :log), legend=:topleft, label="f(x)")

    plot(p1, p2, layout=(2,1), dpi=200)

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
@testset "Allocations" begin
    T = Float64
    sys = Swalbe.SysConst{T}(Lx=25, Ly=26)
    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py, fthermalx, fthermaly = Swalbe.Sys(sys, "CPU", true, T)
    # Struct should do what I expect, hopefully
    @test sys.Lx == 25
    @test sys.Ly == 26
    @test sys.Tmax == 1000
    @test sys.tdump == 100
    @test isa(fout, Array{Float64, 3})
    @test isa(ftemp, Array{Float64, 3})
    @test isa(feq, Array{Float64, 3})
    @test isa(height, Array{Float64, 2})
    @test isa(velx, Array{Float64, 2})
    @test isa(vely, Array{Float64, 2})
    @test isa(vsq, Array{Float64, 2})
    @test isa(pressure, Array{Float64, 2})
    @test isa(dgrad, Array{Float64, 3})
    @test isa(Fx, Array{Float64, 2})
    @test isa(Fy, Array{Float64, 2})
    @test isa(slipx, Array{Float64, 2})
    @test isa(slipy, Array{Float64, 2})
    @test isa(h∇px, Array{Float64, 2})
    @test isa(h∇py, Array{Float64, 2})
    @test isa(fthermalx, Array{Float64, 2})
    @test isa(fthermaly, Array{Float64, 2})

    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py, = Swalbe.Sys(sys, "CPU", false, T)
    # Struct should do what I expect, hopefully
    @test sys.Lx == 25
    @test sys.Ly == 26
    @test isa(fout, Array{Float64, 3})
    @test isa(ftemp, Array{Float64, 3})
    @test isa(feq, Array{Float64, 3})
    @test isa(height, Array{Float64, 2})
    @test isa(velx, Array{Float64, 2})
    @test isa(vely, Array{Float64, 2})
    @test isa(vsq, Array{Float64, 2})
    @test isa(pressure, Array{Float64, 2})
    @test isa(dgrad, Array{Float64, 3})
    @test isa(Fx, Array{Float64, 2})
    @test isa(Fy, Array{Float64, 2})
    @test isa(slipx, Array{Float64, 2})
    @test isa(slipy, Array{Float64, 2})
    @test isa(h∇px, Array{Float64, 2})
    @test isa(h∇py, Array{Float64, 2})

    dyn = Swalbe.Sys(sys, "CPU")
    # Struct should do what I expect, hopefully
    @test sys.Lx == 25
    @test sys.Ly == 26
    @test isa(dyn.fout, Array{Float64, 3})
    @test isa(dyn.ftemp, Array{Float64, 3})
    @test isa(dyn.feq, Array{Float64, 3})
    @test isa(dyn.height, Array{Float64, 2})
    @test isa(dyn.velx, Array{Float64, 2})
    @test isa(dyn.vely, Array{Float64, 2})
    @test isa(dyn.vsq, Array{Float64, 2})
    @test isa(dyn.pressure, Array{Float64, 2})
    @test isa(dyn.dgrad, Array{Float64, 3})
    @test isa(dyn.Fx, Array{Float64, 2})
    @test isa(dyn.Fy, Array{Float64, 2})
    @test isa(dyn.slipx, Array{Float64, 2})
    @test isa(dyn.slipy, Array{Float64, 2})
    @test isa(dyn.h∇px, Array{Float64, 2})
    @test isa(dyn.h∇py, Array{Float64, 2})

    # 1D case
    sys1d = Swalbe.SysConst_1D{T}(L=25)
    fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p, fthermal = Swalbe.Sys(sys1d, true, T)
    # Struct should do what I expect, hopefully
    @test sys1d.L == 25
    @test sys1d.Tmax == 1000
    @test sys1d.tdump == 100
    @test isa(fout, Matrix{Float64})
    @test isa(ftemp, Matrix{Float64})
    @test isa(feq, Matrix{Float64})
    @test isa(height, Vector{Float64})
    @test isa(vel, Vector{Float64})
    @test isa(pressure, Vector{Float64})
    @test isa(dgrad, Matrix{Float64})
    @test isa(F, Vector{Float64})
    @test isa(slip, Vector{Float64})
    @test isa(h∇p, Vector{Float64})
    @test isa(fthermal, Vector{Float64})

    dyn = Swalbe.Sys(sys1d)
    # Struct should do what I expect, hopefully
    @test sys1d.L == 25
    @test sys1d.Tmax == 1000
    @test sys1d.tdump == 100
    @test isa(dyn.fout, Matrix{Float64})
    @test isa(dyn.ftemp, Matrix{Float64})
    @test isa(dyn.feq, Matrix{Float64})
    @test isa(dyn.height, Vector{Float64})
    @test isa(dyn.vel, Vector{Float64})
    @test isa(dyn.pressure, Vector{Float64})
    @test isa(dyn.dgrad, Matrix{Float64})
    @test isa(dyn.F, Vector{Float64})
    @test isa(dyn.slip, Vector{Float64})
    @test isa(dyn.h∇p, Vector{Float64})

    fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p = Swalbe.Sys(sys1d, false, T)
    # Struct should do what I expect, hopefully
    @test sys1d.L == 25
    @test sys1d.Tmax == 1000
    @test sys1d.tdump == 100
    @test isa(fout, Matrix{Float64})
    @test isa(ftemp, Matrix{Float64})
    @test isa(feq, Matrix{Float64})
    @test isa(height, Vector{Float64})
    @test isa(vel, Vector{Float64})
    @test isa(pressure, Vector{Float64})
    @test isa(dgrad, Matrix{Float64})
    @test isa(F, Vector{Float64})
    @test isa(slip, Vector{Float64})
    @test isa(h∇p, Vector{Float64})
end


text:
abstract type AbstractPattern end

indices(p::AbstractPattern) = p.idx

struct AllToAll <: AbstractPattern
    n::Int								# number of vertices
    idx::Vector{Tuple{Int, Int}}		# edges indices
    function AllToAll(n)
        n > 1 ? new(n, vcat([[(i, j) for j in filter(pos -> pos != i, collect(1:n))] for i in collect(1:n)]...)) : error("incorrect length: $n")
    end
end

struct Cycle <: AbstractPattern
    n::Int								# cycle length
    idx::Vector{Tuple{Int, Int}}		# edges indices
    function Cycle(n)
        n >= 1 ? new(n, [(i, (i%n)+1) for i in 1:n]) : error("cannot make a cycle of length $n")
    end
end

struct Path <: AbstractPattern
    n::Int
    idx::Vector{Tuple{Int, Int}}
    function Path(n)
        n > 0 ? new(n, [(i, i+1) for i in collect(1:n-1)]) : error("incorrect length: $n (must be positive)")
    end
end


text:
@testset "randomsubset.jl" begin
    N = 30
    test_Zalg(RandomSubset(10))
    test_Zalg(RandomSubset(10), N; weights=rand(N))
end


text:
using Raytracing
using GeometryBasics
using Random
using CoordinateTransformations
using ColorTypes
using ColorVectorSpace
using StaticArrays
using Images
using ImageIO
using FileIO
using ImageMagick
using Test

@testset "Raytracing.jl" begin
    T = Float64
    origin = Point3{T}(0,0,1.5)
    lookat = Point3{T}(0,0,0)
    vup = Vec3{T}(0,1,0)
    aspectratio = T(16/9)
    vfov = T(40)
    camera = Raytracing.make_camera(origin, lookat, vup, vfov, aspectratio)
    genray = Raytracing.generate_ray(camera, 0.5, 0.5)
    raydir = Vec3{T}(0,0,-1)
    @test Raytracing.origin(genray) ≈ origin
    @test Raytracing.direction(genray) ≈ raydir
    @test genray ≈ Raytracing.Ray(origin, raydir)
    @test ≈(genray, Raytracing.Ray(origin, raydir); atol=1e-2)

    origin2 = Point3{T}(0,0.5,1)
    camera2 = Raytracing.make_camera(origin2, lookat, vup, vfov, aspectratio)
    #@show dump(camera2)
    #@show Raytracing.generate_ray(camera2, 0.5, 0.5)
    #@show Raytracing.generate_ray(camera2, 0.0, 0.0)

    spheresdf = Raytracing.SphereSignedDistanceField{T}(0.5)
    @test spheresdf(origin) ≈ 1.0
    @test spheresdf(lookat) ≈ -0.5
    @test spheresdf(Point3{T}(0,0,0.5)) ≈ 0.0

    #@show Raytracing.normal(spheresdf, origin)
    #@show Raytracing.normalzygote(spheresdf, origin)
    #@test Raytracing.normal(spheresdf, origin) ≈ Raytracing.normalzygote(spheresdf, origin)

    normalmat = Raytracing.NormalShading()

    id = IdentityTransformation()
    spherealbedo = Vec{3, T}(0.8, 0.25, 0.25)

    lambertianmat = Raytracing.Lambertian(spherealbedo)

    sphereshape = Raytracing.Shape(id, spheresdf, lambertianmat)

    skyboxspheresdf = Raytracing.SphereSignedDistanceField{T}(100)
    skyboxsphereshape = Raytracing.Shape(id, skyboxspheresdf, normalmat)

    sceneshapes = (skyboxsphereshape, sphereshape)
    #sceneshapes = @SVector [sphereshape]
    scene = Raytracing.Scene(sceneshapes)

    image_height = 10
    samples_per_pixel = 1
    max_bounces_per_ray = 2
    max_steps_per_bounce = 200
    distance_tolerance = 1e-4
    rng = MersenneTwister(1234)
    rat = 2.0 / 0.70
    len = 0.7
    camera3 = Raytracing.make_camera(Point3{T}(0.0, rat * len, 1/rat * len), lookat, vup, vfov, aspectratio)
    image = Raytracing.raytrace_image(scene, camera3, image_height, samples_per_pixel, max_bounces_per_ray, max_steps_per_bounce, distance_tolerance, rng)

    #savename = "lofi_bistar_sphere_skybox.png"
    #tmpdir = mktempdir()
    #saved_image_path = "$(tmpdir)/$(savename)"
    #save(saved_image_path, image)
    #image_reloaded = load(saved_image_path)

    #reference_image = load("$(savename)")
    #Raytracing.saveimage("$(savename)", image)

    # Regression test
    #@test image_reloaded ≈ reference_image

end


text:
module FromFile

export @from

# This replicates Base.__toplevel__
baremodule __toplevel__
    using Base
end

macro from(path::String, ex)
    esc(from_m(__module__, path, ex))
end

function from_m(m::Module, path::String, ex::Expr)
    ex.head === :using || ex.head === :import || error("expect using/import statement")
    path = abspath(normpath(path))
    loading = Expr(ex.head)
    root = root_module(m)

    if root === Main
        toplevel = __toplevel__
    else
        toplevel_symbol = Symbol("#__toplevel__#")
        if isdefined(root, toplevel_symbol) # package
            toplevel = getfield(root, toplevel_symbol)
        else
            toplevel = Base.eval(root, :(baremodule $toplevel_symbol; using Base; end))
        end
    end

    file_module = load_module(toplevel, root, path)

    for each in ex.args
        each isa Expr || continue

        if each.head === :(:) # using/import A: a, b, c
            m_name = each.args[1].args[1]
            m_name === :(.) && error("cannot load relative module from file")
            push!(loading.args, Expr(:(:), Expr(:., fullname(file_module)..., m_name), each.args[2:end]...) )
        elseif each.head === :(.) # using/import A, B, C
            m_name = each.args[1] # module name
            m_name === :(.) && error("cannot load relative module from file")
            push!(loading.args, Expr(:., fullname(file_module)..., m_name))
        else
            error("invalid syntax $ex")
        end
    end
    return loading
end

function load_module(toplevel::Module, root::Module, path::String)
    if root === Main
        file_module_sym = Symbol(path)
    else
        file_module_sym = Symbol(relpath(path, pathof(root)))
    end

    if isdefined(toplevel, file_module_sym)
        file_module = getfield(toplevel, file_module_sym)
    else
        file_module = Base.eval(toplevel, :(module $(file_module_sym); end))
        Base.include(file_module, path)
    end

    return file_module
end

function root_module(m::Module)
    root = parentmodule(m)
    prev = m
    while root !== prev
        prev = root
        root = parentmodule(root)
    end
    return root
end

end


text:
module LegendrePolynomials

using OffsetArrays

export Pl
export collectPl
export collectPl!
export dnPl
export collectdnPl
export collectdnPl!

function checkdomain(x)
	abs(x) > 1 && throw(DomainError(x,
		"Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))
end

assertnonnegative(l) = (l >= 0 || throw(ArgumentError("l must be >= 0, received " * string(l))))

function checklength(arr, minlength)
	length(arr) >= minlength || throw(ArgumentError(
		"array is not large enough to store all values, require a minimum length of " * string(minlength)))
end

@inline function Pl_recursion(::Type{T}, ℓ, Plm1, Plm2, x) where {T}
	# relation is valid from ℓ = 1
	Pl = ((2ℓ-1) * x * Plm1 - (ℓ-1) * Plm2)/ℓ
	convert(T, Pl)
end

@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, P_n_lm2, x) where {T}
	P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1) - (ℓ - 1) * P_n_lm2)/ℓ
	convert(T, P_n_l)
end

# special cases
# case 1 : l == n, in which case there are no (l-1,n) and (l-2,n) terms
@inline function dPl_recursion(::Type{T}, ℓ, n, ::Nothing, P_nm1_lm1, ::Nothing, x) where {T}
	P_n_l = (2ℓ-1) * P_nm1_lm1
	convert(T, P_n_l)
end
# case 1 : l == n + 1, in which case there's no (l-2,n) term
@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, ::Nothing, x) where {T}
	P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1))/ℓ
	convert(T, P_n_l)
end

polytype(x) = typeof(float(x*x))

"""
	LegendrePolynomialIterator(x, [lmax::Integer])

Return an iterator that generates the values of the Legendre polynomials ``P_\\ell(x)`` for the given `x`.
If `lmax` is specified then only the values of ``P_\\ell(x)`` from `0` to `lmax` are returned.

# Examples
```jldoctest
julia> import LegendrePolynomials: LegendrePolynomialIterator

julia> iter = LegendrePolynomialIterator(0.5, 4);

julia> collect(iter)
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
  1.0
  0.5
 -0.125
 -0.4375
 -0.2890625

julia> iter = LegendrePolynomialIterator(0.5);

julia> collect(Iterators.take(iter, 5)) # evaluete 5 elements (l = 0:4)
5-element Array{Float64,1}:
  1.0
  0.5
 -0.125
 -0.4375
 -0.2890625

julia> collect(Iterators.take(Iterators.drop(iter, 100), 5)) # evaluate Pl for l = 100:104
5-element Array{Float64,1}:
 -0.0605180259618612
  0.02196749072249231
  0.08178451892628381
  0.05963329258495025
 -0.021651535258316177
```
"""
struct LegendrePolynomialIterator{T, L <: Union{Integer, Nothing}, V}
	x :: V
	lmax :: L
	function LegendrePolynomialIterator{T,L,V}(x::V, lmax::L) where {T, L <: Union{Integer, Nothing}, V}
		checkdomain(x)
		new{T,L,V}(x, lmax)
	end
end
LegendrePolynomialIterator(x) = LegendrePolynomialIterator{polytype(x), Nothing, typeof(x)}(x, nothing)
function LegendrePolynomialIterator(x, lmax)
	assertnonnegative(lmax)
	LegendrePolynomialIterator{polytype(x), typeof(lmax), typeof(x)}(x, lmax)
end

Base.eltype(::Type{<:LegendrePolynomialIterator{T}}) where {T} = T
Base.IteratorSize(::Type{<:LegendrePolynomialIterator{<:Any,Nothing}}) = Base.IsInfinite()

function Base.iterate(::LegendrePolynomialIterator{T}) where {T}
	return one(T), (1, (zero(T), one(T)))
end
function Base.iterate(iter::LegendrePolynomialIterator{T,Nothing}, state) where {T}
	l, (Plm2, Plm1) = state
	x = iter.x
	Pl = Pl_recursion(T, l, Plm1, Plm2, x)
	return Pl, (l+1, (Plm1, Pl))
end
function Base.iterate(iter::LegendrePolynomialIterator{T}, state) where {T}
	l, (Plm2, Plm1) = state
	l > iter.lmax && return nothing
	x = iter.x
	Pl = Pl_recursion(T, l, Plm1, Plm2, x)
	return Pl, (l+1, (Plm1, Pl))
end

Base.length(iter::LegendrePolynomialIterator{<:Any,<:Integer}) = iter.lmax + 1
Base.size(iter::LegendrePolynomialIterator{<:Any,<:Integer}) = (length(iter),)
Base.axes(iter::LegendrePolynomialIterator{<:Any,<:Integer}) = (0:iter.lmax, )
Base.keys(iter::LegendrePolynomialIterator{<:Any,<:Integer}) = 0:iter.lmax

function Base.collect(iter::LegendrePolynomialIterator{T,<:Integer}) where {T}
	arr = zeros(T, 0:iter.lmax)
	collectPl!(arr, iter.x; lmax = iter.lmax)
end

Base.copy(iter::LegendrePolynomialIterator) = typeof(iter)(iter.x, iter.lmax)

"""
	Pl(x, l::Integer)

Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and the degree `l`

# Examples
```jldoctest
julia> Pl(1, 2)
1.0

julia> Pl(0.5, 20)
-0.04835838106737356
```
"""
function Pl(x, l::Integer)
	iter = LegendrePolynomialIterator(x)
	d = Iterators.drop(iter, l) # 0 to l-1
	first(d)
end

function doublefactorial(T, n)
	p = convert(T, one(n))
	for i in n:-2:1
		p *= convert(T, i)
	end
	convert(T, p)
end

function _checkvalues(x, l, n)
	assertnonnegative(l)
	checkdomain(x)
	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
end

Base.@propagate_inbounds function _unsafednPl!(cache, x, l, n)
	# unsafe, assumes 1-based indexing
	checklength(cache, l - n + 1)
	if n == l # may short-circuit this
		cache[1] = doublefactorial(eltype(cache), 2l-1)
	else
		collectPl!(cache, x, lmax = l - n)

		for ni in 1:n
			# We denote the terms as P_ni_li

			# li == ni
			P_nim1_nim1 = cache[1]
			P_ni_ni = dPl_recursion(eltype(cache), ni, ni, nothing, P_nim1_nim1, nothing, x)
			cache[1] = P_ni_ni

			# li == ni + 1
			P_nim1_ni = cache[2]
			P_ni_nip1 = dPl_recursion(eltype(cache), ni + 1, ni, P_ni_ni, P_nim1_ni, nothing, x)
			cache[2] = P_ni_nip1

			for li in ni+2:min(l, l - n + ni)
				P_ni_lim2 = cache[li - ni - 1]
				P_ni_lim1 = cache[li - ni]
				P_nim1_lim1 = cache[li - ni + 1]
				P_ni_li = dPl_recursion(eltype(cache), li, ni, P_ni_lim1, P_nim1_lim1, P_ni_lim2, x)
				cache[li - ni + 1] = P_ni_li
			end
		end
	end
	nothing
end

"""
	dnPl(x, l::Integer, n::Integer, [cache::AbstractVector])

Compute the ``n``-th derivative ``d^n P_\\ell(x)/dx^n`` of the Legendre polynomial ``P_\\ell(x)``.
Optionally a pre-allocated vector `cache` may be provided, which must have a minimum length of `l - n + 1`
and may be overwritten during the computation.

The order of the derivative `n` must be non-negative. For `n == 0` this function just returns
Legendre polynomials.

# Examples

```jldoctest
julia> dnPl(0.5, 3, 2) # second derivative of P3(x) at x = 0.5
7.5

julia> dnPl(0.5, 4, 0) == Pl(0.5, 4) # zero-th order derivative == Pl(x)
true
```
"""
Base.@propagate_inbounds function dnPl(x, l::Integer, n::Integer,
	A = begin
		_checkvalues(x, l, n)
		# do not allocate A if the value is trivially zero
		if l < n
			return zero(polytype(x))
		end
		zeros(polytype(x), l - n + 1)
	end
	)

	_checkvalues(x, l, n)
	# check if the value is trivially zero in case A is provided in the function call
	if l < n
		return zero(eltype(A))
	end

	cache = OffsetArrays.no_offset_view(A)
	# function barrier, as no_offset_view may be type-unstable
	_unsafednPl!(cache, x, l, n)

	return cache[l - n + 1]
end

"""
	collectPl!(v::AbstractVector, x; [lmax::Integer = length(v) - 1])

Compute the Legendre Polynomials ``P_\\ell(x)`` for the argument `x` and all degrees `l` in `0:lmax`,
and store them in `v`.

At output `v[firstindex(v) + l] == Pl(x,l)`.

# Examples
```jldoctest
julia> v = zeros(4);

julia> collectPl!(v, 0.5)
4-element Array{Float64,1}:
  1.0
  0.5
 -0.125
 -0.4375

julia> v = zeros(0:4);

julia> collectPl!(v, 0.5, lmax = 3) # only l from 0 to 3 are populated
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
  1.0
  0.5
 -0.125
 -0.4375
  0.0
```
"""
function collectPl!(v::AbstractVector, x; lmax::Integer  = length(v) - 1)
	checklength(v, lmax + 1)

	iter = LegendrePolynomialIterator(x, lmax)

	for (l, Pl) in pairs(iter)
		v[l + firstindex(v)] = Pl
	end

	v
end

"""
	collectPl(x; lmax::Integer)

Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l` in `0:lmax`.
Return `P` with indices `0:lmax`, where `P[l] == Pl(x,l)`

# Examples
```jldoctest
julia> collectPl(0.5, lmax = 4)
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
  1.0
  0.5
 -0.125
 -0.4375
 -0.2890625
```
"""
collectPl(x; lmax::Integer) = collect(LegendrePolynomialIterator(x, lmax))

"""
	collectdnPl(x; n::Integer, lmax::Integer)

Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`.

The order of the derivative `n` must be greater than or equal to zero.

Returns `v` with indices `0:lmax`, where `v[l] == dnPl(x, l, n)`.

# Examples

```jldoctest
julia> collectdnPl(0.5, lmax = 3, n = 2)
4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
 0.0
 0.0
 3.0
 7.5
```
"""
function collectdnPl(x; lmax::Integer, n::Integer)
	assertnonnegative(lmax)
	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
	v = zeros(polytype(x), 0:lmax)
	if lmax >= n
		collectdnPl!(parent(v), x; lmax = lmax, n = n)
	end
	v
end

"""
	collectdnPl!(v::AbstractVector, x; lmax::Integer, n::Integer)

Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`,
and store the result in `v`.

The order of the derivative `n` must be greater than or equal to zero.

At output, `v[l + firstindex(v)] == dnPl(x, l, n)` for `l = 0:lmax`.

# Examples

```jldoctest
julia> v = zeros(4);

julia> collectdnPl!(v, 0.5, lmax = 3, n = 2)
4-element Array{Float64,1}:
 0.0
 0.0
 3.0
 7.5
```
"""
function collectdnPl!(v, x; lmax::Integer, n::Integer)
	assertnonnegative(lmax)
	n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
	checklength(v, lmax + 1)

	# trivially zero for l < n
	fill!((@view v[(0:n-1) .+ firstindex(v)]), zero(eltype(v)))
	# populate the other elements
	@inbounds dnPl(x, lmax, n, @view v[(n:lmax) .+ firstindex(v)])

	v
end

end