text:
#####
##### SHARP
#####

"""
    sharp(
        f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)},
        mask::AbstractArray{Bool, 3},
        vsz::NTuple{3, Real};
        r::Real = 18*minimum(vsz),
        thr::Real = 0.05,
    ) -> Tuple{typeof(similar(f)), typeof(similar(mask))}

Sophisticated harmonic artifact reduction for phase data (SHARP) [1].

### Arguments
- `f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)}`: unwrapped (multi-echo) field/phase
- `mask::AbstractArray{Bool, 3}`: binary mask of region of interest
- `vsz::NTuple{3, Real}`: voxel size for smv kernel

### Keywords
- `r::Real = 18*minimum(vsz)`: radius of smv kernel in units of `vsz`
- `thr::Real = 0.05`: threshold for high pass filter

### Returns
- `typeof(similar(f))`: background corrected local field/phase
- `typeof(similar(mask))`: eroded binary mask

### References
[1] Schweser F, Deistung A, Lehr BW, Reichenbach JR. Quantitative imaging of
    intrinsic magnetic tissue properties using MRI signal phase: an approach to
    in vivo brain iron metabolism?. Neuroimage. 2011 Feb 14;54(4):2789-807.
"""
function sharp(
    f::AbstractArray{<:AbstractFloat, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real};
    r::Real = 18*minimum(vsz),
    thr::Real = 0.05,
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return _sharp!(tzero(f), tzero(mask), f, mask, vsz, r, thr)
end

function _sharp!(
    fl::AbstractArray{<:AbstractFloat, N},
    smask::AbstractArray{Bool, 3},
    f::AbstractArray{T, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real},
    r::Real,
    thr::Real
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))

    size(fl) == size(f)        || throw(DimensionMismatch())
    size(smask) == size(mask)  || throw(DimensionMismatch())
    size(mask) == size(f)[1:3] || throw(DimensionMismatch())

    # crop image and pad for convolution
    Rc = crop_indices(mask)
    pad = ntuple(_ -> Int(cld(r, minimum(vsz))), Val(3))

    fp = padfastfft(@view(f[Rc,1]), pad, rfft=true)
    m  = padfastfft(@view(mask[Rc]), pad, rfft=true)

    # init vars and fft plans
    sz = size(m)
    sz_ = (sz[1]>>1 + 1, sz[2], sz[3])

    s = Array{T}(undef, sz)
    S = Array{T}(undef, sz_)
    F̂ = Array{complex(T)}(undef, sz_)

    FFTW.set_num_threads(num_cores())
    P = plan_rfft(fp)
    iP = inv(P)

    # get smv kernel
    S = _smv_kernel!(S, F̂, s, vsz, r, P)

    # constants
    _thr  = convert(eltype(S), thr)
    _one  = one(eltype(S))
    _zero = zero(eltype(F̂))
    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))

    # erode mask
    s = _tcopyto!(s, m) # in-place type conversion, reuse smv var

    F̂ = mul!(F̂, P, s)
    @inbounds @batch for I in eachindex(F̂)
        F̂[I] *= S[I]
        S[I] = _one - S[I]
    end

    s = mul!(s, iP, F̂)
    @inbounds @batch for I in eachindex(m)
        m[I] = s[I] > _δ
    end

    # SHARP
    @inbounds for t in axes(f, 4)
        if t > 1
            fp = padarray!(fp, @view(f[Rc,t]))
        end

        F̂ = mul!(F̂, P, fp)
        @batch for I in eachindex(F̂)
            F̂[I] *= S[I]
        end

        fp = mul!(fp, iP, F̂)
        @batch for I in eachindex(fp)
            fp[I] *= m[I]
        end

        F̂ = mul!(F̂, P, fp)
        @batch for I in eachindex(F̂)
            F̂[I] = ifelse(abs(S[I]) < _thr, _zero, F̂[I]*inv(S[I]))
        end

        fp = mul!(fp, iP, F̂)
        @batch for I in eachindex(fp)
            fp[I] *= m[I]
        end

        unpadarray!(@view(fl[Rc,t]), fp)
    end

    unpadarray!(@view(smask[Rc]), m)

    return fl, smask
end

#####
##### V-SHARP
#####

"""
    vsharp(
        f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)},
        mask::AbstractArray{Bool, 3},
        vsz::NTuple{3, Real};
        r::AbstractVector{<:Real} = 18*minimum(vsz):-2*maximum(vsz):2*maximum(vsz),
        thr::Real = 0.05,
    ) -> Tuple{typeof(similar(f)), typeof(similar(mask))}

Variable kernels sophisticated harmonic artifact reduction for phase data (V-SHARP) [1].

### Arguments
- `f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)}`: unwrapped (multi-echo) field/phase
- `mask::AbstractArray{Bool, 3}`: binary mask of region of interest
- `vsz::NTuple{3, Real}`: voxel size for smv kernel

### Keywords
- `r::AbstractVector{<:Real} = 18*minimum(vsz):-2*maximum(vsz):2*maximum(vsz)`:
    radii of smv kernels in mm
- `thr::Real = 0.05`: threshold for high pass filter

### Returns
- `typeof(similar(f))`: background corrected local field/phase
- `typeof(similar(mask))`: eroded binary mask

### References
[1] Wu B, Li W, Guidon A, Liu C. Whole brain susceptibility mapping using
    compressed sensing. Magnetic resonance in medicine. 2012 Jan;67(1):137-47.
"""
function vsharp(
    f::AbstractArray{<:AbstractFloat, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real};
    r::AbstractVector{<:Real} = 18*minimum(vsz):-2*maximum(vsz):2*maximum(vsz),
    thr::Real = 0.05,
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return length(r) == 1 ?
        _sharp!(tzero(f), tzero(mask), f, mask, vsz, first(r), thr) :
        _sharp!(tzero(f), tzero(mask), f, mask, vsz, r, thr)
end

function _sharp!(
    fl::AbstractArray{<:AbstractFloat, 3},
    smask::AbstractArray{Bool, 3},
    f::AbstractArray{T, 3},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real},
    r::AbstractVector{<:Real},
    thr::Real,
) where {T<:AbstractFloat}
    size(fl) == size(f)       || throw(DimensionMismatch())
    size(smask) == size(mask) || throw(DimensionMismatch())
    size(mask) == size(f)     || throw(DimensionMismatch())
    length(r) > 0             || throw(ArgumentError("r must not be empty"))

    rs = sort!(unique(collect(r)), rev=true)

    # crop image and pad for convolution
    Rc = crop_indices(mask)
    pad = ntuple(_ -> Int(cld(maximum(rs), minimum(vsz))), Val(3))

    fp = padfastfft(@view(f[Rc]), pad, rfft=true)
    mr = padfastfft(@view(mask[Rc]), pad, rfft=true)

    # init vars and fft plans
    sz = size(mr)
    sz_ = (sz[1]>>1 + 1, sz[2], sz[3])

    s = Array{T}(undef, sz)
    m = Array{Bool}(undef, sz)

    m = tfill!(m, 0)

    FFTW.set_num_threads(num_cores())
    P = plan_rfft(s)
    iP = inv(P)

    S  = Array{T}(undef, sz_)
    iS = Array{T}(undef, sz_)
    F̂  = Array{complex(T)}(undef, sz_)
    M̂  = Array{complex(T)}(undef, sz_)
    F̂p = Array{complex(T)}(undef, sz_)

    F̂p = mul!(F̂p, P, fp)
    fp = tfill!(fp, 0)

    # constants
    _thr  = convert(eltype(S), thr)
    _one  = one(eltype(S))
    _zero = zero(eltype(F̂))
    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))

    # fft of original mask
    s = _tcopyto!(s, mr) # in-place type conversion
    M̂ = mul!(M̂, P, s)

    @inbounds for (i, r) in enumerate(rs)
        # get smv kernel
        S = _smv_kernel!(S, F̂, s, vsz, r, P)

        # erode mask
        @batch for I in eachindex(F̂)
            F̂[I] = S[I] * M̂[I]
            S[I] = _one - S[I]
        end

        s = mul!(s, iP, F̂)
        @batch for I in eachindex(mr)
            mr[I] = s[I] > _δ
        end

        # high-pass filter first (largest) kernel for deconvolution
        if i == 1
            @batch for I in eachindex(iS)
                iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I]))
            end
        end

        # SHARP
        @batch for I in eachindex(F̂)
            F̂[I] = S[I] * F̂p[I]
        end

        s = mul!(s, iP, F̂)
        @batch for I in eachindex(fp)
            if mr[I] && !m[I]
                fp[I] = s[I]
            end
        end

        m, mr = mr, m
    end

    # deconvolution + high-pass filter
    F̂p = mul!(F̂p, P, fp)
    @batch for I in eachindex(F̂p)
        F̂p[I] *= iS[I]
    end

    fp = mul!(fp, iP, F̂p)
    @batch for I in eachindex(fp)
        fp[I] *= m[I]
    end

    unpadarray!(@view(fl[Rc]), fp)
    unpadarray!(@view(smask[Rc]), m)

    return fl, smask
end

function _sharp!(
    fl::AbstractArray{<:AbstractFloat, 4},
    smask::AbstractArray{Bool, 3},
    f::AbstractArray{T, 4},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real},
    r::AbstractVector{<:Real},
    thr::Real,
) where {T<:AbstractFloat}
    size(fl) == size(f)        || throw(DimensionMismatch())
    size(smask) == size(mask)  || throw(DimensionMismatch())
    size(mask) == size(f)[1:3] || throw(DimensionMismatch())
    length(r) > 0              || throw(ArgumentError("r must not be empty"))

    rs = sort!(unique(collect(r)), rev=true)

    # crop image and pad for convolution
    Rc = crop_indices(mask)
    pad = ntuple(_ -> Int(cld(maximum(rs), minimum(vsz))), Val(3))

    fp = padfastfft(@view(f[Rc,:]), pad, rfft=true)
    mr = padfastfft(@view(mask[Rc]), pad, rfft=true)

    # init vars and fft plans
    sz = size(mr)
    sz_ = (sz[1]>>1 + 1, sz[2], sz[3])

    s = Array{T}(undef, sz)
    m = Array{Bool}(undef, sz)

    m = tfill!(m, 0)

    FFTW.set_num_threads(num_cores())
    P = plan_rfft(s)
    iP = inv(P)

    P4 = plan_rfft(fp, 1:3)
    iP4 = inv(P4)

    S  = Array{T}(undef, sz_)
    iS = Array{T}(undef, sz_)
    F̂  = Array{complex(T)}(undef, sz_)
    M̂  = Array{complex(T)}(undef, sz_)
    F̂p = Array{complex(T)}(undef, (sz_..., size(fp, 4)))

    F̂p = mul!(F̂p, P4, fp)
    fp = tfill!(fp, 0)

    # constants
    _thr  = convert(eltype(S), thr)
    _one  = one(eltype(S))
    _zero = zero(eltype(F̂))
    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))

    # fft of original mask
    s = _tcopyto!(s, mr) # in-place type conversion
    M̂ = mul!(M̂, P, s)

    @inbounds for (i, r) in enumerate(rs)
        # get smv kernel
        S = _smv_kernel!(S, F̂, s, vsz, r, P)

        # erode mask
        @batch for I in eachindex(F̂)
            F̂[I] = S[I] * M̂[I]
            S[I] = _one - S[I]
        end

        s = mul!(s, iP, F̂)
        @batch for I in eachindex(mr)
            mr[I] = s[I] > _δ
        end

        # high-pass filter first (largest) kernel for deconvolution
        if i == 1
            @batch for I in eachindex(iS)
                iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I]))
            end
        end

        # SHARP
        for t in axes(fp, 4)
            F̂t = @view(F̂p[:,:,:,t])
            ft = @view(fp[:,:,:,t])

            @batch for I in eachindex(F̂)
                F̂[I] = S[I] * F̂t[I]
            end

            s = mul!(s, iP, F̂)
            @batch for I in eachindex(ft)
                if mr[I] && !m[I]
                    ft[I] = s[I]
                end
            end
        end

        m, mr = mr, m
    end

    # deconvolution + high-pass filter
    F̂p = mul!(F̂p, P4, fp)
    @inbounds for t in axes(F̂p, 4)
        F̂t = @view(F̂p[:,:,:,t])
        @batch for I in eachindex(F̂t)
            F̂t[I] *= iS[I]
        end
    end

    fp = mul!(fp, iP4, F̂p)
    @inbounds for t in axes(fp, 4)
        ft = @view(fp[:,:,:,t])
        @batch for I in eachindex(ft)
            ft[I] *= m[I]
        end
    end

    unpadarray!(@view(fl[Rc,:]), fp)
    unpadarray!(@view(smask[Rc]), m)

    return fl, smask
end


text:
using LinearAlgebra, Test
using SumOfSquares
using DynamicPolynomials

"""
    simple_matrix(optimizer, config::MOIT.Config)

Example 3.77 and 3.79 of
Blekherman, G., Parrilo, P. A., & Thomas, R. R. (Eds.).
Semidefinite optimization and convex algebraic geometry SIAM 2013
"""
function simple_matrix_test(optimizer, config::MOIT.Config)
    @polyvar x
    P = [x^2-2x+2 x; x x^2]

    # Example 3.77
    mat_model = _model(optimizer)
    PolyJuMP.setpolymodule!(mat_model, SumOfSquares)
    mat_cref = @constraint(mat_model, Symmetric(P) in PSDCone())
    JuMP.optimize!(mat_model)
    @test JuMP.termination_status(mat_model) == MOI.OPTIMAL
    @test JuMP.primal_status(mat_model) == MOI.FEASIBLE_POINT
    @test length(gram_matrix(mat_cref).basis.monomials) == 3

    # Example 3.79
    @polyvar y[1:2]
    model = _model(optimizer)
    PolyJuMP.setpolymodule!(model, SumOfSquares)
    cref = @constraint(model, dot(vec(y), P * vec(y)) >= 0)
    JuMP.optimize!(model)
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test length(gram_matrix(cref).basis.monomials) == 3
end
sd_tests["simple_matrix"] = simple_matrix_test

ext: jl

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Langevin

include("transform.jl")


text:
#!/usr/bin/env julia

using DrWatson
quickactivate(@__DIR__)
@show projectdir()
#-----------------------------------------------------------------------------
# load packages
using CovidSurvey
using ArgParse
using Plots
using Dates

# using StatsPlots
using PrettyTables
using DataFrames
using BSON

using Turing
using Memoization
using ReverseDiff
using InlineStrings
using WeakRefStrings
Turing.setrdcache(true)
setadbackend(:reversediff)

plotlyjs()
#-----------------------------------------------------------------------------
# read input
argtable = ArgParseSettings(
    description="This script applies post-processing to Turing.Chains."
)
@add_arg_table! argtable begin
    "--fname", "-f"
        help = "file path to Turing.Chains"
    "--exclude", "-e"
        help = "exclude chains by id. Expected format: 1,2,3 "
        arg_type = String
        default = ""
    "--plot-results", "-p"
        help = "plot and display with firefox"
        arg_type = Bool
        default = false
    "--warmup", "-w"
        help = "set nr. of warmup samples"
        arg_type = Int64
        default = nothing
end

parsed_args = parse_args(ARGS, argtable)
exclude = isempty(parsed_args["exclude"]) ? [] : [parse(Int64,i) for i in split(parsed_args["exclude"], ",")]
#------------------------------------------------------------------------------
Regional.postprocessing(
    parsed_args["fname"];
    plot_results = parsed_args["plot-results"],
    exclude,
    warmup = parsed_args["warmup"]
)

fname = "/home/ankoh/dev/CovidSurvey/reports/regional_effect_estimation/mobility/apple_google_totalcontacts/2021-10-06/CHAIN_chains=10_model=hospit_observ=2021-02-01_preds=CS,MG,MA_regional=true_rwstep=7_semipar=true_std=false_steps=1000_warmup=2000.jls"
isfile(fname)
Regional.postprocessing(
    fname;
    plot_results = false,
    exclude = [],
)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module ShaleDrillingLikelihood_MakeModelsTest

using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters
using Test

@testset "Make Test Datset" begin
    num_i = 50
    (datasm, thetasm, ), (datalg, thetalg,) = MakeTestData(;num_i=num_i)

    @test length.(datasm) == (num_i, 0, 0,)
    @show length.(datalg) == (num_i, num_i, num_i)

end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((1, 1, 1, 2, true, true, false), nothing, relaxed_tols),
    ((1, 1, 1, 2, true, false, true),),
    ((1, 1, 1, 2, false, true, false),),
    ((1, 1, 1, 2, false, false, true),),
    ((1, 1, 1, 2, false, false, false),),
    ]
insts["fast"] = [
    ((2, 2, 2, 2, true, true, false),),
    ((2, 2, 2, 2, true, false, true),),
    ((2, 2, 2, 2, false, true, false),),
    ((2, 2, 2, 2, false, false, true),),
    ((2, 2, 2, 2, false, false, false),),
    ((4, 2, 3, 5, true, true, false),),
    ((4, 2, 3, 5, true, false, true),),
    ]
insts["various"] = vcat(insts["fast"], [
    ((4, 2, 3, 5, false, true, false),),
    ((4, 2, 3, 5, false, false, true),),
    ((4, 2, 3, 5, false, false, false),),
    ])
return (PolyNormJuMP, insts)


text:
function trainfunction1d(func,x)
    if func == 1
        u = sin(4*pi*x)
    elseif func == 2
        u = sin(2*pi*x)
    elseif func == 3
        u = sin(pi*x)
    elseif func == 4
        u = sin(3*pi*x)
    elseif func == 5
        u = sin(2.5*pi*x)
    elseif func == 6
        u = sin(-3*pi*x)
    elseif func == 7
        u = sin(-4*pi*x)
    elseif func == 8
        u = sin(-1*pi*x)
    elseif func == 9
        u = sin(1.5*pi*x)
    elseif func == 10
        u = sin(0.5*pi*x)
    elseif func == 11
        u =  2*x
    elseif func == 12
        u =  6*x
    elseif func == 13
        u =  -4*x
    elseif func == 14
        u =  20*x
    elseif func == 15
        u =  -10*x
    elseif func == 16
        u =  0.2*x
    elseif func == 17
        u =  -0.6*x
    elseif func == 18
        u =  -2*x
    elseif func == 19
        u = 5 * abs(x)
    elseif func == 20
        u = -2 * abs(x)
    elseif func == 21
        u = 10 * abs(x)
    elseif func == 22
        u = 3 * abs(x)
    elseif func == 23
        u = -5 * abs(x)
    elseif func == 24
        u = 0.8 * abs(x)
    elseif func == 25
        u = 0.1 * abs(x)
    elseif func == 26
        u = -4 * abs(x)
    end
    return u
end

function Output1d(func, xi, h)
    if func <= 18
      Y = [0; 1]
    elseif func > 18 && func <=26
      if 0 >= xi-(1/2)*h && 0 <= xi+(1/2)*h
        Y = [1; 0]
      else
        Y = [0; 1]
      end
    elseif func == 27
      Y = [1; 0]
    end
    return Y
end

function troubledcellfunctionstep1d(x, ul, ur, x0)
    u = ul*(x<=x0) + ur*(x>x0)
    return u
end

function good_cell1d(node_coord, h, x0)
    x = node_coord[1] + h/2
    if x0 < x+(1/2)*h && x0 > x-(1/2)*h
        return false
    else
        return true
    end

end

function validfunction1d(func,x)
    if func == 1
        u = sin(2*pi*x)
    elseif func == 2
        u = sin(1*pi*x)
    elseif func == 3
        u = sin(3*pi*x)
    elseif func == 4
        u = sin(1*pi*x) +sin(2*pi*x)+ sin(3*pi*x)
    elseif func == 5
        u = sin(1*pi*x)+sin(2*pi*x)+sin(2*pi*x)+sin(4*pi*x)+sin(5*pi*x)
    elseif func == 6
        u = sin(4*pi*x)
    elseif func == 7
        u = sin(3*pi*x)*sin(2*pi*x)*cos(5*pi*x)
    elseif func == 8
        u = sin(1*pi*x)*sin(4*pi*x)*cos(2*pi*x)
    elseif func == 9
        u = sin(4*pi*x)*sin(0.2*pi*x)*cos(6*pi*x)
    elseif func == 10
        u = sin(pi*x)*sin(3*pi*x)*cos(0.5*pi*x)
    elseif func == 11
        u =  sin(0.2pi*x)*sin(2*pi*x)*cos(3.5*pi*x)
    elseif func == 12
        u =  sin(pi*x) + exp(x)
    elseif func == 13
        u =  sin(2*pi*x) + exp(x)
    elseif func == 14
        u =  sin(0.2*pi*x) + exp(x)
    elseif func == 15
        u =  sin(3.2*pi*x) + exp(x)
    elseif func == 16
        u = sin(5*pi*x) + exp(x)
    end
    return u
end


text:
const UNCONVERGED_STATUS = :unconverged
const CONVERGED_STATUS = :converged
const UNSET_STATUS = :unset

is_converged(out::EPOut) = (out.status == CONVERGED_STATUS)


text:
include("../hopfield.jl")
using Gadfly, Colors

function circLayer(x::Vector{Float64};r::Float64=1.0,c="red")
  line = 0:0.05:2pi
  N = length(line)
  X, Y = zeros(N), zeros(N)
  for i in 1:N
    X[i] = x[1] + r*cos(line[i])
    Y[i] = x[2] + r*sin(line[i])
  end
  return Gadfly.layer(x=X,y=Y, Geom.path(), Gadfly.Theme(default_color=parse(Colorant,c)))
end

function drawLandmPtsLayer(pts::Dict{String, Vector{Float64}}; c::String="blue")
  N = length(pts)
  xy = zeros(2,N)
  lbls = String[]
  i = 0
  lbl = String[]
  for p in pts
    i += 1
    push!(lbls, p[1])
    [xy[j,i] = p[2][j] for j in 1:2]
    push!(lbl, p[1])
  end
  return layer(x=xy[1,:],y=xy[2,:], label=lbl, Geom.point, Geom.label, Gadfly.Theme(default_color=parse(Colorant,c)))
end

incirc(dx,dy,r) = dx^2 + dy^2 <= r^2

function findInRadius(pts::Dict{String, Vector{Float64}}, from::String; r::Float64=1.0)
  rd = Dict{String, Vector{Float64}}()
  xy = pts[from]
  rd[from] = pts[from]
  for p in pts
    if incirc(xy[1]-p[2][1], xy[2]-p[2][2], r)
      rd[p[1]] = p[2]
    end
  end
  return rd
end

function drawLineLayer(from::Vector{Float64}, to::Vector{Float64}; c="magenta")
  return Gadfly.layer(x=[from[1];to[1]],y=[from[2];to[2]], Geom.path(), Gadfly.Theme(default_color=parse(Colorant,c)))
end

function drawDictLinesLayers(pl, pts::Dict{String, Vector{Float64}}, from::String)
  x0 = pts[from][1:2]
  for p in pts
    p[1] == from ? continue : nothing
    l = drawLineLayer(x0, p[2][1:2])
    push!(pl.layers, l[1])
  end
  return pl
end

calcang(x1::Vector{Float64},x2::Vector{Float64}) = atan2(x2[2]-x1[2], x2[1]-x1[1])
calcrange(x1::Vector{Float64},x2::Vector{Float64}) = norm(x1-x2)

function orderedBinIndex(x::Float64; start::Float64=0.0, stop::Float64=1.0,N::Int=3)
  refSpace= linspace(start, stop, N+1)
  if N == 1
    return 1
  end
  for i in 2:(N+1)
    if refSpace[i-1] <= x < refSpace[i]
      return i-1
    end
  end
  return N
end

function constructShapeContexctMat(range::Dict{String, Float64}, ang::Dict{String, Float64}, me::String;
                                   angBins::Int=4, rangeBins::Int=3, r::Float64=1.0)
  sc = zeros(rangeBins, angBins)
  for a in ang
    if String(a[1]) == me  continue; end
    # @show a[1], a[2], range[a[1]]
    aIdx = orderedBinIndex(a[2], start=Float64(-pi), stop=Float64(pi), N=angBins)
    rIdx = orderedBinIndex(range[a[1]], start=0.0, stop=r, N=rangeBins)
    # @show rIdx, aIdx
    sc[rIdx, aIdx] = sc[rIdx, aIdx] + 1
  end
  return sc
end

function myShapeContext(pts::Dict{String, Vector{Float64}}, me::String, refpt::String;
                        r::Float64=1.0, angBins::Int=4, rangeBins::Int=3)
  rd = findInRadius(pts, me, r=r)
  ang = Dict{String, Float64}()
  rang = Dict{String, Float64}()
  # scVec = Dict{String}
  refang = calcang(rd[me], pts[refpt])
  for p in rd
    ang[p[1]] = calcang(rd[me], p[2]) - refang
    rang[p[1]] = calcrange(rd[me], p[2])
  end
  # @show ang
  # @show rang
  mat = constructShapeContexctMat(rang, ang, me, r=r, angBins=angBins, rangeBins=rangeBins)
  return mat
end

function calcShapeContext(pts::Dict{String, Vector{Float64}}, refpt::String;
                          r::Float64=1.0, angBins::Int=4, rangeBins::Int=3)
  rd = findInRadius(pts,refpt,r=r)
  mat = zeros(rangeBins, angBins)
  for doone in keys(rd)
    if String(doone) == refpt  continue; end
    # @show doone
    mat = mat + myShapeContext(pts, doone, refpt, r=r, angBins=angBins, rangeBins=rangeBins)
  end
  return mat
end

function shapeContMatToVec(mat::Array{Float64,2}; duplicity::Int=2)
  v = [map(Int,mat[:].==1)'; map(Int,mat[:].>=2)']#; map(Int,mat[:].>2)']
  # v = map(Int,mat[:].>=1)'
  r,c = size(v)
  V = reshape(v,floor(Int,r*c),1)
  return sign(V[:]-0.5)
end

function shapeContPattern(pts::Dict{String, Vector{Float64}}, refpt::String;
                          r::Float64=1.0, angBins::Int=4, rangeBins::Int=3)
    mat = calcShapeContext(pts, refpt, r=r, angBins=angBins, rangeBins=rangeBins)
    return shapeContMatToVec(mat)
end

# we want a way to select the top 3 descriptors from current list of landmarks
# parameter
function bestDescriptors(d::Dict{String, Vector{Float64}};
                          r::Float64=1.0,Mf::Int=4,angBins::Int=4,rangeBins::Int=3,k::Int=3)

  allDescr = Dict{String, Vector{Int}}()
  skey = collect(keys(d))
  for lm in skey
    if length(findInRadius(d,lm,r=r)) > Mf
      allDescr[lm] = shapeContPattern(d, lm, r=r, angBins=angBins, rangeBins=rangeBins)
    end
  end
  @show ks = collect(keys(allDescr))

  #construct matrix of hamming distances
  n = length(ks)
  H = 99999*eye(Int, n, n)
  for i in 1:(n-1), j = (i+1):n
    H[i,j] = hamming(allDescr[ks[i]],allDescr[ks[j]])
  end
  H = H + H'

  # sort for max of min distances
  Hmin = vec(minimum(H,2))
  p = sortperm(Hmin,rev=true)
  # @show Hmin[p]
  # @show ks[p]

  # take top half and redo
  H1a = H[:,p]
  H1b = H1a[p,:]

  k = n >= k ? k : n
  # @show nn = k #floor(Int, n/2+0.5)
  H1 = H1b[1:k, 1:k]
  @show kss = (ks[p])[1:k]

  H1min = vec(minimum(H1,2))
  p1 = sortperm(H1min,rev=true)
  # @show H1min[p1]
  # @show kss[p1]
  bestlbls = kss[p1]#[1:k]]
  bestDescr = Dict{String, Vector{Int}}()
  for lb in bestlbls
    bestDescr[lb] = allDescr[lb]
  end
  return bestlbls, bestDescr
end

function findMatches(dd::Dict{String, Vector{Float64}}, Patt::Array{Float64,2}, bLbls::Vector{String};
  R::Float64=1.0,angBins::Int=4,rangeBins::Int=3, Mf::Int=4)
  SV = Dict{String, Vector{Int}}()
  Psv = Dict{String, Int}()
  for i in 1:length(dd)
    l = length(findInRadius(dd,"l$(i)",r=R))
    if l > Mf
      val = shapeContPattern(dd, "l$(i)", r=R, angBins=angBins, rangeBins=rangeBins)
      SV[String("l$(i)")] = val
      Psv[String("l$(i)")] = searchHopfield(Patt, val, maxiter=4)
    end
  end

  # display matching results
  mtchs = Dict{Int, Any}()
  fm = Dict{Int, Any}()
  # @show bLbls
  i = 0
  ord = Int[]
  for psv in Psv
    if psv[2] != -1
      i+=1
      hd = hamming(SV[psv[1]], Patt[:,psv[2]])
      mtchs[i] = (psv[1], bLbls[psv[2]], hd)
      push!(ord, hd)
    end
  end
  p = sortperm(ord)
  i = 0
  for j in p
    i+=1
    fm[i] = mtchs[j]
  end
  return fm
end

function convertToMMLst(fm::Dict{Int,Any}, distLim::Int)
  mmlst = Dict{String, Array{Tuple{String, Int},1}}()
  for i in 1:length(fm)
    if fm[i][3] <= distLim
      if haskey(mmlst,fm[i][2])
          push!(mmlst[fm[i][2]], (fm[i][1],fm[i][3]) )
      else
          mmlst[fm[i][2]] =  [ (fm[i][1],fm[i][3])  ]
      end
    end
  end
  return mmlst
end

# start some practicatical tests
if false
d = Dict{String, Vector{Float64}}()
for i in 1:9
  d[String("l$(i)")] = rand(2)
end
for i in 10:20
  d[String("l$(i)")] = (3.0.*rand(2))-1
end

# draw points
pl=plot(drawLandmPtsLayer(d))

R = 0.4
Mf = 4
angBins=6
bLbls, bestDescr = bestDescriptors(d, r=R, Mf=Mf, k=4 )

# draw so we can see
for lm in bLbls
  push!(pl,circLayer(d[lm],r=R)[1])
end
pl

# Training pattern vectors
# @show Plbl = [bLbls[1];bLbls[2];bLbls[3]]#; bLbls[4]]
@show bLbls
v1 = bestDescr[bLbls[1]]
v2 = bestDescr[bLbls[2]]
v3 = bestDescr[bLbls[3]]
v4 = bestDescr[bLbls[4]]

P = [v1';v2';v3';v4']'
P = map(Float64,P)

hamming(v1,v2)
hamming(v1,v3)
hamming(v2,v3)

hamming(v1,v4)
hamming(v2,v4)
hamming(v3,v4)

# how well do these patterns work

dd = deepcopy(d)
@show length(d)
for i in 1:length(d)
  dd[String("l$(i)")] = dd[String("l$(i)")] + 0.0*randn(2)
end
for i in length(d):floor(Int,length(d)*1.0)
  dd[String("l$(i)")] = (3.0.*rand(2))-1
end

findMatches(dd, P, bLbls, R=R, angBins=angBins, Mf=Mf)

l2 = drawLandmPtsLayer(dd,c="green")
push!(pl.layers, l2[1])
pl

# plot(drawLandmPtsLayer(dd))

# random drawing
refpt = "l8"
rd = findInRadius(d,refpt,r=R)
pl=plot(drawLandmPtsLayer(d),
     circLayer(d[refpt],r=R))
drawDictLinesLayers(pl, rd, refpt)

# sanity check

val1 = shapeContPattern(d, "l3", r=R, angBins=angBins)
val2 = shapeContPattern(dd, "l3", r=R, angBins=angBins)
@show hamming(val1,val2)

end

# some training and recognition tests
# ks = collect(keys(rd))
# doone = ks[1] != "l1" ? ks[1] : ks[2]
# doone="l6"
# push!(pl,circLayer(rd[doone],r=0.4)[1])
#
# mat = zeros(3, 4)
# for doone in ["l6";"l8";"l10"]
#   @show doone
#   @show mat = mat + myShapeContext(d, doone, "l1", r=0.4)
# end
#
# v1 = shapeContPattern(d, "l2", r=0.4)
# v2 = shapeContPattern(d, "l4", r=0.4)
# v3 = shapeContPattern(d, "l7", r=0.4)
# v4 = shapeContPattern(d, "l10", r=0.4)
#
# P = [v1';v2';v3']'
#
# a1 = deepcopy(v10)
# a1[1] = -1*v10[1]
# a1[10] = -1*v10[10]
# @show searchHopfield(P, a1)

# trying sorting

# @show A = rand(4,4)
#
# A = 16*eye(4,4)
# for i in 1:3, j in (i+1):4
# A[i,j] = i*j
# end
# A = A + A'
# @show A
#
# [A[i,i]=24.0 for i in 1:4]
# @show Amin = vec(minimum(A,2))
# @show p = sortperm(Amin,rev=true)
# A1 = A[:,p]
# @show A2 = A1[p,:]


text:
using CommutativeFunctions
using Test

@testset "CommutativeFunctions.jl" begin
    @commutative f(a::Int, b::Float64) = a + b
    @test f(2, 3.0) == 5
    @test f(3.0, 2) == 5

    @test_throws Exception @eval (@commutative g(x) = x)
    @test_throws Exception @eval (@commutative g(x, y, z) = x)
    @test_throws Exception @eval @commutative :+
end


text:
# Decision Tree Module

module Tree
export DecisionTree, create_decision_tree, predict

"""
Data structure for Decision Tree\\
One of `children` and `target` will be `nothing`
"""
mutable struct DecisionTree
    col_id::Integer
    children::Union{AbstractDict{Function,DecisionTree},Nothing}
    target::Union{AbstractFloat,Nothing}
end

"""
C4.5 implementation for Decision Tree\\
`X_data` has 2 dims, `Y_data` has 1 dims\\
Returns a computed `DecisionTree` structure
"""
function create_decision_tree(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number; max_depth::Integer=10)::DecisionTree
    @assert ndims(X_data) == ndims(Y_data) + 1 == 2
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert max_depth >= 1

    """
    Output the most frequent value in the vector
    """
    function majority(Y_vec::Array)::Number
        unique_vals = Dict{Number, Integer}()
        for Y_val in Y_vec
            if !haskey(unique_vals, Y_val)
                unique_vals[Y_val] = 1
            else
                unique_vals[Y_val] += 1
            end
        end
        result = sort(collect(unique_vals), by=m->m[2])
        return result[end][1]
    end

    """
    Check if vector has same value
    """
    function is_same_class(Y_vec::Array)::Bool
        val = Y_vec[1]
        for m in Y_vec[2:end]
            if m != val
                return false
            end
        end
        return true
    end

    """
    Get an array of unique values from a vector
    """
    function get_unique_array(X_vec::Array)::Array
        X_max = maximum(X_vec)
        X_min = minimum(X_vec)
        round_digits = -Integer(trunc(log10(X_max-X_min)-2))
        result = []
        for m in X_vec
            if !(trunc(m, digits=round_digits) in result)
                push!(result, m)
            end
        end
        return result
    end

    """
    The entropy function
    """
    function info(Y_data::Array)::AbstractFloat
        portion = Dict{Number,Integer}()
        for m in Y_data
            if !haskey(portion, m)
                portion[m] = 1
            else
                portion[m] += 1
            end
        end
        T = size(Y_data)[1]
        result = 0.0
        for m in keys(portion)
            p = portion[m] / T
            result += -p * log2(p)
        end
        return result
    end

    """
    Detect types for each column of `X_data`\\
    `true` if continuous\\
    `false` if categorical
    """
    function detect_types(X_data::Array)::Array
        result = Array{Bool}(undef, size(X_data)[2])
        for i in 1:size(X_data)[2]
            X_unique = get_unique_array(X_data[:, i])
            if length(X_unique) > 10
                result[i] = true
            else
                result[i] = false
            end
        end
        return result
    end

    function create_decision_tree_recursive(X_data::Array, Y_data::Array, X_type::Array, max_depth::Integer)::DecisionTree
        if is_same_class(Y_data)
            return DecisionTree(-1, nothing, Y_data[1])
        end
        if max_depth <= 0
            return DecisionTree(-1, nothing, majority(Y_data))
        end
        max_index = -1
        max_ratio = 0.0
        max_fn = nothing
        for i in 1:size(X_data)[2]
            X_vec_unique = get_unique_array(X_data[:, i])
            if X_type[i]
                # continous data, then choose binary threshold
                fn_max = nothing
                ratio_max = 0
                for j in 1:(length(X_vec_unique)-1)
                    # choose average as threshold
                    threshold = (X_vec_unique[j] + X_vec_unique[j+1]) / 2.0
                    fn = (m) -> m <= threshold
                    # compute information gain
                    Ii = info(Y_data)
                    Y_1 = Y_data[fn.(X_data[:, i])]
                    Y_2 = Y_data[(!fn).(X_data[:, i])]
                    I1 = info(Y_1)
                    I2 = info(Y_2)
                    gain = Ii - (length(Y_1)/length(Y_data)*I1 + length(Y_2)/length(Y_data)*I2)
                    splitInfo = -(length(Y_1)/length(Y_data))*log2(length(Y_1)/length(Y_data)) - (length(Y_2)/length(Y_data))*log2(length(Y_2)/length(Y_data))
                    ratio = gain / splitInfo
                    if fn_max === nothing
                        fn_max = fn
                        ratio_max = ratio
                    else
                        if ratio > ratio_max
                            fn_max = fn
                            ratio_max = ratio
                        end
                    end
                end
                if ratio_max > max_ratio
                    max_index = i
                    max_fn = fn_max
                    max_ratio = ratio_max
                end
            else
                # skip visited categorical feature (homogeneous)
                if length(X_vec_unique) <= 1
                    continue
                end
                # categorical data, then compute by splitting on each category
                Ii = info(Y_data)
                gain = Ii
                splitInfo = 0.0
                for val in X_vec_unique
                    fn = (m) -> m == val
                    Y0 = Y_data[fn.(X_data[:, i])]
                    I0 = info(Y0)
                    gain -= (length(Y0)/length(Y_data))*I0
                    splitInfo -= (length(Y0)/length(Y_data))*log2(length(Y0)/length(Y_data))
                end
                ratio = gain / splitInfo
                if ratio > max_ratio
                    max_index = i
                    max_fn = nothing
                    max_ratio = ratio
                end
            end
        end
        if max_index < 0
            return DecisionTree(-1, nothing, majority(Y_data))
        end
        children = Dict{Function,DecisionTree}()
        if max_fn === nothing
            # categorical
            X_vec_unique = get_unique_array(X_data[:, max_index])
            for val in X_vec_unique
                fn = (m) -> m == val
                identify = fn.(X_data[:, max_index])
                X_part = X_data[identify, :]
                Y_part = Y_data[identify]
                children[fn] = create_decision_tree_recursive(X_part, Y_part, X_type, max_depth-1)
            end
        else
            # continuous
            identify_1 = max_fn.(X_data[:, max_index])
            X_part_1 = X_data[identify_1, :]
            Y_part_1 = Y_data[identify_1]
            children[max_fn] = create_decision_tree_recursive(X_part_1, Y_part_1, X_type, max_depth-1)
            identify_2 = (!max_fn).(X_data[:, max_index])
            X_part_2 = X_data[identify_2, :]
            Y_part_2 = Y_data[identify_2]
            children[(!max_fn)] = create_decision_tree_recursive(X_part_2, Y_part_2, X_type, max_depth-1)
        end
        return DecisionTree(max_index, children, nothing)
    end

    return create_decision_tree_recursive(X_data, Y_data, detect_types(X_data), max_depth)
end

"""
Prediction function for Decision Tree\\
`X_data` has ndims of 1 or 2\\
Returns an array of predictions\\
Print log when error occurs (undecided children tree)
"""
function predict(X_data::Array{T} where T<:Number, tree::DecisionTree)::Array
    if ndims(X_data) == 1
        X_data = reshape(X_data, (1, size(X_data)[1]))
    end
    @assert ndims(X_data) == 2
    prediction = []
    for i in 1:size(X_data)[1]
        X_vec = X_data[i, :]
        Y_pred = nothing
        tree_copy = tree
        while true
            if tree_copy.target !== nothing
                Y_pred = tree_copy.target
                break
            end
            @assert tree_copy.children !== nothing
            updated = false
            for key_fn in keys(tree_copy.children)
                if key_fn(X_vec[tree_copy.col_id])
                    tree_copy = tree_copy.children[key_fn]
                    updated = true
                    break
                end
            end
            if !updated
                println("Error occured for $X_vec")
                break
            end
        end
        push!(prediction, Y_pred)
    end
    return prediction
end

end

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using OffsetArrays
using .SetProb

function qinit(parm::CParam, meqn::Int, mbc::Int, mx::Int, xlower::Float64, dx::Float64, q::OffsetArray{Float64}, maux::Int, aux::OffsetArray{Float64})

    # Set initial conditions for the q array.
    # This default version prints an error message since it should
    # not be used directly.  Copy this to an application directory and
    # loop over all grid cells to set values of q(1:meqn, 1:mx).

    # integer, intent(in) :: meqn,mbc,mx,maux
    # real(kind=8), intent(in) :: xlower,dx
    # real(kind=8), intent(in) :: aux(maux,1-mbc:mx+mbc)
    # real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc)

    # integer :: i
    # real(kind=8) :: beta, xcell
    # common /cqinit/ beta

    for i = 1 : mx
        xcell = xlower + (i-0.5)*dx
        q[1,i] = exp(-parm.beta * (xcell-0.3)^2)
        q[2,i] = 0.0
        #println("q[1,$i] = $(q[1,i])")
        #println("q[2,$i] = $(q[2,i])")
    end

end # qinit


text:
include("edit_file_global.jl")
simulations = []
include("edit_file_none.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_tgfb.jl")
push!(simulations,rate_constant_update_array)
@show rate_constant_update_array
include("edit_file_mapkin.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_ss.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_ssecad.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_ss_tgfb3ab.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_ss_dnlef1.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_tgfb3.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_dnsmad.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_tgfb3_dnlef1.jl")
push!(simulations,rate_constant_update_array)
include("edit_file_tgfb_tgfb3ab.jl")
push!(simulations,rate_constant_update_array)


text:
module TestConstructors
    using Base.Test
    using DataFrames, DataFrames.Index

    #
    # DataFrame
    #

    df = DataFrame()
    @test isequal(df.columns, Any[])
    @test isequal(df.colindex, Index())

    df = DataFrame(Any[NullableCategoricalVector(zeros(3)),
                       NullableCategoricalVector(ones(3))],
                   Index([:x1, :x2]))
    @test size(df, 1) == 3
    @test size(df, 2) == 2

    @test isequal(df, DataFrame(Any[NullableCategoricalVector(zeros(3)),
                                    NullableCategoricalVector(ones(3))]))
    @test isequal(df, DataFrame(x1 = [0.0, 0.0, 0.0],
                                x2 = [1.0, 1.0, 1.0]))

    df2 = convert(DataFrame, [0.0 1.0;
                              0.0 1.0;
                              0.0 1.0])
    names!(df2, [:x1, :x2])
    @test isequal(df[:x1], NullableArray(df2[:x1]))
    @test isequal(df[:x2], NullableArray(df2[:x2]))

    @test isequal(df, DataFrame(x1 = [0.0, 0.0, 0.0],
                                x2 = [1.0, 1.0, 1.0]))
    @test isequal(df, DataFrame(x1 = [0.0, 0.0, 0.0],
                                x2 = [1.0, 1.0, 1.0],
                                x3 = [2.0, 2.0, 2.0])[[:x1, :x2]])

    df = DataFrame(Int, 2, 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Nullable{Int}, Nullable{Int}]

    df = DataFrame([Int, Float64], [:x1, :x2], 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Nullable{Int}, Nullable{Float64}]

    @test isequal(df, DataFrame([Int, Float64], 2))

end


text:
# Julia wrapper for header: IviDmm.h
# Automatically generated using Clang.jl

function IviDmm_init(logicalName, idQuery, resetDevice, vi)
    ccall((:IviDmm_init, libnidmm), ViStatus, (ViRsrc, ViBoolean, ViBoolean, Ptr{ViSession}), logicalName, idQuery, resetDevice, vi)
end

function IviDmm_close(vi)
    ccall((:IviDmm_close, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_reset(vi)
    ccall((:IviDmm_reset, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_self_test(vi, selfTestResult, selfTestMessage)
    ccall((:IviDmm_self_test, libnidmm), ViStatus, (ViSession, Ptr{ViInt16}, Ptr{ViChar}), vi, selfTestResult, selfTestMessage)
end

function IviDmm_error_query(vi, errorCode, errorMessage)
    ccall((:IviDmm_error_query, libnidmm), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViChar}), vi, errorCode, errorMessage)
end

function IviDmm_error_message(vi, statusCode, message)
    ccall((:IviDmm_error_message, libnidmm), ViStatus, (ViSession, ViStatus, Ptr{ViChar}), vi, statusCode, message)
end

function IviDmm_revision_query(vi, driverRev, instrRev)
    ccall((:IviDmm_revision_query, libnidmm), ViStatus, (ViSession, Ptr{ViChar}, Ptr{ViChar}), vi, driverRev, instrRev)
end

function IviDmm_InvalidateAllAttributes(vi)
    ccall((:IviDmm_InvalidateAllAttributes, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_ResetWithDefaults(vi)
    ccall((:IviDmm_ResetWithDefaults, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_Disable(vi)
    ccall((:IviDmm_Disable, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_InitWithOptions(logicalName, IDQuery, resetDevice, optionString, vi)
    ccall((:IviDmm_InitWithOptions, libnidmm), ViStatus, (ViRsrc, ViBoolean, ViBoolean, ViConstString, Ptr{ViSession}), logicalName, IDQuery, resetDevice, optionString, vi)
end

function IviDmm_GetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviDmm_GetAttributeViInt32, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt32}), vi, channelName, attributeId, value)
end

function IviDmm_SetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViInt32, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViInt32, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviDmm_GetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviDmm_GetAttributeViInt64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt64}), vi, channelName, attributeId, value)
end

function IviDmm_SetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViInt64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViInt64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviDmm_GetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviDmm_GetAttributeViReal64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViReal64}), vi, channelName, attributeId, value)
end

function IviDmm_SetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViReal64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViReal64, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviDmm_GetAttributeViString(vi, channelName, attributeId, bufferSize, value)
    ccall((:IviDmm_GetAttributeViString, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32, Ptr{ViChar}), vi, channelName, attributeId, bufferSize, value)
end

function IviDmm_SetAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViString, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViString, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviDmm_GetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviDmm_GetAttributeViBoolean, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViBoolean}), vi, channelName, attributeId, value)
end

function IviDmm_SetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViBoolean, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViBoolean, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviDmm_GetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviDmm_GetAttributeViSession, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViSession}), vi, channelName, attributeId, value)
end

function IviDmm_SetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviDmm_SetAttributeViSession, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviDmm_CheckAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviDmm_CheckAttributeViSession, libnidmm), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviDmm_LockSession(vi, callerHasLock)
    ccall((:IviDmm_LockSession, libnidmm), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviDmm_UnlockSession(vi, callerHasLock)
    ccall((:IviDmm_UnlockSession, libnidmm), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviDmm_GetError(vi, errorCode, bufferSize, description)
    ccall((:IviDmm_GetError, libnidmm), ViStatus, (ViSession, Ptr{ViStatus}, ViInt32, Ptr{ViChar}), vi, errorCode, bufferSize, description)
end

function IviDmm_ClearError(vi)
    ccall((:IviDmm_ClearError, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_GetNextInterchangeWarning(vi, bufferSize, warning)
    ccall((:IviDmm_GetNextInterchangeWarning, libnidmm), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, warning)
end

function IviDmm_ClearInterchangeWarnings(vi)
    ccall((:IviDmm_ClearInterchangeWarnings, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_ResetInterchangeCheck(vi)
    ccall((:IviDmm_ResetInterchangeCheck, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_GetNextCoercionRecord(vi, bufferSize, record)
    ccall((:IviDmm_GetNextCoercionRecord, libnidmm), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, record)
end

function IviDmm_GetSpecificDriverCHandle(vi, specificDriverCHandle)
    ccall((:IviDmm_GetSpecificDriverCHandle, libnidmm), ViStatus, (ViSession, Ptr{ViSession}), vi, specificDriverCHandle)
end

function IviDmm_GetSpecificDriverIUnknownPtr(vi, specificDriverIUnknownPtr)
    ccall((:IviDmm_GetSpecificDriverIUnknownPtr, libnidmm), ViStatus, (ViSession, Ptr{Cvoid}), vi, specificDriverIUnknownPtr)
end

function IviDmm_ConfigureMeasurement(vi, _function, range, resolution)
    ccall((:IviDmm_ConfigureMeasurement, libnidmm), ViStatus, (ViSession, ViInt32, ViReal64, ViReal64), vi, _function, range, resolution)
end

function IviDmm_ConfigureTrigger(vi, triggerSource, triggerDelay)
    ccall((:IviDmm_ConfigureTrigger, libnidmm), ViStatus, (ViSession, ViInt32, ViReal64), vi, triggerSource, triggerDelay)
end

function IviDmm_Read(vi, maxTime, reading)
    ccall((:IviDmm_Read, libnidmm), ViStatus, (ViSession, ViInt32, Ptr{ViReal64}), vi, maxTime, reading)
end

function IviDmm_Fetch(vi, maxTime, reading)
    ccall((:IviDmm_Fetch, libnidmm), ViStatus, (ViSession, ViInt32, Ptr{ViReal64}), vi, maxTime, reading)
end

function IviDmm_Abort(vi)
    ccall((:IviDmm_Abort, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_Initiate(vi)
    ccall((:IviDmm_Initiate, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_IsOverRange(vi, measurementValue, isOverRange)
    ccall((:IviDmm_IsOverRange, libnidmm), ViStatus, (ViSession, ViReal64, Ptr{ViBoolean}), vi, measurementValue, isOverRange)
end

function IviDmm_ConfigureACBandwidth(vi, minFreq, maxFreq)
    ccall((:IviDmm_ConfigureACBandwidth, libnidmm), ViStatus, (ViSession, ViReal64, ViReal64), vi, minFreq, maxFreq)
end

function IviDmm_ConfigureFrequencyVoltageRange(vi, frequencyVoltageRange)
    ccall((:IviDmm_ConfigureFrequencyVoltageRange, libnidmm), ViStatus, (ViSession, ViReal64), vi, frequencyVoltageRange)
end

function IviDmm_ConfigureTransducerType(vi, transducerType)
    ccall((:IviDmm_ConfigureTransducerType, libnidmm), ViStatus, (ViSession, ViInt32), vi, transducerType)
end

function IviDmm_ConfigureFixedRefJunction(vi, fixedRefJunction)
    ccall((:IviDmm_ConfigureFixedRefJunction, libnidmm), ViStatus, (ViSession, ViReal64), vi, fixedRefJunction)
end

function IviDmm_ConfigureThermocouple(vi, thermocoupleType, refJunctionType)
    ccall((:IviDmm_ConfigureThermocouple, libnidmm), ViStatus, (ViSession, ViInt32, ViInt32), vi, thermocoupleType, refJunctionType)
end

function IviDmm_ConfigureRTD(vi, alpha, resistance)
    ccall((:IviDmm_ConfigureRTD, libnidmm), ViStatus, (ViSession, ViReal64, ViReal64), vi, alpha, resistance)
end

function IviDmm_ConfigureThermistor(vi, resistance)
    ccall((:IviDmm_ConfigureThermistor, libnidmm), ViStatus, (ViSession, ViReal64), vi, resistance)
end

function IviDmm_ConfigureMeasCompleteDest(vi, measCompleteDest)
    ccall((:IviDmm_ConfigureMeasCompleteDest, libnidmm), ViStatus, (ViSession, ViInt32), vi, measCompleteDest)
end

function IviDmm_ConfigureMultiPoint(vi, triggerCount, sampleCount, sampleTrigger, sampleInterval)
    ccall((:IviDmm_ConfigureMultiPoint, libnidmm), ViStatus, (ViSession, ViInt32, ViInt32, ViInt32, ViReal64), vi, triggerCount, sampleCount, sampleTrigger, sampleInterval)
end

function IviDmm_ReadMultiPoint(vi, maxTime, arraySize, readingArray, actualPts)
    ccall((:IviDmm_ReadMultiPoint, libnidmm), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViReal64}, Ptr{ViInt32}), vi, maxTime, arraySize, readingArray, actualPts)
end

function IviDmm_FetchMultiPoint(vi, maxTime, arraySize, readingArray, actualPts)
    ccall((:IviDmm_FetchMultiPoint, libnidmm), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViReal64}, Ptr{ViInt32}), vi, maxTime, arraySize, readingArray, actualPts)
end

function IviDmm_ConfigureTriggerSlope(vi, polarity)
    ccall((:IviDmm_ConfigureTriggerSlope, libnidmm), ViStatus, (ViSession, ViInt32), vi, polarity)
end

function IviDmm_SendSoftwareTrigger(vi)
    ccall((:IviDmm_SendSoftwareTrigger, libnidmm), ViStatus, (ViSession,), vi)
end

function IviDmm_GetApertureTimeInfo(vi, apertureTime, apertureTimeUnits)
    ccall((:IviDmm_GetApertureTimeInfo, libnidmm), ViStatus, (ViSession, Ptr{ViReal64}, Ptr{ViInt32}), vi, apertureTime, apertureTimeUnits)
end

function IviDmm_GetAutoRangeValue(vi, autoRangeValue)
    ccall((:IviDmm_GetAutoRangeValue, libnidmm), ViStatus, (ViSession, Ptr{ViReal64}), vi, autoRangeValue)
end

function IviDmm_ConfigureAutoZeroMode(vi, autoZeroMode)
    ccall((:IviDmm_ConfigureAutoZeroMode, libnidmm), ViStatus, (ViSession, ViInt32), vi, autoZeroMode)
end

function IviDmm_ConfigurePowerLineFrequency(vi, powerLineFreq)
    ccall((:IviDmm_ConfigurePowerLineFrequency, libnidmm), ViStatus, (ViSession, ViReal64), vi, powerLineFreq)
end


text:
"""
    _init_gibbs_sample(bn, evidence)

A random sample of non-evidence nodes uniformly over their domain
"""
@inline function _init_gibbs_sample(bn::BayesNet, evidence::Assignment=Assignment())
    sample = Assignment()

    for cpd in bn.cpds
        nn = name(cpd)
        if haskey(evidence, nn)
            sample[nn] = evidence[nn]
        else
            # pick a random variable from that distribution
            sample[nn] = rand(1:ncategories(bn, nn))
        end
    end

    return sample
end

"""Get the cpd's of a node and its children"""
get_mb_cpds(bn::BayesNet, node::Symbol) = [get(bn, n) for n in vcat(children(bn, node), node)]

"""
    eval_mb_cpd(node, ncategories, assignment, mb_cpds)

Return the potential of all instances of a node given its markove blanket
as a WeightVec:
    P(node | pa_node) * Prod (c in children) P(c | pa_c)

Trys out all possible values of node (assumes categorical)
Assignment should have values for all in the Markov blanket, including the
variable itself.
"""
@inline function eval_mb_cpd(node::Symbol, ncategories::Int,
        assignment::Assignment, mb_cpds::Vector)
    # pdf accepts Assignment only, so swap out current for new value
    old_value = assignment[node]
    p = zeros(ncategories)

    for x in 1:ncategories
        assignment[node] = x

        p[x] = foldl((*), (pdf(cpd, assignment) for cpd in mb_cpds))
    end

    assignment[node] = old_value
    return Weights(p)
end

"""
    infer(GibbsSampling, state::Assignment, InferenceState)

Run Gibbs sampling for `N` iterations. Each iteration changes one node.

Discareds first `burn_in` samples and keeps only the `thin`-th sample.
Ex, if `thin=3`, will discard the first two samples and keep the third.
"""
@with_kw type GibbsSamplingNodewise <: InferenceMethod
    nsamples::Int=2E3
    burn_in::Int=500
    thin::Int=3
    state::Assignment = Assignment()
end
function infer{BN<:DiscreteBayesNet}(im::GibbsSamplingNodewise, inf::InferenceState{BN})

    nsamples, burn_in, thin, x = im.nsamples, im.burn_in, im.thin, im.state

    burn_in < nsamples || throw(ArgumentError("Burn in ($burn_in) " *
                "must be less than number of samples ($nsamples)"))

    # if the math doesn't work out correctly, loop a couple more times ...
    total_num_samples = Int(ceil((nsamples-burn_in) / thin))

    bn = inf.pgm
    nodes = names(bn)
    query = inf.query
    evidence = inf.evidence
    non_evidence = setdiff(nodes, keys(evidence))

    # the current state
    if isempty(im.state)
        im.state = _init_gibbs_sample(bn, inf.evidence)
    else
        _ensure_query_nodes_in_bn_and_not_in_evidence(query, names(bn), evidence)
    end
    x = im.state

    n_cats = Dict(n => ncategories(bn, n) for n in non_evidence)
    # if each non-evidence node is a query variable, and its order as a query
    q_loc = indexin(non_evidence, query)

    ϕ = Factor(query, [n_cats[q] for q in query])
    # manual index into factor.potential
    q_ind = [x[q] for q in query]

    # Markov blankets of each node to sample from
    mb_cpds = Dict(n => get_mb_cpds(bn, n) for n = non_evidence)

    finished = false
    after_burn = false
    num_smpls = 0
    num_iters = 0

    # a sample is samples after sampling just one node
    while !finished
        for (i, n) in enumerate(non_evidence)
            # potential of each instance of n
            wv = eval_mb_cpd(n, n_cats[n], x, mb_cpds[n])
            # sample x_n ~ P(X_n|mb(X))
            x[n] = sample(wv)

            # changing one variable at a time, so can just update that node
            qi = q_loc[i]
            # if it is a query variable, update the index
            qi != 0 && (q_ind[qi] = x[n])

            num_iters += 1

            # kick in the afterburners
            if !after_burn && num_iters > burn_in
                after_burn = true
            end

            # start collecting after the burn in and on the `thin`-th iteration
            if after_burn && ( ((num_iters - burn_in) % thin) == 0 )
                num_smpls += 1

                # sample
                @inbounds ϕ.potential[q_ind...] += 1
            end

            # doubly nested loops!! yay!!
            if num_smpls >= total_num_samples
                finished = true
                break
            end
        end
    end

    normalize!(ϕ)
    return ϕ
end

"""
    infer(im, inf)

Run Gibbs sampling for `N` iterations. Each iteration changes all
nodes.
Discareds first `burn_in` samples and keeps only the `thin`-th sample.
Ex, if `thin=3`, will discard the first two samples and keep the third.
"""
@with_kw type GibbsSamplingFull <: InferenceMethod
    nsamples::Int=2E3
    burn_in::Int=500
    thin::Int=3
    state::Assignment=Assignment()
end
function infer{BN<:DiscreteBayesNet}(im::GibbsSamplingFull, inf::InferenceState{BN})

    nsamples, burn_in, thin = im.nsamples, im.burn_in, im.thin

    burn_in < nsamples || throw(ArgumentError("Burn in ($burn_in) " *
                "must be less than number of samples ($nsamples)"))

    # if the math doesn't work out correctly, loop a couple more times ...
    total_num_samples = Int(ceil((nsamples-burn_in) / thin))

    bn = inf.pgm
    nodes = names(bn)
    query = inf.query
    evidence = inf.evidence
    non_evidence = setdiff(nodes, keys(evidence))

    n_cats = Dict(n => ncategories(bn, n) for n in non_evidence)
    # if each non-evidence node is a query variable, and its order as a query
    q_loc = indexin(non_evidence, query)

    # the current state
    if isempty(im.state)
        im.state = _init_gibbs_sample(bn, inf.evidence)
    else
        _ensure_query_nodes_in_bn_and_not_in_evidence(query, names(bn), evidence)
    end
    x = im.state

    ϕ = Factor(query, [n_cats[q] for q in query])
    # manual index into factor.potential
    q_ind = [x[q] for q in query]

    # Markov blankets of each node to sample from
    mb_cpds = Dict(n => get_mb_cpds(bn, n) for n = non_evidence)

    finished = false
    after_burn = false
    num_smpls = 0
    num_iters = 0

    while !finished
        # generate a new sample
        for (i, n) in enumerate(non_evidence)
            # potential of each instance of n
            wv = eval_mb_cpd(n, n_cats[n], x, mb_cpds[n])
            # sample x_n ~ P(X_n|mb(X))
            x[n] = sample(wv)

            # check if `n` is a query node, and update the index
            qi = q_loc[i]
            qi != 0 && (q_ind[qi] = x[n])
        end

        num_iters += 1

        # kick in the afterburners
        if !after_burn && num_iters > burn_in
            after_burn = true
        end

        # start collecting after the burn in and on the `thin`-th iteration
        if after_burn && ( ((num_iters - burn_in) % thin) == 0 )
            num_smpls += 1
            @inbounds ϕ.potential[q_ind...] += 1
        end

        if num_smpls >= total_num_samples
            finished = true
        end
    end

    normalize!(ϕ)
    return ϕ
end


text:
# | model                  | build  | fitted_params | report | metadata | tests 1 | tests 2 |
# | ---------------------- | ------ | ------------- | ------ | -------- | ------- | ------- |
# | GPRegressor            | ✓      | ✓             | ✗      | ✓        |  ✓      | ✓       |
# | GPClassif              | ✗      | ✗             | ✗      | ✗        |  ✗      | ✗       |

GaussianProcessRegressor_ = SKGP.GaussianProcessRegressor
@sk_reg mutable struct GaussianProcessRegressor <: MLJBase.Deterministic
    kernel::Any        = nothing
    alpha::Union{Float64,AbstractArray} = 1.0e-10
    optimizer::Any     = "fmin_l_bfgs_b"
    n_restarts_optimizer::Int = 0
    normalize_y::Bool  = false
    copy_X_train::Bool = true
    random_state::Any  = nothing
end
MLJBase.fitted_params(model::GaussianProcessRegressor, (fitresult, _, _)) = (
    X_train = fitresult.X_train_,
    y_train = fitresult.y_train_,
    kernel  = fitresult.kernel_,
    L       = fitresult.L_,
    alpha   = fitresult.alpha_,
    log_marginal_likelihood_value = fitresult.log_marginal_likelihood_value_
    )

metadata_model(GaussianProcessRegressor,
    input=MLJBase.Table(MLJBase.Continuous),
    target=AbstractVector{MLJBase.Continuous},
    weights=false,
    descr="Gaussian process regressor."
    )

# ============================================================================
GaussianProcessClassifier_ = SKGP.GaussianProcessClassifier
@sk_clf mutable struct GaussianProcessClassifier <: MLJBase.Probabilistic
    kernel::Any        = nothing
    optimizer::Any     = "fmin_l_bfgs_b"
    n_restarts_optimizer::Int = 0
    copy_X_train::Bool = true
    random_state::Any  = nothing
    max_iter_predict::Int = 100::(_ > 0)
    warm_start::Bool   = false
    multi_class::String = "one_vs_rest"::(_ in ("one_vs_one", "one_vs_rest"))
end
MLJBase.fitted_params(m::GaussianProcessClassifier, (f, _, _)) = (
    kernel    = f.kernel_,
    log_marginal_likelihood_value = f.log_marginal_likelihood_value_,
    classes   = f.classes_,
    n_classes = f.n_classes_
    )
metadata_model(GaussianProcessClassifier,
    input=MLJBase.Table(MLJBase.Continuous),
    target=AbstractVector{<:MLJBase.Finite},
    weights=false,
    descr="Gaussian process classifier."
    )


text:
function add_input(obj::AbstractListCrawler, arg0::Object)
    return jcall(obj, "addInput", void, (Object,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function feed(obj::AbstractListCrawler, arg0::Pattern, arg1::DataLoader, arg2::DataProvidersManager)
    return jcall(obj, "feed", jboolean, (Pattern, DataLoader, DataProvidersManager), arg0, arg1, arg2)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_inputs(obj::AbstractListCrawler)
    return jcall(obj, "getInputs", List, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


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

# %% [markdown] tags=[]
# Optimization by specialization to argument types comes with the tradeoff of increased compilation time. An extreme case of this is shown in the following example.

# %%
f(x) = x^2 + 3x + 2
g(::Val{x}) where x = f(x)
f(10^6), g(Val(10^6))

# %%
using BenchmarkTools
@btime f(10^6)
@btime g($(Val(10^6)));

# %% [markdown]
# `g(Val(10^6))` is ultra fast because it is specialized to the argument type `Val{10^6}` and compiled to `return 1000003000002`.

# %%
@code_typed debuginfo=:none g(Val(10^6))

# %% [markdown]
# However, `g(Val(k))` is compiled separately for each different `k`.  So if you run `g(Val(k))` for a large number of different `k`'s,  it will perform a large number of compilations and will be very slow. (After compilation, though, it will be explosively fast.)

# %%
F(n) = [f(k) for k in 1:n]
G(n) = [g(Val(k)) for k in 1:n]

@time F(10^4)
@time G(10^4);

# %% [markdown]
# The first execution of `G(10^4)` is very slow.

# %%
F(10^4) == G(10^4)

# %% [markdown]
# Thus, it is not reasonable to try to optimize by specialization to the argument types by different large numbers of `Val{k}` types.
#
# On the other hand, the native code specialized to argument types is very fast, so if compilation time is not an issue, optimization by specialization to argument types should be done aggressively.
#
# In short, it is a matter of trade-off.

# %%

max_stars_repo_path: 0010/dispatch by Val.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0010/dispatch by Val.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0010/dispatch by Val.jl

max_forks_repo_name: genkuroki/public


text:
module MeshesRecipes

using Meshes
using AbstractPlotting

import AbstractPlotting: Plot, default_theme, plot!

include("polytypes.jl")

end


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

const BitCount = Union{Signed, Unsigned}

@inline fbc(::Type{T}, x) where {T} = x%unsigned(T)
@inline fbc(::Type{T}, bits::UnitRange{U}) where {T,U} = fbc(T, bits.start):fbc(T, bits.stop)

function bsizeof end
export bsizeof

"""
    bsizeof(x)

Returns the data size of x in bits.
"""
@inline bsizeof(x) = sizeof(x) << 3

function bmask end
export bmask

"""
    bmask(::Type{T}, bit::Integer)::T where {T<:Integer}

Generates a bit mask of type `T` with bit `bit` set to one.
"""
@inline bmask(::Type{T}, bit::BitCount) where {T<:Integer} = one(T) << fbc(T, bit)

"""
    bmask(::Type{T}, bits::UnitRange{<:Integer})::T where {T<:Integer}

Generates a bit mask of type `T` with bit-range `bits` set to one.
"""
@inline bmask(::Type{T}, bits::UnitRange{<:BitCount}) where {T<:Integer} = begin
    fbits = fbc(T, bits)
    #@assert fbits.stop >= fbits.start "Bitmask range of fbits can't be reverse"
    ((one(T) << (fbits.stop - fbits.start + 1)) - one(T)) << fbits.start
end

function lsbmask end
export lsbmask

"""
    lsbmask(::Type{T})::T where {T<:Integer}

Generates a bit mask with only the least significant bit set.
"""
@inline lsbmask(::Type{T}) where {T<:Integer} = one(T)

"""
    lsbmask(::Type{T}, nbits::Integer)::T where {T<:Integer}

Generates a bit mask with only the `nbits` least significant bits set.
"""
@inline lsbmask(::Type{T}, nbits::BitCount) where {T<:Integer} = ~(~zero(T) << fbc(T, nbits))

function msbmask end
export msbmask

"""
    msbmask(::Type{T})::T where {T<:Integer}

Generates a bit mask with only the most significant bit set.
"""
@inline msbmask(::Type{T}) where {T<:Integer} = one(T) << fbc(T, bsizeof(T) - 1)

"""
    msbmask(::Type{T}, nbits::Integer)::T where {T<:Integer}

Generates a bit mask with only the `nbits` most significant bits set.
"""
@inline msbmask(::Type{T}, nbits::BitCount) where {T<:Integer} = ~(~zero(T) >>> fbc(T, nbits))

function bget end
export bget

"""
    bget(x::T, bit::Integer)::Bool where {T<:Integer}

Get the value of bit `bit` of x.
"""
@inline bget(x::T, bit::BitCount) where {T<:Integer} = x & bmask(T, fbc(T, bit)) != zero(T)

"""
    bget(x::T, bits::UnitRange{<:Integer})::T where {T<:Integer}

Get the value of the bit range `bits` of x.
"""
@inline bget(x::T, bits::UnitRange{<:BitCount}) where {T<:Integer} = begin
    fbits = fbc(T, bits)
    (x & bmask(T, fbits)) >>> fbits.start
end

function bset end
export bset

"""
    bset(x::T, bit::Integer)::T where {T<:Integer}

Returns a modified copy of x, with bit `bit` set (to one).
"""
@inline bset(x::T, bit::BitCount) where {T<:Integer} = x | bmask(T, fbc(T, bit))

"""
    bset(x::T, bit::Integer, y::Bool)::T where {T<:Integer}

Returns a modified copy of x, with bit `bit` set to `y`.
"""
@inline bset(x::T, bit::BitCount, y::Bool) where {T<:Integer} = y ? bset(x, fbc(T, bit)) : bclear(x, fbc(T, bit))

"""
    bset(x::T, bits::UnitRange{<:Integer}, y::Integer)::T where {T<:Integer}

Returns a modified copy of x, with bit range `bits` set to `y`.
"""
@inline bset(x::T, bits::UnitRange{<:BitCount}, y::Integer) where {T<:Integer} = begin
    fbits = fbc(T, bits)
    local bm = bmask(T, fbits)
    (x & ~bm) | ((convert(T, y) << fbits.start) & bm)
end

function bclear end
export bclear

"""
    bclear(x::T, bit::Integer)::T where {T<:Integer}

Returns a modified copy of x, with bit `bit` cleared (set to zero).
"""
@inline bclear(x::T, bit::BitCount) where {T<:Integer} = x & ~bmask(T, fbc(T, bit))

function bflip end
export bflip

"""
    bflip(x::T, bit::Integer)::T where {T<:Integer}

Returns a modified copy of x, with bit `bit` flipped.
"""
@inline bflip(x::T, bit::BitCount) where {T<:Integer} = xor(x, bmask(T, fbc(T, bit)))

"""
    bflip(x::T, bits::UnitRange{<:Integer})::T where {T<:Integer}

Returns a modified copy of x, with all bits in bit range `bits` flipped.
"""
@inline bflip(x::T, bits::UnitRange{<:BitCount}) where {T<:Integer} = xor(x, bmask(T, fbc(T, bits)))

function lsbget end
export lsbget

"""
    lsbget(x::Integer)::Bool

Returns the value of the least significant bit of x.
"""
@inline lsbget(x::T) where {T<:Integer} =
    x & lsbmask(T) != zero(T)

"""
    lsbget(x::T)::T where {T <: Integer}

Returns the value of the `nbits` least significant bits of x.
"""
@inline lsbget(x::T, nbits::BitCount) where {T<:Integer} =
    x & lsbmask(T, fbc(T, nbits))

function msbget end
export msbget

"""
    msbget(x::T)::Bool where {T<:Integer}

Returns the value of the most significant bit of x.
"""
@inline msbget(x::T) where {T<:Integer} =
    x & msbmask(T) != zero(T)

"""
    msbget(x::T)::T where {T<:Integer}

Returns the value of the `nbits` most significant bits of x.
"""
@inline msbget(x::T, nbits::BitCount) where {T<:Integer} =
    x >>> (bsizeof(x) - fbc(T, nbits))


text:
using LinearAlgebra

@enum Orientation::Int8 Colinear = 0 Clockwise = -1 CounterClockwise = 1

function mid_point(p1, p2)
    return (p1 .+ p2) ./ 2
end

function orientation(p1, p2, p3) # orientation
    val = (p2[2] - p1[2]) * (p3[1] - p2[1]) - (p2[1] - p1[1]) * (p3[2] - p2[2])
    return val == 0.0 ? Colinear : (val < 0 ? CounterClockwise : Clockwise)
end

# make edges and dual edges always be in CounterClockwise orientation
function rectify_orientation!(de, dv, v)
    p1 = [0, 0]
    p2 = v[2] .- v[1]
    p3 = dv[2] .- dv[1]
    or = orientation(p1, p2, p3)
    if or == Clockwise
        reverse!(de)
    elseif or == Colinear
        error("Found three colinear vertices trying to sneak as a triangle! 😠")
    end
end

function rectify_orientation!(triangles, vertices)
    for tri in triangles
        or = orientation(vertices[tri[1:3]]...)
        if or == Clockwise
            reverse!(tri)
        elseif or == Colinear
            error("Found three colinear vertices trying to sneak as a triangle! 😠")
        end
    end
end

function circumcenter(p1, p2, p3) # circumcenter
    mat = SMatrix{2}((p1[1] - p2[1]), (p1[1] - p3[1]), (p1[2] - p2[2]), (p1[2] - p3[2]))
    b = -SVector((p2' * p2 - p1' * p1) / 2, (p3' * p3 - p1' * p1) / 2)
    x = mat \ b
end

function longitude(p1, p2) # measure of a line
    return norm(p1 .- p2)
end

function area(p) # measure of a polygon
    A = 0.0
    for i = 2:length(p)
        A += p[i - 1][1] * p[i][2] - p[i - 1][2] * p[i][1]
    end
    A += p[end][1] * p[1][2] - p[end][2] * p[1][1]
    return abs(A) / 2
end


text:
n = 8 # rand(2: 10)

# test Hessenberg

A = matrixdepot("rohess", n)

B = tril(A) - diagm(0 => diag(A)) - diagm(-1 => diag(A, -1))

@test B ≈ zeros(n,n) atol=1.0e-7

# test orithogonality

@test A'*A ≈ Matrix(1.0I, n, n) atol=1.0e-7

println("'rohess' passed test...")


text:
module Intersect

using ..Geometry
using ..Diagram

mutable struct Rectangle
   w::Real
   h::Real
end

#mutable struct CPolygonVertex
#    origin::Tuple{Real, Real}
#    next::Union{PolygonVertex, Nothing}
#end
#
## curvilinear polygon
#mutable struct CPolygon
#    head::Union{PolygonVertex, Nothing}
#end

function inbounds(box::Rectangle, p::Tuple{Real, Real})
   return (0 <= p[1] <= box.w) && (0 <= p[2] <= box.h)
end

function generateValidPoints(agen::Tuple{Real, Real}, bgen::Tuple{Real, Real})
   center = ((agen[1] + bgen[1])/2, (agen[2] + bgen[2])/2)

   dir_gen = Geometry.rotateVectorCCW(Geometry.subVector(bgen, agen))

   b = Geometry.addVector(center, dir_gen)
   a = Geometry.subVector(center, dir_gen)

   return a, b
end

function intersectLine(box::Rectangle, agen::Tuple{Real, Real}, bgen::Tuple{Real, Real})
    a, b = generateValidPoints(agen, bgen)

   dir = Geometry.subVector(b, a)

   if dir[1] == 0
      if !(0 <= a[1] <= box.w) # no intersection
         return nothing, nothing
      end

      if a[2] < b[2]
         return (a[1], 0), (a[1], box.h)
      else
         return (a[1], box.h), (a[1], 0)
      end
   elseif dir[2] == 0
      if !(0 <= a[2] <= box.h) # no intersection
         return nothing, nothing
      end

      if a[1] < b[1]
         return (0, a[2]), (box.w, a[2])
      else
         return (box.w, a[2]), (0, a[2])
      end
   else
      # intersection points with each line that delimits the box
      l = (0, a[2] - (a[1])*(dir[2]/dir[1]))
      d = (a[1] - (a[2])*(dir[1]/dir[2]), 0)
      r = (box.w, a[2] - (a[1]-box.w)*(dir[2]/dir[1]))
      u = (a[1] - (a[2]-box.h)*(dir[1]/dir[2]), box.h)

      p = Array{Union{Tuple{Real, Real}, Nothing}, 1}([nothing, nothing])
      i = 1
      for q in [l, r, d, u]
         if inbounds(box, q)
            p[i] = q
            i += 1

            if i > 2
                break
            end
         end
      end

      if p[1] == nothing
          if p[2] == nothing
              return nothing, nothing
          else
              return p[2], p[2]
          end
      elseif p[2] == nothing
          return p[1], p[1]
      end

      if (a[1] < b[1] && p[1] < p[2]) || (a[1] > b[1] && p[1] > p[2])
         return p[1], p[2]
      else
         return p[2], p[1]
      end
   end
end

function norm2(v)
    return v[1]^2 + v[2]^2
end

function intersectEdge(box::Rectangle, he::Diagram.HalfEdge)
   if he.origin == nothing
      return nothing, nothing # no intersection
   end

   cp1, cp2 = intersectLine(box, he.generator, he.twin.generator)

   if cp1 == nothing
      return nothing, nothing # no intersection
   end

   if !inbounds(box, he.origin) && !inbounds(box, he.twin.origin)
       if he.isFixed || he.twin.isFixed
           if (he.origin[1] < 0 && he.twin.origin[1] < 0) ||
                (he.origin[1] > box.w && he.twin.origin[1] > box.w) ||
                (he.origin[2] < 0 && he.twin.origin[2] < 0) ||
                (he.origin[2] > box.h && he.twin.origin[2] > box.h) # same side

                if he.isFixed && he.twin.isFixed
                    return nothing, nothing
                end

                if (he.isFixed && !he.twin.isFixed)
                    fixedOne = he
                    unfixedOne = he.twin
                else
                    fixedOne = he.twin
                    unfixedOne = he
                end

                fixedOneDistance = norm2(Geometry.subVector(fixedOne.origin, cp1))
                unfixedOneDistance = norm2(Geometry.subVector(unfixedOne.origin, cp1))

                if fixedOneDistance < unfixedOneDistance
                    return nothing, nothing
                end
            end
       end
   end

   if he.isFixed && inbounds(box, he.origin)
       ret1 = he.origin
   else
       ret1 = cp1
   end

   if he.twin.isFixed && inbounds(box, he.twin.origin)
       ret2 = he.twin.origin
   else
       ret2 = cp2
   end

   return ret1, ret2
end

function nextCorner(box::Rectangle, p::Tuple{Real, Real})
   if p[1] == 0 && p[2] > 0
      return (0, 0)
   elseif p[2] == 0 && p[1] < box.w
      return (box.w, 0)
   elseif p[1] == box.w && p[2] < box.h
      return (box.w, box.h)
   else
      return (0, box.h)
   end
end

function joinAround!(array::Array{Diagram.HalfEdge}, box::Rectangle, a::Diagram.HalfEdge, b::Diagram.HalfEdge)
   ao = a.origin
   bo = b.origin

   if (ao[1] == bo[1] == 0 && ao[2] > bo[2]) ||
      (ao[1] == bo[1] == box.w && ao[2] < bo[2]) ||
      (ao[2] == bo[2] == 0 && ao[1] < bo[1]) ||
      (ao[2] == bo[2] == box.h && ao[1] > bo[1])
      Diagram.concat(a, b)
   else
      new = Diagram.HalfEdge(nextCorner(box, ao), false, nothing, nothing, nothing, a.origin)
      push!(array, new)
      Diagram.concat(a, new)
      return joinAround!(array, box, new, b)
   end
end

function intersect(V::Diagram.DCEL, box::Rectangle)
   for region in V.regions
      he = region.borderHead

      # find the "leftmost" edge in the region, or determine that it's closed
      while he.isFixed && he.prev != region.borderHead
         he = he.prev
      end

      if he.prev == region.borderHead
         he = he.prev
         he = he.prev
      end

      region.borderHead = he

      while he.next != nothing && he.next != region.borderHead
          he = he.next
      end
      he = region.borderHead

      if !he.isFixed && !he.twin.isFixed # unlimited line
         p1, p2 = intersectLine(box, he.generator, he.twin.generator)
         he.origin = p1
         he.twin.origin = p2
         new = Diagram.HalfEdge(p2, false, nothing, nothing, nothing, he.origin)
         push!(V.halfEdges, new)
         Diagram.concat(he, new)
         joinAround!(V.halfEdges, box, new, he)
      else
         while true
            he.origin, he.twin.origin = intersectEdge(box, he)

            if he.next == nothing
               he.next = region.borderHead
            end

            he = he.next
            if he == region.borderHead
               break
            end
         end

         i = 0
         cycleStart = he
         while he.origin == nothing
            he = he.next
            i += 1

            if he == cycleStart
                readline(stdin)
            end
         end
         region.borderHead = he

         while true
            while he.next.origin == nothing
               he.next = he.next.next
            end

            next = he.next

            if he.twin.origin != he.next.origin
               new = Diagram.HalfEdge(he.twin.origin, false, nothing, nothing, nothing, he.origin)
               push!(V.halfEdges, new)
               Diagram.concat(he, new)
               joinAround!(V.halfEdges, box, new, next)
            end

            he = next
            if he == region.borderHead
               break
            end
         end
      end
   end
end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
using Documenter: DocMeta, HTML, MathJax3, deploydocs, makedocs
using Manopt, Manifolds, Literate, Pluto, PlutoStaticHTML, Pkg
# Load an unregistered package (for now) to update exports of Pluto notebooks

TutorialMenu = Array{Pair{String,String},1}()

#
# Generate Pluto Tutorial HTMLs

# First tutorial with AD
pluto_src_folder = joinpath(@__DIR__, "..", "pluto/")
pluto_output_folder = joinpath(@__DIR__, "src/", "pluto/")
pluto_relative_path = "pluto/"
mkpath(pluto_output_folder)
#
#
# Please do not use the same name as for a(n old) literate Tutorial
pluto_files = [
    "Optimize!",
    "AutomaticDifferentiation",
    "HowToRecord",
    "GeodesicRegression",
    "Bezier",
    "SecondOrderDifference",
    "StochasticGradientDescent",
    "Benchmark",
    "JacobiFields",
]
pluto_titles = [
    "Get started: Optimize!",
    "Use AD in Manopt",
    "How to record values",
    "Do Geodesic regression",
    "Use Bezier Curves",
    "Compute a second order difference",
    "Do stochastic gradient descent",
    "speed up! using `gradF!`",
    "Illustrate Jacobi Fields",
]
# build menu and write files myself - tp set edit url correctly.
for (title, file) in zip(pluto_titles, pluto_files)
    global TutorialMenu
    rendered = build_notebooks( #though not really parallel here
        BuildOptions(
            pluto_src_folder;
            output_format=documenter_output,
            write_files=false,
            use_distributed=false,
        ),
        ["$(file).jl"],
    )
    write(
        pluto_output_folder * file * ".md",
        """
        ```@meta
        EditURL = "$(pluto_src_folder)$(file).jl"
        ```

        $(rendered[1])
        """,
    )
    push!(TutorialMenu, title => joinpath(pluto_relative_path, file * ".md"))
end

generated_path = joinpath(@__DIR__, "src")
base_url = "https://github.com/JuliaManifolds/Manopt.jl/blob/master/"
isdir(generated_path) || mkdir(generated_path)
open(joinpath(generated_path, "contributing.md"), "w") do io
    # Point to source license file
    println(
        io,
        """
        ```@meta
        EditURL = "$(base_url)CONTRIBUTING.md"
        ```
        """,
    )
    # Write the contents out below the meta block
    for line in eachline(joinpath(dirname(@__DIR__), "CONTRIBUTING.md"))
        println(io, line)
    end
end

makedocs(;
    format=HTML(; mathengine=MathJax3(), prettyurls=get(ENV, "CI", nothing) == "true"),
    modules=[Manopt],
    sitename="Manopt.jl",
    pages=[
        "Home" => "index.md",
        "About" => "about.md",
        "How to..." => TutorialMenu,
        "Plans" => [
            "Specify a Solver" => "plans/index.md",
            "Problem" => "plans/problem.md",
            "Options" => "plans/options.md",
            "Stepsize" => "plans/stepsize.md",
        ],
        "Solvers" => [
            "Introduction" => "solvers/index.md",
            "Alternating Gradient Descent" => "solvers/alternating_gradient_descent.md",
            "Chambolle-Pock" => "solvers/ChambollePock.md",
            "Conjugate gradient descent" => "solvers/conjugate_gradient_descent.md",
            "Cyclic Proximal Point" => "solvers/cyclic_proximal_point.md",
            "Douglas–Rachford" => "solvers/DouglasRachford.md",
            "Gradient Descent" => "solvers/gradient_descent.md",
            "Nelder–Mead" => "solvers/NelderMead.md",
            "Particle Swarm Optimization" => "solvers/particle_swarm.md",
            "Quasi-Newton" => "solvers/quasi_Newton.md",
            "Stochastic Gradient Descent" => "solvers/stochastic_gradient_descent.md",
            "Subgradient method" => "solvers/subgradient.md",
            "Steihaug-Toint TCG Method" => "solvers/truncated_conjugate_gradient_descent.md",
            "Trust-Regions Solver" => "solvers/trust_regions.md",
        ],
        "Functions" => [
            "Introduction" => "functions/index.md",
            "Bézier curves" => "functions/bezier.md",
            "Cost functions" => "functions/costs.md",
            "Differentials" => "functions/differentials.md",
            "Adjoint Differentials" => "functions/adjointdifferentials.md",
            "Gradients" => "functions/gradients.md",
            "Jacobi Fields" => "functions/Jacobi_fields.md",
            "Proximal Maps" => "functions/proximal_maps.md",
            "Specific Manifold Functions" => "functions/manifold.md",
        ],
        "Helpers" => [
            "Checks" => "helpers/checks.md",
            "Data" => "helpers/data.md",
            "Error Measures" => "helpers/errorMeasures.md",
            "Exports" => "helpers/exports.md",
        ],
        "Contributing to Manopt.jl" => "contributing.md",
        "Notation" => "notation.md",
        "Function Index" => "list.md",
    ],
)
deploydocs(; repo="github.com/JuliaManifolds/Manopt.jl", push_preview=true)