text:
#CppSimData: A Julia wrapper for CppSimData tools
#-------------------------------------------------------------------------------
__precompile__(true)

module CppSimData

const rootpath = realpath(joinpath(dirname(realpath(@__FILE__)),"../."))
const libroot = joinpath(rootpath, "core/lib")

const objfilename = @static (
if Sys.iswindows()
	joinpath(:x86_64 == Sys.ARCH ? "win64" : "win32", "cppsimdata_lib.dll")
elseif Sys.isapple()
	joinpath("macosx", "cppsimdata_lib.so")
elseif Sys.islinux()
	joinpath(:x86_64 == Sys.ARCH ? "glnxa64" : "glnx86", "cppsimdata_lib.so")
else
	error("No compiled sources for $(Sys.KERNEL) ($(Sys.ARCH))")
end
)
const objfile = joinpath(libroot, objfilename)

include("base.jl")
include("show.jl")

initialize() #Initialize CppSimData module

export loadsig, lssig, evalsig

end #CppSimData
#Last line


text:
using CurveProximityQueries
using ConvexBodyProximityQueries
using Test
using LinearAlgebra
using StaticArrays
using IntervalArithmetic
using Random: seed!

import CurveProximityQueries: differentiate, integrate

@testset "CurveProximityQueries" begin
    @testset "Bernstein Polynomials" begin
        B2 = rand(Bernstein{2,8})
        B3 = rand(Bernstein{3,5})
        @test typeof(B2) <: Bernstein{2,8}
        @test typeof(B3) <: Bernstein{3,5}
        @test sum(norm.((differentiate(integrate(B2.f)) - B2.f).control_points)) ≤ 1e-10
        @test sum(norm.((differentiate(integrate(B3.f)) - B3.f).control_points)) ≤ 1e-10
    end
    @testset "Curve - Polygon" begin
        seed!(1)
        obs = @point zeros(2)
        c = rand(Bernstein{2,7})
        @test abs(minimum_distance(obs, c) - 0.3683522584741768) ≤ 1e-5
        @test tolerance_verification(obs, c, 0.3) == true
        @test tolerance_verification(obs, c, 0.4) == false
        @test collision_detection(obs, c) == false
        @test abs(minimum_distance(c, obs) - 0.3683522584741768) ≤ 1e-5
        @test tolerance_verification(c, obs, 0.3) == true
        @test tolerance_verification(c, obs, 0.4) == false
        @test collision_detection(c, obs) == false

        seed!(1)
        obs = @point zeros(3)
        c = rand(Bernstein{3,11})
        @test abs(minimum_distance(obs, c) - 0.511627527056288) ≤ 1e-5
        @test tolerance_verification(obs, c, 0.5) == true
        @test tolerance_verification(obs, c, 0.6) == false
        @test collision_detection(obs, c) == false
        @test abs(minimum_distance(c, obs) - 0.511627527056288) ≤ 1e-5
        @test tolerance_verification(c, obs, 0.5) == true
        @test tolerance_verification(c, obs, 0.6) == false
        @test collision_detection(c, obs) == false
    end
    @testset "Curve - Curve" begin
        seed!(1)
        c = rand(Bernstein{2,7})
        d = rand(Bernstein{2,3})
        @test abs(minimum_distance(c, d)) ≤ 1e-5
        @test tolerance_verification(c, d, 0.3) == false
        @test collision_detection(c, d) == true

        seed!(1)
        obs = @point zeros(3)
        c = rand(Bernstein{3,11})
        d = rand(Bernstein{3,5})
        @test abs(minimum_distance(c, d) - 0.137873546917387) ≤ 1e-5
        @test tolerance_verification(c, d, 0.1) == true
        @test tolerance_verification(c, d, 0.3) == false
        @test collision_detection(c, d) == false
    end
end

ext: jl


text:
function makecinfo(ecut,Ω,basis="Hermite")
cinfo = Cinfo()
if basis=="Hermite"
  dim=length(Ω)
  if dim==1
    ωx = Ω[1]
    e0 = 0.5*ωx
    ecut < e0 && error("ecut must exceed the zero point energy.")
    Mx,nx,en = nenergy(ecut,e0,ωx,basis)
    Mmax = Mx
    P  = en .< ecut; Mult = sum(P)
    en = P.*en
    M  = [Mx]
  elseif dim==2
    ωx,ωy = Ω
    e0 = 0.5(ωx+ωy)
    ecut < e0 && error("ecut must exceed the zero point energy.")
    Mx,nx,ex = nenergy(ecut,e0,ωx,basis)
    My,ny,ey = nenergy(ecut,e0,ωy,basis)
    M  = [Mx,My]; Mmax = maximum(M)
    en = [ex[i+1]+ey[j+1] for i=nx, j=ny]
    P  = en .< ecut; Mult = sum(P)
    en = P.*en
    M  = [Mx,My]
  elseif dim==3
    ωx,ωy,ωz = Ω
    e0 = 0.5(ωx+ωy+ωz)
    ecut < e0 && error("ecut must exceed the zero point energy.")
    Mx,nx,ex = nenergy(ecut,e0,ωx,basis)
    My,ny,ey = nenergy(ecut,e0,ωy,basis)
    Mz,nz,ez = nenergy(ecut,e0,ωz,basis)
    M  = [Mx,My,Mz]; Mmax = maximum(M)
    en = [ex[i+1]+ey[j+1]+ez[k+1] for i=nx, j=ny, k=nz]
    P  = en .< ecut; Mult = sum(P)
    en = P.*en
else error("Spatial dimension must be 1, 2, or 3.")
  end
end
  @pack cinfo = basis, Ω, ecut, e0, Mmax, M, en, P, Mult
  return cinfo
end


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "quadruped", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

## quadruped on piecewise surface

s = get_simulation("quadruped", "piecewise1_2D_lc", "piecewise", approx = true)
model = s.model
env = s.env
nq = model.dim.q
nc = model.dim.c

ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt))

ref_traj = deepcopy(ref_traj_)

# time
H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 4000 #3000

# barrier parameter
κ_mpc = 1.0e-4

obj = TrackingObjective(model, env, H_mpc,
	q = [Diagonal(1e-2 * [5; 0.02; 0.10; 0.25 * ones(nq-3)]) for t = 1:H_mpc],
    # q = [Diagonal(1e-2 * [10; 0.02; 0.25; 0.25 * ones(nq-3)]) for t = 1:H_mpc],
    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	# mode = :configuration,
	ip_type = :mehrotra,
    n_opts = NewtonOptions(
		solver = :lu_solver,
		r_tol = 3e-4,
		max_iter = 5,
		# max_time = ref_traj.h, # HARD REAL TIME
		),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        altitude_update = true,
        altitude_impact_threshold = 0.05,
        # altitude_verbose = true,
        ),
	ip_opts = InteriorPointOptions(
		max_iter = 100,
		verbose = false,
		r_tol = 1.0e-4,
		κ_tol = 1.0e-4,
		diff_sol = true,
		# κ_reg = 1e-3,
		# γ_reg = 1e-1,
		solver = :empty_solver,
		),
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))

sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = ContactImplicitMPC.InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_tol = 2.0e-6,
        diff_sol = true),
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

@time status = ContactImplicitMPC.simulate!(sim, verbose = true)

# plt = plot(layout=(3,1), legend=false)
# plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
# plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[2,1], hcat(Vector.([u[1:nu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

plot_lines!(vis, model, sim.traj.q[1:1:end])
plot_surface!(vis, env, ylims=[0.3, -0.05])
anim = visualize_meshrobot!(vis, model, sim.traj, sample=5)
# anim = visualize_robot!(vis, model, sim.traj, anim=anim)
anim = visualize_force!(vis, model, sim.traj, anim=anim, h=h_sim)

# # Display ghosts
# t_ghosts = [1, 1333, 2666]
# mvis_ghosts = []
# for (i,t) in enumerate(t_ghosts)
#     α = i/(length(t_ghosts)+1)
#     name = Symbol("ghost$i")
#     mvis = build_meshrobot!(vis, model, name=name, α=α)
#     push!(mvis_ghosts, mvis)
# end
#
# for (i,t) in enumerate(t_ghosts)
#     name = Symbol("ghost$i")
#     set_meshrobot!(vis, mvis_ghosts[i], model, sim.traj.q[t], name=name)
# end

# filename = "quadruped_forces"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)

# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)

# const ContactImplicitMPC = Main


text:
"""
    ScalarBatchAccumulator(net, state, batch_size)

A ScalarBatchAccumulator is used to evaluate the contribution to
local observables ⟨σ|Ô|σ'⟩ψ(σ')/ψ(σ) , but by computing the neural
network ψ(σ) in batches of size `batch_size`. This is essential
to extrract a speedup when using GPUs.

This is an internal implementation detail of NeuralQuantum, and should
not be relied upon.

Once constructed, a ScalarBatchAccumulator is supposed to be used as follows:
- if `isfull(sba) == true` you should not push new elements to it (an error
will be throw otherwise)
- data is pushed as `sba(⟨σ|Ô|σ'⟩, σ', ψ(σ))`.
The configuration should be passed as a vector (if ket) or as a tuple
of two vectors (if density matrix).
"""
mutable struct AccumulatorLogGradPsi{N,A,B,C} <: AbstractMachineGradAccumulator
    bnet::N         # A batched version of the cached neural network

    in_buf::A       # the matrix of Nsites x batchsz used as input
    out_buf::B      # The row vector of outputs
    ∇out_buf::C

    buf_n::Int      # Counter for elements in buffer
    batch_sz::Int   # batch size
end

function AccumulatorLogGradPsi(net::NeuralNetwork, v, batch_sz)
    #bnet     = cached(net, batch_sz)
    bnet     = net

    w        = trainable_first(net)
    RT       = real(eltype(w))
    in_buf   = deepcopy(v)
    out_buf  = out_similar(net)
    ∇out_buf = grad_cache(net, batch_sz)

    return AccumulatorLogGradPsi(
        bnet,
        in_buf, out_buf, ∇out_buf,
        0, batch_sz)
end

Base.length(a::AccumulatorLogGradPsi) = a.batch_sz
Base.count(a::AccumulatorLogGradPsi)  = a.buf_n
isfull(a::AccumulatorLogGradPsi) = count(a) == length(a)

"""
    init!(c::ScalarBatchAccumulator)

Resets the internal counter of the accumulator, deleting
all previously accumulated (but not computed) values.
"""
init!(c::AccumulatorLogGradPsi) = c.buf_n = 0
function init!(c::AccumulatorLogGradPsi, σ)
    init!(c)

    # Reset the state so that later is faster to apply the changes
    state_copy!(c.in_buf, σ)
    return nothing
end

function (c::AccumulatorLogGradPsi)(v::Union{AState, ADoubleState})
    @assert !isfull(c) "Pushed data to a full accumulator."

    # Increase the step in our internal buffer
    # this should be guaranteed to always be < max_capacity
    c.buf_n = c.buf_n + 1

    buf_i = NeuralQuantum.unsafe_get_batch(c.in_buf, c.buf_n)
    state_copy!(buf_i, v)

    return nothing
end

function (c::AccumulatorLogGradPsi)(cngs)
    @assert !isfull(c) "Pushed data to a full accumulator."
    c.buf_n = c.buf_n + 1

    σp = unsafe_get_batch(c.in_buf, c.buf_n)
    @inbounds apply!(σp, cngs)
end

function (c::AccumulatorLogGradPsi)(cngs_l, cngs_r)
    @assert !isfull(c) "Pushed data to a full accumulator."
    c.buf_n = c.buf_n + 1

    σp = unsafe_get_batch(c.in_buf, c.buf_n)
    @inbounds apply!(σp, cngs_l, cngs_r)
end

"""
    process_accumulator!(c)

Processes all states stored in the accumulator, by computing their
relative local contribution.

It is safe to call this even if the accumulator is not full. In this
case all data beyond the count should be disregarded as it was
not initialized.

The output will be returned. You should not assume ownership of
the output, as it is preallocated and will be used for further
computations of the accumulator.
"""
function process_accumulator!(c::AccumulatorLogGradPsi)
    out  = c.out_buf
    ∇out = c.∇out_buf

    # Compute the batch of logψ neural networks
    logψ_and_∇logψ!(∇out, out, c.bnet, c.in_buf)

    # Reset the counter of the batch accumulator
    init!(c)

    return ∇out, out
end


text:
using NonconvexMMA
function convert_common_kwargs(opt::Union{NonconvexMMA.MMA02, NonconvexMMA.MMA87}, opt_kwargs;
    callback=nothing,
    maxiters=nothing,
    maxtime=nothing,
    abstol=nothing,
    reltol=nothing)

    conv_opt_kwargs = (; opt_kwargs...)

    if !isnothing(callback)
        @warn "common callback argument is currently not used by $(opt)"
    end

    if !isnothing(maxiters)
        conv_opt_kwargs = (; conv_opt_kwargs..., outer_maxiter=maxiters)
    end

    if !isnothing(maxtime)
        @warn "common maxtime argument is currently not used by $(opt)"
    end

    if !isnothing(abstol)
        tol_tmp = :tol .∈ Ref(keys(conv_opt_kwargs)) ? conv_opt_kwargs[:tol] : Nonconvex.NonconvexCore.Tolerance()
        tol_tmp_args = (; zip(propertynames(tol_tmp), [getproperty(tol_tmp, j) for j in propertynames(tol_tmp)])...)
        tol_tmp_args = (; tol_tmp_args..., fabs = abstol)

        conv_opt_kwargs = (; conv_opt_kwargs..., tol = Nonconvex.NonconvexCore.Tolerance(tol_tmp_args...))
    end

    if !isnothing(reltol)
        tol_tmp = :tol .∈ Ref(keys(conv_opt_kwargs)) ? conv_opt_kwargs[:tol] : Nonconvex.NonconvexCore.Tolerance()
        tol_tmp_args = (; zip(propertynames(tol_tmp), [getproperty(tol_tmp, j) for j in propertynames(tol_tmp)])...)
        tol_tmp_args = (; tol_tmp_args..., frel= reltol)

        conv_opt_kwargs = (; conv_opt_kwargs..., tol = Nonconvex.NonconvexCore.Tolerance(tol_tmp_args...))
    end

    return conv_opt_kwargs
end

function __create_options(opt::Union{NonconvexMMA.MMA02, NonconvexMMA.MMA87};
    opt_kwargs=nothing)

    options = !isnothing(opt_kwargs) ? NonconvexMMA.MMAOptions(;opt_kwargs...) : NonconvexMMA.MMAOptions()

    return options
end

function _create_options(opt::Union{NonconvexMMA.MMA02, NonconvexMMA.MMA87};
    opt_kwargs=nothing,
    sub_options=nothing,
    convergence_criteria=nothing)

    if !isnothing(convergence_criteria)
        options = (; options =  __create_options(opt, opt_kwargs=opt_kwargs), convcriteria=convergence_criteria)
    else
        options = (; options =  __create_options(opt, opt_kwargs=opt_kwargs))
    end

    return options
end

check_optimizer_backend(opt::Union{NonconvexMMA.MMA02, NonconvexMMA.MMA87}) = false


text:
"""
    MatrixBeta(p, n1, n2)
```julia
p::Int    dimension
n1::Real  degrees of freedom (greater than p - 1)
n2::Real  degrees of freedom (greater than p - 1)
```
The [matrix beta distribution](https://en.wikipedia.org/wiki/Matrix_variate_beta_distribution)
generalizes the beta distribution to ``p\\times p`` real matrices ``\\mathbf{U}``
for which ``\\mathbf{U}`` and ``\\mathbf{I}_p-\\mathbf{U}`` are both positive definite.
If ``\\mathbf{U}\\sim \\textrm{MB}_p(n_1/2, n_2/2)``, then its probability density function is

```math
f(\\mathbf{U}; n_1,n_2) = \\frac{\\Gamma_p(\\frac{n_1+n_2}{2})}{\\Gamma_p(\\frac{n_1}{2})\\Gamma_p(\\frac{n_2}{2})}
|\\mathbf{U}|^{(n_1-p-1)/2}\\left|\\mathbf{I}_p-\\mathbf{U}\\right|^{(n_2-p-1)/2}.
```

If ``\\mathbf{S}_1\\sim \\textrm{W}_p(n_1,\\mathbf{I}_p)`` and
``\\mathbf{S}_2\\sim \\textrm{W}_p(n_2,\\mathbf{I}_p)``
are independent, and we use ``\\mathcal{L}(\\cdot)`` to denote the lower Cholesky factor, then

```math
\\mathbf{U}=\\mathcal{L}(\\mathbf{S}_1+\\mathbf{S}_2)^{-1}\\mathbf{S}_1\\mathcal{L}(\\mathbf{S}_1+\\mathbf{S}_2)^{-\\rm{T}}
```

has ``\\mathbf{U}\\sim \\textrm{MB}_p(n_1/2, n_2/2)``.
"""
struct MatrixBeta{T <: Real, TW <: Wishart} <: ContinuousMatrixDistribution
    W1::TW
    W2::TW
    logc0::T
end

#  -----------------------------------------------------------------------------
#  Constructors
#  -----------------------------------------------------------------------------

function MatrixBeta(p::Int, n1::Real, n2::Real)
    p > 0 || throw(ArgumentError("dim must be positive: got $(p)."))
    logc0 = matrixbeta_logc0(p, n1, n2)
    T = Base.promote_eltype(n1, n2, logc0)
    Ip = ScalMat(p, one(T))
    W1 = Wishart(T(n1), Ip)
    W2 = Wishart(T(n2), Ip)
    MatrixBeta{T, typeof(W1)}(W1, W2, T(logc0))
end

#  -----------------------------------------------------------------------------
#  REPL display
#  -----------------------------------------------------------------------------

 show(io::IO, d::MatrixBeta) = show_multline(io, d, [(:n1, d.W1.df), (:n2, d.W2.df)])

#  -----------------------------------------------------------------------------
#  Conversion
#  -----------------------------------------------------------------------------

function convert(::Type{MatrixBeta{T}}, d::MatrixBeta) where T <: Real
    W1 = convert(Wishart{T}, d.W1)
    W2 = convert(Wishart{T}, d.W2)
    MatrixBeta{T, typeof(W1)}(W1, W2, T(d.logc0))
end

function convert(::Type{MatrixBeta{T}}, W1::Wishart, W2::Wishart, logc0) where T <: Real
    WW1 = convert(Wishart{T}, W1)
    WW2 = convert(Wishart{T}, W2)
    MatrixBeta{T, typeof(WW1)}(WW1, WW2, T(logc0))
end

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

dim(d::MatrixBeta) = dim(d.W1)

size(d::MatrixBeta) = size(d.W1)

rank(d::MatrixBeta) = dim(d)

insupport(d::MatrixBeta, U::AbstractMatrix) = isreal(U) && size(U) == size(d) && isposdef(U) && isposdef(I - U)

params(d::MatrixBeta) = (dim(d), d.W1.df, d.W2.df)

mean(d::MatrixBeta) = ((p, n1, n2) = params(d); Matrix((n1 / (n1 + n2)) * I, p, p))

@inline partype(d::MatrixBeta{T}) where {T <: Real} = T

#  Konno (1988 JJSS) Corollary 3.3.i
function cov(d::MatrixBeta, i::Integer, j::Integer, k::Integer, l::Integer)
    p, n1, n2 = params(d)
    n = n1 + n2
    Ω = Matrix{partype(d)}(I, p, p)
    n1*n2*inv(n*(n - 1)*(n + 2))*(-(2/n)*Ω[i,j]*Ω[k,l] + Ω[j,l]*Ω[i,k] + Ω[i,l]*Ω[k,j])
end

function var(d::MatrixBeta, i::Integer, j::Integer)
    p, n1, n2 = params(d)
    n = n1 + n2
    Ω = Matrix{partype(d)}(I, p, p)
    n1*n2*inv(n*(n - 1)*(n + 2))*((1 - (2/n))*Ω[i,j]^2 + Ω[j,j]*Ω[i,i])
end

#  -----------------------------------------------------------------------------
#  Evaluation
#  -----------------------------------------------------------------------------

function matrixbeta_logc0(p::Int, n1::Real, n2::Real)
    #  returns the natural log of the normalizing constant for the pdf
    return -logmvbeta(p, n1 / 2, n2 / 2)
end

function logkernel(d::MatrixBeta, U::AbstractMatrix)
    p, n1, n2 = params(d)
    ((n1 - p - 1) / 2) * logdet(U) + ((n2 - p - 1) / 2) * logdet(I - U)
end

#  -----------------------------------------------------------------------------
#  Sampling
#  -----------------------------------------------------------------------------

#  Mitra (1970 Sankhyā)

function _rand!(rng::AbstractRNG, d::MatrixBeta, A::AbstractMatrix)
    S1   = PDMat( rand(rng, d.W1) )
    S2   = PDMat( rand(rng, d.W2) )
    S    = S1 + S2
    invL = Matrix( inv(S.chol.L) )
    A .= X_A_Xt(S1, invL)
end

#  -----------------------------------------------------------------------------
#  Test utils
#  -----------------------------------------------------------------------------

function _univariate(d::MatrixBeta)
    check_univariate(d)
    p, n1, n2 = params(d)
    return Beta(n1 / 2, n2 / 2)
end

function _rand_params(::Type{MatrixBeta}, elty, n::Int, p::Int)
    n == p || throw(ArgumentError("dims must be equal for MatrixBeta"))
    n1 = elty( n + 1 + abs(10randn()) )
    n2 = elty( n + 1 + abs(10randn()) )
    return n, n1, n2
end


text:
using MDCDL
using Images, TestImages
using ImageFiltering
using Plots
using LinearAlgebra

function main()
    ### configurations ###
    # dirNsolt = joinpath(Pkg.dir(),"MDCDL","examples","design","sample.json")
    dirNsolt = joinpath(pathof(MDCDL), "..","..","examples", "design", "sample.json")
    dirNsolt = abspath(dirNsolt)
    σ_noise = 1e-2 # variance of AWGN
    λ = 1e-2 # FISTA parameter
    lv = 3 # tree level of NSOLT
    ######################

    ### setup observed image
    # orgImg = testimage("cameraman")
    # u = Array{Float64}(imresize(orgImg,(256,256)))
    orgImg = testimage("lena")
    u = Array{RGB{Float64}}(orgImg)

    # AWGN
    w = mapc.((v)->σ_noise * randn(), u)

    # Degradation (pixel loss)
    P = map((v)-> v > 0.01 ? 1 : 0, rand(Float64,size(u)))

    # Observation
    x = P .* u + w

    ### setup NSOLT
    nsolt = MDCDL.load(dirNsolt)
    println("NSOLT configurations:")
    println(" - Type = $(typeof(nsolt))")
    println(" - Decimation Factor = $(nsolt.decimationFactor)")
    println(" - Number of Channels = $(nsolt.nChannels)")
    println(" - Polyphase Order = $(nsolt.polyphaseOrder)")
    ### show analysis filters

    # setup Multiscale NSOLT
    analyzer = createMultiscaleAnalyzer(nsolt, size(x); level=lv, shape=:vector)
    synthesizer = createMultiscaleSynthesizer(nsolt, size(x); level=lv, shape=:vector)

    # setup FISTA
    y0 = analyze(analyzer, x)

    ∇f(ty) = begin
        - analyze(analyzer, P.*(x - synthesize(synthesizer, ty)))
    end
    prox(ty, η) = MDCDL.softshrink(ty, λ*η)
    viewFcn(itrs, ty, err) = begin
        println("# Iterations=$itrs, error=$err")
    end

    hy = MDCDL.fista(∇f, prox, y0; η=1.0, iterations=400, verboseFunction=viewFcn, absTol=eps())

    # restored image
    ru = synthesize(synthesizer, hy)

    psnr = (a, ref) -> 10*log10(1/(norm(a-ref)^2/length(a)))
    errx = norm(ru - u)

    println("error: $errx")

    plotOrg = plot(u ; xlabel=string("Original Image"), aspect_ratio=:equal);
    plotRes = plot(ru; xlabel=string("Restored Image\nPSNR=",trunc(psnr(ru,u),digits=3)), aspect_ratio=:equal);
    plot(plotOrg, plotObs, plotRes; aspect_ratio=:equal, layout=(1,3))

end

main()

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
abstract type AbstractAggregate end

struct Reduce{O, C, A} <: AbstractAggregate
    op::O
    condf::C
    adjust::A
    checkempty::Bool
end
Reduce(f, condf=nothing, adjust=nothing) = Reduce(f, condf, adjust, false)

check_aggregate(f::Any, ::AbstractVector) = f
check_aggregate(f::typeof(sum), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.add_sum)
check_aggregate(f::typeof(sum∘skipmissing), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.add_sum, !ismissing)
check_aggregate(f::typeof(prod), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.mul_prod)
check_aggregate(f::typeof(prod∘skipmissing), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.mul_prod, !ismissing)
check_aggregate(f::typeof(maximum),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(maximum), v::AbstractVector{<:Union{Missing, Real}}) =
    eltype(v) === Any ? f : Reduce(max)
check_aggregate(f::typeof(maximum∘skipmissing),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(maximum∘skipmissing), v::AbstractVector{<:Union{Missing, Real}}) =
    eltype(v) === Any ? f : Reduce(max, !ismissing, nothing, true)
check_aggregate(f::typeof(minimum),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(minimum), v::AbstractVector{<:Union{Missing, Real}}) =
    eltype(v) === Any ? f : Reduce(min)
check_aggregate(f::typeof(minimum∘skipmissing),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(minimum∘skipmissing), v::AbstractVector{<:Union{Missing, Real}}) =
    eltype(v) === Any ? f : Reduce(min, !ismissing, nothing, true)
check_aggregate(f::typeof(mean), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.add_sum, nothing, /)
check_aggregate(f::typeof(mean∘skipmissing), ::AbstractVector{<:Union{Missing, Number}}) =
    Reduce(Base.add_sum, !ismissing, /)

# Other aggregate functions which are not strictly reductions
struct Aggregate{F, C} <: AbstractAggregate
    f::F
    condf::C
end
Aggregate(f) = Aggregate(f, nothing)

check_aggregate(f::typeof(var), ::AbstractVector{<:Union{Missing, Number}}) =
    Aggregate(var)
check_aggregate(f::typeof(var∘skipmissing), ::AbstractVector{<:Union{Missing, Number}}) =
    Aggregate(var, !ismissing)
check_aggregate(f::typeof(std), ::AbstractVector{<:Union{Missing, Number}}) =
    Aggregate(std)
check_aggregate(f::typeof(std∘skipmissing), ::AbstractVector{<:Union{Missing, Number}}) =
    Aggregate(std, !ismissing)
check_aggregate(f::typeof(first), v::AbstractVector) =
    eltype(v) === Any ? f : Aggregate(first)
check_aggregate(f::typeof(first),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(first∘skipmissing), v::AbstractVector) =
    eltype(v) === Any ? f : Aggregate(first, !ismissing)
check_aggregate(f::typeof(first∘skipmissing),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(last), v::AbstractVector) =
    eltype(v) === Any ? f : Aggregate(last)
check_aggregate(f::typeof(last),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(last∘skipmissing), v::AbstractVector) =
    eltype(v) === Any ? f : Aggregate(last, !ismissing)
check_aggregate(f::typeof(last∘skipmissing),
                ::AbstractVector{<:Union{Missing, MULTI_COLS_TYPE, AbstractVector}}) = f
check_aggregate(f::typeof(length), ::AbstractVector) = Aggregate(length)

# SkipMissing does not support length

# Use a strategy similar to reducedim_init from Base to get the vector of the right type
function groupreduce_init(op, condf, adjust,
                          incol::AbstractVector{U}, gd::GroupedDataFrame) where U
    T = Base.promote_union(U)

    if op === Base.add_sum
        initf = zero
    elseif op === Base.mul_prod
        initf = one
    else
        throw(ErrorException("Unrecognized op $op"))
    end

    Tnm = nonmissingtype(T)
    if isconcretetype(Tnm) && applicable(initf, Tnm)
        tmpv = initf(Tnm)
        initv = op(tmpv, tmpv)
        if adjust isa Nothing
            x = Tnm <: AbstractIrrational ? float(initv) : initv
        else
            x = adjust(initv, 1)
        end
        if condf === !ismissing
            V = typeof(x)
        else
            V = U >: Missing ? Union{typeof(x), Missing} : typeof(x)
        end
        v = similar(incol, V, length(gd))
        fill!(v, x)
        return v
    else
        # do not try to determine the narrowest possible type nor starting value
        # as this is not possible to do correctly in general without processing
        # groups; it will get fixed later in groupreduce!; later we
        # will make use of the fact that this vector is filled with #undef
        # while above the vector is filled with a concrete value
        return Vector{Any}(undef, length(gd))
    end
end

for (op, initf) in ((:max, :typemin), (:min, :typemax))
    @eval begin
        function groupreduce_init(::typeof($op), condf, adjust,
                                  incol::AbstractVector{T}, gd::GroupedDataFrame) where T
            @assert isnothing(adjust)
            S = nonmissingtype(T)
            # !ismissing check is purely an optimization to avoid a copy later
            outcol = similar(incol, condf === !ismissing ? S : T, length(gd))
            # Comparison is possible only between CatValues from the same pool
            outcolT = typeof(outcol).name
            if outcolT.name === :CategoricalArray &&
                nameof(outcolT.module) === :CategoricalArrays
                # we know that CategoricalArray has `pool` field
                outcol.pool = incol.pool
            end
            # It is safe to use a non-missing init value
            # since missing will poison the result if present
            # we assume here that groups are non-empty (current design assures this)
            # + workaround for https://github.com/JuliaLang/julia/issues/36978
            if isconcretetype(S) && hasmethod($initf, Tuple{S}) && !(S <: Irrational)
                fill!(outcol, $initf(S))
            else
                fillfirst!(condf, outcol, incol, gd)
            end
            return outcol
        end
    end
end

function copyto_widen!(res::AbstractVector{T}, x::AbstractVector) where T
    @inbounds for i in eachindex(res, x)
        val = x[i]
        S = typeof(val)
        if S <: T || promote_type(S, T) <: T
            res[i] = val
        else
            newres = Tables.allocatecolumn(promote_type(S, T), length(x))
            return copyto_widen!(newres, x)
        end
    end
    return res
end

function groupreduce!(res::AbstractVector, f, op, condf, adjust, checkempty::Bool,
                      incol::AbstractVector, gd::GroupedDataFrame)
    n = length(gd)
    if adjust !== nothing || checkempty
        counts = zeros(Int, n)
    end
    groups = gd.groups
    @inbounds for i in eachindex(incol, groups)
        gix = groups[i]
        x = incol[i]
        if gix > 0 && (condf === nothing || condf(x))
            # this check should be optimized out if U is not Any
            if eltype(res) === Any && !isassigned(res, gix)
                res[gix] = f(x, gix)
            else
                res[gix] = op(res[gix], f(x, gix))
            end
            if adjust !== nothing || checkempty
                counts[gix] += 1
            end
        end
    end
    # handle the case of an unitialized reduction
    if eltype(res) === Any
        if op === Base.add_sum
            initf = zero
        elseif op === Base.mul_prod
            initf = one
        else
            initf = x -> throw(ErrorException("Unrecognized op $op"))
        end
        @inbounds for gix in eachindex(res)
            if !isassigned(res, gix)
                res[gix] = initf(nonmissingtype(eltype(incol)))
            end
        end
    end
    if adjust !== nothing
        res .= adjust.(res, counts)
    end
    if checkempty && any(iszero, counts)
        throw(ArgumentError("some groups contain only missing values"))
    end
    # Reallocate Vector created in groupreduce_init with min or max
    # for CategoricalVector
    resT = typeof(res).name
    if resT.name === :CategoricalArray &&
        nameof(resT.module) === :CategoricalArrays
        @assert op === min || op === max
        # we know that CategoricalArray has `pool` field
        @assert res.pool === incol.pool
        res.pool = copy(incol.pool)
    end
    if isconcretetype(eltype(res))
        return res
    else
        return copyto_widen!(Tables.allocatecolumn(typeof(first(res)), n), res)
    end
end

# function barrier works around type instability of groupreduce_init due to applicable
groupreduce(f, op, condf, adjust, checkempty::Bool,
            incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce!(groupreduce_init(op, condf, adjust, incol, gd),
                 f, op, condf, adjust, checkempty, incol, gd)
# Avoids the overhead due to Missing when computing reduction
groupreduce(f, op, condf::typeof(!ismissing), adjust, checkempty::Bool,
            incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce!(disallowmissing(groupreduce_init(op, condf, adjust, incol, gd)),
                 f, op, condf, adjust, checkempty, incol, gd)

(r::Reduce)(incol::AbstractVector, gd::GroupedDataFrame) =
    groupreduce((x, i) -> x, r.op, r.condf, r.adjust, r.checkempty, incol, gd)

# this definition is missing in Julia 1.0 LTS and is required by aggregation for var
# TODO: remove this when we drop 1.0 support
if VERSION < v"1.1"
    Base.zero(::Type{Missing}) = missing
end

function (agg::Aggregate{typeof(var)})(incol::AbstractVector, gd::GroupedDataFrame)
    means = groupreduce((x, i) -> x, Base.add_sum, agg.condf, /, false, incol, gd)
    # !ismissing check is purely an optimization to avoid a copy later
    if eltype(means) >: Missing && agg.condf !== !ismissing
        T = Union{Missing, real(eltype(means))}
    else
        T = real(eltype(means))
    end
    res = zeros(T, length(gd))
    return groupreduce!(res, (x, i) -> @inbounds(abs2(x - means[i])), +, agg.condf,
                        (x, l) -> l <= 1 ? oftype(x / (l-1), NaN) : x / (l-1),
                        false, incol, gd)
end

function (agg::Aggregate{typeof(std)})(incol::AbstractVector, gd::GroupedDataFrame)
    outcol = Aggregate(var, agg.condf)(incol, gd)
    if eltype(outcol) <: Union{Missing, Rational}
        return sqrt.(outcol)
    else
        return map!(sqrt, outcol, outcol)
    end
end

for f in (first, last)
    function (agg::Aggregate{typeof(f)})(incol::AbstractVector, gd::GroupedDataFrame)
        n = length(gd)
        outcol = similar(incol, n)
        fillfirst!(agg.condf, outcol, incol, gd, rev=agg.f === last)
        if isconcretetype(eltype(outcol))
            return outcol
        else
            return copyto_widen!(Tables.allocatecolumn(typeof(first(outcol)), n), outcol)
        end
    end
end

function (agg::Aggregate{typeof(length)})(incol::AbstractVector, gd::GroupedDataFrame)
    if getfield(gd, :idx) === nothing
        lens = zeros(Int, length(gd))
        @inbounds for gix in gd.groups
            gix > 0 && (lens[gix] += 1)
        end
        return lens
    else
        return gd.ends .- gd.starts .+ 1
    end
end

isagg((col, (fun, outcol))::Pair{<:ColumnIndex, <:Pair{<:Any, <:SymbolOrString}}, gdf::GroupedDataFrame) =
    check_aggregate(fun, parent(gdf)[!, col]) isa AbstractAggregate
isagg(::Any, gdf::GroupedDataFrame) = false


text:
# Compute order-m efficiency score
function orderm{T<:AbstractDEA}(A::T,M::Int,B::Int)
  Data = getdata(A)
  K = getnrdmu(Data)
  # getdata() returns a shallow copy of the data in A, but we need a deep copy!
  Datacopy = deepcopy(Data)
  if(M > K)
    error("M cannot be larger than the number of DMUs!")
  end
  theta = Array(Float64,K,B)
  @sync @parallel for i=1:B
    setindexes!(Data,rand(1:K,M))
    for j in eachindex(Datacopy)
      theta[j,i] = geteff(A(Datacopy[j]...))
    end
  end
  setindexes!(Data,collect(1:K))
  return mean(theta,2)
end

# Compute order-m efficiency score for FDH
function orderm{T<:RS}(A::FDH{T},M::Int,B::Int)
  Data = getdata(A)
  K = getnrdmu(Data)
  # getdata() returns a shallow copy of the data in A, but we need a deep copy!
  Datacopy = deepcopy(Data)
  if(M > K)
    error("M cannot be larger than the number of DMUs!")
  end
  X,Y = Data[1:end]
  theta = Array(Float64,K,B)
  for j in eachindex(Datacopy)
    # Condition on Y >= y for input-oriented and X <= x for output-oriented
    Xk,Yk = Datacopy[j]
    if(A.input)
      dom = find(Y .>= Yk)
    else
      dom = find(X .<= Xk)
    end
    @sync @parallel for i=1:B
      setindexes!(Data,rand(1:K,M))
      theta[j,i] = geteff(A(Xk,Yk))
    end
  end
  setindexes!(Data,collect(1:K))
  return mean(theta,2)
end

max_stars_repo_name: kepiej/NEPA


text:
using GLMakie
using MixedModels
using MixedModels: dataset
using MixedModelsMakie

# we scale weirdly here so that we have nice scaling in the 3d plot
m = fit(MixedModel,
        @formula(rt_trunc / 100 ~ 1 + spkr*prec*load + (1|subj) + (1+spkr+prec|item)),
        dataset(:kb07);
        contrasts=Dict(:load => HelmertCoding(),
                        :spkr => HelmertCoding(),
                        :prec => HelmertCoding()))

shrinkageplot(m, :item)

let gf = :item, θref = 10000*m.optsum.initial
    reind = findfirst(==(gf), fnames(m))  # convert the symbol gf to an index
    if isnothing(reind)
        throw(ArgumentError("gf=$gf is not one of the grouping factor names, $(fnames(m))"))
    end
    reest = ranef(m)[reind]               # random effects conditional means at estimated θ
    reref = ranef(updateL!(setθ!(m, θref)))[reind]  # same at θref
    updateL!(setθ!(m, m.optsum.final))    # restore parameter estimates and update m

    tails = [Point3f(reref[:, j]) for j in 1:size(reref,2)]
    lengths = [Vec3f(reest[:, j] .- reref[:, j]) for j in 1:size(reref,2)]
    f = Figure()
    ax3 = Axis3(f[1,1])
    arrows!(ax3, tails, lengths)
    ax3.xlabel = "Intercept"
    ax3.ylabel = "spkr"
    ax3.zlabel = "prec"
    f
end


text:
# testing subviews

using ArrayViews
using Base.Test

#### Test Subviews

## tools to facilitate array view testing
function _test_arrview_contents(v, r)
    siz_r = size(r)
    siz_v = size(v)

    if siz_r != siz_v
        error("Incorrect size: get $(siz_v), but expect $(siz_r)")
    end

    for i = 1 : length(v)
        if v[i] != r[i]
            println("v = ")
            println(v)
            println("r = ")
            println(r)
            error("Incorrect content.")
        end
    end
end

function _test_arrview(a, r, subs...)
    _test_arrview_contents(view(a, subs...), r)
    _test_arrview_contents(unsafe_view(a, subs...), r)
end

if VERSION >= v"0.4"
    macro test_arrview(a_, subs...)
        esc(:(_test_arrview($a_, sub($a_, $(subs...)), $(subs...))))
    end
else
    macro test_arrview(a_, subs...)
        esc(:(_test_arrview($a_, $a_[$(subs...)], $(subs...))))
    end
end

#### test views from arrays

avparent = copy(reshape(1.:1680., (8, 7, 6, 5)))

# 1D
@test_arrview(avparent, 3)
@test_arrview(avparent, :)
@test_arrview(avparent, 1:12)
@test_arrview(avparent, 3:2:36)

# 2D
@test_arrview(avparent, 4, 2)
@test_arrview(avparent, 4, :)
@test_arrview(avparent, 4, 3:10)
@test_arrview(avparent, 4, 2:2:10)

@test_arrview(avparent, :, 2)
@test_arrview(avparent, :, :)
@test_arrview(avparent, :, 3:10)
@test_arrview(avparent, :, 2:2:10)

@test_arrview(avparent, 1:6, 2)
@test_arrview(avparent, 1:6, :)
@test_arrview(avparent, 1:6, 3:10)
@test_arrview(avparent, 1:6, 2:2:10)

@test_arrview(avparent, 1:2:8, 2)
@test_arrview(avparent, 1:2:8, :)
@test_arrview(avparent, 1:2:8, 3:10)
@test_arrview(avparent, 1:2:8, 2:2:10)

# 3D
@test_arrview(avparent, 4, 3, 2)
@test_arrview(avparent, 4, 3, :)
@test_arrview(avparent, 4, 3, 2:5)
@test_arrview(avparent, 4, 3, 1:2:6)

@test_arrview(avparent, 4, :, 2)
@test_arrview(avparent, 4, :, :)
@test_arrview(avparent, 4, :, 2:5)
@test_arrview(avparent, 4, :, 1:2:6)

@test_arrview(avparent, 4, 3:7, 2)
@test_arrview(avparent, 4, 3:7, :)
@test_arrview(avparent, 4, 3:7, 2:5)
@test_arrview(avparent, 4, 3:7, 1:2:6)

@test_arrview(avparent, 4, 1:2:5, 2)
@test_arrview(avparent, 4, 1:2:5, :)
@test_arrview(avparent, 4, 1:2:5, 2:5)
@test_arrview(avparent, 4, 1:2:5, 1:2:6)

@test_arrview(avparent, :, 3, 2)
@test_arrview(avparent, :, 3, :)
@test_arrview(avparent, :, 3, 2:5)
@test_arrview(avparent, :, 3, 1:2:6)

@test_arrview(avparent, :, :, 2)
@test_arrview(avparent, :, :, :)
@test_arrview(avparent, :, :, 2:5)
@test_arrview(avparent, :, :, 1:2:6)

@test_arrview(avparent, :, 3:7, 2)
@test_arrview(avparent, :, 3:7, :)
@test_arrview(avparent, :, 3:7, 2:5)
@test_arrview(avparent, :, 3:7, 1:2:6)

@test_arrview(avparent, :, 1:2:5, 2)
@test_arrview(avparent, :, 1:2:5, :)
@test_arrview(avparent, :, 1:2:5, 2:5)
@test_arrview(avparent, :, 1:2:5, 1:2:6)

@test_arrview(avparent, 2:7, 3, 2)
@test_arrview(avparent, 2:7, 3, :)
@test_arrview(avparent, 2:7, 3, 2:5)
@test_arrview(avparent, 2:7, 3, 1:2:6)

@test_arrview(avparent, 2:7, :, 2)
@test_arrview(avparent, 2:7, :, :)
@test_arrview(avparent, 2:7, :, 2:5)
@test_arrview(avparent, 2:7, :, 1:2:6)

@test_arrview(avparent, 2:7, 3:7, 2)
@test_arrview(avparent, 2:7, 3:7, :)
@test_arrview(avparent, 2:7, 3:7, 2:5)
@test_arrview(avparent, 2:7, 3:7, 1:2:6)

@test_arrview(avparent, 2:7, 1:2:5, 2)
@test_arrview(avparent, 2:7, 1:2:5, :)
@test_arrview(avparent, 2:7, 1:2:5, 2:5)
@test_arrview(avparent, 2:7, 1:2:5, 1:2:6)

@test_arrview(avparent, 1:2:7, 3, 2)
@test_arrview(avparent, 1:2:7, 3, :)
@test_arrview(avparent, 1:2:7, 3, 2:5)
@test_arrview(avparent, 1:2:7, 3, 1:2:6)

@test_arrview(avparent, 1:2:7, :, 2)
@test_arrview(avparent, 1:2:7, :, :)
@test_arrview(avparent, 1:2:7, :, 2:5)
@test_arrview(avparent, 1:2:7, :, 1:2:6)

@test_arrview(avparent, 1:2:7, 3:7, 2)
@test_arrview(avparent, 1:2:7, 3:7, :)
@test_arrview(avparent, 1:2:7, 3:7, 2:5)
@test_arrview(avparent, 1:2:7, 3:7, 1:2:6)

@test_arrview(avparent, 1:2:7, 1:2:5, 2)
@test_arrview(avparent, 1:2:7, 1:2:5, :)
@test_arrview(avparent, 1:2:7, 1:2:5, 2:5)
@test_arrview(avparent, 1:2:7, 1:2:5, 1:2:6)

# Some 4D Tests
@test_arrview(avparent, 4, :,     3, 4)
@test_arrview(avparent, 4, :,     :, 4)
@test_arrview(avparent, 4, :,   3:5, 4)
@test_arrview(avparent, 4, :, 1:2:5, 4)

@test_arrview(avparent, :, :,     3, 4)
@test_arrview(avparent, :, :,     :, 4)
@test_arrview(avparent, :, :,   3:5, 4)
@test_arrview(avparent, :, :, 1:2:5, 4)

@test_arrview(avparent, 2:7, :,     3, 4)
@test_arrview(avparent, 2:7, :,     :, 4)
@test_arrview(avparent, 2:7, :,   3:5, 4)
@test_arrview(avparent, 2:7, :, 1:2:5, 4)

@test_arrview(avparent, 1:2:7, :,     3, 4)
@test_arrview(avparent, 1:2:7, :,     :, 4)
@test_arrview(avparent, 1:2:7, :,   3:5, 4)
@test_arrview(avparent, 1:2:7, :, 1:2:5, 4)

@test_arrview(avparent, 4, :,     3, :)
@test_arrview(avparent, 4, :,     :, :)
@test_arrview(avparent, 4, :,   3:5, :)
@test_arrview(avparent, 4, :, 1:2:5, :)

@test_arrview(avparent, :, :,     3, :)
@test_arrview(avparent, :, :,     :, :)
@test_arrview(avparent, :, :,   3:5, :)
@test_arrview(avparent, :, :, 1:2:5, :)

@test_arrview(avparent, 2:7, :,     3, :)
@test_arrview(avparent, 2:7, :,     :, :)
@test_arrview(avparent, 2:7, :,   3:5, :)
@test_arrview(avparent, 2:7, :, 1:2:5, :)

@test_arrview(avparent, 1:2:7, :,     3, :)
@test_arrview(avparent, 1:2:7, :,     :, :)
@test_arrview(avparent, 1:2:7, :,   3:5, :)
@test_arrview(avparent, 1:2:7, :, 1:2:5, :)

@test_arrview(avparent, 4, :,     3, 2:5)
@test_arrview(avparent, 4, :,     :, 2:5)
@test_arrview(avparent, 4, :,   3:5, 2:5)
@test_arrview(avparent, 4, :, 1:2:5, 2:5)

@test_arrview(avparent, :, :,     3, 2:5)
@test_arrview(avparent, :, :,     :, 2:5)
@test_arrview(avparent, :, :,   3:5, 2:5)
@test_arrview(avparent, :, :, 1:2:5, 2:5)

@test_arrview(avparent, 2:7, :,     3, 2:5)
@test_arrview(avparent, 2:7, :,     :, 2:5)
@test_arrview(avparent, 2:7, :,   3:5, 2:5)
@test_arrview(avparent, 2:7, :, 1:2:5, 2:5)

@test_arrview(avparent, 1:2:7, :,     3, 2:5)
@test_arrview(avparent, 1:2:7, :,     :, 2:5)
@test_arrview(avparent, 1:2:7, :,   3:5, 2:5)
@test_arrview(avparent, 1:2:7, :, 1:2:5, 2:5)

@test_arrview(avparent, :, 1, 1, 1)
@test_arrview(avparent, :, 1, 2, 1)
@test_arrview(avparent, :, :, 2, 2)
@test_arrview(avparent, :, :, 1:2, 3:4)

# Some 5D tests
avparent2 = copy(reshape(1.:6720., (8, 7, 6, 5, 4)))

@test_arrview(avparent2, :, 1, 1, 1, 1)
@test_arrview(avparent2, :, 1, 1, 3, 1)
@test_arrview(avparent2, :, 1, 2, 3, 3:4)
@test_arrview(avparent2, 2:7, 1, 1, 1, 3:4)

@test_arrview(avparent2, 1, :, 1, 3, 4)
@test_arrview(avparent2, 2, 3, :, 3, 4)
@test_arrview(avparent2, 3, 4, 5, :, 4)
@test_arrview(avparent2, 2, 3, 4, 4, :)

@test_arrview(avparent2, 1, 2:3, 3, 4, 2)
@test_arrview(avparent2, 2, 1, 3:4, 1, 1)
@test_arrview(avparent2, 3, 1, 1, 4:5, 2)
@test_arrview(avparent2, 4, 1, 1, 3:4, 2:4)

#### Test Subviews of Views

function print_subscripts(subs1, subs2)
    println("Error happens on: ")

    subs1s = join([repr(i) for i in subs1], ", ")
    subs2s = join([repr(i) for i in subs2], ", ")

    println("  subs1 = [$subs1s]")
    println("  subs2 = [$subs2s]")
end

function test_arrview2(a, subs1, subs2)
    v = view(a, subs1...)
    v2r = view(copy(v), subs2...)
    _test_arrview_contents(view(v, subs2...), v2r)
    uv = unsafe_view(a, subs1...)
    _test_arrview_contents(unsafe_view(uv, subs2...), v2r)
end

avparent = copy(reshape(1:6912, (12, 12, 8, 6)))

# 1D --> 1D
for sa in Any[Colon(), 1:36, 2:2:36]
    v1 = view(avparent, sa)
    for sb in Any[4, Colon(), 1:length(v1), 3:2:length(v1)]
        test_arrview2(avparent, (sa,), (sb,))
    end
end

# 2D --> 2D
for sa1 in Any[Colon(), 1:10, 2:2:12], sa2 = Any[Colon(), 1:12, 2:2:16]
    v1 = view(avparent, sa1, sa2)
    for sb1 in Any[4, Colon(), 2:size(v1,1), 2:2:size(v1,1)],
        sb2 in Any[4, Colon(), 2:size(v1,2), 2:2:size(v1,2)]
        test_arrview2(avparent, (sa1, sa2), (sb1, sb2))
    end
end

# 3D --> 3D
for sa1 in Any[Colon(), 1:10, 2:2:12],
    sa2 in Any[Colon(), 1:10, 2:2:12],
    sa3 in Any[Colon(), 1:7, 2:2:8]
    v1 = view(avparent, sa1, sa2, sa3)
    (d1, d2, d3) = size(v1)
    for sb1 in Any[Colon(), 2:d1, 2:2:d1],
        sb2 in Any[Colon(), 2:d2, 2:2:d2],
        sb3 in Any[Colon(), 2:d3, 2:2:d3]
        test_arrview2(avparent, (sa1, sa2, sa3), (sb1, sb2, sb3))
    end
end


text:
################################################################################
#
#   Roots
#
################################################################################

function roots(f::Generic.Poly{T}) where T <: Union{padic, qadic, LocalFieldElem}
  K = base_ring(f)
  e = absolute_ramification_index(K)
  k, mk = ResidueField(K)
  fk = map_coefficients(mk, f)
  #TODO: We don't need a full Hensel factorization.
  lH = Hensel_factorization(f)
  rt = elem_type(K)[]
  Kx = parent(f)
  for (phi, g) in lH
    if isone(degree(phi))
      if isone(degree(g))
        push!(rt, -constant_coefficient(g)//leading_coefficient(g))
      else
        #TODO: We don't need a full slope factorization.
        lS = slope_factorization(g)
        for (h, mh) in lS
          if isone(degree(h))
            r = -constant_coefficient(h)//leading_coefficient(h)
            for j = 1:mh
              push!(rt, r)
            end
          elseif iszero(constant_coefficient(h)) || divides(numerator(e*valuation(constant_coefficient(h))), degree(h))[1]
            rth = _roots(h)
            for j = 1:mh
              append!(rt, rth)
            end
          end
        end
      end
    end
  end
  #the roots need to be refined.
  #rt = refine_roots(f, rt)
  return rt
end

function refine_roots(f::Generic.Poly{T}, rt::Vector{T}) where T <: Union{padic, qadic, LocalFieldElem}
  Rx = parent(f)
  x = gen(Rx)
  factors = typeof(f)[x-y for y in rt]
  push!(factors, divexact(f, prod(factors)))
  factors = lift_factorization(f, factors)
  res = Vector{T}(undef, length(rt))
  for i = 1:length(rt)
    res[i] = -constant_coefficient(factors[i])
  end
  return res
end

function refine_roots1(f::Generic.Poly{T}, rt::Vector{T}) where T <: Union{padic, qadic, LocalFieldElem}
  K = base_ring(f)
  v = numerator(absolute_ramification_index(K)*valuation(reduced_discriminant(f)))
  target_prec = precision(f)
  starting = minimum(Int[precision(x) for x in rt])
  chain = [target_prec]
  i = target_prec
  while i > starting
    i = div(i+1, 2)
    pushfirst!(chain, i)
  end
  der = derivative(f)
  wvect = [inv(der(rt[i])) for i = 1:length(rt)]
  rtnew = copy(rt)
  for i in 1:length(chain)
    for j = 1:length(rtnew)
      wvect[j]*f(rtnew[j])
      rtnew[j] = rtnew[j] - wvect[j]*f(rtnew[j])
      wvect[j] = wvect[j]*(2-wvect[j]*der(rtnew[j]))
    end
  end
  return rtnew
end

function _roots(f::Generic.Poly{T}) where T <: Union{padic, qadic, LocalFieldElem}
  K = base_ring(f)
  k, mk = ResidueField(K)
  fk = map_coefficients(mk, f)
  rts = roots(fk)
  x = gen(parent(f))
  #TODO: Is this setprecision call ok?
  r = setprecision(lift(rts[1], K), precision(f))
  pi = uniformizer(K)
  g = f(pi*x+r)
  g = divexact(g, _content(g))
  rtg = roots(g)
  rts = elem_type(K)[setprecision(r, precision(y)) + pi*y for y in rtg]
  return rts
end

function automorphisms(K::T) where T <: Union{LocalField, FlintQadicField}
  rt = roots(defining_polynomial(K), K)
  return morphism_type(K)[hom(K, K, x) for x in rt]
end

function automorphisms(K::LocalField, L::T) where T <: Union{LocalField, FlintQadicField, FlintPadicField}
  return _automorphisms(K, K, L)
end

function absolute_automorphisms(K::LocalField{qadic, S}) where S
  autsk = small_generating_set(automorphisms(base_field(K)))
  auts = morphism_type(K)[]
  for f in autsk
    fnew = map_coefficients(f, defining_polynomial(K))
    rt = roots(fnew, K)
    for x in rt
      push!(auts, hom(K, K, f, x))
    end
  end
  return closure(auts, *)
end

function absolute_automorphisms(K::LocalField)
  return _automorphisms(K, K, absolute_base_field(K))
end

function absolute_automorphisms(K::FlintQadicField)
  return automorphisms(K)
end

function _automorphisms(K::S, F::T, L::U) where {S <: Union{LocalField, FlintQadicField, FlintPadicField}, T <: Union{LocalField, FlintQadicField, FlintPadicField}, U <: Union{LocalField, FlintQadicField, FlintPadicField}}
  if absolute_degree(K) < absolute_degree(L)
    error("The base field is not naturally a subfield!")
  end
  if K == L
    return morphism_type(K, F)[hom(K, F, F(gen(K)))]
  end
  autsk = _automorphisms(base_field(K), F, L)
  auts = morphism_type(K, F)[]
  for f in autsk
    rt = roots(map_coefficients(f, defining_polynomial(K)))
    for x in rt
      push!(auts, hom(K, F, f, x))
    end
  end
  return auts
end

function small_generating_set(auts::Vector{T}) where T <: LocalFieldMor
  @assert length(auts) >= 1
  @assert domain(auts[1]) == codomain(auts[1])
  for i = 1:length(auts)
    @assert domain(auts[i]) == codomain(auts[i])
    @assert domain(auts[1]) == domain(auts[i])
  end
  if length(auts) == 1
    return eltype(auts)[x for x in auts]
  end
  return small_generating_set(auts, *, id_hom(domain(auts[1])))
end

################################################################################
#
#   Automorphism group
#
################################################################################

function automorphism_group(K::LocalField)
  aut = automorphisms(K)
  mult_table = Matrix{Int}(undef, length(aut), length(aut))
  for s = 1:length(aut)
    for i = 1:length(aut)
      mult_table[s, i] = findfirst(isequal(aut[s]*aut[i]), aut)
    end
  end
  G = GrpGen(mult_table)
  return G, GrpGenToNfMorSet(G, aut, K)
end

@doc Markdown.doc"""
    automorphism_group(L::NumField, K::NumField) -> GenGrp, GrpGenToNfMorSet

Given the number field extension $L$ and $K$, this function returns a group $G$
and a map from $G$ to the automorphisms of $L$ that fix $K$.
"""
function automorphism_group(L::LocalField, K::LocalField)
  aut = automorphisms(L, K)
  mult_table = Matrix{Int}(undef, length(aut), length(aut))
  for s = 1:length(aut)
    for i = 1:length(aut)
      mult_table[s, i] = findfirst(isequal(aut[s]*aut[i]), aut)
    end
  end
  G = GrpGen(mult_table)
  return G, GrpGenToNfMorSet(G, aut, L)
end

@doc Markdown.doc"""
    absolute_automorphism_group(L::LocalField) -> GenGrp, GrpGenToNfMorSet

Given the local field $L$, this function returns a group $G$
and a map from $G$ to the automorphisms of $L$ over the padics.
"""
function absolute_automorphism_group(L::LocalField)
  aut = absolute_automorphisms(L)
  mult_table = Matrix{Int}(undef, length(aut), length(aut))
  for s = 1:length(aut)
    for i = 1:length(aut)
      mult_table[s, i] = findfirst(isequal(aut[s]*aut[i]), aut)
    end
  end
  G = GrpGen(mult_table)
  return G, GrpGenToNfMorSet(G, aut, L)
end


text:
# file contains default parameters used throughout Iris

# PML/cPML parameters
const EXTINCTION = 1e-10 # extinction in PML layer
const SCALING_ANGLE = .15 # phase in conductivity to accelerate evanscent decay
const BL_DEPTH = .1 # default depth of cPML or PML layer

const DEFAULT_BC = :noBC # default boundary condition (be specified as symbol or string)
const DEFAULT_BL = :noBL # default boundary layer (specified as symbol or string)
const LUPACK = :USolver # default LU package, must be one of USolver, MSolver, PSolver (specified as symbol or string)
const INDEX_OFFSET = 7 # position index offset from midpoint where field is real (SPECTRAL), or real and non-vanishing (LASING/SCPA)
const EQUIVALENT_SOURCE_RELATIVE_CUTOFF = 1e-8 # relative tolerance condition for dropping components of equivalent current
const DEFAULT_CFL_NUMBER = 0.9
const DEFAULT_LINEAR_EIGENSOLVER = :Arpack
const DEFAULT_NONLINEAR_EIGENSOLVER = :NonlinearEigenproblems
const NUM_SUBPIXELS = 5
const ORTHOGONALIZE_OVERLAP_THRESHOLD = .1
const SPA_ACCEPTANCE_PHASE = 5*π/180 # largest phase deviation from real axis considered as "real" for single pole approximation

# 2D unsymmetric shape defaults
const NL_NORMAL_ALGORITHM = :LD_MMA # NLopt algorithm used to find the normal, tangent, and distance to shape perimeter
const NL_NORMAL_XTOL = 1e-5 # NLopt's xtol_rel
const NL_NORMAL_FTOL = 1e-5 # NLopt's ftol_rel
const NL_NORMAL_MAXEVAL = 400 # NLopt's maxeval

# for plotting
const SHAPE_FILL_ALPHA = .2
const SHAPE_COLOR = :red
const BOUNDARY_COLOR = :gray
const BOUNDARY_PML_COLOR = :blue
const BOUNDARY_CPML_COLOR = :red
const BOUNDARY_DIRICHLET_LINETYPE = :solid
const BOUNDARY_NEUMANN_LINETYPE = :dash
const BOUNDARY_FLOQUET_MATCHED_LINETYPE = :dashdot
const BOUNDARY_BC_COLOR = :black
const MARKERSHAPE = :rect
const LINK_SCALE_REDUCTION = 1/2
const MARKERSIZE_SCALE = 135
const PLOT_QUANTILE = .95 # optimize the dynamic range of plotting wavefunctions to cover this percentile
const PLOT_SCALE_FUDGE = 1.3 # multiply plot_quantile result by this factor to get colorbar scale

# Pretty Printing Colors
const PRINTED_COLOR_NUMBER = :light_cyan # numbers
const PRINTED_COLOR_VARIABLE = :cyan # varialbes and functions
const PRINTED_COLOR_WARN = :light_yellow # to grab attention but not necessarily b/c of an error
const PRINTED_COLOR_GOOD = :green # successful solution
const PRINTED_COLOR_DARK = 63 # names of structures that are part of simulation, but do not themselves contain calculations
const PRINTED_COLOR_LIGHT = 171 # names of structures used to calculate and/or containing solutions
const PRINTED_COLOR_BAD = :light_red # indicates error
const PRINTED_COLOR_INSTRUCTION = 244 # parenthetical instructions

"""
PML extinction factor: `$EXTINCTION`
"""
EXTINCTION

"""
PML scaling angle: `$SCALING_ANGLE`
"""
SCALING_ANGLE

"""
default PML depth: `$BL_DEPTH`
"""
BL_DEPTH

"""
default boundary condition: `$DEFAULT_BC`
"""
DEFAULT_BC

"""
default boundary layer: `$DEFAULT_BL`
"""
DEFAULT_BL

"""
default LU package (one of `USolver`, `PSolver`, `MSolver` for UMFPACK, Pardiso, MUMPS): `$LUPACK`
"""
LUPACK

"""
index offset from midpoint where field is real (SPECTRAL), or real and nonzero (LASING/SaturableCPA): `$INDEX_OFFSET`
"""
INDEX_OFFSET

"""
relative tolerance condition for dropping components of equivalent current: `$EQUIVALENT_SOURCE_RELATIVE_CUTOFF`
"""
EQUIVALENT_SOURCE_RELATIVE_CUTOFF

"""
Courant-Friedrichs-Lewy number (N*space-to-time step ratio in time domain simulations in `N` dimensions): `$DEFAULT_CFL_NUMBER`
"""
DEFAULT_CFL_NUMBER

"""
`$DEFAULT_LINEAR_EIGENSOLVER`
"""
DEFAULT_LINEAR_EIGENSOLVER

"""
`$DEFAULT_NONLINEAR_EIGENSOLVER`
"""
DEFAULT_NONLINEAR_EIGENSOLVER

"""
subsampling rate for simulation smoothing: `$NUM_SUBPIXELS`
"""
NUM_SUBPIXELS

#
# # Tessellation defaults
# const TESSELLATION_EXPANSION_FACTOR = 1.025 # to make accessing the tessellation a little more stable
#
# # for interpolating in Simulations
# const N_CUBIC = 10
# const N_QUADRATIC = 6
# const N_BILINEAR = 4
# const N_LINEAR = 3
#
# const NNN_BULK = N_CUBIC # number of nearest neighbors used in bulk
# const NNN_SURFACE = N_LINEAR # number of nearest neighbors for surface
# const NNN_CORNER = N_LINEAR # number of nearest neighbors at corners


text:
using Documenter, Causal
# using DocumenterLaTeX

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

makedocs(
    modules=[Causal],
    sitename="Causal",
    pages=[
        "Home" => "index.md",
        "Modeling and Simulation in Causal" => [
            "modeling_and_simulation/modeling.md",
            "modeling_and_simulation/simulation.md",
        ],
        "Tutorials" => [
            "Model Construction" => "tutorials/model_construction.md",
            "Model Simulation" => "tutorials/model_simulation.md",
            "Algebraic Loops" => "tutorials/algebraic_loops.md",
            "Extending Component Library" => "tutorials/defining_new_components.md",
            "Coupled Systems" => "tutorials/coupled_systems.md",
        ],
        "Manual" => [
            "Utilities" => [
                "manual/utilities/callback.md",
                "manual/utilities/buffers.md"
            ],
            "Connections" => [
                "manual/connections/link.md",
                "manual/connections/pin.md",
                "manual/connections/port.md",
            ],
            "Components" => [
                "ComponentsBase" => [
                    "manual/components/componentsbase/hierarchy.md",
                    "manual/components/componentsbase/evolution.md",
                    "manual/components/componentsbase/interpolation.md",
                ],
                "Sources" => [
                    "manual/components/sources/clock.md",
                    "manual/components/sources/generators.md",
                ],
                "Sinks" => [
                    "manual/components/sinks/sinks.md",
                    "manual/components/sinks/writer.md",
                    "manual/components/sinks/printer.md",
                    "manual/components/sinks/scope.md",
                ],
                "Systems" => [
                    "StaticSystems" => [
                        "StaticSystems" => "manual/components/systems/staticsystems/staticsystems.md",
                        # "Subsystem" => "manual/components/systems/staticsystems/subsystem.md",
                    ],
                    "DynamicSystems" => [
                        "DiscreteSystem" => "manual/components/systems/dynamicsystems/discretesystem.md",
                        "ODESystem" => "manual/components/systems/dynamicsystems/odesystem.md",
                        "DAESystem" => "manual/components/systems/dynamicsystems/daesystem.md",
                        "RODESystem" => "manual/components/systems/dynamicsystems/rodesystem.md",
                        "SDESystem" => "manual/components/systems/dynamicsystems/sdesystem.md",
                        "DDESystem" => "manual/components/systems/dynamicsystems/ddesystem.md",
                    ],
                ]
            ],
            "Models" => [
                "manual/models/taskmanager.md",
                "manual/models/simulation.md",
                "manual/models/model.md",
            ],
            "Plugins" => "manual/plugins/plugins.md",
        ]
    ]
)

deploydocs(;
    repo="github.com/zekeriyasari/Causal.jl.git",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "v#.#.#"]
)


text:
using CSV, DataFrames
using TextAnalysis: NaiveBayesClassifier, fit!, predict

spamdata = DataFrame(CSV.read("/Users/kimfung/Documents/Julia/TextAnalysis/spam.csv"; allowmissing=:none))
global m = NaiveBayesClassifier([:ham, :spam])
for row in eachrow(spamdata)
    if row.v1 == "ham"
        fit!(m, filter(isvalid, row.v2), :ham)
    elseif row.v1 == "spam"
        fit!(m, filter(isvalid, row.v2), :spam)
    end
end

function checkspam(msg::String)
    prediction = predict(m, msg)
    if prediction[:spam] > prediction[:ham]
        println("spam")
    else
        println("ham (not spam)")
    end
end

max_stars_repo_name: fungktt/JuliaTextAnalysis

max_issues_repo_name: fungktt/JuliaTextAnalysis

max_forks_repo_name: fungktt/JuliaTextAnalysis


text:
mutable struct CImGuiToggleButton
    on_name
	off_name
	name
    state
    off_color
    on_color
    hover_color
end
# Outer constructor
CImGuiToggleButton(on_name, off_name) = CImGuiToggleButton(on_name, off_name, off_name, false,
                        (0.1, 0.7, 0.7), (0.3, 0.9, 0.9), (0.2, 0.7, 0.7))
function toggle!(x::CImGuiToggleButton)
    x.state = !x.state
	if x.state
		x.name = x.on_name
	else
		x.name = x.off_name
	end
    return nothing
end

mutable struct PSUChannel
	volt_meas
	curr_meas
	volt_set
	curr_set
	ovp_set
	ocp_set
end
# Outer constructor
PSUChannel() = PSUChannel(repeat(["00.000"], 5)..., "off")

mutable struct PST3201Conf
	active
	C1::PSUChannel
	C2::PSUChannel
	C3::PSUChannel
	output_state
	crt_chan
    crt_func
    crt_value
    ocp_state
	shift_btn
	output_btn
	refresh_btn
end
# Outer constructor
PST3201Conf() = PST3201Conf(false,
					PSUChannel(), PSUChannel(), PSUChannel(),
					"off",
					"C1",
					"voltage",
					0.0,
					"na",
					CImGuiToggleButton(" SHIFT ", " SHIFT "),
					CImGuiToggleButton("OUTPUT", "OUTPUT"),
					CImGuiToggleButton("REFRESH", "REFRESH")
					)

# DMM
mutable struct GDM8246Conf
	active
	primary
	secondary
	range
	crt_chan
    crt_func
    crt_value
	shift_btn
	auto_btn
	refresh_btn
end
# Outer constructor
GDM8246Conf() = GDM8246Conf(false, "000.000", "----.---", "na", "C1", "na", 0.0,
					CImGuiToggleButton(" SHIFT ", " SHIFT "),
					CImGuiToggleButton("  AUTO  ", "  AUTO  "),
					CImGuiToggleButton("REFRESH", "REFRESH")
					)

# FGEN
mutable struct GFG3015Conf
	active
	idx_sel
	crt_chan
    crt_func
    crt_value
	amplit_unit
	freq
	amplit
	offs
	duty
	shift_btn
	refresh_btn
end
# Outer constructor
GFG3015Conf() = GFG3015Conf(false, "1", "C1", "na", 0.0, "Vpp", "na","na","na","na",
					CImGuiToggleButton(" SHIFT ", " SHIFT "),
					CImGuiToggleButton("REFRESH", "REFRESH")
					)

#scope
mutable struct TDS2002BConf
	active
	CH1_Volt_div
	CH1_Offset
	CH2_Volt_div
	CH2_Offset
	Time_div
	Trigger_source
	Trigger_level
	Trigger_mode           #triggered or imediate
	CH1_Volt_div_new
	CH1_Offset_new
	CH2_Volt_div_new
	CH2_Offset_new
	Time_div_new
	Trigger_source_new
	Trigger_level_new
	Trigger_mode_new
	Acquire
	Acquire_btn
	t
	y1
	y2
	Measurement_source1
	Measurement_Type1
	Measurement_Value1
	Measurement_source2
	Measurement_Type2
	Measurement_Value2
	Measurement_source3
	Measurement_Type3
	Measurement_Value3
	Measurement_source4
	Measurement_Type4
	Measurement_Value4
	Measurement_source5
	Measurement_Type5
	Measurement_Value5
	end
					# Outer constructor
	TDS2002BConf() = TDS2002BConf(false,
				"1.00",
				"0.00",
				"1.00",
				"0.00",
				"0.001",    #in seconds
				"CH1",
				"01.00",
				"AUTO",
				00.00,   #CH1_Volt_div_new
				00.00,   #CH1_Offset_new
				00.00,   #CH2_Volt_div_new
				00.00,   #CH2_Offset_new
				00.000,    #Time_div_new
				"CH1",
				00.00,
				"AUTO" ,
				"off",
				CImGuiToggleButton("ACQUIRE", "ACQUIRE"),
				collect(1:300),
				randn(300),
				randn(300),
				"CH1",
				"Maximum",
				"00.00",
				"CH1",
				"Maximum",
				"00.00",
				"CH1",
				"Maximum",
				"00.00",
				"CH1",
				"Maximum",
				"00.00",
				"CH1",
				"Maximum",
				"00.00"
				)

# RELAYS matrix
mutable struct RelaysConf
	active
	C1
	C2
	C3
	C4
	C5
	C6
	C7
	C8
	C9
    C10
    C11
    C12
    C13
    C14
    C15
    C16
	Refresh
end
# Outer constructor
RelaysConf() = RelaysConf(false,
					CImGuiToggleButton("C1 ##C1", "C1 ##C1"),
					CImGuiToggleButton("C2 ##C2", "C2 ##C2"),
					CImGuiToggleButton("C3 ##C3", "C3 ##C3"),
					CImGuiToggleButton("C4 ##C4", "C4 ##C4"),
					CImGuiToggleButton("C5 ##C5", "C5 ##C5"),
					CImGuiToggleButton("C6 ##C6", "C6 ##C6"),
					CImGuiToggleButton("C7 ##C7", "C7 ##C7"),
					CImGuiToggleButton("C8 ##C8", "C8 ##C8"),
					CImGuiToggleButton("C9 ##C9", "C9 ##C9"),
                    CImGuiToggleButton("C10 ##C10", "C10 ##C10"),
                    CImGuiToggleButton("C11 ##C11", "C11 ##C11"),
                    CImGuiToggleButton("C12 ##C12", "C12 ##C12"),
                    CImGuiToggleButton("C13 ##C13", "C13 ##C13"),
                    CImGuiToggleButton("C14 ##C14", "C14 ##C14"),
                    CImGuiToggleButton("C15 ##C15", "C15 ##C15"),
                    CImGuiToggleButton("C16 ##C16", "C16 ##C16"),
					CImGuiToggleButton("ON ##Refresh", "OFF##Refresh")
					)

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
"""
    blockparse(str::String)

Initial Meta.parse() on entire input. return funcs and initialize Dicts of argtypes and imports.

# Details:
- funcs: a dict with function,expression as key,values
- argtypes: a dict with with function,argtypes as key,values
- imports: a dict with with function,importstring as key,values
"""
function blockparse(str::String)
    block = pruneLineNumberNode(Meta.parse("begin $str end").args)
    funcs = Dict(blockfuncnames(block) .=> blockfuncexpressions(block))
    argtypes = Dict()
    imports = Dict()

    # fill in top level argtypes for function calls
    for ex in block
        if ex.head == :(call)
            name = ex.args[1]
            args = ex.args[2:end]
            # typeof(eval(arg)) instead of typeof(arg) allows top level call with rand() etc.
            types = typeof.([Base.eval(Evalscope, arg) for arg in args])
            argtypes[name] = types
        end
    end
    return funcs, argtypes, imports
end

pruneLineNumberNode(item) = item
pruneLineNumberNode(items::Array) = [pruneLineNumberNode(item) for item in items if !isa(item, LineNumberNode)]
pruneLineNumberNode(item::Expr) = item.head == :(block) ? pruneLineNumberNode(item.args) : Expr(item.head, pruneLineNumberNode(item.args)...)

function isfunction(ex::Expr)
    if ex.head == Symbol("function") || ex.head == :(=)
        return true
    else
        return false
    end
end
isfunction(ex) = false

blockfuncnames(block) = [ex.args[1].args[1] for ex in block if isfunction(ex)]
blockfuncexpressions(block) = [ex.args[2] for ex in block if isfunction(ex)]

"""
    codeinfo(func::Symbol, argtypes::Array)

Essentially code_typed() with optimize=false

# Details:
Not optimizing gives a much simpler ast without phinodes and boundchecks.
"""
function codeinfo(func, argtypes::Array)
    ct = code_typed(Base.eval(Evalscope, func), Tuple{argtypes...}; optimize=false, debuginfo=:none)[1] # none, source
    cinfo = ct[1]
    Rtype = ct[2]

    for (i, st) in enumerate(cinfo.slottypes)
        if isa(st, Union)
            # iterator variables are union of [nothing, [value,index]]
            # getfield(st,2) gets the [value,index] tuple
            cinfo.slottypes[i] = getfield(st, 2).parameters[1]
            push!(cinfo.slotnames, Symbol("_$(i)i"))
            push!(cinfo.slottypes, getfield(st, 2).parameters[2])
        elseif !isa(st, Compiler.Const) && length(st.parameters) > 1 && istuple(st)
            # wasm can now return multiple values, but tuples dont exist
            # so create individual variables representing the tuple parameters
            # (but keep the "incorrect" tuple type of the original tuple slot)

            for j = 2:length(st.parameters)
                push!(cinfo.slotnames, Symbol("$(cinfo.slotnames[i])$(j)"))
                push!(cinfo.slottypes, st.parameters[j])
            end
            #cinfo.slottypes[i] = st.parameters[1] #actually dont keep the "tuple" type of the first
        end

        # make sure slots have names
        if string(cinfo.slotnames[i]) == ""
            cinfo.slotnames[i] = Symbol("_$i")
        end
    end

    for (i, vt) in enumerate(cinfo.ssavaluetypes)
        if isa(vt, Union)
            cinfo.ssavaluetypes[i] = Tuple{getfield(vt, 2).parameters...}
        elseif isa(vt, Compiler.Const)
            cinfo.ssavaluetypes[i] = typeof(vt.val)
        elseif isa(vt, PartialStruct)
            cinfo.ssavaluetypes[i] = getfield(vt, 1)
        end
    end
    return cinfo, Rtype
end

istuple(slottype) = all([typeof(p) <: DataType for p in slottype.parameters])


text:
# TODO: type for dm::DemandMatrix? or TrafficMatrix?

basic_routing_model_unitary(m::Model, rd::RoutingData) =
    basic_routing_model_unitary(m, rd, rd.model_type.type)
basic_routing_model_unitary(m::Model, rd::RoutingData, ::Any...) =
    error("Not implemented: $(rd.model_type)")
basic_routing_model_unscaled(m::Model, rd::RoutingData, ::Any...) =
    error("Not implemented: $(rd.model_type)")

function total_flow_in_edge(rm::RoutingModel, e::Edge)
    @assert rm.scaled_flows == UnscaledFlows
    return total_flow_in_edge(rm, e, rm.data.model_type.type, Val(UnscaledFlows))
end

function total_flow_in_edge(rm::RoutingModel, e::Edge, dm::Union{Dict{Edge{Int}, Float64}, Nothing})
    # Switch on the flow scale: don't pass a traffic matrix when it's not needed.
    if rm.scaled_flows == UnitaryFlows
        if dm === nothing
            dm = rm.data.traffic_matrix
        end
        @assert dm !== nothing
        return total_flow_in_edge(rm, e, dm, rm.data.model_type.type, Val(UnitaryFlows))
    elseif rm.scaled_flows == UnscaledFlows
        # dm not needed.
        return total_flow_in_edge(rm, e, rm.data.model_type.type, Val(UnscaledFlows))
    else
        error("Not implemented: $(rm.data.model_type)")
    end
end
total_flow_in_edge(rm::RoutingModel, edge::Edge, ::Any...) =
    error("Not implemented: $(rm.data.model_type)")
total_flow_in_edge(rm::RoutingModel, edge::Edge, dm::Union{Dict{Edge{Int}, Float64}, Nothing}, ::Any...) =
    error("Not implemented: $(rm.data.model_type)")

# Interpretation of the dm argument: demand -> value

function capacity_constraints(rm::RoutingModel, dm::Union{Dict{Edge{Int}, Float64}, Nothing}=nothing)
    # Capacity constraints, one per edge.
    @assert isempty(rm.constraints_capacity)
    rm.constraints_capacity = Dict{Edge{Int}, Any}()
    for e in edges(rm.data)
        flow = total_flow_in_edge(rm, e, dm)
        if ! iszero(flow)
            rm.constraints_capacity[e] = @constraint(rm.model, flow <= capacity(rm, e))
        end
    end
end

function mu_capacity_constraints(rm::RoutingModel, dm::Union{Dict{Edge{Int}, Float64}, Nothing}=nothing)
    # Difference with capacity_constraints: <= mu * capacity, instead of simply capacity.

    # Create the mu variable if it does not exist yet.
    if rm.mu === nothing
        rm.mu = @variable(rm.model, lower_bound=0.0)
        set_name(rm.mu, "mu")
    end

    # Capacity constraints, one per edge.
    @assert isempty(rm.constraints_capacity)
    rm.constraints_capacity = Dict{Edge{Int}, Any}()
    for e in edges(rm.data)
        flow = total_flow_in_edge(rm, e, dm)
        if ! iszero(flow)
            rm.constraints_capacity[e] = @constraint(rm.model, flow <= rm.mu * capacity(rm, e))
        end
    end
end


text:
#=
Note: This file was taken from the TimeWarp.jl package

The TimeWarp.jl package is licensed under the MIT "Expat" License:

> Copyright (c) 2016: Alex Williams and Joseph Fowler
>
> 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.
=#

using Distances
include("sequence.jl")
include("utils.jl")

#####################################
#     Basic interface functions     #
#####################################

"""
    cost,i1,i2 = dtw(seq1, seq2, [dist=SqEuclidean])

Find a set of indices (`i1`,`i2`) that align two time series (`seq1`,`seq2`) by
dynamic time warping. Also returns the distance (after warping) according to
the SemiMetric `dist`, which defaults to squared Euclidean distance (see
Distances.jl). If `seq1` and `seq2` are matrices, each column is considered
an observation.
"""
function dtw(
        seq1::Sequence,
        seq2::Sequence,
        dist::SemiMetric=SqEuclidean()
    )

    D = dtw_cost_matrix(seq1, seq2, dist)
    return trackback(D)
end

"""
    cost,i1,i2 = dtw(seq1,seq2,i2min,i2max,[dist=SqEuclidean])

Do DTW to align `seq1` and `seq2` confined to a window. Vectors `i2min` and
`i2max` specify (inclusive) lower and upper bounds for `seq2` for each index in
`seq1`. Thus, `i2min` and `i2max` are required to be the same length as `seq1`.
"""
function dtw(
        seq1::Sequence,
        seq2::Sequence,
        i2min::AbstractVector,
        i2max::AbstractVector,
        dist::SemiMetric = SqEuclidean()
    )

    D = dtw_cost_matrix(seq1, seq2, i2min, i2max, dist)
    return trackback(D)
end

##############################
#  Cost matrix computations  #
##############################
function dtw_cost_matrix{N,T}(
        seq1::Sequence{N,T},
        seq2::Sequence{N,T},
        dist::SemiMetric=SqEuclidean()
    )
    # Build the cost matrix
    m = length(seq2)
    n = length(seq1)
    D = zeros(T, m, n)

    # Initialize first column and first row
    D[1,1] = evaluate(dist, seq1[1], seq2[1])
    for r=2:m
        D[r,1] = D[r-1,1] + evaluate(dist, seq1[1], seq2[r])
    end
    for c=2:n
        D[1,c] = D[1,c-1] + evaluate(dist, seq1[c], seq2[1])
    end

    # Complete the cost matrix
    for c=2:n
        for r=2:m
            best_neighbor_cost = min(D[r-1,c], D[r-1,c-1], D[r,c-1])
            D[r,c] = best_neighbor_cost + evaluate(dist, seq1[c], seq2[r])
        end
    end

    return D
end

function dtw_cost_matrix{N,T,U<:Integer}(
        seq1::Sequence{N,T},
        seq2::Sequence{N,T},
        i2min::AbstractVector{U},
        i2max::AbstractVector{U},
        dist::SemiMetric = SqEuclidean()
    )
    m = length(seq2) # of rows in cost matrix
    n = length(seq1) # of columns in cost matrix
    n == length(i2min) || throw(ArgumentError("i2min does not match length of seq1."))
    n == length(i2max) || throw(ArgumentError("i2max does not match length of seq1."))
    1 == i2min[1] || throw(ArgumentError("i2min must start at 1."))
    m == i2max[end] || throw(ArgumentError("i2max must end at length(seq2)."))

    # Build the (n x m) cost matrix into a WindowedMatrix, because it's ragged.
    # That type gives efficient storage with convenient [r,c] indexing and returns
    # Inf when accessed outside the window.
    D = WindowedMatrix(i2min, i2max, Inf)

    # First column first
    D[1,1] = evaluate(dist, seq1[1], seq2[1])
    for r=2:i2max[1]
        D[r,1] = D[r-1,1]  + evaluate(dist, seq1[1], seq2[r])
    end

    # Complete the cost matrix from columns 2 to m.
    for c=2:n
        for r=i2min[c]:i2max[c]
            best_neighbor_cost = min(D[r-1,c], D[r-1,c-1], D[r,c-1])
            D[r,c] = best_neighbor_cost + evaluate(dist, seq1[c], seq2[r])
        end
    end

    return D
end

########################################
#  Find Best Path through Cost Matrix  #
########################################

"""
    cost,cols,rows = trackback(D::Matrix)

Given the cost matrix `D`, computes the optimal track from end to beginning.
Returns `cols` and `rows` which are vectors respectively holding the track.
"""
function trackback{T<:Number}(D::AbstractMatrix{T})

    # initialize trackback throught rows/columns
    r,c = size(D)
    rows,cols = Int[r],Int[c]

    # estimate that we'll need N⋅logN elements
    N = max(r,c)
    sz = 2*N
    sizehint!(rows,sz)
    sizehint!(cols,sz)

    # do trackback
    while r > 1 && c > 1
        tb,r,c = indmin3(D[r-1,c-1], D[r-1,c], D[r,c-1], r, c)
        push!(rows,r)
        push!(cols,c)
    end
    # Possibly either r>1 or c>1 at this point (but not both).
    # Add the unfinished part of the track to reach [1,1]
    for r=r-1:-1:1
        push!(rows,r)
        push!(cols,1)
    end
    for c=c-1:-1:1
        push!(rows,1)
        push!(cols,c)
    end
    return D[end,end], reverse(cols), reverse(rows)
end

###############################################
# Wrapper functions for abstract array inputs #
###############################################

function dtw(
        seq1::AbstractArray,
        seq2::AbstractArray,
        dist::SemiMetric=SqEuclidean()
    )
    dtw(Sequence(seq1),Sequence(seq2),dist)
end

function dtw(
        seq1::AbstractArray,
        seq2::AbstractArray,
        i2min::AbstractVector,
        i2max::AbstractVector,
        dist::SemiMetric = SqEuclidean()
    )
    dtw(Sequence(seq1),Sequence(seq2),i2min,i2max,dist)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(maybe_evaluate_builtin), JuliaStackFrame, Expr})
    precompile(Tuple{typeof(getargs), Vector{Any}, JuliaStackFrame})
    precompile(Tuple{typeof(get_call_framecode), Vector{Any}, JuliaFrameCode, Int})
    for f in (evaluate_call!,
              evaluate_foreigncall!,
              evaluate_methoddef!,
              evaluate_structtype!,
              evaluate_abstracttype!,
              evaluate_primitivetype!)
        precompile(Tuple{typeof(f), Compiled, JuliaStackFrame, Expr, JuliaProgramCounter})
        precompile(Tuple{typeof(f), Vector{JuliaStackFrame}, JuliaStackFrame, Expr, JuliaProgramCounter})
    end
    precompile(Tuple{typeof(lookup_global_refs!), Expr})
    precompile(Tuple{typeof(lookup_or_eval), Compiled, JuliaStackFrame, Any, JuliaProgramCounter})
    precompile(Tuple{typeof(lookup_or_eval), Vector{JuliaStackFrame}, JuliaStackFrame, Any, JuliaProgramCounter})
    precompile(Tuple{typeof(eval_rhs), Compiled, JuliaStackFrame, Expr, JuliaProgramCounter})
    precompile(Tuple{typeof(eval_rhs), Vector{JuliaStackFrame}, JuliaStackFrame, Expr, JuliaProgramCounter})
    precompile(Tuple{typeof(_step_expr!), Compiled, JuliaStackFrame, Any, JuliaProgramCounter, Bool})
    precompile(Tuple{typeof(_step_expr!), Vector{JuliaStackFrame}, JuliaStackFrame, Any, JuliaProgramCounter, Bool})
    for f in (finish!, finish_and_return!)
        precompile(Tuple{typeof(f), Compiled, JuliaStackFrame, JuliaProgramCounter, Bool})
        precompile(Tuple{typeof(f), Vector{JuliaStackFrame}, JuliaStackFrame, JuliaProgramCounter, Bool})
    end
    precompile(Tuple{typeof(split_expressions), Module, Expr})
    precompile(Tuple{typeof(prepare_thunk), Module, Expr})
    precompile(Tuple{typeof(prepare_locals), JuliaFrameCode, Vector{Any}})
    precompile(Tuple{typeof(prepare_args), Any, Vector{Any}, Vector{Any}})
    precompile(Tuple{typeof(prepare_call), Any, Vector{Any}})
    precompile(Tuple{typeof(Core.kwfunc(prepare_call)), NamedTuple{(:enter_generated,),Tuple{Bool}}, typeof(JuliaInterpreter.prepare_call), Function, Vector{Any}})
    precompile(Tuple{typeof(build_frame), JuliaFrameCode, Vector{Any}, Core.SimpleVector})
    precompile(Tuple{typeof(extract_args), Module, Expr})
    precompile(Tuple{typeof(enter_call), Int, Int})
    precompile(Tuple{typeof(enter_call_expr), Expr})
    precompile(Tuple{typeof(copy_codeinfo), Core.CodeInfo})
    precompile(Tuple{typeof(optimize!), Core.CodeInfo, Module})
    precompile(Tuple{typeof(set_structtype_const), Module, Symbol})
    precompile(Tuple{typeof(namedtuple), Vector{Any}})
    precompile(Tuple{typeof(resolvefc), Any})
    precompile(Tuple{typeof(check_isdefined), JuliaStackFrame, Expr})
    precompile(Tuple{typeof(check_isdefined), JuliaStackFrame, Symbol})
    precompile(Tuple{typeof(find_used), Core.CodeInfo})
    precompile(Tuple{typeof(do_assignment!), JuliaStackFrame, Int, Int})
    precompile(Tuple{typeof(pc_expr), JuliaStackFrame, Nothing})
end


text:
module GaussianFermions

using MatrixProductStateTools
using SymTensors
using LinearAlgebra

include("gmps.jl")
include("gmps2mps.jl")

export GaussianMPS
export corrmat2gmps
export gmps2mps

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
function root_sum_of_squares(signals::AbstractArray{<: Number}, dim::Val{N}=Val(1)) where N
	result = dropdims(sum(abs2, signals; dims=N); dims=N)
	@. result = sqrt(result)
end

function roemer(
	signals::AbstractArray{<: Number},
	sensitivity::AbstractArray{<: Number},
	dim::Val{N}=Val(1)
) where N
	dropdims(
		sum(conj.(sensitivity) .* signals; dims=N) ./ sum(abs2, sensitivity; dims=N);
		dims=N
	)
end

"""
	plan_sensitivities(sensitivities::Union{AbstractMatrix{<: Number}, AbstractArray{<: Number, 3}}) = S

sensitivities[spatial dimensions..., channels]
shape = (spatial dimensions, other dimension)

"""
function plan_sensitivities(
	sensitivities::AbstractArray{<: Number, N},
	shape::NTuple{N, <: Integer}
) where N

	# Get dimensions
	# TODO: check spatial dims
	other_dims = shape[N]
	shape = shape[1:N-1]
	shape_s = size(sensitivities)
	@assert shape == shape_s[1:N-1]
	channels = shape_s[N]
	spatial_dimensions = prod(shape)
	input_dimension = spatial_dimensions * other_dims
	output_dimension = input_dimension * channels

	# Reshape
	sensitivities = reshape(sensitivities, spatial_dimensions, channels, 1)
	conj_sensitivities = conj.(sensitivities)

	S = LinearMap{ComplexF64}(
		x::AbstractVector{<: Complex} -> begin
			Sx = sensitivities .* reshape(x, spatial_dimensions, 1, other_dims)
			vec(Sx)
		end,
		y::AbstractVector{<: Complex} -> begin
			y = reshape(y, spatial_dimensions, channels, other_dims)
			SHy = sum(conj_sensitivities .* y; dims=2)
			vec(SHy)
		end,
		output_dimension, input_dimension
	)
	return S
end

max_stars_repo_name: felixhorger/MRIRecon

max_issues_repo_name: felixhorger/MRIRecon

max_forks_repo_name: felixhorger/MRIRecon


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

using ForwardDiff
using Printf
using LinearAlgebra

const Vec = Vector{Float64}
const Mat = Matrix{Float64}
const iMat = Matrix{Int}

N_CAM_PARAMS = 11
ROT_IDX = 1
C_IDX = 4
F_IDX = 7
X0_IDX = 8
RAD_IDX = 10

##################### objective #############################

function rodrigues_rotate_point(rot :: Vector{T}, X :: Vector{T}) where T
  sqtheta = sum(rot.*rot)
  if sqtheta > 1e-10
    theta = sqrt(sqtheta)
    costheta = cos(theta)
    sintheta = sin(theta)
    theta_inverse = 1. / theta

    w = theta_inverse * rot
    w_cross_X = cross(w,X)
    tmp = dot(w,X) * (1. - costheta)

    X*costheta + w_cross_X * sintheta + w * tmp
  else
    X + cross(rot,X)
  end
end

function radial_distort(rad_params,proj)
  rsq = sum(proj.*proj)
  L = 1. + rad_params[1]*rsq + rad_params[2]*rsq*rsq
  proj*L
end

function project(cam,X)
  Xcam = rodrigues_rotate_point(cam[ROT_IDX:ROT_IDX+2], X - cam[C_IDX:C_IDX+2])
  distorted = radial_distort(cam[RAD_IDX:RAD_IDX+1], Xcam[1:2]/Xcam[3])
  distorted*cam[F_IDX] + cam[X0_IDX:X0_IDX+1]
end

function compute_reproj_err(cam,X,w,feat)
  return w*(project(cam,X) - feat)
end

function ba_objective(cams,X,w,obs,feats)
  reproj_err = similar(feats)
  for i in 1:size(feats,2)
    reproj_err[:,i] = compute_reproj_err(cams[:,obs[1,i]],X[:,obs[2,i]],w[i],feats[:,i])
  end
  w_err = 1.0 .- w.*w
  (reproj_err, w_err)
end

#################### derivatives extra ##########################

function pack(cam,X,w)
  [cam[:];X[:];w]
end

function unpack(packed)
  packed[1:end-4],packed[end-3:end-1],packed[end]
end

function compute_w_err(w)
    1.0 - w*w
end
compute_w_err_d = x -> ForwardDiff.derivative(compute_w_err, x)

function compute_reproj_err_d(params, feat)
  cam, X, w = unpack(params)
  compute_reproj_err(cam,X,w,feat)
end

function compute_ba_J(cams,X,w,obs,feats)
  p = size(obs,2)
  reproj_err_d = zeros(2*p, N_CAM_PARAMS + 3 + 1)
  for i in 1:p
    compute_reproj_err_d_i = x -> compute_reproj_err_d(x, feats[:,i])
    idx = (2*(i-1))+1
    reproj_err_d[idx:idx+1,:] = ForwardDiff.jacobian(compute_reproj_err_d_i,
              pack(cams[:,obs[1,i]],X[:,obs[2,i]],w[i]))
  end
  w_err_d = zeros(1,p)
  for i in 1:p
    w_err_d[i] = compute_w_err_d(w[i])
  end
  (reproj_err_d, w_err_d)
end

# precompile as much as possible
precompile(rodrigues_rotate_point,(Vec,Vec))
precompile(radial_distort,(Vec,Vec))
precompile(project,(Vec,Vec))
precompile(compute_reproj_err,(Vec,Vec,Float64,Vec))
precompile(ba_objective,(Mat,Mat,Vec,iMat,Mat))

precompile(compute_ba_J,(Mat,Mat,Vec,iMat,Mat))
precompile(compute_reproj_err_d,(Vec,Vec))
precompile(compute_w_err_d,(Vec,))

#################### IO ############################
function read_ba_instance(fn)
  fid = open(fn)
  lines = readlines(fid)
  close(fid)
  line=split(lines[1]," ")
  n = parse(Int,line[1])
  m = parse(Int,line[2])
  p = parse(Int,line[3])
  off = 2

  one_cam = zeros(Float64,N_CAM_PARAMS,1)
  line=split(lines[off]," ")
  for i in 1:N_CAM_PARAMS
    one_cam[i] = parse(Float64,line[i])
  end
  cams = repeat(one_cam,1,n)
  off += 1

  one_X = zeros(Float64,3,1)
  line=split(lines[off]," ")
  for i in 1:3
    one_X[i] = parse(Float64,line[i])
  end
  X = repeat(one_X,1,m)
  off += 1

  one_w = parse(Float64,lines[off])
  w = repeat([one_w],1,p)
  off += 1

  one_feat = zeros(Float64,2,1)
  line=split(lines[off]," ")
  for i in 1:2
    one_feat[i] = parse(Float64,line[i])
  end
  feats = repeat(one_feat,1,p)

  camIdx = 1
  ptIdx = 1
  obs = zeros(Int,2,p)
  for i in 1:p
    obs[1,i] = camIdx
    obs[2,i] = ptIdx
    camIdx = (camIdx%n) + 1
    ptIdx = (ptIdx%m) + 1
  end

  (cams,X,w,obs,feats)
end

function write_times(fn,tf,tJ)
  @printf "Writing to %s\n" fn
  fid = open(fn,"w")
  @printf fid "%f %f\r\n" tf tJ
  @printf fid "tf tJ\r\n"
  close(fid)
end

##################### run it ###########################

# Read instance
function main(ARGS)
  dir_in = ARGS[1]
  dir_out = ARGS[2]
  fn = ARGS[3]
  nruns_f = parse(Int,ARGS[4])
  nruns_J = parse(Int,ARGS[5])

  fn_in = string(dir_in, fn)
  fn_out = string(dir_out, fn)

  cams,X,w,obs,feats = read_ba_instance(string(fn_in,".txt"))

  # compute once for precompilation
  reproj_err, w_err = ba_objective(cams,X,w,obs,feats)
  # Time runs
  tf = @elapsed for i in 1:nruns_f
    reproj_err, w_err = ba_objective(cams,X,w,obs,feats)
  end
  tf /= nruns_f
  @printf "tf: %g\n" tf
  #println(reproj_err)
  #println(w_err)

  # compute once for precompilation
  J = compute_ba_J(cams,X,w,obs,feats)
  # Time runs
  tJ = @elapsed for i in 1:nruns_J
    J = compute_ba_J(cams,X,w,obs,feats)
  end
  tJ /= nruns_J;
  @printf "tJ: %g\n" tJ
  #println("J:")
  #println(J)

  name = "Julia"

  #write_J(string(fn_out,"_J_",name,".txt"),J)
  write_times(string(fn_out,"_times_",name,".txt"),tf,tJ)
end

if isinteractive()
  main([raw"C:\dev\github\autodiff/data/ba/", raw"C:\dev\github\autodiff/tmp/Release/ba/Julia/", "ba1", "1", "1", "60"])
else
  main(ARGS)
end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
# chart patterns functions
@testset "Chart Patterns" begin
    Random.seed!(SEED)
    @testset "Array" begin
        X = cumsum(randn(N, 3), dims=1)
        tmp = renko(X, use_atr=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = renko(X, use_atr=false)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    planaritytest, outerplanaritytest,
    isplanar, isouterplanar

PlanarityTestDS = Vector{Vector{Vector{Int}}}

struct PlanarityTestState
    graph::PlainGraph

    rank::Dict{Int,Int} # tree node n, dfsindex(n)
    inedge::Dict{Int,Int} # tree node n, inedge(n)
    isthick::Dict{Int,Bool} # edge e, isthick(e)

    source::Dict{Int,Int} # edge e, source node dfsindex
    loworder::Dict{Int,Vector{Int}} # tree edge e, list of outedge o in low(o) order
    cotree::Dict{Int,Int} # cotree edge e, low(e)
    treeedge::Vector{Int} # tree edges in dfs order

    function PlanarityTestState(graph)
        new(graph, Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), [])
    end
end

function dfs!(state::PlanarityTestState, u::Int)
    state.rank[u] = length(state.rank) + 1
    buckets = Dict{Int,Vector{Int}}() # low(e), edges
    lows = Dict{Int,Int}() # edge, lownode(e)
    for (inc, adj) in neighbors(state.graph, u)
        inc == get(state.inedge, u, nothing) && continue # Predecessor
        haskey(state.cotree, inc) && continue # Visited
        state.source[inc] = state.rank[u]
        if !haskey(state.rank, adj)
            # Tree node
            push!(state.treeedge, inc)
            state.inedge[adj] = inc
            low = dfs!(state, adj)
        else
            # Cotree node
            low = adj
            state.cotree[inc] = state.rank[adj]
            state.isthick[inc] = false
        end
        low === nothing && continue
        lows[inc] = low
        rank = state.rank[low]
        haskey(buckets, rank) || (buckets[rank] = Int[])
        push!(buckets[rank], inc)
    end
    u in keys(state.inedge) || return # Root
    # low(e) ordering
    inedge = state.inedge[u]
    fringes = [i for i in values(state.cotree) if i < state.rank[u]]
    state.isthick[inedge] = length(unique(fringes)) > 1
    sorted = Int[]
    for k in sort(collect(keys(buckets)))
        append!(sorted, sort(buckets[k], by=x->state.isthick[x]))
    end
    state.loworder[inedge] = sorted
    isempty(lows) && return # No cycles
    return lows[state.loworder[inedge][1]]
end

function merge!(ds1::PlanarityTestDS, ds2::PlanarityTestDS, cotree::Dict{Int,Int})
    if isempty(ds1)
        append!(ds1, ds2)
        @debug "Trunk"
        return true
    end
    @assert !isempty(ds2)
    b1 = cotree[ds1[1][1][1]]
    b2 = cotree[ds2[1][1][1]]
    newds = PlanarityTestDS()
    lowtop = 0
    while lowtop <= b2 && !isempty(ds1)
        pair = popfirst!(ds1)
        topl = isempty(pair[1]) ? 0 : cotree[pair[1][end]]
        topr = isempty(pair[2]) ? 0 : cotree[pair[2][end]]
        if topl > b2 && topr > b2
            @debug "Not planer: j+1 case"
            return false
        end
        lowtop = max(topl, topr)
        @assert lowtop != 0
        newl = topl >= topr ? 1 : 2
        newr = topl >= topr ? 2 : 1
        push!(newds, [pair[newl], pair[newr]])
    end
    # newds[end] is now j + 1
    if lowtop <= b2
        push!(newds, [[], []])
    end

    fused1 = Int[]
    while !isempty(ds1)
        pair = popfirst!(ds1)
        if !isempty(pair[1]) && !isempty(pair[2])
            @debug "Not planer: should be trichromatic"
            return false
        end
        append!(fused1, pair[1], pair[2])
    end
    append!(newds[end][1], fused1)

    fused2 = Int[]
    for (i, pair) in enumerate(ds2)
        if !isempty(pair[1]) && !isempty(pair[2])
            @debug "Not planer: dichromatic branch cannot be merged"
            return false
        end
        cell = i == 1 && b1 == b2 ? newds[1][1] : fused2
        append!(cell, pair[1], pair[2])
    end
    whichtofuse = isempty(newds[end][1]) ? 1 : 2
    append!(newds[end][whichtofuse], fused2)
    if isempty(newds[end][1]) && isempty(newds[end][2])
        # ds2 has only a b1 == b2 edge
        pop!(newds)
    end
    # Already ds1 has used up
    append!(ds1, newds)
    @debug "Merged: ds $(ds1)"
    return true
end

function planaritytest(graph::OrderedGraph)
    # Do DFS to determine treeedge, cotree, loworder, source
    state = PlanarityTestState(plaingraph(graph))
    nodes = nodeset(graph)
    while !isempty(nodes)
        n = pop!(nodes)
        dfs!(state, n)
        setdiff!(nodes, keys(state.rank))
    end
    @debug state

    # L-R check
    ds = Dict(c => Vector{Vector{Int}}[[[c], []]] for c in keys(state.cotree))
    while !isempty(state.treeedge)
        e = pop!(state.treeedge)
        ds[e] = PlanarityTestDS()
        @debug "inedge: $(e)"
        for i in state.loworder[e]
            @debug "stem: $(i), ds: $(ds[i])"
            if !merge!(ds[e], ds[i], state.cotree)
                return false
            end
        end
        finished = filter(x -> state.cotree[x] == state.source[e], keys(state.cotree))
        @debug "removed: $(finished)"
        if !isempty(ds[e])
            setdiff!(ds[e][end][1], finished)
            setdiff!(ds[e][end][2], finished)
            if isempty(ds[e][end][1]) && isempty(ds[e][end][2])
                pop!(ds[e])
            end
        end
    end
    return true
end

function outerplanaritytest(graph::UndirectedGraph)
    # Add a node connects all other nodes
    newg = plaingraph(graph)
    newnode = addnode!(newg)
    for n in nodeset(graph)
        addedge!(newg, newnode, n)
    end
    return planaritytest(newg)
end

"""
    isplanar(graph::UndirectedGraph) -> Bool

Return whether the graph is planar.

# Reference

1. de Fraysseix, H., & Ossona de Mendez, P. (2012). Trémaux trees and planarity.
   European Journal of Combinatorics, 33(3), 279–293.
   https://doi.org/10.1016/j.ejc.2011.09.012
1. Trémaux trees and planarity. https://arxiv.org/abs/math/0610935
"""
function isplanar(graph::UndirectedGraph)
    if isdefined(graph, :cache) && haskey(graph.cache, :isouterplanar)
        graph.cache[:isouterplanar] && return true
    end
    edgecount(graph) < 9 && return true
    edgecount(graph) > nodecount(graph) * 3 - 6 && return false
    return planaritytest(graph)
end

"""
    isouterplanar(graph::UndirectedGraph) -> Bool

Return whether the graph is outerplanar. The outerplanarity test is based on
a planarity test (see [`isplanar`](@ref)).
"""
function isouterplanar(graph::UndirectedGraph)
    if isdefined(graph, :cache) && haskey(graph.cache, :isplanar)
        graph.cache[:isplanar] || return false
    end
    edgecount(graph) < 6 && return true
    edgecount(graph) > nodecount(graph) * 2 - 6 && return false
    return outerplanaritytest(graph)
end


text:
using Distributions, Optim, LinearAlgebra, Random, Plots

# Define Thinning Operators
import Base.∘

# Standard Binomial Thinning
"""
    ∘
Thinning operator for p∘X, d∘X or p∘d.

* `p`: Thinning probability (∈ [0, 1])
* `X`: Positive integer
* `d`: Discrete distribution

Random numbers are generated according to (binomial) thinning.

`p∘X` computes ``\\sum_{i = 1}^X{Z_i}`` with ``Z_i\\sim\\text{Bin}(1, p)``

`d∘X` computes ``\\sum_{i = 1}^X{Z_i}`` with ``Z_i\\sim`` `d`

`p∘d` computes ``\\sum_{i = 1}^X{Z_i}`` with ``Z_i\\sim\\text{Bin}(1, n)`` and ``X\\sim`` `d`
"""
p::Float64 ∘ X::Integer = begin
    if (0 <= p <= 1) & (X >= 0)
        rand(Binomial(X, p))
    else
        if !(0 <= p <= 1)
            error("Invalid thinning probability.")
        end
        if X < 0
            error("X must be a positive integer.")
        end
    end
end

# Generic Thinning for any distribution as summands
∘(d::DiscreteDistribution, X::Integer) = sum(rand(d, X))

# Binomial Thinning with distribution for sum limit
p::Float64 ∘ d::DiscreteDistribution = begin
    if (0 <= p <= 1)
        rand(Binomial(rand(d), p))
    else
        error("Invalid thinning probability.")
    end
end
# Generic thinning with distribution for sum limit
dInner::Distribution ∘ dOuter::DiscreteDistribution = sum(rand(dInner, rand(dOuter)))

# CRN thinning (Binomial)
function ∘(p::Float64, X::Integer, u::Float64)
    if !(0 <= p <= 1)
        error("Invalid thinning probability.")
    end

    if X < 0
        error("X must be a positive integer.")
    end

    if !(0 <= u <= 1)
        error("Invalid CRN.")
    end

    quantile(Binomial(X, p), u)
end

# CRN thinning (generic)
function ∘(d::Distribution, X::Integer, u::Vector{Float64})
    if !all(0 .<= u .<= 1)
        error("Invalid CRN.")
    end

    if X < 0
        error("X must be a positive integer.")
    end

    if length(u) < X
        error("Too few CRNs.")
    end

    sum(quantile.(d, u[1:X]))
end


text:
# E. Dynamic Parallelism

## streams

export CuDeviceStream

struct CuDeviceStream
    handle::cudaStream_t

    function CuDeviceStream(flags=cudaStreamNonBlocking)
        handle_ref = Ref{cudaStream_t}()
        cudaStreamCreateWithFlags(handle_ref, flags)
        return new(handle_ref[])
    end
end

Base.unsafe_convert(::Type{cudaStream_t}, s::CuDeviceStream) = s.handle

function unsafe_destroy!(s::CuDeviceStream)
    cudaStreamDestroy(s)
    return
end

## execution

# device-side counterpart of launch
@inline function device_launch(f, blocks, threads, shmem, stream, args...)
    blockdim = CuDim3(blocks)
    threaddim = CuDim3(threads)

    buf = parameter_buffer(f, blockdim, threaddim, shmem, args...)
    cudaLaunchDeviceV2(buf, stream)

    return
end

@generated function parameter_buffer(f, blocks, threads, shmem, args...)
    # allocate a buffer
    ex = quote
        Base.@_inline_meta
        buf = cudaGetParameterBufferV2(f, blocks, threads, shmem)
        ptr = Base.unsafe_convert(Ptr{UInt32}, buf)
    end

    # store the parameters
    #
    # D.3.2.2. Parameter Buffer Layout
    # > Each individual parameter placed in the parameter buffer is required to be aligned.
    # > That is, each parameter must be placed at the n-th byte in the parameter buffer,
    # > where n is the smallest multiple of the parameter size that is greater than the
    # > offset of the last byte taken by the preceding parameter. The maximum size of the
    # > parameter buffer is 4KB.
    #
    # NOTE: the above seems wrong, and we should use the parameter alignment, not its size.
    last_offset = 0
    for i in 1:length(args)
        T = args[i]
        align = Base.datatype_alignment(T)
        offset = Base.cld(last_offset, align) * align
        push!(ex.args, :(
            Base.pointerset(convert(Ptr{$T}, ptr+$offset), args[$i], 1, $align)
        ))
        last_offset = offset + sizeof(T)
    end

    push!(ex.args, :(return buf))

    return ex
end

## synchronization

"""
    device_synchronize()

Wait for the device to finish. This is the device side version,
and should not be called from the host.

`device_synchronize` acts as a synchronization point for
child grids in the context of dynamic parallelism.
"""
device_synchronize() = cudaDeviceSynchronize()


text:
using Documenter, TangentMps

makedocs(;
    modules=[TangentMps],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ho-oto/TangentMps.jl/blob/{commit}{path}#L{line}",
    sitename="TangentMps.jl",
    authors="Hayate Nakano",
)

deploydocs(;
    repo="github.com/ho-oto/TangentMps.jl",
)