text:
@testset "rand" begin
    # Some arbitrary alignment
    x = iterdates()
    expected_times = _eval_fast(x).times

    @testset "scalar" begin
        n = rand(x)
        # calling multiple times should give different nodes, since we use a different rng.
        @test n != rand(x)

        # If we specify an explicit rng, then we should get identical nodes back.
        rng = MersenneTwister()
        n = rand(rng, x)
        @test n === rand(rng, x)

        @test value_type(n) == Float64
        block = _eval(n)
        @test block.times == expected_times
        @test all(0 .<= block.values .< 1)

        n = rand(rng, x, Int32)
        @test n === rand(rng, x, Int32)
        @test value_type(n) == Int32
        block = _eval(n)
        @test block.times == expected_times
        @test all(typemin(Int32) .<= block.values .<= typemax(Int32))
    end

    @testset "array" begin
        n = rand(x, ())
        @test value_type(n) == Array{Float64,0}

        # If we specify an explicit rng, then we should get identical nodes back.
        rng = MersenneTwister()
        n = rand(rng, x, (2,))
        @test n === rand(rng, x, (2,))
        @test n != rand(rng, x, (3,))
        @test n != rand(rng, x, (2, 3))

        @test rand(rng, x, (2,)) === rand(rng, x, 2)
        @test rand(rng, x, Float64, (2,)) === rand(rng, x, 2)
        @test rand(rng, x, (2,)) === rand(rng, x, Float64, 2)
        @test rand(rng, x, (2, 3)) === rand(rng, x, 2, 3)
        @test rand(rng, x, Float64, (2, 3)) === rand(rng, x, 2, 3)
        @test rand(rng, x, (2, 3)) === rand(rng, x, Float64, 2, 3)

        @test value_type(n) == Vector{Float64}
        block = _eval(n)
        @test block.times == expected_times
        @test all(map(value -> all(0 .<= value .< 1), block.values))
    end
end


text:
using DataDeps

register(DataDep(
    "splice",
    """
    Dataset: splice
    Donors: G. Towell, M. Noordewier, and J. Shavlik
    Website (Original): https://archive.ics.uci.edu/ml/datasets/Molecular+Biology+(Splice-junction+Gene+Sequences)
    Obtained from: https://www.openml.org/d/40670
    See also: https://www.rdocumentation.org/packages/mlbench/versions/2.1-1/topics/DNA

    Observations: 3186
    Features:     180
    Classes:      3
    """,
    "https://www.openml.org/data/get_csv/4965245/dna.arff",
    "2775c83180277226d445e66acdc6370f298e5d4811ef20b8ce9d111485600b5c",
    post_fetch_method = (path -> begin
        MVDA.process_dataset(path,
            header=true,
            target_index=181,
            feature_indices=1:180,
            ext=".csv",
        )
    end),
))

ext: jl


text:
"""
tupling of two or more discrete functions
if we have
    f1 : D → R1
    f2 : D → R2
        …
    fk : D → Rk
then FunTupling((f1, …, fk)) is F
    F : D → R1 × … × Rk
"""
struct FunTupling{N, T, D, R} <: CompositeFunction
    factors :: NTuple{N, T}
    dom :: D
    rng :: R
end

const FTuple = FunTupling

function FTuple(factors)
    dom = domain(first(factors))
    rng = DirectProduct(map(codomain, factors))
    return FTuple(factors, dom, rng)
end

# creating functional tupling from two tuples
# not necessarily equal length !

function FTuple(from :: NTuple{N, Int}, to :: NTuple{M, Int}) where N where M
    funs = tuple(map(x -> ExtendedResidueFunction(from, x), to)...)
    return FTuple(funs)
end

# creating functional product from dom and codom

function FTuple(dom :: AbstractProduct, codom :: AbstractProduct)
    #fromto = zip(from, to)
    #funs = tuple(map(ResidueFunction, fromto)...)
    #return FProd(funs)
    from = size(dom)
    to = size(codom)
    return FTuple(from, to)
end

function (ft :: FTuple{N, T, D, R})(v) where {N,T,D,R}
    funs = factors(ft)
    n = length(funs)
    return NTuple{N}(funs[i](v) for i in 1 : n)
end

function Base.setindex!(ft :: FTuple, v :: NTuple, k)
    fs = factors(ft)
    for (i, val) in enumerate(v)
        fs[i][k] = val
    end
end

(fs :: FTuple)(v...) = fs(tuple(v...))

FTuple(factors...) = FTuple(tuple(factors...))

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
function single_solve(state, solve!, problem)
    t0 = time_ns()
    state = solve!(state, problem)
    t1 = time_ns()
    t_exec = to_micros(t1-t0)
    return state, t_exec
end

to_micros(t) = Float64(t/10^3)

function construct_solver(A, N)
    return (s, p) -> algo = A() ; return algo(s, N, p) end
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
@doc raw"""
Axisymmetric tokamak equilibrium in (r,θ,ϕ) coordinates with covariant
components of the vector potential given by
```math
A (r, \theta, \phi) = B_0 \, \bigg( 0 , \, \frac{r R_0}{\cos (\theta)} - \bigg( \frac{R_0}{\cos (\theta)} \bigg)^2 \, \ln \bigg( \frac{R}{R_0} \bigg) , \, - \frac{r^2}{2 q_0} \bigg)^T ,
```
resulting in the magnetic field with covariant components
```math
B (r, \theta, \phi) = \frac{B_0}{q_0} \, \bigg( 0 , \, \frac{r^2}{R}, \, q_0 R_0 \bigg)^T ,
```
where $R = R_0 + r \cos \theta$.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `q₀`: safety factor at magnetic axis
"""
module AxisymmetricTokamakToroidal

    using RecipesBase

    import ..ElectromagneticFields
    import ..ElectromagneticFields: AnalyticEquilibrium, code

    export  AxisymmetricTokamakToroidalEquilibrium

    const DEFAULT_R₀ = 1.0
    const DEFAULT_B₀ = 1.0
    const DEFAULT_q₀ = 2.0

    struct AxisymmetricTokamakToroidalEquilibrium{T <: Number} <: AnalyticEquilibrium
        name::String
        R₀::T
        B₀::T
        q₀::T

        function AxisymmetricTokamakToroidalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number
            new("AxisymmetricTokamakEquilibriumToroidal", R₀, B₀, q₀)
        end
    end

    AxisymmetricTokamakToroidalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakToroidalEquilibrium{T}(R₀, B₀, q₀)

    function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        AxisymmetricTokamakToroidalEquilibrium(R₀, B₀, q₀)
    end

    function ITER()
        AxisymmetricTokamakToroidalEquilibrium(ITER_R₀, ITER_B₀, ITER_q₀)
    end

    macro code(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        code(init(R₀, B₀, q₀); escape=true)
    end

    macro code_iter()
        code(ITER(); escape=true)
    end

    function Base.show(io::IO, equ::AxisymmetricTokamakToroidalEquilibrium)
        print(io, "Axisymmetric Tokamak Equilibrium in Toroidal Coordinates with\n")
        print(io, "  R₀ = ", equ.R₀, "\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  q₀ = ", equ.q₀)
    end

    r(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[1]
    θ(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[2]
    ϕ(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[3]
    R(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = equ.R₀ + r(x,equ) * cos(θ(x,equ))
    X(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x,equ) * cos(ϕ(x,equ))
    Y(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x,equ) * sin(ϕ(x,equ))
    Z(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x,equ) * sin(θ(x,equ))

    ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x,equ) * R(x,equ)

    ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = + equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * cos(θ(x,equ)) - log(R(x,equ) / equ.R₀) * sin(θ(x,equ)) ) / 2
    ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = - equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * sin(θ(x,equ)) + log(R(x,equ) / equ.R₀) * cos(θ(x,equ)) ) * r(x,equ) / 2
    ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = + equ.B₀ * r(x,equ)^2 / equ.q₀ / 2

    ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = X(ξ,equ)
    ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Y(ξ,equ)
    ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Z(ξ,equ)

    ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = sqrt((sqrt(x[1]^2 + x[2]^2)-equ.R₀)^2 + x[3]^2)
    ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = atan(x[3], sqrt(x[1]^2 + x[2]^2)-equ.R₀)
    ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = atan(x[2], x[1])

    ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = one(eltype(x))
    ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x, equ)^2
    ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x, equ)^2

    ElectromagneticFields.get_functions(::AxisymmetricTokamakToroidalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ)

    function ElectromagneticFields.periodicity(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium)
        p = zero(x)
        p[2] = 2π
        p[3] = 2π
        return p
    end

    @recipe function f(equ::AxisymmetricTokamakToroidalEquilibrium;
                       nx = 100, ny = 120, levels = 50, size = (400,400),
                       xlims = (  0.5 * equ.R₀,   1.5 * equ.R₀),
                       ylims = (- 0.5 * equ.R₀, + 0.5 * equ.R₀))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        zgrid = LinRange(ylims[1], ylims[2], ny)
        rgrid = [ElectromagneticFields.ξ¹([xgrid[i], 0.0, zgrid[j]], equ) for i in eachindex(xgrid), j in eachindex(zgrid)]
        θgrid = [ElectromagneticFields.ξ²([xgrid[i], 0.0, zgrid[j]], equ) for i in eachindex(xgrid), j in eachindex(zgrid)]
        pot   = [ElectromagneticFields.A₃([rgrid[i,j], θgrid[i,j], 0.0], equ) / xgrid[i] for i in eachindex(xgrid), j in eachindex(zgrid)]

        seriestype   := :contour
        aspect_ratio := :equal
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        (xgrid, zgrid, pot')
    end

end


text:
import StaticArrays:
    MVector
import ..AbstractTransferOperator
import ..TransferOperatorTriangulationExact

import StateSpaceReconstruction:
    MutableSimplex,
    orientation

export transferoperator_triangulation_exact

"""
    transferoperator_triangulation_exact(invariant_pts) -> TransferOperatorTriangulationExact

Compute an approximation of the transfer operator over a triangulated partition of the
provided points (`invariant_pts`). The estimator assumes that the points are invariant
under the forward projection of the simplices in the triangulation, which means that
the last point must be contained within the convex hull of all previous points.

*Tip: The `invariantize` function from `StateSpaceReconstruction` can be used to move the
last point toward the center of the embedding until invariance is achieved.*

This function computes *exact* simplex intersections and is thus *slow*. It should
only be used on small datasets.

## Example

```julia
# Create some random points and make sure they are invariant
pts = [rand(3) for i = 1:30]
invariant_pts = invariantize(pts)

# Compute the transfer operator
to = transferoperator_triangulation_exact(invariant_pts)
```

## References
The idea of using a triangulation estimator to compute the transfer operator was
introduced in [1], where it is used to compute physical invariant measures. To compute
a physical invariant measure from a transfer operator `to`, simply call
`invariantmeasure(to)`.

1. Froyland, Gary. "Estimating physical invariant measures and space averages of dynamical systems indicators." Bulletin of the Australian Mathematical Society 56.1 (1997): 157-159.
"""
function transferoperator_triangulation_exact(invariant_pts)
    DT = DelaunayTriangulation(invariant_pts[1:end-1])
    transferoperator_triangulation_exact(invariant_pts, DT)
end

"""
    transferoperator_triangulation_exact(invariant_pts, DT::DelaunayTriangulation)

Approximate the transfer operator over a triangulation of `invariant_pts`. Does the same
as the single-argument version of the function, but requires that the triangulation
is precomputed on `invariant_pts[1:(end - 1)]` (excluding the last point).

## Example
```julia
# Create some random points and make sure they are invariant
pts = [rand(3) for i = 1:30]
invariant_pts = invariantize(pts)

# Triangulate all but the last point (exclude the last, so that an image simplex
# always exists)
DT = DelaunayTriangulation(invariant_pts[1:(end - 1)])

# Compute the transfer operator
to = transferoperator_triangulation_exact(invariant_pts, DT)
```

## References
The idea of using a triangulation estimator to compute the transfer operator was
introduced in [1], where it is used to compute physical invariant measures. To compute
a physical invariant measure from a transfer operator `to`, simply call
`invariantmeasure(to)`.

1. Froyland, Gary. "Estimating physical invariant measures and space averages of dynamical systems indicators." Bulletin of the Australian Mathematical Society 56.1 (1997): 157-159.
"""
function transferoperator_triangulation_exact(invariant_pts,
        DT::DelaunayTriangulation)

    dim = length(invariant_pts[1])
    n_simplices = length(DT)

    #=
    Tolerance for simplex size when filling the transfer matrix
    =#
    ϵ::Float64 = 1e-8 / n_simplices

    # Pre-allocate the simplices
    image_simplex = MutableSimplex(zeros(Float64, dim, dim + 1))
    simplex = MutableSimplex(zeros(Float64, dim, dim + 1))

    # Pre-allocate transfer matrix
    TO = zeros(Float64, n_simplices, n_simplices)

    for j in 1:n_simplices
        # Fill the image simplex
        for k = 1:(dim + 1)
            image_simplex[k] .= 0.0
            image_simplex[k] = invariant_pts[DT[j][k] + 1]
        end

        imvol = abs(orientation(image_simplex))

        for i = 1:n_simplices
            # Fill which we're testing if arrives in the image simplex,
            # and if it does, how much of it overlaps with the image
            # simplex.
            for k = 1:(dim + 1)
                simplex[k] .= 0.0
                simplex[k] = invariant_pts[DT[i][k]]
            end

            # Only compute the entry of the transfer matrix
            # if simplices are of sufficient size.
            vol = abs(orientation(simplex))

            if vol * imvol > 0 && (vol/imvol) > ϵ
                TO[j, i] = (simplex ∩ image_simplex) / imvol
            end
        end
    end
    TransferOperatorTriangulationExact(TO)
end


text:
X = [6.91 3.34516 0.00916154 0.0306394 0.0217702; 7.88 3.39959 0.034967 0.00916154 0.0209481; 7.56 3.38934 0.0247539 0.034967 0.022901; 6.88 3.28777 0.0239149 0.0247539 0.0233713; 6.56 3.23396 0.178572 0.0239149 0.0243629; 6.62 3.17776 0.0436451 0.178572 0.0377708; 7.15 3.26011 0.045801 0.0436451 0.0401421; 7.38 3.30292 0.0284987 0.045801 0.0400468; 7.0 3.34005 0.0328766 0.0284987 0.0405198; 6.58 3.32693 0.0232238 0.0328766 0.0413081; 6.84 3.38406 0.0228948 0.0232238 0.0410688; 7.25 3.44035 0.0113077 0.0228948 0.0415791; 7.06 3.50324 0.0183375 0.0113077 0.0399681; 7.16 3.5064 0.014676 0.0183375 0.0407328; 7.61 3.54386 0.0091747 0.014676 0.0390419; 7.28 3.6063 0.00775049 0.0091747 0.0377436; 6.81 3.53619 0.011547 0.00775049 0.0363966; 6.32 3.53367 0.016315 0.011547 0.0224778; 6.13 3.5053 0.0101825 0.016315 0.0202003; 6.23 3.53123 0.037272 0.0101825 0.0172321; 6.15 3.51983 0.0111909 0.037272 0.0179632; 6.05 3.48445 0.0135623 0.0111909 0.016156; 6.31 3.48418 0.0285729 0.0135623 0.0153509; 7.86 3.56537 0.0165191 0.0285729 0.0158241; 7.82 3.60818 0.0103206 0.0165191 0.0162584; 7.64 3.63117 0.00279467 0.0103206 0.0155903; 8.96 3.6298 0.0129423 0.00279467 0.0146002; 8.91 3.65772 0.00874122 0.0129423 0.0149142; 8.66 3.6484 0.0131621 0.00874122 0.0149967; 7.74 3.62345 0.00623531 0.0131621 0.0151313; 8.0 3.63318 0.0106876 0.00623531 0.0142913; 8.59 3.66264 0.00957714 0.0106876 0.0143334; 8.19 3.64083 0.00748615 0.00957714 0.0120255; 8.12 3.65437 0.0141145 0.00748615 0.0117168; 8.13 3.6768 0.00796359 0.0141145 0.0117628; 7.91 3.67883 0.0137971 0.00796359 0.0100454; 7.45 3.63124 0.012523 0.0137971 0.00981853; 7.06 3.59675 0.0269155 0.012523 0.0100021; 6.73 3.59296 0.022324 0.0269155 0.0120121; 6.37 3.52516 0.0328174 0.022324 0.0127939; 6.27 3.50922 0.066186 0.0328174 0.0148003; 6.27 3.49457 0.0436384 0.066186 0.0192189; 6.78 3.53243 0.0168895 0.0436384 0.0223359; 6.68 3.50432 0.0147592 0.0168895 0.0228527; 6.76 3.5616 0.0142414 0.0147592 0.0232845; 6.75 3.52881 0.0120729 0.0142414 0.0238475; 6.24 3.40511 0.0353162 0.0120729 0.0236773; 6.08 3.40797 0.0502831 0.0353162 0.0259567; 6.08 3.41142 0.0279779 0.0502831 0.0289972; 6.1 3.39677 0.0149894 0.0279779 0.0302851; 6.16 3.36188 0.0129986 0.0149894 0.0292913; 5.92 3.31185 0.0167614 0.0129986 0.0285142; 5.32 3.25197 0.0207611 0.0167614 0.0271762; 4.21 3.21546 0.035183 0.0207611 0.0233908; 3.94 3.14431 0.0522645 0.035183 0.0226862; 3.54 3.02392 0.0487077 0.0522645 0.0256341; 3.21 2.90006 0.0951968 0.0487077 0.0284631; 4.05 3.04712 0.0663837 0.0951968 0.0352094; 4.04 2.99893 0.0356827 0.0663837 0.0397353; 3.72 2.96024 0.0327594 0.0356827 0.0397658; 4.53 3.0946 0.0341305 0.0327594 0.0383055; 4.36 3.13052 0.0160746 0.0341305 0.0388182; 4.55 3.14737 0.0189406 0.0160746 0.0389087; 4.87 3.19572 0.0137644 0.0189406 0.0394038; 4.96 3.23776 0.0202806 0.0137644 0.0391541; 5.16 3.28317 0.00843024 0.0202806 0.039114; 4.79 3.21636 0.00979102 0.00843024 0.0368846; 4.49 3.19312 0.0246488 0.00979102 0.0333452; 3.91 3.14332 0.0227786 0.0246488 0.0313403; 3.98 3.20656 0.017014 0.0227786 0.0253054; 4.2 3.22841 0.00925353 0.017014 0.0211913; 4.12 3.21194 0.0177494 0.00925353 0.0189888; 4.65 3.32641 0.0178187 0.0177494 0.017738; 4.73 3.32821 0.0130985 0.0178187 0.0163787; 5.04 3.34891 0.00772945 0.0130985 0.0161307; 4.93 3.35116 0.0159885 0.00772945 0.0151964; 4.85 3.30655 0.00913078 0.0159885 0.0153818; 5.12 3.33897 0.00717949 0.00913078 0.0144526; 4.97 3.3392 0.00504268 0.00717949 0.0143484; 4.74 3.30255 0.00862218 0.00504268 0.0139527; 4.59 3.31912 0.00806807 0.00862218 0.0126171; 4.66 3.29984 0.0104511 0.00806807 0.0113913; 4.57 3.24836 0.0104141 0.0104511 0.0108444; 4.8 3.3019 0.0034504 0.0104141 0.0109411; 4.54 3.26805 0.00512096 0.0034504 0.00974949; 4.36 3.21178 0.00334708 0.00512096 0.00869135; 4.23 3.19157 0.00390479 0.00334708 0.00787873; 4.18 3.19161 0.0113716 0.00390479 0.00756001; 4.58 3.14776 0.00576937 0.0113716 0.00717527; 4.68 3.18543 0.00303635 0.00576937 0.00689515; 4.7 3.17076 0.00423945 0.00303635 0.00654989; 4.63 3.12922 0.00326957 0.00423945 0.00648295; 4.34 3.11883 0.00372302 0.00326957 0.0060369; 4.03 3.07498 0.00883822 0.00372302 0.00567482; 4.23 3.08759 0.0124923 0.00883822 0.00554041; 3.89 3.08416 0.00620245 0.0124923 0.0057136; 3.62 3.02044 0.0101739 0.00620245 0.00594293; 3.54 2.98736 0.00613571 0.0101739 0.00636401; 3.68 3.00996 0.00346752 0.00613571 0.0065964; 4.1 3.08448 0.0126186 0.00346752 0.00655996; 4.12 3.08578 0.00759755 0.0126186 0.00666387; 3.79 3.06113 0.00889366 0.00759755 0.00681622; 4.04 3.11789 0.00673254 0.00889366 0.00730433; 4.27 3.1335 0.00764436 0.00673254 0.00751209; 4.08 3.1146 0.00967736 0.00764436 0.00787666; 3.66 2.98809 0.0221603 0.00967736 0.00837285; 3.67 3.0037 0.0279763 0.0221603 0.00948302; 3.46 3.00753 0.0144274 0.0279763 0.0107734; 3.54 3.04931 0.00786153 0.0144274 0.0114588; 3.37 2.99988 0.00841317 0.00786153 0.0112661; 3.51 3.05306 0.00369175 0.00841317 0.0114559; 3.55 3.055 0.00534635 0.00369175 0.0114745; 3.55 3.01038 0.00871432 0.00534635 0.0108685; 3.61 3.04535 0.00409256 0.00871432 0.0109616; 3.69 3.05586 0.00729858 0.00409256 0.0105615; 3.82 3.09307 0.00269815 0.00729858 0.0106087; 3.49 3.08286 0.00971329 0.00269815 0.0101965; 3.11 3.01164 0.0391782 0.00971329 0.0101995; 3.26 3.03684 0.0113255 0.0391782 0.0116176; 3.23 3.05318 0.00315729 0.0113255 0.0102301; 3.2 3.11622 0.00859721 0.00315729 0.00929089; 2.9 3.08902 0.0185874 0.00859721 0.0093522; 2.81 2.94537 0.0465228 0.0185874 0.0102; 3.03 2.99738 0.0258757 0.0465228 0.0137693; 3.2 3.03413 0.0154738 0.0258757 0.0154801; 3.19 3.05545 0.0106913 0.0154738 0.0160434; 3.24 3.12024 0.0113258 0.0106913 0.0165933; 3.23 3.12448 0.0112618 0.0113258 0.0169289; 3.13 3.13807 0.022712 0.0112618 0.0176425; 3.1 3.17185 0.0100923 0.022712 0.0187257; 3.23 3.23345 0.019462 0.0100923 0.0163019; 3.51 3.18895 0.0336651 0.019462 0.0169799; 3.4 3.14736 0.0141817 0.0336651 0.0195223; 3.53 3.14055 0.0148984 0.0141817 0.0199876; 3.52 3.17614 0.0107487 0.0148984 0.0196802; 3.57 3.16231 0.00869592 0.0107487 0.016699; 4.09 3.1553 0.00844733 0.00869592 0.0152674; 3.93 3.13451 0.0111374 0.00844733 0.0146818; 3.88 3.13557 0.0114623 0.0111374 0.014719; 3.72 3.04799 0.0112517 0.0114623 0.0147304; 3.54 2.97676 0.023295 0.0112517 0.0147296; 3.77 3.00593 0.0128169 0.023295 0.0147781; 3.89 3.06416 0.00964972 0.0128169 0.0150052; 3.59 3.02532 0.00462449 0.00964972 0.0141875; 3.56 2.99966 0.0246003 0.00462449 0.0117675; 3.52 2.9272 0.00919473 0.0246003 0.0126357; 3.38 2.90702 0.0134035 0.00919473 0.0121604; 3.48 2.95349 0.00758963 0.0134035 0.0123816; 3.15 2.90091 0.00714592 0.00758963 0.0122894; 3.01 2.8697 0.00952303 0.00714592 0.012181; 2.91 2.84458 0.00821384 0.00952303 0.0120464; 3.27 2.94364 0.0342597 0.00821384 0.0117757; 3.33 2.95187 0.0103447 0.0342597 0.0136931; 4.16 3.06757 0.044261 0.0103447 0.0126139; 4.49 3.10277 0.0388504 0.044261 0.0152342; 4.06 3.1185 0.0221918 0.0388504 0.0176676; 4.2 3.16569 0.0243738 0.0221918 0.0191315; 4.52 3.18055 0.0125656 0.0243738 0.0191127; 4.94 3.20116 0.00985972 0.0125656 0.0193936; 5.51 3.28033 0.00730012 0.00985972 0.0190983; 5.38 3.29387 0.0126653 0.00730012 0.0190741; 5.39 3.31238 0.0103151 0.0126653 0.0195341; 4.98 3.28944 0.0224266 0.0103151 0.0196001; 4.94 3.287 0.0104573 0.0224266 0.0207845; 4.66 3.28721 0.0107423 0.0104573 0.018801; 4.3 3.25894 0.00751638 0.0107423 0.0188341; 4.67 3.28505 0.00444809 0.00751638 0.015772; 4.53 3.24763 0.00283246 0.00444809 0.0129052; 4.58 3.23609 0.00880856 0.00283246 0.0112919; 4.33 3.2036 0.0192514 0.00880856 0.00999479; 4.24 3.20037 0.00807262 0.0192514 0.0105519; 4.11 3.19511 0.00607678 0.00807262 0.010403; 3.61 3.13225 0.00761034 0.00607678 0.0103011; 3.61 3.14333 0.0106308 0.00761034 0.00987983; 3.47 3.1236 0.00531365 0.0106308 0.00990613; 4.04 3.21632 0.0262053 0.00531365 0.00848005; 4.18 3.20479 0.00781393 0.0262053 0.00979239; 4.19 3.20406 0.00983917 0.00781393 0.00954836; 4.12 3.1823 0.00681722 0.00983917 0.00974192; 4.35 3.19097 0.0107634 0.00681722 0.00993935; 4.89 3.27162 0.00763859 0.0107634 0.0106003; 4.86 3.27631 0.00816225 0.00763859 0.0105028; 4.79 3.2586 0.00615828 0.00816225 0.00957867; 4.91 3.25602 0.00446743 0.00615828 0.00941915; 5.21 3.29826 0.00667639 0.00446743 0.00928503; 5.09 3.30691 0.00559602 0.00667639 0.0092072; 4.88 3.29817 0.0106564 0.00559602 0.00878764; 4.68 3.27826 0.00607752 0.0106564 0.00923287; 4.19 3.23006 0.0087431 0.00607752 0.00755555; 4.63 3.2754 0.0051831 0.0087431 0.00763298; 4.97 3.32 0.00711929 0.0051831 0.00724498; 5.02 3.36029 0.00934003 0.00711929 0.00727015; 5.14 3.36592 0.0139622 0.00934003 0.00715153; 5.72 3.42134 0.00909857 0.0139622 0.0076785; 6.12 3.46731 0.0119759 0.00909857 0.00775652; 5.83 3.44969 0.020886 0.0119759 0.00824133; 6.04 3.49307 0.0123068 0.020886 0.00960954; 5.82 3.502 0.0130381 0.0123068 0.0100787; 5.39 3.43079 0.0184731 0.0130381 0.0106989; 5.89 3.49017 0.00967079 0.0184731 0.0113503; 5.04 3.4105 0.0380501 0.00967079 0.0116497; 5.35 3.43992 0.00631833 0.0380501 0.014092; 5.29 3.42948 0.0150208 0.00631833 0.0141866; 5.03 3.39508 0.010825 0.0150208 0.0148451; 5.5 3.51604 0.0165197 0.010825 0.0149688; 5.5 3.56035 0.00961305 0.0165197 0.0151819; 5.3 3.57257 0.0139911 0.00961305 0.0152248; 5.29 3.53802 0.0348094 0.0139911 0.0153927; 5.28 3.53243 0.0242302 0.0348094 0.016553; 5.44 3.56934 0.0107992 0.0242302 0.0175466; 5.54 3.59675 0.00392679 0.0107992 0.0173601; 5.67 3.62055 0.0173006 0.00392679 0.0161479; 5.54 3.58289 0.0165517 0.0173006 0.0167837; 4.0 3.31877 0.449884 0.0165517 0.0149921; 3.59 3.251 0.0502355 0.449884 0.051956; 3.88 3.30317 0.0423798 0.0502355 0.0548905; 4.2 3.33764 0.0562727 0.0423798 0.0575201; 4.3 3.34078 0.0102739 0.0562727 0.0608328; 3.98 3.31169 0.0123757 0.0102739 0.0608879; 4.08 3.33196 0.0220608 0.0123757 0.0607533; 4.31 3.31787 0.0175004 0.0220608 0.0596909; 4.63 3.35675 0.0171355 0.0175004 0.0591301; 4.66 3.32944 0.00673564 0.0171355 0.0596581; 4.62 3.28449 0.00640583 0.00673564 0.0598922; 4.71 3.31424 0.00657673 0.00640583 0.0589843; 4.68 3.31156 0.0080261 0.00657673 0.058153; 4.63 3.2675 0.00631955 0.0080261 0.0213315; 4.54 3.26114 0.00396134 0.00631955 0.0176719; 4.94 3.29906 0.00368396 0.00396134 0.0144703; 4.79 3.3027 0.00857479 0.00368396 0.0100879; 5.07 3.3072 0.00594052 0.00857479 0.00994636; 5.12 3.34335 0.00445941 0.00594052 0.00941009; 5.33 3.3501 0.00628534 0.00445941 0.00794331; 5.08 3.33707 0.00976656 0.00628534 0.00700872; 5.68 3.41238 0.00224568 0.00976656 0.00639465; 5.67 3.41702 0.00368053 0.00224568 0.00602049; 5.74 3.39544 0.00247236 0.00368053 0.00579338; 5.34 3.38719 0.0301903 0.00247236 0.00545135; 5.27 3.40847 0.00342608 0.0301903 0.00729836; 5.26 3.42122 0.0183523 0.00342608 0.00705724; 4.7 3.36182 0.0111556 0.0183523 0.00825649; 4.84 3.36178 0.0069524 0.0111556 0.00887913; 4.69 3.38359 0.00729703 0.0069524 0.00874393; 4.51 3.33798 0.00803578 0.00729703 0.00885697; 4.98 3.42249 0.0076521 0.00803578 0.009155; 4.84 3.40811 0.012491 0.0076521 0.0092689; 4.51 3.40122 0.00645773 0.012491 0.00949594; 4.03 3.30269 0.0401628 0.00645773 0.00984694; 3.5 3.25093 0.0329481 0.0401628 0.0128871; 3.28 3.23414 0.0486498 0.0329481 0.0154268; 3.71 3.28873 0.0311593 0.0486498 0.0169651; 3.94 3.30842 0.0305317 0.0311593 0.0192762; 4.17 3.35548 0.0386024 0.0305317 0.0202911; 4.64 3.4199 0.0522188 0.0386024 0.0225784; 4.72 3.44391 0.0258859 0.0522188 0.0263506; 4.83 3.453 0.0285783 0.0258859 0.0278996; 5.03 3.4934 0.0154173 0.0285783 0.0296115; 4.57 3.44752 0.00780699 0.0154173 0.0302586; 4.95 3.50342 0.00726545 0.00780699 0.0298683; 5.2 3.52335 0.027933 0.00726545 0.0299356; 5.09 3.50929 0.00749819 0.027933 0.0289164; 5.0 3.53384 0.0276792 0.00749819 0.0267956; 4.51 3.48732 0.0265612 0.0276792 0.025048; 4.94 3.58907 0.0251771 0.0265612 0.0246649; 5.32 3.6124 0.0091169 0.0222993 0.0228601; 5.64 3.5815 0.00673357 0.0091169 0.0192682; 5.72 3.59794 0.0283169 0.00673357 0.0176722; 5.84 3.60186 0.008784 0.0283169 0.0176504; 5.85 3.56991 0.0107662 0.008784 0.0170977; 5.84 3.61464 0.00845005 0.0107662 0.0173443; 5.47 3.59109 0.00552929 0.00845005 0.017443; 5.51 3.59216 0.00645006 0.00552929 0.015576; 5.63 3.60315 0.00780406 0.00645006 0.0154887; 6.08 3.644 0.00636316 0.00780406 0.0138324; 6.68 3.6524 0.00330784 0.00636316 0.0121492; 6.86 3.66341 0.00908019 0.00330784 0.0103268; 7.05 3.66873 0.0252781 0.00908019 0.00922519; 7.03 3.68041 0.00464459 0.0252781 0.010572; 6.72 3.65512 0.0245012 0.00464459 0.0103979; 6.79 3.68486 0.0128932 0.0245012 0.0100799; 7.2 3.67772 0.00897489 0.0128932 0.0104223; 7.36 3.67424 0.00474715 0.00897489 0.0102731; 7.64 3.71526 0.00281754 0.00474715 0.00996449; 7.76 3.70474 0.0115949 0.00281754 0.00973851; 7.28 3.71766 0.0152335 0.0115949 0.0101672; 7.16 3.70428 0.0193429 0.0152335 0.0107864; 7.29 3.7247 0.00548308 0.0193429 0.011868; 7.81 3.75003 0.0119648 0.00548308 0.0120493; 7.3 3.72102 0.0143266 0.0119648 0.0122897; 6.45 3.67827 0.0154827 0.0143266 0.011377; 6.62 3.66743 0.0165288 0.0154827 0.0122802; 6.86 3.67639 0.00836789 0.0165288 0.0116158; 6.58 3.64312 0.00915724 0.00836789 0.0112387; 6.78 3.67028 0.00431793 0.00915724 0.0112539; 6.83 3.71263 0.0045763 0.00431793 0.0112182; 6.39 3.68896 0.00647471 0.0045763 0.0113647; 6.53 3.69806 0.0131606 0.00647471 0.010938; 6.14 3.65788 0.0107468 0.0131606 0.0107653; 6.08 3.65094 0.0147703 0.0107468 0.0100489; 6.32 3.66313 0.0110439 0.0147703 0.0108229; 6.7 3.70241 0.0158727 0.0110439 0.0107461; 6.37 3.7217 0.0067033 0.0158727 0.010875; 6.66 3.73835 0.00259856 0.0067033 0.0101434; 7.13 3.7672 0.0106309 0.00259856 0.00898252; 6.99 3.80036 0.0165413 0.0106309 0.0091711; 7.39 3.80775 0.00958021 0.0165413 0.00978644; 7.69 3.80609 0.0047203 0.00958021 0.010225; 8.1 3.83439 0.00852898 0.0047203 0.010237; 7.43 3.82408 0.0145607 0.00852898 0.0104081; 7.81 3.85115 0.0100367 0.0145607 0.0105248; 7.64 3.86011 0.0107774 0.0100367 0.0104656; 8.11 3.88359 0.0144297 0.0107774 0.0101329; 8.33 3.89565 0.00860898 0.0144297 0.0104151; 8.36 3.89407 0.0191316 0.00860898 0.00980974; 8.27 3.92599 0.015852 0.0191316 0.0108454; 8.4 3.93866 0.00817759 0.015852 0.0119499; 8.13 3.91944 0.0074614 0.00817759 0.0117454; 8.13 3.87432 0.029429 0.0074614 0.0109888; 8.36 3.90538 0.0122298 0.029429 0.0126429; 8.7 3.94505 0.00574302 0.0122298 0.0132686; 9.22 3.9589 0.00840977 0.00574302 0.0130365; 9.82 4.01297 0.0074137 0.00840977 0.0125239; 9.45 3.99599 0.0122658 0.0074137 0.0123053; 10.27 4.05142 0.0111815 0.0122658 0.0124294; 10.52 4.04107 0.0226492 0.0111815 0.0121587; 9.38 3.97932 0.0257261 0.0226492 0.0133287; 9.99 4.02011 0.0367676 0.0257261 0.0138782; 10.04 4.08863 0.0136816 0.0367676 0.0156212; 10.32 4.12882 0.0139546 0.0136816 0.0160799; 11.64 4.19258 0.0189652 0.0139546 0.016621; 10.86 4.1522 0.021 0.0189652 0.015749; 11.53 4.19955 0.0271049 0.021 0.0164798; 11.2 4.16416 0.130715 0.0271049 0.01826; 11.86 4.18913 0.067412 0.130715 0.0284521; 11.76 4.19618 0.0382175 0.067412 0.033452; 11.74 4.19439 0.0646634 0.0382175 0.0356146; 12.72 4.25714 0.0166831 0.0646634 0.0400714; 12.8 4.30706 0.0290318 0.0166831 0.0395743; 12.92 4.3061 0.0786414 0.0290318 0.0398497; 12.94 4.27883 0.0286844 0.0786414 0.0433392; 13.15 4.30207 0.0458306 0.0286844 0.0445895; 13.4 4.28265 0.0285645 0.0458306 0.0472458; 9.97 4.1 0.189989 0.0285645 0.0480457; 8.59 4.15447 0.28989 0.189989 0.0621281; 11.96 4.2269 0.275574 0.28989 0.0840269; 11.8 4.28339 0.129115 0.275574 0.0960984; 11.6 4.34213 0.0997603 0.129115 0.10124; 11.95 4.37331 0.1358 0.0997603 0.106369; 12.19 4.33384 0.131824 0.1358 0.112297; 12.74 4.36179 0.0868624 0.131824 0.121892; 13.22 4.41355 0.262455 0.0868624 0.126711; 12.41 4.34772 0.180371 0.262455 0.142029; 12.97 4.39403 0.2056 0.180371 0.154669; 12.56 4.35731 0.0415169 0.2056 0.167984; 11.82 4.34541 0.112041 0.0415169 0.169063; 11.15 4.3134 0.060023 0.112041 0.162567; 12.22 4.37393 0.182118 0.060023 0.143412; 12.05 4.40506 0.142499 0.182118 0.135624; 11.83 4.48345 0.0717098 0.142499 0.136739; 10.78 4.44236 0.129219 0.0717098 0.134401; 10.24 4.46961 0.0791773 0.129219 0.133853; 11.07 4.52347 0.199013 0.0791773 0.129466; 10.52 4.44625 0.47681 0.199013 0.138812; 10.86 4.45183 0.181941 0.47681 0.156675; 10.59 4.50674 0.154481 0.181941 0.156806; 11.62 4.48882 0.120288 0.154481 0.152546; 12.94 4.56382 0.0823393 0.120288 0.15911; 12.06 4.52088 0.127238 0.0823393 0.156635; 12.04 4.49382 0.279058 0.127238 0.162236; 10.91 4.39027 0.108493 0.279058 0.170314; 11.4 4.40858 0.289244 0.108493 0.167481; 11.64 4.4445 0.368885 0.289244 0.185608; 10.89 4.34369 0.139849 0.368885 0.205581; 10.67 4.2863 0.265159 0.139849 0.210637; 11.09 4.37095 0.457359 0.265159 0.216149; 11.4 4.37866 0.0964968 0.457359 0.214528; 11.0 4.36738 0.0744541 0.0964968 0.207407; 10.85 4.34446 0.07582 0.0744541 0.200738; 10.41 4.28599 0.0691933 0.07582 0.197033; 8.98 4.15481 0.371706 0.0691933 0.195937; 8.65 4.18935 0.131877 0.371706 0.21631; 9.29 4.2647 0.0804986 0.131877 0.204045; 9.01 4.27409 0.0720876 0.0804986 0.201712; 8.99 4.25772 0.101816 0.0720876 0.183615; 8.69 4.22984 0.115279 0.101816 0.16136; 9.39 4.28624 0.0686651 0.115279 0.159312; 9.22 4.24027 0.0564201 0.0686651 0.142938; 9.77 4.22566 0.0918655 0.0564201 0.109526; 8.94 4.13781 0.0961874 0.0918655 0.10914; 8.16 4.05331 0.204122 0.0961874 0.110951; 8.37 4.06602 0.164481 0.204122 0.121643; 7.23 3.97475 0.159079 0.164481 0.129584; 7.6 4.03632 0.350299 0.159079 0.111865; 8.19 4.09455 0.147462 0.350299 0.130067; 7.57 4.03106 0.0721512 0.147462 0.135647; 7.11 4.00721 0.122622 0.0721512 0.135652; 7.05 3.98658 0.0673919 0.122622 0.137386; 6.85 3.98509 0.15058 0.0673919 0.133395; 7.66 4.0415 0.0671606 0.15058 0.140222; 7.98 4.10112 0.0601661 0.0671606 0.141117; 8.02 4.10536 0.0712229 0.0601661 0.138475; 8.37 4.11672 0.05785 0.0712229 0.136395; 8.38 4.13307 0.0305843 0.05785 0.124205; 8.35 4.11118 0.0578079 0.0305843 0.113047; 8.93 4.14733 0.0378894 0.0578079 0.104608; 8.88 4.1475 0.0400608 0.0378894 0.0785739; 9.15 4.17566 0.0187112 0.0400608 0.0696239; 8.78 4.19044 0.0315006 0.0187112 0.0651705; 8.99 4.19352 0.0446757 0.0315006 0.0575771; 8.72 4.16501 0.0435281 0.0446757 0.0556841; 8.71 4.13255 0.035298 0.0435281 0.0467631; 8.63 4.11758 0.0243698 0.035298 0.0441079; 8.6 4.12153 0.0275137 0.0243698 0.0411249; 8.57 4.07439 0.0348378 0.0275137 0.0374825; 8.76 4.06591 0.0350053 0.0348378 0.0355648; 8.61 4.07835 0.0408418 0.0350053 0.0359332; 8.8 4.08572 0.0298743 0.0408418 0.0345194; 8.83 4.00317 0.0177237 0.0298743 0.0338514; 8.95 4.02169 0.0139754 0.0177237 0.03199; 7.68 3.99784 0.022127 0.0139754 0.0315953; 7.72 3.9955 0.0194192 0.022127 0.0308142; 7.48 3.95248 0.0235278 0.0194192 0.0287095; 7.43 3.92234 0.0387397 0.0235278 0.0270428; 7.49 3.94802 0.0377844 0.0387397 0.0273296; 7.61 3.95308 0.010477 0.0377844 0.0284475; 7.75 3.98572 0.0105807 0.010477 0.0270278; 7.65 3.96797 0.0148468 0.0105807 0.0250063; 7.93 3.95684 0.0156236 0.0148468 0.0233265; 7.92 3.93125 0.0425062 0.0156236 0.021225; 8.23 4.05059 0.0248735 0.0425062 0.0222776; 8.33 4.04234 0.00959146 0.0248735 0.0228734; 8.09 4.06415 0.0240059 0.00959146 0.0225081; 8.49 4.0444 0.0221307 0.0240059 0.0226647; 8.46 4.06501 0.022206 0.0221307 0.0228906; 8.61 4.05974 0.0249786 0.022206 0.0227805; 8.17 4.00087 0.0718094 0.0249786 0.0216337; 7.97 3.99773 0.0984261 0.0718094 0.0244692; 7.77 3.99294 0.0585436 0.0984261 0.0317982; 7.9 3.99801 0.022539 0.0585436 0.0357951; 8.14 4.01733 0.0219371 0.022539 0.0364362; 8.29 4.05283 0.0181744 0.0219371 0.0369623; 8.07 4.05087 0.0450168 0.0181744 0.0349346; 8.18 4.04149 0.017293 0.0450168 0.0366132; 7.6 4.05922 0.031752 0.017293 0.037255; 7.28 4.0272 0.0537384 0.031752 0.0379005; 7.03 4.0271 0.0691728 0.0537384 0.0405345; 7.34 4.06588 0.0216218 0.0691728 0.0444484; 7.46 4.06005 0.0215354 0.0216218 0.0441687; 7.04 4.04495 0.0356381 0.0215354 0.0399792; 6.73 4.00501 0.0663055 0.0356381 0.0347469; 6.41 4.00674 0.152702 0.0663055 0.0353937; 6.44 4.02339 0.0734082 0.152702 0.0462406; 6.53 4.03715 0.068253 0.0734082 0.0505299; 6.05 3.99065 0.203762 0.068253 0.0547031; 5.68 3.95808 0.0932746 0.203762 0.0679318; 5.16 3.88403 0.182583 0.0932746 0.0742636; 4.77 3.86691 0.121 0.182583 0.0868329; 4.46 3.84246 0.552488 0.121 0.092438; 4.7 3.89913 0.134563 0.552488 0.132714; 4.47 3.92449 0.0818618 0.134563 0.142126; 3.97 3.82291 0.159035 0.0818618 0.147153; 4.35 3.80592 0.423966 0.159035 0.157436; 4.38 3.81696 0.160222 0.423966 0.187241; 4.47 3.71217 1.03051 0.160222 0.187868; 3.57 3.51219 1.19969 1.03051 0.267626; 2.77 3.42492 1.22939 1.19969 0.361912; 2.6 3.46834 0.80178 1.22939 0.447381; 2.49 3.39047 0.709717 0.80178 0.506423; 2.23 3.28251 0.441693 0.709717 0.550351; 2.65 3.3649 0.859499 0.441693 0.577076; 3.4 3.47591 0.363713 0.859499 0.60266; 3.85 3.56191 0.252403 0.363713 0.621756; 4.01 3.5967 0.0898218 0.252403 0.635967; 4.6 3.69113 0.106431 0.0898218 0.6302; 5.04 3.73584 0.0508081 0.106431 0.603738; 5.2 3.80446 0.0644148 0.0508081 0.594621; 4.88 3.79134 0.092285 0.0644148 0.514113; 4.95 3.8476 0.0530035 0.092285 0.421829; 5.05 3.89209 0.0227195 0.0530035 0.323797; 4.67 3.85633 0.0640625 0.0227195 0.258876; 5.07 3.90217 0.055555 0.0640625 0.205071; 5.44 3.98438 0.0160577 0.055555 0.172893; 5.35 4.00866 0.0840481 0.0160577 0.102606; 4.9 3.91423 0.146865 0.0840481 0.0793008; 4.6 3.855 0.079108 0.146865 0.070506; 5.08 3.89936 0.0414568 0.079108 0.0696132; 4.75 3.84393 0.0480049 0.0414568 0.0641987; 4.8 3.91494 0.0637759 0.0480049 0.0639651; 4.82 3.93236 0.0233433 0.0637759 0.0639118; 4.66 3.93095 0.0454893 0.0233433 0.0581667; 4.95 3.96794 0.0226372 0.0454893 0.0575405; 5.04 3.98323 0.0359067 0.0226372 0.0575336; 5.34 4.00765 0.0313676 0.0359067 0.0551873; 5.07 4.00176 0.037243 0.0313676 0.0531717; 5.06 4.01534 0.0157077 0.037243 0.0549371; 4.99 3.98971 0.0221936 0.0157077 0.0492421; 4.8 3.95249 0.0482368 0.0221936 0.0388528; 4.51 3.9315 0.0643981 0.0482368 0.0362802; 4.02 3.85815 0.381819 0.0643981 0.038192; 3.54 3.76478 0.203012 0.381819 0.0660099; 3.87 3.87409 0.329498 0.203012 0.0776129; 3.52 3.84395 0.243452 0.329498 0.103126; 3.63 3.83922 0.104327 0.243452 0.119623; 4.08 3.88328 0.0429844 0.104327 0.12643; 4.36 3.90645 0.0305972 0.0429844 0.12702; 4.56 3.91139 0.0574196 0.0305972 0.126956; 4.3 3.89668 0.0457753 0.0574196 0.128637; 3.79 3.81411 0.0263635 0.0457753 0.131143; 3.98 3.83369 0.0787785 0.0263635 0.13149; 3.9 3.83949 0.0417438 0.0787785 0.134035; 4.08 3.84185 0.0246239 0.0417438 0.132147; 4.33 3.84118 0.0386762 0.0246239 0.102381; 4.55 3.81186 0.021289 0.0386762 0.0886866; 4.63 3.77378 0.0395661 0.021289 0.0630025; 4.9 3.71841 0.0191646 0.0395661 0.046012; 4.91029 3.76446 0.0223779 0.0191646 0.0389152; 4.8286 3.75276 0.0459764 0.0223779 0.037198; 4.80042 3.77553 0.0181368 0.0459764 0.0384796; 5.06895 3.80082 0.0473042 0.0181368 0.035206; 5.31549 3.78916 0.0360042 0.0473042 0.0353334; 5.12041 3.75204 0.0720758 0.0360042 0.0361368; 5.57379 3.80159 0.0151387 0.0720758 0.0355783; 5.3427 3.75831 0.0310281 0.0151387 0.0333612; 5.54734 3.7994 0.0175181 0.0310281 0.0338948; 5.65045 3.84005 0.0270681 0.0175181 0.0321317; 5.86609 3.88612 0.0194642 0.0270681 0.0326133; 6.02478 3.94835 0.0213394 0.0194642 0.0309381; 5.82015 3.90927 0.0370949 0.0213394 0.0311193; 6.06535 3.93439 0.0392782 0.0370949 0.0323457; 5.98842 3.92786 0.0266122 0.0392782 0.0317876; 5.86352 3.9153 0.0490829 0.0266122 0.0324938; 5.88856 3.93099 0.0242753 0.0490829 0.0326421; 5.86827 3.94484 0.0125888 0.0242753 0.0316647; 5.90605 3.91651 0.0214852 0.0125888 0.0267074; 6.0181 3.94649 0.0125183 0.0214852 0.0272363; 6.06792 3.91148 0.0174854 0.0125183 0.0256938; 6.07831 3.92527 0.0473901 0.0174854 0.0256911; 6.08681 3.92808 0.00232388 0.0473901 0.0273846; 6.16194 3.91007 0.0405431 0.00232388 0.0259562; 5.49262 3.87884 0.0505182 0.0405431 0.0275565; 5.92662 3.91669 0.0196693 0.0505182 0.0286751; 5.85899 3.8938 0.0334556 0.0196693 0.0270411; 6.20164 3.89778 0.00868597 0.0334556 0.0276114; 6.23238 3.89925 0.0129965 0.00868597 0.0242449; 6.47019 3.87353 0.0220447 0.0129965 0.023305; 6.62655 3.88593 0.0188653 0.0220447 0.024093; 6.13074 3.8132 0.079256 0.0188653 0.0238747; 5.82759 3.76922 0.0750067 0.079256 0.0294362; 6.13214 3.83488 0.0460227 0.0750067 0.0342296; 6.14162 3.8313 0.0254232 0.0460227 0.0341157; 6.27796 3.80376 0.0716164 0.0254232 0.0360406; 5.51006 3.7403 0.111635 0.0716164 0.03863; 5.13294 3.74316 0.084967 0.111635 0.0437231];

y = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];

offset = [3.17805, 3.49651, 3.68888, 3.68888, 3.73767, 3.29584, 3.52636, 3.49651, 3.46574, 3.46574, 3.66356, 3.71357, 3.68888, 3.46574, 3.43399, 3.49651, 3.4012, 3.49651, 3.63759, 3.43399, 3.49651, 3.61092, 3.46574, 3.68888, 3.73767, 3.4012, 3.52636, 3.43399, 3.4012, 3.46574, 3.66356, 3.71357, 3.63759, 3.43399, 3.66356, 3.43399, 3.49651, 3.63759, 3.68888, 3.63759, 3.7612, 3.55535, 3.52636, 3.71357, 3.58352, 3.68888, 3.55535, 3.49651, 3.71357, 3.55535, 3.52636, 3.55535, 3.82864, 3.29584, 3.58352, 3.55535, 3.66356, 3.58352, 3.58352, 3.43399, 3.71357, 3.21888, 3.52636, 3.52636, 3.46574, 3.43399, 3.49651, 3.55535, 3.71357, 3.7612, 3.43399, 3.52636, 3.43399, 3.63759, 3.95124, 3.66356, 3.49651, 3.66356, 3.58352, 3.55535, 3.71357, 3.68888, 3.4012, 3.78419, 3.21888, 3.29584, 3.63759, 3.3673, 3.55535, 3.66356, 3.46574, 3.68888, 3.68888, 3.55535, 3.63759, 3.66356, 3.7612, 3.63759, 3.66356, 3.7612, 3.52636, 3.8712, 3.4012, 3.55535, 3.58352, 3.68888, 3.95124, 3.52636, 3.55535, 3.3322, 3.61092, 3.52636, 3.66356, 3.98898, 3.8712, 3.55535, 3.49651, 3.98898, 3.7612, 3.3673, 3.68888, 3.3673, 3.7612, 3.52636, 3.55535, 3.3322, 3.58352, 3.3673, 3.52636, 3.68888, 3.4012, 3.66356, 3.43399, 3.52636, 3.55535, 3.71357, 3.68888, 3.58352, 3.63759, 3.71357, 3.58352, 3.78419, 3.55535, 3.71357, 3.52636, 3.4012, 3.3322, 3.52636, 3.55535, 3.63759, 3.71357, 3.71357, 3.61092, 3.52636, 3.71357, 3.68888, 3.58352, 3.49651, 3.68888, 3.49651, 3.49651, 3.61092, 3.52636, 3.43399, 3.4012, 3.21888, 3.17805, 3.58352, 3.7612, 3.63759, 3.63759, 3.55535, 3.68888, 3.68888, 3.43399, 3.66356, 3.3322, 3.52636, 3.73767, 3.68888, 3.82864, 3.61092, 3.68888, 3.89182, 3.55535, 3.52636, 3.85015, 3.61092, 3.78419, 3.71357, 3.49651, 3.82864, 3.58352, 3.4012, 3.43399, 3.55535, 3.66356, 3.66356, 3.71357, 3.66356, 3.63759, 3.68888, 3.55535, 3.68888, 3.29584, 3.61092, 3.63759, 3.49651, 3.52636, 3.66356, 3.3673, 3.43399, 3.3673, 4.09434, 3.78419, 3.63759, 3.3673, 3.43399, 3.3322, 3.43399, 3.2581, 3.46574, 3.09104, 3.43399, 3.3322, 3.29584, 3.3322, 3.46574, 3.85015, 3.68888, 3.71357, 3.46574, 3.73767, 3.97029, 3.82864, 3.91202, 3.91202, 4.14313, 3.7612, 3.55535, 4.09434, 4.02535, 3.63759, 3.52636, 3.89182, 3.63759, 3.71357, 4.17439, 3.93183, 3.93183, 3.80666, 3.63759, 3.21888, 3.55535, 3.17805, 3.17805, 3.80666, 3.21888, 3.3673, 3.63759, 3.2581, 3.49651, 3.29584, 2.89037, 3.17805, 3.49651, 3.49651, 3.55535, 3.3322, 3.46574, 3.2581, 3.2581, 3.43399, 3.17805, 3.17805, 5.61677, 5.42935, 5.70711, 5.69036, 5.50533, 5.7301, 5.743, 5.49306, 5.58725, 5.53733, 5.53733, 5.53339, 5.69709, 5.43808, 5.66988, 5.43808, 5.65599, 6.00635, 5.96615, 5.86647, 5.90536, 5.96101, 5.85793, 6.01127, 5.90263, 5.76205, 5.98394, 5.9162, 5.53339, 6.01616, 5.63121, 5.87212, 5.33754, 5.84644, 5.90536, 5.70378, 5.9108, 5.59471, 5.68017, 5.50939, 5.73979, 5.94803, 5.91889, 6.17794, 5.98141, 6.10702, 5.81711, 5.66296, 5.8999, 5.88053, 5.95584, 5.92158, 5.93754, 6.01127, 5.8522, 5.273, 5.22575, 5.42935, 4.76217, 5.0876, 5.0876, 5.1299, 5.0876, 5.22575, 5.17615, 5.11199, 5.19296, 5.24175, 5.33272, 5.28827, 5.11199, 5.17615, 5.39816, 5.03044, 5.1299, 5.0876, 4.90527, 4.99043, 5.09375, 5.0814, 4.92725, 5.09987, 4.82028, 4.96981, 5.19296, 4.96284, 5.47646, 5.15329, 5.04986, 5.25227, 5.15906, 5.14166, 5.09987, 5.51343, 5.12396, 5.37528, 5.35186, 5.40268, 5.45104, 5.39363, 5.48894, 5.1985, 5.42053, 5.48894, 5.41165, 5.54126, 5.3799, 5.26786, 5.5835, 5.41165, 5.40717, 5.59471, 5.59099, 5.45532, 5.66643, 5.57595, 5.56068, 5.75257, 5.53733, 5.70044, 5.78383, 5.49306, 5.56452, 5.50126, 5.80212, 5.97635, 5.6058, 5.2832, 5.50126, 5.55296, 5.273, 5.54518, 4.95583, 4.8752, 5.24175, 5.02388, 5.14166, 5.17048, 5.14166, 5.20401, 4.91998, 5.38907, 5.27811, 5.20401, 5.18739, 4.8752, 5.0689, 5.18739, 5.1985, 5.0876, 5.21494, 5.03695, 5.273, 5.35659, 5.15329, 5.11799, 5.12396, 4.99721, 5.14166, 4.97673, 5.0689, 5.22036, 5.16479, 5.29832, 5.42053, 5.273, 4.96284, 5.10595, 4.92725, 4.90527, 5.19296, 4.98361, 5.09987, 5.29832, 4.7185, 5.31321, 5.24175, 5.273, 4.85203, 4.54329, 4.55388, 4.52179, 4.44265, 4.41884, 4.47734, 4.5326, 4.57471, 4.20469, 4.59512, 4.7362, 4.60517, 4.78749, 4.55388, 4.41884, 4.38203, 4.47734, 4.38203, 4.20469, 4.39445, 4.44265, 3.80666, 4.34381, 4.06044, 4.12713, 4.84419, 4.60517, 4.82028, 4.81218, 4.38203, 4.70953, 4.63473, 4.51086, 4.41884, 4.64439, 4.66344, 4.61512, 4.52179, 4.21951, 4.40672, 4.57471, 4.65396, 4.49981, 4.44265, 4.51086, 4.80402, 4.48864, 4.40672, 4.51086, 4.45435, 4.57471, 4.58497, 4.49981, 4.64439, 4.47734, 4.47734, 4.7185, 4.7185, 4.79579, 4.56435, 4.78749, 4.85981, 4.68213, 4.60517, 4.60517, 4.54329, 4.79579, 5.01064, 4.74493, 4.82028, 5.03695, 4.67283, 4.61512, 4.70048, 4.44265, 4.20469, 4.54329, 4.40672, 4.68213, 4.48864, 4.11087, 4.62497, 4.65396, 4.75359, 4.7362, 4.56435, 4.54329, 4.77068, 4.67283, 4.34381, 4.94164, 4.59512, 4.72739, 4.65396, 4.70048, 4.39445, 4.74493, 4.7362, 4.56435];


text:
@testset "UI" begin
@testset "ArrayOp" begin
    ctx = Context()
    log = Dagger.LocalEventLog()
    ctx.log_sink = log

    sz = 2^4
    bsz = 2^2
    X = rand(Float32, sz, sz)
    XD = Distribute(Blocks(bsz, bsz), X)

    dag = XD*XD
    collect(ctx, dag)
    logs = Dagger.get_logs!(log)
    @test logs isa Vector{Dagger.Timespan}
    @test length(logs) > 0
    plan = Dagger.show_plan(logs)
    @test plan isa String
    @test occursin("digraph {", plan)
    @test occursin("Comm:", plan)
    @test occursin("Move:", plan)
    @test occursin("Compute:", plan)
    @test endswith(plan, "}\n")
end

@testset "Argument Merging" begin
    ctx = Context()
    log = Dagger.LocalEventLog()
    ctx.log_sink = log

    X = rand(Float32, 4, 3)
    a = delayed(sum)(X)
    b = delayed(sum)(X)
    c = delayed(+)(a,b)

    X1 = rand(Float32, 5, 10)
    X2 = rand(Float32, 5, 10)
    f = delayed(sum)(X1)
    g = delayed(sum)(X2)
    h = delayed(+)(f,g)

    j = delayed(+)(c,h)
    collect(ctx, j)
    logs = Dagger.get_logs!(log)
    plan = Dagger.show_plan(logs, j)
end
end


text:
__precompile__()

module Pardiso

if !isfile(joinpath(@__DIR__, "..", "deps", "deps.jl"))
    error("""please run Pkg.build("Pardiso") before loading the package""")
end

include("../deps/deps.jl")

function show_build_log()
    logfile = joinpath(@__DIR__, "..", "deps", "build.log")
    if !isfile(logfile)
        error("no log file found")
    else
        println(read(logfile, String))
    end
end

using Libdl
using SparseArrays
using LinearAlgebra
import Base.show

if !LOCAL_MKL_FOUND
    import MKL_jll
end

if LinearAlgebra.BLAS.vendor() === :mkl && LinearAlgebra.BlasInt == Int64
    const MklInt = Int64
    const PARDISO_FUNC = :pardiso_64
else
    const MklInt = Int32
    const PARDISO_FUNC = :pardiso
end

if Sys.iswindows()
    const libmkl_rt = "mkl_rt"
elseif Sys.isapple()
    const libmkl_rt = "@rpath/libmkl_rt.dylib"
else
    const libmkl_rt = "libmkl_rt"
end

export PardisoSolver, MKLPardisoSolver
export set_iparm!, set_dparm!, set_matrixtype!, set_solver!, set_phase!, set_msglvl!, set_nprocs!
export get_iparm, get_iparms, get_dparm, get_dparms
export get_matrixtype, get_solver, get_phase, get_msglvl, get_nprocs
export set_maxfct!, set_perm!, set_mnum!
export get_maxfct, get_perm, get_mnum
export checkmatrix, checkvec, printstats, pardisoinit, pardiso
export solve, solve!
export get_matrix
export schur_complement, pardisogetschur
export fix_iparm!

struct PardisoException <: Exception
    info::String
end

struct PardisoPosDefException <: Exception
    info::String
end

Base.showerror(io::IO, e::Union{PardisoException,PardisoPosDefException}) = print(io, e.info);

const PardisoNumTypes = Union{Float64,ComplexF64}

abstract type AbstractPardisoSolver end

function collect_gfortran_lib_candidates(main_ver)
    candidates = String[]
    if Sys.isapple()
        homebrew_gcc = "/usr/local/Cellar/gcc/"
        isdir(homebrew_gcc) || return String[]
        vers = readdir(homebrew_gcc)
        filter!(x -> startswith(x, "$main_ver."), vers)
        for v in vers
            push!(candidates, joinpath(homebrew_gcc, v, "lib/gcc/$main_ver/"))
        end
    elseif Sys.islinux()
        gcc_path = "/usr/lib/gcc/x86_64-linux-gnu/"
        isdir(gcc_path) || return String[]
        vers = readdir(gcc_path)
        filter!(x -> startswith(x, "$main_ver.") || isequal(x, "$main_ver"), vers)
        for v in vers
            push!(candidates, joinpath(gcc_path, v))
        end
    end
    return candidates
end

load_lib_fortran(lib::String, v::Int) = load_lib_fortran(lib, [v])
function load_lib_fortran(lib::String, vs::Vector{Int})
    candidates = String[]
    for v in vs
        append!(candidates, collect_gfortran_lib_candidates(v))
    end
    path = Libdl.find_library(lib, candidates)
    isempty(path) && (path = lib)
    Libdl.dlopen(path * "." * Libdl.dlext, Libdl.RTLD_GLOBAL)
end

# Pardiso
const init = Ref{Ptr}()
const pardiso_f = Ref{Ptr}()
const pardiso_chkmatrix = Ref{Ptr}()
const pardiso_chkmatrix_z = Ref{Ptr}()
const pardiso_printstats = Ref{Ptr}()
const pardiso_printstats_z = Ref{Ptr}()
const pardiso_chkvec = Ref{Ptr}()
const pardiso_chkvec_z = Ref{Ptr}()
const pardiso_get_schur_f = Ref{Ptr}()
const PARDISO_LOADED = Ref(false)

function __init__()
    if !haskey(ENV, "PARDISOLICMESSAGE")
        ENV["PARDISOLICMESSAGE"] = 1
    end

    if LOCAL_MKL_FOUND && !haskey(ENV, "MKLROOT")
        @warn "MKLROOT not set, MKL Pardiso solver will not be functional"
    end

    # This is apparently needed for MKL to not get stuck on 1 thread when
    # libpardiso is loaded in the block below...
    get_nprocs_mkl()

    if PARDISO_LIB_FOUND
        try
            libpardiso = Libdl.dlopen(PARDISO_PATH)
            init[] = Libdl.dlsym(libpardiso, "pardisoinit")
            pardiso_f[] = Libdl.dlsym(libpardiso, "pardiso")
            pardiso_chkmatrix[] = Libdl.dlsym(libpardiso, "pardiso_chkmatrix")
            pardiso_chkmatrix_z[] = Libdl.dlsym(libpardiso, "pardiso_chkmatrix_z")
            pardiso_printstats[] = Libdl.dlsym(libpardiso, "pardiso_printstats")
            pardiso_printstats_z[] = Libdl.dlsym(libpardiso, "pardiso_printstats_z")
            pardiso_chkvec[] = Libdl.dlsym(libpardiso, "pardiso_chkvec")
            pardiso_chkvec_z[] = Libdl.dlsym(libpardiso, "pardiso_chkvec_z")
            pardiso_get_schur_f[] = Libdl.dlsym(libpardiso, "pardiso_get_schur")

            if Sys.isunix()
                gfortran_v = [8, 9]
                for lib in ("libgfortran", "libgomp")
                    load_lib_fortran(lib, gfortran_v)
                end
            end

            # Windows Pardiso lib comes with BLAS + LAPACK prebaked but not on UNIX so we open them here
            # if not MKL is loaded
            if Sys.isunix()
                ptr = C_NULL
                for l in ("libblas", "libblas.so.3")
                    ptr = Libdl.dlopen_e(l, Libdl.RTLD_GLOBAL)
                    if ptr !== C_NULL
                        break
                    end
                end
                if ptr == C_NULL
                    error("could not load blas library")
                end
            end
            PARDISO_LOADED[] = true
        catch e
            @error("Pardiso did not manage to load, error thrown was: $(sprint(showerror, e))")
        end
    end
end
include("enums.jl")
include("project_pardiso.jl")
include("mkl_pardiso.jl")

# Getters and setters
set_matrixtype!(ps::AbstractPardisoSolver, v::Int) = set_matrixtype!(ps, MatrixType(v))
function set_matrixtype!(ps::AbstractPardisoSolver, v::MatrixType)
    ps.mtype = v
end

get_matrixtype(ps::AbstractPardisoSolver) = ps.mtype
get_iparm(ps::AbstractPardisoSolver, i::Integer) = ps.iparm[i]
get_iparms(ps::AbstractPardisoSolver) = ps.iparm
set_iparm!(ps::AbstractPardisoSolver, i::Integer, v::Integer) = ps.iparm[i] = v

get_mnum(ps::AbstractPardisoSolver) = ps.mnum
set_mnum!(ps::AbstractPardisoSolver, mnum::Integer) = ps.mnum = mnum

get_maxfct(ps::AbstractPardisoSolver) = ps.maxfct
set_maxfct!(ps::AbstractPardisoSolver, maxfct::Integer) = ps.maxfct = maxfct

get_perm(ps::AbstractPardisoSolver) = ps.perm
set_perm!(ps::AbstractPardisoSolver, perm::Vector{T}) where {T <: Integer} = ps.perm = convert(Vector{Int32}, perm)

get_phase(ps::AbstractPardisoSolver) = ps.phase

set_phase!(ps::AbstractPardisoSolver, v::Int) = set_phase!(ps, Phase(v))
function set_phase!(ps::AbstractPardisoSolver, v::Phase)
    ps.phase = v
end

get_msglvl(ps::AbstractPardisoSolver) = ps.msglvl

set_msglvl!(ps::AbstractPardisoSolver, v::Integer) = set_msglvl!(ps, MessageLevel(v))
function set_msglvl!(ps::AbstractPardisoSolver, v::MessageLevel)
    ps.msglvl = v
end

function pardisoinit(ps::AbstractPardisoSolver)
    ccall_pardisoinit(ps)
    return
end

function solve(ps::AbstractPardisoSolver, A::SparseMatrixCSC{Tv,Ti},
               B::StridedVecOrMat{Tv}, T::Symbol=:N) where {Ti, Tv <: PardisoNumTypes}
  X = copy(B)
  solve!(ps, X, A, B, T)
  return X
end

function fix_iparm!(ps::AbstractPardisoSolver, T::Symbol)
    # We need to set the transpose flag in PARDISO when we DON'T want
    # a transpose in Julia because we are passing a CSC formatted
    # matrix to PARDISO which expects a CSR matrix.
    if T === :N
        if isa(ps, PardisoSolver)
            set_iparm!(ps, 12, 1)
        else
            # iparm[12] = 1 is complex conjugate in MKL
            set_iparm!(ps, 12, 2)
        end
    elseif T === :C || T === :T
        set_iparm!(ps, 12, 0)
    else
        throw(ArgumentError("only :T, :N and :C, are valid transpose symbols"))
    end
end

function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
                A::SparseMatrixCSC{Tv,Ti}, B::StridedVecOrMat{Tv},
                T::Symbol=:N) where {Ti, Tv <: PardisoNumTypes}
    set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)

    # This is the heuristics for choosing what matrix type to use
    ##################################################################
    # - If hermitian try to solve with symmetric positive definite.
    #   - On pos def exception, solve instead with symmetric indefinite.
    # - If complex and symmetric, solve with symmetric complex solver
    # - Else solve as unsymmetric.
    if ishermitian(A)
        eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_POSDEF) : set_matrixtype!(ps, COMPLEX_HERM_POSDEF)
        pardisoinit(ps)
        fix_iparm!(ps, T)
        try
            pardiso(ps, X, get_matrix(ps, A, T), B)
        catch e
            set_phase!(ps, RELEASE_ALL)
            pardiso(ps, X, A, B)
            set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
            if !isa(e, PardisoPosDefException)
                rethrow()
            end
            eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_INDEF) : set_matrixtype!(ps, COMPLEX_HERM_INDEF)
            pardisoinit(ps)
            fix_iparm!(ps, T)
            pardiso(ps, X, get_matrix(ps, A, T), B)
        end
    elseif issymmetric(A)
        set_matrixtype!(ps, COMPLEX_SYM)
        pardisoinit(ps)
        fix_iparm!(ps, T)
        pardiso(ps, X, get_matrix(ps, A, T), B)
    else
        eltype(A) == Float64 ? set_matrixtype!(ps, REAL_NONSYM) : set_matrixtype!(ps, COMPLEX_NONSYM)
        pardisoinit(ps)
        fix_iparm!(ps, T)
        pardiso(ps, X, get_matrix(ps, A, T), B)
    end

    # Release memory, TODO: We are running the convert on IA and JA here
    # again which is unnecessary.
    set_phase!(ps, RELEASE_ALL)
    pardiso(ps, X, A, B)

    set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
    return X
end

function get_matrix(ps::AbstractPardisoSolver, A, T)
    mtype = get_matrixtype(ps)

    if isposornegdef(mtype)
        if ps isa MKLPardisoSolver
            T == :C && return conj(tril(A))
            return tril(A)
        elseif ps isa PardisoSolver
            T == :T && return tril(A)
            return conj(tril(A))
        end
    end

    if !issymmetric(mtype)
        T == :C && return conj(A)
        return A
    end

    if mtype == COMPLEX_SYM
        T == :C && return conj(tril(A))
        return tril(A)
    end

    error("Unhandled matrix type")
end

function pardiso(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv}, A::SparseMatrixCSC{Tv,Ti},
                 B::StridedVecOrMat{Tv}) where {Ti, Tv <: PardisoNumTypes}
    if length(X) != 0
        dim_check(X, A, B)
    end

    if Tv <: Complex && isreal(get_matrixtype(ps))
        throw(ErrorException(string("input matrix is complex while PardisoSolver ",
                                    "has a real matrix type set: $(get_matrixtype(ps))")))
    end

    if Tv <: Real && !isreal(get_matrixtype(ps))
        throw(ErrorException(string("input matrix is real while PardisoSolver ",
                                    "has a complex matrix type set: $(get_matrixtype(ps))")))
    end

    N = size(A, 2)

    AA = A.nzval
    IA = convert(Vector{Int32}, A.colptr)
    JA = convert(Vector{Int32}, A.rowval)

    resize!(ps.perm, size(B, 1))

    NRHS = size(B, 2)

    ccall_pardiso(ps, N, A.nzval, A.colptr, A.rowval, NRHS, B, X)
end

pardiso(ps::AbstractPardisoSolver) = ccall_pardiso(ps, Int32(0), Float64[], Int32[], Int32[], Int32(0), Float64[], Float64[])
function pardiso(ps::AbstractPardisoSolver, A::SparseMatrixCSC{Tv,Ti}, B::StridedVecOrMat{Tv}) where {Ti, Tv <: PardisoNumTypes}
    pardiso(ps, Tv[], A, B)
end

# populated rows of S determine schur complment block
"""
    schur_complement(ps,A,x) -> S

Schur complement `S` of the submatrix defined by the nonzero entries of `x` in matrix `A`.
If `n=nnz(x)`, then `S` is `n`-by-`n`.

WARNING: for complex `M`, seems to be unstable, made worse as number of nonzero elements in `M` decreases
"""
schur_complement(ps::AbstractPardisoSolver,A,x::SparseVector,T::Symbol=:N) = _schur_complement_permuted(ps,A,x.nzind,T)
schur_complement(ps::AbstractPardisoSolver,A,x::SparseMatrixCSC,T::Symbol=:N) = _schur_complement_permuted(ps,A,unique!(sort!(x.rowval)),T)

# permute A and then compute complement of lower right-hand `n`-by-`n` block
function _schur_complement_permuted(ps,A,rows,T::Symbol)
    P = sparse(vcat(setdiff(1:A.n,rows),rows),1:A.n,1,size(A)...)
    schur_complement(ps,P'*A*P,length(rows),T)
end
# or integer gives last n rows and columns as schur complement block
"""
    schur_complement(ps,M,n) -> S

Schur complement `S` of upper-left block in `M`, where `n` is the size of lower-right block (and therefore also of Schur complement)

WARNING: for complex `M`, seems to be unstable, made worse as number of nonzero elements in `M` decreases
"""
function schur_complement(ps::AbstractPardisoSolver,A::SparseMatrixCSC{Tv},n::Integer,T::Symbol=:N) where Tv <: PardisoNumTypes

    n ≥ size(A,1) ? throw(ErrorException("complement block size n=$n≥A.m=$(A.m)")) : nothing
    # Tv<:Complex ? (@warn "unstable for complex types, unknown why") : nothing

    pardisoinit(ps)
    original_phase = get_phase(ps)
    original_iparms = get_iparms(ps)
    set_iparm!(ps,1,1) # use custom IPARM
    set_iparm!(ps,38,n) # set Schur complement block size to n
    set_phase!(ps,12) # analyze and factorize
    B = Matrix{Tv}(undef,size(A,1),0) # dummy array to feed to pardiso

    # transpose via IPARM(12) doesn't work at factorize step (only on entry to solve step)
    if T==:N
        M = permutedims(A)
        set_iparm!(ps, 12, 1)
    elseif T == :C
        M = conj(permutedims(A))
        set_iparm!(ps, 12, 0)
    elseif T == :T
        M = A
        set_iparm!(ps, 12, 0)
    else
        throw(ArgumentError("only :T, :N and :C, are valid transpose symbols"))
    end

    pardiso(ps,B,M,B)
    S = pardisogetschur(ps) # get schur complement matrix

    set_phase!(ps, RELEASE_ALL)
    pardiso(ps, B, M, B)
    set_phase!(ps, original_phase) # reset phase to user setting
    for i ∈ eachindex(original_iparms)
        set_iparm!(ps,i,original_iparms[i])
    end

    return S
end

"""
    pardisogetschur(ps) -> S

retrieve schur complement from PardisoSolver `ps`.
"""
function pardisogetschur(ps::AbstractPardisoSolver)
    nnzschur = get_iparm(ps, 39)
    nschur = get_iparm(ps,38)
    T = isreal(get_matrixtype(ps)) ? Float64 : ComplexF64
    if nnzschur==0
        return spzeros(T,nschur,nschur)
    else
        S = Vector{T}(undef,nnzschur)
        IS = Vector{Int32}(undef,nschur)
        JS = Vector{Int32}(undef,nnzschur)
        ccall_pardiso_get_schur(ps,S,IS,JS)
        IS = pushfirst!(IS,Int32(1)) # some issue with IS (nschur+1 doesn't seem to work)
        S = permutedims(SparseMatrixCSC(nschur,nschur,IS,JS,S)) # really constructing CSR and then transposing
        return S
    end
end

function dim_check(X, A, B)
    size(X) == size(B) || throw(DimensionMismatch(string("solution has $(size(X)), ",
                                                         "RHS has size as $(size(B)).")))
    size(A, 1) == size(B, 1) || throw(DimensionMismatch(string("matrix has $(size(A,1)) ",
                                                               "rows, RHS has $(size(B,1)) rows.")))
    size(B, 1) == stride(B, 2) || throw(DimensionMismatch(
                                            string("Only memory-contiguous RHS supported")))
end

end # module


text:
export Filt

"""
`Filt([domainType=Float64::Type,] dim_in::Tuple, b::AbstractVector, [a::AbstractVector,])`

`Filt(x::AbstractVector, b::AbstractVector, [a::AbstractVector,])`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractVector`, returns a vector `y` filtered by an IIR filter of coefficients `b` and `a`. If only `b` is provided a FIR is used to comute `y` instead.

"""
struct Filt{T,N} <: LinearOperator
	dim_in::NTuple{N,Int}
	b::AbstractVector{T}
	a::AbstractVector{T}
	si::AbstractVector{T}
	function Filt{T,N}(dim_in,b,a) where {T,N}
		isempty(b) && throw(ArgumentError("filter vector b must be non-empty"))
		isempty(a) && throw(ArgumentError("filter vector a must be non-empty"))
		a[1] == 0  && throw(ArgumentError("filter vector a[1] must be nonzero"))
		as = length(a)
		bs = length(b)
		sz = max(as, bs)
		silen = sz - 1

		# Filter coefficient normalization
		if a[1] != 1
			norml = a[1]
			a ./= norml
			b ./= norml
		end
		# Pad the coefficients with zeros if needed
		bs<sz   && (b = copyto!(zeros(eltype(b), sz), b))
		1<as<sz && (a = copyto!(zeros(eltype(a), sz), a))

		si = zeros(promote_type(eltype(b), eltype(a)), max(length(a), length(b))-1)
		new(dim_in,b,a,si)
	end
end

# Constructors

#default constructor
function Filt(T::Type, dim_in::NTuple{N,Int}, b::AbstractVector, a::AbstractVector) where {N}
	eltype(b) != T && error("eltype of b is $(eltype(b)), should be $T")
	eltype(a) != T && error("eltype of a is $(eltype(b)), should be $T")
	Filt{T,N}(dim_in,b,a)
end

function Filt(T::Type, dim_in::NTuple{N,Int}, b::AbstractVector) where {N}
	eltype(b) != T && error("eltype of b is $(eltype(b)), should be $T")
	Filt{T,N}(dim_in,b,[convert(T,1.0)])
end

Filt(dim_in::Int,  b::AbstractVector, a::AbstractVector) =
Filt(eltype(b),(dim_in,), b, a)

Filt(dim_in::Tuple,  b::AbstractVector, a::AbstractVector) =
Filt(eltype(b), dim_in, b, a)

Filt(dim_in::Int,  b::AbstractVector) =
Filt(eltype(b),(dim_in,), b)

Filt(dim_in::Tuple,  b::AbstractVector) =
Filt(eltype(b), dim_in, b)

Filt(x::AbstractArray, b::AbstractVector, a::AbstractVector) =
Filt(eltype(x),size(x), b, a)

Filt(x::AbstractArray, b::AbstractVector) =
Filt(size(x), b)

# Mappings

function mul!(y::AbstractArray{T},L::Filt{T,N},x::AbstractArray{T}) where {T,N}
	for col = 1:size(x,2)
		length(L.a) != 1 ? iir!(y,L.b,L.a,x,L.si,col,col) : fir!(y,L.b,x,L.si,col,col)
	end
end

function mul!(y::AbstractArray{T},L::AdjointOperator{Filt{T,N}},x::AbstractArray{T}) where {T,N}
    A = L.A
	for col = 1:size(x,2)
		length(A.a) != 1 ? iir_rev!(y,A.b,A.a,x,A.si,col,col) : fir_rev!(y,A.b,x,A.si,col,col)
	end
end

function iir!(y, b, a, x, si, coly, colx)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,colx]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi - a[j+1]*val
        end
        si[silen] = b[silen+1]*xi - a[silen+1]*val
        y[i,coly] = val
    end
    si .= 0. #reset state
end

# Utilities

function iir_rev!(y, b, a, x, si, coly, colx)
    silen = length(si)
    @inbounds for i=size(x, 1):-1:1
        xi = x[i,colx]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi - a[j+1]*val
        end
        si[silen] = b[silen+1]*xi - a[silen+1]*val
        y[i,coly] = val
    end
    si .= 0.
end

function fir!(y, b, x, si, coly, colx)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,colx]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi
        end
        si[silen] = b[silen+1]*xi
        y[i,coly] = val
    end
    si .= 0.
end

function fir_rev!(y, b, x, si, coly, colx)
    silen = length(si)
    @inbounds for i=size(x, 1):-1:1
        xi = x[i,colx]
        val = si[1] + b[1]*xi
        for j=1:(silen-1)
            si[j] = si[j+1] + b[j+1]*xi
        end
        si[silen] = b[silen+1]*xi
        y[i,coly] = val
    end
    si .= 0.
end

# Properties

domainType(L::Filt{T, N}) where {T, N} = T
codomainType(L::Filt{T, N}) where {T, N} = T

size(L::Filt) = L.dim_in, L.dim_in

fun_name(L::Filt)  = size(L.a,1) != 1 ? "IIR" : "FIR"

#TODO find out a way to verify this,
# probably for IIR it means zeros inside unit circle
is_full_row_rank(L::Filt)    = true
is_full_column_rank(L::Filt) = true


text:
#function DemoLinRegCUDA4()

    # note that this linear regression code is not the fastest way to do this since it is wasteful in terms of storage. We would need to code an operation Ax+beta*b where A,x,b are CudaArrays and beta is a Float64. This would mean we can then avoid the separate A*x (which stores this on the graph) and then Ax-y. Hence we can cut the storage significantly this way.

    #include("setup.jl")

    #AtransB(A,B)=A'*B

    StartCode()

    A=ADvariable()
#B=ADvariable()
#    loss=sum(AtransB(A,A))
#    loss=meanSquare(A')
#     loss=sum(A'-B')
    # loss=sumSquare(diagm(B)*A*diagm(B))
#    loss=sumSquare(diagm(A.*A))
    loss=sumSquare(abs(A))

    net=EndCode() # defines the graph

    # instantiate parameter nodes and inputs:
    #value=Array(Any,NodeCounter()) # function values on the nodes

    net.value[A]=randn(2,2)
    #net.value[B]=randn(2,1)

    net=compile(net;debug=true) # compile and preallocate memory

    @gpu CUDArt.init([0])
@gpu net=convert(net,"GPU")

gradcheck(net;showgrad=true) # use a small number of datapoints and small network to check the gradient, otherwise this will be very slow
 #   @cpu gradcheck(net;showgrad=true) # use a small number of datapoints and small network to check the gradient, otherwise this will be very slow
#@gpu gradcheckGPU(net;showgrad=false)

#end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
The anchor for the logging system. It maintains a list of all the loggers
requested by applications and a reference to the configuration.
It will be atomically updated whenever a reconfigure occurs.
"""
type LoggerContext <: LifeCycle.Object
    name::AbstractString
    loggers::LOGGERS
    configLocation::AbstractString
    config::Configuration
    state::LifeCycle.State
end
LoggerContext(name::AbstractString, cfgloc::AbstractString, cfg::Configuration) =
    LoggerContext(name, LOGGERS(), cfgloc, cfg, LifeCycle.INITIALIZED)
LoggerContext(name::AbstractString, cfgloc::AbstractString="") =
    LoggerContext(name, cfgloc, DefaultConfiguration())

show(io::IO, ctx::LoggerContext) = print(io, "LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]")

"Setup context with a new configuration and return previous one"
function setconfig!(ctx::LoggerContext, cfg::Configuration)
    prevcfg = ctx.config

    if isdefined(cfg, :properties)
        !haskey(cfg.properties, "host") && setindex!(cfg.properties, "host", gethostname())
        !haskey(cfg.properties, "context") && setindex!(cfg.properties, "context", ctx.name)
    end
    start(cfg)
    ctx.config = cfg
    ctx.configLocation = source(cfg)
    prevcfg !== nothing && stop(prevcfg)

    return prevcfg
end

"Setup context with a configuration located at the provided location"
function setconfig!(ctx::LoggerContext, cfgloc::AbstractString)
    ctx.configLocation = cfgloc
    reconfigure!(ctx)
end

"Reconfigure context"
function reconfigure!(ctx::LoggerContext)
    debug(LOGGER, "Reconfiguration started for context[name=$(ctx.name)] at $(ctx.configLocation)")
    cfg = getconfig(ctx.configLocation, ctx.name, current_module())
    @assert cfg !== nothing "No configuration found"
    setconfig!(ctx, cfg)
    debug(LOGGER, "Reconfiguration complete for context[name=$(ctx.name)] at $(ctx.configLocation)")
end

"Start context."
function start(ctx::LoggerContext)
    debug(LOGGER, "Starting LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]...")
    if state(ctx) == LifeCycle.INITIALIZED || state(ctx) == LifeCycle.STOPPED
        state!(ctx, LifeCycle.STARTING)
        reconfigure!(ctx)
        # add shutdown hook
        atexit(()->begin
                    debug(LOGGER, symbol("SHUTDOWN HOOK"), "Stopping LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]")
                    stop(ctx)
                   end)
        state!(ctx, LifeCycle.STARTED)
    end
    debug(LOGGER, "LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))] started OK.")
end

"Start context with a specific configuration."
function start(ctx::LoggerContext, cfg::Configuration)
    debug(LOGGER, "Starting LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]...")
    if state(ctx) == LifeCycle.INITIALIZED || state(ctx) == LifeCycle.STOPPED
        state!(ctx, LifeCycle.STARTING)
        # add shutdown hook
        atexit(()->begin
                    debug(LOGGER, symbol("SHUTDOWN HOOK"), "Stopping LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]")
                    stop(ctx)
                   end)
        state!(ctx, LifeCycle.STARTED)
    end
    setconfig!(ctx, cfg)
    debug(LOGGER, "LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))] started OK.")
end

"Stop context."
function stop(ctx::LoggerContext)
    debug(LOGGER, "Stopping LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))]...")
    state(ctx) == LifeCycle.STOPPED && return
    state!(ctx, LifeCycle.STOPPING)

    prevcfg = ctx.config
    ctx.config = NullConfiguration()
    stop(prevcfg)

    state!(ctx, LifeCycle.STOPPED)
    debug(LOGGER, "Stopped LoggerContext[name=$(ctx.name), state=$(string(state(ctx)))].")
end

"Checks if a logger with the specified name exists."
in(lname::AbstractString, ctx::LoggerContext) = haskey(ctx.loggers, lname)

"Returns a logger from a logger context"
function logger(ctx::LoggerContext, lname::AbstractString,
                msgfactory::DataType=LOG4JL_DEFAULT_MESSAGE)
    # return logger in exists
    lname in ctx && return ctx.loggers[lname]

    # otherwise create new logger and return it
    lgr = Logger(lname, msgfactory, logger(ctx.config, lname), filter(ctx.config))
    ctx.loggers[lname] = lgr
    return lgr
end

"Returns all loggers."
loggers(ctx::LoggerContext) = values(ctx.loggers)


text:
using LinearAlgebra

export makeLinearCubicPP, makeCubicPP, C2, C2Hyman89, C2HymanNonNegative, C2MP, Bessel, HuynRational, VanAlbada, VanLeer, FritschButland, Brodlie
export evaluateDerivative, evaluateSecondDerivative

abstract type DerivativeKind end
struct C2 <: DerivativeKind end
struct C2Hyman89 <: DerivativeKind end
struct C2HymanNonNegative <: DerivativeKind end
struct C2MP <: DerivativeKind end
struct C2MP2 <: DerivativeKind end
@enum PPBoundary NOT_A_KNOT = 0 FIRST_DERIVATIVE = 1 SECOND_DERIVATIVE = 2 FIRST_DIFFERENCE = 3

abstract type LimiterDerivative <: DerivativeKind end
struct Hermite <: DerivativeKind end
struct Bessel <: LimiterDerivative end
struct HuynRational <: LimiterDerivative end
struct VanLeer <: LimiterDerivative end
struct VanAlbada <: LimiterDerivative end
struct FritschButland <: LimiterDerivative end #Fritsch Butland 1980
struct Brodlie <: LimiterDerivative end #Fritch Butland 1984

struct PP{N,T<:Real,TX}
    a::AbstractArray{T}
    b::AbstractArray{T}
    c::AbstractMatrix{T} #c[i,:] = coeff of x^{i+1}
    x::AbstractArray{TX}
    PP(N::Int, T, TX, n::Int) = new{N,T,TX}(zeros(T, n), zeros(T, n), zeros(T, (N - 1, n - 1)), zeros(TX, n))
    PP(N::Int, a::AbstractArray{T}, b::AbstractArray{T}, c::AbstractMatrix{T}, x::AbstractArray{TX}) where {T<:Real,TX} =
        new{N,T,TX}(a, b, c, x)
end

Base.length(p::PP) = Base.length(p.x)
Base.size(p::PP) = Base.size(p.x)
Base.broadcastable(p::PP) = Ref(p)

function makeLinearPP(x::AbstractArray{TX}, y::AbstractArray{T}) where {T,TX}
    pp = PP(2, T, TX, length(y))
    computeLinearPP(pp, x, y)
    return pp
end

function computeLinearPP(pp::PP{N,T,TX}, x::AbstractArray{TX}, y::AbstractArray{T}) where {N,T,TX}
    n = length(x)
    if n <= 1
        pp.a[1:end] = y
        pp.b[1:end] = zeros(n)
        pp.c[1:end, 1:end] = zeros((N - 1), n - 1)
        pp.x[1:end] = x
    elseif n == 2
        t = y[2] - y[1]
        b = zeros(n)
        if abs(x[2] - x[1]) > eps()
            b[1] = t / (x[2] - x[1])
        end
        b[2] = b[1]
        pp.a[1:end] = y
        pp.b[1:end] = b
        pp.c[1:end, 1:end] = zeros((N - 1), n - 1)
        pp.x[1:end] = x
    else
        # on xi, xi+1, f(x)= yi (xi+1-x) + yi (x-xi) = A + B (x-xi) => B = (yi-yi+1)/(xi-xi+1)
        b = zeros(n)
        for i = 1:n-1
            b[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
        end
        pp.a[1:end] = y
        pp.b[1:end] = b
        pp.c[1:end, 1:end] = zeros((N - 1), n - 1)
        pp.x[1:end] = x
    end
end

function makeCubicPP(
    x::AbstractArray{TX},
    y::AbstractArray{T},
    leftBoundary::PPBoundary,
    leftValue::T,
    rightBoundary::PPBoundary,
    rightValue::T,
    kind::DerivativeKind,
) where {T,TX}
    pp = PP(3, T, TX, length(y))
    computePP(pp, x, y, leftBoundary, leftValue, rightBoundary, rightValue, kind)
    return pp
end

function computePP(
    pp::PP{3,T,TX},
    x::AbstractArray{TX},
    y::AbstractArray{T},
    leftBoundary::PPBoundary,
    leftValue::T,
    rightBoundary::PPBoundary,
    rightValue::T,
    kind::Union{C2,C2Hyman89,C2HymanNonNegative,C2MP,C2MP2},
) where {T,TX}
    n = length(y)
    if n <= 2
        computeLinearPP(pp, x, y)
        return
    end

    dx = x[2:end] - x[1:end-1]
    S = (y[2:end] - y[1:end-1]) ./ dx
    middle = zeros(TX, n)
    alpha = zeros(T, n)
    lower = zeros(TX, n - 1)
    upper = zeros(TX, n - 1)
    for i = 2:n-1
        lower[i-1] = dx[i]
        upper[i] = dx[i-1]
        middle[i] = 2 * (dx[i] + dx[i-1])
        alpha[i] = 3 * (dx[i] * S[i-1] + dx[i-1] * S[i])
    end
    #middle[2:n-1] = 3 * (dx[2:n-1] + dx[1:n-2])
    #alpha[2:n-1] = 3 * (dx[2:n-1] .* S[1:n-2] + dx[1:n-2] .* S[2:n-1])
    if leftBoundary == NOT_A_KNOT
        middle[1] = dx[2] * (dx[2] + dx[1])
        upper[1] = (dx[2] + dx[1]) * (dx[2] + dx[1])
        alpha[1] = S[1] * dx[2] * (2 * dx[2] + 3 * dx[1]) + S[2] * dx[1]^2
    elseif leftBoundary == FIRST_DERIVATIVE
        middle[1] = one(TX)
        upper[1] = zero(TX)
        alpha[1] = leftValue
    elseif leftBoundary == FIRST_DIFFERENCE
        middle[1] = one(TX)
        upper[1] = zero(TX)
        alpha[1] = S[1]
    elseif leftBoundary == SECOND_DERIVATIVE
        middle[1] = 2 * one(TX)
        upper[1] = one(TX)
        alpha[1] = 3 * S[1] - leftValue * dx[1] / 2
    end
    if rightBoundary == NOT_A_KNOT
        lower[n-1] = -(dx[n-1] + dx[n-2]) * (dx[n-1] + dx[n-2])
        middle[n] = -dx[n-2] * (dx[n-2] + dx[n-1])
        alpha[n] = -S[n-2] * dx[n-1]^2 - S[n-1] * dx[n-2] * (3 * dx[n-1] + 2 * dx[n-2])
    elseif rightBoundary == FIRST_DERIVATIVE
        middle[n] = one(TX)
        lower[n-1] = zero(TX)
        alpha[n] = rightValue
    elseif rightBoundary == FIRST_DIFFERENCE
        middle[n] = one(TX)
        lower[n-1] = zero(TX)
        alpha[n] = S[n-1]
    elseif rightBoundary == SECOND_DERIVATIVE
        middle[n] = 2 * one(TX)
        lower[n-1] = zero(TX)
        alpha[n] = 3 * S[n-1] - rightValue * dx[n-1] / 2
    end
    tri = LinearAlgebra.Tridiagonal(lower, middle, upper)
    fPrime = tri \ alpha
    filterSlope(kind, y, fPrime, dx, S)
    c = pp.c
    for i = 1:n-1
        c[1, i] = (3 * S[i] - fPrime[i+1] - 2 * fPrime[i]) / dx[i]
        c[2, i] = (fPrime[i+1] + fPrime[i] - 2 * S[i]) / (dx[i]^2)
    end
    pp.a[1:end] = y
    pp.b[1:end] = fPrime
    pp.x[1:end] = x
end

function computePP(
    pp::PP{3,T,TX},
    x::AbstractArray{TX},
    y::AbstractArray{T},
    leftBoundary::PPBoundary,
    leftValue::T,
    rightBoundary::PPBoundary,
    rightValue::T,
    limiter::LimiterDerivative,
) where {T,TX}
    n = length(y)
    if n <= 2
        computeLinearPP(pp, x, y)
        return
    end

    dx = x[2:end] - x[1:end-1]
    S = (y[2:end] - y[1:end-1]) ./ dx
    b = pp.b
    fillDerivativeEstimate(limiter, dx, S, b)
    if leftBoundary == NOT_A_KNOT
        b[1] = S[2] * dx[1] / (dx[2] * (dx[2] + dx[1])) - S[1] * ((dx[2] / dx[1] + 2) / (dx[2] + dx[1]))
    elseif leftBoundary == FIRST_DIFFERENCE
        b[1] = S[1]
    elseif leftBoundary == FIRST_DERIVATIVE
        b[1] = leftValue
    elseif leftBoundary == SECOND_DERIVATIVE
        #c[1] = leftValue * 0.5
        b[1] = (-leftValue / 2 * dx[1] - b[2] + 3 * S[1]) / 2
    end
    if rightBoundary == NOT_A_KNOT
        b[n] =
            S[n-2] * dx[n-1] / (dx[n-2] * (dx[n-2] + dx[n-1])) -
            S[n-1] * ((dx[n-2] / dx[n-1] + 2) / (dx[n-2] + dx[n-1]))
    elseif rightBoundary == FIRST_DERIVATIVE
        b[n] = rightValue
    elseif rightBoundary == FIRST_DIFFERENCE
        b[n] = S[n-1]
    elseif rightBoundary == SECOND_DERIVATIVE
        b[n] = (rightValue * dx[n-1] + 6 * S[n-1] - 2 * b[n-1]) / 4
    end
    c = pp.c
    for i = 1:n-1
        c[1, i] = (3 * S[i] - b[i+1] - 2 * b[i]) / dx[i]
        c[2, i] = (b[i+1] + b[i] - 2 * S[i]) / (dx[i]^2)
    end
    pp.a[1:end] = y
    pp.x[1:end] = x
end

function evaluate(self::PP{3,T,TX}, z::TZ) where {T,TX,TZ}
    if z <= self.x[1]
        return self.b[1] * (z - self.x[1]) + self.a[1]
    elseif z >= self.x[end]
        return self.b[end] * (z - self.x[end]) + self.a[end]
    end
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if z != self.x[i] && i > 1
        i -= 1
    end
    h = z - self.x[i]
    return self.a[i] + h * (self.b[i] + h * (self.c[1, i] + h * (self.c[2, i])))
end

function evaluateDerivative(self::PP{3,T,TX}, z::TZ) where {T,TX,TZ}
    if z <= self.x[1]
        return self.b[1]
    elseif z >= self.x[end]
        rightSlope = self.b[end]
        return rightSlope
    end
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if z != self.x[i] && i > 1
        i -= 1
    end
    h = z - self.x[i]
    return self.b[i] + h * (2 * self.c[1, i] + h * (3 * self.c[2, i]))
end
function evaluateSecondDerivative(self::PP{3,T,TX}, z::TZ) where {T,TX,TZ}
    if z <= self.x[1]
        return self.b[1]
    elseif z >= self.x[end]
        rightSlope = self.b[end]
        return rightSlope
    end
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if z != self.x[i] && i > 1
        i -= 1
    end
    h = z - self.x[i]
    return 2 * self.c[1, i] + h * (3 * 2 * self.c[2, i])
end

(spl::PP{N,T,TX})(x::TZ) where {N,T,TX,TZ} = evaluate(spl, x)
function (spl::PP{N,T,TX})(x::AbstractArray) where {N,T,TX}
    evaluate.(spl, x)
end

function filterSlope(kind::C2, y::AbstractArray{T}, b::AbstractArray{T}, dx::AbstractArray{TX}, S::AbstractArray{T}) where {T,TX}
    #do nothing
end

function filterSlope(kind::C2MP, y::AbstractArray{T}, b::AbstractArray{T}, dx::AbstractArray{TX}, S::AbstractArray{T}) where {T,TX}
    n = length(y)
    if S[1] > 0
        b[1] = min(max(0, b[1]), 3 * S[1])
    else
        b[1] = max(min(0, b[1]), 3 * S[1])
    end
    if S[n-1] > 0
        b[n] = min(max(0, b[n]), 3 * S[n-1])
    else
        b[n] = max(min(0, b[n]), 3 * S[n-1])
    end

    for i = 2:n-1
        Sim = S[i-1]
        Sip = S[i]
        if Sim * Sip <= 0
            b[i] = 0
        elseif Sim > 0 && Sip > 0
            b[i] = min(max(0, b[i]), 3 * min(Sim, Sip))
        else
            b[i] = max(min(0, b[i]), 3 * max(Sim, Sip))
        end
    end
end

minmod(s::T, t::T) where {T} = s * t <= 0 ? zero(T) : sign(s) * min(abs(s), abs(t))

function filterSlope(kind::C2MP2, y::AbstractArray{T}, b::AbstractArray{T}, dx::AbstractArray{TX}, S::AbstractArray{T}) where {T,TX}
    #r = s/t and G = t*g(r)
    n = length(y)
    b[1] = minmod(b[1], 3 * S[1])
    b[n] = minmod(b[n], 3 * S[n-1])

    for i = 2:n-1
        Sim = S[i-1]
        Sip = S[i]
        lowerBound = minmod(Sim, Sip)
        upperBound = (sign(Sim) + sign(Sip)) / 2 * min(max(abs(Sim), abs(Sip)), 3 * min(abs(Sim), abs(Sip)))
        b[i] = min(max(lowerBound, b[i]), upperBound)
    end
end

function filterSlope(kind::C2HymanNonNegative, y::AbstractArray{T}, b::AbstractArray{T}, dx::AbstractArray{TX}, S::AbstractArray{T}) where {T,TX}
    n = length(y)
    tau0 = sign(y[1])
    #Warning the paper (3.3) is wrong as it does not obey (3.1)
    b[1] = tau0 * max(-3 * tau0 * y[1] / dx[1], tau0 * b[1])
    for i = 2:n-1
        taui = sign(y[i])
        b[i] = taui * min(3 * taui * y[i] / dx[i-1], max(-3 * taui * y[i] / dx[i], taui * b[i]))
    end
    taun = sign(y[n])
    b[n] = taun * min(3 * taun * y[n] / dx[n-1], taun * b[n])
end

function filterSlope(kind::C2Hyman89, y::AbstractArray{T}, b::AbstractArray{T}, dx::AbstractArray{TX}, S::AbstractArray{T}) where {T,TX}
    n = length(y)
    tmp = b

    local correction, pm, pu, pd, M
    if (tmp[1] * S[1]) > 0
        correction = sign(tmp[1]) * min(abs(tmp[1]), abs(3 * S[1]))
    else
        correction = zero(T)
    end
    if correction != tmp[1]
        tmp[1] = correction
    end
    for i = 2:n-1
        pm = ((S[i-1] * dx[i]) + (S[i] * dx[i-1])) / (dx[i-1] + dx[i])
        M = 3 * min(min(abs(S[i]), abs(S[i-1])), abs(pm))
        if i > 2
            if ((S[i-1] - S[i-2]) * (S[i] - S[i-1])) > 0
                pd = ((S[i-1] * ((2 * dx[i-1]) + dx[i-2])) - (S[i-2] * dx[i-1])) / (dx[i-2] + dx[i-1])
                if (pm * pd) > 0 && (pm * (S[i-1] - S[i-2])) > 0
                    M = max(M, 3 * min(abs(pm), abs(pd)) / 2)
                end
            end
        end
        if i < (n - 1)
            if ((S[i] - S[i-1]) * (S[i+1] - S[i])) > 0
                pu = ((S[i] * ((2 * dx[i]) + dx[i+1])) - (S[i+1] * dx[i])) / (dx[i] + dx[i+1])
                if ((pm * pu) > 0) && ((-pm * (S[i] - S[i-1])) > 0)
                    M = max(M, 3 * min(abs(pm), abs(pu)) / 2)
                end
            end
        end
        if (tmp[i] * pm) > 0
            correction = sign(tmp[i]) * min(abs(tmp[i]), M)
        else
            correction = zero(T)
        end
        if correction != tmp[i]
            tmp[i] = correction
        end
    end
    if (tmp[n] * S[n-1]) > 0
        correction = sign(tmp[n]) * min(abs(tmp[n]), abs(3 * S[n-1]))
    else
        correction = zero(T)
    end
    if correction != tmp[n]
        tmp[n] = correction
    end
end

function estimateDerivativeParabolic(x::AbstractArray{TX}, y::AbstractArray{T}) where {T,TX}
    n = length(x)
    b = zeros(n)
    dx = x[2:end] - x[1:end-1]
    S = (y[2:end] - y[1:end-1]) ./ dx
    for i = 2:n-1
        b[i] = (dx[i-1] * S[i] + dx[i] * S[i-1]) / (dx[i-1] + dx[i])
    end
    return b
end

function limit(::HuynRational, s::T, t::T) where {T}
    st = s * t
    if st <= 0
        return zero(T)
    end
    return st * 3 * (s + t) / (s^2 + 4 * st + t^2)
end

function limit(::VanAlbada, s::T, t::T) where {T}
    st = s * t
    if st == 0
        return s
    end
    return st * (s + t) / (s^2 + t^2)
end

function limit(::VanLeer, s::T, t::T) where {T}
    st = s * t
    if st <= 0
        return zero(T)
    end
    return 2 * st / (s + t) #warning, the product s*t can be zero even when s or t are not 0, this rewrite helps saving accuracy
end

function limit(::FritschButland, s::T, t::T) where {T}
    st = s * t
    if st <= 0
        return zero(T)
    end
    if abs(s) <= abs(t)
        return 3 * st / (2 * s + t)
    end
    return 3 * st / (s + 2 * t)   # (1+dxp / (dxp+dx)) * t + (2 - dxp / (dxp+dx)) * s
end

function fillDerivativeEstimate(limiter::LimiterDerivative, dx::AbstractArray{TX}, S::AbstractArray{T}, b::AbstractArray{T}) where {T,TX}
    n = length(S)
    for i = 2:n
        b[i] = limit(limiter, S[i-1], S[i])
    end
end

function fillDerivativeEstimate(limiter::Bessel, dx::AbstractArray{TX}, S::AbstractArray{T}, b::AbstractArray{T}) where {T,TX}
    # @. b[2:end-1] = (dx[1:end-1] * S[2:end] + dx[2:end] * S[1:end-1]) / (dx[1:end-1] + dx[2:end])
    n = length(S)
    for i = 2:n
        b[i] = (dx[i-1] * S[i] + dx[i] * S[i-1]) / (dx[i-1] + dx[i])
    end
end

function fillDerivativeEstimate(limiter::Brodlie, dx::AbstractArray{TX}, S::AbstractArray{T}, b::AbstractArray{T}) where {T,TX}
    n = length(S)
    for i = 2:n
        s, t = S[i-1], S[i]
        st = s * t
        if st == 0
            b[i] = s
        else
            α = (dx[i-1] + 2 * dx[i]) / (3 * (dx[i-1] + dx[i]))
            b[i] = (st) / (α * t + (1 - α) * s)
        end
    end
end


text:
module Door

using ..Ahorn, Maple

const textures = ["wood", "metal"]
const placements = Ahorn.PlacementDict(
    "Door ($(uppercasefirst(texture)))" => Ahorn.EntityPlacement(
        Maple.Door,
        "point",
        Dict{String, Any}(
            "type" => texture
        )
    ) for texture in textures
)

function doorSprite(entity::Maple.Door)
    variant = get(entity.data, "type", "wood")

    return variant == "wood" ? "objects/door/door00.png" : "objects/door/metaldoor00.png"
end

Ahorn.editingOptions(entity::Maple.Door) = Dict{String, Any}(
    "type" => textures
)

function Ahorn.selection(entity::Maple.Door)
    x, y = Ahorn.position(entity)
    sprite = doorSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Door, room::Maple.Room)
    sprite = doorSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
using Distributions
using ProgressMeter
using Base.LinAlg.BLAS
using Compat
using Devectorize
using HDF5
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np

import StatsBase.fit

function update_weights!(rbm, h_pos, v_pos, h_neg, v_neg, lr; approx="CD")
    # dW = zeros(size(rbm.W))
    # dW = lr * ( (h_pos * v_pos') - (h_neg * v_neg') )
    gemm!('N', 'T', lr, h_neg, v_neg, 0.0, rbm.dW)
    gemm!('N', 'T', lr, h_pos, v_pos, -1.0, rbm.dW)

    if contains(approx,"tap")
        buf2 = gemm('N', 'T', h_neg-abs2(h_neg), v_neg-abs2(v_neg)) .* rbm.W
        axpy!(-lr, buf2, rbm.dW)
    end

    if approx == "tap3"
        buf3 = gemm('N','T', (h_neg-abs2(h_neg)) .* (0.5-h_neg), (v_neg-abs2(v_neg)) .* (0.5-v_neg)) .* rbm.W2
        axpy!(-2.0*lr, buf3, rbm.dW)
    end

    # rbm.dW += rbm.momentum * rbm.dW_prev
    axpy!(rbm.momentum, rbm.dW_prev, rbm.dW)
    # rbm.W +=  dW
    axpy!(1.0, rbm.dW, rbm.W)
    rbm.W2=rbm.W.*rbm.W

    if approx == "tap3"
        rbm.W3=rbm.W2.*rbm.W
    end
    # save current dW
    copy!(rbm.dW_prev, rbm.dW)
end

function update_weights_QuadraticPenalty!(rbm, h_pos, v_pos, h_neg, v_neg, lr, decay_mag; approx="CD")
    # dW = zeros(size(rbm.W))

    gemm!('N', 'T', lr, h_neg, v_neg, 0.0, rbm.dW)
    gemm!('N', 'T', lr, h_pos, v_pos, -1.0, rbm.dW)

    if contains(approx,"tap")
        buf2 = gemm('N', 'T', h_neg-abs2(h_neg), v_neg-abs2(v_neg)) .* rbm.W
        axpy!(-lr, buf2, rbm.dW)
    end

    if approx == "tap3"
        buf3 = gemm('N','T', (h_neg-abs2(h_neg)) .* (0.5-h_neg), (v_neg-abs2(v_neg)) .* (0.5-v_neg)) .* rbm.W2
        axpy!(-2.0*lr, buf3, rbm.dW)
    end

    # rbm.W += rbm.momentum * rbm.dW_prev
    axpy!(rbm.momentum, rbm.dW_prev, rbm.dW)

    # Apply Weight-Decay Penalty
    # rbm.W += -lr * L2-Penalty-Gradient
    axpy!(-lr*decay_mag,rbm.W,rbm.dW)

    # rbm.W +=  dW
    axpy!(1.0, rbm.dW, rbm.W)

    rbm.W2=rbm.W.*rbm.W

    if approx == "tap3"
        rbm.W3=rbm.W2.*rbm.W
    end
    # save current dW
    copy!(rbm.dW_prev, rbm.dW)
end

function update_weights_LinearPenalty!(rbm, h_pos, v_pos, h_neg, v_neg, lr, decay_mag ; approx="CD")
    # dW = zeros(size(rbm.W))
    # dW = (h_pos * v_pos') - (h_neg * v_neg')
    gemm!('N', 'T', lr, h_neg, v_neg, 0.0, rbm.dW)          # Not flushing rbm.dW since we multiply w/ 0.0
    gemm!('N', 'T', lr, h_pos, v_pos, -1.0,rbm.dW)

    if contains(approx,"tap")
        buf2 = gemm('N', 'T', h_neg-abs2(h_neg), v_neg-abs2(v_neg)) .* rbm.W
        axpy!(-lr, buf2, rbm.dW)
    end

    if approx == "tap3"
        buf3 = gemm('N','T', (h_neg-abs2(h_neg)) .* (0.5-h_neg), (v_neg-abs2(v_neg)) .* (0.5-v_neg)) .* rbm.W2
        axpy!(-2.0*lr, buf3, rbm.dW)
    end

    # rbm.W += rbm.momentum * rbm.dW_prev
    axpy!(rbm.momentum, rbm.dW_prev, rbm.dW)

    # Apply Weight-Decay Penalty
    # rbm.W += -lr * L1-Penalty-Gradient
    axpy!(-lr*decay_mag,sign(rbm.W),rbm.dW)

    # rbm.W += lr * dW
    axpy!(1.0, rbm.dW, rbm.W)

    rbm.W2=rbm.W.*rbm.W

    if approx == "tap3"
        rbm.W3=rbm.W2.*rbm.W
    end
    # save current dW
    copy!(rbm.dW_prev, rbm.dW)
end

function contdiv(rbm::RBM, vis::Mat{Float64}, n_gibbs::Int; approx="CD")
    if approx == "CD"
        v_pos, h_pos, v_neg, h_neg = gibbs(rbm, vis; n_times=n_gibbs)
    else
        v_pos, h_pos, v_neg, h_neg = iter_mag(rbm, vis; n_times=n_gibbs, approx=approx)
    end
    return v_pos, h_pos, v_neg, h_neg
end

function persistent_contdiv(rbm::RBM, vis::Mat{Float64}, n_gibbs::Int; approx="CD")
    if size(rbm.persistent_chain_vis) != size(vis)
        # persistent_chain not initialized or batch size changed, re-initialize
        rbm.persistent_chain_vis = vis
        rbm.persistent_chain_hid = hid_means(rbm, vis)
    end

    if approx == "CD"
        # take positive samples from real data
        v_pos, h_pos, _, _ = gibbs(rbm, vis; n_times=1)
        # take negative samples from "fantasy particles"
        _, _, v_neg, h_neg = gibbs(rbm, rbm.persistent_chain_vis; n_times=n_gibbs)
        copy!(rbm.persistent_chain_vis,v_neg)
    else
        v_pos, h_pos, v_neg, h_neg = iter_mag_persist!(rbm, vis; n_times=n_gibbs, approx=approx)
    end

    return v_pos, h_pos, v_neg, h_neg
end

function fit_batch!(rbm::RBM, vis::Mat{Float64};
                    persistent=true, lr=0.1, n_gibbs=1,
                    weight_decay="none",decay_magnitude=0.01, approx="CD")

    sampler = persistent ? persistent_contdiv : contdiv
    v_pos, h_pos, v_neg, h_neg = sampler(rbm, vis, n_gibbs; approx=approx)

    # Gradient Update on Weights
    if weight_decay=="l2"
        update_weights_QuadraticPenalty!(rbm, h_pos, v_pos, h_neg, v_neg, lr, decay_magnitude, approx=approx)
    elseif weight_decay=="l1"
        update_weights_LinearPenalty!(rbm, h_pos, v_pos, h_neg, v_neg, lr, decay_magnitude, approx=approx)
    else
        update_weights!(rbm, h_pos, v_pos, h_neg, v_neg, lr, approx=approx)
    end

    rbm.hbias += vec(lr * (sum(h_pos, 2) - sum(h_neg, 2)))
    rbm.vbias += vec(lr * (sum(v_pos, 2) - sum(v_neg, 2)))

    return rbm
end

function fit(rbm::RBM, X::Mat{Float64};
             persistent=true, lr=0.1, n_iter=10, batch_size=100, n_gibbs=1,
             weight_decay="none",decay_magnitude=0.01,validation=[],
             monitor_every=5,monitor_vis=false, approx="CD",
             persistent_start=1)
#=
The core RBM training function. Learns the weights and biasings using
either standard Contrastive Divergence (CD) or Persistent CD, depending on
the user options.

### Required Inputs
- *rbm:* RBM object, initialized by `RBM()`/`GRBM()`
- *X:* Set of training vectors

### Optional Inputs
 - *persistent:* Whether or not to use persistent-CD [default=true]
 - *persistent_start:* At which epoch to start using the persistent chains. Only
                       applicable for the case that `persistent=true`.
                       [default=1]
 - *lr:* Learning rate [default=0.1]
 - *n_iter:* Number of training epochs [default=10]
 - *batch_size:* Minibatch size [default=100]
 - *n_gibbs:* Number of Gibbs sampling steps on the Markov Chain [default=1]
 - *weight_decay:* A string value representing the regularization to add to apply to the
                   weight magnitude during training {"none","l1","l2"}. [default="none"]
 - *decay_magnitude:* Relative importance assigned to the weight regularization. Smaller
                      values represent less regularization. Should be in range (0,1).
                      [default=0.01]
 - *validation:* An array of validation samples, e.g. a held out set of training data.
                 If passed, `fit` will also track generalization progress during training.
                 [default=empty-set]
 - *score_every:* Controls at which epoch the progress of the fit is monitored. Useful to
                  speed up the fit procedure if detailed progress monitoring is not required.
                  [default=5]
=#
    @assert minimum(X) >= 0 && maximum(X) <= 1

    n_valid=0
    n_features = size(X, 1)
    n_samples = size(X, 2)
    n_hidden = size(rbm.W,1)
    n_batches = @compat Int(ceil(n_samples / batch_size))
    N = n_hidden+n_features

    # Check for the existence of a validation set
    flag_use_validation=false
    if length(validation)!=0
        flag_use_validation=true
        n_valid=size(validation,2)
    end

    # Create the historical monitor
    ProgressMonitor = Monitor(n_iter,monitor_every;monitor_vis=monitor_vis,
                                                   validation=flag_use_validation)

    # Print info to user
    m_ = rbm.momentum
    info("=====================================")
    info("RBM Training")
    info("=====================================")
    info("  + Training Samples:   $n_samples")
    info("  + Features:           $n_features")
    info("  + Hidden Units:       $n_hidden")
    info("  + Epochs to run:      $n_iter")
    info("  + Persistent ?:       $persistent")
    info("  + Training approx:    $approx")
    info("  + Momentum:           $m_")
    info("  + Learning rate:      $lr")
    info("  + Gibbs Steps:        $n_gibbs")
    info("  + Weight Decay?:      $weight_decay")
    info("  + Weight Decay Mag.:  $decay_magnitude")
    info("  + Validation Set?:    $flag_use_validation")
    info("  + Validation Samples: $n_valid")
    info("=====================================")

    # Scale the learning rate by the batch size
    lr=lr/batch_size

    # Random initialization of the persistent chain
    # It is okay if it isn't used in the actual training procedure.
    p = shuffle!(collect(1:n_samples))[1:batch_size]
    rbm.persistent_chain_vis = Array(Float64,n_features,batch_size)
    for i=1:batch_size
        rbm.persistent_chain_vis[:,i] = X[:,p[i]]
    end
    rbm.persistent_chain_hid = hid_means(rbm, rbm.persistent_chain_vis)

    use_persistent = false
    for itr=1:n_iter
        # Check to see if we can use persistence at this epoch
        use_persistent = itr>=persistent_start ? persistent : false

        tic()
        @showprogress 1 "Fitting Batches..." for i=1:n_batches
            batch = X[:, ((i-1)*batch_size + 1):min(i*batch_size, end)]
            batch = full(batch)

            fit_batch!(rbm, batch; persistent=use_persistent,
                                   n_gibbs=n_gibbs,
                                   weight_decay=weight_decay,
                                   decay_magnitude=decay_magnitude,
                                   lr=lr, approx=approx)

        end
        walltime_µs=(toq()/n_batches/N)*1e6

        UpdateMonitor!(rbm,ProgressMonitor,X,itr;bt=walltime_µs,validation=validation)
        ShowMonitor(rbm,ProgressMonitor,itr)
    end

    return rbm, ProgressMonitor
end


text:
# fbp/z-test.jl

using Test

@test ellipse_im(:test)
@test embed(:test)
@test image_geom_add_unitv(:test)
@test image_geom(:test)


text:
__precompile__()
module QDXML

using ParserCombinator

export XML

abstract XML

include("types.jl")
include("lexer.jl")
include("utils.jl")

parse_string(s) = parse_one(s, all) |> Document
parse_file(f) = readall(f) |> parse_string

"""
Parse the given filename, open file or string as XML. If the document doesn't start with doctype or a tag, please use `parse_string(s)`
"""
XML(s) = isa(s, AbstractString) && first(s) == '<' ? parse_string(s) : parse_file(s)

end # module