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

using Behavior.Gherkin: parsescenario!, issuccessful, Given, When, Then, ByLineParser, ScenarioStep

@testset "Scenario descriptions" begin
    @testset "Scenario long description, one line; Long description is available" begin
        text = """
        Scenario: Some description
            This is a longer description
            Given some precondition

            When some action

            Then some postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description"
    end

    @testset "Scenario long description, three lines; Long description is available" begin
        text = """
        Scenario: Some description
            This is a longer description
            This is another line.
            This is a third line.
            Given some precondition

            When some action

            Then some postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\nThis is another line.\nThis is a third line."
    end

    @testset "Scenario with blank lines between description and steps; Trailing blank lines are ignored" begin
        text = """
        Scenario: Some description
            This is a longer description
            This is another line.
            This is a third line.

            Given some precondition

            When some action

            Then some postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\nThis is another line.\nThis is a third line."
    end

    @testset "Scenario with blank lines in description; Blank lines in description are included" begin
        text = """
        Scenario: Some description
            This is a longer description

            This is another line.
            This is a third line.

            Given some precondition

            When some action

            Then some postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\n\nThis is another line.\nThis is a third line."
    end

    @testset "Long description without steps; Zero steps and a long description" begin
        text = """
        Scenario: Some description
            NotAStep some more text
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)
        @test issuccessful(result)
        scenario = result.value

        @test isempty(scenario.steps)
        @test scenario.long_description == "NotAStep some more text"
    end
end

@testset "Outline descriptions " begin
    @testset "Scenario long description, one line; Long description is available" begin
        text = """
        Scenario Outline: Some description
            This is a longer description
            Given some precondition <foo>

            When some action

            Then some postcondition

        Examples:
            | foo |
            |  1  |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description"
    end

    @testset "Scenario long description, three lines; Long description is available" begin
        text = """
        Scenario Outline: Some description
            This is a longer description
            This is another line.
            This is a third line.
            Given some precondition <foo>

            When some action

            Then some postcondition

        Examples:
            | foo |
            |  1  |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\nThis is another line.\nThis is a third line."
    end

    @testset "Scenario with blank lines between description and steps; Trailing blank lines are ignored" begin
        text = """
        Scenario Outline: Some description
            This is a longer description
            This is another line.
            This is a third line.

            Given some precondition <foo>

            When some action

            Then some postcondition

        Examples:
            | foo |
            |  1  |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\nThis is another line.\nThis is a third line."
    end

    @testset "Scenario with blank lines in description; Blank lines in description are included" begin
        text = """
        Scenario Outline: Some description
            This is a longer description

            This is another line.
            This is a third line.

            Given some precondition <foo>

            When some action

            Then some postcondition

        Examples:
            | foo |
            |  1  |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.long_description == "This is a longer description\n\nThis is another line.\nThis is a third line."
    end

    # @testset "Long description without steps; Undefined" begin
    #     # We have problems supporting empty Scenario Outlines because I don't want
    #     # to make "Examples" an unallowed word in scenario descriptions.
    #     # If this turns out to be a problem, we'll have to come up with a solution.
    #     text = """
    #     Scenario Outline: Some description
    #         NotAStep some more text
    #
    #     Examples:
    #         | foo |
    #         |  1  |
    #     """
    #
    #     byline = ByLineParser(text)
    #     result = parsescenario!(byline)
    # end
end


text:
using Vulkan
using BenchmarkTools

# Benchmarks:
# - Creation of minimalist create info structs:
#   - directly
#   - from the high-level struct passed as a literal
#   - via a conversion from an existing high-level struct
# - Creation of high-level structs

println("InstanceCreateInfo")
@btime _InstanceCreateInfo([], []; application_info = _ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))
@btime _InstanceCreateInfo(InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine")))
ci = InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))
@btime _InstanceCreateInfo($ci)
@btime InstanceCreateInfo([], []; application_info = ApplicationInfo(v"1.0", v"1.5", v"1.2"; application_name="MyApp", engine_name="MyEngine"))

println("DescriptorPoolCreateInfo")
@btime _DescriptorPoolCreateInfo(3, _DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))
@btime _DescriptorPoolCreateInfo(DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6])))
ci = DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))
@btime _DescriptorPoolCreateInfo($dp)
@btime DescriptorPoolCreateInfo(3, DescriptorPoolSize.([DESCRIPTOR_TYPE_SAMPLER, DESCRIPTOR_TYPE_STORAGE_BUFFER], [4, 6]))

println("DeviceCreateInfo")
@btime _DeviceCreateInfo([_DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])
@btime _DeviceCreateInfo(DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], []))
ci = DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])
@btime _DeviceCreateInfo($ci)
@btime DeviceCreateInfo([DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])], [], [])

ext: jl


text:
# given a model and a new set of bounds check
# how they will affect a given obj reaction
function check_newbounds(S, b, lb, ub, newlb, newub,
        check_obj::Int, idxs = eachindex(lb);
        check_obj_atol = 1e-4,
        verbose = true,
        batchlen = 50,
        solver = Clp.Optimizer
    )

    lp_model = build_lp_model(S, b, lb, ub, solver)
    fbaout = fba(lp_model, check_obj)
    isempty(fbaout) && error("fba returns is empty")
    ref_obj_val = objval(fbaout)

    # working copy
    wlb, wub = deepcopy.([lb, ub])

    icount = length(idxs)
    batchlen = max(1, min(batchlen, length(idxs)))
    batches = [idxs[i0:(min(i0 + batchlen - 1, icount))] for i0 in 1:batchlen:icount]
    verbose && (prog = Progress(length(batches); desc = "Checking bounds  "))
    for batch in batches

        # Test whole batch (this use the heuristic that only a few rxns will affect the biomass)
        wlb[batch] = @view newlb[batch]
        wub[batch] = @view newub[batch]

        # (TODO: This can be improved by updating only the relevant lp_model bounds)
        # Update bounds
        set_lb_con!(lp_model, wlb)
        set_ub_con!(lp_model, wub)

        # fba
        fbaout = fba(lp_model, check_obj)
        new_obj_val = objval(fbaout)
        if !isapprox(new_obj_val, ref_obj_val; atol = check_obj_atol)

            # reset batch
            wlb[batch] = @view lb[batch]
            wub[batch] = @view ub[batch]

            # If whole batch fail, test independent idxs
            for idx in batch

                # check both first
                wlb[idx], wub[idx] = newlb[idx], newub[idx]

                # Update bounds
                set_lb_con!(lp_model, wlb)
                set_ub_con!(lp_model, wub)

                # fba
                fbaout = fba(lp_model, check_obj)
                new_obj_val = objval(fbaout)

                if !isapprox(new_obj_val, ref_obj_val; atol = check_obj_atol)
                    wlb[idx], wub[idx] = lb[idx], ub[idx] # reset both

                    # if obj_val changed I check each bound
                    for (wcol, newb, oldb) in [
                            (wlb, newlb[idx], lb[idx]),
                            (wub, newub[idx], ub[idx])
                        ]
                        wcol[idx] = newb

                        # Update bounds
                        set_lb_con!(lp_model, wlb)
                        set_ub_con!(lp_model, wub)

                        # fba
                        fbaout = fba(lp_model, check_obj)
                        new_obj_val = objval(fbaout)
                        if !isapprox(new_obj_val, ref_obj_val; atol = check_obj_atol)
                            wcol[idx] = oldb # reset
                        end
                    end
                end
            end
        end

        verbose && next!(prog; showvalues = [
                (:ref_obj_val, ref_obj_val),
                (:curr_obj_val, new_obj_val)
            ]
        )

    end # for batch in batches

    verbose && finish!(prog)

    return wlb, wub
end


text:
using Oceananigans
using Oceananigans.Units
using Printf

@info "Building a grid..."

Ny = 128
Ly = 0.3

Nx = Int(Ny/2)
Lx = Ly/2

Nz = 2Ny # number of points in the vertical direction
Lz = Ly # domain depth

refinement = 1.5 # controls spacing near surface (higher means finer spaced)
stretching = 8   # controls rate of stretching at bottom
h(k) = (k - 1) / Nz
ζ₀(k) = 1 + (h(k) - 1) / refinement
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

arch = GPU()

grid = VerticallyStretchedRectilinearGrid(architecture = arch,
                                          size = (Nx, Ny, Nz),
                                          halo = (3, 3, 3),
                                          x = (0, Lx),
                                          y = (0, Ly),
                                          z_faces = k -> Lz * (ζ₀(k) * Σ(k) - 1),
                                          topology = (Periodic, Bounded, Bounded))
#=
grid = RegularRectilinearGrid(size = (Nx, Ny, Nz),
                              halo = (3, 3, 3),
                              x = (0, Lx),  # longer in streamwise direction
                              y = (0, Ly),  # wide enough to avoid finite-width effects?
                              z = (-Lz, 0), # not full depth, but deep enough?
                              topology = (Periodic, Bounded, Bounded))
=#

@show grid

#####
##### Parameters
#####

a = 0.000      # m, reference wave amplitude at t ≈ τᵃ
τ = 1e-4       # m² s⁻², surface stress
k = 2π / 0.02  # m⁻¹, constant wavenumber
ν = 1.05e-6    # m² s⁻¹, kinematic viscosity
κ = 1e-6       # tracer diffusivity

struct ConstantStokesShear{T}
    a :: T
    k :: T
    ω :: T
end

ConstantStokesShear(a, k; g = 9.81, T = 7.2e-5) =
    ConstantStokesShear{Float64}(a, k, sqrt(g * k + T * k^3))

(sh::ConstantStokesShear)(z, t) = 2 * sh.a^2 * sh.k^2 * sh.ω * exp(2 * sh.k * z)

∂z_uˢ = ConstantStokesShear(a, k)

# Calculations
ω = ∂z_uˢ.ω
ϵ = a * k
u★ = sqrt(abs(τ))
La = k * ν^(3/2) / (a * u★ * √(ω))

@info """

    Wave parameters | Values
    =============== | ======
                  a | $a
                  k | $k
             2π / k | $(2π / k)
                  ϵ | $ϵ
                 La | $La
               La⁻¹ | $(1 / La)
"""

#u_top_bc = FluxBoundaryCondition((x, y, t) -> - 1e-5 * sqrt(t))
u_top_bc = FluxBoundaryCondition(-τ)
u_bcs = FieldBoundaryConditions(top = u_top_bc)
boundary_conditions = (; u = u_bcs)

@info "Modeling..."

model = NonhydrostaticModel(architecture = arch,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            grid = grid,
                            tracers = :c,
                            boundary_conditions = boundary_conditions,
                            closure = IsotropicDiffusivity(ν=ν, κ=κ),
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            coriolis = nothing,
                            buoyancy = nothing)

Δz = try
    grid.Δz # for a regular grid
catch
    minimum(parent(grid.Δzᵃᵃᶜ)) # for a stretched grid
end

uᵢ(x, y, z) = 1e-9 * randn()
set!(model, u=uᵢ, v=uᵢ, w=uᵢ)

c = model.tracers.c
view(interior(c), :, :, grid.Nz) .= 1

@info "Revvving up a simulation..."

simulation = Simulation(model, Δt=0.01, stop_time=1minutes)

wizard = TimeStepWizard(cfl=0.5, max_Δt=1.0)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

function progress(sim)

    umax = maximum(abs, sim.model.velocities.u)
    vmax = maximum(abs, sim.model.velocities.v)
    wmax = maximum(abs, sim.model.velocities.w)

    t = time(sim)
    h = √(ν * t)
    Re = umax * h / ν

    @info @sprintf("Time: %s, iteration: %d, next Δt: %s, max|U|: (%.2e, %.2e, %.2e)  m s⁻¹, hk: %.2f, Re: %.1f",
                   prettytime(t),
                   iteration(sim),
                   prettytime(sim.Δt),
                   umax, vmax, wmax,
                   h * k,
                   Re)

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

#####
##### Set up output
#####

prefix = @sprintf("simple_instability_%d_%d_%d_τ%.1e_a%.1e", grid.Nx, grid.Ny, grid.Nz, τ, a)

outputs = merge(model.velocities, model.tracers)

simulation.output_writers[:yz] = JLD2OutputWriter(model, outputs,
                                                    schedule = TimeInterval(0.1),
                                                    force = true,
                                                    prefix = prefix * "_yz",
                                                    field_slicer = FieldSlicer(i = round(Int, grid.Nx/2)))

simulation.output_writers[:xz] = JLD2OutputWriter(model, outputs,
                                                  schedule = TimeInterval(0.1),
                                                  force = true,
                                                  field_slicer = FieldSlicer(j = round(Int, grid.Ny/2)),
                                                  prefix = prefix * "_xz")

simulation.output_writers[:netcdf_xz] = NetCDFOutputWriter(model, outputs,
                                                           schedule = TimeInterval(0.1),
                                                           mode = "c",
                                                           field_slicer = FieldSlicer(j = round(Int, grid.Ny/2)),
                                                           filepath = prefix * "_xz.nc")

simulation.output_writers[:z] = NetCDFOutputWriter(model, outputs,
                                                   schedule = TimeInterval(0.01),
                                                   mode = "c",
                                                   field_slicer = FieldSlicer(i = round(Int, grid.Nx/2),
                                                                              j = round(Int, grid.Ny/2)),
                                                   filepath = prefix * "_z.nc")

C = AveragedField(model.tracers.c, dims=(1, 2))
U = AveragedField(model.velocities.u, dims=(1, 2))

simulation.output_writers[:averages] = NetCDFOutputWriter(model, (c=C, u=U),
                                                          schedule = TimeInterval(0.1),
                                                          mode = "c",
                                                          filepath = prefix * "_averages.nc")

@time run!(simulation)

@info "Simulation complete: $simulation. Output:"

for (name, writer) in simulation.output_writers
    absfilepath = abspath(writer.filepath)
    @info "OutputWriter $name, $absfilepath:\n $writer"
end


text:
__precompile__()

module GSW

deps = joinpath(dirname(dirname(pathof(GSW))), "deps", "deps.jl")
isfile(deps) ? include(deps) : error("GSW is not properly installed")

include("gen_gswteos_h.jl")
include("gen_gswteos10.jl")
# package code goes here

end # module


text:
add DataFrames Distributions StatsBase GLM CategoricalArrays Dates NLopt
add ForwardDiff Formatting TexTables DataStructures Atom Juno
add LinearAlgebra Gadfly HCubature HTTP Measures Formatting ODBC
add CodecZlib CSV GZip NLopt CategoricalArrays Serialization Compose
add Distributed Random Cairo Fontconfig Compose uCSV Revise StaticArrays

add DataFrames Distributions StatsBase GLM CategoricalArrays Dates NLopt ForwardDiff Formatting TexTables DataStructures Atom Juno LinearAlgebra Gadfly HCubature HTTP Measures Formatting ODBC CodecZlib CSV GZip NLopt CategoricalArrays Serialization Compose Distributed Random Cairo Fontconfig Compose uCSV Revise StaticArrays

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module _DistEn
export DistEn
using StatsBase: Histogram, fit, skewness
    """
        Dist, Ppi = DistEn(Sig)

    Returns the distribution entropy estimate (`Dist`) and the
    corresponding distribution probabilities (`Ppi`) estimated from
    the data sequence (`Sig`) using the default  parameters:
    embedding dimension = 2, time delay = 1, binning method = 'Sturges',
    logarithm = base 2, normalisation = w.r.t # of histogram bins

        Dist, Ppi = DistEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, Bins::Union{Int,String}="Sturges", Logx::Real=2, Norm::Bool=true)

    Returns the distribution entropy estimate (`Dist`) estimated from
    the data sequence (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer  \n
    `tau`   - Time Delay, a positive integer  \n
    `Bins`  - Histogram bin selection method for distance distribution,
              one of the following:  \n
              an integer > 1 indicating the number of bins, or one of the
              following strings {'sturges','sqrt','rice','doanes'}
              [default: 'sturges']
    `Logx`  - Logarithm base, a positive scalar (enter 0 for natural log)   \n
    `Norm`  - Normalisation of DistEn value:  \n
              [false]  no normalisation.
              [true]   normalises w.r.t # of histogram bins - default

    # See also `XDistEn`, `DistEn2D`, `MSEn`, `K2En`

    # References:
        [1] Li, Peng, et al.,
            "Assessing the complexity of short-term heartbeat interval
            series by distribution entropy."
            Medical & biological engineering & computing
            53.1 (2015): 77-87.

    """
    function DistEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1,
        Bins::Union{Int,String}="Sturges", Logx::Real=2, Norm::Bool=true)

    Logx == 0  ? Logx = exp(1) : nothing

    N = size(Sig)[1]
    (N>10) ? nothing : error("Sig:  must be a numeric vector")
    (m > 0) ? nothing : error("m:   must be an integer > 0")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing :  error("Logx:  must be a positive number > 0")
    if typeof(Bins)<:Int
        (Bins>1) ? nothing :
        error("Bins:  must be an integer > 1 (or name of binning method)")
    elseif typeof(Bins)<:String
        (lowercase(Bins) in ["sturges","sqrt","rice","doanes"]) ? nothing :
        error("Bins:    must be one of the following strings
                'sturges', 'sqrt', 'rice', 'doanes' (or an integer >1)")
    end

    Nx = size(Sig)[1] - ((m-1)*tau)
    Zm = zeros(Nx,m)
    for n = 1:m
        Zm[:,n] = Sig[(n-1)*tau + 1:Nx+(n-1)*tau]
    end

    DistMat = zeros(Int(Nx*(Nx-1)/2))
    for k = 1:Nx-1
        Ix = [Int((k-1)*(Nx - k/2)+1), Int(k*(Nx-((k+1)/2)))]
        DistMat[Ix[1]:Ix[2]] = maximum(abs.(transpose(Zm[k,:]) .- Zm[k+1:end,:]),dims=2)
        # DistMat[Ix[1]:Ix[2]] = maximum(abs.(repeat(Zm[k,:],outer=(1,Nx-k)) .- Zm[k+1:end,:]'),dims=2)
    end

    Ny = size(DistMat)[1]
    if eltype(Bins)<:Char
        if lowercase(Bins) == "sturges"
            Bx = ceil(log2(Ny) + 1)
        elseif lowercase(Bins) == "rice"
            Bx = ceil(2*(Ny^(1/3)))
        elseif lowercase(Bins) == "sqrt"
            Bx = ceil(sqrt(Ny))
        elseif lowercase(Bins) == "doanes"
            sigma = sqrt(6*(Ny-2)/((Ny+1)*(Ny+3)))
            Bx = ceil(1+log2(Ny)+log2(1+abs(skewness(DistMat)/sigma)))
        else
            error("Please enter a valid binning method")
        end
    else
        Bx = Bins
    end
    By = collect(range(minimum(DistMat),maximum(DistMat),length=Int(Bx+1)))
    By[end] += 1; By[1]-= 1
    Ppi = fit(Histogram, DistMat, By).weights/Ny

    if round(sum(Ppi),digits=6) != 1
        @warn("Potential error estimating probabilities (p = $(sum(Ppi))")
        Ppi = Ppi[Ppi.!=0]
    elseif any(Ppi.==0)
        print("Note: $(sum(Ppi.==0))/$(length(Ppi)) bins were empty \n")
        Ppi = Ppi[Ppi.!=0]
    end
    Dist = -sum(Ppi.*log.(Logx, Ppi))
    if Norm
        Dist = Dist/(log(Logx, Bx));
    end

    return Dist, Ppi
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

        http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
matplotlib.use("agg") # or try "macosx"
using SparseArrays

include("chip_unstructured_geometry.jl")
include("chip_unstructured_solver.jl")

# physical constants
k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

k_chip = k_chip_ref * ones(length(chip_fem_idx))

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp; 0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

# unit quantities for nondimensionalization
u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 14    # number of iterations for Newton's method

# k_fem = k_air * constant(ones(ndof))
# k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
# k_fem = scatter_update(k_fem, chip_fem_idx, k_chip * ones(length(chip_fem_idx)))
# kgauss = dof_to_gauss_points(k_fem, mesh)
# LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, mesh))

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

B = constant(compute_interaction_matrix(mesh))
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
     bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

S0 = zeros(nelem+3*ndof)
S = solve_navier_stokes(S0, NT, k_chip)

sess = Session(); init(sess)
output = run(sess, S)

matwrite("xz_chip_unstructured_data.mat",
    Dict(
        "V"=>output[end, :]
    ))

u_out, v_out, p_out, T_out = output[NT+1,1:nnode], output[NT+1,ndof+1:ndof+nnode],
                             output[NT+1,2*ndof+1:2*ndof+nelem],output[NT+1,2*ndof+nelem+1:2*ndof+nelem+nnode]

figure(figsize=(10,10))
subplot(221)
title("u velocity")
visualize_scalar_on_fem_points(u_out .* u_std, mesh)
subplot(222)
title("v velocity")
visualize_scalar_on_fem_points(v_out .* u_std, mesh)
subplot(223)
visualize_scalar_on_fvm_points(p_out .* p_std, mesh)
title("pressure")
subplot(224)
title("temperature")
visualize_scalar_on_fem_points(T_out.* T_infty .+ T_infty, mesh)
tight_layout()
savefig("forward_solution_unstructured.pdf")

print("Solution range:",
    "\n [u velocity] \t min:", minimum(u_out .* u_std), ",\t max:", maximum(u_out .* u_std),
    "\n [v velocity] \t min:", minimum(v_out .* u_std), ",\t max:", maximum(v_out .* u_std),
    "\n [pressure]   \t min:", minimum(p_out .* p_std), ",\t max:", maximum(p_out .* p_std),
    "\n [temperature]\t min:", minimum(T_out.* T_infty .+ T_infty), ",\t\t\t max:", maximum(T_out.* T_infty .+ T_infty))


text:
module LASindex

    using FileIO
    using RegionTrees
    using StaticArrays

    export
    quadtree

    include("fileio.jl")
    include("quadtree.jl")
    include("spatial.jl")
    include("util.jl")

    function __init__()
      add_format(format"LAX", "LASX", ".lax", [:LASindex])
    end

end # module


text:
export deformulaMinusOneToOne, deformulaZeroToInf, deint

struct Formula{T <: Real}
    range::Tuple{T,T}
    phi
    phidash
end

const deformulaZeroToInf = Formula(
    (-6.8, 6.8),
    t -> exp(pi * sinh(t) / 2),
    t -> pi * cosh(t) * exp(pi * sinh(t) / 2) / 2
)

const deformulaMinusOneToOne = Formula(
    (-3.0, 3.0),
    t -> tanh(pi * sinh(t) / 2),
    t -> pi * cosh(t) * (1 / cosh(pi * sinh(t) / 2)) * (1 / cosh(pi * sinh(t) / 2)) / 2
)

struct DeformulaResult{T <: Real}
    s::T
    t::Vector{T}
    x::Vector{T}
    w::Vector{T}
    h::T
end

function _calcWeight!(data::Vector{Tuple{T,T,T}}, t::T, f, phi, phidash; abstol::T = eps(T)) where {T <: Real}
    local xtmp::T = phi(t)
    local wtmp::T = phidash(t) * f(xtmp)
    if !isnan(wtmp) && abs(wtmp) > abstol
        if !isfinite(wtmp)
            error("Error: weight becomes NaN")
        end
        push!(data, (t, xtmp, wtmp))
    end
end

"""
_deint(f, formula; reltol::T = 1.0e-8, abstol::T = eps(T), d = 8, maxiter = 16)

Compute the numerical integration for f with double exponential formula.

    int f(x) dx = h * sum w_i

Parameters:
- f: integrand function
- formula: double expoential formula.
- reltol: tolerance for relative errors
- abstol: tolerance for absolute errors
- d: the initial number of divides
- maxiter: the maximum number of iterations to increase the number of divides twice.
Return value (tuple):
- s: the value of integration
- t: a sequence for divides on the transformed domain
- x: a sequence for divides; x_i
- w: a sequence of weights (unscaled); w_i
- h: a scale parameter for w_i
"""

function _deint(f, formula::Formula{T};
        reltol::T = 1.0e-9, abstol::T = eps(T), d = 8, maxiter = 16) where {T <: Real}
    local lower::T = formula.range[1]
    local upper::T = formula.range[2]
    local h::T = (upper - lower) / d

    data = Vector{Tuple{T,T,T}}()
    for t = LinRange(lower, upper, d+1)
        _calcWeight!(data, t, f, formula.phi, formula.phidash, abstol=abstol)
    end
    local s::T = sum([x[3] for x in data]) * h

    local iter = 1
    while true
        iter += 1
        prev = s + 1.0
        if iter > maxiter
            error("Warning: The number of iteration attains maxtier $(maxiter)")
        end
        h /= 2
        for t = LinRange(lower+h, upper-h, d)
            _calcWeight!(data, t, f, formula.phi, formula.phidash, abstol=abstol)
        end
        d *= 2
        s = sum([x[3] for x in data]) * h
        aerror = s + 1.0 - prev
        rerror = aerror / prev
        if abs(rerror) < reltol
            break
        end
    end
    sort!(data, by=first)
    DeformulaResult(s, map(x->x[1], data), map(x->x[2], data), map(x->x[3], data), h)
end

"""
deint(f, lower, upper; reltol::T = 1.0e-8, abstol::T = eps(T), d = 8, maxiter = 16)

Compute the numerical integration for f on [lower, upper] with double exponential formula.

    int f(x) dx = h * sum w_i

Parameters:
- f: integrand function
- lower: lower value of range. This should be finite.
- upper: upper value of range. This is allowed to take infinite.
- reltol: tolerance for relative errors
- abstol: tolerance for absolute errors
- d: the initial number of divides
- maxiter: the maximum number of iterations to increase the number of divides twice.
Return value (tuple):
- the value of integration
- a sequence for divides on the transformed domain
- a sequence for divides; t_i
- a sequence of weights (unscaled); w_i
- a scale parameter for w_i; h
"""

function deint(f, lower::Float64, upper::Float64;
    reltol::Float64 = 1.0e-8, abstol::Float64 = eps(Float64), d = 8, maxiter = 16)
    if lower == 0.0 && upper == Inf64
        _deint(f, deformulaZeroToInf, reltol=reltol, abstol=abstol, d=d, maxiter=maxiter)
    elseif lower != 0.0
        result = _deint(deformulaZeroToInf, reltol=reltol, abstol=abstol, d=d, maxiter=maxiter) do x
            f(x + lower)
        end
        x = result.x .+ lower
        DeformulaResult(result.s, result.t, x, result.w, result.h)
    elseif lower == -1.0 && upper == 1.0
        _deint(f, deformulaMinusOneToOne, reltol=reltol, abstol=abstol, d=d, maxiter=maxiter)
    elseif upper != Inf64
        result = _deint(deformulaMinusOneToOne, reltol=reltol, abstol=abstol, d=d, maxiter=maxiter) do x
            d = (upper - lower)/2
            f(d*(x + 1) + lower) * d
        end
        x = ((upper - lower)/2) .* (result.x .+ 1.0) .+ lower
        DeformulaResult(result.s, result.t, x, result.w, result.h)
    else
        error("lower and upper should be finite or semi-infinite.")
    end
end


text:
using Test, YaoBase

@test_throws ErrorException @assert_locs_inbounds 4 (1, 2, 3, 4, 5)
@test_throws ErrorException @assert_locs_inbounds 4 collect(1:5)

@test_throws LocationConflictError @assert_locs_safe 4 (1, 2, 2, 1)
@test_throws LocationConflictError @assert_locs_safe 4 [1, 3, 2, 2]


text:
using Test, SpecialFunctions
import TruncatedNormal as TN

@test TN.tnmom1c(0, -Inf, Inf) == 0
@test TN.tnmom1c(0, -Inf, 0) ≈ TN.tnmean(-Inf, 0)
@test TN.tnmom1c(0, 0, +Inf) ≈ TN.tnmean(0, +Inf)

@test TN.tnmom1c(0, +Inf, +Inf) == +Inf
@test TN.tnmom1c(0, -Inf, -Inf) == -Inf
@test TN.tnmom1c(1, +Inf, +Inf) == +Inf
@test TN.tnmom1c(1, -Inf, -Inf) == -Inf
@test TN.tnmom1c(-Inf, +Inf, +Inf) == +Inf
@test TN.tnmom1c(+Inf, -Inf, -Inf) == -Inf

@test isnan(TN.tnmom1c(+Inf, +Inf, +Inf))
@test isnan(TN.tnmom1c(-Inf, -Inf, -Inf))

@test isnan(TN.tnmom1c(0, 1, 0))
@test isnan(TN.tnmom1c(0, +Inf, 0))
@test isnan(TN.tnmom1c(0, 0, -Inf))
@test isnan(TN.tnmom1c(0, +Inf, -Inf))
@test isnan(TN.tnmom1c(1, 1, 0))
@test isnan(TN.tnmom1c(1, +Inf, 0))
@test isnan(TN.tnmom1c(1, 0, -Inf))

for c = (0, 1), x = (0, NaN), y = (0, NaN)
    if isnan(x) || isnan(y)
        @test isnan(TN.tnmom1c(c, x, y))
    end
end

for x = exp.(-100:100), c = -100:100
    @test TN.tnmom1c(c, -x, x) ≈ -c
end

for x = -10:10, c = -10:10
    @test TN.tnmom1c(c, x, +Inf) ≈ TN.tnmean(x, +Inf) - c
    @test TN.tnmom1c(c, -Inf, x) ≈ TN.tnmean(-Inf, x) - c
end

#= TN.tnmom1c(c, a, b) = 0 when c == tnmean(a,b) =#
for a = exp.(-10:100), b = exp.(-10:100)
    a ≤ b || continue
    c = TN.tnmean(a, b)
    @test abs(TN.tnmom1c(c, a, b)) ≤ max(1e-10, c * 1e-10)
    @test TN.tnmom1c(a, a, b) ≥ 0
    @test TN.tnmom1c(b, a, b) ≤ 0
end

@test TN.tnmom1c(0, -Inf, Inf, 0, 1) ≈ 0

for x = -10:10
    @test TN.tnmom1c(0, x, +Inf, 0, 1) ≈ TN.tnmean(x, +Inf)
    @test TN.tnmom1c(0, -Inf, x, 0, 1) ≈ TN.tnmean(-Inf, x)
end

@test TN.tnmom1c(10, 0, 10) ≈ -9.20211543919713464412024961258
@test TN.tnmom1c(1, 3, 4) ≈ 2.26045428559002115390313889722
@test TN.tnmom1c(-1, 3, 4) ≈ 4.26045428559002115390313889722
@test TN.tnmom1c(2, 3, 4) ≈ 1.26045428559002115390313889722
@test TN.tnmom1c(-2, 3, 4) ≈ 5.26045428559002115390313889722
@test TN.tnmom1c(0, -1e200, 1e200) ≈ 0
@test TN.tnmom1c(1, -1e200, 1e200) ≈ -1
@test TN.tnmom1c(1e200, -1e200, 1e200) ≈ -1e200
@test TN.tnmom1c(0, -1e6, -999000) ≈ TN.tnmean(-1e6, -999000)
@test TN.tnmom1c(0, +1e6, +Inf) ≈ TN.tnmean(+1e6, +Inf)
@test TN.tnmom1c(0, -Inf, -1e6) ≈ TN.tnmean(-Inf, -1e6)

@test TN.tnmom1c(999000, 999000, 1e6) ≈ 1.00100100099899498898098100910e-6
@test TN.tnmom1c(999500, 999000, 1e6) ≈ -499.999998998998999001005011019
@test TN.tnmom1c(1e6, 999000, 1e6) ≈ -999.999998998998999001005011019


text:
@testset "Game" begin
    @test all(hasfield.(WordleGame, [:target, :number, :guesses]))

    @testset "Contructors" begin
        @testset "Valid Inputs" begin
            # target word must be in wordle list or in list of valid guesses
            @test_throws ErrorException WordleGame("aaaaa")
            @test_throws ErrorException WordleGame("AAAAA")

            # valid target word can be specified in capitalization agnostic form
            @test typeof(WordleGame("words")) == WordleGame
            @test typeof(WordleGame("Words")) == WordleGame
            @test typeof(WordleGame("WORDS")) == WordleGame
            @test typeof(WordleGame("wOrDs")) == WordleGame

            # target word must be five letters long
            @test typeof(WordleGame("words")) == WordleGame
            @test_throws ErrorException WordleGame("a")
            @test_throws ErrorException WordleGame("aaaaaaa")
            @test_throws ErrorException WordleGame("")

            # number can be specified
            @test typeof(WordleGame("words", 1, Wordle.WordleGuess[])) == WordleGame

            # 0 < number < LATEST_WORDLE_NUMBER
            @test_throws ErrorException WordleGame("words", 0, Wordle.WordleGuess[])
            @test_throws ErrorException WordleGame("words", -1, Wordle.WordleGuess[])
            @test begin
                typeof(WordleGame("words", Wordle.LATEST_WORDLE_NUMBER,
                                  Wordle.WordleGuess[])) == WordleGame
            end
            @test_throws ErrorException begin
                WordleGame("words", Wordle.LATEST_WORDLE_NUMBER + 1,
                           Wordle.WordleGuess[])
            end

            # number of guesses cannot be more than six
            for i = 1:7
                🤔 = Wordle.WordleGuess("words", fill(Wordle.CORRECT, 5))

                if i > 6
                    @test_throws ErrorException WordleGame("words", nothing, fill(🤔, i))
                else
                    @test @suppress typeof(WordleGame("words", nothing, fill(🤔, i))) == WordleGame
                end
            end
        end

        @testset "Convenience Constructors" begin
            @testset "Specify WordleGame by Wordle Number" begin
                @test typeof(WordleGame(1)) == WordleGame

                game = WordleGame(1)
                @test game.target == first(Wordle.WORDLE_LIST)
                @test game.number == 1
            end

            @testset "Daily Wordle" begin
                @test typeof(WordleGame()) == WordleGame

                game = WordleGame()
                also_game = WordleGame(Wordle.LATEST_WORDLE_NUMBER)

                @test game.target == also_game.target
                @test game.number == also_game.number
            end
        end

        @testset "Display" begin
            game = WordleGame("words")

            @test (@capture_out print(game)) == "Wordle 0/6"

            # make a really bad guess
            guess!(game, "quiet")

            @test begin
                (@capture_out print(game)) == strip("""
                Wordle 1/6

                ⬛⬛⬛⬛⬛
                """)
            end

            # make a guess a present letter
            guess!(game, "dodge")

            @test begin
                (@capture_out print(game)) == strip("""
                Wordle 2/6

                ⬛⬛⬛⬛⬛
                🟨🟩⬛⬛⬛
                """)
            end

            # make a guess with two of the same, present letters. Should only
            # get feedback for one of the two.
            guess!(game, "looks")

            @test begin
                (@capture_out print(game)) == strip("""
                Wordle 3/6

                ⬛⬛⬛⬛⬛
                🟨🟩⬛⬛⬛
                ⬛🟩⬛⬛🟩
                """)
            end

            # make a guess that changes a known letter to a correct letter
            guess!(game, "lords")

            @test begin
                (@capture_out print(game)) == strip("""
                Wordle 4/6

                ⬛⬛⬛⬛⬛
                🟨🟩⬛⬛⬛
                ⬛🟩⬛⬛🟩
                ⬛🟩🟩🟩🟩
                """)
            end

            # make a guess the word!
            guess!(game, "words")

            @test begin
                (@capture_out print(game)) == strip("""
                Wordle 5/6

                ⬛⬛⬛⬛⬛
                🟨🟩⬛⬛⬛
                ⬛🟩⬛⬛🟩
                ⬛🟩🟩🟩🟩
                🟩🟩🟩🟩🟩
                """)
            end
        end

        @testset "Guesses" begin
            game = WordleGame("words")

            for i = 1:6
                @test guess!(game, "bored") !== nothing
            end

            # can't guess more than 6 times
            @test_throws ErrorException guess!(game, "bored")

            # can't guess after target has been guessed
            game = WordleGame("words")
            @test @suppress guess!(game, "words") !== nothing
            @test_throws ErrorException guess!(game, "bored")

            game = WordleGame("seams")
            # if n is the number of times a letter appears in a wordle and a
            # guess has >n occurrences of that letter, only the first n occurrences
            # have feedback
            @test @capture_out(print(guess(game, "sissy"))) == "🟩⬛🟨⬛⬛"
            @test @capture_out(print(guess(game, "beams"))) == "⬛🟩🟩🟩🟩"
        end

        @testset "Hard Mode" begin
            game = WordleGame("words"; hard = true)
            guess!(game, "bored")

            # if in hard mode, your next guess must contain the letters you know
            # are in the wordle
            @test_throws ErrorException guess!(game, "cubic")
            @test typeof(guess!(game, "cords")) == WordleGame

            game = WordleGame("words")
            guess!(game, "bored")

            # but this sequences of guesses should be allowed in regular mode.
            @test typeof(guess!(game, "cubic")) == WordleGame
        end

        @testset "Available Letters" begin
            game = WordleGame("words")
            @test available_letters(game) == fill(Wordle.UNGUESSED, 26)
            @test (@capture_out print_available_letters(game)) == """
            a b c d e f g h i j k l m n o p q r s t u v w x y z
            ⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜⬜
            """

            # printing available letters uses available letters and happens to
            # be a bit more compact so we will use its outputs as a surrogate
            guess!(game, "quiet")
            @test (@capture_out print_available_letters(game)) == """
            a b c d e f g h i j k l m n o p q r s t u v w x y z
            ⬜⬜⬜⬜⬛⬜⬜⬜⬛⬜⬜⬜⬜⬜⬜⬜⬛⬜⬜⬛⬛⬜⬜⬜⬜⬜
            """

            guess!(game, "dodge")
            @test (@capture_out print_available_letters(game)) == """
            a b c d e f g h i j k l m n o p q r s t u v w x y z
            ⬜⬜⬜🟨⬛⬜⬛⬜⬛⬜⬜⬜⬜⬜🟩⬜⬛⬜⬜⬛⬛⬜⬜⬜⬜⬜
            """

            # if the outcome for a letter in a given guess is present but it
            # was previously correct, then it should still be marked as
            # correct in this list
            guess!(game, "quiet")
            @test (@capture_out print_available_letters(game)) == """
            a b c d e f g h i j k l m n o p q r s t u v w x y z
            ⬜⬜⬜🟨⬛⬜⬛⬜⬛⬜⬜⬜⬜⬜🟩⬜⬛⬜⬜⬛⬛⬜⬜⬜⬜⬜
            """

            # let's test going from "present" to "correct"
            guess!(game, "woods")
            @test (@capture_out print_available_letters(game)) == """
            a b c d e f g h i j k l m n o p q r s t u v w x y z
            ⬜⬜⬜🟩⬛⬜⬛⬜⬛⬜⬜⬜⬜⬜🟩⬜⬛⬜🟩⬛⬛⬜🟩⬜⬜⬜
            """
        end

        @testset "Utility Functions" begin
            game = WordleGame("words")

            @test Wordle.target(game) == "words"
            @test nguess(game) == 0

            for i = 1:6
                guess!(game, "bored")
                @test nguess(game) == i
                @test "bored" ∈ game
            end
        end

    end
end


text:
using Serialization

mutable struct TextureAtlas
    rectangle_packer::RectanglePacker
    mapping         ::Dict{Any, Int} # styled glyph to index in sprite_attributes
    index           ::Int
    data            ::Matrix{Float16}
    attributes      ::Vector{Vec4f0}
    scale           ::Vector{Vec2f0}
    extent          ::Vector{FontExtent{Float64}}
end

function TextureAtlas(initial_size = (2048, 2048))
    TextureAtlas(
        RectanglePacker(SimpleRectangle(0, 0, initial_size...)),
        Dict{Any, Int}(),
        1,
        zeros(Float16, initial_size...),
        Vec4f0[],
        Vec2f0[],
        FontExtent{Float64}[]
    )
end
assetpath(files...) = joinpath(@__DIR__, "..", "..", "assets", files...)

begin #basically a singleton for the textureatlas
    # random list of chars we cache
    # basically to make runtests fast, until we figure out a better way to cache
    # newly rendered chars.
    const _tobe_cached = [
        'π','∮','⋅','→','∞','∑','∏','∀','∈','ℝ','⌈','⌉','−','⌊','⌋','α','∧','β','∨','ℕ','⊆','₀',
        '⊂','ℤ','ℚ','ℂ','⊥','≠','≡','≤','≪','⊤','⇒','⇔','₂','⇌','Ω','⌀',
    ]
    const _cache_path = joinpath(dirname(@__FILE__), "..", ".cache", "texture_atlas.jls")
    const _default_font = Vector{Ptr{FreeType.FT_FaceRec}}[]
    const _alternative_fonts = Vector{Ptr{FreeType.FT_FaceRec}}[]

    function defaultfont()
        if isempty(_default_font)
            push!(_default_font, newface(assetpath("fonts", "DejaVuSans.ttf")))
        end
        _default_font[]
    end

    function alternativefonts()
        if isempty(_alternative_fonts)
            alternatives = [
                "DejaVuSans.ttf",
                "NotoSansCJKkr-Regular.otf",
                "NotoSansCuneiform-Regular.ttf",
                "NotoSansSymbols-Regular.ttf",
                "FiraMono-Medium.ttf"
            ]
            for font in alternatives
                push!(_alternative_fonts, newface(assetpath("fonts", font)))
            end
        end
        _alternative_fonts
    end

    function cached_load()
        if isfile(_cache_path)
            try
                return open(_cache_path) do io
                    dict = Serialization.deserialize(io)
                    fields = map(fieldnames(TextureAtlas)) do n
                        v = dict[n]
                        isa(v, Vector) ? copy(v) : v # otherwise there seems to be a problem with resizing
                    end
                    TextureAtlas(fields...)
                end
            catch e
                @info("You can likely ignore the following warning, if you just switched Julia versions for GLVisualize")
                @warn(e)
                rm(_cache_path)
            end
        end
        atlas = TextureAtlas()
        @info("Caching fonts, this may take a while. Needed only on first run!")
        for c in '\u0000':'\u00ff' #make sure all ascii is mapped linearly
            insert_glyph!(atlas, c, defaultfont())
        end
        for c in _tobe_cached
            insert_glyph!(atlas, c, defaultfont())
        end
        to_cache(atlas) # cache it
        return atlas
    end

    function to_cache(atlas)
        if !ispath(dirname(_cache_path))
            mkdir(dirname(_cache_path))
        end
        open(_cache_path, "w") do io
            dict = Dict(map(fieldnames(typeof(atlas))) do name
                name => getfield(atlas, name)
            end)
            Serialization.serialize(io, dict)
        end
    end
    const global_texture_atlas = RefValue{TextureAtlas}()
    function get_texture_atlas()
        if isassigned(global_texture_atlas)
            global_texture_atlas[]
        else
            global_texture_atlas[] = cached_load() # initialize only on demand
            global_texture_atlas[]
        end
    end

end

function glyph_index!(atlas::TextureAtlas, c::Char, font)
    if FT_Get_Char_Index(font[], c) == 0
        for afont in alternativefonts()
            if FT_Get_Char_Index(afont[], c) != 0
                font = afont
            end
        end
    end
    if c < '\u00ff' && font == defaultfont() # characters up to '\u00ff'(255), are directly mapped for default font
        Int(c)+1
    else #others must be looked up, since they're inserted when used first
        return insert_glyph!(atlas, c, font)
    end
end

glyph_scale!(c::Char, scale) = glyph_scale!(get_texture_atlas(), c, defaultfont(), scale)
glyph_uv_width!(c::Char) = glyph_uv_width!(get_texture_atlas(), c, defaultfont())

function glyph_uv_width!(atlas::TextureAtlas, c::Char, font)
    atlas.attributes[glyph_index!(atlas, c, font)]
end
function glyph_scale!(atlas::TextureAtlas, c::Char, font, scale)
    atlas.scale[glyph_index!(atlas, c, font)] .* (scale * 0.02)
end
function glyph_extent!(atlas::TextureAtlas, c::Char, font)
    atlas.extent[glyph_index!(atlas, c, font)]
end

function bearing(extent)
    Point2f0(
        extent.horizontal_bearing[1],
        -(extent.scale[2] - extent.horizontal_bearing[2])
    )
end

function glyph_bearing!(atlas::TextureAtlas, c::Char, font, scale)
    bearing(atlas.extent[glyph_index!(atlas, c, font)]) .* Point2f0(scale * 0.02)
end

function glyph_advance!(atlas::TextureAtlas, c::Char, font, scale)
    atlas.extent[glyph_index!(atlas, c, font)].advance .* (scale * 0.02)
end

insert_glyph!(atlas::TextureAtlas, glyph::Char, font) = get!(atlas.mapping, (glyph, font)) do
    uv, extent, width_nopadd, pad = render(atlas, glyph, font)
    tex_size       = Vec2f0(size(atlas.data))
    uv_start       = Vec2f0(uv.x, uv.y)
    uv_width       = Vec2f0(uv.w, uv.h)
    real_start     = uv_start + pad - 1 # include padding
    # padd one additional pixel
    relative_start = real_start ./ tex_size # use normalized texture coordinates
    relative_width = (real_start + width_nopadd + 2) ./ tex_size

    uv_offset_width = Vec4f0(relative_start..., relative_width...)
    i = atlas.index
    push!(atlas.attributes, uv_offset_width)
    push!(atlas.scale, Vec2f0(width_nopadd + 2))
    push!(atlas.extent, extent)
    atlas.index = i + 1
    return i
end

function sdistancefield(img, downsample = 8, pad = 8*downsample)
    w, h = size(img)
    wpad = 0; hpad = 0;
    while w % downsample != 0
        w += 1
    end
    while h % downsample != 0
        h += 1
    end
    w, h = w + 2pad, h + 2pad #pad this, to avoid cuttoffs

    in_or_out = Matrix{Bool}(undef, w, h)
    @inbounds for i=1:w, j=1:h
        x, y = i-pad, j-pad
        in_or_out[i,j] = checkbounds(Bool, img, x, y) && img[x,y] > 0.5 * 255
    end
    yres, xres = div(w, downsample), div(h, downsample)
    sd = sdf(in_or_out, xres, yres)
    Float16.(sd)
end

const font_render_callbacks = Function[]

function font_render_callback!(f)
    push!(font_render_callbacks, f)
end

function remove_font_render_callback!(f)
    filter!(f2-> f2 != f, font_render_callbacks)
end

function render(atlas::TextureAtlas, glyph::Char, font, downsample = 5, pad = 8)
    #select_font_face(cc, font)
    if glyph == '\n' # don't render  newline
        glyph = ' '
    end
    bitmap, extent = renderface(font, glyph, (50*downsample, 50*downsample))
    sd = sdistancefield(bitmap, downsample, downsample*pad)
    extent = extent ./ Vec2f0(downsample)
    rect = SimpleRectangle(0, 0, size(sd)...)
    uv = push!(atlas.rectangle_packer, rect) #find out where to place the rectangle
    uv == nothing && error("texture atlas is too small. Resizing not implemented yet. Please file an issue at GLVisualize if you encounter this") #TODO resize surface
    atlas.data[uv.area] = sd
    for f in font_render_callbacks
        f(sd, uv.area)
    end
    uv.area, extent, Vec2f0(size(bitmap)) ./ downsample, pad
end

make_iter(x) = repeated(x)
make_iter(x::AbstractArray) = x

function get_iter(defaultfunc, dictlike, key)
    make_iter(get(defaultfunc, dictlike, key))
end

using Base.Iterators: Repeated

function getposition(text, text2, fonts, scales, start_pos)
    calc_position(text2, start_pos, scales, fonts, text.text.atlas)
end
function getoffsets(text, text2, fonts, scales)
    calc_offset(text2, scales, fonts, text.text.atlas)
end

iter_or_array(x) = repeated(x)
iter_or_array(x::Repeated) = x
iter_or_array(x::Array) = x
iter_or_array(x::Vector{Ptr{FreeType.FT_FaceRec}}) = repeated(x)

function calc_position(
        last_pos, start_pos,
        atlas, glyph, font,
        scale, lineheight = 1.5
    )
    advance_x, advance_y = glyph_advance!(atlas, glyph, font, scale)
    if isnewline(glyph)
        return Point2f0(start_pos[1], last_pos[2] - advance_y * lineheight) #reset to startx
    else
        return last_pos + Point2f0(advance_x, 0)
    end
end
function calc_position(glyphs, start_pos, scales, fonts, atlas)
    positions = zeros(Point2f0, length(glyphs))
    last_pos  = Point2f0(start_pos)
    s, f = iter_or_array(scales), iter_or_array(fonts)
    c1 = first(glyphs)
    for (i, (c2, scale, font)) in enumerate(zip(glyphs, s, f))
        c2 == '\r' && continue # stupid windows!
        b = glyph_bearing!(atlas, c2, font, scale)
        positions[i] = last_pos .+ b
        last_pos = calc_position(last_pos, start_pos, atlas, c2, font, scale)
    end
    positions
end

function calc_offset(glyphs, scales, fonts, atlas)
    offsets = fill(Point2f0(0.0), length(glyphs))
    s, f = iter_or_array(scales), iter_or_array(fonts)
    c1 = first(glyphs)
    for (i, (c2, scale, font)) in enumerate(zip(glyphs, s, f))
        c2 == '\r' && continue # stupid windows!
        offsets[i] = Point2f0(glyph_bearing!(atlas, c2, font, scale))
        c1 = c2
    end
    offsets # bearing is the glyph offset
end

isnewline(x) = x == '\n'


text:
# Methods for porting Indicators.jl functions to TS objects from Temporal.jl package
function close_fun(X::TS, f::Function, flds::Vector{Symbol}; args...)
    if size(X,2) == 1
        return ts(f(X.values; args...), X.index, flds)
    elseif size(X,2) > 1 && has_close(X)
        return ts(f(cl(X).values; args...), X.index, flds)
    else
        error("Must be univariate or contain Close/Settle/Last.")
    end
end
function hlc_fun(X::TS, f::Function, flds::Vector{Symbol}; args...)
    if size(X,2) == 3
        return ts(f(X.values; args...), X.index, flds)
    elseif size(X,2) > 3 && has_high(X) && has_low(X) && has_close(X)
        return ts(f(hlc(X).values; args...), X.index, flds)
    else
        error("Argument must have 3 columns or have High, Low, and Close/Settle/Last fields.")
    end
end
function hl_fun(X::TS, f::Function, flds::Vector{Symbol}; args...)
    if size(X,2) == 2
        return ts(f(X.values; args...), X.index, flds)
    elseif size(X,2) > 2 && has_high(X) && has_low(X)
        return ts(f(hl(X).values; args...), X.index, flds)
    else
        error("Argument must have 2 columns or have High and Low fields.")
    end
end

###### run.jl ######
function runcov(x::TS{V,T}, y::TS{V,T}; args...) where {V,T}
    @assert size(x,2) == 1 && size(y,2) == 1 "Arguments x and y must both be univariate (have only one column)."
    z = [x y].values
    return ts(runcov(z[:,1], z[:,2]; args...), x.index, :RunCov)
end
function runcor(x::TS{V,T}, y::TS{V,T}; args...) where {V,T}
    @assert size(x,2) == 1 && size(y,2) == 1 "Arguments x and y must both be univariate (have only one column)."
    z = [x y].values
    ts(runcor(z[:,1], z[:,2]; args...), x.index, :RunCor)
end
mode(X::TS{V,T}) where {V,T} = mode(X.values)
runmean(X::TS{V,T}; args...) where {V,T} = close_fun(X, runmean, [:RunMean]; args...)
runsum(X::TS{V,T}; args...) where {V,T} = close_fun(X, runsum, [:RunSum]; args...)
runmad(X::TS{V,T}; args...) where {V,T} = close_fun(X, runmad, [:RunMAD]; args...)
runvar(X::TS{V,T}; args...) where {V,T} = close_fun(X, runvar, [:RunVar]; args...)
runmax(X::TS{V,T}; args...) where {V,T} = close_fun(X, runmax, [:RunMax]; args...)
runmin(X::TS{V,T}; args...) where {V,T} = close_fun(X, runmin, [:RunMin]; args...)
runsd(X::TS{V,T}; args...) where {V,T} = close_fun(X, runsd, [:RunSD]; args...)
runquantile(X::TS{V,T}; args...) where {V,T} = close_fun(X, runquantile, [:RunQuantile]; args...)
wilder_sum(X::TS{V,T}; args...) where {V,T} = close_fun(X, wilder_sum, [:WilderSum]; args...)
runacf(X::TS{V,T}; n::Int=10, maxlag::Int=n-3, lags::AbstractArray{Int,1}=0:maxlag, cumulative::Bool=true) where {V,T} =
    close_fun(X, runacf, [Symbol(i) for i in lags]; n=n, maxlag=maxlag, lags=lags, cumulative=cumulative)

##### ma.jl ######
sma(X::TS{V,T}; args...) where {V,T} = close_fun(X, sma, [:SMA]; args...)
hma(X::TS{V,T}; args...) where {V,T} = close_fun(X, hma, [:HMA]; args...)
mma(X::TS{V,T}; args...) where {V,T} = close_fun(X, mma, [:MMA]; args...)
swma(X::TS{V,T}; args...) where {V,T} = close_fun(X, swma, [:SWMA]; args...)
kama(X::TS{V,T}; args...) where {V,T} = close_fun(X, kama, [:KAMA]; args...)
alma(X::TS{V,T}; args...) where {V,T} = close_fun(X, alma, [:ALMA]; args...)
trima(X::TS{V,T}; args...) where {V,T} = close_fun(X, trima, [:TRIMA]; args...)
wma(X::TS{V,T}; args...) where {V,T} = close_fun(X, wma, [:WMA]; args...)
ema(X::TS{V,T}; args...) where {V,T} = close_fun(X, ema, [:EMA]; args...)
dema(X::TS{V,T}; args...) where {V,T} = close_fun(X, dema, [:DEMA]; args...)
tema(X::TS{V,T}; args...) where {V,T} = close_fun(X, tema, [:TEMA]; args...)
zlema(X::TS{V,T}; args...) where {V,T} = close_fun(X, zlema, [:ZLEMA]; args...)
mama(X::TS{V,T}; args...) where {V,T} = close_fun(X, mama, [:MAMA,:FAMA]; args...)

##### reg.jl ######
mlr_beta(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_beta, [:Intercept,:Slope]; args...)
mlr_slope(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_slope, [:Slope]; args...)
mlr_intercept(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_intercept, [:Intercept]; args...)
mlr(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr, [:MLR]; args...)
mlr_se(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_se, [:StdErr]; args...)
mlr_ub(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_ub, [:MLRUB]; args...)
mlr_lb(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_lb, [:MLRLB]; args...)
mlr_bands(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_bands, [:MLRLB,:MLR,:MLRUB]; args...)
mlr_rsq(X::TS{V,T}; args...) where {V,T} = close_fun(X, mlr_rsq, [:RSquared]; args...)

##### mom.jl ######
momentum(X::TS{V,T}; args...) where {V,T} = close_fun(X, momentum, [:Momentum]; args...)
roc(X::TS{V,T}; args...) where {V,T} = close_fun(X, roc, [:ROC]; args...)
macd(X::TS{V,T}; args...) where {V,T} = close_fun(X, macd, [:MACD,:Signal,:Histogram]; args...)
rsi(X::TS{V,T}; args...) where {V,T} = close_fun(X, rsi, [:RSI]; args...)
psar(X::TS{V,T}; args...) where {V,T} = hl_fun(X, psar, [:PSAR]; args...)
kst(X::TS{V,T}; args...) where {V,T} = close_fun(X, kst, [:KST]; args...)
wpr(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, wpr, [:WPR]; args...)
adx(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, adx, [:DiPlus,:DiMinus,:ADX]; args...)
cci(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, cci, [:CCI]; args...)
stoch(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, stoch, [:Stochastic,:Signal]; args...)
smi(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, smi, [:SMI,:Signal]; args...)
donch(X::TS{V,T}; args...) where {V,T} = hl_fun(X, donch, [:Low,:Mid,:High]; args...)
aroon(X::TS{V,T}; args...) where {V,T} = hl_fun(X, aroon, [:AroonUp,:AroonDn,:AroonOsc]; args...)

##### vol.jl ######
bbands(X::TS{V,T}; args...) where {V,T} = close_fun(X, bbands, [:LB,:MA,:UB]; args...)
tr(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, tr, [:TR]; args...)
atr(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, atr, [:ATR]; args...)
keltner(X::TS{V,T}; args...) where {V,T} = hlc_fun(X, keltner, [:KeltnerLower,:KeltnerMiddle,:KeltnerUpper]; args...)

##### trendy.jl #####
maxima(X::TS{V,T}; args...) where {V,T} = close_fun(X, maxima, [:Maxima]; args...)
minima(X::TS{V,T}; args...) where {V,T} = close_fun(X, minima, [:Minima]; args...)
support(X::TS{V,T}; args...) where {V,T} = close_fun(X, support, [:Support]; args...)
resistance(X::TS{V,T}; args...) where {V,T} = close_fun(X, resistance, [:Resistance]; args...)

#### utils.jl ####
crossover(x::TS, y::TS) = ts(crossover(x.values, y.values), x.index, [:CrossOver])
crossunder(x::TS, y::TS) = ts(crossunder(x.values, y.values), x.index, [:CrossUnder])


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.02"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	mu, sigma = 10, 4
	n, N = 10, 10^6

	sMeans = Array{Float64}(undef, N)
	sVars  = Array{Float64}(undef, N)
	tStats = Array{Float64}(undef, N)

	for i in 1:N
		data       = rand(Normal(mu,sigma),n)
		sampleMean = mean(data)
		sampleVars = var(data)
		sMeans[i]  = sampleMean
		sVars[i]   = sampleVars
		tStats[i]  = (sampleMean - mu)/(sqrt(sampleVars/n))
	end

	xRangeMean = 5:0.1:15
	xRangeVar = 0:0.1:60
	xRangeTStat = -5:0.1:5
end

# ╔═╡ 94ec1808-2529-11eb-35ba-e9565c6cdbc1
begin
	p1 = stephist(sMeans, bins=50, c=:blue, normed=true, legend=false)
	p1 = plot!(xRangeMean, pdf.(Normal(mu,sigma/sqrt(n)), xRangeMean),
		c=:red, xlims=(5,15), ylims=(0,0.35), xlabel="Sample mean",ylabel="Density")
	p2 = stephist(sVars, bins=50, c=:blue, normed=true, label="Simulated")
	p2 = plot!(xRangeVar, (n-1)/sigma^2*pdf.(Chisq(n-1), xRangeVar*(n-1)/sigma^2),
		c=:red, label="Analytic", xlims=(0,60), ylims=(0,0.06),
		xlabel="Sample Variance",ylabel="Density")
	p3 = stephist(tStats, bins=100, c=:blue, normed=true, legend=false)
	p3 = plot!(xRangeTStat, pdf.(TDist(n-1), xRangeTStat),
		c=:red, xlims=(-5,5), ylims=(0,0.4), xlabel="t-statistic",ylabel="Density")

	plot(p1, p2, p3, layout = (1,3), size=(1200, 600))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.02"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═94ec1808-2529-11eb-35ba-e9565c6cdbc1
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
"""
    tquotas(x::Union{Accumulated, Incremental}

Returns a vector containing the quotas development
for the incremental or accumulated triangle `x`.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> tquotas(ti)
3-element Vector{Float64}:
 0.4
 0.8
 1.0

 julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> tquotas(ta)
3-element Vector{Float64}:
 0.4
 0.8
 1.0
```
"""
function tquotas(x::Union{Accumulated, Incremental})
    return reverse(1 ./ cumprod(reverse(tfactors(x))))
end

"""
    expositionbf(primes::Vector, lossratios::Vector)

Returns a vector containing the estimated losses when
having a vector of `primes` and a vector of `lossratios`.

# Example
```
julia> primes = [15, 7, 3];

julia> lossratios = [0.85, 0.7, 0.79];

julia> expositionbf(primes, lossratios)
3-element Vector{Float64}:
 12.75
  4.8999999999999995
  2.37
```
"""
function expositionbf(primes::Vector, lossratios::Vector)
    return primes .* lossratios
end

"""
    fillbf(x::Accumulated, primes::Vector, lossratios::Vector)

Returns an accumulated completed triangle, containing both the known accumulated claims
and the Bornhuetter-Ferguson accumulated estimates, when having a `primes` vector and
a `lossratios` vector.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> primes = [15, 7, 3];

julia> lossratios = [0.85, 0.7, 0.79];

julia> fillbf(ta, primes, lossratios)
Accumulated Completed Triangle {Float64}:
 3.0  8.0    10.0
 4.0  6.0     6.98
 1.0  1.948   2.422
```
"""
function fillbf(x::Accumulated, primes::Vector, lossratios::Vector)
    cuo = tquotas(x)
    exposicion = expositionbf(primes, lossratios)
    matriz_completa = convert(Matrix{Float64}, deepcopy(x.claims))
    for i in 1:size(matriz_completa, 2)
        for j in 1:size(matriz_completa, 1)
            if (i + j) > (size(matriz_completa, 1) + 1)
                matriz_completa[i, j] = matriz_completa[i,size(matriz_completa, 2)+1-i] + (cuo[j] - cuo[size(matriz_completa, 2)+1-i]) * exposicion[i]
            end
        end
    end
    return AccumulatedCompleted(matriz_completa)
end

"""
    fillbf(x::Incremental, primes::Vector, lossratios::Vector)

Returns an incremental completed triangle, containing both the known incremental claims
and the Bornhuetter-Ferguson incremental estimates, when having a `primes` vector and
a `lossratios` vector.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> primes = [15, 7, 3];

julia> lossratios = [0.85, 0.7, 0.79];

julia> fillbf(ti, primes, lossratios)
Incremental Completed Triangle {Float64}:
 3.0  5.0    2.0
 4.0  2.0    0.98
 1.0  0.948  0.474
```
"""
function fillbf(x::Incremental, primes::Vector, lossratios::Vector)
    y = taccum(x)
    z = fillbf(y, primes, lossratios)
    des = zeros(Float64, size(z.claims, 1), size(z.claims, 1))
    des[:,1] = z.claims[:,1]
    for j in 2:size(z.claims, 1)
        for i in 1:size(z.claims, 1)
            des[i, j] = z.claims[i, j] - z.claims[i, j - 1]
        end
    end
    return IncrementalCompleted(des)
end

"""
    ibnrbf(x::Accumulated, primes::Vector, lossratios::Vector)

Returns the estimated incurred but not reported reserve under the Bornhuetter-Ferguson
methodology, for an accumulated triangle `x`, a `primes` vector and a `lossratios` vector.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> primes = [15, 7, 3];

julia> lossratios = [0.85, 0.7, 0.79];

julia> ibnrbf(ta, primes, lossratios)
2.4019999999999992
```
"""
function ibnrbf(x::Accumulated, primes::Vector, lossratios::Vector)
    tri_completado = fillbf(x, primes, lossratios)
    ib = 0
    for i in 1:size(tri_completado.claims, 2)
        j = (tri_completado.claims[i,end] - tri_completado.claims[i, end - i + 1])
        ib = ib + j
    end
    return ib
end

"""
    ibnrbf(x::Incremental, primes::Vector, lossratios::Vector)

Returns the estimated incurred but not reported reserve under the Bornhuetter-Ferguson
methodology, for an incremental triangle `x`, a `primes` vector and a `lossratios` vector.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> primes = [15, 7, 3];

julia> lossratios = [0.85, 0.7, 0.79];

julia> ibnrbf(ti, primes, lossratios)
2.4019999999999992
```
"""
function ibnrbf(x::Incremental, primes::Vector, lossratios::Vector)
    return ibnrbf(taccum(x), primes, lossratios)
end

"""
    lossratiosbf(x::Accumulated, primes::Vector, periods::Int)

Returns a vector of loss ratios, obtained from an accumulated triangle `x` and
a `primes` vector.
The loss ratios are obtained by summing over the last `periods` number of
incurred claims and primes.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> primes = [15, 7, 3];

julia> periods = 2;

julia> lossratiosbf(ta, primes, periods)
3-element Vector{Float64}:
 0.7272727272727273
 0.7272727272727273
 0.7
```
"""
function lossratiosbf(x::Accumulated, primes::Vector, periods::Int)
    xi = incurred(x)
    tmp = suma_movil_vector(xi, periods)./suma_movil_vector(primes, periods)
    lrbf= zeros(Float64, length(xi) - length(tmp))
    fill!(lrbf, tmp[1])
    return vcat(lrbf, tmp)
end


text:
@testset "SymbolicRegression" begin
    Random.seed!(1223)
    # Generate a multivariate function for OccamNet
    X = rand(2,20)
    f(x) = [sin(x[1]); exp(x[2])]
    Y = hcat(map(f, eachcol(X))...)
    # Define the options
    opts = EQSearch([+, *, sin, exp], maxdepth = 1)
    # Define the problem
    prob = DirectDataDrivenProblem(X, Y)
    # Solve the problem
    res = solve(prob, opts, numprocs = 0, multithreading = false)
    sys = result(res)

    x = states(sys)
    m = metrics(res)

    @test m.Complexity == 4
    @test m.Error <= eps()
    @test m.AICC == Inf
    @test isequal([x.rhs for x in equations(sys)], [sin(x[1]); exp(x[2])])
end


text:
"""
Fitting functions for DPP and k-DPP.
"""

# TODO


text:
module TestMatrixEqn

using QuantEcon
using Base.Test
using FactCheck

rough_kwargs = Dict(:atol => 1e-7, :rtol => 1e-7)

facts("Testing matrix_eqn.jl") do
    context("simple test where X is all zero") do
        A = eye(4) .* .95
        B = zeros(4, 4)

        X = solve_discrete_lyapunov(A, B)

        @fact X --> roughly(B)
    end

    context("simple test where X is same as B") do
        A = fill(0.5, 2, 2)
        B = [.5 -.5; -.5 .5]

        X = solve_discrete_lyapunov(A, B)

        @fact X --> roughly(B)
    end

    context("simple_ones lyap") do
        A = zeros(4, 4)
        B = ones(4, 4)

        sol = solve_discrete_lyapunov(A, B)

        @fact sol --> roughly(B)
    end

    context("scalars lyap") do
        a, b = 0.5, 0.75
        sol = solve_discrete_lyapunov(a, b)

        @fact sol --> ones(1, 1)
    end

    context("testing ricatti golden_num_float") do
        val = solve_discrete_riccati(1.0, 1.0, 1.0, 1.0)
        gold_ratio = (1 + sqrt(5)) / 2.
        @fact val[1] --> roughly(gold_ratio)
    end

    context("testing ricatti golden_num_2d") do
        A, B, R, Q = eye(2), eye(2), eye(2), eye(2)
        gold_diag = eye(2) .* (1 + sqrt(5)) ./ 2.
        val = solve_discrete_riccati(A, B, Q, R)
        @fact val --> roughly(gold_diag)
    end

    context("test tjm 1") do
        A = [0.0 0.1 0.0
             0.0 0.0 0.1
             0.0 0.0 0.0]
        B = [1.0 0.0
             0.0 0.0
             0.0 1.0]
        Q = [10^5 0.0 0.0
             0.0 10^3 0.0
             0.0 0.0 -10.0]
        R = [0.0 0.0
             0.0 1.0]
        X = solve_discrete_riccati(A, B, Q, R)
        Y = diagm([1e5, 1e3, 0.0])

        @fact X --> roughly(Y; rough_kwargs...)
    end

    context("test tjm 2") do
        A = [0.0 -1.0
             0.0 2.0]
        B = [1.0 0.0
             1.0 1.0]
        Q = [1.0 0.0
             0.0 0.0]
        R = [4.0 2.0
             2.0 1.0]
        X = solve_discrete_riccati(A, B, Q, R)
        Y = zeros(2, 2)
        Y[1, 1] = 1
        @fact X --> roughly(Y; rough_kwargs...)
    end

    context("test tjm 3") do
        r = 0.5
        I = eye(2)
        A = [2.0+r^2 0.0
             0.0     0.0]
        B = I
        R = [1.0 r
             r   r*r]
        Q = I - A' * A + A' * ((R + I) \ A)
        X = solve_discrete_riccati(A, B, Q, R)
        Y = eye(2)
        @fact X --> roughly(Y; rough_kwargs...)
    end

end  # facts
end  # module


text:
"""
_phunif(ph, ufact)

Get a uniformed PH.
"""

function _phunif(ph::GPH{Tv,MatT}, ufact::Tv = 1.01) where {Tv,MatT}
    P, qv = unif(ph.T, ufact)
    return ph.dim, ph.alpha, P, ph.tau, qv
end


text:
# rect_sino.jl

using MIRT: rect_sino, sino_geom, jim
import MIRT: trapezoid

using Plots: plot
using Test: @test, @test_throws, @inferred

"""
`rect_sino_show()`
show examples
"""
function rect_sino_show( ;
	down::Int = 4,
	rect::Matrix = [70 20 70 30 0 1],
	orbit::Real = 360,
	na::Int = 400,
	oversample::Int = 2,
)

#=
	ig = image_geom(nx=512, ny=504, fov=500)
	ig = ig.down(down)
	xtrue = rect_im(ig, rect; oversample=4)
=#

	arg = (na=na, down=down, orbit=orbit, offset=0.25)
	geoms = (
		sino_geom(:par, nb = 444, d = 1 ; arg...),
		sino_geom(:fan, nb = 888, d = 1 ; arg..., dsd = 949, dod = 408),
		sino_geom(:fan, nb = 888, d = 1 ; arg...,
			dsd = 949, dod = 408, dfs = Inf), # flat fan
		sino_geom(:moj, nb = 666, d = 1 ; arg...),
	)

	ngeom = length(geoms)
	pl = Array{Any}(undef, ngeom)

	for ii=1:ngeom
		sg = geoms[ii]
		over = (sg isa SinoMoj) ? 1 : oversample
		sino = rect_sino(sg, rect; oversample=over)
		pl[ii] = jim(sg.r, sg.ad, sino, title="$(typeof(sg))")
	end
	plot(pl...)
end

param = [60 10.5 38 27 0 1]
sg = sino_geom(:ge1, down=8)
#@inferred # todo
rect_sino(sg, param; xscale=-1, yscale=-1) # test scale
@inferred trapezoid(2.5, 1, 2, 3, 4.)
t = LinRange(-2, 5, 101)
f = trapezoid.(t, 1, 2, 3, 4.)

rect_sino_show()


text:
using DifferentialEquations
using LinearAlgebra
using Plots
using ForwardDiff
using ProgressMeter

# Thermal Model

function generate1DStencil!(A, N, Δx)
    # Generate the centered difference matrix manually
    for i in 1:N, j in 1:N

    if abs(i-j)<=1
        (A[i,j]+=1)
    end
    if i==j
        (A[i,j]-=3)
    end

    end

    A = A/(Δx^2)
end

function double_partial(x)
    dx = zeros(size(x))
    for i in 2:length(dx)-1
      dx[i] = -1*x[i-1] + 2*x[i] + -1*x[i+1]
    end
    dx[1] = 2*x[1] + -1*x[2]
    dx[end] = -1*x[end-1] + 2*x[end]
    dx
end

# Discretize in time and space
Δx = 1
x = Δx:Δx:100-Δx # Solve only for the interior: the endpoints are known to be zero!

Δt = 0.1
t = Δt:Δt:1000-Δt # Solve only for the interior: the endpoints are known to be zero!

N = length(x)
A = zeros(N,N)
generate1DStencil!(A, N, Δx)

# f(u, p, t) = zeros(size(u))
f(u, p, t) = gen_forcing(p)
# u0_func(x) = sin.(2π * x)

# Given x coord of components, give initial conditions that are equiv to a gaussian around that coordinate
function gen_forcing(positions)
    out_vec = zeros(size(x))
    out_vec = convert.(eltype(positions),out_vec)
    for pos in positions
        out_vec += gaussian.(x, pos, 5.0)
        # out_vec += box.(x, pos)
    end
    out_vec
end

function gaussian(x, mu, sig)
    return exp(-(x - mu)^2.0 / (2 * sig^2.0))
end

function box(x, pos)
    sz = 5.0
    if pos-sz <= x <= pos+sz
        return 1.0
    else
        return 0.0
    end
end

function heat_eq(u, p, t)
    du = A * u + f(u,p,t)
    du
end

# Predict
function predict(positions)
    # Simply roll this forward in time and we get our solution
    tspan = (0.0, 1.0)
    # tspan_dual = convert.(eltype(positions),tspan)
    # @show typeof(positions)
    # u₀ = gen_forcing(positions)
    u₀ = zeros(size(x))
    u_dual = convert.(eltype(positions),u₀)
    # @show typeof(u₀)
    prob = ODEProblem(heat_eq, u_dual, tspan, positions)
    sol = solve(prob)
    # @show typeof(sol[end])

    # display("Finished ODE Solve")
    sol
end

# Loss: Given a list of positions, what is the cost / loss of the
function loss(positions)
    simulation = predict(positions)
    # display(simulation)
    final_state = simulation[end]
    # display(size(final_state)[1])
    # display(final_state)
    return norm(final_state, 2)
end

# Gradient Descent: Optimize initial positions
function optimize(positions, η=0.1; plt=false)

    @showprogress for idx in 1:1000
        grads = ForwardDiff.gradient(s -> loss(s), positions) # Magic
        positions .-= η*grads
        # display(positions)
        # display(loss(positions))

        if loss(positions) < 0.1
            break
        end

        if plt
            sol = predict(positions)
            plt_sol(sol)
        end
    end

    positions
end

function plot_loss_surf()
    loss_surf = zeros((90, 90))
    for x_1 in 5.0:94.0
        for x_2 in 5.0:94.0
            loss_surf[Int(x_1) - 4, Int(x_2) - 4] = loss([x_1, x_2])
        end
    end

    # sol_f(x, y) = loss_surf[x, y]
    plot(loss_surf, st=:heatmap, color=:viridis)

end

# optimize([25.0, 30.0])

function plt_sol(sol)
    sol_f(t, x) = sol[t, x]
    display(plot(1:length(sol[end]),
            1:length(sol.t),
            sol_f,
            st=:surface,
            color=:viridis,
            xlabel="Position",
            ylabel="Time (sec)",
            title="1D Temperature Over Time"))
end

function plt_1D(sol)
    display(plot(sol[end]))
end

function plt_1D_w_pos(sol, positions)
    display(plot(sol[end],
                 linewidth=5,
                 xlabel="Position (meters)",
                 ylabel="Temperature",
                 label="Component positions at " * first(string(positions[1]),4) * " and " * first(string(positions[2]),4)))
end

function optimize_and_save_thermal(positions, η=0.1)
    solns = []
    append!(solns, [deepcopy(positions)])
    @showprogress for idx in 1:1000
        grads = ForwardDiff.gradient(s -> loss(s), positions) # Magic
        positions .-= η*grads
        # display(positions)
        # display(loss(positions))

        if loss(positions) < 0.1
            break
        end

        if idx in [500,1000]
            append!(solns, [deepcopy(positions)])
        end
    end

    solns
end

function test_and_save_thermal()
    p = [25.0, 30.0]
    solns = optimize_and_save_thermal(p)
    display(plot())
    i = 1
    idxs = [1,500,1000]
    for pos in solns
        sol = predict(pos)
        plt_1D_w_pos(sol, pos)
        plot!(title="1D Thermal Optimizer: Iteration " * string(idxs[i]), ylims = (0.0,2.0))
        savefig("plots/thermal_plot_"*string(i))
        i += 1
    end
end


text:
	rootbase = pwd()

	function webserverpath(folder::AbstractString)
	  global rootbase = joinpath(pwd(),folder)
	end

	"""
	File(file::String)

	Return a file content located in "rootbase"
	"""
	function File(file::String)
	    return File("", file)
	end

	"""
	File(folder::String, file::String)

	Return a file content located in "joinpath(rootbase, folder)"
	"""
	function File(folder::String, file::String)
	  path = joinpath(rootbase, folder, file)
	  f = open(path)
	  return read(f, String)
	end

	"""
	files(roop::String, file::String)

	Create a GET route to file
	"""
	function files(folder::String, file::String)
	    extension="text/plain"
	    ext= split(file,".")
	    if(length(ext)>1)
	      my_extension = ext[end]
	      if (haskey(mimetypes, my_extension))
	        extension = mimetypes[my_extension]
	      end
	    end
	    data = File(folder, file)
	    routefile = replace(joinpath(folder,file),"\\" => "/")
	    Get( string("/",routefile) , (request,HTTP)->(begin
	      HTTP.Response(200
          , HTTP.mkheaders(["Content-Type" => extension])
          , body=data)
	    end) )
	end

	"""
	Create routes to files inside "rootbase"
	"""
	function WebServer(folder::String, exten::String)
	  path = joinpath(rootbase, folder)
	  ls= readdir(path)
	  for i = 1:length(ls)
	    if isfile( joinpath(path, ls[i]) )
	      ext = split(ls[i], ".")
	      if length(ext) > 1 && ext[1] != "" && occursin(exten, ext[end] )
	        files(folder,ls[i])
	      end
	    elseif isdir( joinpath(path, ls[i]) )
	      WebServer(joinpath(folder, ls[i]) ,exten)
	    end
	  end
	end

	function webserverfiles(load::AbstractString)
	  if load == "*"
	    WebServer("", "")
	  else
	    WebServer("", load)
	  end
	end


text:
"Integration rule for approximating the continuous integrals for the equations of motion"
abstract type QuadratureRule end

"Integration rules of the form `x′ = f(x,u)`, where `x′` is the next state"
abstract type Explicit <: QuadratureRule end

"Integration rules of the form x′ = f(x,u,x′,u′), where x′,u′ are the states and controls at the next time step."
abstract type Implicit <: QuadratureRule end

(::Type{Q})(model::AbstractModel) where {Q<:QuadratureRule} = Q(state_dim(model), control_dim(model))

@autodiff struct DiscretizedDynamics{L,Q} <: DiscreteDynamics
    continuous_dynamics::L
    integrator::Q
    function DiscretizedDynamics(
        continuous_dynamics::L,
        integrator::Q,
    ) where {L<:ContinuousDynamics,Q<:QuadratureRule}
        new{L,Q}(continuous_dynamics, integrator)
    end
end
function DiscretizedDynamics{Q}(
    continuous_dynamics::L,
) where {L<:ContinuousDynamics} where {Q}
    n, m = size(continuous_dynamics)
    DiscretizedDynamics(continuous_dynamics, Q(n, m))
end

const ImplicitDynamicsModel{L,Q} = DiscretizedDynamics{L,Q} where {L,Q<:Implicit}

@inline state_dim(model::DiscretizedDynamics) = state_dim(model.continuous_dynamics)
@inline control_dim(model::DiscretizedDynamics) = control_dim(model.continuous_dynamics)
@inline errstate_dim(model::DiscretizedDynamics) = errstate_dim(model.continuous_dynamics)
@inline state_diff(model::DiscretizedDynamics, x, x0) = state_diff(model.continuous_dynamics, x, x0)
statevectortype(::Type{<:DiscretizedDynamics{L}}) where L = statevectortype(L)
default_diffmethod(::DiscretizedDynamics) = ForwardAD()

@inline integration(model::DiscretizedDynamics) = model.integrator
discrete_dynamics(model::DiscretizedDynamics, x, u, t, dt) =
    integrate(integration(model), model.continuous_dynamics, x, u, t, dt)
discrete_dynamics!(model::DiscretizedDynamics, xn, x, u, t, dt) =
    integrate!(integration(model), model.continuous_dynamics, xn, x, u, t, dt)

jacobian!(sig::FunctionSignature, ::UserDefined, model::DiscretizedDynamics, J, xn, z) =
    jacobian!(
        integration(model),
        sig,
        model.continuous_dynamics,
        J,
        xn,
        state(z),
        control(z),
        time(z),
        timestep(z),
    )

########################################
# Implicit Dynamics
########################################

# TODO: overwrite `evaluate` to solve for the next state using Newton
# TODO: overwrite `jacobian` to provide A,B using implicit function theorem

dynamics_error(
    model::ImplicitDynamicsModel,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error(integration(model), model.continuous_dynamics, z2, z1)

dynamics_error!(
    model::ImplicitDynamicsModel,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error!(integration(model), model.continuous_dynamics, y2, y1, z2, z1)

dynamics_error_jacobian!(
    sig::FunctionSignature,
    ::UserDefined,
    model::ImplicitDynamicsModel,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error_jacobian!(
    integration(model),
    sig,
    model.continuous_dynamics,
    J2,
    J1,
    y2,
    y1,
    z2,
    z1,
)


text:
include("../../src/Constructors.jl")
include("../../src/SelfPlay.jl")
include("../../src/ReplayBuffer.jl")
include("../../src/Learning.jl")
include("game.jl")
include("params.jl")

buffer = Dict{Int, GameHistory}()

latest_training_step=1000 #Change this to the iteration number of latest saved networks

# representation= deserialize(joinpath(conf.networks_path,"$(latest_training_step)_representation.bin"))
# prediction= deserialize(joinpath(conf.networks_path,"$(latest_training_step)_prediction.bin"))
# dynamics= deserialize(joinpath(conf.networks_path,"$(latest_training_step)_dynamics.bin"))

representation= init_representation(hyper)
prediction= init_prediction(hyper)
dynamics= init_dynamics(hyper)

NNs = (representation=representation, prediction=prediction, dynamics=dynamics)

env=TicTacToe()

competitive_play!(buffer_to_disk=false)


text:
struct ReshapeHandler end

is_hardware_primitive(::Type{typeof(reshape)}, args...) = true
gethandler(broadcasted, ::Type{typeof(reshape)}, args...) =
    !broadcasted ? ReshapeHandler() : error("Cannot generate hardware for broadcasted reshape.")
init_state(::ReshapeHandler) = nothing

function (handler::ReshapeHandler)(buffer, netlist, state, inputs, outputs)
    for input in inputs[2:end]
        delete!(netlist, input)
    end

    if issigned(inputs[1])
        setsigned!(netlist, outputs[1], true)
        write(buffer, """
            $stdcomment
            // BEGIN reshape
            assign $(name(outputs[1]))_p = $(name(inputs[1]))_p;
            assign $(name(outputs[1]))_m = $(name(inputs[1]))_m;
            // END reshape
            \n""")
    else
        write(buffer, """
            $stdcomment
            // BEGIN reshape
            assign $(name(outputs[1])) = $(name(inputs[1]));
            // END reshape
            \n""")
    end

    return buffer, state
end


text:
module CxxROS

using Cxx
include("std.jl")

try
    global const path_to_ros = normpath(joinpath(ENV["ROS_ROOT"], "../.."))
    isdir(path_to_ros)
catch
    throw(ErrorException("please set environment variavle: ROS_ROOT"))
end

addHeaderDir(path_to_ros * "include", kind=C_System)

cxx"""
#include <ros/ros.h>
"""
include("boost.jl")

# importing objects...

Libdl.dlopen(path_to_ros * "lib/libroscpp.so", Libdl.RTLD_GLOBAL)

include("time.jl")

# importing datatypes...

include("msgs/std_msgs.jl")
include("msgs/geometry_msgs.jl")
include("msgs/sensor_msgs.jl")
include("msgs/nav_msgs.jl")

include("tf.jl")

end #module

# no show interface for elements of boost::array cause exception
function Base.show(io::IO, ptr::Cxx.CppRef{Cxx.CxxArrayType{T},Cxx.NullCVR}) where T
    println(io, "Array<", Base.typename(T), ",1>")
end

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
module Turing

##############
# Dependency #
########################################################################
# NOTE: when using anything from external packages,                    #
#       let's keep the practice of explictly writing Package.something #
#       to indicate that's not implemented inside Turing.jl            #
########################################################################

using Requires, Reexport, ForwardDiff
using Bijectors, StatsFuns, SpecialFunctions
using Statistics, LinearAlgebra, ProgressMeter
using Markdown, Libtask, MacroTools
@reexport using Distributions, MCMCChains, Libtask
using Tracker: Tracker

import Base: ~, ==, convert, hash, promote_rule, rand, getindex, setindex!
import MCMCChains: AbstractChains, Chains

const PROGRESS = Ref(true)
function turnprogress(switch::Bool)
    @info("[Turing]: global PROGRESS is set as $switch")
    PROGRESS[] = switch
end

# Constants for caching
const CACHERESET  = 0b00
const CACHEIDCS   = 0b10
const CACHERANGES = 0b01

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_TURING", "0")))

# Include the interface. Temporary until the interface is moved
# to MCMCChains. CSP 2019-05-12
include("interface/Interface.jl")
using .Interface
import .Interface: AbstractSampler

"""
    struct Model{pvars, dvars, F, TData, TDefaults}
        f::F
        data::TData
        defaults::TDefaults
    end

A `Model` struct with parameter variables `pvars`, data variables `dvars`, inner
function `f`, `data::NamedTuple` and `defaults::NamedTuple`.
"""
struct Model{pvars,
    dvars,
    F,
    TData,
    TDefaults
} <: Sampleable{VariateForm,ValueSupport} # May need to find better types
    f::F
    data::TData
    defaults::TDefaults
end
function Model{pvars, dvars}(f::F, data::TD, defaults::TDefaults) where {pvars, dvars, F, TD, TDefaults}
    return Model{pvars, dvars, F, TD, TDefaults}(f, data, defaults)
end
get_pvars(m::Model{params}) where {params} = Tuple(params.types)
get_dvars(m::Model{params, data}) where {params, data} = Tuple(data.types)
get_defaults(m::Model) = m.defaults
@generated function in_pvars(::Val{sym}, ::Model{params}) where {sym, params}
    return sym in params.types ? :(true) : :(false)
end
@generated function in_dvars(::Val{sym}, ::Model{params, data}) where {sym, params, data}
    return sym in data.types ? :(true) : :(false)
end
(model::Model)(args...; kwargs...) = model.f(args..., model; kwargs...)
function runmodel! end
function getspace end

struct Selector
    gid :: UInt64
    tag :: Symbol # :default, :invalid, :Gibbs, :HMC, etc.
end
Selector() = Selector(time_ns(), :default)
Selector(tag::Symbol) = Selector(time_ns(), tag)
hash(s::Selector) = hash(s.gid)
==(s1::Selector, s2::Selector) = s1.gid == s2.gid

"""
Robust initialization method for model parameters in Hamiltonian samplers.
"""
struct SampleFromUniform <: AbstractSampler end
struct SampleFromPrior <: AbstractSampler end

getspace(::Union{SampleFromPrior, SampleFromUniform}) = ()
getspace(::Type{<:SampleFromPrior}) = ()
getspace(::Type{<:SampleFromUniform}) = ()

"""
An abstract type that mutable sampler state structs inherit from.
"""
abstract type AbstractSamplerState end

"""
    Sampler{T}

Generic interface for implementing inference algorithms.
An implementation of an algorithm should include the following:

1. A type specifying the algorithm and its parameters, derived from InferenceAlgorithm
2. A method of `sample` function that produces results of inference, which is where actual inference happens.

Turing translates models to chunks that call the modelling functions at specified points.
The dispatch is based on the value of a `sampler` variable.
To include a new inference algorithm implements the requirements mentioned above in a separate file,
then include that file at the end of this one.
"""
mutable struct Sampler{T, S<:AbstractSamplerState} <: AbstractSampler
    alg      ::  T
    info     ::  Dict{Symbol, Any} # sampler infomation
    selector ::  Selector
    state    ::  S
end
Sampler(alg) = Sampler(alg, Selector())
Sampler(alg, model::Model) = Sampler(alg, model, Selector())
Sampler(alg, model::Model, s::Selector) = Sampler(alg, model, s)

include("utilities/Utilities.jl")
using .Utilities
include("core/Core.jl")
using .Core
include("inference/Inference.jl")  # inference algorithms
using .Inference
include("variational/VariationalInference.jl")
using .Variational

# TODO: re-design `sample` interface in MCMCChains, which unify CmdStan and Turing.
#   Related: https://github.com/TuringLang/Turing.jl/issues/746
#@init @require CmdStan="593b3428-ca2f-500c-ae53-031589ec8ddd" @eval begin
#     @eval Utilities begin
#         using ..Turing.CmdStan: CmdStan, Adapt, Hmc
#         using ..Turing: HMC, HMCDA, NUTS
#         include("utilities/stan-interface.jl")
#     end
# end

@init @require LogDensityProblems="6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c" @eval Inference begin
    using ..Turing.LogDensityProblems: LogDensityProblems, AbstractLogDensityProblem, ValueGradient
    struct FunctionLogDensity{F} <: AbstractLogDensityProblem
        dimension::Int
        f::F
    end

    LogDensityProblems.dimension(ℓ::FunctionLogDensity) = ℓ.dimension

    function LogDensityProblems.logdensity(
        ::Type{ValueGradient},
        ℓ::FunctionLogDensity,
        x::AbstractVector,
    )
        return ℓ.f(x)::ValueGradient
    end
end
@init @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" @eval Inference begin
    using ..Turing.DynamicHMC: DynamicHMC, NUTS_init_tune_mcmc
    include("contrib/inference/dynamichmc.jl")
end

# Random probability measures.
include("stdlib/distributions.jl")
include("stdlib/RandomMeasures.jl")

###########
# Exports #
###########

# Turing essentials - modelling macros and inference algorithms
export  @model,                 # modelling
        @VarName,

        MH,                     # classic sampling
        Gibbs,

        HMC,                    # Hamiltonian-like sampling
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,
        DynamicNUTS,
        ANUTS,

        IS,                     # particle-based sampling
        SMC,
        CSMC,
        PG,
        PIMH,
        PMMH,
        IPMCMC,

        vi,                    # variational inference
        ADVI,

        sample,                 # inference
        psample,
        setchunksize,
        resume,

        auto_tune_chunk_size!,  # helper
        setadbackend,
        setadsafe,

        turnprogress,           # debugging

        Flat,
        FlatPos,
        BinomialLogit,
        VecBinomialLogit,
        OrderedLogistic

end


text:
using GSL
using LinearAlgebra

@testset "spherical bessel functions" begin
    x = 1.1 + 1.1im
    @test sbesselj(1, x) ≈ sin(x)/x^2 - cos(x)/x
    @test shankelh1(1, x) ≈ - exp(im*x) * (x + im) / (x^2)

    x = 3.1 + 1.1im
    n = 1
    @test 2 * diffsbessel(shankelh1,n,x) ≈ shankelh1(n-1, x) - (shankelh1(n, x) + x*shankelh1(n+1, x))/x

    x = 3.1 - 2.1im
    n = 2
    @test 2 * diffsbessel(shankelh1,n,x) ≈ shankelh1(n-1, x) - (shankelh1(n, x) + x*shankelh1(n+1, x))/x

    xs = rand(20)*10;
    n = 15; ns = 0:n

    # Compare with spherical bessel from GSL, which can only accept real arguments..
    errs = map(xs) do x
        maximum(abs.(sf_bessel_jl_array(n,x) - [sbesselj(n,x) for n in ns]))
    end;
    @test maximum(errs) < 20*eps(Float64)

    errs = map(xs) do x
        maximum(
            abs.(
                im .* sf_bessel_yl_array(n,x) -
                [shankelh1(n,x) for n in ns] +
                sf_bessel_jl_array(n,x)
            ) ./ abs.(sf_bessel_yl_array(n,x))
        )
    end;
    @test maximum(errs) < 1e4*eps(Float64)

end

@testset "Associated legendre functions" begin
    x = rand(1)[1] * 0.99
    l_max = 2

    # a few associated legendre functions without the Condon-Shortly factor
    Plm_arr = [1,x,sqrt(1-x^2),(3x^2-1)/2, 3x*sqrt(1-x^2),3*(1-x^2)]

    ls, ms = associated_legendre_indices(l_max)

    @test sf_legendre_array(GSL_SF_LEGENDRE_NONE, l_max, x)[1:length(ls)] ≈ Plm_arr

    sph_factors = map(eachindex(ls)) do i
       (-1)^ms[i] * sqrt((2ls[i] + 1)/(4pi) * factorial(ls[i]-ms[i]) / factorial(ls[i]+ms[i]))
    end

    condon_phase = (-1).^ms

    @test condon_phase .* sf_legendre_array(GSL_SF_LEGENDRE_SPHARM, l_max, x)[1:length(ls)] ≈ sph_factors .* Plm_arr

end

@testset "Spherical harmonics" begin
    θ = rand(1)[1] * 0.99
    φ = rand(1)[1] * 0.99

    l_max = 8 # small l_max due to factorial formula below

    ls, ms = spherical_harmonics_indices(l_max)
    sphs = spherical_harmonics(l_max, θ, φ)

    # check special case l == abs(m)
    inds = findall(ls .== abs.(ms))

    for i in inds
        @test sphs[i] ≈ (sign(ms[i]))^ls[i] / (2^ls[i] * factorial(ls[i])) *
            sqrt(factorial(2*ls[i] + 1) / (4pi)) * sin(θ)^ls[i] * exp(im * ms[i] * φ)
    end

    l_max = 30
    ls, ms = spherical_harmonics_indices(l_max)
    sphs = spherical_harmonics(l_max, θ, φ)

    # special case m == 0, reduce to just Legendre polynomials
    Ps = sf_legendre_Pl_array(l_max, cos(θ))
    inds = findall(ms .== 0)

    for l in 0:l_max
        @test sphs[inds][l+1] ≈ sqrt((2l+1)/(4pi)) * Ps[l+1]
    end

    #sphs[inds] .≈ sqrt.((2 .* (0:l_max) .+ 1) ./ (4pi)) .* Ps

    # special case, north pole
    θ = 0.0
    sphs = spherical_harmonics(l_max, θ, φ)

    for i in eachindex(sphs)
        if ms[i] == 0
            @test sphs[i] ≈ sqrt((2*ls[i] + 1)/(4pi))
        else
            @test sphs[i] ≈ 0.0
        end
        i += 1
    end

end

@testset "Gaunt and Wigner symbols" begin
    l1 = rand(1:100)
    l2 = rand(1:100)
    l3 = rand(abs(l1-l2):2:(l1+l2)) # guarantees that iseven(l1+l2+l3)

    @test iseven(l1+l2+l3)

    m1 = rand(-l1:l1)
    m2 = rand(-l2:l2)
    m3 = m1-m2

    @test_throws(DomainError,gaunt_coefficients(l1,2*l1,l2,m2,l3,m3))
    @test_throws(MethodError,gaunt_coefficients(l1,m1,l2,m2,0.1,m3))

    # the spherical harmonics linearisation formula
    θ, φ = rand(2) .* pi

    l_small = 6
    l_max = 2*l_small # needs to be larger than l_small

    ls, ms = spherical_harmonics_indices(l_max);
    Ys = spherical_harmonics(l_max, θ, φ);

    cs = reshape(
        [
            gaunt_coefficients(l1,m1,l2,m2,l3,m3)
        for l3 = 0:l_max for m3 = -l3:l3 for l2 = 0:l_small for m2 = -l2:l2 for l1 = 0:l_small for m1 = -l1:l1]
    , ((l_small + 1)^2,(l_small + 1)^2,(l_max+1)^2));

    # cs = reshape(
    #     [
    #         (l1,m1,l2,m2,l3,m3)
    #     for l3 = 0:l_small for m3 = -l3:l3 for l2 = 0:l_small for m2 = -l2:l2 for l1 = 0:l_max for m1 = -l1:l1]
    # , ((l_small + 1)^2,(l_small + 1)^2,(l_max+1)^2));

    for n2 in 1:(l_small + 1)^2, n3 in 1:(l_small + 1)^2
        @test 4pi * (-1)^(ms[n3]+ms[n2]) * (1.0im)^(ls[n3]-ls[n2]) * Ys[n3] * conj(Ys[n2]) ≈
            sum( (1.0im).^(-ls) .* (-1.0).^ms .* conj.(Ys) .* cs[n2,n3,:]) atol = 1e-12
    end

end


text:
#####
##### Utilities to batch oracles
#####

"""
Utilities for batchifying oracle evaluation.

The typical workflow should be the following:

  1. Launch an inference server using `launch_server` to obtain a request channel.
  2. Create several `BatchedOracle` that send queries along this channel.
  3. When a client to the server is done, it calls `client_done!` so that the server
     does not have to wait for its next query anymore.
"""
module Batchifier

using ..AlphaZero: MCTS, Util

export BatchedOracle

"""
    launch_server(f; num_workers, batch_size)

Launch an inference requests server.

# Arguments

  - `f` is the function to be evaluated (e.g. a neural network).
    It takes a batch of inputs and returns a batch of results
    of similar size.
  - `num_workers` is the number of workers that are expected to
    query the server.
  - `batch_size` corresponds to the batch size that is used to evaluate `f`.
    Note that one must have `batch_size <= num_workers`

# How to use

This function returns a channel along which workers can send queries.
A query can be either:

    - `:none` when a worker is done and about to terminate
    - a named tuple with fields `query` (the input to be given to `f`) and
      `answer_channel` (the channel the sever must use to return its answer).

The server stops automatically after all workers send `:none`.
"""
function launch_server(f; num_workers, batch_size)
  @assert batch_size <= num_workers
  channel = Channel(num_workers)
  # The server is spawned on the main thread for maximal responsiveness
  Util.@tspawn_main Util.@printing_errors begin
    num_active = num_workers
    pending = []
    while num_active > 0
      req = take!(channel)
      if req == :done
        num_active -= 1
        if num_active < batch_size
          batch_size = num_active
        end
      else
        push!(pending, req)
      end
      @assert length(pending) <= num_active
      @assert batch_size <= num_active
      if length(pending) >= batch_size && length(pending) > 0
        batch = [p.query for p in pending]
        results = f(batch)
        for i in eachindex(pending)
          put!(pending[i].answer_channel, results[i])
        end
        empty!(pending)
      end
    end
  end
  return channel
end

"""
    client_done!(reqc)

This function is to be called every time a client to the inference server identified
by request channel `reqc` is done and won't send queries anymore.

It should be called ``n`` times in total, where ``n`` is the number of workers that was
indicated to [`launch_server`](@ref).

!!! note

    This function does not take a `BatchedOracle` as an argument as there can be several
    oracles per worker (e.g. a worker can simulate games between two players that each
    rely on a neural network).
"""
client_done!(reqc) = put!(reqc, :done)

"""
    BatchedOracle(reqc, preprocess=(x->x))

Create an oracle that delegates its job to an inference server.

- When called on a state, this oracle sends a query to the server identified by
  request channel `reqc`. This call is blocking until every other active worker also
  sends its query.
- A `preprocess` function can be provided to ransform the passed state before it is
  sent to the server as a query.
"""
struct BatchedOracle{F}
  preprocess :: F
  reqchan :: Channel
  anschan :: Channel
  function BatchedOracle(reqchan, preprocess=(x->x))
    return new{typeof(preprocess)}(preprocess, reqchan, Channel(1))
  end
end

function (oracle::BatchedOracle)(state)
  query = oracle.preprocess(state)
  put!(oracle.reqchan, (query=query, answer_channel=oracle.anschan))
  answer = take!(oracle.anschan)
  return answer
end

end