text:
struct Firm{T <: NTuple}
  ζ::NTuple #inverse demand function
  τ::T #tax
  α::Float64
  γ::Float64 #production share
  δ::Float64 #depreciation
  Θ::T #technology per state

  ψ::Vector{NTuple} #kroniker product of all possible state values
  π::Matrix{CartesianIndex} #policy states
  v::Matrix{Float64} #value states

  R::Array{Float64,3}
  M::Array{T,3}
end

function getyl(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  kern::Float64 = (k^(F.γ*F.α)/pl) - (ph/pl)* (yh/F.Θ[1])^F.α
  yl::Float64 = (kern ≥ 0.0) ? (kern)^(1/F.α)*F.Θ[2] : NaN

  return yl
end

function getylSIMD(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  kern::Float64 = (k^(F.γ*F.α)/pl) - (ph/pl)* (yh/F.Θ[1])^F.α
  yl::Float64 = (kern ≥ 0.0) ? (kern)^(1/F.α)*F.Θ[2] : NaN

  return yl
end

#pass the following parameters to get the pricing kernel
function getkernel(s::Int, yₜ::Float64, yₜ₁h::Float64, yₜ₁l::Float64, p::Matrix{Float64},
  F::T where T<:Firm)::NTuple{3,Float64}

  #derive the remaining parameters
  #hopefully this will aid in debugging
  ζ₀::Float64 = F.ζ[1]
  ζ₁::Float64 = F.ζ[2]
  Pₜ::Float64 = ζ₀ - ζ₁*yₜ

  pₜ₁h::Float64 = p[1,s]
  pₜ₁l::Float64 = 1 - pₜ₁h

  Θₜ::Float64 = F.Θ[s]
  Θh::Float64 = F.Θ[1]
  Θl::Float64 = F.Θ[2]

  τh::Float64 = F.τ[1]
  τl::Float64 = F.τ[2]

  #1/(F.ζ[1]-F.ζ[2]*yₜ₁h) * (yₜ₁h/yₜ)^(F.α-1) * (F.Θ[1]/F.Θ[s])^(-F.α)
  #now compute the intermediate quantities
  ξʰ::Float64 = (1/(ζ₀ - ζ₁*yₜ₁h)) * (yₜ₁h/yₜ)^(F.α-1) * (Θh/Θₜ)^(-F.α)
  ξˡ::Float64 = (1/(ζ₀ - ζ₁*yₜ₁l)) * (yₜ₁l/yₜ)^(F.α-1) * (Θl/Θₜ)^(-F.α)
  ξ::Float64 = ξʰ*pₜ₁h + ξˡ*pₜ₁l

  κʰ::Float64 = (ζ₀ - 2*ζ₁*yₜ₁h)*τh
  κˡ::Float64 = (ζ₀ - 2*ζ₁*yₜ₁l)*τl

  Mdenom::Float64 = ξ*(1-τh)*(1-τl) + pₜ₁l*κˡ*ξˡ*(1-τh) + pₜ₁h*κʰ*ξʰ*(1-τl)

  Mₜ₁h::Float64 = Pₜ*ξʰ*(1-τl)/(Mdenom)
  Mₜ₁l::Float64 = Pₜ*ξˡ*(1-τh)/(Mdenom)

  return (Mₜ₁h, Mₜ₁l, Mₜ₁l)

end

function fillRCol!(F::T, p::Matrix{Float64}, s::Int, c::Int)::T where T<:Firm
  N::Int = length(F.ψ)
  @assert s ≤ size(p,1)
  @assert c ≤ N

  @inbounds @simd for r ∈ 1:N #2x performance gain or more
  #make some type assertions to help with debugging and optimization
  local τₜ::Float64
  local yₜ::Float64 #will hold the true produciton, not just the high state production
  local yₜ₁l::Float64 #future low stae productivity (fully specified by high state prod and state)
  local yₜl::Float64 #current low state productivitiy

  local pₜ₁h::Float64 #probability of transitioning to a high state
  local pₜ₁l::Float64 #probability of transitioning to a low state
  local pₜh::Float64 #probability of having transitioned from a high state
  local pₜl::Float64 #probability of having transitioned from a low state

  local dₜ::Float64 #the dividend payment
  local Mₜ₁::NTuple{3, Float64} #pricing kernel given the current state

  local Pₜ::Float64 #current price
  local Pₜ₁h::Float64 #future price in high state
  local Pₜ₁l::Float64 #future price in low state

  local allocationₜ::Float64 #used to check the feasability of the current allocaiton
  local allocationₜ₁::Float64 #used to check the feasability of the future allocaiton

  (kₜ₁::Float64, yₜ₁h::Float64) = F. ψ[r]
  (kₜ::Float64, yₜh::Float64) =  F.ψ[c]
  τₜ = F.τ[s]

  # compute the future low state production from the high state production
  pₜ₁h = p[1,s]
  pₜ₁l = 1.0 - p[1,s]
  yₜ₁l = getyl(yₜ₁h, kₜ₁, pₜ₁h, pₜ₁l, F::Firm)

  #compute the current low income state
  pₜh = (s == 2) ? p[1,1] : p[1,2]
  pₜl = 1.0 - pₜh
  yₜl = getyl(yₜh, kₜ, pₜh, 1 - pₜh, F::Firm)

  if s == 1 #need to recover the previously selected low productivity
    yₜ = yₜh
  else
    yₜ = yₜl
  end

  yₜ = max(yₜ, N_TOL)
  yₜ₁l = max(yₜ₁l, N_TOL)
  yₜ₁h = max(yₜ₁h, N_TOL)

  #prices
  Pₜ = F.ζ[1] - F.ζ[2] * yₜ
  Pₜ₁h = F.ζ[1] - F.ζ[2] * yₜ₁h
  Pₜ₁l = F.ζ[1] - F.ζ[2] * yₜ₁l

  if isnan(yₜl) || isnan(yₜ₁l) || (Pₜ<0) || (Pₜ₁h<0) || (Pₜ₁l<0)
    Mₜ₁ = (0.0,0.0,0.0)
    dₜ = 0.0
  else
    dₜ = Pₜ*yₜ*(1-τₜ)-(kₜ₁-(1-F.δ)*kₜ)+τₜ*F.δ*kₜ
    Mₜ₁ = getkernel(s, yₜ, yₜ₁h, yₜ₁l, p, F)
  end

  #if sum((isnan).([Mₜ₁...;dₜ])) ≠ 0
    #println("ERROR: Mₜ₁: $Mₜ₁, dₜ: $dₜ")
  #  println("(s, yₜ, yₜ₁h, yₜ₁l, p):", "$((s, yₜ, yₜ₁h, yₜ₁l, p))")
  #end
  #@assert sum((isnan).([Mₜ₁...;dₜ])) == 0 #no NaNs allowed
  #RFR::Float64 = (sum(Mₜ₁) > 0.01) ?  1/(Mₜ₁[1]*pₜ₁h + Mₜ₁[2]*pₜ₁l) : NaN
  #println("RFR: $RFR Mₜ₁: $Mₜ₁, dₜ: $dₜ")

  F.R[r,c,s] = dₜ
  F.M[r,c,s] = Mₜ₁

  end

  return F
end

function fillR!(F::T, p::Matrix{Float64})::T where T<:Firm
  N::Int = length(F.ψ)
  S::Int = length(F.Θ)

  @inbounds for s ∈ 1:S, c ∈ 1:N
    fillRCol!(F, p, s, c)
  end

  return F
end

function Firm(;  ζ::NTuple= DEF_ζ,
  τ::T = DEF_τ,
  δ::Float64 = DEF_δ,
  Θ::T = DEF_Θ,
  α::Float64=DEF_α,
  γ::Float64 = DEF_γ,
  kstates::Vector{Float64} = DEF_KSTATES,
  yhstates::Vector{Float64} = DEF_YSTATES,
  p::Matrix{Float64} = DEF_p)::Firm{T} where T<:NTuple #for performance reasons, don't use vectors

  #rules that need to be followed for this approach to work
  @assert Θ[2] == Θ[3] #tech level has to be the same between states
  @assert τ[2] == τ[3] #tax is identical between the two states
  @assert sum(p[:,2] .≠ p[:,3]) == 0 #transition probabilities are the same if originating from low
  @assert sum(p) == 3.0 #probabilities add to 1
  @assert ζ[2] ≥ 0.0 #demand curve is weakly downward sloping

  ψ::Vector{NTuple} = combine2tuples(kstates, yhstates)
  N::Int = length(ψ)
  S::Int = length(Θ)

  π::Matrix{CartesianIndex} = Matrix{CartesianIndex}(undef, N, S)
  v::Matrix{Float64} = zeros(N,S)
  R::Array{Float64,3} = Array{Float64,3}(undef, N, N, S)
  M::Array{T,3} = fill(T(zeros(S)), N, N, S) #preallocate

  F::Firm{T} = Firm{T}(ζ, τ, α, γ, δ, Θ, ψ, π, v, R, M)

  #build up R
  fillR!(F, p)

  #println("getYl: ",getyl(1.4, 2.3, .3, .9, F)) #NOTE: Checks this formula
  #println("getYlcheck: ",getylcheck(1.4, 2.3, .3, .9, F))

  #println("getkernel: ",getkernel(1, 2.5, 2.7, 3.4, p, F)) #NOTE: Checks this formula
  #println("getkernel: ",getkernelcheck(1, 2.5, 2.7, 3.4, p, F), "\n")
  #Mdisp::Matrix{Float64} = Matrix{Float64}(undef,N,N)
  #for r ∈ 1:N, c∈ 1:N
  #  Mdisp[r,c] = F.M[r,c,2][2]
  #end

  #display(Mdisp)
  return F
end

#=
function getkernelcheck(s::Int, yₜ::Float64, yₜ₁h::Float64, yₜ₁l::Float64, p::Matrix{Float64},
  F::T where T<:Firm)::NTuple{3,Float64}

  #define derived parameters
  pₜ₁h = p[1,s]
  pₜ₁l = 1 - pₜ₁h
  Pₜ::Float64 = F.ζ[1] - F.ζ[2]*yₜ

  ξʰ::Float64 = 1/(F.ζ[1]-F.ζ[2]*yₜ₁h) * (yₜ₁h/yₜ)^(F.α-1) * (F.Θ[1]/F.Θ[s])^(-F.α)
  ξˡ::Float64 = 1/(F.ζ[1]-F.ζ[2]*yₜ₁l) * (yₜ₁l/yₜ)^(F.α-1) * (F.Θ[2]/F.Θ[s])^(-F.α)
  ξ::Float64 = pₜ₁h*ξʰ + pₜ₁l*ξˡ
  κʰ::Float64 = (F.ζ[1]-2*F.ζ[2]*yₜ₁h)*F.τ[1]
  κˡ::Float64 = (F.ζ[1]-2*F.ζ[2]*yₜ₁l)*F.τ[2]

  Mₜ₁l::Float64 = (Pₜ*ξˡ*(1-F.τ[1])/
    (ξ*(1-F.τ[1])*(1-F.τ[2]) + pₜ₁l*ξˡ*κˡ*(1-F.τ[1]) + pₜ₁h*ξʰ*κʰ*(1-F.τ[2])))
  Mₜ₁h::Float64 = ((Pₜ-pₜ₁l*κˡ*Mₜ₁l)*ξʰ /
    (ξ*(1-F.τ[1]) + pₜ₁h*ξʰ*κʰ))

  #println("Mdenom: $((ξ*(1-F.τ[1])*(1-F.τ[2]) + pₜ₁l*ξˡ*κˡ*(1-F.τ[1]) + pₜ₁h*ξʰ*κʰ*(1-F.τ[2])))")
  return (Mₜ₁h, Mₜ₁l, Mₜ₁l)
end=#

#=
function getylcheck(yh::Float64, k::Float64, ph::Float64, pl::Float64, F::Firm)
  Yl::Float64 = ((1/pl)*(k^(F.γ*F.α)-ph*(yh/F.Θ[1])^F.α))^(1/F.α)*F.Θ[2]
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using StaticArrays
using LinearAlgebra
using Printf
using Images
using Manifolds

using SLAM

include("kitty.jl")

function kitty_slam(kitty_dir, save_dir, sequence, n_frames, stereo = true)
    dataset = KittyDataset(kitty_dir, sequence; stereo)
    println(dataset)

    save_dir = joinpath(save_dir, "kitty-$sequence")
    isdir(save_dir) || mkdir(save_dir)

    fx, fy = dataset.K[1, 1], dataset.K[2, 2]
    cx, cy = dataset.K[1:2, 3]
    # TODO KITTY dataset should do that
    height, width = 376, 1241
    # height, width = 370, 1226

    camera = SLAM.Camera(;fx, fy, cx, cy, height, width)
    right_camera = SLAM.Camera(;fx, fy, cx, cy, height, width, Ti0=dataset.Ti0)
    params = Params(;
        stereo, do_local_bundle_adjustment=true, map_filtering=true)

    saver = ReplaySaver()
    slam_manager = SlamManager(params, camera; right_camera, slam_io=saver)
    slam_manager_thread = Threads.@spawn run!(slam_manager)

    t1 = time()
    for i in 1:n_frames
        timestamp = dataset.timestamps[i]
        left_frame, right_frame = dataset[i]
        left_frame = Gray{Float64}.(left_frame)

        if params.stereo
            right_frame = Gray{Float64}.(right_frame)
            add_stereo_image!(slam_manager, left_frame, right_frame, timestamp)
        else
            add_image!(slam_manager, left_frame, timestamp)
        end

        q_size = get_queue_size(slam_manager)
        f_size = length(slam_manager.mapper.estimator.frame_queue)
        m_size = length(slam_manager.mapper.keyframe_queue)
        while q_size > 0 || f_size > 0 || m_size > 0
            sleep(1e-2)
            q_size = get_queue_size(slam_manager)
            f_size = length(slam_manager.mapper.estimator.frame_queue)
            m_size = length(slam_manager.mapper.keyframe_queue)
        end

        sleep(1e-2)
    end

    slam_manager.exit_required = true
    wait(slam_manager_thread)

    t2 = time()
    @info "SLAM took: $(t2 - t1) seconds."

    SLAM.save(saver, save_dir)
    @info "Saved ReplaySaver."
    slam_manager
end

ext: jl


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    # 2x + y
    @test my_f(2,1) == 5
    @test my_f(2,2) == 6
    @test my_f(2,3) == 7
end

@testset "Derivative test" begin
    # 2x + y
    @test derivative_of_my_f(2,1) == 2
end

max_stars_repo_name: filippigianluca/MyExample

max_issues_repo_name: filippigianluca/MyExample


text:
# This file is a part of ChangesOfVariables.jl, licensed under the MIT License (MIT).

using ChangesOfVariables
using Test

using LinearAlgebra

using ChangesOfVariables
using ChangesOfVariables: test_with_logabsdet_jacobian
using ChainRulesTestUtils

include("getjacobian.jl")

@testset "with_logabsdet_jacobian" begin
    @test with_logabsdet_jacobian(sum, rand(5)) == NoLogAbsDetJacobian{typeof(sum),Vector{Float64}}()
    @test_throws MethodError _, _ = with_logabsdet_jacobian(sum, rand(5))

    function ChangesOfVariables.with_logabsdet_jacobian(::typeof(foo), x)
        y = foo(x)
        ladj = -x + 2 * log(y)
        (y, ladj)
    end

    x = 4.2
    X = rand(10)
    A = rand(5, 5)
    CA = Complex.(rand(5, 5), rand(5, 5))

    isaprx(a, b; kwargs...) = isapprox(a,b; kwargs...)
    isaprx(a::NTuple{N,Any}, b::NTuple{N,Any}; kwargs...) where N = all(map((a,b) -> isaprx(a, b; kwargs...), a, b))

    test_with_logabsdet_jacobian(foo, x, getjacobian)

    @static if VERSION >= v"1.6"
        test_with_logabsdet_jacobian(log ∘ foo, x, getjacobian)
    end

    @testset "with_logabsdet_jacobian on mapped and broadcasted" begin
        for f in (Base.Fix1(map, foo), Base.Fix1(broadcast, foo))
            for arg in (x, fill(x,), Ref(x), (x,), X)
                test_with_logabsdet_jacobian(f, arg, getjacobian, compare = isaprx)
            end
        end
    end

    @testset "with_logabsdet_jacobian on identity, adjoint and transpose" begin
        for f in (identity, adjoint, transpose)
            for arg in (x, A)
                test_with_logabsdet_jacobian(f, arg, getjacobian)
            end
        end
    end

    @testset "with_logabsdet_jacobian on inv" begin
        for arg in (x, A, CA)
            test_with_logabsdet_jacobian(inv, arg, getjacobian)
        end
    end

    @testset "with_logabsdet_jacobian on log and exp functions" begin
        for f in (+, -, exp, log, exp2, log2, exp10, log10, expm1, log1p)
            test_with_logabsdet_jacobian(f, x, getjacobian)
        end
    end

    @testset "rrules" begin
        for map_or_bc in (map, broadcast)
            test_rrule(ChangesOfVariables._with_ladj_on_mapped, map_or_bc, [(1.0, 2.0), (3.0, 4.0), (5.0, 6.0)])
        end
    end
end


text:
export devicelist

Base.@kwdef mutable struct DWFDevice
    enumid::Cint
    name::String
    serial::String
    hdwf = Ref(HDWF(0))
end
display(dev::DWFDevice) = println("ID: $(dev.enumid)\nName: $(dev.name)\nSerial Number: $(dev.serial)")

"""
    enumdev(i<:Integer)

Return device by enumeration id.
"""
function enumdev(i::T) where {T<:Integer}
    # Digilent WaveForms SDK uses 0 based indexing
    szDeviceName = zeros(Cchar, 32)
    szSN = zeros(Cchar, 32)
    FDwfEnumDeviceName(i, szDeviceName)
    FDwfEnumSN(i, szSN)
    return DWFDevice(enumid=i, name=str(szDeviceName), serial=name=str(szSN))
end

"""
    devicelist(;info::Bool=false, silent::Bool=false)

Show device list and return vector of device objects.
"""
function devicelist(;info::Bool=false, silent::Bool=false)
    cDevice = Ref(Cint(0))
    szError = zeros(Cchar, 512)

    # detect connected all supported devices
    if FDwfEnum(enumfilterAll, cDevice) == 0
        FDwfGetLastErrorMsg(szError)
        error("FDwfEnum: $szError")
    end
    # list information about each device
    !silent && @info "Found $(cDevice[]) devices:"
    devices = DWFDevice[]
    for i in 0:cDevice[]-1

        dev = enumdev(i)
        push!(devices, dev)

        !silent && println("Device: $(i+1): $(dev)")
        # before opening, check if the device isn�t already opened by other application, like: WaveForms

        if info
            if !open(dev)
                open!(dev)
                infostring = devinfo(dev)
                println(infostring)
                close!(dev)
            else
                @warn "Device already open. Couldn't collect info"
            end
        end

    end
    # before application exit make sure to close all opened devices by this process
    closeall!()
    return devices
end

"""
    open(dev::DWFDevice)

Return whether device is open
"""
function open(dev::DWFDevice)
    fIsInUse = Ref(Cint(0))
    FDwfEnumDeviceIsOpened(dev.enumid, fIsInUse)
    return fIsInUse[] == 1
end

"""
    checkopen(dev::DWFDevice)

Check & error whether device is open
"""
function checkopen(dev::DWFDevice)
    #TODO: Check not working
    #!open(dev) && error("Device isn't open")
    return nothing
end

"""
    open!(dev::DWFDevice)

Open device.
"""
function open!(dev::DWFDevice)
    if FDwfDeviceOpen(dev.enumid, dev.hdwf) != 1
        szError = zeros(Cchar, 512)
        FDwfGetLastErrorMsg(szError)
        @error "Device in use. FDwfDeviceOpen error: $(str(szError))"
    end
    return nothing
end

"""
    close!(dev::DWFDevice)

Close device.
"""
function close!(dev::DWFDevice)
    FDwfDeviceClose(dev.hdwf[])
    dev.hdwf = Ref(hdwfNone)
    return nothing
end

"""
    closeall!()

Close all devices.
"""
function closeall!()
    FDwfDeviceCloseAll()
    return nothing
end

"""
    devinfo(dev::DWFDevice)

Return some device information.
"""
function devinfo(dev::DWFDevice)
    numchannels = AnalogInChannelCount(dev)
    maxfreq = AnalogInFrequencyInfo(dev)
    return "Number of analog input channels: $(numchannels). Maximum freq.: $(maxfreq) Hz"
end


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

# This file tests the new Julia 1.6+ extension interface of TerminalMenus
# To trigger the new interface, at least one configuration keyword argument must be supplied.

# Check to make sure types are imported properly
@test RadioMenu{TerminalMenus.Config} <: TerminalMenus.ConfiguredMenu  # TODO Julia 2.0: delete parameter

# Constructor
@test RadioMenu(["one", "two", "three"]; charset=:ascii).pagesize == 3
@test RadioMenu(string.(1:30), pagesize=-1, charset=:ascii).pagesize == 30
@test RadioMenu(string.(1:4), pagesize=10, charset=:ascii).pagesize == 4
@test RadioMenu(string.(1:100); charset=:ascii).pagesize == 10

radio_menu = RadioMenu(string.(1:20); charset=:ascii)
@test TerminalMenus.options(radio_menu) == string.(1:20)
radio_menu.selected = 2
TerminalMenus.cancel(radio_menu)
@test radio_menu.selected == -1
@test TerminalMenus.header(radio_menu) == ""

# Output
for kws in ((charset=:ascii,),
            (charset=:unicode,),
            (cursor='@',))
    local radio_menu
    radio_menu = RadioMenu(string.(1:10); kws...)
    c = isdefined(kws, :cursor) ? kws.cursor : (kws.charset === :ascii ? '>' : '→')
    buf = IOBuffer()
    TerminalMenus.writeline(buf, radio_menu, 1, true)
    @test String(take!(buf)) == "1"
    TerminalMenus.printmenu(buf, radio_menu, 1; init=true)
    @test startswith(String(take!(buf)), "\e[2K $c 1")
    TerminalMenus.printmenu(buf, radio_menu, 2; init=true)
    @test startswith(String(take!(buf)), string("\e[2K   1\r\n\e[2K $c 2"))
end
@testset begin "cursor page"
    radio_menu = RadioMenu(string.(1:20); charset=:ascii)
    buf = IOBuffer()
    TerminalMenus.printmenu(buf, radio_menu, 19; init=true)
    @test String(take!(buf)) == "\e[2K^  11\r\n\e[2K   12\r\n\e[2K   13\r\n\e[2K   14\r\n\e[2K   15\r\n\e[2K   16\r\n\e[2K   17\r\n\e[2K   18\r\n\e[2K > 19\r\n\e[2K   20"
    TerminalMenus.printmenu(buf, radio_menu, 8; init=true)
    @test String(take!(buf)) == "\e[2K^  4\r\n\e[2K   5\r\n\e[2K   6\r\n\e[2K   7\r\n\e[2K > 8\r\n\e[2K   9\r\n\e[2K   10\r\n\e[2K   11\r\n\e[2K   12\r\n\e[2Kv  13"
end

# Test using stdin
radio_menu = RadioMenu(string.(1:10); charset=:ascii)
@test simulate_input(3, radio_menu, :down, :down, :enter)
radio_menu = RadioMenu(["single option"], charset=:ascii)
@test simulate_input(1, radio_menu, :up, :up, :down, :up, :enter)
radio_menu = RadioMenu(string.(1:3), pagesize=1, charset=:ascii)
@test simulate_input(3, radio_menu, :down, :down, :down, :down, :enter)
radio_menu = RadioMenu(["apple", "banana", "cherry"]; keybindings=collect('a':'c'), charset=:ascii)
@test simulate_input(2, radio_menu, 'b')

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


text:
import GivensBoson.find_offdiagonal_maximun
@testset "Functional Test" begin
    value,index,flag = find_offdiagonal_maximun([1 0.9;0.9 1])
    @test (value,index,flag)==(0.9,[1,2],1)

    value,index,flag = find_offdiagonal_maximun([1 1.0;1.0 1])
    @test (value,index,flag)==(1.0,[1,2],1)

    value,index,flag = find_offdiagonal_maximun([1 1.0;1.0 1],zeromode=true)
    @test (value,index,flag)==(1.0,[1,2],2)

    value,index,flag = find_offdiagonal_maximun([1 .0;0 1])
    @test (value,index,flag)==(0.0,[1,1],3)
end


text:
#workspace()

using IntervalArithmetic
using Clp
using EAGODomainReduction

type t_solver
    LP_solver
end

type test_opt
    numVar
    numConstr
    f
    g
    solver::t_solver
    gL
    gU
    gL_loc
    gU_loc
end

f(x) = x[1]
g(x) = [x[1]+x[2];
        x[1]-x[2]]

in_opt = t_solver(ClpSolver())
opt = test_opt(2,2,f,g,in_opt,[2 -1],[3 0],[1,2],[1,2])
X = [Interval(0.0,2.0) for i=1:2]
STD_LP_Probe!(X,[opt],1.0)
println("Standard Probing Reduced X: ", X)


text:
using Pardiso: Pardiso, set_solver!, pardisoinit, set_phase!, set_iparm!, set_msglvl!, get_iparm, pardiso, PardisoSolver, MKLPardisoSolver, get_nprocs, set_nprocs!, set_matrixtype!
export PardisoDirectKKTSolver, PardisoIndirectKKTSolver, MKLPardisoSolver

abstract type AbstractPardisoKKTSolver <: AbstractKKTSolver end

#---------------------------
#Pardiso 5.0 Direct Solver Configuration
#---------------------------

struct PardisoDirectKKTSolver{Tv, Ti} <: AbstractPardisoKKTSolver

    ps::PardisoSolver
    K::SparseMatrixCSC{Tv, Ti}
    m::Integer
    n::Integer
    work::Vector{Tv}  #working memory for Pardiso

    function PardisoDirectKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma, rho;
        iparm::Dict{Int64, Int64} = Dict{Int64, Int64}(), msg_level_on::Bool = false) where {Tv, Ti}

        m, n, K, ps, work = _pardiso_common_init(P, A, sigma, rho, PardisoSolver, iparm, msg_level_on)

        set_solver!(ps,Pardiso.DIRECT_SOLVER)
        _kktutils_direct_solver_init(ps, K, work, m, n)

        return new{Tv, Ti}(ps, K, m, n, work)
    end
end

#---------------------------
#Pardiso 5.0 Indirect Solver Configuration
#---------------------------

struct PardisoIndirectKKTSolver{Tv, Ti} <: AbstractPardisoKKTSolver

    ps::PardisoSolver
    K::SparseMatrixCSC{Tv, Ti}
    m::Integer
    n::Integer
    work::Vector{Tv}  #working memory for Pardiso

    function PardisoIndirectKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma, rho;
        iparm::Dict{Int64, Int64} = Dict{Int64, Int64}(), msg_level_on::Bool = false) where {Tv, Ti}

        m, n, K, ps, work = _pardiso_common_init(P, A, sigma, rho, PardisoSolver, iparm, msg_level_on)

        set_solver!(ps, Pardiso.ITERATIVE_SOLVER)
        pardisoinit(ps)

        ## set phase to solving for iterations
        set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)

        return new{Tv, Ti}(ps, K, m, n, work)
    end
end

#---------------------------
#MKL Pardiso Solver Configuration
#---------------------------

struct MKLPardisoKKTSolver{Tv, Ti} <: AbstractPardisoKKTSolver

    ps::MKLPardisoSolver
    K::SparseMatrixCSC{Tv, Ti}
    m::Integer
    n::Integer
    work::Vector{Tv}  #working memory for Pardiso

    function MKLPardisoKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma, rho;
        iparm::Dict{Int64, Int64} = Dict{Int64, Int64}(), msg_level_on::Bool = false, num_threads::Int64 = 1) where {Tv, Ti}

        m, n, K, ps, work = _pardiso_common_init(P, A, sigma, rho, MKLPardisoSolver, iparm, msg_level_on)

        _kktutils_direct_solver_init(ps, K, work, m, n)
        if num_threads >= 1
            set_nprocs!(ps, num_threads)
        else
            error("The number of CPU threads has to be greater than 0.")
        end

        return new{Tv, Ti}(ps ,K, m, n, work)
    end
end

#-------------------
#Pardiso common utils

function _pardiso_common_init(P, A, sigma, rho, Solver::Type, iparm::Dict{Int64, Int64}, msg_level_on::Bool)

    #allow Pardiso solvers do this first

    m,n  = _kktutils_check_dims(P, A, sigma, rho)
    K    = _kktutils_make_kkt(P, A, sigma, rho, :L)
    work = zeros(eltype(A), m + n)
    ps   = Solver()

    #set to symmetric indefinite
    set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF)

    if msg_level_on
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    else
        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_OFF)
    end

    # set to allow non-default iparams
    set_iparm!(ps, 1, 1)

    # set user specific parameters
    for (i, v) in iparm
        set_iparm!(ps, i, v)
    end

    return m, n, K, ps, work
end

function _kktutils_direct_solver_init(ps, K, work, m, n)

    # set the transpose flag (Pardiso: CSR, Julia CSC)
    set_iparm!(ps, 12, 1)

    # Analyze the matrix and compute a symbolic factorization.
    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, K, work)

    # Compute the numeric factorization.
     set_phase!(ps, Pardiso.NUM_FACT)
     pardiso(ps, K, work)
    _pardiso_check_inertia(ps, m, n)

    ## set phase to solving for iterations
    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
end

solve!(s::AbstractPardisoKKTSolver, lhs, rhs) = pardiso(s.ps, lhs, s.K, rhs)

function update_rho!(s::AbstractPardisoKKTSolver, rho::Union{Vector,AbstractFloat})
    update_kkt_matrix(s.K, s.n, s.m, rho)

    #only refactor for the direct solver
    if s.ps.solver == Pardiso.DIRECT_SOLVER
        # Compute the numeric factorization again,
        #but skipping the analysis phase
        set_phase!(s.ps, Pardiso.NUM_FACT)
        pardiso(s.ps, s.K, s.work)
        _pardiso_check_inertia(s.ps, s.m, s.n)

        ## set back to solving phase
        set_phase!(s.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    end
end

function _pardiso_check_inertia(ps, m, n)

    num_pos_eigenvalues = get_iparm(ps, 22)
    num_neg_eigenvalues = get_iparm(ps, 23)

    num_neg_eigenvalues == m || throw("P matrix appears to have negative eigenvalues")
end

function free_memory!(s::AbstractPardisoKKTSolver)
    set_phase!(s.ps, Pardiso.RELEASE_ALL)
    pardiso(s.ps)
end

#number of solver threads active
get_number_of_threads(s::AbstractPardisoKKTSolver) = get_nprocs(s.ps)

#MKL can set number of threads directly
set_number_of_threads(s::MKLPardisoKKTSolver, i) = set_nprocs!(s.ps, i)

#Non-MKL can only set number of threads via an environment variable
set_number_of_threads(s::AbstractPardisoKKTSolver, i) =
    error("Pardiso 5.0 (non-MKL) must set thread count via environment variable, e.g. ENV[\"OMP_NUM_THREADS\"] = 4")


text:
using Lux, NNlib, Random, Test, Zygote

rng = Random.default_rng()
Random.seed!(rng, 0)

include("../utils.jl")

@testset "Pooling" begin
    x = randn(rng, Float32, 10, 10, 3, 2)
    y = randn(rng, Float32, 20, 20, 3, 2)

    layer = AdaptiveMaxPool((5, 5))
    ps, st = Lux.setup(rng, layer)

    @test layer(x, ps, st)[1] == maxpool(x, PoolDims(x, 2))
    run_JET_tests(layer, x, ps, st)

    layer = AdaptiveMeanPool((5, 5))
    ps, st = Lux.setup(rng, layer)

    @test layer(x, ps, st)[1] == meanpool(x, PoolDims(x, 2))
    run_JET_tests(layer, x, ps, st)

    layer = AdaptiveMaxPool((10, 5))
    ps, st = Lux.setup(rng, layer)

    @test layer(y, ps, st)[1] == maxpool(y, PoolDims(y, (2, 4)))
    run_JET_tests(layer, x, ps, st)

    layer = AdaptiveMeanPool((10, 5))
    ps, st = Lux.setup(rng, layer)

    @test layer(y, ps, st)[1] == meanpool(y, PoolDims(y, (2, 4)))
    run_JET_tests(layer, x, ps, st)

    layer = GlobalMaxPool()
    ps, st = Lux.setup(rng, layer)

    @test size(layer(x, ps, st)[1]) == (1, 1, 3, 2)
    run_JET_tests(layer, x, ps, st)

    layer = GlobalMeanPool()
    ps, st = Lux.setup(rng, layer)

    @test size(layer(x, ps, st)[1]) == (1, 1, 3, 2)
    run_JET_tests(layer, x, ps, st)

    layer = MaxPool((2, 2))
    ps, st = Lux.setup(rng, layer)

    @test layer(x, ps, st)[1] == maxpool(x, PoolDims(x, 2))
    run_JET_tests(layer, x, ps, st)

    layer = MeanPool((2, 2))
    ps, st = Lux.setup(rng, layer)

    @test layer(x, ps, st)[1] == meanpool(x, PoolDims(x, 2))
    run_JET_tests(layer, x, ps, st)

    @testset "$ltype SamePad windowsize $k" for ltype in (MeanPool, MaxPool),
                                                k in ((1,), (2,), (3,), (4, 5), (6, 7, 8))

        x = ones(Float32, (k .+ 3)..., 1, 1)

        layer = ltype(k; pad=Lux.SamePad())
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1])[1:(end - 2)] == cld.(size(x)[1:(end - 2)], k)
        run_JET_tests(layer, x, ps, st; call_broken=length(k) == 1)
    end
end

@testset "CNN" begin
    @testset "Grouped Conv" begin
        x = rand(rng, Float32, 4, 6, 1)
        layer = Conv((3,), 6 => 2; groups=2)
        ps, st = Lux.setup(rng, layer)

        @test size(ps.weight) == (3, 3, 2)
        @test size(layer(x, ps, st)[1]) == (2, 2, 1)
        run_JET_tests(layer, x, ps, st; call_broken=true)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        x = rand(rng, Float32, 4, 4, 6, 1)
        layer = Conv((3, 3), 6 => 2; groups=2)
        ps, st = Lux.setup(rng, layer)

        @test size(ps.weight) == (3, 3, 3, 2)
        @test size(layer(x, ps, st)[1]) == (2, 2, 2, 1)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        x = rand(rng, Float32, 4, 4, 4, 6, 1)
        layer = Conv((3, 3, 3), 6 => 2; groups=2)
        ps, st = Lux.setup(rng, layer)

        @test size(ps.weight) == (3, 3, 3, 3, 2)
        @test size(layer(x, ps, st)[1]) == (2, 2, 2, 2, 1)
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        # Test that we cannot ask for non-integer multiplication factors
        layer = Conv((2, 2), 3 => 10, groups=2)
        @test_throws AssertionError Lux.setup(rng, layer)
        layer = Conv((2, 2), 2 => 9, groups=2)
        @test_throws AssertionError Lux.setup(rng, layer)
    end

    @testset "Asymmetric Padding" begin
        layer = Conv((3, 3), 1 => 1, relu; pad=(0, 1, 1, 2))
        x = ones(Float32, 28, 28, 1, 1)
        ps, st = Lux.setup(rng, layer)

        ps.weight .= 1.0
        ps.bias .= 0.0

        y_hat = layer(x, ps, st)[1][:, :, 1, 1]
        @test size(y_hat) == (27, 29)
        @test y_hat[1, 1] ≈ 6.0
        @test y_hat[2, 2] ≈ 9.0
        @test y_hat[end, 1] ≈ 4.0
        @test y_hat[1, end] ≈ 3.0
        @test y_hat[1, end - 1] ≈ 6.0
        @test y_hat[end, end] ≈ 2.0

        run_JET_tests(layer, x, ps, st)
    end

    @testset "Variable BitWidth Parameters" begin
        # https://github.com/FluxML/Flux.jl/issues/1421
        layer = Conv((5, 5), 10 => 20, identity; init_weight=Base.randn)
        ps, st = Lux.setup(rng, layer)
        @test ps.bias isa Array{Float64, 4}
    end

    @testset "Depthwise Conv" begin
        x = randn(rng, Float32, 4, 4, 3, 2)

        layer = Conv((2, 2), 3 => 15; groups=3)
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1], 3) == 15
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        layer = Conv((2, 2), 3 => 9; groups=3)
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1], 3) == 9
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        layer = Conv((2, 2), 3 => 9; groups=3, bias=false)
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1], 3) == 9
        run_JET_tests(layer, x, ps, st)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        # Test that we cannot ask for non-integer multiplication factors
        layer = Conv((2, 2), 3 => 10; groups=3)
        @test_throws AssertionError Lux.setup(rng, layer)
    end

    @testset "Conv SamePad kernelsize $k" for k in ((1,), (2,), (3,), (2, 3), (1, 2, 3))
        x = ones(Float32, (k .+ 3)..., 1, 1)

        layer = Conv(k, 1 => 1; pad=Lux.SamePad())
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1]) == size(x)
        run_JET_tests(layer, x, ps, st; call_broken=length(k) == 1)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        layer = Conv(k, 1 => 1; pad=Lux.SamePad(), dilation=k .÷ 2)
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1]) == size(x)
        run_JET_tests(layer, x, ps, st; call_broken=length(k) == 1)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)

        stride = 3
        layer = Conv(k, 1 => 1; pad=Lux.SamePad(), stride=stride)
        ps, st = Lux.setup(rng, layer)

        @test size(layer(x, ps, st)[1])[1:(end - 2)] == cld.(size(x)[1:(end - 2)], stride)
        run_JET_tests(layer, x, ps, st; call_broken=length(k) == 1)
        test_gradient_correctness_fdm((x, ps) -> sum(layer(x, ps, st)[1]), x, ps;
                                      atol=1.0f-3, rtol=1.0f-3)
    end

    @testset "Conv with non quadratic window #700" begin
        x = zeros(Float32, 7, 7, 1, 1)
        x[4, 4, 1, 1] = 1

        layer = Conv((3, 3), 1 => 1)
        ps, st = Lux.setup(rng, layer)

        y = zeros(eltype(ps.weight), 5, 5, 1, 1)
        y[2:(end - 1), 2:(end - 1), 1, 1] = ps.weight
        @test y ≈ layer(x, ps, st)[1]
        run_JET_tests(layer, x, ps, st)

        layer = Conv((3, 1), 1 => 1)
        ps, st = Lux.setup(rng, layer)

        y = zeros(eltype(ps.weight), 5, 7, 1, 1)
        y[2:(end - 1), 4, 1, 1] = ps.weight
        @test y ≈ layer(x, ps, st)[1]
        run_JET_tests(layer, x, ps, st)

        layer = Conv((1, 3), 1 => 1)
        ps, st = Lux.setup(rng, layer)

        y = zeros(eltype(ps.weight), 7, 5, 1, 1)
        y[4, 2:(end - 1), 1, 1] = ps.weight
        @test y ≈ layer(x, ps, st)[1]
        run_JET_tests(layer, x, ps, st)
    end
end


text:
module TreeView

using LightGraphs, TikzGraphs
using MacroTools
using CommonSubexpressions

export LabelledTree, walk_tree, walk_tree!, draw, @tree, @tree_with_call,
        tikz_representation

export make_dag, @dag, @dag_cse

abstract type LabelledDiGraph
end

struct LabelledTree <: LabelledDiGraph
    g::DiGraph
    labels::Vector{Any}
end

add_numbered_vertex!(g) = (add_vertex!(g); top = nv(g))  # returns the number of the new vertex

include("tree.jl")
include("dag.jl")
include("display.jl")

end # module


text:
#=ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер.
 Робот - в произвольной клетке поля.
   РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.=#
function find_marker(r)
    steps=1
    while ismarker(r)==false
        up_left_moves!(r,Nord,steps)
        steps+=1
        down_right_moves!(r,Sud,steps)
        steps+=1
    end
end
function up_left_moves!(r::Robot,Nord::HorizonSide,steps::Int)
    for _ in 1:steps
        if ismarker(r)==false
            move!(r,Nord)
        else
            break
        end
    end
    for _ in 1:steps
        if ismarker(r)==false
             move!(r,West)
        end
    end
end
function down_right_moves!(r::Robot,Sud::HorizonSide,steps::Int)
    for _ in 1:steps
        if ismarker(r)==false
            move!(r,Sud)
        else
            break
        end
    end
    for _ in 1:steps
        if ismarker(r)==false
             move!(r,Ost)
        end
    end
end
moves!(r::Robot, side::HorizonSide, num_steps::Int) =
for _ in 1:num_steps
    move!(r,side)
end

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
using Knet.Layers21: Embed, Linear, Block, 𝑵
using Knet.Ops20: RNN

@isdefined(trmor_data) || include("data.jl")

surface_forms = first.(collect(trmor_data))
surface_vocab = trmor_data.src.vocab
rperm = randperm(MersenneTwister(1234), length(surface_forms))
tst = surface_forms[rperm[1:80]]
val = surface_forms[rperm[81:160]]
trn = surface_forms[rperm[161:end]]

# Best model: (i = 10745, trn = 1.2661633f0, val = 1.4611613f0, tst = 1.4005132f0)
# Saved in charlm.jld2: could look at hidden states at morpheme boundaries

function charlm_model(;
                      vocab=length(surface_vocab.i2w),
                      hidden=350,
                      embed=64,
                      init1=1,
                      drop1=0.2,
                      drop2=0.3,
                      o...)
    Block(Embed(vocab,embed; init=𝑵(init1)),      # B,1+T -> E,B,1+T
          x->dropout(x,drop1),
          RNN(embed,hidden; o...), # -> H,B,T+1
          x->dropout(x,drop2),
          Linear(hidden,vocab),    # -> V,B,T+1
          )
end

charlm_input(x) = hcat(fill(surface_vocab.eos, size(x,1)), x)
charlm_output(x) = hcat(x, fill(surface_vocab.eos, size(x,1)))
charlm_loss(model, batch::Matrix{Int64}) = nll(model(charlm_input(batch)), charlm_output(batch))

function charlm_loss(model, data::Vector{Matrix{Int64}}; average=true)
    sum = cnt = 0
    for batch in data
        (z,n) = nll(model(charlm_input(batch)), charlm_output(batch); average=false)
        sum += z; cnt += n
    end
    average ? sum / cnt : (sum, cnt)
end

function train(model; epochs=35)
    progress!(adam(x->charlm_loss(model,x), ncycle(tuple.(trn),epochs))) do p
        println()
        (trn=charlm_loss(model,trn), val=charlm_loss(model,val), tst=charlm_loss(model,tst))
    end
end

function fasttrain(model; epochs=30)
    adam!(x->charlm_loss(model,x), ncycle(tuple.(trn), epochs))
    charlm_loss(model,trn)
end

# goldensection(optimize_hidden, 1; verbose=true)

function optimize_hidden(logh)  # around 350 gives best train loss in 30 epochs ~1.235, gru around ~1.25
    @show hidden = round(Int, 128*exp(logh[1]/3))
    if hidden > 1000; return NaN; end
    model = charlm_model(; hidden)
    fasttrain(model)
end

function optimize_init1(x)      # doesn't make much difference embed init N(1) is ok
    @show init1 = exp(x[1])
    model = charlm_model(; init1)
    fasttrain(model)
end

function optimize_hidden2(logh) # 2-layer rnn gives optimum around 256 hidden with 1.25 training loss
    @show hidden = round(Int, 128*exp(logh[1]/3))
    if hidden > 1000; return NaN; end
    model = charlm_model(; hidden, numLayers=2)
    fasttrain(model)
end

function optimize_embed(logh)   # embed=64 ok
    @show embed = round(Int, 64*exp(logh[1]/3))
    if embed > 1000; return NaN; end
    model = charlm_model(; embed)
    fasttrain(model)
end

function fastval(model; epochs=30)
    adam!(x->charlm_loss(model,x), ncycle(tuple.(trn), epochs))
    charlm_loss(model,val)
end

function optimize_drop(logp)    # drop1=0.2 drop2=0.3 is optimal with 1.26 trn, 1.46 val, 1.40 tst loss
    @show drop1,drop2 = sigm.(logp)
    model = charlm_model(; drop1, drop2)
    fastval(model)
end

# goldensection(optimize_drop, 2; verbose=true)

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
#Heatmap plots

function PlotHeatmap(label_id, x::AbstractArray{T}, rows::Integer, cols::Integer, scale_min::Real = 0., scale_max::Real = 0., label_fmt::String = "&.1f", bounds_min = LibCImPlot.ImPlotPoint(0.0,0.0),
                 bounds_max = LibCImPlot.ImPlotPoint(1.0,1.0)) where {T<:ImPlotData}
     LibCImPlot.PlotHeatmap(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotHeatmap(label_id, x::AbstractArray{T}, rows::Integer, cols::Integer, scale_min::Real = 0., scale_max::Real = 0., label_fmt::String = "&.1f", bounds_min = LibCImPlot.ImPlotPoint(0.0,0.0),
                 bounds_max = LibCImPlot.ImPlotPoint(1.0,1.0)) where {T<:Float64}
    LibCImPlot.PlotHeatmap(label_id, Float64.(x), rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotHeatmap(x::AbstractArray{T}, rows, cols, scale_min = 0.0, scale_max = 1.0;
                 label_id::String = "", label_fmt="%.1f", bounds_min = LibCImPlot.ImPlotPoint(0.0,0.0),
                 bounds_max = LibCImPlot.ImPlotPoint(1.0,1.0)) where {T<:ImPlotData}

    LibCImPlot.PlotHeatmap(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end


text:
struct TransmissionNetwork <: AbstractTN
  external::BitArray{1}
  internal::BitArray{2}

  function TransmissionNetwork(individuals::Integer)
    return new(fill(0, individuals),
               fill(0, (individuals, individuals)))
  end

  function TransmissionNetwork(external::BitArray{1},
                               internal::BitArray{2})
    if !(length(external) == size(internal, 1) == size(internal, 2))
      error("Mismatched BitArray sizes")
    end
    multiple_exposures = findall((sum(internal, dims=1)[:] .+ external) .> 1)
    if length(multiple_exposures) > 0
      error("Multiple exposures detected for individual(s): $multiple_exposures")
    end
    return new(external, internal)
  end

  function TransmissionNetwork(starting_states::DiseaseStates)
    return new(starting_states .!= Ref(State_S), fill(0, (length(starting_states), length(starting_states))))
  end
end

function innerjoin(x::TN, y::TN) where {TN <: TransmissionNetwork}
  return TN(x.external .* y.external, x.internal .* y.internal)
end

function leftantijoin(x::TN, y::TN) where {TN <: TransmissionNetwork}
  return TN(x.external .* .!y.external, x.internal .* .!y.internal)
end


text:
function MSDE.dZ!(u, dz::Tuple{Int64, Float64, <:SVector}, Z, P::MSDE.SDEKernel)
    i = u[1]
    (Z[i+1][1] - Z[i][1], Z[i+1][2] - Z[i][2],  Z[i+1][3] -  Z[i][3])
end

wrapdiagonal(x) = x
wrapdiagonal(x::Vector) = Diagonal(x)

function MSDE.tangent!(du::Tuple{Int64, Float64, <:SVector}, u, dz, P::MSDE.SDEKernel)
    MSDE.@unpack f, g, p, noise_rate_prototype = P
    k1 = f(u[3],p,u[2])
    g1 = wrapdiagonal(g(u[3],p,u[2]))
    du3 = k1*dz[2] + g1*dz[3]
    (dz[1], dz[2], du3)
end

function MSDE.exponential_map!(u::Tuple{Int64, Float64, <:SVector}, du, P::MSDE.SDEKernel)
    (u[1] + du[1], u[2] + du[2], u[3] + du[3])
end

function myNoiseGrid(t,W,Z=nothing;reset=true)
  val = W[1]
  curt = t[1]
  dt = t[1]
  curW = copy(val)
  dW = copy(val)
  if Z==nothing
    curZ = nothing
    dZ = nothing
  else
    curZ = copy(Z[1])
    dZ = copy(Z[1])
  end
  DiffEqNoiseProcess.NoiseGrid{typeof(val),ndims(val),typeof(dt),typeof(dW),typeof(dZ),typeof(Z),false}(
            t,W,W,Z,curt,curW,curZ,dt,dW,dZ,true,reset)
end

"""
  myinnov(t)

t:  either an array of Numbers or StepRangeLen{Number}
returns NoiseGrid, with a Wiener process on the grid t.
fac: factor that allows to scale the noise values (e.g., set them to zero)
"""
function myinnov(t, 𝕏_ ,fac=1)
    dt = diff(t)
    w = [fac*sqrt(dt[i])*randn(𝕏_) for i in 1:length(t)-1]
    brownian_values = cumsum(pushfirst!(w, zero(𝕏_)))
    myNoiseGrid(t,brownian_values)
end

"""
  pcn_innov(Z, ρ)

update NoiseGrid Z by pCN-step with parameter ρ. Taking ρ=1 just returns a
NoiseGrid with the same values as in Z
"""
function pcn_innov(Z, ρ, 𝕏_, fac=1)
    Znew = myinnov(Z.t, 𝕏_, fac)
    a = cumsum(pushfirst!(ρ * diff(Z.W) + sqrt(1. - ρ^2) * diff(Znew.W), zero(𝕏_)))
    myNoiseGrid(Z.t, a)
end

"""
    forwardsample(tree::Tree, rootval::𝕏, θ, dt0, f, g)

Forward sample a diffusiion process with (drift, diffusivity) = (f,g) and parameter θ on a grid with target mesh width dt0.
Starting value is rootval.

Returns Xd::Vector{𝕏} and segs, which is a vector of ODE-sols. Note that segs[1] should be skipped (not defined)
"""
#function forwardsample(tree::Tree, rootval::𝕏, θ, dt0, f, g)
function forwardsample(tree::Tree, rootval, θ, dt0, f, g, SDEalg)
    Xd = [rootval]    # save endpoints
    segs = Vector{Any}(undef, tree.n) # save all guided segments
    for i in eachindex(tree.T)
        #local κ, x, xT, dt, u
        i == 1 && continue  # skip root-node  (has no parent)

        # define timegrid on segment
        i′ = tree.Par[i]
        t = tree.T[i′]
        m = round(Int, (tree.T[i]-t)/dt0)
        dt = (tree.T[i] - t)/m

        u = Xd[i′]   # starting value
        # set noise_rate_prototype for comparison with StochasticDiffEq package
        κ = MSDE.SDEKernel(f, g, t:dt:tree.T[i], θ, Diagonal(θ[2]))
        # choose SDEalg =  MSDE.EulerMaruyama!() for new fast solver...
        NG = myinnov(t:dt:tree.T[i], 𝕏)
        xT, seg  = MSDE.sample(κ, u, SDEalg, NG)
        push!(Xd, xT)
        segs[i] = seg
    end
    Xd, segs
end

"""
    bwfiltertree!(Q, tree::Tree, θlin, dt0; apply_time_change=false)

Backward filter on the tree, using the lienar process specified by θlin, targetting mesh-width dt0
Q must have been initialised on all leaves using a WGaussian

Returns updated Q and messages
"""
function bwfiltertree!(Q, tree::Tree, θlin, dt0; apply_time_change=false, alg=Tsit5())
    T = tree.T
    messages = Vector{Any}(undef, tree.n)
    for i in reverse(eachindex(T))
        i == 1 && continue  # skip root node  (has no parent)
        ipar = tree.Par[i]

        δ = T[i]-T[ipar]
        dt = δ/round(Int, δ/dt0)
        if apply_time_change
            tvals = MSDE.timechange(T[ipar]:dt:T[i])
        else
            tvals = T[ipar]:dt:T[i]
        end

        κ̃ = MSDE.SDEKernel(MS.AffineMap(θlin[1], θlin[2]), MS.ConstantMap(θlin[3]), tvals, θlin)

        # message, u = MSDE.backwardfilter(κ̃, Q[i], alg=OrdinaryDiffEq.Tsit5(), abstol=1e-12, reltol=1e-12, apply_timechange=apply_time_change)
        message, u = MSDE.backwardfilter(κ̃, Q[i], apply_timechange=apply_time_change, alg=alg)
        messages[i] = message
        if tree.lastone[i] # last child is encountered first backwards
            Q[ipar] = u
        else
            Q[ipar] = MS.fuse(u, Q[ipar]; unfused=false)[2]
        end
    end
    Q, messages
end

"""
    fwguidtree!(X, guidedsegs, messages, tree::Tree, f, g, θ, Z; apply_time_change=false)

Returns updated `X` and `guidedsegs`, as also the loglikelihood vector `ll` (at all vertices of the tree) and
`𝐋` the loglikelihood summed over all leaf-indices.
"""

function fwguidtree!(X, guidedsegs, Q, messages, tree::Tree, f, g, θ, Z, SDEalg; apply_time_change=false)
    ll = zeros(tree.n)
    for i in eachindex(tree.T)
        i == 1 && continue  # skip root-node (has no parent)
        κ = MSDE.SDEKernel(f, g, messages[i].ts, θ, zeros(d,d))
        ipar = tree.Par[i]
        (solend, llnew), res = MSDE.forwardguiding(κ, messages[i], (X[ipar], 0.0), SDEalg, Z[i-1],
                                                            inplace=false, apply_timechange=apply_time_change)
        ll[i] = llnew + ll[ipar] * tree.lastone[i]
        X[i] = solend
        guidedsegs[i] = res
    end
    𝐋 = sum(ll[tree.lids]) + logdensity(Q[1], X[1]) #logdensity(convert(WGaussian{(:F,:Γ,:c)},Q[1]), X[1])
    X, guidedsegs, ll, 𝐋
end


text:
module DemoCards

import Base: basename
using Mustache
using ImageCore
using FileIO
using JSON
using YAML

const config_filename = "config.json"
# directly copy these folders without processing
const ignored_dirnames = ["assets"]

include("compat.jl")

include("types/card.jl")
include("types/section.jl")
include("types/page.jl")

include("utils.jl")
include("show.jl")
include("generate.jl")
include("cardthemes.jl")

export makedemos, cardtheme

"""

This package breaks the rendering of the whole demo page into three types:

* `DemoPage` contains serveral `DemoSection`s;
* `DemoSection` contains either serveral `DemoSection`s or serveral `DemoCard`s;
* `DemoCard` consists of cover image, title and other necessary information.
"""
DemoCards

end # module


text:
#=
---------------------------------------------------------------------
This function computes the element mass and Laplacian matrices.

The strategy follows Algorithm 12.7 in F.X. Giraldo's Introduction to Element-based
Galerkin Methods using Tensor-Product Bases: Analysis, Algorithms, and Applications.

Written by F.X. Giraldo on 5/2019
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function element_matrices(ψ,dψ,ξ_x,ξ_y,η_x,η_y,jac,ωq,Ne,Np,Nq,DFloat)

    #Initialize Matrices
    M=zeros(DFloat,Np,Np,Np,Np,Ne)
    L=zeros(DFloat,Np,Np,Np,Np,Ne)

    #Construct Mass and Differentiation Matrices
    for e=1:Ne
        for l=1:Nq, k=1:Nq
            wq=ωq[k]*ωq[l]*jac[k,l,e]
            for j=1:Np, i=1:Np
                Ψ_JK=ψ[i,k]*ψ[j,l] #h_ik*h_jl
                dΨdx_JK=dψ[i,k]*ψ[j,l]*ξ_x[k,l,e] + ψ[i,k]*dψ[j,l]*η_x[k,l,e]
                dΨdy_JK=dψ[i,k]*ψ[j,l]*ξ_y[k,l,e] + ψ[i,k]*dψ[j,l]*η_y[k,l,e]
                for n=1:Np, m=1:Np
                    Ψ_IK=ψ[m,k]*ψ[n,l] #h_ik*h_jl
                    dΨdx_IK=dψ[m,k]*ψ[n,l]*ξ_x[k,l,e] + ψ[m,k]*dψ[n,l]*η_x[k,l,e]
                    dΨdy_IK=dψ[m,k]*ψ[n,l]*ξ_y[k,l,e] + ψ[m,k]*dψ[n,l]*η_y[k,l,e]
                    M[m,n,i,j,e]+=wq*Ψ_IK*Ψ_JK
                    L[m,n,i,j,e]+=wq*( dΨdx_IK*dΨdx_JK + dΨdy_IK*dΨdy_JK )
                end #m,n
            end #j,i
        end #k,l
    end #e

    return (M,L)
end


text:
abstract type AbstractProblem{dType <: Number, tType <: Real, arrayType <: AbstractArray{dType}} end

"Returns the parent equation object of the problem."
equation(prob::AbstractProblem) = error("equation() not implemented for ", typeof(prob), ".")

"Returns a NamedTuple containing all functions (e.g. vector fields) provided by the equation."
functions(prob::AbstractProblem) = functions(equation(prob))

"Returns a NamedTuple containing all solutions provided by the equation."
solutions(prob::AbstractProblem) = solutions(equation(prob))

"Returns a NamedTuple containing all invariants provided by the equation."
invariants(prob::AbstractProblem) = invariants(equation(prob))

tspan(prob::AbstractProblem) = error("tspan() not implemented for ", typeof(prob), ".")
tstep(prob::AbstractProblem) = error("tstep() not implemented for ", typeof(prob), ".")
tbegin(prob::AbstractProblem) = tspan(prob)[begin]
tend(prob::AbstractProblem) = tspan(prob)[end]

GeometricBase.parameters(prob::AbstractProblem) = error("parameters() not implemented for ", typeof(prob), ".")
GeometricBase.periodicity(prob::AbstractProblem) = periodicity(equation(prob))

hassolution(prob::AbstractProblem) = hassolution(equation(prob))
hasvectorfield(prob::AbstractProblem) = hasvectorfield(equation(prob))
hasprimary(prob::AbstractProblem) = hasprimary(equation(prob))
hassecondary(prob::AbstractProblem) = hassecondary(equation(prob))

hasinvariants(prob::AbstractProblem) = hasinvariants(equation(prob))
hasparameters(prob::AbstractProblem) = hasparameters(equation(prob))
hasperiodicity(prob::AbstractProblem) = hasperiodicity(equation(prob))

hashamiltonian(prob::AbstractProblem) = hashamiltonian(equation(prob))
haslagrangian(prob::AbstractProblem) = haslagrangian(equation(prob))


text:
"""Test if a value is numeric."""
is_numeric(value) = isa(value, Real)


text:
using NLSolvers, Test

function fourth(∇²f, ∇f, x)
    if !(∇²f == nothing)
        ∇²f = 12x^2 - sin(x)
    end
    if !(∇f == nothing)
        ∇f = 4x^3 + cos(x)
    end

    fx = x^4 + sin(x)
    if ∇f == nothing && ∇²f == nothing
        return fx
    elseif ∇²f == nothing
        return fx, ∇f
    else
        return fx, ∇f, ∇²f
    end
end

@testset "scalar no-alloc" begin
scalar_obj = TwiceDiffed(fourth; infer=true)
_alloc = @allocated minimize(scalar_obj, 4.0, SR1(Direct()))
_alloc = @allocated minimize(scalar_obj, 4.0, SR1(Direct()))
@test _alloc == 0

_alloc = @allocated minimize(scalar_obj, 4.0, BFGS(Direct()))
_alloc = @allocated minimize(scalar_obj, 4.0, BFGS(Direct()))
@test _alloc == 0

_alloc = @allocated minimize(scalar_obj, 4.0, DFP(Direct()))
_alloc = @allocated minimize(scalar_obj, 4.0, DFP(Direct()))
@test _alloc == 0

_alloc = @allocated minimize(scalar_obj, 4.0, Newton())
_alloc = @allocated minimize(scalar_obj, 4.0, Newton())
@test _alloc == 0

_alloc = @allocated minimize(scalar_obj, 4.0, Newton())
_alloc = @allocated minimize(scalar_obj, 4.0, Newton())
@test _alloc == 0
end


text:
module JuliaAPI

using Compat
using Compat.REPL
using Compat.Dates
using Compat.Pkg

import PyCall
using PyCall: PyObject, pyjlwrap_new

using Requires

function eval_str(code::AbstractString;
                  scope::Module = Main,
                  filename::AbstractString = "string",
                  auto_jlwrap = true,
                  force_jlwrap = false)
    result = include_string(scope, code, filename)
    if force_jlwrap
        return pyjlwrap_new(result)
    elseif auto_jlwrap
        return _wrap(result)
    end
    return result
end

const NpyNumber = Union{
    (t for t in values(PyCall.npy_typestrs) if t <: Number)...
}

_wrap(obj::Any) = pyjlwrap_new(obj)
# Wrap the object if it is desirable to _not_ invoke PyCall's
# automatic conversion.

_wrap(obj::Union{
    Nothing,
    Integer,
    NpyNumber,            # to exclude ForwardDiff.Dual etc.
    Array{<: NpyNumber},  # ditto
    AbstractString,  # should it be just String?
    Dates.AbstractTime,
    IO,
}) = obj
# It's OK to include some types that are not supported PyCall.  In
# that case, those objects are passed through pyjlwrap_new anyway.
# What should be avoided here is to include some types that would be
# converted by PyCall in an irreversible manner (e.g., Symbol,
# BitArray, etc.).

"""
    wrapcall(f, args...; kwargs...)

Wrap what `f(args...; kwargs...)` returns.
"""
wrapcall(f::Base.Callable, args...; kwargs...) = _wrap(f(args...; kwargs...))

struct WrappingCallback
    o::PyObject
    force::Bool
end

function (f::WrappingCallback)(args...; kwargs...)
    wrap = f.force ? pyjlwrap_new : _wrap
    return f.o(wrap.(args)...; (k => wrap(v) for (k, v) in kwargs)...)
end

set_var(name::String, value) = set_var(Symbol(name), value)

function set_var(name::Symbol, value)
    Base.eval(Main, :($name = $value))
    nothing
end

@static if VERSION < v"0.7-"
    getproperty_str(obj, name) = getfield(obj, Symbol(name))
else
    getproperty_str(obj, name) = getproperty(obj, Symbol(name))
end

@static if VERSION < v"0.7-"
    dir(m::Module; all=true, imported=false) = String.(names(m, all, imported))
    dir(::T; _...) where T = String.(fieldnames(T))
else
    dir(m::Module; kwargs...) = String.(names(m; all=true, kwargs...))
    dir(m; all=true) = String.(propertynames(m, all))
end

struct _jlwrap_type end  # a type that would be wrapped as jlwrap by PyCall

get_jlwrap_prototype() = _jlwrap_type()

pybroadcast(op::String, args...) =
    _pybroadcast(eval(Meta.parse(strip(op))), args...)
pybroadcast(op, args...) = _pybroadcast(op, args...)
_pybroadcast(op, args::AbstractArray...) = op.(args...)
_pybroadcast(op, args...) = op(args...)

struct PySlice
    start
    stop
    step
end

struct PyEllipsis end

pygetindex(x, indices...) =
    getindex(x, convert_pyindices(x, indices)...)
pysetindex!(x, value, indices...) =
    setindex!(x, value, convert_pyindices(x, indices)...)

function convert_pyindices(x, indices)
    if applicable(firstindex, x, 1)
        # arrays
        return convert_pyindex.((x,),
                                process_pyindices(ndims(x), indices),
                                1:length(indices))
    elseif indices isa Tuple{Integer} && applicable(firstindex, x)
        # tuples, etc.
        return (firstindex(x) + indices[1],)
    else
        # dictionaries, etc.
        return indices
    end
end

function process_pyindices(nd, indices)
    colons = repeat([:], inner=nd - length(indices))
    i = findfirst(i -> i isa PyEllipsis, indices)
    if i === nothing
        return (indices..., colons...)
    else
        return (indices[1:i - 1]..., colons..., ndices[i + 1:end]...)
    end
end

convert_pyindex(x, i, ::Int) = i

convert_pyindex(x, i::Integer, d::Int) = firstindex(x, d) + i

function convert_pyindex(x, slice::PySlice, d::Int)
    start = (slice.start === nothing ? 0 : slice.start) + firstindex(x, d)
    stop = slice.stop === nothing ? lastindex(x, d) : slice.stop + 1
    step = slice.step === nothing ? 1 : slice.step
    if step == 1
        return start:stop
    else
        return start:step:stop
    end
end

@static if VERSION < v"0.7-"
    completions(_a...; __k...) = String[]
else
    function completions(string, pos, context_module = Main)
        ret, _, should_complete =
            REPL.completions(string, pos, context_module)
        if should_complete
            return map(REPL.completion_text, ret)
        else
            return String[]
        end
    end
end

"""
    start_repl(; <keyword arguments>)

Start Julia REPL.

# Keyword Arguments
- `interactive::Bool = true`
- `quiet::Bool = true`
- `banner::Bool = false`
- `history_file::Bool = true`
- `color_set::Bool = false`: "color (configuration is already) set"
"""
function start_repl(;
        interactive::Bool = true,
        quiet::Bool = true,
        banner::Bool = false,
        history_file::Bool = true,
        color_set::Bool = false,
        )
    if !(stdout isa Base.TTY)
        error("stdout is not a TTY")
    elseif !(stdin isa Base.TTY)
        error("stdin is not a TTY")
    end
    was_interactive = Base.is_interactive
    try
        # Required for Pkg.__init__ to setup the REPL mode:
        Base.eval(:(is_interactive = $interactive))

        Base.run_main_repl(interactive, quiet, banner, history_file, color_set)
    finally
        Base.eval(:(is_interactive = $was_interactive))
    end
end

function __init__()
    # Don't wrap JuliaPy wrappers
    @require Pandas="eadc2687-ae89-51f9-a5d9-86b5a6373a9c" @eval begin
        _wrap(obj::Pandas.PandasWrapped) = obj
    end
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" @eval begin
        _wrap(obj::SymPy.SymbolicObject) = obj
    end
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" @eval begin
        _wrap(obj::PyPlot.Figure) = obj
    end
end

module IOPiper

const orig_stdin  = Ref{IO}()
const orig_stdout = Ref{IO}()
const orig_stderr = Ref{IO}()

function __init__()
    orig_stdin[]  = stdin
    orig_stdout[] = stdout
    orig_stderr[] = stderr
end

"""
If `d` ends with an incomplete UTF8-encoded character, return the number of trailing incomplete bytes.
Otherwise, return `0`.

Taken from IJulia.jl.
"""
function num_utf8_trailing(d::Vector{UInt8})
    i = length(d)
    # find last non-continuation byte in d:
    while i >= 1 && ((d[i] & 0xc0) == 0x80)
        i -= 1
    end
    i < 1 && return 0
    c = d[i]
    # compute number of expected UTF-8 bytes starting at i:
    n = c <= 0x7f ? 1 : c < 0xe0 ? 2 : c < 0xf0 ? 3 : 4
    nend = length(d) + 1 - i # num bytes from i to end
    return nend == n ? 0 : nend
end

function pipe_stream(sender::IO, receiver, buf::IO = IOBuffer())
    try
        while !eof(sender)
            nb = bytesavailable(sender)
            write(buf, read(sender, nb))

            # Taken from IJulia.send_stream:
            d = take!(buf)
            n = num_utf8_trailing(d)
            dextra = d[end-(n-1):end]
            resize!(d, length(d) - n)
            s = String(copy(d))

            write(buf, dextra)
            receiver(s)  # check isvalid(String, s)?
        end
    catch e
        if !isa(e, InterruptException)
            rethrow()
        end
        pipe_stream(sender, receiver, buf)
    end
end

const read_stdout = Ref{Base.PipeEndpoint}()
const read_stderr = Ref{Base.PipeEndpoint}()

function pipe_std_outputs(out_receiver, err_receiver)
    global readout_task
    global readerr_task
    read_stdout[], = redirect_stdout()
    readout_task = @async pipe_stream(read_stdout[], out_receiver)
    read_stderr[], = redirect_stderr()
    readerr_task = @async pipe_stream(read_stderr[], err_receiver)
end

end  # module

end  # module


text:
using JuMP, Ipopt

function Aiyagari.get_optimum(states, agg_state, 𝔼V, params, a_grid)
  @unpack p, r = agg_state
  @unpack β, θ, δ = params

  model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

  MOI.set(model, MOI.RawParameter("print_level"), 0)
  #model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

  register(model, :u, 2, u, autodiff=true)
  register(model, :𝔼V, 1, w -> 𝔼V(w), autodiff=true)

  guess = sum(states)/2

  @variable(model, c >= eps(), start=guess)
  @variable(model, h >= eps(), start=guess/2)

  @NLparameter(model, w == states.a)
  @NLparameter(model, y == states.z)

  w_next = @NLexpression(model, w + y - c - p * h * (r + δ))
  m = @NLexpression(model, p * h + c - y - w)
  @NLconstraint(model, (1+r) * m <= p * h * (1-δ) * θ)
  @NLobjective(model, Max, u(c,h) + β * 𝔼V(w_next) )

  MOI.set(model, MOI.Silent(), true)

  JuMP.optimize!(model)

  val = objective_value(model)

  pol_full = (c=value(c), h=value(h), m=value(m), w_next=value(w_next), conv=termination_status(model))
  pol = pol_full.w_next, pol_full.h

  conv = termination_status(model) == 4

  (pol=pol, pol_full=pol_full, val=val, conv=conv)
end

# constrained optimization is provided in Optim!
#https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/#nonlinear-constrained-optimization


text:
struct CrossValidation
    indices
    partition::Integer
    samples::Array
    folds::Integer
    use_all::Bool
end

function CrossValidation(folds::Integer, samples::AbstractArray;shuffle::Bool=true,use_all::Bool=true)
    ids = LinearIndices(samples) |> collect
    shuffle!(ids)
    partition = length(ids) ÷ folds
    CrossValidation(ids, partition, samples, folds, use_all)
end

function (x::CrossValidation)(fold::Integer;)
    (fold < 1 || fold > x.folds) && throw(ArgumentError("Please check 1 ≤ fold ≤ $(x.folds)"))
    last = x.use_all ? length(x.indices) : x.folds * x.partition
    select_start = (fold - 1) * x.partition + 1
    select_end = fold * x.partition
    train = x.samples[x.indices[1:select_start - 1]]
    append!(train, x.samples[x.indices[select_end + 1:last]])
    test = x.samples[x.indices[select_start:select_end]]
    return (train = train, test = test)
end

struct ValidateResult{T <: AbstractFloat}
    dir::String
    scores::Matrix{T}
    labels::Matrix{Bool}
    mask::Matrix{Bool}
    drugs::Array{String}
    diseases::Array{String}
end

show(io::IO,vr::ValidateResult) = print(io, "Validate Results($(vr.dir))")

function ValidateResult(dir::String)
    labels = readdlm(joinpath(dir, "RD.csv"), ',', Bool)
    drugs = readdlm(joinpath(dir, "drugs.csv"), ',', String)
    diseases = readdlm(joinpath(dir, "diseases.csv"), ',', String)
    mask = zeros(size(labels))
    scores = zeros(Float32, size(labels))
    folds = map(x->joinpath(abspath(dir),x),readdir(dir))
    filter!(isdir,folds)
    for f in folds
        mask = mask .+ readdlm(joinpath(dir, f, "mask.csv"), ',', Bool)
        scores = scores .+ readdlm(joinpath(dir, f, "result.csv"), ',', Float32)
    end
    scores = scores ./ mask
    ValidateResult(dir, scores, labels, map(>(0), mask), drugs, diseases)
end

function analyze(res::ValidateResult;union_all::Bool=false)
    if union_all
        cms = threshold_scan(res.labels[res.mask], res.scores[res.mask])
    else
        cms_all = []
        for (name, l, s, m) in zip(res.drugs, eachrow(res.labels), eachrow(res.scores), eachrow(res.mask))
            sum(m)==0 && continue
            tmp = threshold_scan(l[m], s[m])
            push!(cms_all, tmp)
        end
        cms=[]
        for i=1:length(cms_all[1])
            cm=ConfusionMatrix(zeros(eltype(cms_all[1][1]),size(cms_all[1][1])))
            for each in cms_all
                length(each)<i && continue
                cm=cm+each[i]
            end
            push!(cms,cm)
        end
    end
    cms
end

function roc(r::ValidateResult;union_all::Bool=false)
    cms=analyze(r,union_all=union_all)
    TPR = tpr.(cms, 2)
    FPR = fpr.(cms, 2)
    FPR,TPR
end

function pr(r::ValidateResult,union_all::Bool=false)
    cms=analyze(r,union_all=union_all)
    TPR = tpr.(cms, 2)
    FPR = fpr.(cms, 2)
    FPR,TPR
end

function analyze_by_drugs(result::ValidateResult)
    df = DataFrame("Drug Name" => String[], "AUPR" => Float64[], "AUROC" => Float64[])
    for (name, l, s, m) in zip(result.drugs, eachrow(result.labels), eachrow(result.scores), eachrow(result.mask))
        sum(m)==0 && continue
        cms = threshold_scan(l[m], s[m])
        TPR = tpr.(cms, 2)
        FPR = fpr.(cms, 2)
        P = precision.(cms, 2)
        R = recall.(cms, 2)
        aupr = auc(R, P)
        auroc = auc(FPR, TPR)
        push!(df, (name, aupr, auroc))
    end
    push!(df, ("Average of $(size(df, 1)) drugs", mean(df.AUPR), mean(df.AUROC)))
    df
end

function pr_roc(res::ValidateResult;union_all::Bool=false)
    cms=analyze(res;union_all=union_all)
    TPR = tpr.(cms, 2)
    FPR = fpr.(cms, 2)
    P = precision.(cms, 2)
    R = recall.(cms, 2)
    R,P,FPR,TPR
end

function plot(res::ValidateResult;union_all::Bool=false,)
    R,P,FPR,TPR=pr_roc(res,union_all=union_all)
    aupr = auc(R, P)
    auroc = auc(FPR, TPR)
    left = plot(R, P, title="P-R(auc=$(round(aupr, digits=3)))")
    right = plot(FPR, TPR, title="ROC(auc=$(round(auroc, digits=3)))")
    plot(left, right, layout=(1, 2)) |> display
end

export CrossValidation,ValidateResult,analyze_by_drugs,analyze,pr,roc,pr_roc

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
include("defs.jl")

bd = billiard_sinai(0.25)
ps = particlebeam(0.2, 0.75, π + π/4 + 0.414235, 100, 0.002)

fig, = interactive_billiard(
    bd, ps;
    tail = 500, backgroundcolor = RGBf0(1,1,1),
    vr = 0.05
)

# record_interaction(fig, joinpath(@__DIR__, "sinai_billiard.mp4"))

# %% Ergodicity in the Sinai billiard
bd = billiard_sinai(0.25)

fig, = interactive_billiard_bmap(bd)

# record_interaction(fig, joinpath(@__DIR__, "sinai_ergodicity.mp4"); total_time = 20)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Revise
using PyPlot
using LinearAlgebra
using Infiltrator
using Bem2d

function stylesubplots_local()
    gca().set_aspect("equal")
    gca().set_xticks([])
    gca().set_yticks([])
    return nothing
end

function plot18(els, x, y, disp1, stress1, string1, disp2, stress2, string2, title_string)
    # Set contour levels for displacements and stresses
    contourvecdispy = collect(LinRange(-5e-11, 5e-11, 51))
    contourvecdispx = collect(LinRange(-1e-12, 1e-12, 51))
    contourvecstress = collect(LinRange(-1e-2, 1e-2, 51))
    cmap = PyPlot.get_cmap("seismic")
    fontsize = 30
    PyPlot.figure(figsize=(30, 20))

    # Fields from first model
    PyPlot.subplot(3, 6, 1)
    PyPlot.quiver(x[:], y[:], disp1[:, 1], disp1[:, 2], units="width", color="b")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("u", fontsize=fontsize)

    PyPlot.subplot(3, 6, 2)
    PyPlot.contourf(x, y, reshape(disp1[:, 1], size(x)), contourvecdispx, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp1[:, 1], size(x)), contourvecdispx, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("ux (" * string1 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 3)
    PyPlot.contourf(x, y, reshape(disp1[:, 2], size(x)), contourvecdispy, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp1[:, 2], size(x)), contourvecdispy, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("uy (" * string1 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 4)
    PyPlot.contourf(x, y, reshape(stress1[:, 1], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 1], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxx (" * string1 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 5)
    PyPlot.contourf(x, y, reshape(stress1[:, 2], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 2], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("syy (" * string1 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 6)
    PyPlot.contourf(x, y, reshape(stress1[:, 3], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 3], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxy (" * string1 * ")", fontsize=fontsize)

    # Fields from second model
    PyPlot.subplot(3, 6, 7)
    PyPlot.quiver(x[:], y[:], disp2[:, 1], disp2[:, 2], units="width", color="b")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("u", fontsize=fontsize)

    PyPlot.subplot(3, 6, 8)
    PyPlot.contourf(x, y, reshape(disp2[:, 1], size(x)), contourvecdispx, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp2[:, 1], size(x)), contourvecdispx, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("ux (" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 9)
    PyPlot.contourf(x, y, reshape(disp2[:, 2], size(x)), contourvecdispy, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp2[:, 2], size(x)), contourvecdispy, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("uy (" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 10)
    PyPlot.contourf(x, y, reshape(stress2[:, 1], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress2[:, 1], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxx (" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 11)
    PyPlot.contourf(x, y, reshape(stress2[:, 2], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress2[:, 2], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("syy (" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 12)
    PyPlot.contourf(x, y, reshape(stress2[:, 3], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress2[:, 3], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxy (" * string2 * ")", fontsize=fontsize)

    # Residuals
    PyPlot.subplot(3, 6, 13)
    PyPlot.quiver(x[:], y[:], disp1[:, 1]-disp2[:, 1], disp1[:, 2]-disp2[:, 2], units="width", color="b")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("u", fontsize=fontsize)

    PyPlot.subplot(3, 6, 14)
    PyPlot.contourf(x, y, reshape(disp1[:, 1]-disp2[:, 1], size(x)), contourvecdispx, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp1[:, 1]-disp2[:, 1], size(x)), contourvecdispx, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("uy (" * string1 * "-" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 15)
    PyPlot.contourf(x, y, reshape(disp1[:, 2]-disp2[:, 2], size(x)), contourvecdispy, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(disp1[:, 2]-disp2[:, 2], size(x)), contourvecdispy, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("uy (" * string1 * "-" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 16)
    PyPlot.contourf(x, y, reshape(stress1[:, 1]-stress2[:, 1], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 1]-stress2[:, 1], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxx (" * string1 * "-" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 17)
    PyPlot.contourf(x, y, reshape(stress1[:, 2]-stress2[:, 2], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 2]-stress2[:, 2], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("syy (" * string1 * "-" * string2 * ")", fontsize=fontsize)

    PyPlot.subplot(3, 6, 18)
    PyPlot.contourf(x, y, reshape(stress1[:, 3]-stress2[:, 3], size(x)), contourvecstress, cmap=cmap)
    PyPlot.colorbar(fraction=0.020, pad=0.05, extend="both")
    PyPlot.contour(x, y, reshape(stress1[:, 3]-stress2[:, 3], size(x)), contourvecstress, linewidths=0.25, colors="k")
    stylesubplots_local()
    Bem2d.plotelements(els)
    PyPlot.title("sxy (" * string1 * "-" * string2 * ")", fontsize=fontsize)

    PyPlot.suptitle(title_string, fontsize=fontsize)
    PyPlot.show()
end

# function local_subplot(x, y, mat, npts, title_string)
#     fontsize = 20
#     contour_levels = 50
#     contour_color = "white"
#     contour_linewidth = 0.5
#     color_scale = 1e6
#     PyPlot.contourf(reshape(x, npts, npts), reshape(y, npts, npts), reshape(mat, npts, npts), levels=contour_levels)
#     cbar = PyPlot.colorbar(fraction=0.020, pad=0.05, extend = "both")
#     PyPlot.contour(reshape(x, npts, npts), reshape(y, npts, npts), reshape(mat, npts, npts), levels=contour_levels, colors=contour_color, linewidths=contour_linewidth)
#     PyPlot.title(title_string, fontsize=fontsize)
#     PyPlot.xlim([-1000, 1000])
#     PyPlot.ylim([-1000, 1000])
#     PyPlot.xticks([])
#     PyPlot.yticks([])
#     PyPlot.gca().set_aspect("equal")
#     PyPlot.gca().tick_params(labelsize=fontsize)
# end

"""
   kelvinUD(x, y, xoffset, yoffset, fx, fy, mu, nu)

   Analytic point source Kelvin solution from Crouch and Starfield
   section 4.2.
"""
function kelvinUD(x, y, xoffset, yoffset, fx, fy, mu, nu)
    x = @. x - xoffset
    y = @. y - yoffset
    Ukelvin = zeros(length(x), 2)
    Skelvin = zeros(length(x), 3)
    C = 1/(4*pi*(1-nu))
    r = @. sqrt(x^2+y^2)
    g = @. -C * log(r)
    gx = @. -C * x/(x^2+y^2)
    gy = @. -C * y/(x^2+y^2)
    gxy = @. C * 2*x*y/(x^2+y^2)^2
    gxx = @. C * (x^2-y^2)/(x^2+y^2)^2
    gyy = -gxx
    Ukelvin[:, 1] = @. fx/(2*mu)*((3-4*nu)*g-x*gx) + fy/(2*mu)*(-y*gx)
    Ukelvin[:, 2] = @. fx/(2*mu)*(-x*gy) + fy/(2*mu)*((3-4*nu)*g-y*gy)
    Skelvin[:, 1] = @. fx*(2*(1-nu)*gx-x*gxx) + fy*(2*nu*gy-y*gxx)
    Skelvin[:, 2] = @. fx*(2*nu*gx-x*gyy) + fy*(2*(1-nu)*gy-y*gyy)
    Skelvin[:, 3] = @. fx*((1-2*nu)*gy-x*gxy) + fy*((1-2*nu)*gx-y*gxy)
    return Ukelvin, Skelvin
end

function kelvin()
    mu = 3e10
    nu = 0.25

    # Observation coordinates
    npts = 50
    obswidth = 5
    x, y = Bem2d.obsgrid(-obswidth, -obswidth, obswidth, obswidth, npts)
    fx = 0.0
    fy = 1.0

    # Offset from origin
    xoffset = 0.0
    yoffset = 0.0

    # Analytic Kelvin point source solution
    Ukelvin, Skelvin = kelvinUD(x, y, xoffset, yoffset, fx, fy, mu, nu)

    # BEM solution
    a = 0.00005
    els = Elements(Int(2))
    els.x1[1] = -a + xoffset
    els.y1[1] = 0.0 + yoffset
    els.x2[1] = a + xoffset
    els.y2[1] = 0.0 + yoffset
    els.name[1] = "point"
    standardize_elements!(els)
    idx = getidxdict(els)
    Ubem, Sbem = constdispstress(trac2dispstress, x, y, els, idx["point"], fx/(2*a), fy/(2*a), mu, nu)

    # 18 panel plot
    close("all")
    fontsize = 20
    xobs = reshape(x, npts, npts)
    yobs = reshape(y, npts, npts)
    plot18(els, xobs, yobs, Ubem, Sbem, "BEM", Ukelvin, Skelvin, "Kelvin point", "fy (BEM vs. Kelvin point)")
    tight_layout()
    show()
end
kelvin()


text:
using GeometryBasics, RPRMakie
using Colors, FileIO, ImageShow
using Colors: N0f8

RPRMakie.activate!()
fig = Figure(; resolution=(1200, 800), fontsize=26)
radiance = 10000
lights = [EnvironmentLight(0.5, load(RPR.assetpath("studio026.exr"))),
          PointLight(Vec3f(0, 0, 20), RGBf(radiance, radiance, radiance))]

ax = LScene(fig[1, 1]; scenekw=(lights=lights,))
Makie.axis3d!(ax.scene, Rect3f(Vec3f(0), Vec3f(1, 1, 1.2)))
screen = RPRScreen(ax.scene; iterations=200, plugin=RPR.Northstar)
matsys = screen.matsys
rectMesh = FRect3D(Vec3f0(-0.5, -0.5, 0), Vec3f0(1))
recmesh = GeometryBasics.normal_mesh(rectMesh)
n = 100
pos = [Point3f(i, j, 0) ./ 10 for i in 1:n for j in 1:n]
z = rand(n, n)
mat = RPR.Chrome(matsys)
mat.roughness = Vec4f(0.3)
meshscatter!(ax, pos; marker=recmesh, markersize=Vec3f.(0.1, 0.1, z[:]), material=mat, color=vec(z))

cam = cameracontrols(ax.scene)
cam.eyeposition[] = Float32[5, 22, 12]
cam.lookat[] = Float32[5, 5, -0.5]
cam.upvector[] = Float32[0.0, 0.0, 1.0]
cam.fov[] = 14.0

# show the screen directly in e.g. VSCode (display unecessary, )
# display(screen)
# alternatively, get the Image directly from the screen

image = colorbuffer(screen)::Matrix{RGB{N0f8}}

save("bars.png", image)


text:
@lexer RLexer let
    keywords = [
        "if", "else", "for", "while", "repeat", "in", "next", "break", "return", "switch",
        "function"
    ]
    types = [
        "array", "category", "character", "complex", "double", "function", "integer",
        "list", "logical", "matrix", "numeric", "vector", "data.frame", "c"
    ]
    builtin_symbols = [
        "NULL", "NA(_(integer|real|complex|character)_)?", "letters", "LETTERS", "Inf",
        "TRUE", "FALSE", "NaN", "pi", "\\.\\.(\\.|[0-9]+)"
    ]
    Dict(
        :name => "R",
        :description => "A lexer for the R language for statistical computing.",
        :aliases => ["r", "rlang"],
        :filenames => ["*.r"],
        :comments => "Based on the Pygments lexer.",
        :tokens => Dict(
            :comments => [
                (r"#.*$"m, COMMENT_SINGLE),
            ],
            :valid_name => [
                (r"[a-zA-Z][\w.]*", NAME),
                (r"\.[a-zA-Z_][\w.]*", NAME),
            ],
            :punctuation => [
                (r"(\[{1,2}|\]{1,2}|\(|\)|;|,)", PUNCTUATION),
            ],
            :keywords => [
                (words(keywords; suffix = "\\b(?![\\w.])"), KEYWORD_RESERVED),
                (words(types; suffix = "\\b(?![\\w.])"), KEYWORD_TYPE),
                (r"(library|require|attach|detach|source)\b(?![\\w.])", KEYWORD_NAMESPACE),
            ],
            :operators => [
                (r"(<<?-|->>?|-|==|<=|>=|<|>|&&?|!=|\|\|?|\?)", OPERATOR),
                (r"(\*|\+|\^|/|!|%[^%]*%|=|~|\$|@|:{1,3})", OPERATOR)
            ],
            :builtin_symbols => [
                (words(builtin_symbols; suffix = "(?![\\w.])"), KEYWORD_CONSTANT),
                (r"(T|F)\b", NAME_BUILTIN_PSEUDO),
            ],
            :numbers => [
                (r"0[xX][a-fA-F0-9]+([pP][0-9]+)?[Li]?", NUMBER_HEX),
                (r"[+-]?([0-9]+(\.[0-9]+)?|\.[0-9]+|\.)([eE][+-]?[0-9]+)?[Li]?", NUMBER),
            ],
            :statements => [
                :comments,
                # whitespaces
                (r"\s+", TEXT),
                (r"`.*?`", STRING_BACKTICK),
                (r"'", STRING, :string_squote),
                (r"\"", STRING, :string_dquote),
                :builtin_symbols,
                :numbers,
                :keywords,
                :valid_name,
                :operators,
                :punctuation,
            ],
            :root => [
                :statements,
                (r"\{|\}", PUNCTUATION),
                (r".", TEXT),
            ],
            :string_squote => [
                (r"([^'\\]|\\.)*'", STRING, :__pop__),
            ],
            :string_dquote => [
                (r"([^\"\\]|\\.)*\"", STRING, :__pop__),
            ],
        ),
    )
end


text:
include("FrequencyAnalysis.jl")
include("SubFrames.jl")
include("Reconstruction.jl")


text:
module TestPkg

export f1

"""
    f1(a)

Calculate f1.

# Examples

```jldoctest
julia> f1(3)
6
```

```jldoctest
julia> f1(4)
8
```
"""
f1(a) = 2a

end


text:
#
# Using LabelledArrays source code and types to
# "pass through" all properties to an underlying
# LArray or SLArray. This allows us to write
# a new abstract type that *functions*
# like a LabelledArray type!
#
# All code in this file was copied, and modified
# from LabelledArrays.jl source code. Their license is shown below.
# All credit goes to LabelledArrays.jl developers.
#

const __LABELLED_ARRAYS_LICENSE = """
The LabelledArrays.jl package is licensed under the MIT "Expat" License:

> Copyright (c) 2017: Christopher Rackauckas.
>
>
> Permission is hereby granted, free of charge, to any person obtaining a copy
>
> of this software and associated documentation files (the "Software"), to deal
>
> in the Software without restriction, including without limitation the rights
>
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
>
> copies of the Software, and to permit persons to whom the Software is
>
> furnished to do so, subject to the following conditions:
>
>
>
> The above copyright notice and this permission notice shall be included in all
>
> copies or substantial portions of the Software.
>
>
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
>
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
>
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
>
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
>
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
>
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
>
> SOFTWARE.
>
>
"""

const __LABELLED_ARRAYS_CREDITS = """
This source code which provides this functionality was copied directly from LabelledArrays.jl source code. The LabelledArrays.jl license text is provided in Julia's **Extended Help** section (accessible via `@doc`, or `??` in Julia's REPL).

# Extended Help

__LabelledArrays.jl License__

$__LABELLED_ARRAYS_LICENSE
"""

"""
$(TYPEDEF)

A supertype for types that *function* like
`LabelledArray.LArray` or `LabelledArray.SLArray`
instances, but are under a new type tree. This is
used in `GeneralAstrodynamics` for parameterizing
astrodynamics state vectors by physical units.

!!! note
    All subtypes __must__ have only one field:
    a `LabelledArrays.LArray` or `LabelledArrays.SLArray`
    field called `__rawdata`. All methods on this abstract
    type require this field to be called `__rawdata`!

Nearly all code which acts on this type is copied
and / or modified from LabelledArrays.jl source code.
All credit goes to LabelledArrays.jl developers.
The LabelledArrays.jl LICENSE file is provided
in this docstring under Julia's __Extended Help__
docstring section.

# Extended help

__LabelledArrays.jl License__

$__LABELLED_ARRAYS_LICENSE
"""
abstract type ParameterizedLabelledArray{F,N,T,L} <: DenseArray{F,N} end

"""
$(SIGNATURES)

Returns dot-accessible property names for a `ParameterizedLabelledArray`.
"""
Base.propertynames(::ParameterizedLabelledArray{F,N,T}) where {F,N,T} = T isa NamedTuple ? keys(T) : T

"""
$(SIGNATURES)

Overrides `Base.getproperty` for all `ParameterizedLabelledArray` instances.

$__LABELLED_ARRAYS_CREDITS
"""
Base.@propagate_inbounds function Base.getproperty(x::ParameterizedLabelledArray, s::Symbol)
    if s ∈ propertynames(x)
        return getproperty(getfield(x, :__rawdata), s)
    end
    return getfield(x, s) # will throw an error if s is not :__rawdata!
end

"""
$(SIGNATURES)

Sets indices of a `ParameterizedLabelledArray` via label.

$__LABELLED_ARRAYS_CREDITS
"""
Base.@propagate_inbounds function Base.setproperty!(x::ParameterizedLabelledArray, s::Symbol,y)
    if s ∈ propertynames(x)
        return setproperty!(getfield(x, :__rawdata), s, y)
    end
    setfield!(x, s, y)
end

"""
$(SIGNATURES)

Overrides `similar` for `ParameterizedLabelledArray` instances.

$__LABELLED_ARRAYS_CREDITS
"""
function Base.similar(x::ParameterizedLabelledArray, ::Type{S}, dims::NTuple{N,Int}) where {S,N}
    tmp = similar(x.__rawdata, S, dims)
    return typeof(x)(tmp)
end

"""
$(SIGNATURES)

Shallow copies a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.copy(x::ParameterizedLabelledArray) = typeof(x)(copy(getfield(x,:__rawdata)))

"""
$(SIGNATURES)

Deep copies a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.deepcopy(x::ParameterizedLabelledArray) = typeof(x)(deepcopy(getfield(x, :__rawdata)))

"""
$(SIGNATURES)

Copies one `ParameterizedLabelledArray` to another.

$__LABELLED_ARRAYS_CREDITS
"""
Base.copyto!(x::C,y::C) where C <: ParameterizedLabelledArray = copyto!(getfield(x,:__rawdata),getfield(y,:__rawdata))

"""
$(SIGNATURES)

Provides `unsafe_convert` for `ParameterizedLabelledArray` types for use with LAPACK.

$__LABELLED_ARRAYS_CREDITS
"""
Base.unsafe_convert(::Type{Ptr{T}}, a::ParameterizedLabelledArray{F}) where {T, F} = Base.unsafe_convert(Ptr{T}, getfield(a,:__rawdata))

"""
$(SIGNATURES)

Converts the underlying floating point type for a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.convert(::Type{T}, x) where {T<:ParameterizedLabelledArray} = T(x)

"""
$(SIGNATURES)

Converts the underlying floating point type for a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.convert(::Type{T}, x::T) where {T<:ParameterizedLabelledArray} = x

"""
$(SIGNATURES)

Converts the underlying floating point type for a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.convert(::Type{<:Array},x::ParameterizedLabelledArray) = convert(Array, getfield(x,:__rawdata))

"""
$(SIGNATURES)

Reshapes a `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
ArrayInterface.restructure(x::ParameterizedLabelledArray{F1}, y::ParameterizedLabelledArray{F2}) where {F1, F2} = reshape(y, size(x)...)

"""
$(SIGNATURES)

Implements `dataids` for a `ParameterizedLabelledArray` instance.

$__LABELLED_ARRAYS_CREDITS
"""
Base.dataids(A::ParameterizedLabelledArray) = Base.dataids(A.__rawdata)

"""
$(SIGNATURES)

Returns the index of the `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.getindex(state::ParameterizedLabelledArray, args...) = Base.getindex(state.__rawdata, args...)

"""
$(SIGNATURES)

Sets the index of the `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.setindex!(state::ParameterizedLabelledArray, args...) = Base.setindex!(state.__rawdata, args...)

"""
$(SIGNATURES)

Returns the memory stride for any `ParameterizedLabelledArray`.

$__LABELLED_ARRAYS_CREDITS
"""
Base.elsize(::ParameterizedLabelledArray{F}) where F = sizeof(F)


text:
# Translation of pylikwid example from https://github.com/RRZE-HPC/pylikwid
using Test
using Printf
using LIKWID

list = Int[]
cpus = [0,1]
@test LIKWID.PerfMon.init(cpus)
group = LIKWID.PerfMon.add_event_set("INSTR_RETIRED_ANY:FIXC0")
@test group == 0
@test LIKWID.PerfMon.setup_counters(group)
@test LIKWID.PerfMon.start_counters()
for i in 1:1_000_000
    push!(list,i)
end
@test LIKWID.PerfMon.stop_counters()
for thread in 1:length(cpus)
    tidx = thread-1
    @printf("Result CPU %d : %f\n", cpus[thread], LIKWID.PerfMon.get_result(group,0,tidx))
end


text:
using LinearAlgebra, Random, ExactD, Test

@testset "exponential operators" begin
Random.seed!(0)
L = 6
N = 3
basis = make_LN_basis(L, N)
state = rand(length(basis))

@testset "exponential J2" begin
    α = rand(L, L)
    cstate = deepcopy(state)
    exp_J2!(cstate, basis, α)
    @test cstate[1] ≈ state[1]*exp(α[1,2]+α[2,1]+α[1,3]+α[3,1]+α[2,3]+α[3,2]+α[1,1]+α[2,2]+α[3,3])
    @test cstate[4] ≈ state[4]*exp(α[4,2]+α[2,4]+α[4,3]+α[3,4]+α[2,3]+α[3,2]+α[4,4]+α[2,2]+α[3,3])
    @test cstate[8] ≈ state[8]*exp(α[1,4]+α[4,1]+α[1,5]+α[5,1]+α[4,5]+α[5,4]+α[1,1]+α[4,4]+α[5,5])
    @test cstate[16] ≈ state[16]*exp(α[4,6]+α[6,4]+α[4,3]+α[3,4]+α[6,3]+α[3,6]+α[4,4]+α[6,6]+α[3,3])
    @test cstate[20] ≈ state[20]*exp(α[6,5]+α[5,6]+α[6,4]+α[4,6]+α[5,4]+α[4,5]+α[6,6]+α[5,5]+α[4,4])
end

@testset "exponential J1" begin
    α = rand(L)
    cstate = deepcopy(state)
    exp_J1!(cstate, basis, α)
    @test cstate[1] ≈ state[1]*exp(α[1]+α[2]+α[3])
    @test cstate[4] ≈ state[4]*exp(α[4]+α[2]+α[3])
    @test cstate[8] ≈ state[8]*exp(α[1]+α[4]+α[5])
    @test cstate[16] ≈ state[16]*exp(α[4]+α[6]+α[3])
    @test cstate[20] ≈ state[20]*exp(α[6]+α[5]+α[4])
end

end


text:
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using MTHE493
using Plots

#%%

image_dir = "./data/small_dataset/"
dims = (112,92)
images = load_images(image_dir; dims=dims)

image_matrix = images_to_image_matrix(map(x->x[1], images))
labels = map(x->x[2], images)
label_map = Dict(String.(unique(labels)) .=> collect(1:length(unique(labels))))
proj, mean_vec = get_projection(image_matrix, 400)
tile_images(get_eigenfaces(proj, 25, dims), 5)

img = images[1][1]
[img eigenface_rep_to_image(image_to_eigenface_rep(img, proj, mean_vec), proj, mean_vec, dims)]

gr()
performance = [test_model(proj, images, d, 3, 1) for d in 1:15]
plot(performance)


text:
using Metatheory
using Metatheory.EGraphs
using Metatheory.Util
# Description here:
# https://www.philipzucker.com/metatheory-progress/

# https://github.com/AlgebraicJulia/Catlab.jl/blob/ce2fde9c63a8aab65cf2a7697f43cd24e5e00b3a/src/theories/Monoidal.jl#L127

cat_rules = @theory begin
    f ⋅ id(b) => f
    id(a) ⋅ f => f

    #f => f ⋅ id(dom(type(f)))
    #f => id(cod(type(f))) ⋅ f

    a ⊗ₒ munit() => a
    munit() ⊗ₒ a => a

    #a => a ⊗ₒ munit()
    #a => munit() ⊗ₒ a

    f ⋅ (g ⋅ h) == (f ⋅ g) ⋅ h

end

monoidal_rules = @theory begin
    id(munit()) ⊗ₘ f => f
    f ⊗ₘ id(munit()) => f
    a ⊗ₒ (b ⊗ₒ c) == (a ⊗ₒ b) ⊗ₒ c
    f ⊗ₘ (h ⊗ₘ j) == (f ⊗ₘ h) ⊗ₘ j
    id(a ⊗ₒ b) == id(a) ⊗ₘ id(b)

    (f ⊗ₘ p) ⋅ (g ⊗ₘ q) |>
    # future metatheory macros will clean this up
    begin
        fcod = addexpr!(_egraph, :(cod(type($f))))
        gdom = addexpr!(_egraph, :(dom(type($g))))
        pcod = addexpr!(_egraph, :(cod(type($p))))
        qdom = addexpr!(_egraph, :(dom(type($q))))
        if (fcod == gdom && pcod == qdom)
            :(($f ⋅ $g) ⊗ₘ ($p ⋅ $q))
        else
            :(($f ⊗ₘ $p) ⋅ ($g ⊗ₘ $q))
        end
    end

    (f ⋅ g) ⊗ₘ (p ⋅ q) => (f ⊗ₘ p) ⋅ (g ⊗ₘ q)
end

sym_rules = @theory begin
    σ(a, b) ⋅ σ(b, a) == id(a ⊗ₒ b)
    (σ(a, b) ⊗ₘ id(c)) ⋅ (id(b) ⊗ₘ σ(a, c)) == σ(a, b ⊗ₒ c)
    (id(a) ⊗ₘ σ(b, c)) ⋅ (σ(a, c) ⊗ₘ id(b)) == σ(a ⊗ₒ b, c)

    (f ⊗ₘ h) ⋅ σ(a, b) |> begin
        fcod = addexpr!(_egraph, :(cod(type($f)))).id
        hcod = addexpr!(_egraph, :(cod(type($h)))).id
        if fcod == a && hcod == b
            :(σ(dom(type($f)), dom(type($h))) ⋅ ($h ⊗ₘ $f))
        else
            :(($f ⊗ₘ $h) ⋅ σ($a, $b))
        end
    end

    σ(c, d) ⋅ (h ⊗ₘ f) |> begin
        fdom = addexpr!(_egraph, :(dom(type($f)))).id
        hdom = addexpr!(_egraph, :(dom(type($h)))).id
        if fdom == c && hdom == d
            :(($f ⊗ₘ $h) ⋅ σ(cod(type($f)), cod(type($h))))
        else
            :(σ($c, $d) ⋅ ($h ⊗ₘ $f))
        end
    end

    # these rules arer not catlab
    σ(a, munit()) == id(a)
    σ(munit(), a) == id(a)
    σ(munit(), munit()) => id(munit() ⊗ₒ munit())

end

diag_rules = @theory begin

    Δ(a) ⋅ (⋄(a) ⊗ₘ id(a)) == id(a)
    Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a)) == id(a)
    Δ(a) ⋅ σ(a, a) == Δ(a)

    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ σ(a, b) ⊗ₘ id(b)) == Δ(a ⊗ₒ b)

    Δ(a) ⋅ (Δ(a) ⊗ₘ id(a)) == Δ(a) ⋅ (id(a) ⊗ₘ Δ(a))
    ⋄(a ⊗ₒ b) == ⋄(a) ⊗ₘ ⋄(b)

    Δ(munit()) == id(munit())
    ⋄(munit()) == id(munit())
end

cart_rules = @theory begin

    Δ(a) ⋅ (f ⊗ₘ k) |> begin
        a_id = find(_egraph, a)
        if (
            addexpr!(_egraph, :(dom(type($f)))).id == a_id &&
            addexpr!(_egraph, :(dom(type($k)))).id == a_id
        )
            :(pair($f, $k))
        else
            :(Δ($a) ⋅ ($f ⊗ₘ $k))
        end
    end

    pair(f, k) => Δ(dom(type(f))) ⋅ (f ⊗ₘ k)
    proj1(a, b) == id(a) ⊗ₘ ⋄(b)
    proj2(a, b) == ⋄(a) ⊗ₘ id(b)
    f ⋅ ⋄(b) => ⋄(dom(type(f)))
    # Has to invent f. Hard to fix.
    # ⋄(b) => f ⋅ ⋄(b)

    f ⋅ Δ(b) => Δ(dom(type(f))) ⋅ (f ⊗ₘ f)
    Δ(a) ⋅ (f ⊗ₘ f) => f ⋅ Δ(cod(type(f)))
end

typing_rules = @theory begin
    dom(hom(a, b)) => a
    cod(hom(a, b)) => b
    type(id(a)) => hom(a, a)
    type(f ⋅ g) => hom(dom(type(f)), cod(type(g)))
    type(f ⊗ₘ g) => hom(dom(type(f)) ⊗ₒ dom(type(g)), cod(type(f)) ⊗ₒ cod(type(g)))
    type(a ⊗ₒ b) => :ob
    type(munit()) => :ob
    type(σ(a, b)) => hom(a ⊗ₒ b, b ⊗ₒ a)
    type(⋄(a)) => hom(a, munit())
    type(Δ(a)) => hom(a, a ⊗ₒ a)
    type(pair(f, g)) => hom(dom(type(f)), cod(type(f)) ⊗ₒ cod(type(g)))
    type(proj1(a, b)) => hom(a ⊗ₒ b, a)
    type(proj2(a, b)) => hom(a ⊗ₒ b, b)
end

rules = typing_rules ∪ cat_rules ∪ monoidal_rules ∪ sym_rules ∪ diag_rules ∪ cart_rules ∪ typing_rules

# A goofy little helper macro
# Taking inspiration from Lean/Dafny/Agda
macro calc(e...)
    theory = eval(e[1])
    e = rmlines(e[2])
    @assert e.head == :block
    for (a, b) in zip(e.args[1:end-1], e.args[2:end])
        println(a, " =? ", b)
        @time println(areequal(theory, a, b; timeout = 40))
    end
end

# Get the Julia motor hummin'
@calc rules begin

    ((⋄(a) ⊗ₘ ⋄(b)) ⋅ σ(munit(), munit()))
    (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a)))

end

@calc rules begin
    id(a ⊗ₒ b)
    id(a) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ (Δ(b) ⋅ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ ⋄(a)) ⊗ₘ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ (⋄(a) ⊗ₘ ⋄(b)) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ ((⋄(a) ⊗ₘ ⋄(b)) ⋅ σ(munit(), munit())) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⋅ id(a)) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ σ(a, b) ⊗ₘ id(b)) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b)))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (proj1(a, b) ⊗ₘ proj2(a, b))
    pair(proj1(a, b), proj2(a, b))
end

# shorter proof also accepted
@calc rules begin
    id(a ⊗ₒ b)
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    pair(proj1(a, b), proj2(a, b))
end


text:
using MATLAB

function patch(Γ)

    V = vertexarray(Γ)
    F = cellarray(Γ)

    @mput V F
    mat"""
        hold('on')
        patch('Vertices', V, 'Faces', F, 'FaceColor', 'red')
    """

    mat"axis equal"
    mat"colorbar"
    mat"view(3)"

end

function patch(Γ, C)

    V = vertexarray(Γ)
    F = cellarray(Γ)

    @mput V F C
    mat"""
        hold('on')
        patch('Vertices', V, 'Faces', F, 'FaceVertexCData', C, 'FaceColor', 'flat')
    """

    mat"axis equal"
    mat"colorbar"
    mat"view(3)"
end

function jmatlab_quiver(mesh, facecurs)
    els = elements(mesh)
    cts = center.(els)
    C = [ct[i] for ct in cts, i in 1:3]
    N = [fc[i] for fc in facecurs, i in 1:3]
    @mput C N
    mat"quiver3(C(:,1),C(:,2),C(:,3), N(:,1),N(:,2),N(:,3))"
end

function jmatlab_quiver(m)
    els = [chart(m, cells(m,i)) for i in 1:numcells(m)]
    C = zeros(numcells(m),3)
    N = zeros(numcells(m),3)
    for i in 1:size(C,1)
        s = els[i]
        p = neighborhood(s, [1,1]/3)
        n = normal(p)
        c = cartesian(p)
        C[i,:] = [c[1],c[2],c[3]]
        N[i,:] = [n[1],n[2],n[3]]
    end
    @mput C N
    mat"quiver3(C(:,1),C(:,2),C(:,3), N(:,1),N(:,2),N(:,3))"
end

function jmatlab_plot(x,y)
    @mput x
    @mput y
    mat"plot(x,y)"
end

# function quiver(Γ, fcr)
#
#   D = length(eltype(fcr))
#   V = vertexarray(Γ)
#   F = cellarray(Γ)
#
#   C = [cartesian(center(chart(Γ,c))) for c in cells(Γ)]
#   C = [C[i][j] for i in eachindex(C), j in 1:3]
#
#   Q = [fcr[i][j] for i in eachindex(fcr), j in 1:3]
#
#   @show size(C)
#   @show size(Q)
#
#   @mput V F C Q
#   mat"""
#   patch('Vertices', V, 'Faces', F, 'FaceColor', 'red')
#   quiver3(C(:,1),C(:,2),C(:,3),Q(:,1),Q(:,2),Q(:,3))
#   """
# end


text:
using EngEconomics

# Given
# Option 1: 1,500,000 install; 200,000 every 10 years; i=0.05; lifespan = infinite
installation = 1500000
A = 200000
i = 0.05
ic = (1 + i)^(1 * 10)-1
N = 10
P = installation + A / ic


text:
module MeasureToolbox

# Import the necessary packages.
import FastGaussQuadrature
import JuMP
import Distributions
const JuMPC = JuMP.Containers
using ..InfiniteOpt

# include jl files here
include("integrals.jl")
include("expectations.jl")
include("support_sums.jl")

# Export DataTypes
export Automatic, UniTrapezoid, UniMCSampling, UniIndepMCSampling, Quadrature,
GaussRadau, GaussHermite, GaussLegendre, GaussRadau, GaussLobatto, FEGaussLobatto,
GaussJacobi, GaussChebyshev, GaussLaguerre, MultiMCSampling, MultiIndepMCSampling,
InternalGaussLobatto

# Export methods
export uni_integral_defaults, set_uni_integral_defaults,
clear_uni_integral_defaults, integral, multi_integral_defaults,
set_multi_integral_defaults, clear_multi_integral_defaults, expect, support_sum,
generate_integral_data, 𝔼, ∫, generate_expect_data

# Export macros
export @expect, @𝔼, @support_sum, @integral, @∫

end # end module


text:
@testset "uniform.jl" begin
    # affine transformations
    test_affine_transformations(Uniform, rand(), 4 + rand())
end


text:
const chebyquad = let res_init=zeros(8), jac_init=zeros(8,8), x_init=[1:8;]/9

    function res(x, r)
        n = length(r)
        r[:] .= 0.0
        for j = 1:n
            tim1 = 1.0e0
            ti = 2.0e0 * x[j] - 1.0e0
            z = ti + ti
            for i = 1:n
                r[i] = r[i] + ti
                tip1 = z * ti - tim1
                tim1 = ti
                ti = tip1
            end
        end
        floatn = Float64(n)
        for i = 1:n
            ti = 0.0e0
            if mod(i,2) == 0
                ti = -1.0e0 / Float64(i*i - 1)
            end
            r[i] = ti - r[i]/floatn
        end
        return r
    end

    function jac(x, j)
        n = size(j)[1]
        for k = 1:n
            tim1 = -1.e0 / Float64(n)
            z = 2.e0 * x[k] - 1.e0
            ti = z * tim1
            tpim1 = 0.e0
            tpi = 2.e0 * tim1
            z = z + z
            for i = 1:n
                j[i, k] = tpi
                tpip1 = 4.e0 * ti + z * tpi - tpim1
                tpim1 = tpi
                tpi = tpip1
                tip1 = z * ti - tim1
                tim1 = ti
                ti = tip1
            end
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("chebyquad", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end