text:
using LinearAlgebra

export fitEllipse

function fitEllipse(x::Array{Float64, 1}, y::Array{Float64, 1})
    # outputs [A; B; C; D; E; F], a Vector containing the parameters of the equation
    # Ax² + Bxy + Cy² + Dx + Ey + F = 0
    #
    # Based on Direct Least Squares Fitting of Ellipses
    # by A. W. Fitzgibbon, M. Pilu, R. B. Fisher
    # http://cseweb.ucsd.edu/~mdailey/Face-Coord/ellipse-specific-fitting.pdf

    # Build design matrix
    D = [x.*x x.*y y.*y x y ones(Float64, size(x))]

    # Build scatter matrix
    S = transpose(D) * D

    # Build 6x6 constraint matrix
    C = zeros(Float64, (6, 6))
    C[1, 3] = C[3, 1] = -2
    C[2, 2] = 1

    # Solve generalised eigensystem
    geval, gevec = eigen(S, C)
    # geval - eigenvalues
    # gevec - eigenvectors

    # Find the only negative eigenvalue index
    negindex = findall(x->x<0, geval)

    # Get fitted parameters
    return a = gevec[:, negindex]

end


text:
#=
# 2D DG poisson
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("dgpoisson3d");

@useLog("dgpoisson3dlog")

@domain(3)
@solver(DG)
@functionSpace(LEGENDRE, 1)

n = 12;
@mesh(HEXMESH, n)

@variable(u)
@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)

@coefficient(f, "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@coefficient(beta, 100000*n)
@weakForm(u, "-dot(grad(u), grad(v)) - f*v + surface(dot(ave(grad(u)), jump(v)) + dot(ave(grad(v)), jump(u)) - beta*dot(jump(u), jump(v)))")

solve(u);

# # solution is stored in the variable's "values"
# using Plots
# pyplot();
# display(plot(Finch.grid_data.allnodes[1,:], Finch.grid_data.allnodes[2,:], u.values[:], st=:surface))

# exact solution is sin(pi*x)*sin(pi*y)*sin(pi*z)
# check error
maxerr = 0;
exact(x,y,z) = sin(pi*x)*sin(pi*y)*sin(pi*z);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    z = Finch.grid_data.allnodes[3,i];
    err = abs(u.values[i] - exact(x,y,z));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

@finalize()

ext: jl

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using BinaryBuilder

# Collection of sources required to build Bonmin and Couenne
name = "Bonmin-Couenne"
version = v"0.0"
sources = [
    "https://github.com/coin-or/Couenne/archive/releases/0.5.6.tar.gz" =>
    "",

    "https://github.com/coin-or/Bonmin/archive/releases/1.8.6.tar.gz" =>
    ""
]

# Bash recipe for building across all platforms
script = raw"""
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:i686),
    Windows(:x86_64),
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),
    MacOS(),
    FreeBSD(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "couenne", :couenne),
    ExecutableProduct(prefix, "bonmin", :bonmin)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jgoldfar/BonminBuilder

max_issues_repo_name: jgoldfar/BonminBuilder

max_forks_repo_name: jgoldfar/BonminBuilder


text:
module Perf
using ThreadingTools
using ThreadingTools: tname
using BenchmarkTools

macro race(f, args...)
    tf = tname(f)
    tt_call   = :(($tf)($(args...)))
    base_call = :((Base.$f)($(args...)))
    call_str = string(:($f($(args...))))
    quote
        println("Benchmark: ", $call_str)
        print("Base: ")
        @btime $(base_call)
        print("TT  : ")
        @btime $(tt_call)
        println("#"^80)
    end |> esc
end

@info "Running benchmarks on $(Threads.nthreads()) threads."
data = randn(10^6)
@race(sum,     sin, data)
@race(prod,    sin, data)
@race(minimum, sin, data)
@race(maximum, sin, data)
@race(mapreduce, sin, +, data)
@race(map, sin, data)

end#module


text:
"""
scipy.spatial module

- [Spatial algorithms and data structures (scipy.spatial) Reference Guide](https://docs.scipy.org/doc/scipy/reference/spatial.html)

# Examples

You can calculate several rotation representations:

```julia-repl
julia> R = SciPy.spatial.transform.Rotation;

julia> r = R.from_quat([0, 0, sin(π/4.0), cos(π/4)]);

julia> r.as_matrix()
3×3 Array{Float64,2}:
 2.22045e-16  -1.0          0.0
 1.0           2.22045e-16  0.0
 0.0           0.0          1.0

julia> r.as_rotvec()
3-element Array{Float64,1}:
 0.0
 0.0
 1.5707963267948963

julia> r.as_euler("zyx", degrees=true)
3-element Array{Float64,1}:
 90.0
  0.0
  0.0
```
"""
module spatial

using PyCall
import PyCall: hasproperty # Base.hasproperty in Julia 1.2

import ..pyspatial
import .._generate_docstring
import ..LazyHelp

const _ignore_funcs = ["spatial"]

# This list can be generated by print_scipy_api_list.py script
all_properties = ["ConvexHull", "Delaunay", "HalfspaceIntersection", "KDTree", "Rectangle", "SphericalVoronoi", "Voronoi", "cKDTree", "ckdtree", "convex_hull_plo
t_2d", "delaunay_plot_2d", "distance", "distance_matrix", "geometric_slerp", "kdtree", "minkowski_distance", "minkowski_distance_p", "procrustes
", "qhull", "tsearch", "voronoi_plot_2d", "distance", "transform"]

for f in all_properties
    f in _ignore_funcs && continue

    sf = Symbol(f)
    @eval @doc LazyHelp(pyspatial, $f) $sf(args...; kws...) = pycall(pyspatial.$f, PyAny, args...; kws...)
end

function __init__()
    copy!(pyspatial, pyimport_conda("scipy.spatial", "scipy"))
end

end # module


text:
function run_simulation(filter::AbstractFilter, x0::AbstractVector, action_seq::AbstractArray)
    """ Simulate N=length(action_seq) points correspding to a
    Dynamic/Observation Model encapsulated in a KF struct.
    Start from an initial state belief s0"""
    meas = Vector{typeof(x0)}()
    states = [x0]
    for u in action_seq
        xp = dynamic(filter, states[end], u)
        y = observation(filter, xp)
        push!(meas, y)
        push!(states, xp)
    end
    return states, meas
end

function run_filter(filter::AbstractFilter, s0::State, action_history::AbstractArray,
    measurement_history::AbstractArray)
    """
    Run Filter on a measurement_history points, for a given action_history
    Start from an initial state belief s0
    """
    @assert length(action_history) == length(measurement_history)
    states = [s0]
    for (u, y) in zip(action_history, measurement_history)
        # filter.A += rand(Normal(0, 0.005), (3, 3))
        sp = prediction(filter, states[end], u)
        sn = correction(filter, sp, y)
        push!(states, sn)
    end
    return states
end

function run_param_filter(θ, param_filter::AbstractParamFilter, s0::State, action_history::AbstractArray,
    measurement_history::AbstractArray)
    """
    Run parametrized Klaman Filter on a measurement_history points, for a given action_history
    Start from an initial state belief s0
    """
    @assert length(action_history) == length(measurement_history)
    states = [s0]
    for (u, y) in zip(action_history, measurement_history)
        sp = param_prediction(θ, param_filter, states[end], u)
        sn = param_correction(θ, param_filter, sp, y)
        push!(states, sn)
    end
    return states
end

function run_kf_gradient(θ, param_kf::ParamKalmanFilter, s0::State, action_history::AbstractArray, measurement_history::AbstractArray,
    opt, epochs)
    """
    Run gradient descent on unknown parameters of a linear on non linear state space model (parametrized abstract filter).
    returns the found paramters after n epochs, and the associated loss function
    """
    @assert length(action_history)==length(measurement_history)
    # Compute initial state esimates
    loss = []
    for i in ProgressBar(1:epochs)
        states = run_param_filter(θ, param_kf, s0, action_history, measurement_history)
        ∇, = gradient(theta -> kf_likelihood(theta, param_kf, states, action_history, measurement_history), θ)
        update!(opt, θ, ∇)
        l = kf_likelihood(θ, param_kf, states, action_history, measurement_history)
        push!(loss, l)
    end
    return θ, loss
end

function run_ekf_gradient(θ, param_ekf::ExtendedParamKalmanFilter, s0::State, action_history::AbstractArray, measurement_history::AbstractArray,
    opt, epochs)
    """
    Run gradient descent on unknown parameters of a linear on non linear state space model (parametrized abstract filter).
    returns the found paramters after n epochs, and the associated loss function
    """
    @assert length(action_history)==length(measurement_history)
    # Compute initial state esimates
    loss = []
    for i in ProgressBar(1:epochs)
        states = run_param_filter(θ, param_ekf, s0, action_history, measurement_history)
        ∇, = gradient(theta -> ekf_likelihood(theta, param_ekf, states, action_history, measurement_history), θ)
        update!(opt, θ, ∇)
        l = ekf_likelihood(θ, param_ekf, states, action_history, measurement_history)
        push!(loss, l)
    end
    return θ, loss
end


text:
"""
    isweekend(dt)

Returns `true` for Saturdays or Sundays.
Returns `false` otherwise.
"""
function isweekend(dt::Date)
    return dayofweek(dt) in [6, 7]
end

"""
    isbday(calendar, dt)

Returns `true` for weekends or holidays.
Returns `false` otherwise.
"""
function isbday(hc::HolidayCalendar, dt::Date)
    if _getcachestate(hc)
        hcc :: HolidayCalendarCache = _getholidaycalendarcache(hc)
        return isbday(hcc, dt)
    else
        return !(isweekend(dt) || isholiday(hc, dt))
    end
end

isbday(calendar, dt) = isbday(convert(HolidayCalendar, calendar), dt)

"""
    tobday(calendar, dt; [forward=true])

Adjusts `dt` to next Business Day if it's not a Business Day.
If `isbday(dt)`, returns `dt`.
"""
function tobday(hc::HolidayCalendar, dt::Date; forward::Bool = true)
    if isbday(hc, dt)
        return dt
    else
        local next::Date
        local increment::Int = forward ? 1 : -1
        next = dt + Dates.Day(increment)

        while !isbday(hc, next)
            next += Dates.Day(increment)
        end
    end

    return next
end

tobday(calendar, dt; forward::Bool = true) = tobday(convert(HolidayCalendar,calendar), dt; forward=forward)

"""
    advancebdays(calendar, dt, bdays_count)

Increments given date `dt` by `bdays_count`.
Decrements it if `bdays_count` is negative.
`bdays_count` can be a `Int`, `Vector{Int}` or a `UnitRange`.

Computation starts by next Business Day if `dt` is not a Business Day.
"""
function advancebdays(hc::HolidayCalendar, dt::Date, bdays_count::Int)
    local result::Date = tobday(hc, dt)

    # does nothing
    if bdays_count == 0
        return result
    end

    # if bdays_count is positive, goes forward. Otherwise, goes backwards.
    local increment::Int = bdays_count > 0 ? +1 : -1

    local num_iterations::Int = abs(bdays_count)

    while num_iterations > 0
        result += Dates.Day(increment)

        # Looks for previous / next Business Day
        while !isbday(hc, result)
            result += Dates.Day(increment)
        end

        num_iterations += -1
    end

    return result
end

advancebdays(calendar, dt, bdays_count) = advancebdays(convert(HolidayCalendar,calendar), dt, bdays_count)

"""
    bdays(calendar, dt0, dt1)

Counts the number of Business Days between `dt0` and `dt1`.
Returns instances of `Dates.Day`.

Computation is always based on next Business Day if given dates are not Business Days.
"""
function bdays(hc::HolidayCalendar, dt0::Date, dt1::Date)
    if _getcachestate(hc)
        hcc::HolidayCalendarCache = _getholidaycalendarcache(hc)
        return bdays(hcc, dt0, dt1)
    else
        dt0 = tobday(hc, dt0)
        dt1 = tobday(hc, dt1)
        inc::Int = dt0 <= dt1 ? +1 : -1

        local result::Int = 0
        while dt0 != dt1
            dt0 += Dates.Day(inc)

            # Looks for next/last Business Day
            while !isbday(hc, dt0)
                dt0 += Dates.Day(inc)
            end

            result += inc
        end

        return Dates.Day(result)
    end
end

bdays(calendar, dt0, dt1) = bdays(convert(HolidayCalendar, calendar), dt0, dt1)


text:
# Gaussian mixture model
struct GMM{T<:AbstractVector{<:Real}, M<:Normal, S<:InverseGamma, W<:Dirichlet} <: MCMC.Model
  y::T
  K::Int
  mu::M
  sigmasq::S
  eta::W
end
function GMM(y::AbstractVector{<:Real}, K::Int; mu=Normal(0, 1),
             sigmasq=InverseGamma(3, 2), eta=Dirichlet(K, 1/K))
  return GMM(y, K, mu, sigmasq, eta)
end
function MCMC.make_init_state(m::GMM)
  nobs = length(m.y)
  eta = rand(m.eta)
  return (mu=rand(m.mu, m.K), sigmasq=rand(m.sigmasq, m.K),
          eta=eta, lambda=rand(Categorical(eta), nobs))
end
function MCMC.update(m::GMM, s::T) where T
  s = setproperty!!(s, :mu, update_mu(m, s))
  s = setproperty!!(s, :sigmasq, update_sigmasq(m, s))
  s = setproperty!!(s, :eta, update_eta(m, s))
  s = setproperty!!(s, :lambda, update_lambda(m, s))

  return s
end
function update_mu(m::GMM, s::T) where T
  m_mu, s_mu = params(m.mu)

  vkernels = zero.(s.mu)
  mkernels = zero.(s.mu)

  for n in eachindex(m.y)
    k = s.lambda[n]
    vkernels[k] += 1
    mkernels[k] += m.y[n]
  end

  vnew = 1 ./ (s_mu^-2 .+ vkernels ./ s.sigmasq)
  mnew = vnew .* (m_mu/(s_mu^2) .+ mkernels ./ s.sigmasq)

  return randn(m.K) .* sqrt.(vnew) + mnew
end
function update_sigmasq(m::GMM, s::T) where T
  a, b = params(m.sigmasq)
  akernel = zero.(s.sigmasq)
  bkernel = zero.(s.sigmasq)

  for n in eachindex(m.y)
    k = s.lambda[n]
    akernel[k] += 1
    bkernel[k] += (m.y[n] - s.mu[k]) ^ 2
  end

  anew = a .+ akernel / 2
  bnew = b .+ bkernel / 2
  s.sigmasq .= rand.(InverseGamma.(anew, bnew))

  return s.sigmasq
end
function update_eta(m::GMM, s::T) where T
  anew = copy(collect(m.eta.alpha))
  for lam in s.lambda
    anew[lam] += 1
  end
  return rand(Dirichlet(anew))
end
function update_lambda(m::GMM, s::T) where T
  logeta = s.eta
  lambda = [let
              logmix = normlogpdf.(s.mu, sqrt.(s.sigmasq), yn) + logeta
              wsample_logprob(logmix)
            end for yn in m.y]
  return lambda
end
@testset "GMM" begin
  true_model = MixtureModel([Normal(3, .6), Normal(-2, 1)], [.3, .7])

  Random.seed!(0)

  nobs = 1000
  y = rand(true_model, nobs)
  model = GMM(y, 2)
  nsamps = 1000
  chain, _ = mcmc(model, nsamps, nburn=1000)

  mu = hcat(map(c -> c.mu, chain)...)
  sigma = hcat(map(c -> sqrt.(c.sigmasq), chain)...)
  eta = hcat(map(c -> c.eta, chain)...)
  lambda = hcat(map(c -> c.lambda, chain)...)

  @test size(mu) == (model.K, nsamps)
  @test size(sigma) == (model.K, nsamps)
  @test size(eta) == (model.K, nsamps)
  @test size(lambda) == (nobs, nsamps)

  function match_components(postmean::AbstractVector{<:Real}, truth::AbstractVector{<:Real}, atol)
    cond1 = abs2(postmean[1] - truth[1]) < atol
    cond2 = abs2(postmean[2] - truth[2]) < atol
    cond3 = abs2(postmean[1] - truth[2]) < atol
    cond4 = abs2(postmean[2] - truth[1]) < atol
    return (cond1 && cond2) || (cond3 && cond4)
  end
  function match_components(postmean, sym::Symbol, atol)
    truth = [getfield(true_model.components[1], sym),
             getfield(true_model.components[2], sym)]
    return match_components(postmean, truth, atol)
  end
  mean_mu = vec(mean(mu, dims=2))
  mean_sigma = vec(mean(sigma, dims=2))
  mean_eta = vec(mean(eta, dims=2))
  # println(mean_mu)
  # println(mean_sigma)
  # println(mean_eta)
  @test match_components(mean_mu, :μ, 1e-1)
  @test match_components(mean_sigma, :σ, 1e-2)
  @test match_components(mean_eta, true_model.prior.p, 1e-2)
end


text:
# ambiguity.jl
# tests needed only because of code added to resolve method ambiguities

using LinearAlgebra: Diagonal, UpperTriangular
using LinearAlgebra: Adjoint, Transpose, Symmetric
using LinearAlgebra: TransposeAbsVec, AdjointAbsVec
using LinearAlgebra: givens
import LinearAlgebra
using LinearMapsAA
using Test: @test, @testset, @test_throws

M = rand(2,2)
A = LinearMapAA(M)

@testset "Diagonal" begin
    D = Diagonal(1:2)
    @test Matrix(A * D) == M * D
    @test Matrix(D * A) == D * M
end

@testset "UpperTri" begin
    U = UpperTriangular(M)
    @test Matrix(A * U) == M * U
    @test Matrix(U * A) == U * M
end

@testset "TransposeVec" begin
    xt = Transpose(1:2)
    @test xt isa TransposeAbsVec
    @test xt * M == xt * A
end

@testset "AdjointVec" begin
    r = Adjoint(1:2)
    @test r isa AdjointAbsVec
    @test r * M == r * A
end

@testset "Transpose" begin
# work-around per
# https://github.com/Jutho/LinearMaps.jl/issues/147
    LinearAlgebra.isposdef(A::Transpose) = LinearAlgebra.isposdef(parent(A))
    S = LinearAlgebra.Symmetric(M)
    T = LinearAlgebra.Transpose(S)
    @test Matrix(A * T) == M * T # failed prior to isposdef overload
    @test Matrix(T * A) == T * M
end

@testset "Adjoint" begin
    C = rand(ComplexF32, 2, 2)
    H = LinearAlgebra.Hermitian(C)
    J = Adjoint(H)
    LinearAlgebra.isposdef(A::Adjoint) = LinearAlgebra.isposdef(parent(A))
    @test Matrix(A * J) == M * J
    @test Matrix(J * A) == J * M
end

@testset "AbsRot" begin
    G, _ = givens(1., 2., 3, 4)
    R = LinearAlgebra.Rotation([G, G])
    AR = adjoint(R); # semicolon needed, otherwise "show" error!
    # size(AR) # fails because AR is size-less and not an AbstractMatrix
    AR isa Adjoint{<:Any,<:LinearAlgebra.AbstractRotation} # true
    # ones(4,4) * AR # works
    @test_throws String A * AR
end


text:
#=Copyright (c) 2017 Gabriel Goh

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.=#

csize = 5;
y = [-1*ones(csize); 1*ones(csize); -1*ones(csize);1*ones(csize);-1*ones(csize);]
y = y + 0.3*cat(1, [[0,1,0,1,0] for i in 1:5]...)
x = cat(1,[collect(c + linspace(-0.03,0.03, csize)) for c in linspace(-1,1,5)]...)

pdegree = 25

# Construct Vander* Matrix
A = zeros(length(y),pdegree+1)
for i = 1:size(A,1), j = 1:size(A,2)
	A[i,j] = x[i]^(j-1)
end

Σ,U = eig(A'A)

p(w,x) = sum([w[i]*x^(i-1) for i = 1:length(w)])

Q,R = qr(A)
ws = (R'R)\(A'*y)
plot(r,[p(ws,i) for i in r])
plot!(ylim = [-2,2])
scatter!(x,y)

r = linspace(x[1]-0.1,1+0.1,500)

w = 100*randn(size(w))
α = 1/max(eig(A'A)[1]...)

# for ind = 1:2
#   plot!(r,[(U'*ws)[ind]*p(U[:,ind],i) for i in r])
# end
# plot!(ylim = [0,15])
# scatter!(x,y)

@gif for i = 1:100
  println(i)
  w = getiterk(A'A, A'y, α)(1.4^i)
  println(norm(A'*(A*w - y)))
  plot(size = (880,180))
  scatter!(x,y)
  plot!(r,[p(w,i) for i in r])
  plot!(ylim = [-2,2])
end

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
rationalize(0.9054054)
rationalize(0.518518)
rationalize(0.75)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module JungleHelperCockatiel

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Cockatiel" Cockatiel(x::Integer, y::Integer, facingLeft::Bool=true, sprite::String="")

const placements = Ahorn.PlacementDict(
    "Cockatiel (Jungle Helper)" => Ahorn.EntityPlacement(
        Cockatiel
    )
)

function Ahorn.selection(entity::Cockatiel)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 6, y - 4, 12, 12)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Cockatiel, room::Maple.Room)
    scaleX = (get(entity.data, "facingLeft", false) ? 1 : -1)
    Ahorn.drawSprite(ctx, "JungleHelper/Cockatiel/idle00", 0, 2, sx = scaleX)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
# This file contains the tests for the `CSV` module

module TestCSV

using PDFHighlights
using Suppressor
using SyntaxTree
using Test

# Print the header
println("\e[1;32mRUNNING\e[0m: TestCSV.jl")

const pdf = joinpath(@__DIR__, "..", "pdf", "TestPDF.pdf")
const HEADER = PDFHighlights.Internal.Constants.HEADER

macro initialize()
    return esc(
        quote
            csv = mktemp((path, _) -> string(path, ".csv"))
            initialize(csv)
        end
    )
end

@testset "initialize" begin

    @initialize

    @test readlines(csv) == String[HEADER,]
    @test_nowarn initialize(csv)

    csv = mktemp((path, _) -> string(path, ".csv"))
    touch(csv)
    initialize(csv)

    @test readlines(csv) == String[HEADER,]
    @test_nowarn initialize(csv)

    @test_throws PDFHighlights.Internal.Exceptions.NotCSV("oof") initialize("oof")

end

@testset "_check" begin

    @initialize

    @test_nowarn PDFHighlights.Internal.CSV._check(csv)

    line = "\"Highlight with \"\"quotes\"\"\",\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test_nowarn PDFHighlights.Internal.CSV._check(csv)

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",\"Note 1\",\"Love\""

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test_throws(
        PDFHighlights.Internal.Exceptions.LastElementIsNotAnInteger(csv, 2),
        PDFHighlights.Internal.CSV._check(csv),
    )

    line = "Wrong header"

    open(csv, "w") do io
        println(io, line)
    end

    @test_throws(
        PDFHighlights.Internal.Exceptions.IncorrectHeader(csv),
        PDFHighlights.Internal.CSV._check(csv),
    )

    line = "Column 1, Column 2, Column 3, Column 4, Column 5, 6"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test_throws(
        PDFHighlights.Internal.Exceptions.NotFiveColumns(csv, 2),
        PDFHighlights.Internal.CSV._check(csv),
    )

end

@testset "import_highlights" begin

    @initialize

    @suppress_out @test_nowarn import_highlights(csv, pdf)
    @suppress_out @test_nowarn import_highlights(csv, pdf)

    @initialize

    @suppress_out @test_nowarn import_highlights(csv, pdf; quiet=true)
    @suppress_out @test_nowarn import_highlights(csv, pdf; quiet=true)

    @initialize

    s1 = @capture_out import_highlights(csv, pdf)
    s2 = """

        CSV: "$(basename(csv))"
        PDF: "$(basename(pdf))"
        Highlights (found / added): 7 / 7

    """

    @test s1 == s2

    s1 = @capture_out import_highlights(csv, pdf)
    s2 = """

        CSV: "$(basename(csv))"
        PDF: "$(basename(pdf))"
        Highlights (found / added): 7 / 0

    """

    @test s1 == s2

    @test @capture_out(import_highlights(csv, pdf; quiet=true)) == ""

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        import_highlights("oof", pdf),
    )

    @initialize

    dir = joinpath(@__DIR__, "..")

    s1 = @capture_out import_highlights(csv, dir)
    s2 = """

        CSV: "$(basename(csv))"
        Directory: "$(basename(dirname(normpath(dir))))"
        Highlights (found / added): 7 / 7

    """

    @test s1 == s2

    s1 = @capture_out import_highlights(csv, dir)
    s2 = """

        CSV: "$(basename(csv))"
        Directory: "$(basename(dirname(normpath(dir))))"
        Highlights (found / added): 7 / 0

    """

    @test s1 == s2

    @initialize

    dir = dirname(normpath(@__DIR__, ".."))

    s1 = @capture_out import_highlights(csv, dir)
    s2 = """

        CSV: "$(basename(csv))"
        Directory: "$(basename(dir))"
        Highlights (found / added): 7 / 7

    """

    @test s1 == s2

    s1 = @capture_out import_highlights(csv, dir)
    s2 = """

        CSV: "$(basename(csv))"
        Directory: "$(basename(dir))"
        Highlights (found / added): 7 / 0

    """

    @test s1 == s2

end

@testset "_get_authors_from_CSV" begin

    @initialize

    @test PDFHighlights.Internal.CSV._get_authors_from_CSV(csv) == String[]

    line = "\"Highlight 1\",\"Title 1\",\"Author with \"\"quotes\"\"\",\"Note 1\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test PDFHighlights.Internal.CSV._get_authors_from_CSV(csv) ==
    String["Author with \"quotes\""]

    line = "\"Highlight 1\",\"Title 1\",Author without quotes,\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_authors_from_CSV(csv) ==
    String["Author without quotes"]

    line = "\"Highlight 1\",\"Title 1\",\"\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_authors_from_CSV(csv) ==
    String[""]

    line = "\"Highlight 1\",\"Title 1\",,\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_authors_from_CSV(csv) ==
    String[""]

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        PDFHighlights.Internal.CSV._get_authors_from_CSV("oof"),
    )

end

@testset "_get_highlights_from_CSV" begin

    @initialize

    @test PDFHighlights.Internal.CSV._get_highlights_from_CSV(csv) == String[]

    line = "\"Highlight with \"\"quotes\"\"\",\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test PDFHighlights.Internal.CSV._get_highlights_from_CSV(csv) ==
    String["Highlight with \"quotes\""]

    line = "Highlight without quotes,\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_highlights_from_CSV(csv) ==
    String["Highlight without quotes"]

    line = "\"\",\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_highlights_from_CSV(csv) ==
    String[""]

    line = ",\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_highlights_from_CSV(csv) ==
    String[""]

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        PDFHighlights.Internal.CSV._get_highlights_from_CSV("oof"),
    )

end

@testset "_get_titles_from_CSV" begin

    @initialize

    @test PDFHighlights.Internal.CSV._get_titles_from_CSV(csv) == String[]

    line = "\"Highlight 1\",\"Title with \"\"quotes\"\"\",\"Author 1\",\"Note 1\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test PDFHighlights.Internal.CSV._get_titles_from_CSV(csv) ==
    String["Title with \"quotes\""]

    line = "\"Highlight 1\",Title without quotes,\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_titles_from_CSV(csv) ==
    String["Title without quotes"]

    line = "\"Highlight 1\",\"\",\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_titles_from_CSV(csv) ==
    String[""]

    line = "\"Highlight 1\",,\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test PDFHighlights.Internal.CSV._get_titles_from_CSV(csv) ==
    String[""]

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        PDFHighlights.Internal.CSV._get_titles_from_CSV("oof"),
    )

end

@testset "@extract" begin

    @test @macroexpand(PDFHighlights.Internal.CSV.@extract(highlights)) |> linefilter! ==
    quote
        if current_comma_index == 1
            highlights[line_index] = ""
        else
            piece = line[1:(current_comma_index - 1)]
            if piece == "\"\""
                highlights[line_index] = ""
            elseif startswith(piece, "\"")
                highlights[line_index] = chop(piece; head = 1, tail = 1)
            else
                highlights[line_index] = piece
            end
        end
    end |> linefilter!

    @test @macroexpand(PDFHighlights.Internal.CSV.@extract(locations)) |> linefilter! ==
    quote
        if current_comma_index == lastindex(line)
            locations[line_index] = 0
        else
            locations[line_index] = parse(
                Int32,
                line[(current_comma_index + 1):end],
            )
        end
    end |> linefilter!

    array = :titles

    @test @macroexpand(PDFHighlights.Internal.CSV.@extract(titles)) |> linefilter! ==
    quote
        if current_comma_index == previous_comma_index + 1
            $(array)[line_index] = ""
        else
            piece = line[(previous_comma_index + 1):(current_comma_index - 1)]
            if piece == "\"\""
                $(array)[line_index] = ""
            elseif startswith(piece, "\"")
                $(array)[line_index] = chop(piece; head = 1, tail = 1)
            else
                $(array)[line_index] = piece
            end
        end
    end |> linefilter!

end

@testset "get_all" begin

    @initialize

    @test get_all(csv) == (String[], String[], String[], String[], Int32[])

    line = "\"Highlight 1\",\"Title with \"\"quotes\"\"\",\"Author 1\",\"Note 1\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test get_all(csv) == (
        String["Highlight 1"],
        String["Title with \"quotes\""],
        String["Author 1"],
        String["Note 1"],
        Int32[1],
    )

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        get_all("oof"),
    )

end

@testset "get_locations" begin

    @initialize

    @test get_locations(csv) == Int32[]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",\"Note 1\","

    open(csv, "a") do io
        println(io, line)
    end

    @test get_locations(csv) ==
    Int32[0]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",\"Note 1\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test get_locations(csv) ==
    Int32[1]

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        get_locations("oof"),
    )

end

@testset "get_notes" begin

    @initialize

    @test get_notes(csv) == String[]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",\"Note with \"\"quotes\"\"\",1"

    open(csv, "a") do io
        println(io, line)
    end

    @test get_notes(csv) ==
    String["Note with \"quotes\""]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",Note without quotes,1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test get_notes(csv) ==
    String["Note without quotes"]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",\"\",1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test get_notes(csv) ==
    String[""]

    line = "\"Highlight 1\",\"Title 1\",\"Author 1\",,1"

    open(csv, "w") do io
        println(io, HEADER, '\n', line)
    end

    @test get_notes(csv) ==
    String[""]

    @test_throws(
        PDFHighlights.Internal.Exceptions.IntegrityCheckFailed("oof"),
        get_notes("oof"),
    )

end

end


text:
using OPFLearn
using Test

# Distributed tests will be performed if a string of an integer > 3, representing the number
# of processors to use is passed to the package test function as follows,
# Pkg.test("OPFLearn"; test_args = ["4"])
if !isempty(ARGS)
	try
		global nproc = parse(Int, ARGS[1])
	catch
		global nproc = nothing
	end
else
	global nproc = nothing
end

@testset "OPFLearn.jl" begin
	include("sample.jl")
	include("create_samples.jl")
	include("io.jl")
end

if !isnothing(nproc) && nproc > 2
	@testset "Distributed" begin
		include("distributed.jl")
	end
end


text:
using HPAT
using ParallelAccelerator
ParallelAccelerator.ParallelIR.PIRHoistParfors(1)
ParallelAccelerator.CGen.setRecreateLoops(true)
ParallelAccelerator.CGen.setRecreateConds(true)

@acc hpat function kernelscore_test(n)
    X = rand(n)
    points = [-1.0, 2.0, 5.0]
    N = size(points,1)
    b = 0.5
    exps::Float64 = 0.0
    @par exps(+) for i in 1:length(X)
        d = -(X[i]-points).^2./(2*b^2)
        m = minimum(d)
        exps += m-log(b*N)+log(sum(exp(d-m)))
    end
    return exps
end

println(kernelscore_test(128))

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
abstract type ShepardType <: InterpolationMethod end

export Shepard

"""
    Shepard(P = 2)

Standard Shepard interpolation with power parameter `P`.
"""
struct Shepard{T} <: ShepardType where T <: Real
    P::T
end
Shepard() = Shepard(2)

struct ShepardInterpolant{T1, T2, F, M} <: ScatteredInterpolant where {T1 <: AbstractArray, T2 <: AbstractMatrix{<:Real}}

    data::T1
    points::T2
    idw::F
    metric::M
end

# No need to compute anything here, everything is done in the evaluation step.
function interpolate(idw::ShepardType,
                     points::AbstractArray{<:Real,2},
                     samples::AbstractArray{<:Number,N};
                     metric = Euclidean()) where {N}

    return ShepardInterpolant(samples, points, idw, metric)
end

function evaluate(itp::ShepardInterpolant, points::AbstractArray{<:Real,2})

    # Compute distances between sample points and interpolation points
    d = pairwise(itp.metric, itp.points, points;dims=2)

    # Evaluate point by point
    m = size(points, 2)
    n = size(itp.data, 2)
    values = zeros(eltype(itp.data), m, n)
    for i = 1:m

        d_col = d[:,i]

        # If an interpolation point coincide with a sampling point, just return the
        # original data. Otherwise, compute distance-weighted sum
        if !all(r > 0 for r in d_col)
            ind = findfirst(x -> x ≈ 0.0, d_col)
            values[i,:] = itp.data[ind, :]
        else
            values[i,:] = evaluatePoint(itp.idw, itp.points, itp.data, d_col)
        end
    end

    return values
end

# Original Shepard
function evaluatePoint(idw::Shepard,
                       dataPoints::AbstractArray{<:Real,2},
                       data::AbstractArray{<:Number,N},
                       d::AbstractVector) where {N}

    # Compute weigths and return the weighted sum
    w = 1.0./(d.^idw.P)
    value = sum(w.*data, dims = 1)./sum(w)
end


text:
function testdump()
    println("\n in testdump")
    dim = 10
    hnsw = createHnswApi(Float32, 8, 16, "DistL1")
    # block // insertion
    nbinsert = 500
    parallel = true
    if parallel
        datas = rand(Float32, (dim, nbinsert));
        data_insert = map(i -> ( rand(Float32, dim) , UInt(i) ), 1:size(datas)[2])
        insert(hnsw, data_insert)
    else
        for i in 1:nbinsert
            data = rand(Float32, dim);
            insert(hnsw, data, UInt64(i))
        end
    end
    #
    v1 = [Float32(0.5) for i in 1:dim]
    neighbours = search(hnsw, v1, 8, 16)
    @info " 0.5 vector neighbours list : " neighbours
    #
    res = fileDump(hnsw, "testdumpfromjulia")
    if res > 0
        true
    else
        @warn "file dump failed : " "testdumpfromjulia"
        false
    end
end

function testreload()
    println("\n in testreload")
    # we know type and distance we reload from so we do not need to call
    # description = getDescription("testdumpfromjulia")
    res = loadHnsw("testdumpfromjulia", Float32, "DistL1")
    if res != nothing
        println("loadHnsw from file succeeded")
        description , hnswapi = res
    else
        return false
    end
    dim = description.data_dimension
    v1 = [Float32(0.5) for i in 1:dim]
    neighbours = search(hnswapi, v1, 8, 16)
    @info " 0.5 vector neighbours list : " neighbours
    true
end


text:
export linpath
export momentumpath

"""
    momentumpath

Generate a list of momenta

# Arguments
* `anchorpoints`
* (Optional) `nseg` - number of points in each segment
"""
function linpath(anchors::AbstractVector...; nseg::Integer=100)
    n = length(anchors)
    if n <= 1
        throw(ArgumentError("Number of anchor points should be more than one"))
    end
    d = length(anchors[1])
    for i in 2:n
        if length(anchors[i]) != d
            throw(DimensionMismatch("All anchor points should have the same dimension"))
        end
    end
    out = Matrix{Float64}(undef, (d, (n-1) * nseg + 1))
    for i in 1:n-1
        from = anchors[i]
        to = anchors[i+1]
        dvec = (to - from) ./ nseg
        for j in 0:(nseg-1)
            out[:, (i-1)*nseg + j + 1] = from + dvec * j
        end
    end
    out[:, end] = anchors[end]
    return out
end

"""
    momentumpath

The anchorpoints are given in units of the reciprocal lattice vectors.
"""
function momentumpath(
    unitcell::UnitCell,
    anchorpoints::AbstractVector{<:AbstractVector{<:Number}},
)
    real_anchorpoints = [unitcell.reciprocallatticevectors * ap for ap in anchorpoints]
    return linpath(real_anchorpoints...)
end


text:
include("deHoop.jl");

# Precompute some guass quadrature coefficieints
function precompute_quad!(mesh,gauss_order)
    # define Gauss-Legendre quadrature points and weights
    local pt,wg;
    if (gauss_order==3)
	# points
	pt = [2/3 1/6 1/6;
	      1/6 2/3 1/6;
	      1/6 1/6 2/3];
	# weights
        wg = [1/3;
              1/3;
              1/3];
    else
        error("The Gauss quadrature order you specified is not implemented")
    end
    # pre-compute guassP and guassW for all triangles
    for tri = mesh
	VK = [tri.P1'; tri.P2'; tri.P3'];
	tri.aux.gaussP = transpose(pt * VK);
	tri.aux.gaussW = tri.area * wg;
    end
end

function assemble_BEM(lambda, mesh, gauss_order,der=0)

    # Number of triangles = size of matrix
    n = size(mesh,1);

    # Storage for the computed matrix
    T = zeros(ComplexF64,n, n);

    # prepare auxiliary vectors
    rowind = vec(repeat(1:gauss_order,1,gauss_order)'); # [1 1 1 2 2 2 3 3 3 ]
    colind = vec(repeat(1:gauss_order,1,gauss_order))  # [1 2 3 1 2 3 1 2 3]

    for row = 1:n
	@inbounds for col = row:n # Skip because of symmetry

	    # compute (row, col)-entry of T
            trir=mesh[row];
            tric=mesh[col]

            # Compute a distance vector where zero elements
            # are set to one and marked in the idx vector
            S1=trir.aux.gaussP[:, rowind];
            S2=tric.aux.gaussP[:, colind];
            s0=(S1 - S2).^2
	    dist = sqrt.(sum(s0, dims=1));
	    idx = (dist .== 0);
	    dist[idx] .= 1;

            # E = derivative of greens function
            local E
            if (der==0) # Zero'th derivative
	        E = exp.((1im * lambda) .* dist) .- 1;
	        E[vec(idx)] .= 1im * lambda; # Modify in dist-zero points
            elseif (der==1)
	        E = (1im*dist).*(exp.((1im * lambda) .* dist));
	        E[vec(idx)] .= 1im;
            else
	        E = ((1im*dist).^der).*(exp.((1im * lambda) .* dist));
	        E[vec(idx)] .= 0;
            end

            # Compute the matrix element
            mgrW=trir.aux.gaussW[rowind]
            mgcW=tric.aux.gaussW[colind]
            mgrW2=trir.aux.gaussW
            local aa=mgrW .* (mgcW ./ dist');
            T[row,col]=(E * aa)[1] / (4 * pi);
            if (der==0) # Note: this constant term disappears when we differentiate
                sing_kern=deHoop(trir.aux.gaussP, tric); # Compute singular integral
	        T[row, col] += sing_kern' * mgrW2 / (4 * pi);
            end

	    # Set element on other side of diagonal
	    T[col, row] = T[row, col];
	end
    end

    return T

end


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

meshgrid(v::AbstractVector) = meshgrid(v, v)

"""
This function performs the same operation as MATLAB or numpy's meshgrid()
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repeat(vx, m, 1), repeat(vy, 1, n))
end

"""
This function performs the same operation as MATLAB or numpy's meshgrid()
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}, vz::AbstractVector{T}) where T
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, 1, n, 1)
    vy = reshape(vy, m, 1, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end


text:
#!/usr/bin/env julia
using Pkg
push!(LOAD_PATH,"../src/Retriever.jl")
include("../src/Retriever.jl")

Pkg.add("SQLite")
Pkg.add("MySQL")

using Test

# Run tests

@time include("test_retriever.jl")


text:
import Tulip

"""
    QuantileRegression

Quantile regression representation

## Fields

* `τ`: the quantile value
* `X`: the model matrix
* `β`: the coefficients
* `y`: the response vector
* `wts`: the weights
* `wrkres`: the working residuals
* `fitdispersion`: if true, the dispersion is estimated otherwise it is kept fixed
* `fitted`: if true, the model was already fitted
"""
mutable struct QuantileRegression{
    T<:AbstractFloat,
    M<:AbstractMatrix{T},
    V<:AbstractVector{T},
} <: AbstractRobustModel{T}
    τ::T
    X::M
    β::V
    y::V
    wts::V
    wrkres::V
    fitdispersion::Bool
    fitted::Bool

    function QuantileRegression(
        τ::T,
        X::M,
        y::V,
        wts::V,
    ) where {V<:AbstractVector{T},M<:AbstractMatrix{T}} where {T<:AbstractFloat}
        n = length(y)
        m, p = size(X)
        ll = length(wts)
        ll == 0 || ll == n || error("length of wts is $ll, must be $n or 0")
        m == n || error("X has $m rows, must be like y, $n")
        new{T,M,V}(τ, X, zeros(T, p), y, wts, similar(y), false, false)
    end
end

function show(io::IO, obj::QuantileRegression)
    println(
        io,
        "Quantile regression for quantile: $(obj.τ)\n\nCoefficients:\n",
        coeftable(obj),
    )
end

function show(io::IO, obj::TableRegressionModel{M,T}) where {T,M<:QuantileRegression}
    println(
        io,
        "Quantile regression for quantile: $(obj.model.τ)\n\n$(obj.mf.f)\n\nCoefficients:\n",
        coeftable(obj),
    )
end

"""
    quantreg(X, y, args...; kwargs...)

An alias for `fit(QuantileRegression, X, y; kwargs...)`.

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and a `DataFrame`.
"""
quantreg(X, y, args...; kwargs...) = fit(QuantileRegression, X, y, args...; kwargs...)

"""
    fit(::Type{M}, X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
        y::AbstractVector{T}; quantile::AbstractFloat=0.5,
        dofit::Bool          = true,
        wts::FPVector        = similar(y, 0),
        fitdispersion::Bool  = false,
        fitargs...) where {M<:QuantileRegression, T<:AbstractFloat}

Fit a quantile regression model with the model matrix (or formula) X and response vector (or dataframe) y.

It is solved using the exact interior method.

# Arguments

- `X`: the model matrix (it can be dense or sparse) or a formula
- `y`: the response vector or a dataframe.

# Keywords

- `quantile::AbstractFloat=0.5`: the quantile value for the regression, between 0 and 1.
- `dofit::Bool = true`: if `false`, return the model object without fitting;
- `wts::Vector = []`: a weight vector, should be empty if no weights are used;
- `fitdispersion::Bool = false`: reevaluate the dispersion;
- `fitargs...`: other keyword arguments like `verbose` to print iteration details.

# Output

the RobustLinearModel object.

"""
function fit(
    ::Type{M},
    X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
    y::AbstractVector{T};
    quantile::AbstractFloat=0.5,
    dofit::Bool=true,
    wts::FPVector=similar(y, 0),
    fitdispersion::Bool=false,
    fitargs...,
) where {M<:QuantileRegression,T<:AbstractFloat}

    # Check that X and y have the same number of observations
    if size(X, 1) != size(y, 1)
        mess = "number of rows in X and y must match: $(size(X, 1)) != $(size(y, 1))"
        throw(DimensionMismatch(mess))
    end

    # Check quantile is an allowed value
    (0 < quantile < 1) ||
        error("quantile should be a number between 0 and 1 excluded: $(quantile)")

    m = QuantileRegression(quantile, X, y, wts)
    return dofit ? fit!(m; fitargs...) : m
end

function fit(
    ::Type{M},
    X::Union{AbstractMatrix,SparseMatrixCSC},
    y::AbstractVector;
    kwargs...,
) where {M<:QuantileRegression}
    fit(M, float(X), float(y); kwargs...)
end

"""
    refit!(m::QuantileRegression,
          [y::FPVector ;
           verbose::Bool=false,
           quantile::Union{Nothing,
           AbstractFloat}=nothing,
          ]
    )

Refit the `QuantileRegression` model with the new values for the response,
weights and quantile. This function assumes that `m` was correctly initialized.
"""
function refit!(m::QuantileRegression, y::FPVector; kwargs...)
    # Update y
    # Check that the old and new y have the same number of rows
    if size(m.y, 1) != size(y, 1)
        mess =
            "the new response vector should have the same dimension: " *
            "$(size(m.y, 1)) != $(size(y, 1))"
        throw(DimensionMismatch(mess))
    end
    copyto!(m.y, y)

    refit!(m; kwargs...)
end

function refit!(
    m::QuantileRegression{T};
    quantile::Union{Nothing,AbstractFloat}=nothing,
    wts::Union{Nothing,FPVector}=nothing,
    kwargs...,
) where {T}

    # Update wts
    n = length(m.y)
    if !isnothing(wts) && (length(wts) in (0, n))
        copy!(m.wts, wts)
    end

    # Update quantile
    if !isnothing(quantile)
        (0 < quantile < 1) ||
            error("quantile should be a number between 0 and 1 excluded: $(quantile)")
        m.τ = quantile
    end

    m.fitted = false
    fit!(m; kwargs...)
end

"""
    fit!(m::QuantileRegression;
         verbose::Bool=false,
         quantile::Union{Nothing, AbstractFloat}=nothing,
         correct_leverage::Bool=false,
         kwargs...)

Optimize the objective of a `QuantileRegression`.  When `verbose` is `true` the values
of the objective and the parameters are printed on stdout at each function evaluation.

This function assumes that `m` was correctly initialized.
This function returns early if the model was already fitted, instead call `refit!`.
"""
function fit!(
    m::QuantileRegression{T};
    verbose::Bool=false,
    quantile::Union{Nothing,AbstractFloat}=nothing,
    correct_leverage::Bool=false,
    kwargs...,
) where {T}

    # Return early if model has the fit flag set
    m.fitted && return m

    if correct_leverage
        wts = m.wts
        copy!(wts, leverage_weights(m))
        ## TODO: maybe multiply by the old wts?
    end

    # Update quantile
    if !isnothing(quantile)
        (0 < quantile < 1) || throw(
            DomainError(quantile, "quantile should be a number between 0 and 1 excluded"),
        )
        m.τ = quantile
    end

    interiormethod!(m.β, m.wrkres, m.X, m.y, m.τ; wts=m.wts, verbose=verbose)

    m.fitted = true
    m
end

function interiormethod(X, y, τ; kwargs...)
    n, p = size(X)
    T = eltype(y)
    interiormethod!(zeros(T, p), zeros(T, n), X, y, τ; kwargs...)
end

function interiormethod!(βout, rout, X, y, τ; wts=[], verbose::Bool=false)
    model = Tulip.Model{Float64}()
    pb = model.pbdata  # Internal problem data

    n, p = size(X)

    # Define variables
    β = Vector{Int}(undef, p)
    u = Vector{Int}(undef, n)
    v = Vector{Int}(undef, n)
    for i in 1:p
        β[i] = Tulip.add_variable!(pb, Int[], Float64[], 0.0 , -Inf , Inf, "β$i")
    end
    for i in 1:n
        u[i] = Tulip.add_variable!(pb, Int[], Float64[], τ   ,  0.0 , Inf, "u$i")
    end
    for i in 1:n
        v[i] = Tulip.add_variable!(pb, Int[], Float64[], 1-τ ,  0.0 , Inf, "v$i")
    end
    #    @variable(model, β[1:p])
    #    @variable(model, u[1:n] >= 0)
    #    @variable(model, v[1:n] >= 0)

    # Define objective
    # Nothing to do, already defined with the variables
    #    e = ones(n)
    #    @objective(model, Min, τ*dot(e, u) + (1-τ)*dot(e, v) )

    Wy, WX = if isempty(wts)
        y, X
    else
        (wts .* y), (wts .* X)
    end

    # Define constraints
    resid = Vector{Int}(undef, n)
    for i in 1:n
        ci = vcat(β, [u[i], v[i]])
        # Call Vector to transform to dense vector
        cci = vcat(Vector(WX[i, :]), [1.0, -1.0])
        vi = Wy[i]
        resid[i] = Tulip.add_constraint!(pb, ci, cci, vi, vi, "resid$i")
    end
    #    @constraint(model, resid, Wy .== WX * β + u - v)

    Tulip.optimize!(model)

    copyto!(βout, model.solution.x[β])
    copyto!(rout, model.solution.x[u] - model.solution.x[v])
    #    copyto!(βout, value.(β))
    #    copyto!(rout, value.(u) - value.(v))

    if verbose
        println("coef: ", βout)
        println("res: ", rout)
    end
    βout, rout
end

_objective(r::AbstractFloat, τ::AbstractFloat) = (τ - (r < 0)) * r

"""
    hall_sheather_bandwidth()
Optimal bandwidth for sparsity estimation according to Hall and Sheather (1988)
"""
function hall_sheather_bandwidth(q::Real, n::Int, α::Real=0.05)
    zα = quantile(Normal(), 1 - α / 2)
    ## Estimate of r=s/s" from a Normal distribution
    zq = quantile(Normal(), q)
    f = pdf(Normal(), zq)
    r = f^2 / (2 * zq^2 + 1)
    (zα^2 * 1.5 * r / n)^(1 / 5)
end

"""
    jones_bandwidth()
Optimal bandwidth for kernel sparsity estimation according to Jones (1992)
"""
function jones_bandwidth(q::Real, n::Int; kernel=:epanechnikov)
    ## From the kernel k(u) with domain [-1, 1]:
    ## kk = ∫k²(u)du / ( ∫u²k(u)du )^2

    kk = if kernel == :epanechnikov
        # For Epanechnikov kernel: kk = 3/5 / (1/5)^2
        15
    elseif kernel == :triangle
        # For triangle kernel (k(u) = 1 - |u|): kk = 2/3 / (1/6)^2
        24
    elseif kernel == :window
        # That is the Bofinger estimate of the bandwidth
        # For window kernel (k(u) = 1/2): kk = 1/2 / (1/3)^2
        4.5
    else
        error("kernel $(kernel) is not defined for sparsity estimation.")
    end

    ## Estimate of r=s/s" from a Normal distribution
    zq = quantile(Normal(), q)
    f = pdf(Normal(), zq)
    r = f^2 / (2 * zq^2 + 1)
    (kk * r^2 / n)^(1 / 5)
end

"""
    bofinger_bandwidth()
Optimal bandwidth for sparsity estimation according to Bofinger (1975)
"""
function bofinger_bandwidth(q::Real, n::Int)
    return jones_bandwidth(q, n; kernel=:window)
end

function epanechnikov_kernel(x::Real)
    if abs(x) < 1
        return 3 / 4 * (1 - x^2)
    else
        return zero(typeof(x))
    end
end

function triangle_kernel(x::Real)
    if abs(x) < 1
        return 1 - abs(x)
    else
        return zero(typeof(x))
    end
end

function window_kernel(x::Real)
    if abs(x) < 1
        return one(typeof(x)) / 2
    else
        return zero(typeof(x))
    end
end

"""
    sparcity(m::QuantileRegression, α::Real=0.05)

Compute the sparcity or quantile density ŝ(0)
using formula from Jones (1992) - Estimating densities, quantiles,
quantile densities and density quantiles.

"""
function sparcity(
    m::QuantileRegression;
    bw_method::Symbol=:jones,
    α::Real=0.05,
    kernel::Symbol=:epanechnikov,
)
    u = m.τ
    n = nobs(m)
    ## Select the optimal bandwidth from different methods
    h = if bw_method == :jones
        jones_bandwidth(u, n; kernel=kernel)
    elseif bw_method == :bofinger
        bofinger_bandwidth(u, n)
    elseif bw_method == :hall_sheather
        hall_sheather_bandwidth(u, n, α)
    else
        error(
            "only :jones, :bofinger and :hall_sheather methods for estimating the"*
            " optimal bandwidth are allowed: $(bw_method)",
        )
    end

    ## Select kernel for density estimation
    k = if kernel == :epanechnikov
        epanechnikov_kernel
    elseif kernel == :triangle
        triangle_kernel
    elseif kernel == :window
        window_kernel
    else
        error("only :epanechnikov, :triangle and :window kernels are allowed: $(kernel)")
    end

    r = sort(residuals(m))

    s0 = 0
    for i in eachindex(r)
        if i == 1
            s0 += r[i] / h * k(u / h)
        else
            s0 += (r[i] - r[i-1]) / h * k((u - (i - 1) / n) / h)
            if i == lastindex(r)
                s0 += r[i] / h * k((u - 1) / h)
            end
        end
    end
    s0
end

function location_variance(
    m::QuantileRegression,
    sqr::Bool=false;
    bw_method::Symbol=:jones,
    α::Real=0.05,
    kernel::Symbol=:epanechnikov,
)
    v = sparcity(m; bw_method=bw_method, α=α, kernel=kernel)^2
    v *= m.τ * (1 - m.τ)
    v *= (nobs(m) / dof_residual(m))
    return sqr ? v : sqrt(v)
end

"""
    nobs(m::QuantileRegression)
For linear and generalized linear models, returns the number of elements of the response.
"""
function nobs(m::QuantileRegression)::Int
    if !isempty(m.wts)
        count(!iszero, m.wts)
    else
        size(m.y, 1)
    end
end

coef(m::QuantileRegression) = m.β

dispersion(m::QuantileRegression) = mean(abs.(residuals(m)))

stderror(m::QuantileRegression) = location_variance(m, false) .* sqrt.(diag(vcov(m)))

function weights(m::QuantileRegression{T}) where {T<:AbstractFloat}
    if isempty(m.wts)
        weights(ones(T, length(m.y)))
    else
        weights(m.wts)
    end
end

workingweights(m::QuantileRegression) = m.wrkres

response(m::QuantileRegression) = m.y

isfitted(m::QuantileRegression) = m.fitted

islinear(m::QuantileRegression) = true

fitted(m::QuantileRegression) = m.y - m.wrkres

residuals(m::QuantileRegression) = m.wrkres

predict(m::QuantileRegression, newX::AbstractMatrix) = newX * coef(m)
predict(m::QuantileRegression) = fitted(m)

function nulldeviance(m::QuantileRegression)
    μ = quantile(m.y, m.τ)
    if isempty(m.wts)
        sum(_objective.(m.wrkres, Ref(m.τ)))
    else
        sum(m.wts .* _objective.(m.wrkres, Ref(m.τ)))
    end
end

deviance(m::QuantileRegression) = sum(_objective.(m.wrkres, Ref(m.τ)))

## TODO: define correctly the loglikelihood of the full model
fullloglikelihood(m::QuantileRegression) = 0

loglikelihood(m::QuantileRegression) = fullloglikelihood(m) - deviance(m) / 2

nullloglikelihood(m::QuantileRegression) = fullloglikelihood(m) - nulldeviance(m) / 2

modelmatrix(m::QuantileRegression) = m.X

vcov(m::QuantileRegression) =
    inv(Hermitian(float(Matrix(modelmatrix(m)' * (weights(m) .* modelmatrix(m))))))

projectionmatrix(m::QuantileRegression) =
    Hermitian(modelmatrix(m) * vcov(m) * modelmatrix(m)') .* weights(m)

leverage_weights(m::QuantileRegression) = sqrt.(1 .- leverage(m))


text:
"""
path related functions.
"""
module PATH
using Libdl
using Pkg
using ..Comonicon

function project(m::Module, xs...)
    path = pathof(m)
    path === nothing && return dirname(Pkg.project().path)
    return joinpath(dirname(dirname(path)), xs...)
end

project(xs...) = project(Comonicon, xs...)

sysimg() = "libcomonicon.$(Libdl.dlext)"
sysimg(name) = "lib$name.$(Libdl.dlext)"

"""
    default_exename()

Default Julia executable name: `joinpath(Sys.BINDIR, Base.julia_exename())`
"""
default_exename() = joinpath(Sys.BINDIR, Base.julia_exename())

"""
    default_julia_bin()

Return the default path to `.julia/bin`.
"""
default_julia_bin() = joinpath(first(DEPOT_PATH), "bin")

"""
    default_julia_fpath()

Return the default path to `.julia/completions`
"""
default_julia_fpath() = joinpath(first(DEPOT_PATH), "completions")
end


text:
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Saw: saw

@testset "Saw in white numbers" begin
    tol = 0.0001
        mat = [
            GreyNumber{Float64}(25.0, 25.0) GreyNumber{Float64}(65.0, 65.0) GreyNumber{Float64}(07.0, 07.0) GreyNumber{Float64}(20.0, 20.0);
            GreyNumber{Float64}(21.0, 21.0) GreyNumber{Float64}(78.0, 78.0) GreyNumber{Float64}(06.0, 06.0) GreyNumber{Float64}(24.0, 24.0);
            GreyNumber{Float64}(19.0, 19.0) GreyNumber{Float64}(53.0, 53.0) GreyNumber{Float64}(05.0, 05.0) GreyNumber{Float64}(33.0, 33.0);
            GreyNumber{Float64}(22.0, 22.0) GreyNumber{Float64}(25.0, 25.0) GreyNumber{Float64}(02.0, 02.0) GreyNumber{Float64}(31.0, 31.0)
        ]
        weights = [0.25, 0.25, 0.25, 0.25]
        fns = [maximum, maximum, minimum, maximum]
        result = saw(mat, weights, fns)
        @test isapprox(result[:scores][1], GreyNumber{Float64}(0.681277, 0.681277), atol = tol)
        @test isapprox(result[:scores][2], GreyNumber{Float64}(0.725151, 0.725151), atol = tol)
        @test isapprox(result[:scores][3], GreyNumber{Float64}(0.709871, 0.709871), atol = tol)
        @test isapprox(result[:scores][4], GreyNumber{Float64}(0.784976, 0.784976), atol = tol)
end


text:
## test_circulant_embedding.jl : test circulant embedding for GRF generation

@testset "circulant embedding      " begin

## 1d Exponential ##
cov = CovarianceFunction(1,Exponential(0.1))
pts = 0:1//1001:1
grf = GaussianRandomField(cov,CirculantEmbedding(),pts)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Exponential)
@test ndims(grf.cov) == 1
@test isa(grf,GaussianRandomField{C,CirculantEmbedding} where {C})
@test length(grf.pts[1]) == length(pts)
@test length(sample(grf)) == length(pts)

## 1d Matern with padding ##
m = CovarianceFunction(1,Matern(1.0,1.0))
grf = GaussianRandomField(m,CirculantEmbedding(),pts,padding=8)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Matern)
@test ndims(grf.cov) == 1
@test isa(grf,GaussianRandomField{C,CirculantEmbedding} where {C})
@test length(grf.pts[1]) == length(pts)
@test length(sample(grf)) == length(pts)

## 2d Anisotropic ##
A = [1000 0; 0 1000]
a = AnisotropicExponential(A)
c = CovarianceFunction(2,a)
pts1 = range(0,stop = 1,length = 128)
pts2 = range(0,stop = 1,length = 256)
@suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2,padding=4)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,AnisotropicExponential)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,CirculantEmbedding} where {C})
@test length(grf.pts) == 2
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)

## 3d AnisotropicExponential ##
A = [100 0 0; 0 100 0; 0 0 100]
a = AnisotropicExponential(A)
c = CovarianceFunction(3,a)
pts = 0:0.025:1
@suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts,pts,pts,padding=1,measure=false)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,AnisotropicExponential)
@test ndims(grf.cov) == 3
@test isa(grf,GaussianRandomField{C,CirculantEmbedding} where {C})
@test length(grf.pts) == 3
@test length(grf.pts[1]) == length(pts)
@test length(grf.pts[2]) == length(pts)
@test length(grf.pts[3]) == length(pts)
@test size(sample(grf),1) == length(pts)
@test size(sample(grf),2) == length(pts)
@test size(sample(grf),3) == length(pts)

## Larger domain ##
A = [1 0.8; 0.8 1]
m = AnisotropicExponential(A)
c = CovarianceFunction(2,m)
pts1 = range(-5,stop = 5,length = 128)
pts2 = range(10,stop = 0,length = 128)
@suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2,padding=2)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,AnisotropicExponential)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,CirculantEmbedding} where {C})
@test length(grf.pts) == 2
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)

## Dedicated plotting commands ##
m = Matern(0.5,2)
c = CovarianceFunction(3,m)
pts1 = range(0,stop = 1,length = 16)
pts2 = range(0,stop = 1,length = 8)
pts3 = range(0,stop = 1,length = 4)
#plot_covariance_matrix(c,pts1,pts2,pts3); close()
m = Matern(0.1,1)
c = CovarianceFunction(2,m)
pts = range(0,stop = 1,length = 512)
g = GaussianRandomField(c,CirculantEmbedding(),pts,pts,padding=4)
#contourf(g); close()
#plot_eigenvalues(g); close()

end


text:
# Create example point data
x = 8 .* rand(100) .- 4
y = 8 .* rand(100) .- 4
z = sin.(x) + cos.(y)
# Use the "hot" colormap
colormap("hot")
contourf(x, y, z)


text:
using BenchmarkTools

phasefile = joinpath("data", "small", "Phase.nii")
magfile = joinpath("data", "small", "Mag.nii")
echo = 3
phase = niread(phasefile).raw[:,:,:,echo]
mag = niread(magfile).raw[:,:,:,echo]

big(I, factor=3) = repeat(I; outer=factor .* [1,1,1])
bigphase = big(phase)
bigmag = big(mag)
@btime ROMEO.calculateweights($bigphase; mag=$bigmag)
@btime unwrap($bigphase; mag=$bigmag)

# akku - high power
 # 450ms
 # 1.4s

 # PC
  # 360ms
  # 1.1s


text:
using Plots, Plots.PlotMeasures, LaTeXStrings

include("SDOF_Model.jl")

(@isdefined TARGET_FOLDER) ? nothing : TARGET_FOLDER = ""

p = sdof()
prob = IVP(p)
asol = analytic_solution(p)
T = period(p)

tmax = 4.0
dom = 0:0.001:tmax

f(α; alg=VREP(δ=α*T, approx_model=StepIntersect(setops=:concrete))) = solve(prob, tspan=(0.0, tmax), alg=alg)

# solution using numerical scheme
q = InitialValueProblem_quad(p)

g(α; alg=Newmark(Δt=α*T, α=0.25, δ=0.5)) = solve(q, alg, NSTEPS=round(Int, tmax / (α*T)));

function plot_x_vs_t(α)
    fig = plot(xlab=L"t", ylab=L"u(t)",
               legendfontsize=12, legend=:none,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([0.0, 1.0, 2.0, 3.0, 4.0], [L"0.0", L"1.0", L"2.0", L"3.0", L"4.0"]),
               ytick=([-1.0, -0.5, 0.0, 0.5, 1.0], [L"-1.0", L"-0.5", L"0.0", L"0.5", L"1.0"]),
               xlims=(0.0, 4.5), ylims=(-1.2, 1.3),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    L1 = LineSegment([3.3, -1.2], [3.7, -1.2])
    L2 = LineSegment([3.3, 1.2], [3.7, 1.2])
    L3 = LineSegment([3.3, -1.2], [3.3, 1.2])
    L4 = LineSegment([3.7, -1.2], [3.7, 1.2])

    #plot!(fig, UnionSetArray([L1, L2, L3, L4]), color=:black, style=:dash, lw=1.5, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L1, color=:black, style=:dash, lw=2.5, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L2, color=:black, style=:dash, lw=2.5, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L3, color=:black, style=:dash, lw=2.5, seriestype=:solid, marker=:none, alpha=1.)
    plot!(fig, L4, color=:black, style=:dash, lw=2.5, seriestype=:solid, marker=:none, alpha=1.)

    # Flowpipe
    println("Flowpipe solution (VREP)")
    fp = f(α) # warm-up
    @time fp = f(α)
    plot!(fig, fp, vars=(0, 1), lab=L"\alpha = %$α~\textrm{Flowpipe}", color=:lightblue, lw=0.5)

    # Numerical solution
    println("Numerical solution (Newmark)")
    y = g(α, alg=Newmark(Δt=α*T, α=0.25, δ=0.5)) # warm-up
    @time y = g(α, alg=Newmark(Δt=α*T, α=0.25, δ=0.5))
    tdom = 0:(α*T):tmax
    plot!(fig, tdom, displacements(y, 1), linestyle=:solid, marker=:circle, lab=L"\alpha = %$α~\textrm{Newmark}", color=:red)

    println("Numerical solution (Bathe)")
    y = g(α, alg=Bathe(Δt=α*T)) # warm-up
    @time y = g(α, alg=Bathe(Δt=α*T))
    plot!(fig, tdom, displacements(y, 1), linestyle=:solid, marker=:utriangle, lab=L"\alpha = %$α~\textrm{Bathe}", color=:green)

    # Analytic solution
    plot!(fig, dom, asol.(dom), lab=L"\textrm{Analytic}", color=:magenta, lw=1.2)

    fig
end

function plot_x_vs_t_zoom(α)
    fig = plot(xlab=L"t", ylab=L"u(t)", legend=:topleft,
               legendfont=font(15, "Times"),
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([3.3, 3.5, 3.7], [L"3.3", L"3.5", L"3.7"]),
               ytick=([-1.0, -0.5, 0.0, 0.5, 1.0], [L"-1.0", L"-0.5", L"0.0", L"0.5", L"1.0"]),
               xlims=(3.28, 3.72), ylims=(-1.2, 1.2),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    a = 3.3; b = 3.7

    # flowpipe solutions
    aux = f(α)(a .. b)
    plot!(fig, aux[1:end-1], vars=(0, 1), lab=L"\textrm{Flowpipe}", color=:lightblue)

    # numerical solutions
    tdom = a:(α*T):b
    ai = round(Int, a / (α * T))
    bi = round(Int, b / (α * T))

    y = g(α, alg=Newmark(Δt=α*T, α=0.25, δ=0.5)) |> displacements
    plot!(fig, tdom, [yi[1] for yi in y[ai:bi]], linestyle=:solid, lw=2.5, marker=:circle, markersize=10.0, lab=L"\textrm{Newmark}", color=:red)

    y = g(α, alg=Bathe(Δt=α*T)) |> displacements
    plot!(fig, tdom, [yi[1] for yi in y[ai:bi]], linestyle=:solid, lw=2.5, marker=:utriangle, markersize=10.0, lab=L"\textrm{Bathe}", color=:green)

    # analytic solution
    domw = a:0.001:b
    plot!(fig, domw, asol.(domw), lab=L"\textrm{Analytic}", color=:magenta, lw=2.0)

    fig
end

# ==========================
# Figures for α = 0.05
# ==========================
println("Solving for α = 0.05")
fig = plot_x_vs_t(0.05)
savefig(fig, joinpath(TARGET_FOLDER, "sdof5a.pdf"))

fig = plot_x_vs_t_zoom(0.05)
savefig(fig, joinpath(TARGET_FOLDER, "sdof5b.pdf"))

# ==========================
# Figures for α = 0.1
# ==========================

println("Solving for α = 0.1")
fig = plot_x_vs_t(0.1)
savefig(fig, joinpath(TARGET_FOLDER, "sdof6a.pdf"))

fig = plot_x_vs_t_zoom(0.1)
savefig(fig, joinpath(TARGET_FOLDER, "sdof6b.pdf"))

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
module Solutions

    using HDF5
    using SharedArrays

    using Base: TwicePrecision
    using Random

    using ..CommonFunctions
    using ..Config
    using ..Equations
    using ..Utils

    export DEFAULT_NSAVE, DEFAULT_NWRITE

    const DEFAULT_NSAVE = 1
    const DEFAULT_NWRITE = 0
    const DEFAULT_SCONV = :strong

    export SolutionVector

    SolutionVector{DT} = Union{Vector{DT}, Vector{TwicePrecision{DT}}}

    export get_data!, set_data!
    export DataSeries, PDataSeries, SDataSeries

    include("solutions/dataseries.jl")

    export TimeSeries, compute_timeseries!

    include("solutions/timeseries.jl")

    export SemiMartingale
    export WienerProcess, generate_wienerprocess!

    include("solutions/wienerprocess.jl")

    export Solution, ParallelSolution, DeterministicSolution, StochasticSolution
    export nsave, nsamples, ntime, timesteps, offset, conv, hdf5
    export get_solution, get_solution!, set_solution!

    include("solutions/solution.jl")

    export AtomicSolution,
           AtomicSolutionODE, AtomicSolutionPODE,
           AtomicSolutionDAE, AtomicSolutionPDAE,
           AtomicSolutionSDE, AtomicSolutionPSDE
    export update!, cut_periodic_solution!,
           get_increments, get_increments!, set_increments!

    include("solutions/atomic_solution.jl")
    include("solutions/atomic_solution_ode.jl")
    include("solutions/atomic_solution_pode.jl")
    include("solutions/atomic_solution_dae.jl")
    include("solutions/atomic_solution_pdae.jl")
    include("solutions/atomic_solution_sde.jl")
    include("solutions/atomic_solution_psde.jl")

    export SolutionODE, SSolutionODE, PSolutionODE, SolutionPODE, SSolutionPODE, PSolutionPODE
    export SolutionDAE, SSolutionDAE, PSolutionDAE, SolutionPDAE, SSolutionPDAE, PSolutionPDAE
    export SolutionSDE, SSolutionSDE, PSolutionSDE, SolutionPSDE, SSolutionPSDE, PSolutionPSDE
    export get_initial_conditions, get_initial_conditions!, set_initial_conditions!,
           create_hdf5, create_hdf5!

    include("solutions/solution_ode.jl")
    include("solutions/solution_pode.jl")
    include("solutions/solution_dae.jl")
    include("solutions/solution_pdae.jl")
    include("solutions/solution_sde.jl")
    include("solutions/solution_psde.jl")

    include("solutions/solution_hdf5.jl")

end


text:
"""
    House

Data stucture that represents a house and its appliances.

# Attributes
- `index`: Index of that house
- `num_timesteps`: Number of time-steps in the time horizon
- `netload_min`: House's minimum net load, for each time-step
- `netload_max`:
- `price`: Electricity price, in dollar per kWh
- `appliances`: List of appliances in that house
"""
struct House <: Resource
    index::Int
    num_timesteps::Int

    netload_min::Vector{Float64}
    netload_max::Vector{Float64}
    price::Vector{Float64}

    appliances::Vector{Resource}
end

function House(;
    index::Integer=0,
    T::Integer=0,
    netload_min::Vector{Float64}=Vector{Float64}(),
    netload_max::Vector{Float64}=Vector{Float64}(),
    price::Vector{Float64}=Vector{Float64}(),
    appliances::Vector{Resource}=Vector{Resource}()
)

    # Dimension checks
    T == size(netload_min, 1) || throw(DimensionMismatch("Invalid netload_min"))
    T == size(netload_max, 1) || throw(DimensionMismatch("Invalid netload_max"))
    T == size(price, 1) || throw(DimensionMismatch("Invalid price"))

    House(
        index, T,
        copy(netload_min), copy(netload_max), copy(price),
        appliances
    )
end

"""
    LindaOracleMIP(h, solver)

Construct a MIP oracle that schedules the House.
"""
function build_model!(h::House, model::MOI.ModelLike)
    # Emtpy model
    # TODO: use an UninstanitatedOptimizer instead
    MOI.empty!(model)

    # Minimize
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    # Initial objective function is empty
    fobj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(fobj)}(), fobj)
    # TODO: set parameters

    T = h.num_timesteps

    # Initialize
    var2idx = Dict{Tuple{Symbol, Int, Symbol, Int}, MOI.VariableIndex}()
    con2idx = Dict{Tuple{Symbol, Int, Symbol, Int}, MOI.ConstraintIndex}()

    # ==========================================
    #    I. Add local variables
    # ==========================================
    pnet = MOI.add_variables(model, T)  # House's net load
    # Variable bounds
    for t in 1:T
        MOI.add_constraint(model,
            MOI.SingleVariable(pnet[t]),
            MOI.Interval(h.netload_min[t], h.netload_max[t])
        )
    end

    # Set objective coefficients (?)
    # append!(obj, h.price)
    for t in 1:T
        MOI.modify(model,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(pnet[t], h.price[t])
        )
    end

    # Update list of variable indices
    for t in 1:T
        var2idx[(:house, h.index, :pnet, t)] = pnet[t]  # net load
    end

    # ==========================================
    #    II. Add local constraints
    # ==========================================
    # Linking constraint for house's net load
    #   pnet[h, t] - sum_{r} pnet[r, t] = 0.0
    for t in 1:T
        cidx = MOI.add_constraint(model,
            MOI.ScalarAffineFunction{Float64}(
                [MOI.ScalarAffineTerm(1.0, pnet[t])],
                0.0
            ),
            MOI.EqualTo(0.0)
        )

        # Record constraint index
        con2idx[(:house, h.index, :link, t)] = cidx
    end

    # ==========================================
    #    III. Add sub-resources to current model
    # ==========================================
    for r in h.appliances
        add_resource_to_model!(model, h, r, var2idx, con2idx)
    end

    return model, var2idx, con2idx
end


text:
using SafeTestsets

@safetestset "Day04" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Day04

    @testset "digits" begin
        using AdventOfCode.Day04: digits

        @test digits(171309) == (-1, 1,7,1,3,0,9, 10)
    end

    @testset "is_valid1" begin
        using AdventOfCode.Day04: is_valid1

        @test  is_valid1((-1, 1,1,1,1,1,1, 10))
        @test !is_valid1((-1, 2,2,3,4,5,0, 10))
        @test !is_valid1((-1, 1,2,3,7,8,9, 10))
    end

    @testset "is_valid2" begin
        using AdventOfCode.Day04: is_valid2

        @test  is_valid2((-1, 1,1,2,2,3,3, 10))
        @test !is_valid2((-1, 1,2,3,4,4,4, 10))
        @test  is_valid2((-1, 1,1,1,1,2,2, 10))
    end

    btest("4.1") do
        part1(171309, 643603) == 1625
    end

    btest("4.2") do
        part2(171309, 643603) == 1111
    end
end


text:
# we write a function that decomposes a square matrix A as A = LU
# where L is lower triangular and U is upper triangular
function LUdecomp(A)

    # we get the size of A
    n, m = size(A)

    # we check if A is indeed a square matrix
    if n != m
        error("Input must be a square matrix")
    end

    # we initialize the output matrices
    L = zeros(n, n)
    U = zeros(n, n)

    # we make a copy of A in order not to modify the original
    B = copy(A)

    for k = 1:n

        # we set the entries of the kth row of U
        for j = k:n
            U[k,j] = B[k,j]
        end

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("** A^(k-1)_{k,k} == 0 in LU decomp")
        end

        # we set the entries of the kth column of L
        for i = k:n
           L[i,k] = B[i,k] / U[k,k]
        end

        # we modify A for the next iteration
        for i = k:n
            for j = k:n
                B[i,j] -= L[i,k] * U[k,j]
            end
        end
    end

    # we return L and U
    return (L, U)
end

max_stars_repo_path: 1a/LU Decomposition/LUdecomp.jl

max_issues_repo_path: 1a/LU Decomposition/LUdecomp.jl

max_forks_repo_path: 1a/LU Decomposition/LUdecomp.jl


text:
using Pkg, Test
Pkg.activate("../../.")
include("../../src/Toolips.jl")
using Main: Toolips
Toolips.new_webapp()

@testset "Interface" begin

    @testset "text/html functions" begin
        @testset "html" begin
            htm = html("
Test
")
            @test typeof(htm) == Function

        end
    end

    @testset "indexing/iter" begin

    end

    @testset "style interface" begin

    end

    @testset "serving/routing" begin

    end

    @testset "requests" begin

    end

end # - Interface

@testset "Server" begin

    @testset "extensions" begin

    end

    @testset "core" begin

    end
end # - Server

@testset "App" begin

end # - App


text:
abstract type CNFLayer <: Function end
Flux.trainable(m::CNFLayer) = (m.p,)

"""
Constructs a continuous-time recurrent neural network, also known as a neural
ordinary differential equation (neural ODE), with fast gradient calculation
via adjoints [1] and specialized for density estimation based on continuous
normalizing flows (CNF) [2]. At a high level this corresponds to the following steps:

1. Parametrize the variable of interest x(t) as a function f(z,θ,t) of a base variable z(t) with known density p_z;
2. Use the transformation of variables formula to predict the density p_x as a function of the density p_z and the trace of the Jacobian of f;
3. Choose the parameter θ to minimize a loss function of p_x (usually the negative likelihood of the data);

After these steps one may use the NN model and the learned θ to predict the density p_x for new values of x.

```julia
NeuralODE(model,basedist=nothing,monte_carlo=false,tspan,args...;kwargs...)
```
Arguments:
- `model`: A Chain neural network that defines the ̇x.
- `basedist`: Distribution of the base variable. Set to the unit normal by default.
- `monte_carlo`: Method for calculating the trace of the Jacobian. The default monte_carlo = false calculates the Jacobian and its trace directly. monte_carlo = true uses the stochastic approach presented in [3] to provide an unbiased estimate for the trace.
- `tspan`: The timespan to be solved on.
- `kwargs`: Additional arguments splatted to the ODE solver. See the
  [Common Solver Arguments](https://docs.sciml.ai/dev/basics/common_solver_opts/)
  documentation for more details.
Ref
[1]L. S. Pontryagin, Mathematical Theory of Optimal Processes. CRC Press, 1987.
[2]R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud. Neural Ordinary Differential Equations. arXiv preprint at arXiv1806.07366, 2019.
[3]W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, D. Duvenaud. FFJORD: Free-Form Continuous Dynamic For Scalable Reversible Generative Models. arXiv preprint at arXiv1810.01367, 2018.

"""
struct FFJORD{M,P,RE,Distribution,Bool,T,A,K} <: CNFLayer
    model::M
    p::P
    re::RE
    basedist::Distribution
    monte_carlo::Bool
    tspan::T
    args::A
    kwargs::K

    function FFJORD(model,tspan,args...;p = nothing,basedist=nothing,monte_carlo=false,kwargs...)
        _p,re = Flux.destructure(model)
        if p === nothing
            p = _p
        end
        if basedist === nothing
            size_input = size(model[1].W)[2]
            basedist = MvNormal(zeros(size_input), I + zeros(size_input,size_input))
        end
        new{typeof(model),typeof(p),typeof(re),typeof(basedist),
            typeof(monte_carlo),typeof(tspan),typeof(args),typeof(kwargs)}(
            model,p,re,basedist,monte_carlo,tspan,args,kwargs)
    end
end

function jacobian_fn(f, x::AbstractVector)
   y::AbstractVector, back = Zygote.pullback(f, x)
   ȳ(i) = [i == j for j = 1:length(y)]
   vcat([transpose(back(ȳ(i))[1]) for i = 1:length(y)]...)
 end

function ffjord(du,u,p,t,re,monte_carlo,e)
    z = @view u[1:end-1]
    m = re(p)
    if monte_carlo
        _, back = Zygote.pullback(m,z)
        eJ = back(e)[1]
        trace_jac = (eJ.*e)[1]
    else
        J = jacobian_fn(m, z)
        trace_jac = length(z) == 1 ? sum(J) : tr(J)
    end
    du[1:end-1] = m(z)
    du[end] = -trace_jac
end

function (n::FFJORD)(x,p=n.p,monte_carlo=n.monte_carlo)
    e = monte_carlo ? randn(Float32,length(x)) : nothing
    ffjord_ = (du,u,p,t)->ffjord(du,u,p,t,n.re,monte_carlo,e)
    prob = ODEProblem{true}(ffjord_,vcat(x,0f0),n.tspan,p)
    pred = solve(prob,n.args...;n.kwargs...)[:,end]
    pz = n.basedist
    z = pred[1:end-1]
    delta_logp = pred[end]
    logpz = log.(pdf(pz, z))
    logpx = logpz .- delta_logp
    return logpx[1]
end


text:
module DynamicBayesianNetworkBehaviors

using AutomotiveDrivingModels
using NLopt
using Distributions
using BayesNets
using SmileExtra

import LightGraphs: topological_sort_by_dfs, in_neighbors, nv, ne
import Discretizers: encode
import AutomotiveDrivingModels:
    ModelTargets,
    AbstractVehicleBehavior,
    AbstractVehicleBehaviorPreallocatedData,
    AbstractVehicleBehaviorTrainParams,

    preallocate_learning_data,
    select_action,
    calc_action_loglikelihood,
    train,
    _reverse_smoothing_sequential_moving_average,
    is_in_fold

export
    DBNModel,
    DynamicBayesianNetworkBehavior,
    DBNSimParams,
    GraphLearningResult,
    ParentFeatures,
    BN_TrainParams,

    DirichletPrior,
    UniformPrior,
    BDeuPrior,

    DEFAULT_INDICATORS,
    DEFAULT_DISCRETIZERS,

    dbnmodel,
    build_bn,

    indexof,
    is_parent,
    parent_indeces,

    get_target_lat,
    get_target_lon,
    get_indicators,
    get_indicators_for_target,
    get_indicators_lat,
    get_indicators_lon,

    find_target_indeces,

    get_num_vertices,
    get_num_edges,

    get_total_sample_count,
    get_bin_counts,
    get_marginal_probability,
    get_counts_for_assignment,
    get_emstats,

    encode,
    sample!,
    sample_unset!,
    sample_and_logP!,
    calc_log_probability_of_assignment,
    calc_probability_distribution_over_assignments!,
    calc_probability_for_uniform_sample_from_bin,

    export_to_text,
    print_structure,

    get_input_acceleration,
    get_input_turnrate,
    # infer_action_lon_from_input_acceleration,
    # infer_action_lat_from_input_turnrate,

    discretize,
    discretize_cleaned,
    drop_invalid_discretization_rows,
    convert_dataset_to_matrix,
    calc_bincounts_array,

    select_action,
    calc_action_loglikelihood,
    train

##############################################################

include("dirichlet_priors.jl")
include("DBNModel.jl")
include("behaviormodel.jl")
include("default_params.jl")
include("graph_learning_result.jl")
include("learning.jl")

include("LinearGaussianBayesianNetBehavior.jl")

end # module


text:
# Autogenerated wrapper script for Conduit_jll for aarch64-apple-darwin
export libconduit, libconduit_blueprint, libconduit_relay

JLLWrappers.@generate_wrapper_header("Conduit")
JLLWrappers.@declare_library_product(libconduit, "@rpath/libconduit.dylib")
JLLWrappers.@declare_library_product(libconduit_blueprint, "@rpath/libconduit_blueprint.dylib")
JLLWrappers.@declare_library_product(libconduit_relay, "@rpath/libconduit_relay.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libconduit,
        "lib/libconduit.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libconduit_blueprint,
        "lib/libconduit_blueprint.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libconduit_relay,
        "lib/libconduit_relay.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()