text:
#  using PortAudio;
using FileIO;
using DSP;
using WAV;
using ProgressMeter;

using AdalmPluto;

"""
FM demodulation
"""
function fmDemod(sig)
    out = zeros(Float32,length(sig));
    @inbounds @simd for n ∈ (1:length(sig)-1)
        out[n+1] = angle(sig[n+1]*conj(sig[n]));
    end
    return out;
end

"""
Records a `seconds` long FM sample in (pwd=.../AdalmPluto/examples)/samples/fm.wav
Returns 0 if no error happens.
"""
function fmRecord(seconds, carrierFreq::Int64 = Int64(96e6), samplingRate::Int64 = Int64(2.8e6))
    if last(split(pwd(), "/")) == "AdalmPluto.jl"; cd("./examples"); end;
    dirpath = joinpath(pwd(), "samples");
    filepath = joinpath(dirpath, "fm.wav");
    if !isdir(dirpath)
        mkdir(dirpath);
    end

    audioRendering = 48e3;
    # initialize audio file with a beep
    t = 0.0:1/audioRendering:prevfloat(1.0);
    f = 1e3;
    y = sin.(2pi * f * t) * 0.1;
    wavwrite(y, joinpath(pwd(), "samples", "fm.wav"), Fs=audioRendering);

    # fm demodulation components
    quadRate = 384e3;
    decim1 = Int(samplingRate ÷ quadRate);
    h1 = digitalfilter(Lowpass(44100, fs=samplingRate), FIRWindow(hamming(128)));
    decim2 = Int(quadRate ÷ audioRendering);
    h2 = digitalfilter(Lowpass(Int(audioRendering ÷ 2), fs=quadRate), FIRWindow(hamming(128)));

    try
        # init pluto
        global pluto = openPluto(carrierFreq, samplingRate, 64);

        AdalmPluto.ad9361_baseband_auto_rate(C_iio_context_find_device(pluto.ctx, "ad9361-phy"), Int(samplingRate));

        # demodulation
        rawSamples = zeros(ComplexF32, 1024*1024);
        wavSmpCount = 0;
        progress = Progress(Int(seconds*audioRendering), barglyphs=BarGlyphs("[=> ]"));
        while(wavSmpCount < seconds*audioRendering)
            n = recv!(rawSamples, pluto);
            samples = filt(h1, rawSamples);
            samples = samples[1+64:decim1:end-64];
            samples = fmDemod(samples);
            samples = filt(h2, samples);
            wavSamples = samples[1+64:decim2:end-64];
            wavSmpCount += length(wavSamples);
            update!(progress, Int(wavSmpCount));
            wavappend(wavSamples, filepath);
        end
        println(""); # progress bar fix
    catch e
        # free allocated C stuff (needed to access pluto again)
        AdalmPluto.close(pluto);
        rethrow(e);
    end

    AdalmPluto.close(pluto);

    return 0;
end

fmRecord(10);


text:
# High-level convenience functions to make standard models
# Note: When adding a function here, also add a method taking an AbstractSystem
#       to external/atomsbase.jl

"""
Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model.
Use `extra_terms` to add additional terms.
"""
function model_atomic(lattice::AbstractMatrix,
                      atoms::Vector{<:Element},
                      positions::Vector{<:AbstractVector};
                      extra_terms=[], kwargs...)
    @assert !(:terms in keys(kwargs))
    terms = [Kinetic(),
             AtomicLocal(),
             AtomicNonlocal(),
             Ewald(),
             PspCorrection(),
             extra_terms...]
    if :temperature in keys(kwargs) && kwargs[:temperature] != 0
        terms = [terms..., Entropy()]
    end
    Model(lattice, atoms, positions; model_name="atomic", terms, kwargs...)
end

"""
Build a DFT model from the specified atoms, with the specified functionals.
"""
function model_DFT(lattice::AbstractMatrix,
                   atoms::Vector{<:Element},
                   positions::Vector{<:AbstractVector},
                   xc::Xc;
                   extra_terms=[], kwargs...)
    model_name = isempty(xc.functionals) ? "rHF" : join(xc.functionals, "+")
    model_atomic(lattice, atoms, positions;
                 extra_terms=[Hartree(), xc, extra_terms...], model_name, kwargs...)
end
function model_DFT(lattice::AbstractMatrix,
                   atoms::Vector{<:Element},
                   positions::Vector{<:AbstractVector},
                   functionals;
                   kwargs...)
    model_DFT(lattice, atoms, positions, Xc(functionals); kwargs...)
end

"""
Build an LDA model (Teter93 parametrization) from the specified atoms.
"""
function model_LDA(lattice::AbstractMatrix, atoms::Vector{<:Element},
                   positions::Vector{<:AbstractVector}; kwargs...)
    model_DFT(lattice, atoms, positions, [:lda_x, :lda_c_pw]; kwargs...)
end

"""
Build an PBE-GGA model from the specified atoms.
DOI:10.1103/PhysRevLett.77.3865
"""
function model_PBE(lattice::AbstractMatrix, atoms::Vector{<:Element},
                   positions::Vector{<:AbstractVector}; kwargs...)
    model_DFT(lattice, atoms, positions, [:gga_x_pbe, :gga_c_pbe]; kwargs...)
end

"""
Build a SCAN meta-GGA model from the specified atoms.
DOI:10.1103/PhysRevLett.115.036402
"""
function model_SCAN(lattice::AbstractMatrix, atoms::Vector{<:Element},
                    positions::Vector{<:AbstractVector}; kwargs...)
    model_DFT(lattice, atoms, positions, [:mgga_x_scan, :mgga_c_scan]; kwargs...)
end

ext: jl


text:
module Test18

using Automa
using Automa.RegExp: @re_str
using Test

@testset "Test18" begin
    machine = Automa.compile(re"\0\a\b\t\n\v\r\x00\xff\xFF[\\][^\\]")
    for generator in (:table, :inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        init_code = Automa.generate_init_code(ctx, machine)
        exec_code = Automa.generate_exec_code(ctx, machine)
        validate = @eval function (data)
            $(init_code)
            p_end = p_eof = lastindex(data)
            $(exec_code)
            return cs == 0
        end
        @test validate(b"abracadabra") == false
        @test validate(b"\0\a\b\t\n\v\r\x00\xff\xFF\\!") == true
        @test validate(b"\0\a\b\t\n\v\r\x00\xff\xFF\\\\") == false
    end
end

end


text:
# This script tests the internal functions of eqcond.jl
using Test, OrderedCollections, BenchmarkTools
include("../../../src/includeall.jl")

time_functions = false
rp = joinpath(dirname(@__FILE__), "../../reference/models/li2020")

# Set up input arguments and expected behavior
inside_input = load(joinpath(rp, "inside_iteration_input.jld2"))
inside_step3 = load(joinpath(rp, "inside_iteration_Step3.jld2"))
inside_output = load(joinpath(rp, "inside_iteration_loop_output.jld2"))
Q̂_input = load(joinpath(rp, "Qhat_calculation_input.jld2"))
prepQ̂_input = load(joinpath(rp, "post_processing_input.jld2"))
aug_out = load(joinpath(rp, "augment_variables_out.jld2"), "results")
Q̂_out   = load(joinpath(rp, "hat_Q_calculation_output.jld2"))

@testset "Inside Iteration" begin
    m = Li2020()
    stategrid, funcvar, derivs, endo = initialize!(m)
    θ = parameters_to_named_tuple(m)
    p₀ = aug_out["p0"]
    p₁ = aug_out["p1"]
    max_jump = (p₁ - p₀) / p₁
    p_fitted = extrapolate(interpolate((stategrid[:w], ), vec(inside_input["p_last"]), get_setting(m, :p_interpolant)), Line())
    Q = get_setting(m, :avg_gdp) * get_setting(m, :gov_bond_gdp_level)
    Φ = get_setting(m, :Φ)
    κp_grid = get_setting(m, :κp_grid)
    κp, xK, xg, a, b, c = inside_iteration_li2020(inside_input["w"], p_fitted, inside_input["p"], inside_input["pw"], inside_input["xK"],
                                                  inside_input["xg"], Q, inside_input["hat_Q"], max_jump, θ, κp_grid = κp_grid,
                                                  nlsolve_tol = 1e-6, xg_tol = 2e-3, verbose = :high)

    # First test one loop
    @test a
    @test b
    @test c
    @test κp ≈ inside_step3["kappap"] atol=1e-7
    @test xK ≈ inside_step3["xK"] atol=1e-5
    @test xg ≈ inside_step3["xg"] atol=1e-7

    # Now test the whole loop
    p = vec(inside_input["p_last"])
    Q̂ = fill(inside_input["hat_Q"], length(p))
    xg = Q ./ (2. .* stategrid[:w] .* p)
    xg[1] = 2 * xg[2]
    xg[end] = θ[:ϵ]
    ∂p∂w = differentiate(stategrid[:w], p)
    ψ  = map((x, y) -> (Φ(x, θ[:χ], θ[:δ]) + θ[:ρ] * (x + y) - θ[:AL]) / (θ[:AH] - θ[:AL]), p, Q̂)
    xK = ψ ./ stategrid[:w] .* p ./ (p + Q̂)

    solved_κp = BitArray(undef, length(p))
    solved_xK = BitArray(undef, length(p))
    solved_xg = BitArray(undef, length(p))
    solved_xK[1] = true
    solved_κp[1] = true
    solved_xg[1] = true
    solved_xK[end] = true
    solved_κp[end] = true
    solved_xg[end] = true
    κp_new = zeros(length(p))
    xg_new = deepcopy(xg)
    p_new = deepcopy(p)
    for i in 2:(length(stategrid) - 1)
        wᵢ    = stategrid[:w][i]
        pᵢ    = p[i]
        xKᵢ   = xK[i]
        xgᵢ   = xg[i]
        ∂p∂wᵢ = ∂p∂w[i]
        Q̂ᵢ    = Q̂[i]

        κp_new[i], xKᵢ, xg_new[i], succeed_κp, succeed_xK, succeed_xg =
            inside_iteration_li2020(wᵢ, p_fitted, pᵢ, ∂p∂wᵢ, xKᵢ, xgᵢ, Q, Q̂ᵢ, max_jump, θ;
                                    κp_grid = κp_grid, nlsolve_tol = 1e-6, xg_tol = 2e-3, nlsolve_iter = 400)

        ψᵢ = xKᵢ * wᵢ / (pᵢ / (pᵢ + Q̂ᵢ))
        p_new[i] = nlsolve(x -> Φ(x[1], θ[:χ], θ[:δ]) .+ θ[:ρ] .* (x .+ Q̂ᵢ) .- (ψᵢ * θ[:AH] + (1. - ψᵢ) * θ[:AL]),
                           [(p₀ + p₁) / 2.]).zero[1]
        solved_κp[i] = succeed_κp
        solved_xK[i] = succeed_xK
        solved_xg[i] = succeed_xg
    end

    problem_inds = [11, 14, 15, 16, 21, 97]
    nonproblem_inds = setdiff(1:length(p), problem_inds)
    @test vec(inside_output["xg_new"])[nonproblem_inds] ≈ xg_new[nonproblem_inds] atol=9e-5
    @test abs(inside_output["xg_new"][11] - xg_new[11]) < 1.32
    @test abs(inside_output["xg_new"][14] - xg_new[14]) < 0.6
    @test abs(inside_output["xg_new"][15] - xg_new[15]) < 0.31
    @test abs(inside_output["xg_new"][16] - xg_new[16]) < 0.03
    @test abs(inside_output["xg_new"][21] - xg_new[21]) < 2e-4
    @test abs(inside_output["xg_new"][97] - xg_new[97]) < 2e-3
    @test vec(inside_output["p_new"]) ≈ p_new atol=8e-5
    @test vec(inside_output["kappap_vec"]) ≈ κp_new atol=5.5e-4
    @test all(solved_κp)
    @test all(solved_xK)
    @test all(solved_xg)
end

@testset "Prepare for and calculate Q̂" begin
    m = Li2020()
    stategrid, funcvar, derivs, endo = initialize!(m)
    θ = parameters_to_named_tuple(m)
    f_μK = get_setting(m, :μK)
    Φ = get_setting(m, :Φ)
    yg_tol = get_setting(m, :yg_tol)
    firesale_bound = get_setting(m, :firesale_bound)
    firesale_interpolant = get_setting(m, :firesale_interpolant)
    Q = get_setting(m, :avg_gdp) * get_setting(m, :gov_bond_gdp_level)
    funcvar[:p] = vec(prepQ̂_input["p_new"])
    funcvar[:Q̂] = vec(prepQ̂_input["hat_Q_last"])
    funcvar[:xg] = vec(prepQ̂_input["xg_new"])
    endo[:κp] = vec(prepQ̂_input["kappap_vec"])
    prepare_Q̂!(stategrid, funcvar, derivs, endo, θ, f_μK, Φ, yg_tol, firesale_bound, firesale_interpolant, Q)

    @test @test_matrix_approx_eq vec(aug_out["muw_vec"]) endo[:μw]
    @test @test_matrix_approx_eq vec(aug_out["sigmaw_vec"]) endo[:σw]
    @test @test_matrix_approx_eq vec(aug_out["kappaw_vec"]) endo[:κw]
    @test @test_matrix_approx_eq vec(aug_out["rg_vec"]) endo[:rg]
    @test @test_matrix_approx_eq vec(aug_out["rh_vec"]) endo[:rh]
    @test @test_matrix_approx_eq vec(aug_out["rf_vec"]) endo[:rf]

    if time_functions
        BenchmarkTools.@btime begin
            Q̂_calculation!(stategrid, zeros(length(stategrid)), vec(Q̂_input["muw_vec"]), vec(Q̂_input["sigmaw_vec"]),
                           vec(Q̂_input["kappaw_vec"]), vec(Q̂_input["rf_vec"]), vec(Q̂_input["rg_vec"]), vec(Q̂_input["rh_vec"]),
                           Q̂_input["Qval"], Q̂_input["lambda"])
        end
    end

    my_Q̂ = Q̂_calculation!(stategrid, zeros(length(stategrid)), vec(Q̂_input["muw_vec"]), vec(Q̂_input["sigmaw_vec"]),
                          vec(Q̂_input["kappaw_vec"]), vec(Q̂_input["rf_vec"]), vec(Q̂_input["rg_vec"]), vec(Q̂_input["rh_vec"]),
                          Q̂_input["Qval"], Q̂_input["lambda"])
    @test maximum(abs.(my_Q̂ - vec(Q̂_out["hat_Q_vec"]))) < 5e-4
end

@testset "eqcond" begin
    m = Li2020()
    Q = get_setting(m, :avg_gdp) * get_setting(m, :gov_bond_gdp_level)
    stategrid, funcvar, derivs, endo = initialize!(m)
    bc = get_setting(m, :boundary_conditions)
    bc[:p] .= [aug_out["p0"], aug_out["p1"]] # Change boundary conditions to match MATLAB output
    indiv_convergence = OrderedDict{Symbol, Vector{eltype(m)}}(:Q̂ => [0., 1e-3])
    p = vec(inside_input["p_last"])
    Q̂ = fill(inside_input["hat_Q"], length(p))
    θ = parameters_to_named_tuple(m)
    xg = Q ./ (2. .* stategrid[:w] .* p)
    xg[1] = 2 * xg[2]
    xg[end] = θ[:ϵ]
    funcvar[:p] .= p
    funcvar[:xg] .= xg
    funcvar[:Q̂] .= Q̂
    θ = parameters_to_named_tuple(m)
    f = eqcond(m)
    @info "The following output is expected; checking high verbosity statements are shown."
    p_new, xg_new, Q̂_new = f(stategrid, funcvar, derivs, endo, θ; individual_convergence = indiv_convergence,
                             verbose = :high)
    # Very rough bounds
    @test p_new ≈ aug_out["p_new"] atol=1e-2
    @test p_new[1:41] ≈ aug_out["p_new"][1:41] atol=2e-3
    @test p_new[42:48] ≈ aug_out["p_new"][42:48] atol=6e-3
    @test p_new[49:end - 1] ≈ aug_out["p_new"][49:end - 1] atol=1e-6
    @test p_new[end] == p_new[end - 1]
    @test Q̂_new ≈ Q̂_out["hat_Q_vec"] atol=1e-2
    @test xg_new[1:12] ≈ inside_output["xg_new"][1:12]
    @test xg_new[13:19] ≈ inside_output["xg_new"][13:19] atol=1.5
    @test xg_new[20:30] ≈ inside_output["xg_new"][20:30] atol=8e-2
    @test xg_new[30:end] ≈ inside_output["xg_new"][30:end] atol=5e-3

    if time_functions
        BenchmarkTools.@btime begin # roughly 3s vs. 30s of MATLAB (2.8s if not starting from guess of zeros for Q̂)
            p = vec(inside_input["p_last"])
            Q̂ = fill(inside_input["hat_Q"], length(p))
            θ = parameters_to_named_tuple(m)
            xg = Q ./ (2. .* stategrid[:w] .* p)
            xg[1] = 2 * xg[2]
            xg[end] = θ[:ϵ]
            funcvar[:p] .= p
            funcvar[:xg] .= xg
            funcvar[:Q̂] .= Q̂
            indiv_convergence = OrderedDict{Symbol, Vector{eltype(m)}}(:Q̂ => [0., 1e-3])
            f(stategrid, funcvar, derivs, endo, θ; individual_convergence = indiv_convergence)
        end
    end
end

nothing


text:
#=

# Vulkan types

The API possesses a few data structures that exhibit a different behavior. We list all of these here along with their properties and features that we hope will free the developer of many patterns and obligations of the Vulkan API.

## Handles

Handles are opaque pointers to internal Vulkan objects. Almost every handle must be created and destroyed with API commands. Some handles have a parent handle (see [Parent handle access](@ref Parent-handle-access) for navigating through the resulting handle hierarchy), which *must not* be destroyed before its children. For this we provide wrappers around creation functions with an automatic finalization feature (see [Automatic finalization](@ref Automatic-finalization)) that uses a simple reference couting system. This alleviates the burden of tracking when a handle can be freed and freeing it, in conformance with the Vulkan specification.

Most handles are typically created with a `*CreateInfo` or `*AllocateInfo` structure, that packs creation parameters to be provided to the API creation function. To allow for nice one-liners that don't involve long create info names, [these create info parameters are exposed](@ref expose-create-info-args) in the creation function, automatically building the create info structure for you.

!!! tip
    Handle types have a constructor defined that wraps around the creation function and automatically unwrap the result (see [Error handling](@ref Error-handling)).

### Automatic finalization

Handles can be created with the API functions `vkCreate*` and `vkAllocate*`, and most of them must be destroyed after use with a `vkDestroy*` or `vkFree*`. More importantly, they must be destroyed with the same allocator and parent handle that created them. To automate this, new mutable handle types were defined to allow for the registration of a [finalizer](https://docs.julialang.org/en/v1/base/base/#Base.finalizer). The `create_*` and `allocate_*` wrappers automatically register the corresponding destructor in a finalizer, so that you don't need to worry about destructors (except for `CommandBuffer`s and `DescriptorSet`s, see below). A handle's finalizer, and therefore its API destructor, will execute when there are no program-accessible references to this handle.

`CommandBuffer`s and `DescriptorSet`s do not register any destructor and are never implicitly freed. You will have to explicitly free those resources yourself with `free_command_buffers` and `free_descriptor_sets` respectively. The reason for that is that they are supposed to be freed in batches for performance considerations.

Therefore, you should **never** explicitly call a destructor (`destroy_*` or `free_*`) except for `CommandBuffer` and `DescriptorSet`. Otherwise objects will be destroyed/freed twice, likely resulting in a crash.

Finalizers [can be run with](https://docs.julialang.org/en/v1/base/base/#Base.finalize) `finalize(handle)`, which may be useful in situations where you want to free an object without waiting for finalizers to be triggered by the Julia runtime. This can be done safely: the finalizers won't run twice if you run them manually.

However, finalizers can be run in arbitrary order, and some handles require to be destroyed only after all their children (such as `VkDevice`s). To avoid crashes related to bad finalization execution order, a simple thread-safe reference counting system is used to make sure that a handle is destroyed **only after all its children are destroyed**.

!!! note
    If you need to construct a handle from an opaque pointer (obtained, for example, via an external library such as a `VkSurfaceKHR` from GLFW), you can use the constructor `(::Type{<:Handle})(ptr::Ptr{Cvoid}, destructor[, parent])` as in

    ```julia
    surface_ptr = GLFW.CreateWindowSurface(instance, window)
    SurfaceKHR(surface_ptr, x -> destroy_surface_khr(instance, x), instance)
    ```

    If you don't need to destroy the surface (for example, if the external library does it for you), then you can just pass in the `identity` function as a destructor.

This introduces a small overhead, since the parent handle and allocator are stored in an anonymous function for each handle at creation. However, it should be minor compared to the execution time of the API destructors.

Because finalization order is the source of many Vulkan bugs, particularly when objects implicitly depend on other objects being alive, there is a [preference](@ref Preferences) `LOG_DESTRUCTION` that allows you to log all destructions if set to `"true"` for debugging purposes.

### [Exposition of create info arguments](@id expose-create-info-args)

Handles that can only be created with a single API constructor possess an additional constructor that wraps around the generated create/allocate\* functions, building the required \*\[Create/Allocate\]Info from exposed arguments. That way, you do not have to explicitly construct this intermediate structure, which reduces boilerplate code.

For example

=#

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

const instance = Instance([], [])
const physical_device = first(unwrap(enumerate_physical_devices(instance)))
const device = Device(
    physical_device, [DeviceQueueCreateInfo(0, [1.])],
    [],
    [],
)

fence = unwrap(create_fence(device, FenceCreateInfo()))
fence_signaled = unwrap(
                   create_fence(
                     device,
                     FenceCreateInfo(flags=FENCE_CREATE_SIGNALED_BIT);
                   )
                 )

# can be replaced with

fence = unwrap(create_fence(device))
fence_signaled = unwrap(
                   create_fence(
                     device;
                     flags=FENCE_CREATE_SIGNALED_BIT,
                     )
                   )

#=

Note that we `unwrap` the result every time, assuming that the `create_fence` function did not return any error. See the [error handling](@ref Error-handling) section for more information.

Furthermore, handle types have a generated constructor that exposes the same arguments as the create/allocate\* functions, but automatically unwrapping the result so you don't have to call it manually. The above can then be further reduced into

=#

fence = Fence(device)
fence_signaled = Fence(
                   device,
                   flags=FENCE_CREATE_SIGNALED_BIT;
                 )

# When multiple info structures are requested, only the main one has its fields exposed as function arguments:

## the ApplicationInfo has to be provided manually
Instance([], []; application_info = ApplicationInfo(v"0.1", v"0.1", v"1.2"))

# Or, if an array of create infos is expected, then you will have to provide it yourself:

## the array of DeviceQueueCreateInfo has to be provided manually
Device(physical_device, [DeviceQueueCreateInfo(0, [1.0])], [], [])

# When multiple handles are constructed at the same time, no additional constructor is defined and you need to call the `create_*` function manually

command_pool = CommandPool(device, 0)
command_buffers = unwrap(
                    allocate_command_buffers(
                      device,
                      CommandBufferAllocateInfo(
                        command_pool,
                        COMMAND_BUFFER_LEVEL_PRIMARY,
                        3,
                      )
                    )
                  )

#=

### Parent handle access

Handles store their parent handle if they have one. For example, it is possible to just carry a `Pipeline` around and access its `device` field whenever you like, and the `physical_device` field of this device and so on until you reach the instance that has no parent. That way, it is not necessary to keep all the other parent handles in scope (e.g. `Device` and `Instance` handles) to make an API call.

=#

command_pool.device.physical_device.instance

# `Base.parent` was extended to work on handles that possess a parent handle.

parent(device) == device.physical_device == physical_device

#=

## Structures

Vulkan structures, such as `Extent2D`, `InstanceCreateInfo` and `PhysicalDeviceFeatures` were wrapped into two different structures each.

First, there are high-level-ish structures that are easier to work with, replacing C types by idiomatic Julia types. The user will almost never need to worry about pointers, and will be able to use types such as Julia arrays, strings and version numbers in a transparent manner. All the data is stored in its fields, and is accessible just like you would expect from a Julia structure. Equality and hashing are implemented with [AutoHashEquals.jl](https://github.com/andrewcooke/AutoHashEquals.jl), to facilitate their use in dictionaries.

However, storing data and converting it to a C-compatible format for API calls introduce some overhead for these high-level structures. It will not matter in non performance-critical sections, but may be a limitation in tight loops. Therefore, lower-level, minimalistic types are available, which only carry the smallest possible amount of information needed to use the API. Although they can be created very similarly to high-level structures, the high-level arguments are directly converted to their C-compatible counterpart. The main drawback is a lack of introspection, since these minimalistic structures carry C types, so these structures are badly suited to keeping data around. Minimalistic structures are not defined for types that are only returned and never requested by Vulkan, such as `ExtensionProperties`.

There is a final family of Vulkan types that you may encounter. Those are the barebones VulkanCore.jl types, which you won't have to worry about in all cases *except when you need to pass functions to the API*. In this case, inputs will not be automatically converted for you, and you will have to define the appropriate signature before obtaining function pointers with `Base.@cfunction`. You can access these types from the (exported) module `Vulkan.core`.

To summarize:

- High-level structs:
  - should be used most of the time.
  - store values in a way that makes it easy to retrieve them later.
  - introduce a small overhead, which may be a concern in some performance-critical sections.

- Low-level structs:
  - offer performance advantages over high-level structs.
  - may be preferred in performance-critical sections.
  - are not meant for introspection capabilities.
  - are not defined for structures not needed by the API.

- VulkanCore structs:
  - should never be used directly, except as argument types for functions intended to be passed to the API.

In general, high-level and low-level structs can be used interchangeably as function arguments to constructors or API functions, at the condition that they are not mixed together.

Using either high-level or low-level structs should be a performance matter, and as such it is encouraged to profile applications before using low-level structs all: they are faster, but can require additional bookkeeping due to a lack of introspection.

Typically, it is easier to use high-level types for create info arguments to handles that are created at a low frequency; this includes `Instance`, `Device` or `SwapchainKHR` handles for example. Their create info structures may contain precious information that needs to be accessed by the application, e.g. to make sure that image formats in a render pass comply with the swapchain image format, or to check instance or device extensions before using extension functionality.

API functions and structures accept either low-level structs or high-level structs. For commands with low-level structs, you currently need to provide typed arrays (i.e. not `[]` which are of type `Vector{Any}`).

In general:
- High-level structs are returned from functions with high-level arguments.
- Low-level structs are returned from functions with low-level arguments.

The only exception currently is for functions that would have the same low-level/high-level argument types, for which only one version is available that returns values in low-level types.

### Conversion

In Vulkan, version numbers are `UInt32` with a special encoding to extract major, minor and patch numbers, booleans are represented as `UInt32` values, `String`s as `Ptr{UInt8}` or `NTuple{N, UInt8}`. Those types have a natural counterpart in Julia, being respectively `VersionNumber`, `Bool` and `String` types. The wrapper makes it possible to work with these idiomatic Julia types, automatically converting them wherever necessary.

### Automatic insertion of inferable arguments

In some places, part of the arguments of a function or of the fields of a structure can only take one logical value. It can be divided into two sets:

1. The structure type `sType` of certain structures
2. Arguments related to the start and length of a pointer which represents an array

The second set is a consequence of using a higher-level language than C. In C, the pointer alone does not provide any information regarding the number of elements it holds. In Julia, array-like values can be constructed in many different ways, being an `Array`, a `NTuple` or other container types which provide a `length` method.

#### Structure type

Many API structures possess a `sType` field which must be set to a unique value. This is done to favor the extendability of the API, but is unnecessary boilerplate for the user. Worse, this is an error-prone process which may lead to crashes. All the constructors of this wrapper do not expose this `sType` argument, and hardcode the expected value.

#### Pointer lengths

The length of array pointers is automatically deduced from the length of the container passed in as argument.

#### Pointer starts

Some API functions require to specify the start of a pointer array as an argument. They have been hardcoded to 0 (first element in C), since it is always possible to pass in a sub-array (e.g. a view).

## Bitmask flags

In Vulkan, the value of some flags carry meaning through a bitmask structure. Bitmasks define bit values which they can be a composition of (using bitwise _and_, _or_, and _xor_ operations). However, the associated flag type is defined as a `UInt32`, which allows any value to be passed in as a flag. This opens up the door to incorrect usage that may be hard to debug. To circumvent that, every bitmask flag now has one associated type which prevents combinations with flags of other bitmask types.

For example, consider the vanilla `VkSampleCountFlags` type (alias for `UInt32`) with bits defined via the enumerated type `VkSampleCountFlagBits`:

```@repl
using Vulkan.core
VK_SAMPLE_COUNT_1_BIT isa VkSampleCountFlagBits
VK_SAMPLE_COUNT_1_BIT === VkSampleCountFlagBits(1)
VK_SAMPLE_COUNT_1_BIT === VkSampleCountFlags(1)
VK_SAMPLE_COUNT_1_BIT | VK_SAMPLE_COUNT_2_BIT === VkSampleCountFlags(3)
VK_SAMPLE_COUNT_1_BIT & VK_SAMPLE_COUNT_2_BIT === VkSampleCountFlags(0)
VK_SAMPLE_COUNT_1_BIT & VK_SURFACE_TRANSFORM_IDENTITY_BIT_KHR === VkSampleCountFlags(1)
```

Those two types are combined into one `SampleCountFlag`:

```@repl
using Vulkan
SampleCountFlag <: BitMask
SurfaceTransformFlagKHR <: BitMask # another bitmask flag
SAMPLE_COUNT_1_BIT | SAMPLE_COUNT_2_BIT === SampleCountFlag(3)
SAMPLE_COUNT_1_BIT & SAMPLE_COUNT_2_BIT === SampleCountFlag(0)
SAMPLE_COUNT_1_BIT & SURFACE_TRANSFORM_IDENTITY_BIT_KHR
UInt32(typemax(SampleCountFlag)) === UInt32(core.VK_SAMPLE_COUNT_FLAG_BITS_MAX_ENUM)
```

All functions that were expecting a `VkSampleCountFlags` (`UInt32`) value will have their wrapped versions expect a value of type `SampleCountFlag`. Furthermore, the `*FLAG_BITS_MAX_ENUM` values are removed. This value is the same for all enums and can be accessed via `typemax(T)` where `T` is a `BitMask` (e.g. `SampleCountFlag`).

=#


text:
# Wrapper function for prediction

"""
    predict(results, h, nChain, Xnew)
Function for forecasting Count Data models.

* `results`: Estimation results
* `h`: Number of steps to forecast
* `nChain`: Number of Chains for simulation based forecast (optional)
* `Xnew`: New values for regressors (only in case of regressors)

# Example
```julia-repl
# 10-step-ahead forecast
predict(results, 10, 10000)
```

The function either returns point forecasts if `nChain` is not specified or
generates multiple time series according to estiamtion results.
The latter is used to compute forecast intervals and is the default for INARMA models.
"""
function predict(results::INGARCHresults,
    h::T where T<:Integer,
    Xnew::Array{T, 2} where T<:AbstractFloat = zeros(0, 0))

    r = length(results.model.external)

    if ndims(Xnew) == 1
        if r > 1
            error("Dimensions of Xnew and external do not match.")
        end

        if length(Xnew) != h
            error("Dimensions of Xnew and h do not match.")
        end

        Xnew = reshape(Xnew, (1, h))
    end

    if length(Xnew) != 0
        if (size(Xnew) == (h, r)) & (h != r)
            Xnew = convert(Array{Float64, 2}, Xnew')
        end

        if size(Xnew)[1] != r
            error("Dimensions of Xnew and external do not match.")
        end

        if size(Xnew)[2] != h
            error("Dimensions of Xnew and external do not match.")
        end
    end

    if ndims(Xnew) != 2
        error("Incorrect type of Xnew.")
    end

    y = results.y
    T = length(y)
    zi = results.model.zi
    nb = results.model.distr == "NegativeBinomial"
    logl = results.model.link == "Log"
    lin = !logl
    pars = results.pars

    if typeof(results.model) != IIDModel
        p = length(results.model.pastObs)
    else
        p = 0
    end

    if typeof(results.model) == INGARCHModel
        q = length(results.model.pastMean)
    else
        q = 0
    end

    if p == 0
        P = 0
    else
        P = maximum(results.model.pastObs)
    end

    if q == 0
        Q = 0
    else
        Q = maximum(results.model.pastMean)
    end

    r = length(results.model.external)

    rI = 0
    rE = 0
    iI = []
    iE = []

    if r > 0
        iE = findall(results.model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    β0 = pars.β0
    α = pars.α
    β = pars.β
    η = pars.η
    ϕ = pars.ϕ
    ω = pars.ω

    λ = zeros(T + h)
    λ[1:T] = results.λ

    if logl
        ν = log.(λ)
    else
        ν = λ
    end

    X = results.model.X
    if rE > 0
        for i = iE
            ν[1:T] = ν[1:T] .- η[i].*X[i, :]
        end
    end

    Y = zeros(T + h)
    Y[1:T] = y

    for t = (T+1):(T+h)
        ν[t] = β0
        if p > 0
            if logl
                ν[t] += sum(α.*log.(Y[t .- results.model.pastObs] .+ 1))
            else
                ν[t] += sum(α.*Y[t .- results.model.pastObs])
            end
        end

        if q > 0
            ν[t] += sum(β.*ν[t .- results.model.pastMean])
        end

        if rI > 0
            ν[t] += sum(η[iI].*Xnew[iI, t - T])
        end

        λ[t] = copy(ν[t])

        if rE > 0
            λ[t] += sum(η[iE].*Xnew[iE, t - T])
        end

        if logl
            λ[t] = exp(λ[t])
        end

        Y[t] = copy(λ[t]*(1 - ω))
    end

    return Y[T+1:end]
end

function predict(results::INGARCHresults,
    h::T where T<:Integer,
    nChain::T where T<:Integer,
    Xnew::Array{T, 2} where T<:AbstractFloat = zeros(0, 0))

    r = length(results.model.external)

    if ndims(Xnew) == 1
        if r > 1
            error("Dimensions of Xnew and external do not match.")
        end

        if length(Xnew) != h
            error("Dimensions of Xnew and h do not match.")
        end

        Xnew = reshape(Xnew, (1, h))
    end

    if length(Xnew) != 0
        if (size(Xnew) == (h, r)) & (h != r)
            Xnew = convert(Array{Float64, 2}, Xnew')
        end

        if size(Xnew)[1] != length(results.model.external)
            error("Dimensions of Xnew and external do not match.")
        end

        if size(Xnew)[2] != h
            error("Dimensions of Xnew and external do not match.")
        end
    end

    if ndims(Xnew) != 2
        error("Incorrect type of Xnew.")
    end

    y = results.y
    T = length(y)
    zi = results.model.zi
    nb = results.model.distr == "NegativeBinomial"
    logl = results.model.link == "Log"
    lin = !logl
    pars = results.pars

    if typeof(results.model) != IIDModel
        p = length(results.model.pastObs)
    else
        p = 0
    end

    if typeof(results.model) == INGARCHModel
        q = length(results.model.pastMean)
    else
        q = 0
    end

    if p == 0
        P = 0
    else
        P = maximum(results.model.pastObs)
    end

    if q == 0
        Q = 0
    else
        Q = maximum(results.model.pastMean)
    end

    M = maximum([P, Q])

    λ = zeros(nChain, M + h)
    ν = zeros(nChain, M + h)

    rI = 0
    rE = 0
    iI = []
    iE = []

    if r > 0
        iE = findall(results.model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    β0 = Float64(pars.β0)
    α = pars.α
    β = pars.β
    if nb
        ϕ = pars.ϕ[1]
    else
        ϕ = 0.0
    end

    η = pars.η
    ω = pars.ω

    Y = zeros(Int64, nChain, M + h)
    λ = zeros(nChain, M + h)
    ν = zeros(nChain, M + h)

    λOld = results.λ

    if logl
        νOld = log.(λOld)
    else
        νOld = λOld
    end

    X = results.model.X
    if rE > 0
        for i = iE
            νOld = νOld .- η[i].*X[i, :]
        end
    end

    for i = 1:nChain
        λ[i, 1:M] = λOld[(end - M + 1):end]
        ν[i, 1:M] = νOld[(end - M + 1):end]
        Y[i, 1:M] = y[(end - M + 1):end]
    end

    @simd for i = 1:nChain
        for t = (M+1):(M+h)
            ν[i, t] = β0
            if p > 0
                if logl
                    ν[i, t] += sum(α.*log.(Y[i, t .- results.model.pastObs] .+ 1))
                else
                    ν[i, t] += sum(α.*Y[i, t .- results.model.pastObs])
                end
            end

            if q > 0
                ν[i, t] += sum(β.*ν[i, t .- results.model.pastMean])
            end

            if rI > 0
                ν[i, t] += sum(η[iI].*Xnew[iI, t - M])
            end

            λ[i, t] = ν[i, t]

            if rE > 0
                λ[i, t] += sum(η[iE].*Xnew[iE, t - M])
            end

            if logl
                λ[i, t] = exp(λ[i, t])
            end

            if nb
                p = ϕ/(ϕ + λ[i, t])
                Y[i, t] = rand(NegativeBinomial(ϕ, p))*(rand() > ω)
            else
                Y[i, t] = rand(Poisson(λ[i, t]))*(rand() > ω)
            end
        end
    end
    Yout = Y[:, M+1:end]

    meanY = mean(Yout, dims = 1)[1, :]
    Q = zeros(2, h)

    for i = 1:h
        Q[:, i] = quantile(Yout[:, i], [0.025, 0.975])
    end

    return meanY, Q, Yout
end

function predict(results::INARMAresults,
    h::T where T<:Integer,
    nChain::T where T<:Integer,
    Xnew::Array{T, 2} where T<:AbstractFloat = zeros(0, 0))

    r = length(results.model.external)

    if ndims(Xnew) == 1
        if length(results.model.external) > 1
            error("Dimensions of Xnew and external do not match.")
        end

        if length(Xnew) != h
            error("Dimensions of Xnew and h do not match.")
        end

        Xnew = reshape(Xnew, (1, h))
    end

    if length(Xnew) != 0
        if (size(Xnew) == (h, r)) & (h != r)
            Xnew = convert(Array{Float64, 2}, Xnew')
        end

        if size(Xnew)[1] != length(results.model.external)
            error("Dimensions of Xnew and external do not match.")
        end

        if size(Xnew)[2] != h
            error("Dimensions of Xnew and external do not match.")
        end
    else
        if length(results.model.external) > 0
            error("Dimensions of Xnew and external do not match.")
        end
    end

    if ndims(Xnew) != 2
        error("Incorrect type of Xnew.")
    end

    y = results.y
    T = length(y)
    zi = results.model.zi
    nb1 = results.model.distr[1] == "NegativeBinomial"
    nb2 = results.model.distr[2] == "NegativeBinomial"
    logl1 = results.model.link[1] == "Log"
    logl2 = results.model.link[2] == "Log"
    lin1 = !logl1
    lin2 = !logl2
    pars = results.pars

    if typeof(results.model) != INMAModel
        p = length(results.model.pastObs)
    else
        p = 0
    end

    if typeof(results.model) != INARModel
        q = length(results.model.pastMean)
    else
        q = 0
    end

    if p == 0
        P = 0
    else
        P = maximum(results.model.pastObs)
    end
    if q == 0
        Q = 0
    else
        Q = maximum(results.model.pastMean)
    end

    M = maximum([P, Q])

    rI = 0
    rE = 0
    iI = []
    iE = []

    if r > 0
        iE = findall(results.model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    β0 = pars.β0
    α = pars.α
    β = pars.β
    if nb1
        ϕ1 = pars.ϕ[1]
    end
    if nb2
        ϕ2 = pars.ϕ[2]
    end
    η = Vector{Float64}(pars.η)
    ω = pars.ω

    Y = fill(0, (nChain, M + h))
    R = fill(0, (nChain, M + h))
    Z = fill(0, (nChain, M + h))

    μ = zeros(M + h)
    λ = fill(β0, M + h)

    # Initialization
    for i = 1:M
        if rE > 0
            for i = iE
                μ[i] += (η[i]*results.model.X[i, T - M + i])
            end
        end

        if rI > 0
            for i = iI
                λ[i] += η[i]*results.model.X[i, T - M + i]
            end
        end

        if nb1
            dR = MixtureModel([NegativeBinomial(ϕ1, ϕ1/(ϕ1 + λ[i])), Poisson(0)], [1 - ω, ω])
        else
            dR = MixtureModel([Poisson(λ[i]), Poisson(0)], [1 - ω, ω])
        end

        if nb2
            dZ = NegativeBinomial(ϕ2, ϕ2/(ϕ2 + μ[i]))
        else
            dZ = Poisson(μ[i])
        end

        for j = 1:nChain
            while true
                Rtemp = rand(dR)
                Ztemp = rand(dZ)
                if Rtemp + Ztemp <= y[T - M + i]
                    R[j, i] = Rtemp
                    Z[j, i] = Ztemp
                    break
                end
            end
            Y[j, i] = y[T - M + i]
        end
    end

    for t = (M+1):(M+h)
        if rE > 0
            for i = iE
                μ[t] += η[i]*Xnew[i, t - M]
            end
        end

        if rI > 0
            for i = iI
                λ[t] += η[i]*Xnew[i, t - M]
            end
        end

        if nb1
            dR = MixtureModel([NegativeBinomial(ϕ1, ϕ1/(ϕ1 + λ[t])), Poisson(0)], [1 - ω, ω])
        else
            dR = MixtureModel([Poisson(λ[t]), Poisson(0)], [1 - ω, ω])
        end

        if nb2
            dZ = NegativeBinomial(ϕ2, ϕ2/(ϕ2 + μ[t]))
        else
            dZ = Poisson(μ[t])
        end

        for i = 1:nChain
            R[i, t] = rand(dR)*(rand() >= ω)
            Z[i, t] = rand(dZ)
            Y[i, t] = R[i, t] + Z[i, t]

            for j = 1:p
                Y[i, t] += α[j]∘Y[i, t .- results.model.pastObs[j]]
            end

            for j = 1:q
                Y[i, t] += β[j]∘R[i, t .- results.model.pastMean[j]]
            end
        end
    end

    Y = Y[:, M+1:end]
    meanY = mean(Y, dims = 1)[1, :]
    Q = zeros(2, h)

    for i = 1:h
        Q[:, i] = quantile(Y[:, i], [0.025, 0.975])
    end

    return meanY, Q, Y
end

function predict(results::INARMAresults,
    h::T where T<:Integer,
    Xnew::Array{T, 2} where T<:AbstractFloat = zeros(0, 0))

    predict(results, h, 10000, Xnew)
end


text:
module FlightMechanicsUtils

using LinearAlgebra
using StaticArrays

export γ_AIR, R_AIR
export gD
export RAD2DEG, DEG2RAD, M2FT, FT2M, KEL2RANK, RANK2KEL, KG2LB, LB2KG, SLUG2KG, KG2SLUG
export PSF2PA, PA2PSF, SLUGFT32KGM3, KGM32SLUGFT3
include("constants.jl")

export atmosphere_isa
include("atmosphere.jl")

export rigid_body_velocity, rigid_body_acceleration
export pqr_2_ψθϕ_dot, pqr_2_quat_dot, ψθϕ_dot_2_pqr
export uvw_to_tasαβ, uvw_dot_to_tasαβ_dot, tasαβ_dot_to_uvw_dot
export rate_of_climb_constrain_no_wind
include("kinematics.jl")

export coordinated_turn_bank, steiner_inertia, translate_forces_moments
include("mechanics.jl")

export euler_angles, quaternions, rotation_matrix_zyx
export body2horizon, body2wind
export wind2body, wind2horizon
export horizon2body, horizon2wind
export horizon2ecef, ecef2horizon
include("rotations.jl")

export Ellipsoid
export Clarke1866, Clarke1880, International, Bessel, Everest, ModifiedEverest,
    AustralianNational, SouthAmerican1969, Airy, ModifiedAiry, Hough, Fischer1960SouthAsia,
    Fischer1960Mercury, Fischer1968, WGS60, WGS66, WGS72, WGS84
export llh2ecef, ecef2llh
include("ellipsoid.jl")

export qc2cas, qc2eas, qc2tas
export cas2eas, cas2tas
export eas2cas, eas2tas
export tas2cas, tas2eas
export compressible_qinf, incompressible_qinf
include("anemometry.jl")

end


text:
module ModelParameters

export AbstractParameter, Parameter, getparameters,
    show, lower, upper, initial, value, constant, convert,
    +, *, -, /, <, isequal, isapprox, one, zero, norm, real, imag,
    isparameter, isfixed, reset!, promote_rule, promote_type

import Base: +, *, -, /, \, ==, <, promote_rule, promote_type, convert,
    push!, isequal, one, oneunit, zero, norm, real, imag, isapprox

abstract type AbstractParameter <: Number end

mutable struct Parameter{T<:Number} <: AbstractParameter
    value::T
    initial::T
    lower::T
    upper::T
    isfixed::Bool
    function Parameter{T}(initial, lower, upper) where T<:Number
        if lower > upper
            throw(ArgumentError(
                "Invalid interval: lower bound is greater than upper bound."))
        elseif initial < lower || initial > upper
            throw(ArgumentError(
                "Initial value is outside the interval."))
        end
        isfixed = lower ≈ upper
        new{T}(initial, initial, lower, upper, isfixed)
    end
end

function Parameter(initial::T, lower::S, upper::R) where
    {T<:Number, S<:Number, R<:Number}
    typ = promote_type(T, promote_type(S, R))
    Parameter{typ}(typ(initial), typ(lower), typ(upper))
end
Parameter(v::T) where {T} = Parameter(v, typemin(T), typemax(T))
Parameter(lower, upper) = Parameter(lower + (upper - lower) / 2, lower, upper)
constant(v) = Parameter(v, v, v)

function push!(prm::Parameter, v)
    if v < lower(prm) || v > upper(prm)
        throw(ArgumentError("$v is outside the interval."))
    end
    prm.value = v
    return prm
end

function reset!(prm::Parameter)
    prm.value = prm.initial
    return prm
end

lower(par::Parameter) = par.lower
upper(par::Parameter) = par.upper
initial(par::Parameter) = par.initial
value(par::Parameter) = par.value
isfixed(prm::Parameter) = prm.isfixed
isparameter(prm::Parameter) = !prm.isfixed
one(::Type{Parameter{T}}) where {T} = constant(one(T))
oneunit(::Type{Parameter{T}}) where {T} = constant(oneunit(T))
zero(::Type{Parameter{T}}) where {T} = constant(zero(T))
one(::Parameter{T}) where {T} = constant(one(T))
oneunit(::Parameter{T}) where {T} = constant(oneunit(T))
zero(::Parameter{T}) where {T} = constant(zero(T))
norm(p::Parameter) = p.value
real(p::Parameter) = p.value
imag(p::Parameter) = 0.0

convert(::Type{Parameter{T}}, v::T) where {T<:Number} = constant(v)
convert(::Type{Parameter{T}}, v::S) where {T<:Number,S<:Number} = constant(T(v))
function convert(::Type{Parameter{T}}, v::Parameter{S}) where
    {T<:Number,S<:Number}
    typ = promote_type(T, S)
    Parameter{promote_type(T,S)}(typ(initial(v)), typ(lower(v)), typ(upper(v)))
end
promote_rule(::Type{Parameter{T}}, ::Type{S}) where
    {T<:Number,S<:Number} = Parameter{promote_type(T,S)}
promote_rule(::Type{T}, ::Type{Parameter{S}}) where
    {T<:Number,S<:Number} = Parameter{promote_type(T,S)}
promote_rule(::Type{Parameter{T}}, ::Type{Parameter{S}}) where
    {T<:Number,S<:Number} = Parameter{promote_type(T,S)}
promote_type(::Type{Parameter{T}}, ::Type{S}) where
    {T<:Number,S<:Number} = Parameter{promote_type(T,S)}
promote_type(::Type{Parameter{T}}, ::Type{Parameter{S}}) where
    {T<:Number,S<:Number} = Parameter{promote_type(T,S)}

Base.show(io::IO, prm::Parameter) = prm.isfixed ? print(io, prm.value) :
    print(io, prm.lower, " ≤ ", prm.value, " ≤ ", prm.upper)

(+)(lhs::Parameter, rhs::Parameter) = lhs.value + rhs.value
(-)(lhs::Parameter, rhs::Parameter) = lhs.value - rhs.value
(*)(lhs::Parameter, rhs::Parameter) = lhs.value * rhs.value
(/)(lhs::Parameter, rhs::Parameter) = lhs.value / rhs.value

(+)(lhs::Parameter, rhs::Number) = lhs.value + rhs
(+)(lhs::Number, rhs::Parameter) = lhs + rhs.value
(-)(lhs::Parameter, rhs::Number) = lhs.value - rhs
(-)(lhs::Number, rhs::Parameter) = lhs - rhs.value
(*)(lhs::Parameter, rhs::Number) = lhs.value * rhs
(*)(lhs::Number, rhs::Parameter) = lhs * rhs.value
(/)(lhs::Parameter, rhs::Number) = lhs.value / rhs
(/)(lhs::Number, rhs::Parameter) = lhs / rhs.value

isequal(lhs::Parameter, rhs::Parameter) = isequal(lhs.value, rhs.value)
isequal(lhs::Number, rhs::Parameter) = isequal(lhs, rhs.value)
isequal(lhs::Parameter, rhs::Number) = isequal(lhs.value, rhs)
==(lhs::Parameter, rhs::Parameter) = ==(lhs.value, rhs.value)
==(lhs::Number, rhs::Parameter) = ==(lhs, rhs.value)
==(lhs::Parameter, rhs::Number) = ==(lhs.value, rhs)
(<)(lhs::Parameter, rhs::Parameter) = lhs.value < rhs.value
(<)(lhs::Number, rhs::Parameter) = lhs < rhs.value
(<)(lhs::Parameter, rhs::Number) = lhs.value < rhs
isapprox(lhs::Parameter, rhs::Parameter) = isapprox(lhs.value, rhs.value)
isapprox(lhs::Number, rhs::Parameter) = isapprox(lhs, rhs.value)
isapprox(lhs::Parameter, rhs::Number) = isapprox(lhs.value, rhs)

const fun = (
    :abs2, :acosh, :acot, :acotd, :acoth, :acsc, :acscd, :acsch, :asec,
    :asecd, :asinh, :atan, :atand, :cbrt, :cos, :cosd, :cosh, :cot, :cotd,
    :coth, :csc, :cscd, :csch, :exp, :exp2, :expm1, :gamma, :inv, :lgamma,
    :log, :log10, :log1p, :log2, :sec, :secd, :sech, :sin, :sind, :sinh,
    :sqrt, :tan, :tand, :tanh, :acos, :acosd, :asech, :asin, :asind, :atanh,
)
for f in fun
    @eval begin
        import Base.$f
        export $f
        $f(par::Parameter) = $f(par.value)
    end
end

getparameters(data) = getparameters(Parameter[], data)
function getparameters(params, prm::Parameter)
    isparameter(prm) && push!(params, prm)
    return params
end
function getparameters(params, arr::AbstractArray)
    foreach(x -> getparameters(params, x), arr)
    return params
end
function getparameters(params, val)
    fields = fieldnames(val)
    foreach(x -> getparameters(params, getfield(val, x)), fields)
    return params
end

include("macros.jl")

end # module


text:
# Compose pseudo-forms for simple symbols, all parameterized by center and size

module Shape

using Measures
using Compose: x_measure, y_measure, circle, rectangle, polygon, line, px

shape_docstr(kinds) = "Draw $kinds at the coordinates specified in `xs` and
`ys` of size `rs`"

"""
    square(xs, ys, rs)

$(shape_docstr("squares"))
"""
function square(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
    n = max(length(xs), length(ys), length(rs))

    s = 1/sqrt(2)
    polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)

    for i in 1:n
        x = x_measure(xs[mod1(i, length(xs))])
        y = y_measure(ys[mod1(i, length(ys))])
        r = rs[mod1(i, length(rs))]
        polys[i] = Tuple{Measure,Measure}[(x-s*r,y-s*r),(x-s*r,y+s*r),(x+s*r,y+s*r),(x+s*r,y-s*r)]
    end

    return polygon(polys)
end

"""
    diamond(xs, ys, rs)

$(shape_docstr("diamonds"))
"""
function diamond(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
    n = max(length(xs), length(ys), length(rs))

    polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
    for i in 1:n
        x = x_measure(xs[1 + i % length(xs)])
        y = y_measure(ys[1 + i % length(ys)])
        r = rs[1 + i % length(rs)]
        polys[i] = Tuple{Measure, Measure}[(x, y - r), (x + r, y), (x, y + r), (x - r, y)]
    end

    return polygon(polys)
end

"""
    cross(xs, ys, rs)

$(shape_docstr("crosses"))
"""
function cross(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    u = 0.4r

    # make a "plus sign"
    polys[i] = Tuple{Measure, Measure}[
      (x-r, y-u), (x-r, y+u), # L edge
      (x-u, y+u),             # BL inside
      (x-u, y+r), (x+u, y+r), # B edge
      (x+u, y+u),             # BR inside
      (x+r, y+u), (x+r, y-u), # R edge
      (x+u, y-u),             # TR inside
      (x+u, y-r), (x-u, y-r), # T edge
      (x-u, y-u)              # TL inside
    ]
  end

  return polygon(polys)
end

"""
    xcross(xs, ys, rs)

$(shape_docstr("rotated crosses"))
"""
function xcross(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
  s = 1/sqrt(5)
  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    u = s*r
    polys[i] = Tuple{Measure, Measure}[
      (x, y - u), (x + u, y - 2u), (x + 2u, y - u),
      (x + u, y), (x + 2u, y + u), (x + u, y + 2u),
      (x, y + u), (x - u, y + 2u), (x - 2u, y + u),
      (x - u, y), (x - 2u, y - u), (x - u, y - 2u)
    ]
  end

  return polygon(polys)
end

"""
    utriangle(xs, ys, rs, scalar=1)

$(shape_docstr("upward-pointing triangles"))
"""
function utriangle(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray, scalar = 1)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)

  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    u = 0.87 * r
    v = scalar * r
    polys[i] = Tuple{Measure, Measure}[
      (x - u, y + 0.5v),
      (x + u, y + 0.5v),
      (x, y - v),  # two of these so the centroid is such that it zooms well
      (x, y - v)
    ]
  end

  return polygon(polys)
end

"""
    dtriangle(xs, ys, rs)

$(shape_docstr("downward-pointing triangles"))
"""
dtriangle(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray) = utriangle(xs, ys, rs, -1)

"""
    star1(xs, ys, rs, scalar=1)

$(shape_docstr("five-pointed stars"))
"""
function star1(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray, scalar = 1)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)

  # some magic scalars
  sx = 0.7
  sy1, sy2 = 1.2, 0.4

  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    polys[i] = Tuple{Measure, Measure}[
      (x-sx*r, y+r),      # BL
      (x,      y-sy1*r),  # T
      (x+sx*r, y+r),      # BR
      (x-r,    y-sy2*r),  # L
      (x+r,    y-sy2*r)   # R
    ]
  end

  return polygon(polys)
end

"""
    star2(xs, ys, rs, scalar=1)

$(shape_docstr("four-pointed stars"))
"""
function star2(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray, scalar = 1)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    u = 0.4r
    polys[i] = Tuple{Measure, Measure}[
      (x-u, y),   (x-r, y-r), # TL
      (x,   y-u), (x+r, y-r), # TR
      (x+u, y),   (x+r, y+r), # BR
      (x,   y+u), (x-r, y+r)  # BL
    ]
  end

  return polygon(polys)
end

"""
    hexagon(xs, ys, rs)

$(shape_docstr("hexagons"))
"""
function hexagon(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))

  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    v = 0.5r
    u = 0.87r
    polys[i] = Tuple{Measure, Measure}[
      (x-u, y-v), (x-u, y+v), # L edge
      (x, y+r),               # B
      (x+u, y+v), (x+u, y-v), # R edge
      (x, y-r)                # T
    ]
  end

  return polygon(polys)
end

"""
    octagon(xs, ys, rs)

$(shape_docstr("octagons"))
"""
function octagon(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))

  polys = Vector{Vector{Tuple{Measure, Measure}}}(undef, n)
  for i in 1:n
    x = x_measure(xs[mod1(i, length(xs))])
    y = y_measure(ys[mod1(i, length(ys))])
    r = rs[mod1(i, length(rs))]
    u = 0.4r

    polys[i] = Tuple{Measure, Measure}[
      (x-r, y-u), (x-r, y+u), # L edge
      (x-u, y+r), (x+u, y+r), # B edge
      (x+r, y+u), (x+r, y-u), # R edge
      (x+u, y-r), (x-u, y-r)  # T edge
    ]
  end

  return polygon(polys)
end

"""
    hline(xs, ys, rs)

$(shape_docstr("horizontal lines"))
"""
function hline(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
    n = max(length(xs), length(ys), length(rs))

    line_ps = Vector{Vector{Tuple{Measure,Measure}}}(undef, n)
    for i in 1:n
        x = x_measure(xs[1 + i % length(xs)])
        y = y_measure(ys[1 + i % length(ys)])
        r = rs[1 + i % length(rs)]
        line_ps[i] = Tuple{Measure, Measure}[(x-r,y-1.5px),(x+r,y-1.5px),(x+r,y+1.5px),(x-r,y+1.5px)]
    end

    return polygon(line_ps)
end

"""
    vline(xs, ys, rs)

$(shape_docstr("vertical lines"))
"""
function vline(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
    n = max(length(xs), length(ys), length(rs))

    line_ps = Vector{Vector{Tuple{Measure,Measure}}}(undef, n)
    for i in 1:n
        x = x_measure(xs[1 + i % length(xs)])
        y = y_measure(ys[1 + i % length(ys)])
        r = rs[1 + i % length(rs)]
        line_ps[i] = Tuple{Measure, Measure}[(x-1.5px,y-r),(x-1.5px,y+r),(x+1.5px,y+r),(x+1.5px,y-r)]
    end

    return polygon(line_ps)
end

end  # module Shape


text:
function formreturns(panel::DataFrame, ms::MeasureSpec;
  returnsname::String = PARAM[:returnsname],
  savepath::String = PARAM[:workingpath],
  returntype::String = PARAM[:crspdatalabel],
  frequencysuffix::String = PARAM[:crspfrequencysuffix],
  Fvol = PARAM[:Fvol],
  volscale::Float64 = PARAM[:volscale]
  #testdate::Date = PARAM[:testdate],
  )::DataFrame

  #housekeeping
  local returnspath::String = "$savepath\\$(returntype)_$(returnsname).csv"
  local returns::DataFrame

  #The below block directs the panel creation
  if PARAM[:refreshreturns]
    #returns = DataFrame(date=sort!(unique(panel.date)))
    returns = vcat(
      DataFrame(date=Date[], Wvolume=MFloat64[], volume=MFloat64[]),
      combine(groupby(panel, :date),
        :date => last => :date,
        :Wdvol => ((v)->sum(skipmissing(v))) => :Wvolume,
        :dvol => ((v)->sum(skipmissing(v)) ./ volscale) => :volume,
        ))

    #pad the dates to get an extra data point of returns
    if frequencysuffix ≡ "m"
      push!(returns, (; date=minimum(returns.date) - Month(1), Wvolume=missing, volume=missing))
    elseif frequencysuffix ≡ "w"
      push!(returns, (; date=minimum(returns.date) - Day(7), Wvolume=missing, volume=missing))
    else
      throw("unrecognized date frequency")
    end

    sret = returns[completecases(returns, [:volume, :Wvolume]), :]
    println("initial: size: $(size(sret)) dates: $(minimum(sret.date)):$(maximum(sret.date))")
    sort!(returns, :date)

    returns = addffreturns!(returns)
    sret = returns[completecases(returns, [:volume, :Wvolume]), :]
    println("afterff: size: $(size(sret)) dates: $(minimum(sret.date)):$(maximum(sret.date))")
    sort!(returns, :date)

    returns = addstrategyreturns!(returns, panel, ms)
    sret = returns[completecases(returns, [:volume, :Wvolume]), :]
    println("afterstrategy: size: $(size(sret)) dates: $(minimum(sret.date)):$(maximum(sret.date))")
    sort!(returns, :date)

    returns |> CSV.write(returnspath)

    sum(completecases(returns)) == nrow(returns) || @warn "WARNING: Some dates in the returns file
      are missing benchmark returns. Make sure you know why before disabling this error- in
      particular, it could be a symptom of stale benchmark data. You can check the returns file
      to see what is missing. On the other hand, if its just a result of using a lag
      then its probably safe to ignore this warning."
  else
    returns = CSV.File(returnspath) |> DataFrame
  end

  return returns
end

getreturns() = CSV.File(
    "$(PARAM[:workingpath])\\$(PARAM[:crspdatalabel])_$(PARAM[:returnsname]).csv") |> DataFrame

#compute the strategy returns
function addstrategyreturns!(returns::DataFrame, panel::DataFrame, ms::MeasureSpec)

  @assert issorted(returns, :date)

  for ξ ∈ ms.ξs
    #acquire the weight field
    #NOTE- not sure whether to use Lw or w here- doesn't mkae a big difference
    Fw = ξ.X[:FLw]
    strategyreturns = combine(groupby(view(panel, completecases(panel, [Fw]),:), :date),
      :date => last => :date, #just record the date
      [Fw, :ret] => weightedreturn => ξ.Fξ,
      [Fw, :ret] => ((w,ret)->weightedreturn(w .|> (x)->max(x,0.0),ret)) => Symbol(ξ.Fξ,:_long),
      [Fw, :ret] => ((w,ret)->weightedreturn(w .|> (x)->min(x,0.0),ret)) => Symbol(ξ.Fξ,:_short),
      )

    #=the below code reverses the sign of the return to the same as the effect on leverage
    #note this is not the same as passive change in portfolio leverage- rather, its the
    #passive change in portfolio leverage caused by net gains and losses
    Basically, profits from the long side are leveraging, while profits from short side
    #are deleveraging
    Logic to create the measure is as follows: start with the leverage effect
      (long profits - short profits)
      Then if part of the leverage effect is due to comovement (e.g. one side
      created gains, the other losses) then subtract out the offsetting part of the positions
      as that is most likely the effect of the market. If both sides created the same
      sign P&L, then the leverage effect is entirely due to the P&L so don't subtract anything.
      (this is actually the logic of the check)
      =#

    Flong, Fshort, Fnet = Symbol(ξ.Fξ,:_long), Symbol(ξ.Fξ,:_short), ξ.Fξ
    transform!(strategyreturns,
      [Flong, Fshort] => ((long, short) ->  long .- short) => Symbol(ξ.Fξ,:_retlevgross),
      [Flong, Fshort, Fnet] => ByRow((long, short, net) ->
        ifelse(sign(long)==sign(short),long - short, ifelse(abs(long)>abs(short),net, -net))) =>
        Symbol(ξ.Fξ,:_retlev),
      [Flong, Fshort, Fnet] => ByRow((long, short, net) -> long - short -
        ifelse(sign(long)==sign(short),0.0, 2.0*sign(long - short)*min(abs(long),abs(short)))) =>
        Symbol(ξ.Fξ,:_retlevcheck))
    @assert strategyreturns[!, Symbol(ξ.Fξ,:_retlev)] ≈ strategyreturns[!, Symbol(ξ.Fξ,:_retlevcheck)] "
      zip(strategyreturns[!, Symbol(ξ.Fξ,:_retlev)],strategyreturns[!, Symbol(ξ.Fξ,:_retlevcheck)])
        $(zip(strategyreturns[!, Symbol(ξ.Fξ,:_retlev)],
          strategyreturns[!, Symbol(ξ.Fξ,:_retlevcheck)]) |> collect)"
    select!(strategyreturns, Not(Symbol(ξ.Fξ,:_retlevcheck)))

    #the below is another choice- using a leftjoin uses all available data, but
    #using an innerjoin keeps the time windows consistent
    #the below is a compromise- use a left join but bound how far it goes back
    startdate = returns.date[
      max(something(findfirst(isequal(minimum(strategyreturns.date)), returns.date),1)-1,1)]
    boundedreturns = view(returns,
      (returns.date .≥ startdate) .&
      (returns.date .≤ maximum(strategyreturns.date)), :)
    returns=leftjoin(boundedreturns, strategyreturns, on=:date)

    if !issorted(returns, :date)
      sort!(returns, :date)
    end
    #below is the old innerjoin code
    #=returns = innerjoin(returns, combine(groupby(strategyreturns, :date),
      :date => last => :date, #just record the date
      [Fw, :ret] => weightedreturn => ξ.Fξ), on=:date)=#
  end

  return returns
end

#note this natually will lead to higher/lower returns depending on the gross leverage employed
#e.g. if the long leg returns +10%, short leg returns -10%, and both sides have weights
#that add to 1.0/-1.0, then the resulting return will be 1.0*0.1 - 1.0*-0.1 = 20%
weightedreturn(weights, returns) = sum(weights .* returns)

#function addstrategy

#merges in the FF returns
function addffreturns!(returns::DataFrame)
  #this is less important when just combining the time series as opposed to creating a panel
  #but we might as well be consistent
  local datefrequency
  if PARAM[:crspfrequencysuffix] == "m"
    datefrequency = :month
  elseif PARAM[:crspfrequencysuffix] == "w"
    datefrequency = :week
  else
    @assert false
  end

  #do I break up mergeff or bootstrap a group into returns?
  returns = mergeff(returns; datefrequency, testlabel="FFreturns")

  return returns
end

#################################################
#the rest of these methods handle processing wide-form returns
#and merging the returns with a panel

#computes a cumulative return skipping missing values in place
function skippingcumret(rets::AbstractVector{T}) where T<:MFloat64
  base::Float64 = 1.0

  out::Vector{T} = rets .+ 1.0
  completeout = view(out,out .!== missing) #want to skip missing values while totalling
  completeout .= cumprod(completeout)

  return out
end

#creates a cumulative index from a return stream
function returns2index(dates::AbstractVector{Date},
  rets::AbstractVector{<:MFloat64})

  issorted(dates) || error("dates must be sorted to form index")

  @assert length(dates) == length(rets)
  index::Vector{Float64} = skippingcumret(rets)

  return :index, DataFrame(date=dates, index=index)
end

#creates a cumulative index from a return stream
function returns2index(df::AbstractDataFrame;
  Fdate::Symbol = :date,
  Frets::Vector{Symbol} = error("Frets is required"))

  issorted(df, Fdate) || error("df must be sorted by date")
  #create the cumulative index
  indexdf = DataFrame(date=df.date)
  Findexes = (Fret->Symbol(Fret,:_idx)).(Frets)
  for (Fret, Findex) ∈ zip(Frets, Findexes)
    indexdf[!, Findex] = skippingcumret(df[!, Fret])
  end
  @assert issetequal(propertynames(indexdf), [:date; Findexes])

  #WARNING the compelete cases is important here
  indexdf = indexdf[completecases(indexdf[!, Findexes]),:]

  return Findexes,indexdf
end

#merge in returns from a file of indexes
#1) merge indexdf with unique dates
#2) generate returns
#3) merge the returns with the provided df

#WARNING WARNING WARNING the join here does not necessitate a complete match between
#the df and the index- callers responsibility to check index dates are comprehensive
#WARNING WARNING WARNING order of rows is undefined
function index2returns!(df::DataFrame, indexdf::AbstractDataFrame;
  Findexes::Vector{Symbol} = throw("Findexes is required"),
  Frets::Vector{Symbol} = throw("Frets is required"),
  datefrequency::Symbol = throw("date frequency is required (:day, :month, :week)"),
  testlabel::String = throw("testlabel is required"),
  Fgroup::NSymbol = nothing, #only used for testing purposes if provided
  validateexactfrequency = false)

  if Fgroup === nothing
    allunique(df.date) || throw("Fgroup must be provided for testing purposes
      if dates in df are not unique")
  end

  @assert allunique(indexdf.date)
  @assert issorted(indexdf, :date)
  @assert all(indexdf.date .!== missing)
  @assert all(df.date .!== missing)
  if :indexdate ∈ propertynames(df)
    throw("Column indexdate should not be in df. Call select!(bm, Not(:indexdate))
      to delete the column prior to calling index2returns")
  end

  #join on the dates, treating months as a special case
  bm = DataFrame(date=sort!(unique(df.date)))
  if datefrequency === :month
    indexdf.dateym = (d->year(d) + month(d) / 100).(indexdf.date)
    indexdf = combine(last, groupby(indexdf, :dateym))
    bm.dateym = (d->year(d) + month(d) / 100).(bm.date)

    #avoid duplicate columns in the join
    rename!(indexdf, :date=>:indexdate)

    #throwing the error here from funds-mf:Merge key(s) in df1 are not unique
    bm = leftjoin(bm, indexdf, on=:dateym, validate=(true, true))
    sort!(bm, [:dateym])

    #the following creates a crosswalk for testing and validation
    indexdf = leftjoin(indexdf, DataFrame(date=bm.date, dateym = bm.dateym), on=:dateym)
    sort!(indexdf, :dateym)
  elseif datefrequency ∈ [:day, :week]
    bm = leftjoin(bm, indexdf, on=:date, validate=(true, true))
    sort!(bm, :date)
  elseif datefrequency === :year
    throw("year date frequency not supported (see month for an idea how to do this)")
  else
    throw("unrecognized date frequency $datefrequency")
  end

  #bm |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_bmlab1.csv")
  #@info "Wrote bm1. datefrequency=$datefrequency"
  #check for missing values
  #=for Findex ∈ Findexes
    indicesmissing = sum(bm[!, Findex] .≡ missing)
    if indicesmissing > 0
      @warn ("BM dates found with no corresponding index. Testlabel: $testlabel, " *
        "Field: $Findex, nummissing=$indicesmissing
        Check testfiles ending with $(testlabel)_ind2ret and $(testlabel)_indmerged to" *
        "get some details.")
    end
  end=#

  #another check on the above join check on the above join
  validatedates(bm.date |> Vector{Date}, frequency=datefrequency; validateexactfrequency)

  #validate dates by comparing the sets of the benchmark and the df,
  #ignoring mismatches beyond the endpoints
  indexdates = indexdf[completecases(indexdf, [Findexes; :date]), :date]
  dfnotindex = setdiff(bm.date, indexdates)
  dfnotindex = (!isempty(dfnotindex))  ? dfnotindex[(minimum(indexdates) .≤ dfnotindex) .&
    (dfnotindex .≤ maximum(indexdates))] : dfnotindex

  #=indexnotdf = setdiff(indexdates, bm.date)
  indexnotdf = (!isempty(indexnotdf)) ? indexnotdf[(minimum(bm.date) .≤ indexnotdf) .&
    (indexnotdf .≤ maximum(bm.date))] : indexnotdf=#

  (!isempty(dfnotindex)) && (@warn("df dates not in index (interior only): $dfnotindex" ))
  #(!isempty(indexnotdf)) && (@warn("index dates not in df (interior only): $indexnotdf") )

  for (Fret, Findex) ∈ zip(Frets, Findexes)
    bm[!, Fret] = [missing; bm[2:end, Findex] ./ bm[1:(end-1), Findex] .- 1]
  end

  #join the index with the panel data
  df = leftjoin(df, bm, on=:date, validate=(false, true))

  #write output for testing purposes
  if PARAM[:testindex2returns]
    bm |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_$(testlabel)_ind2ret.csv")
    local testview::SubDataFrame
    if Fgroup ≡ nothing
      testview = view(df, :, :)
    elseif Fgroup ≡ :permno
      testpermno::Int = PARAM[:testpermnomult]
      testview = view(df, (df[!, Fgroup] .% testpermno) .== 0, :)
    elseif Fgroup ≡ :fundid
      testids = values(PARAM[:fundbetatestfundid])
      testview = view(df, (g->g ∈ testids).(df[!, Fgroup]), :)
    end
    testview |> CSV.write(
      "$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_$(testlabel)_indmerged.csv")
  end

  if :dateym ∈ propertynames(df)
    select!(df, Not([:dateym]))
  end
  select!(df, Not(Findexes))
  return df
end

#merge in returns from the returns file
#core of code is in index2returns- just calls a generic function to generate the indexdf
#WARNING WARNING WARNING the joins here do not necessitate a complete match between
#the df and the index- callers responsibility to check index dates are comprehensive
#WARNING WARNING WARNING order of rows is undefined
function mergereturns(df::DataFrame, returns::AbstractDataFrame;
  Frets::Vector{Symbol} = throw("Frets is required"),
  datefrequency::Symbol = throw("date frequency is required (:day, :month, :week)"),
  testlabel::String = throw("testlabel is required"),
  Fgroup::NSymbol = nothing, #only used for testing purposes if provided
  validateexactfrequency = false)

  #form an index for each return stream
  Findexes, indexdf = returns2index(returns; Frets)

  #returns |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_returnstemp.csv")
  #indexdf |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_indextemp.csv")

  #WARNING WARNING WARNING order of rows is undefined
  return index2returns!(df, indexdf;
    Frets, Findexes, datefrequency, testlabel, Fgroup, validateexactfrequency)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Pkg; Pkg.activate("../")

using StatsPlots
using Turing
using Distributions
using StatsFuns
import Random

@model GMM(y, K) = begin
    mu ~ filldist(Normal(0, .3), K)
    sigmasq ~ filldist(InverseGamma(3, 2), K)
    sigma = sqrt.(sigmasq)
    eta ~ Dirichlet(K, 1/K)
    y .~ UnivariateGMM(mu, sigma, Categorical(eta))

    return (mu=mu, sigma=sigma, eta=eta)
end;

Random.seed!(1);
y = rand(UnivariateGMM([-2, 2], [0.5, 0.8], Categorical([.3, .7])), 1000)
histogram(y, bins=50)
m = GMM(y, 5)
chain = sample(m, NUTS(500, 0.65), 1000)
return_values = generated_quantities(m, chain)

mu = hcat(getfield.(return_values, :mu)...)
sigma = hcat(getfield.(return_values, :sigma)...)
eta = hcat(getfield.(return_values, :eta)...)

plot(mu')
mean(mu, dims=2)

yy = range(-6, 6, length=200)
ypred = [pdf.(UnivariateGMM(r[:mu], r[:sigma], Categorical(r[:eta])), yy)
         for r in return_values]
yp = hcat(ypred...)
plot(yy, mean(yp, dims=2))
histogram!(y, alpha=0.3, normalize=true, bins=50)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# Autogenerated wrapper script for PDFHighlights_jll for x86_64-w64-mingw32
export PDFHighlightsWrapper

using Poppler_jll
JLLWrappers.@generate_wrapper_header("PDFHighlights")
JLLWrappers.@declare_library_product(PDFHighlightsWrapper, "PDFHighlightsWrapper.dll")
function __init__()
    JLLWrappers.@generate_init_header(Poppler_jll)
    JLLWrappers.@init_library_product(
        PDFHighlightsWrapper,
        "bin\\PDFHighlightsWrapper.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Base.position

export AbstractSystem
export BoundaryCondition, DirichletZero, Periodic, infinite_box
export bounding_box, boundary_conditions, periodicity, n_dimensions, species_type
export position, velocity, element, atomic_mass, atomic_number, atomic_symbol

#
# Identifier for boundary conditions per dimension
#
abstract type BoundaryCondition end
struct DirichletZero <: BoundaryCondition end  # Dirichlet zero boundary (i.e. molecular context)
struct Periodic <: BoundaryCondition end  # Periodic BCs

infinite_box(::Val{1}) = [[Inf]]u"bohr"
infinite_box(::Val{2}) = [[Inf, 0], [0, Inf]]u"bohr"
infinite_box(::Val{3}) = [[Inf, 0, 0], [0, Inf, 0], [0, 0, Inf]]u"bohr"
infinite_box(dim::Int) = infinite_box(Val(dim))

#
# Abstract system
#
"""
    AbstractSystem{D}

A `D`-dimensional system.
"""
abstract type AbstractSystem{D} end

"""
    bounding_box(sys::AbstractSystem{D})

Return a vector of length `D` of vectors of length `D` that describe the "box" in which the system `sys` is defined.
"""
function bounding_box end

"""
    boundary_conditions(sys::AbstractSystem{D})

Return a vector of length `D` of `BoundaryCondition` objects, one for each direction described by `bounding_box(sys)`.
"""
function boundary_conditions end

"""
    species_type(::AbstractSystem)

Return the type used to represent a species or atom.
"""
function species_type end

"""Return vector indicating whether the system is periodic along a dimension."""
periodicity(sys::AbstractSystem) = [isa(bc, Periodic) for bc in boundary_conditions(sys)]

"""
    n_dimensions(::AbstractSystem)
    n_dimensions(atom)

Return number of dimensions.
"""
n_dimensions(::AbstractSystem{D}) where {D} = D
# Note: Can't use ndims, because that is ndims(sys) == 1 (because of indexing interface)

# indexing and iteration interface...need to implement getindex and length, here are default dispatches for others
Base.size(s::AbstractSystem) = (length(s),)
Base.firstindex(::AbstractSystem) = 1
Base.lastindex(s::AbstractSystem) = length(s)
# default to 1D indexing
Base.iterate(sys::AbstractSystem, state = firstindex(sys)) =
    (firstindex(sys) <= state <= length(sys)) ? (@inbounds sys[state], state + 1) : nothing

# TODO Support similar, push, ...

#
# Species property accessors from systems and species
#
"""The element corresponding to a species/atom (or missing)."""
function element end

"""
    position(sys::AbstractSystem{D})
    position(sys::AbstractSystem, index)
    position(species)

Return a vector of positions of every particle in the system `sys`. Return type
should be a vector of vectors each containing `D` elements that are
`<:Unitful.Length`. If an index is passed or the action is on a `species`,
return only the position of the referenced `species` / species on that index.
"""
position(sys::AbstractSystem)        = position.(sys)    # in Cartesian coordinates!
position(sys::AbstractSystem, index) = position(sys[index])

"""
    velocity(sys::AbstractSystem{D})
    velocity(sys::AbstractSystem, index)
    velocity(species)

Return a vector of velocities of every particle in the system `sys`. Return
type should be a vector of vectors each containing `D` elements that are
`<:Unitful.Velocity`. If an index is passed or the action is on a `species`,
return only the velocity of the referenced `species`. Returned value of the function
may be `missing`.
"""
velocity(sys::AbstractSystem)        = velocity.(sys)    # in Cartesian coordinates!
velocity(sys::AbstractSystem, index) = velocity(sys[index])

"""
    atomic_mass(sys::AbstractSystem)
    atomic_mass(sys::AbstractSystem, i)
    atomic_mass(species)

Vector of atomic masses in the system `sys` or the atomic mass of a particular `species` /
the `i`th species in `sys`. The elements are `<: Unitful.Mass`.
"""
atomic_mass(sys::AbstractSystem)        = atomic_mass.(sys)
atomic_mass(sys::AbstractSystem, index) = atomic_mass(sys[index])

"""
    atomic_symbol(sys::AbstractSystem)
    atomic_symbol(sys::AbstractSystem, i)
    atomic_symbol(species)

Vector of atomic symbols in the system `sys` or the atomic symbol of a particular `species` /
the `i`th species in `sys`.
"""
atomic_symbol(sys::AbstractSystem)        = atomic_symbol.(sys)
atomic_symbol(sys::AbstractSystem, index) = atomic_symbol(sys[index])

"""
    atomic_number(sys::AbstractSystem)
    atomic_number(sys::AbstractSystem, i)
    atomic_number(species)

Vector of atomic numbers in the system `sys` or the atomic number of a particular `species` /
the `i`th species in `sys`.
"""
atomic_number(sys::AbstractSystem)        = atomic_number.(sys)
atomic_number(sys::AbstractSystem, index) = atomic_number(sys[index])


text:
module TestScheduleUndirectedWiringDiagrams
using Test

using LinearAlgebra: dot, tr

using Catlab.Theories: codom
using Catlab.CategoricalAlgebra.CSets
using Catlab.WiringDiagrams
using Catlab.WiringDiagrams.ScheduleUndirectedWiringDiagrams:
  composites, composite_junction, composite_ports, parent, box_parent
using Catlab.LinearAlgebra.TensorNetworks: @tensor_network

composite_junctions(x::ACSet, c) =
  composite_junction.(Ref(x), composite_ports(x, c))

# Open linear path
##################

d = @tensor_network out[v,z] = A[v,w] * B[w,x] * C[x,y] * D[y,z]
s = schedule(d, order=reverse(boxes(d)))
@test length(composites(s)) == 3
@test parent(s) == [2,3,3]
@test box_parent(s) == [3,2,1,1]

nd = to_nested_diagram(s)
@test parent(nd) == [2,3,3]
@test box_parent(nd) == [3,2,1,1]
@test composite_junctions(nd, 1:3) == [[3,5], [2,5], [1,5]]

matrices = map(randn, [(3,4), (4,5), (5,6), (6,7)])
out = eval_schedule(nd, matrices)
@test out ≈ foldr(*, matrices)

# Closed cycle
##############

d = @tensor_network out[] = A[w,x] * B[x,y] * C[y,z] * D[z,w]
s = schedule(d)
nd = to_nested_diagram(s)
@test parent(nd) == [2,3,3]
@test box_parent(nd) == [1,1,2,3]
@test map(length, composite_ports(nd, 1:3)) == [2,2,0]
@test composite_junctions(nd, 1:2) == [[1,3], [1,4]]

matrices = map(randn, [(10,5), (5,5), (5,5), (5,10)])
out = eval_schedule(nd, matrices)
@test out[] ≈ tr(foldl(*, matrices))

# Tensor product
################

d = @tensor_network out[w,x,y,z] = A[w,x] * B[y,z]
s = schedule(d)
@test parent(s) == [1]
@test box_parent(s) == [1,1]

A, B = randn((3,4)), randn((5,6))
out = eval_schedule(s, [A, B])
@test out ≈ (reshape(A, (3,4,1,1)) .* reshape(B, (1,1,5,6)))

# Frobenius inner product
#########################

d = @tensor_network out[] = A[x,y] * B[x,y]
s = schedule(d)
@test parent(s) == [1]
@test box_parent(s) == [1,1]

A, B = randn((5,5)), randn((5,5))
out = eval_schedule(s, [A, B])
@test out[] ≈ dot(vec(A), vec(B))

end


text:
export BrickletCO2CO2ConcentrationCallbackThreshold
struct BrickletCO2CO2ConcentrationCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletCO2Identity
struct BrickletCO2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletCO2
"""
Measures CO2 concentration in ppm
"""
mutable struct BrickletCO2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletCO2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_co2")
        deviceInternal = package.BrickletCO2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_co2_concentration
"""
    $(SIGNATURES)

Returns the measured CO2 concentration.

If you want to get the CO2 concentration periodically, it is recommended to use
the :cb:`CO2 Concentration` callback and set the period with
:func:`Set CO2 Concentration Callback Period`.
"""
function get_co2_concentration(device::BrickletCO2)
    return device.deviceInternal.get_co2_concentration()
end

export set_co2_concentration_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`CO2 Concentration` callback is
triggered periodically. A value of 0 turns the callback off.

The :cb:`CO2 Concentration` callback is only triggered if the CO2 concentration
has changed since the last triggering.
"""
function set_co2_concentration_callback_period(device::BrickletCO2, period)
    device.deviceInternal.set_co2_concentration_callback_period(period)
end

export get_co2_concentration_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set CO2 Concentration Callback Period`.
"""
function get_co2_concentration_callback_period(device::BrickletCO2)
    return device.deviceInternal.get_co2_concentration_callback_period()
end

export set_co2_concentration_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`CO2 Concentration Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the CO2 concentration is *outside* the min and max values"
 "'i'",    "Callback is triggered when the CO2 concentration is *inside* the min and max values"
 "'<'",    "Callback is triggered when the CO2 concentration is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the CO2 concentration is greater than the min value (max is ignored)"
"""
function set_co2_concentration_callback_threshold(device::BrickletCO2, option, min, max)
    device.deviceInternal.set_co2_concentration_callback_threshold(option, min, max)
end

export get_co2_concentration_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set CO2 Concentration Callback Threshold`.
"""
function get_co2_concentration_callback_threshold(device::BrickletCO2)
    return device.deviceInternal.get_co2_concentration_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`CO2 Concentration Reached`,

are triggered, if the thresholds

* :func:`Set CO2 Concentration Callback Threshold`,

keep being reached.
"""
function set_debounce_period(device::BrickletCO2, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletCO2)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletCO2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletCO2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
Bayesian Poisson Mixture Model
"""
module PoissonMixtureModel
using StatsFuns.logsumexp, SpecialFunctions.digamma
using Distributions

export Gam, BPMM, Poi, PMM
export sample_PMM, sample_data, winner_takes_all
export learn_GS, learn_CGS, learn_VI

####################
## Types
struct Gam
    # Parameters of Gamma distribution
    a::Vector{Float64}
    b::Float64
end

struct BPMM
    # Parameters of Bayesian Poisson Mixture Model
    D::Int
    K::Int
    alpha::Vector{Float64}
    cmp::Vector{Gam}
end

struct Poi
    # Parameters of Poisson Distribution
    lambda::Vector{Float64}
end

struct PMM
    # Parameters of Poisson Mixture Model
    D::Int
    K::Int
    phi::Vector{Float64}
    cmp::Vector{Poi}
end

####################
## Common functions
"""
Sample a PMM given hyperparameters.
"""
function sample_PMM(bpmm::BPMM)
    cmp = Vector{Poi}()
    for c in bpmm.cmp
        lambda = Vector{Float64}()
        for d in 1 : bpmm.D
            push!(lambda, rand(Gamma(c.a[d], 1.0/c.b)))
        end
        push!(cmp, Poi(lambda))
    end
    phi = rand(Dirichlet(bpmm.alpha))
    return PMM(bpmm.D, bpmm.K, phi, cmp)
end

"""
Sample data from a specific PMM model.
"""
function sample_data(pmm::PMM, N::Int)
    X = zeros(pmm.D, N)
    S = categorical_sample(pmm.phi, N)
    for n in 1 : N
        k = indmax(S[:, n])
        for d in 1 : pmm.D
            X[d,n] = rand(Poisson(pmm.cmp[k].lambda[d]))
        end
    end
    return X, S
end

categorical_sample(p::Vector{Float64}) = categorical_sample(p, 1)[:,1]
function categorical_sample(p::Vector{Float64}, N::Int)
    K = length(p)
    S = zeros(K, N)
    S_tmp = rand(Categorical(p), N)
    for k in 1 : K
        S[k,find(S_tmp.==k)] = 1
    end
    return S
end

function init_S(X::Matrix{Float64}, bpmm::BPMM)
    N = size(X, 2)
    K = bpmm.K
    S = categorical_sample(ones(K)/K, N)
    return S
end

function calc_ELBO(X::Matrix{Float64}, pri::BPMM, pos::BPMM)
    ln_expt_S = update_S(pos, X)
    expt_S = exp.(ln_expt_S)
    K, N = size(expt_S)
    D = size(X, 1)

    expt_ln_lambda = zeros(D, K)
    expt_lambda = zeros(D, K)
    expt_ln_lkh = 0
    for k in 1 : K
        expt_ln_lambda[:,k] = digamma.(pos.cmp[k].a) - log.(pos.cmp[k].b)
        expt_lambda[:,k] = pos.cmp[k].a / pos.cmp[k].b
        for n in 1 : N
            expt_ln_lkh += expt_S[k,n] * (X[:, n]' * expt_ln_lambda[:,k]
                                       - sum(expt_lambda[:,k]) - sum(lgamma.(X[:,n]+1)))[1]
        end
    end

    expt_ln_pS = sum(expt_S' * (digamma.(pos.alpha) - digamma.(sum(pos.alpha))))
    expt_ln_qS = sum(expt_S .* ln_expt_S)

    KL_lambda = 0
    for k in 1 : K
        KL_lambda += (sum(pos.cmp[k].a)*log.(pos.cmp[k].b) - sum(pri.cmp[k].a)*log.(pri.cmp[k].b)
                      - sum(lgamma.(pos.cmp[k].a)) + sum(lgamma.(pri.cmp[k].a))
                      + (pos.cmp[k].a - pri.cmp[k].a)' * expt_ln_lambda[:,k]
                      + (pri.cmp[k].b - pos.cmp[k].b) * sum(expt_lambda[:,k])
                      )[1]
    end
    KL_pi = (lgamma.(sum(pos.alpha)) - lgamma.(sum(pri.alpha))
             - sum(lgamma.(pos.alpha)) + sum(lgamma.(pri.alpha))
             + (pos.alpha - pri.alpha)' * (digamma.(pos.alpha) - digamma.(sum(pos.alpha)))
             )[1]

    VB = expt_ln_lkh + expt_ln_pS - expt_ln_qS - (KL_lambda + KL_pi)
    return VB
end

function add_stats(bpmm::BPMM, X::Matrix{Float64}, S::Matrix{Float64})
    D = bpmm.D
    K = bpmm.K
    sum_S = sum(S, 2)
    alpha = [bpmm.alpha[k] + sum_S[k] for k in 1 : K]
    cmp = Vector{Gam}()

    XS = X*S';
    for k in 1 : K
        a = [(bpmm.cmp[k].a[d] + XS[d,k])::Float64 for d in 1 : D]
        b = bpmm.cmp[k].b + sum_S[k]
        push!(cmp, Gam(a, b))
    end
    return BPMM(D, K, alpha, cmp)
end

remove_stats(bpmm::BPMM, X::Matrix{Float64}, S::Matrix{Float64}) = add_stats(bpmm, X, -S)

####################
## used for Variational Inference
function update_S(bpmm::BPMM, X::Matrix{Float64})
    D, N = size(X)
    K = bpmm.K
    ln_expt_S = zeros(K, N)
    tmp = zeros(K)

    sum_digamma_tmp = digamma.(sum(bpmm.alpha))
    for k in 1 : K
        tmp[k] = - sum(bpmm.cmp[k].a) / bpmm.cmp[k].b
        tmp[k] += digamma.(bpmm.alpha[k]) - sum_digamma_tmp
    end
    ln_lambda_X = [X'*(digamma.(bpmm.cmp[k].a) - log.(bpmm.cmp[k].b)) for k in 1 : K]
    for n in 1 : N
        tmp_ln_pi =  [tmp[k] + ln_lambda_X[k][n] for k in 1 : K]
        ln_expt_S[:,n] = tmp_ln_pi - logsumexp(tmp_ln_pi)
    end
    return ln_expt_S
end

"""
Pick single states having a max probability.
"""
function winner_takes_all(S::Matrix{Float64})
    S_ret = zeros(size(S))
    for n in 1 : size(S_ret, 2)
        idx = indmax(S[:,n])
        S_ret[idx,n] = 1
    end
    return S_ret
end

####################
## used for Gibbs Sampling
function sample_S_GS(pmm::PMM, X::Matrix{Float64})
    D, N = size(X)
    K = pmm.K
    S = zeros(K, N)

    tmp = [-sum(pmm.cmp[k].lambda) + log.(pmm.phi[k]) for k in 1 : K]
    ln_lambda_X = [X'*log.(pmm.cmp[k].lambda) for k in 1 : K]
    for n in 1 : N
        tmp_ln_phi = [(tmp[k] + ln_lambda_X[k][n])::Float64 for k in 1 : K]
        tmp_ln_phi = tmp_ln_phi - logsumexp(tmp_ln_phi)
        S[:,n] = categorical_sample(exp.(tmp_ln_phi))
    end
    return S
end

####################
## used for Collapsed Gibbs Sampling
function calc_ln_NB(Xn::Vector{Float64}, gam::Gam)
    ln_lkh = [(gam.a[d]*log.(gam.b)
               - lgamma.(gam.a[d])
               + lgamma.(Xn[d] + gam.a[d])
               - (Xn[d] + gam.a[d])*log.(gam.b + 1)
               )::Float64 for d in 1 : size(Xn, 1)]
    return sum(ln_lkh)
end

function sample_Sn(Xn::Vector{Float64}, bpmm::BPMM)
    ln_tmp = [(calc_ln_NB(Xn, bpmm.cmp[k]) + log.(bpmm.alpha[k])) for k in 1 : bpmm.K]
    ln_tmp = ln_tmp -  logsumexp(ln_tmp)
    Sn = categorical_sample(exp.(ln_tmp))
    return Sn
end

function sample_S_CGS(S::Matrix{Float64}, X::Matrix{Float64}, bpmm::BPMM)
    D, N = size(X)
    K = size(S, 1)
    for n in randperm(N)
        # remove
        bpmm = remove_stats(bpmm, X[:,[n]], S[:,[n]])
        # sample
        S[:,n] = sample_Sn(X[:,n], bpmm)
        # insert
        bpmm = add_stats(bpmm, X[:,[n]], S[:,[n]])
    end
    return S, bpmm
end

####################
## Algorithm main
"""
Compute posterior distribution via variational inference.
"""
function learn_VI(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    # initialisation
    expt_S = init_S(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, expt_S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # E-step
        expt_S = exp.(update_S(bpmm, X))
        # M-step
        bpmm = add_stats(prior_bpmm, X, expt_S)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bpmm, bpmm)
    end

    return expt_S, bpmm, VB
end

"""
Compute posterior distribution via Gibbs sampling.
"""
function learn_GS(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    # initialisation
    S = init_S(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # sample parameters
        pmm = sample_PMM(bpmm)
        # sample latent variables
        S = sample_S_GS(pmm, X)
        # update current model
        bpmm = add_stats(prior_bpmm, X, S)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bpmm, bpmm)
    end

    return S, bpmm, VB
end

"""
Compute posterior distribution via collapsed Gibbs sampling.
"""
function learn_CGS(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    # initialisation
    S = init_S(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, S)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # directly sample S
        S, bpmm = sample_S_CGS(S, X, bpmm)
        # calc VB
        VB[i] = calc_ELBO(X, prior_bpmm, bpmm)
    end

    return S, bpmm, VB
end

end


text:
"""
    Chebyshev( xmin, xmax, nx )

Chebyshev nodes

"""
struct Chebyshev <: NodesType

    nx    :: Int64
    xmin  :: Float64
    xmax  :: Float64
    xk    :: Vector{Float64}

    function Chebyshev( xmin, xmax, nx )

        xk  = zeros(Float64, nx)
        θ = (π:-π/(nx-1):-1e-14)
        xk  .= ((cos.(θ).+1)*(xmax-xmin))/2 .+ xmin
        new( nx, xmin, xmax, xk )

    end

end


text:
using Test
using MDCDL
using FFTW
using Random

@testset "convolutionalOperator" begin
    include("testsetGenerator.jl")

    Random.seed!(3923528829)

    @testset "Analysis" begin
        maxDims = 2
        for d in 1:maxDims, dt in [ Float64, Complex{Float64} ]
            cfgs = vec([ (crdf.I, crord.I .- 1, nch) for crdf in CartesianIndices(tuple(fill(4,d)...)), crord in CartesianIndices(tuple(fill(2+1,d)...)), nch in 2:10 ])

            subcfgs = randsubseq(cfgs, 30 / length(cfgs))

            for (df, ord, nch) in subcfgs
                szFilter = df .* (ord .+ 1)
                afs = [ rand(dt, szFilter) for p in 1:sum(nch) ]
                szx = rand(1:4) .* df .* (ord .+ 1)
                ca = ConvolutionalOperator(ParallelFilters(afs, df, ord, nch), szx)

                x = rand(dt, szx)
                y = analyze(ca, x)

                myfilter(A, h) = begin
                    ha = zero(A)
                    ha[UnitRange.(1, size(h))...] = h
                    if dt <: Real
                        real(ifft(fft(A).*fft(ha)))
                    else
                        ifft(fft(A).*fft(ha))
                    end
                end
                offset = df .- 1

                sys = map(afs) do af
                    fx = myfilter(x, af)
                    dwfx = downsample(fx, df, offset)
                    circshift(dwfx, (-1 .* fld.(ord,2)))
                end
                @test all(y .≈ sys)
            end
        end
    end

    @testset "Synthesis" begin
        maxDims = 2
        for d in 1:maxDims, dt in [ Float64, Complex{Float64} ]
            cfgs = vec([ (crdf.I, crord.I .- 1, nch) for crdf in CartesianIndices(tuple(fill(4,d)...)), crord in CartesianIndices(tuple(fill(2+1,d)...)), nch in 2:10 ])

            subcfgs = randsubseq(cfgs, 30 / length(cfgs))

            for (df, ord, nch) in subcfgs
                szFilter = df .* (ord .+ 1)
                sfs = [ rand(dt, szFilter) for p in 1:sum(nch) ]
                szy = rand(1:4) .* (ord .+ 1)
                # cs = ConvolutionalOperator(sfs, szy, df, ord, nch)
                cs = ConvolutionalOperator(ParallelFilters(sfs, df, ord, nch), szy)

                y = [ rand(dt, szy) for p in 1:sum(nch) ]
                rx = synthesize(cs, y)

                myfilter(A, h) = begin
                    ha = zero(A)
                    ha[UnitRange.(1, size(h))...] = h
                    if dt <: Real
                        real(ifft(fft(A).*fft(ha)))
                    else
                        ifft(fft(A).*fft(ha))
                    end
                end
                # check synthesizer
                fy = map(y, sfs) do  yp, fp
                    upyp = MDCDL.upsample(yp, df)
                    myfilter( upyp, fp )
                end
                subrxs = sum(fy)
                subrx = circshift(subrxs, -1 .* df .* cld.(ord,2))

                @test subrx ≈ rx
            end
        end
    end
end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
import Random
import RestrictedBoltzmannMachines as RBMs

using Test: @test, @testset, @inferred
using Statistics: mean, std
using LinearAlgebra, Statistics
using LogExpFunctions: softmax

@testset "onehot" begin
    q = 10
    A = rand(1:q, 10, 7, 3)
    X = @inferred RBMs.onehot_encode(A, 1:q)
    @test size(X) == (q, size(A)...)
    for i in CartesianIndices(A), a in 1:q
        @test X[a,i] == (A[i] == a)
    end
    @test A == @inferred RBMs.onehot_decode(X)
end

@testset "categorical_rand" begin
    Random.seed!(1)
    ps = [0.2, 0.5, 0.3]
    counts = Dict{Int,Int}()
    for _ = 1:1000000
        s = RBMs.categorical_rand(ps)
        counts[s] = get(counts, s, 0) + 1
    end
    for (s, c) in counts
        @test c ./ sum(values(counts)) ≈ ps[s] atol=1e-2
    end
end

@testset "categorical_sample" begin
    Random.seed!(2)
    q = 10
    P = rand(q, 4, 3, 7)
    P ./= sum(P; dims=1)
    X = mean(RBMs.onehot_encode(RBMs.categorical_sample(P), 1:q) for _ in 1:10^6)
    @test X ≈ P rtol=0.1
end

@testset "gumbel" begin
    Random.seed!(3)
    @test mean(RBMs.randgumbel() for _ = 1:10^6) ≈ MathConstants.γ rtol=0.01
    @test std(RBMs.randgumbel() for _ = 1:10^6) ≈ π / √6 rtol=0.01
end

@testset "categorical_sample_from_logits" begin
    Random.seed!(52)
    q = 3
    logits = randn(q,2)
    p = mean(RBMs.onehot_encode(RBMs.categorical_sample_from_logits(logits), 1:q) for _ in 1:10^6)
    @test p ≈ softmax(logits; dims=1) rtol=0.01
end

@testset "categorical_sample_from_logits_gumbel" begin
    Random.seed!(52)
    q = 3
    logits = randn(q,2)
    p = mean(RBMs.onehot_encode(RBMs.categorical_sample_from_logits_gumbel(logits), 1:q) for _ in 1:10^6)
    @test p ≈ softmax(logits; dims=1) rtol=0.01
end


text:
two = BigInt(2)
ten = BigInt(10)

for i in 1:9
    a = BigInt(i)

    for j in 2:50
        a += i*two^(j-1)*ten^(j-1)

        at = a % ten^j
        s = string(at)
        si = s[end]*s[1:end-1]
        s2 = string(2at)

        if si == s2
            println(s)
        end
    end
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc