text:
cd(@__DIR__)
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
using Pkg
Pkg.activate(@__DIR__)
CI && Pkg.instantiate()

using DynamicalSystems
using Entropies, RecurrenceAnalysis, DelayEmbeddings, ChaosTools, DynamicalSystemsBase
using Neighborhood
using Documenter, PyPlot
using DocumenterTools: Themes
import Downloads

# %%
# download the themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    Downloads.download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

# %% Build docs
# Style for plots inside documentation
include("style.jl")

PyPlot.ioff()
cd(@__DIR__)
ENV["JULIA_DEBUG"] = "Documenter"

makedocs(
modules=[DynamicalSystems, ChaosTools, DynamicalSystemsBase,
         DelayEmbeddings, RecurrenceAnalysis, Entropies, Neighborhood],
doctest=false,
sitename= "DynamicalSystems.jl",
root = @__DIR__,
format = Documenter.HTML(
    prettyurls = CI,
    assets = [
        "assets/logo.ico",
        asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
        ],
    collapselevel = 1,
    ),
pages = [
    "Introduction" => "index.md",
    "Contents" => "contents.md",
    "Dynamical systems" => [
        "Dynamical System Definition" => "ds/general.md",
        "Predefined Dynamical Systems" => "ds/predefined.md",
        "Numerical Data" => "embedding/dataset.md",
    ],
    "Entropies" => [
        "Entropies & Probabilities" => "entropies/api.md",
        "Probabilities Estimators" => "entropies/estimators.md",
    ],
    "DelayEmbeddings" => [
        "Delay Coordinates Embedding" => "embedding/reconstruction.md",
        "Traditional Optimal Embedding" => "embedding/traditional.md",
        "Unified Optimal Embedding" => "embedding/unified.md",
    ],
    "ChaosTools" => [
       "Orbit Diagrams & PSOS" => "chaos/orbitdiagram.md",
       "Lyapunov Exponents" => "chaos/lyapunovs.md",
       "Detecting & Categorizing Chaos" => "chaos/chaos_detection.md",
       "Fractal Dimension" => "chaos/fractaldim.md",
       "Nonlinear Timeseries Analysis" => "chaos/nlts.md",
       "Periodicity & Ergodicity" => "chaos/periodicity.md",
       "Attractor Basins, Tipping Points" => "chaos/basins.md",
    ],
    "RecurrenceAnalysis" => [
        "Recurrence Plots" => "rqa/rplots.md",
        "Recurrence Quantification Analysis" => "rqa/quantification.md",
        "Windowed RQA" => "rqa/windowed.md",
        "Recurrence Networks" => "rqa/networks.md",
    ],
    "Advanced Documentation" => "advanced.md",
    "Contributor Guide" => "contributors_guide.md",
],
expandfirst = ["index.md"],
)

PyPlot.close("all")
PyPlot.ion()

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/DynamicalSystems.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
# Autogenerated wrapper script for SCS_jll for aarch64-apple-darwin
export libscsdir, libscsindir

using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("SCS")
JLLWrappers.@declare_library_product(libscsdir, "@rpath/libscsdir.dylib")
JLLWrappers.@declare_library_product(libscsindir, "@rpath/libscsindir.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libscsdir,
        "lib/libscsdir.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

rgb_matrix = fill(RGB(0., 0., 0.0), (6,))
img_1 = Raytracing.HDRimage(3, 2, rgb_matrix)
img_2 = Raytracing.HDRimage(3, 2)

@test img_1.width == 3
@test img_1.height == 2
@test img_1.rgb_m == img_2.rgb_m

# Controllo della scrittura errata nell'assert
@test_throws AssertionError img = Raytracing.HDRimage(3, 3, rgb_matrix)
@test_throws AssertionError img = Raytracing.HDRimage(1, 3, rgb_matrix)


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_UniformPigment" begin
     color = RGB(1.0, 2.0, 3.0)
     pigment = UniformPigment(color)

     get_color(pigment, Vec2d(0.0, 0.0)) ≈ color
     @test get_color(pigment, Vec2d(1.0, 0.0)) ≈ color
     @test get_color(pigment, Vec2d(0.0, 1.0)) ≈ color
     @test get_color(pigment, Vec2d(1.0, 1.0)) ≈ color
     @test get_color(pigment, Vec2d(0.5, 0.5)) ≈ color
end

@testset "test_ImagePigment" begin
     image = HDRimage(2, 2)
     Raytracing.set_pixel(image, 0, 0, RGB{Float32}(1.0, 2.0, 3.0))
     Raytracing.set_pixel(image, 1, 0, RGB{Float32}(2.0, 3.0, 1.0))
     Raytracing.set_pixel(image, 0, 1, RGB{Float32}(2.0, 1.0, 3.0))
     Raytracing.set_pixel(image, 1, 1, RGB{Float32}(3.0, 2.0, 1.0))

     pigment = ImagePigment(image)
     @test get_color(pigment, Vec2d(0.0, 0.0)) ≈ RGB{Float32}(1.0, 2.0, 3.0)
     @test get_color(pigment, Vec2d(1.0, 0.0)) ≈ RGB{Float32}(2.0, 3.0, 1.0)
     @test get_color(pigment, Vec2d(0.0, 1.0)) ≈ RGB{Float32}(2.0, 1.0, 3.0)
     @test get_color(pigment, Vec2d(1.0, 1.0)) ≈ RGB{Float32}(3.0, 2.0, 1.0)
end

@testset "test_CheckeredPigment" begin
     color1 = RGB{Float32}(1.0, 2.0, 3.0)
     color2 = RGB{Float32}(10.0, 20.0, 30.0)

     pigment = CheckeredPigment(color1, color2, 2)

     # With num_of_steps == 2, the pattern should be the following:
     #
     #              (0.5, 0)
     #   (0, 0) +------+------+ (1, 0)
     #          |      |      |
     #          | col1 | col2 |
     #          |      |      |
     # (0, 0.5) +------+------+ (1, 0.5)
     #          |      |      |
     #          | col2 | col1 |
     #          |      |      |
     #   (0, 1) +------+------+ (1, 1)
     #              (0.5, 1)
     @test get_color(pigment, Vec2d(0.25, 0.25)) ≈ color1
     @test get_color(pigment, Vec2d(0.75, 0.25)) ≈ color2
     @test get_color(pigment, Vec2d(0.25, 0.75)) ≈ color2
     @test get_color(pigment, Vec2d(0.75, 0.75)) ≈ color1
end


text:
"""
A note on the thinking behind this module:

we want to make it as simple as possible to optimise stuff, thats the goal here
in order to do that is it worth adding in some helper structs that store bits
of information that are useful?

It's optional, we use multiple dispatch so that familiar users can do it manually
but if you want something fast then we handle everything internally?

Previously the big problem was trying to do too much. We should make it play
nicely with something like pulse_sim.jl but the two should be separate. Pulse_sim
is for the simulation of pulses and pulse sequences efficiently.

This is for numerical optimal control and not that. Definitely need to have interop though
    so that I can keep writing pulse sequences there but use this to optimise them
    as we see fit in future.
"""

#abstract type Problem end

# abstract type ClosedSystem <: Problem end

"""
Store some pieces of information about your problem in this handy
and barely typed struct. We need to keep things flexible but also pleasant to use
"""
struct ProblemInfo
    # dimensions of your quantum system, useful to know to verify Hams make sense
    dim::Int
    # whether you've tried to optimise it
    opt::Bool
    # result of the optimisation attempt, useful for obvious reasons
    res

    # nervous about including these, maybe its better that the user really defines
    # all this stuff themselves. Don't want to over engineer things
    # drift Hamiltonian/Hamiltonians
    H_drift
    # control Hamiltonian/Hamiltonians
    H_ctrl
    # # I would like to keep a space for a full time dependent Hamiltonian, so that we
    # # can integrate stuff nicely!
    # full_Hamiltonian
    # timestep? would make it possible to use integrator quickly
    timestep
    # total time
    T
    # functional, here we can assume that things work
    functional
end

"""
I'm not sure that this is actually the right way to go, unless we copy the values
out of this struct and into the code directly and just use this to sort of save
    the state of the code as it stands.
"""


text:
# Practical Machine learning
# Reinforcement learning - Q learning example
# Chapter 12

immutable Experience
    s0::NNMatrix
    a0::Int64
    r0::Float64
    s1::NNMatrix
end

type DQN
    matrices::Array{NNMatrix,1}
    numStates::Int64 # length of the state vector (input vector)
    numHidden::Int64 # number of hidden nodes
    numActions::Int64 # number of actions (output vector)
    w1::NNMatrix
    b1::NNMatrix
    w2::NNMatrix
    b2::NNMatrix
    gamma::Float64 # reward discount
    epsilon::Float64 # epsilon-greedy policy
    alpha::Float64 # learning rate
    errorClamp::Float64
    expSize::Int64
    expAddProb::Float64
    expLearn::Int64
    expCnt::Int64
    experiences::Vector{Experience}
    s0::NNMatrix
    a0::Int64
    r0::Float64
    s1::NNMatrix
    a1::Int64
    lastG::Graph
    solver::Solver
    function DQN(numStates::Int64, numHidden::Int64, numActions::Int64;
                 std=0.02, gamma=0.75, epsilon=0.1, alpha=0.00001, errorClamp=5.0,
                 expSize=5000, expAddProb=0.05, expLearn=5)

        matrices = Array(NNMatrix, 0) # reference to matrices used by solver
        w1 = randNNMat(numHidden,  numStates, std);  push!(matrices, w1)
        b1 = randNNMat(numHidden,          1, std);  push!(matrices, b1)
        w2 = randNNMat(numActions, numHidden, std);  push!(matrices, w2)
        b2 = randNNMat(numActions,         1, std);  push!(matrices, b2)

        new(matrices, numStates, numHidden, numActions, w1, b1, w2, b2,
            gamma, epsilon, alpha, errorClamp, expSize, expAddProb, expLearn,0, Array(Experience,0),
            NNMatrix(numStates, 1), 0, typemin(Float64), NNMatrix(numStates, 1), 0, Graph(),Solver())
    end
end

function forward(m::DQN, s::NNMatrix, doBP::Bool=false)
    g = Graph(doBP)
    h0 =  add(g, mul(g, m.w1, s), m.b1)
    hd = tanh(g, h0) # hidden state
    a = add(g, mul(g, m.w2, hd), m.b2) # action vector
    m.lastG = g
    return a
end

function act(m::DQN, s::NNMatrix)
    a = rand()<=m.epsilon? rand(1:m.numActions):indmax(forward(m, s, false).w)
    return a # return selected action
end

function learnFromTuple(m::DQN, s0::NNMatrix, a0::Int64, r0::Float64, s1::NNMatrix)
    tmat = forward(m,s1,false)
    qmax = r0 + m.gamma * tmat.w[indmax(tmat.w)]

    pred = forward(m, s0, true)
    tdError = pred.w[a0,1] - qmax
#     println("tmat=$(tmat.w) max=$(indmax(tmat.w)) q=$qmax r0=$(r0) tMax=$(tmat.w[indmax(tmat.w)]) a0=$a0 pred=$(pred.w) tdError=$tdError")

    tdErrorClamp = minimum([maximum([tdError,-m.errorClamp]),m.errorClamp]) # huber loss to robustify
    pred.dw[a0,1] = tdErrorClamp
    backprop(m.lastG)
    solverstats = step(m.solver, m.matrices, m.alpha, 1e-06, m.errorClamp)

#     for k = 1:length(m.matrices)
#         @inbounds mat = m.matrices[k] # mat ref
#         @inbounds for j = 1:mat.d, i = 1:mat.n
#             mat.w[i,j] += - m.alpha * mat.dw[i]
#             mat.dw[i,j] = 0
#         end
#     end
    return tdErrorClamp
end

function learn(m::DQN, s0::NNMatrix, a0::Int64, r0::Float64, s1::NNMatrix)

    tdError = learnFromTuple(m, s0, a0, r0, s1)

    if rand() <= m.expAddProb
        m.expCnt = m.expCnt >= m.expSize? 1 : m.expCnt + 1
        if length(m.experiences) < m.expSize
            push!(m.experiences,Experience(s0, a0, r0, s1))
        else
            m.experiences[m.expCnt] = Experience(s0, a0, r0, s1)
        end
    end

    #　memory replay
    len = length(m.experiences)
    for i = 1:min(len,m.expLearn)
        r = rand(1:len)
        exp = m.experiences[r]
        err = learnFromTuple(m, exp.s0, exp.a0, exp.r0, exp.s1)
    end

    return tdError
end


text:
# Autogenerated wrapper script for Gdbm_jll for x86_64-apple-darwin
export gdbm_dump, gdbm_load, gdbmtool, libgdbm

using Libiconv_jll
using Readline_jll
JLLWrappers.@generate_wrapper_header("Gdbm")
JLLWrappers.@declare_executable_product(gdbm_dump)
JLLWrappers.@declare_executable_product(gdbm_load)
JLLWrappers.@declare_executable_product(gdbmtool)
JLLWrappers.@declare_library_product(libgdbm, "@rpath/libgdbm.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Readline_jll)
    JLLWrappers.@init_executable_product(
        gdbm_dump,
        "bin/gdbm_dump",
    )

    JLLWrappers.@init_executable_product(
        gdbm_load,
        "bin/gdbm_load",
    )

    JLLWrappers.@init_executable_product(
        gdbmtool,
        "bin/gdbmtool",
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CommutativeRings
using Test
using Polynomials

if !isdefined(Polynomials, :Poly)
    Poly = Polynomials.Polynomial
end

ismonic = CommutativeRings.ismonic

import Base: ==
function ==(a::UnivariatePolynomial, b::Poly)
    va = getproperty.(a.coeff, :val)
    vb = coeffs(b)
    n = length(vb)
    while n > 0 && vb[n] == 0
        n -= 1
    end
    va == vb[1:n]
end

let S = ZZ{Int}, P = UnivariatePolynomial{S,:x}, PQ = UnivariatePolynomial{QQ{Int},:x}

x = P([0, 1])
CP = (Int[], [1], [0, 0, 4], [2, 1], [1,0,30])

@testset "constructors" begin
    @test basetype(P) == S
    @test depth(P) == 2
    @test P([S(0), S(1)]).coeff[2] == S(1)
    @test P([1]).coeff[1] == S(1)
    @test length(P(Int[]).coeff) == 0
    @test length(P(Int[0]).coeff) == 0
    @test eltype(UnivariatePolynomial(:x, S[]).coeff) == S
    @test iscoprime(x^2 + 1, x + 1)
    @test !iscoprime(4x^2 + 4x + 1, 2x + 1)

    @test typeof(UnivariatePolynomial(:x, ZZ(1))) == P
    co = [1,2,3]
    @test copy(P(co)) == P(co)
    @test copy(P(co)).coeff !== P(co).coeff

    @test isunit(P([-1]))
    @test isunit(P([1]))
    @test !isunit(P([0]))
    @test length(zero(P).coeff) == 0
    z = zero(P)
    o = one(P)
    @test length(o.coeff) == 1 && o.coeff[1] == S(1)
    @test iszero(z)
    @test !iszero(o)
    @test isone(o)
    @test !isone(z)

    @test P([1,2,3]) == P([1,2,3])
    @test P([1,2]) != P([1])

    @test P(1) == 1
    @test P(ZZ(Int8(1))) == 1
    @test 1 + P(1) == P(2)
    @test ZZ(1) + P(1) == P(2)
    @test ZZ(1) - P(1) == z
    @test P(1) == UnivariatePolynomial{ZZ{Int8},:y}(1)
    @test P([0,1]) != UnivariatePolynomial{ZZ{Int}, :y}([0,1])
    @test P(1) == UnivariatePolynomial{ZZ{Int},:x}(1)
    @test UnivariatePolynomial{ZZ{Int8},:x}([1]) == P([1])
    @test hash(UnivariatePolynomial{ZZ{Int8},:x}([1])) == hash(P([1]))
end

@testset "varnames and generators" for P in (S[:x],)
    x = monom(P)
    @test [x] == generators(P)
    @test varnames(P) == [:x]
    @test varnames(P(0)) == [:x]

    p = P(S[1,2,3])
    @test deg(p) == 2
    @test multideg(p) == [2]
    @test LC(p) == S(3)
    @test LM(p) == x^2
    @test LT(p) == 3*x^2
end

@testset "promotion of types" begin
    @test promote_type(P, PQ) == PQ
    @test promote_type(P, ZZ{Int}) == P
    @test promote_type(P, Int) == P
    @test promote_type(P, Rational{Int8}) == UnivariatePolynomial{QQ{Int},:x}
    @test promote_rule(P, UnivariatePolynomial{S,:y}) == Base.Bottom
end

@testset "operation $cp $(string(op)) $cq" for op in (+, -, *), cp in CP, cq in CP
    p = P(S.(cp))
    q = P(S.(cq))
    @test op(p, q) == op(Poly(cp), Poly(cq))
end

@testset "operation unary - and ^ $cp" for cp in CP
    p = P(S.(cp))
    @test -p == -Poly(cp)
    @test p^10 == Poly(cp)^10
    @test p^0 == one(P)
    !isunit(p) && @test_throws DomainError p^-1
    @test zero(p)^0 == one(P)
    @test one(p)^-1 == one(P)
    @test zero(p)^0 == one(P)
end

hasunitlead(p::UnivariatePolynomial) = !iszero(p) && isunit(p.coeff[end])

@testset "operation divrem($cp,$cq)" for cp in CP, cq in CP
    p = P(S.(cp))
    q = P(S.(cq))
    if hasunitlead(q) || deg(p) < deg(q) || ( p == q && !iszero(q))
        @test divrem(p, q) == divrem(Poly(cp), Poly(cq))
        @test deg(rem(p, q)) < deg(q)
    elseif !iszero(q)
        @test deg(rem(p, q)) >= deg(q)
    else
        @test_throws DomainError divrem(p, q)
    end
end

@testset "gcd calculations" begin

    P = UnivariatePolynomial{ZZ{Int},:x}
    @test P([0, 1]) == monom(P)
    x = P([0,1])
    @test deg(x + 2x^2 + 1 +0*x^3) == 2
    p = x^8 + x^6 -3x^4 - 3x^3 + 8x^2 +2x - 5
    q = 21 + 3x^6 + 5x^4 - 4*x^2 -9x
    @test deg(pgcd(p, q)) == 0
    r1 = pgcd(p, q)
    s = (ZZ(1) + 3x + x^2)^5
    @test pgcd(p * s, q * s) / r1 == s

    p = P([1, 2, 1])
    q = P([1, 1])
    @test gcd(p, q) == q
    @test gcd([p, p, q]) == q
    @test gcdx(p, q) == (q, P(Int[]), P([1]))
    @test lcm(p, q) == p
    @test lcm([p, q, p]) == p

    @test UnivariatePolynomial{ZZ{Int},:x}(p) === p
    @test UnivariatePolynomial{ZZ{Int32},:x}(p) == p
    @test div(p, q) == q
    @test content(p) == ZZ(1)
    @test primpart(p) == p
    @test content(-p) == ZZ(-1)
    @test primpart(-p) == p

    Q = UnivariatePolynomial{QQ{Int},:x}
    qp = Q(-p) / 12
    @test primpart(qp) == Q(p)
    @test content(qp) == -1//12

    @test iszero(LC(zero(P)))
    @test isone(LC(one(P)))
    @test LC(p) == 1
    @test multideg(one(P)) == [0]
    @test multideg(p) == [2]

    PP = UnivariatePolynomial{P,:z}
    pp = PP([p, q, s])

    # call show
    @test sprint(show, p) !== nothing
    @test sprint(show, zero(P)) !== nothing
    @test sprint(show, one(P)) !== nothing
    @test sprint(show, x^3 + x) !== nothing
end

@testset "pseudo gcd" begin
    S = ZZ{BigInt}
    P = UnivariatePolynomial{S,:x}
    p = P([1, 4, 5, 1, 6, 0, 3])
    q = P([5, 1, 1, 3])
    x = P([0, 1])
    @test content(2p) == S(2)
    @test primpart(12p) == p
    @test LC(p) == S(3)
    a, b, f = pdivrem(p, q)
    @test f * p == a * q + b
    a, b, f = pdivrem(p, P([2]))
    @test f * p == a * P([2]) + b

    @test gcd(p, q) == one(P)
    @test pgcd(p, q) == P([1])
    @test_throws DomainError gcdx(p, q)
    g, u, v, f = pgcdx(p, q)
    @test iszero(p * u + q * v - g * f)
    @test deg(g) == 0
    p = 3x^10 - 3x^9 - 3x^8 + 3x^7 - 2x^5 - 3x^4 - 2x^3 - 1x - 3
    q = 2x^4 + 3x^3 - 1x^2 + x - 2
    g, u, v, f = pgcdx(q, p)
    @test iszero(q * u + p * v - g * f)
    @test deg(g) == 0
    s = 3x^3 - 2x^2 - ZZ(1)
    p *= s
    q *= s
    g, u, v, f = pgcdx(p, q)
    @test iszero(p * u + q * v - g * f)
    @test isdiv(g, s)

    S = ZZmod{31,Int}
    P = UnivariatePolynomial{S,:x}
    x = P([0, 1])
    p = x^3 + 3x^2 - 4
    q = 4 - (-x)
    g, u, v, f = pgcdx(p, q)
    @test iszero(p * u + q * v - g * f)
    p2 = 2p
    @test p2 / 2 == p

    p = 12x^2 + 45x + 60
    q = 5x^2 + 1
    g, u, v, f = pgcdx(p, q)
    @test iszero(p * u + q * v - g * f)

    qq = [q, P([2,1])]
    a, r = divrem(p, qq)
    @test sum(a .* qq) + r == p

    x = UnivariatePolynomial{ZZ{Int},:x}([0, 1])
    @test CommutativeRings.invert(x, x^2 + 1) == -x
    @test_throws DomainError CommutativeRings.invert(x + 1, x^2 + 1)

    @test CommutativeRings.issimple(1.0)
    @test CommutativeRings.issimple(x^3)
    @test !CommutativeRings.issimple(x+1)
end

@testset "polynomials of other structures" begin

    Zp = ZZmod{17, Int8}
    P = Zp[:x]
    @test P([1,2]) != Zp[:y]([1,2])
    x = P([0, 1])
    p = (3x^2 + 2x + 1) * (x - 1)
    q = x^2 - x
    @test p / (x-1) == (3x^2 + 2x + 1)
    @test gcd(p, q) == x - 1

    @test ismonom(P([0,2]))
    @test !ismonom(P([1,2]))
    @test ismonic(P([2,1]))
    @test !ismonic(P([2,-1]))

end

@testset "apply polynomial" begin
    x = monom(ZZ{Int}[:x], 1)
    y = monom(ZZ{Int}[:y], 1)
    @test x(1) == 1
    @test x(0) == 0
    @test (y^2)(x+1) == x^2 + 2x + 1
    @test (0*y)(x^2) == 0
end

@testset "abstract types" begin
    @test QQ[:x]([1//2, 2//3]) == QQ{Int}[:x]([1//2, 2//3])
end

@testset "GF[:x]" begin
    G = GF(7, 3)
    GX = G[:x]
    x = monom(GX)
    p = x^3 + G[(7*2+1)*7+5]x + 1
    q = x^3 + 1 + G[200]*x
    @test isirreducible(p)
    @test factor(q) |> length == 3
    @test deg(p^2) == 6
    @test p * q == q * p
    @test p + q - p == q
    @test q / (x + G[28]) |> deg == 2
end

x = monom(ZZ{BigInt}[:x])
@testset "resultant($a,$b)" for (a, b, c) in ((0,0,0), (0,1,0), (0, x,0), (3,0,0),
                                       (3,2,1), (3,x,3), (3,5x^2+1,3^2), (7,5x^7-1,7^7),
                                       (2x^32+1, 3*x^11+13, 906811932214969750127235270540901107329),
                                       (2x^31+1, 3x^11+1,-617673396281899))

    @test resultant(a, b) == c
    @test resultant(b, a) == c * (-1)^(deg(a)*deg(b)+2)
    @test resultant(11a, 13b) == (a*b == 0 ? 0 : c * big(11)^deg(b) * big(13)^deg(a))

    ref = CommutativeRings.resultant_naive
    @test a*b == 0 || ref(a*x^0, b*x^0) == c
    @test a*b == 0 || ref(b*x^0, a*x^0) == c * (-1)^(deg(a)*deg(b)+2)
end

end


text:
module MuseInference

using AbstractDifferentiation
using Base.Iterators: peel, repeated
using ComponentArrays
using CovarianceEstimation
using Dates
using Distributed
using Distributions
using FileIO
using FiniteDifferences
using ForwardDiff
using LinearAlgebra
using Markdown
using NamedTupleTools
using Optim
using ProgressMeter
using Random
using Requires
using Setfield
using Statistics
using UnPack

export MuseProblem, MuseResult, muse, muse!, get_J!, get_H!

include("util.jl")
include("ad.jl")
include("progress.jl")
include("muse.jl")
include("manual.jl")
@init @require Turing="fce5fe82-541a-59a6-adf8-730c64b5f9a0" include("turing.jl")

end


text:
#=ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля
   +учитываются внутренние прямоугольные перегородки,которые изолированы друг от друга и от внешней перегородки

   РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так:
   нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке,
   следующий - за исключением двух последних клеток на Востоке, и т.д.
=#
function mountain_but_harder(r::Robot)
    steps=0
    num_steps=[]
    while isborder(r,Sud)==false && isborder(r,West)==false
        push!(num_steps, get_num_movements!(r, West))
        push!(num_steps, get_num_movements!(r, Sud))
        steps+=2
    end
    putmarker!(r)
    putmarkers!(r,Ost)
    for_width=get_num_movements!(r, West)#=ширина поля=#
    while isborder(r,Nord)==false && for_width>0
        width=for_width #фиксированная величина(ширина) границ,доступных для движения
        move!(r,Nord)
        while width>0
            putmarker!(r)
            width-=build_trajectory(r,Ost)
        end
        for_marking_trajectory(r,West)
        for_width-=1  #уменьшаем на 1 ширину линии(траектория относительно ширины)
    end
    get_num_movements!(r,Sud)
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for k in 1:num_steps[steps]
            move!(r,side)
        end
        steps-=1
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end
function build_trajectory(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==true && isborder(r,next(side))==false #цикл для схода с линии маркеровки(для обхода перегородок)
        move!(r,next(side))
        num_steps+=1
    end
    count=0
    if isborder(r,side)==false #считаем клетки которые НЕ НУЖНО обходить
        move!(r,side)
        count+=1
    end
    if num_steps !=0 # true если встречались перегородки
        while isborder(r,inverse(next(side)))==true
            move!(r,side)
            count+=1
        end
        for a in 1:num_steps
            move!(r,inverse(next(side)))
        end
    end
    return count
end

function for_marking_trajectory(r::Robot, side::HorizonSide)
    count=1
    while count!=0
        count=build_trajectory(r,side)
    end
end

function get_num_movements!(r::Robot, side::HorizonSide)
    num_steps = 0
    while isborder(r,side)==false
        move!(r,side)
        num_steps += 1
    end
    return num_steps
end

next(side::HorizonSide)=HorizonSide(mod(Int(side)+1,4))

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
#RegressionTypes = (#=glm=Val{:glm}(), =#cholesky=Val{:cholesky}(), qr=Val{:qr}())
reg(X::AbstractMatrix{T},y::Vector{T}, ::Val{:glm}) where T = GLM.coef(lm(X, y))
reg(X::AbstractMatrix{T},y::AbstractArray{T}, ::Val{:cholesky}
  ) where T = cholesky!(Symmetric(X'*X))\ (X' * y)
reg(X::AbstractMatrix{T},y::Vector{T}, ::Val{:qr}) where T = qr(X)\y
reg(X::AbstractMatrix{T}, y::AbstractArray{T}, ::Val{:svd}) where T = svd(X)\y
reg(X::CuMatrix{T}, y::CuArray{T}, ::Val{:cholesky}) where T = hyperreg(X,y) |> CUDA.vec
reg(X::AbstractMatrix{T}, y::AbstractArray{T}, ::Val{:choleskyzygote}) where T = cholesky(X'*X)\ (X' * y)
reg(X::CuMatrix{T}, y::AnyCuVecOrMat{T}, ::Val{:choleskyzygote}) where T = CUDA.cholesky(X'*X)\ (X' * y)
#reg(X::CuMatrix{T}, y::CuVector{T}, ::Val{:none}) where {TV, T} = ones(T, size(X,2)) .* 0.75 |> CuVector{T}
#reg(X::AbstractMatrix{T}, y::AbstractArray{T}, ::Val{:none}) where {TV, T} = ones(T, size(X,2)) .* 0.75
#these methods run with an intercept, which is thrown out and not returned.
reg(X::AbstractMatrix{T}, y::AbstractArray{T},::Val{:intercept}) where T = throw(
  "intercept regressions should not be directly called")
reg(X::AbstractMatrix{T}, y::AbstractArray{T},::Val{:interceptzygote}) where T = throw(
  "intercept regressions should not be directly called")

#=reg(X::AbstractMatrix{T}, y::AbstractArray{T},::Val{:intercept}) where T = reg(
  hcat(ones(T, length(y)), X), y, Val{:cholesky}())[2:end]
reg(X::CuMatrix{T}, y::CuArray{T}, ::Val{:intercept}) where T = reg(
  CUDA.hcat(CUDA.ones(T, length(y)), X), y, Val{:cholesky}())[2:end]
reg(X::AbstractMatrix{T}, y::CuArray{T}, ::Val{:interceptzygote}) where T = reg(
  CUDA.hcat(zygotecuones(T, length(y)), X), y, Val{:choleskyzygote}())[2:end]=#

function reg(X::Matrix{T}, y::Array{T}, ::Val{:fallbackreg}) where T

  #@eval Main X=$X|>deepcopy
  #@eval Main y=$y|>deepcopy
  #= NOTE: Let USVt be the SVD of X, . Note by propertys of svd, Ut*U=Vt*V=I and S  is diagonal
  Then we seek to solve:
    USVt*β = y
    SVt*β = U'y
    β = V*S^-1*U'y = V*(diag(S^-1) .* U' * y)
  =#
  #
  #sX = svd(X)
  #β = sX.V * (sX.S .^-1 .* (sX.U' * y))

  return pinv(X'*X, rtol=sqrt(eps(T)))*(X'*y)
end
#these are just conversion methods for the fallback
reg(X::TX, y::Ty, RegressionType::Val{:fallbackreg}, ::Type{T} = eltype(X) #fallback gpu method
  ) where {TX,Ty,T} = reg(X|>Matrix, y|>Array{T}, RegressionType) |> Ty
reg(X::Matrix{T}, y::Array{T}, RegressionType::Val{:zygotefallbackreg}
  ) where T = reg(X,y,Val{:fallbackreg}())
reg(X::TX, y::Ty, ::Val{:zygotefallbackreg}, ::Type{T} = eltype(X) #fallback gpu method
  ) where {TX,Ty,T} = reg(X|>zygotedevice2host,y|>zygotedevice2host,Val{:fallbackreg}()) |> Ty

#reg(X::AbstractMatrix, y::AbstractArray, ::Val{:qrzygote}) = qr(X)\y

#not sure if there is a potential gpu benefit here- probably would need a very custom kernal
#NOTE: fairly sure these are not needed. If they are needed, just uncomment
#const MIN_MAG = .01
#const MAX_MAG = 3.0
#nonzero(x::T,y) where T = ifelse(isfinite(x) && abs(x) < MIN_MAG, y, x)

#this sets A₁ to 1.0 and all other values to the amount implied from G
#of course, this destroys any current values of A
function initializeÃandG!(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T}) where {TM,TV,T}
  Θ.Ã[:,1] .= T(1.0)

  #below is better but it needs to sync with Zygote
  #someones::TM = ones(size(Θ.A,1),1) |> TM

  #copyto!(Θ.Ã, cumprod(hcat(someones, Θ.G),dims=2))

  #below is a bandaid to sync up the initialization with Zygote
  Θ.Ã[:,2:end] .= cumprodbyrow(Θ.G)

  return nothing
end

#this factors X. NOTE assumes we are solving for A₁ AND that A₁ .== 1
function initializeÃandfactor(Θ::AbstractVolumePartsIter{TM, TV,T},
  Xv::AbstractXYIter) where {TM, TV,T}

  @unpack Ã, LÃ = Θ.expand
  @unpack ws, Rtws, RLws= Xv

  initializeÃandG!(Θ, Xv)
  @assert all(Θ.Ã[:,1] .≈ T(1.0)) #always true if Ã is initialized

  factored = genabsμₖ.(Ã', LÃ', ws, Rtws, RLws)

  return factored
end

#this factors X but assumes we are solving for A₁ AND that A₁ .== 1
function initializeÃandfactor(Θ::VolumePartsIter,
  Xv::AbstractXYIter{TM,TV,T}) where {TM<:CuMatrix, TV<:CuVector, T}

  @unpack Ã, LÃ = Θ.expand
  @unpack ws, Rtws, RLws= Xv

  initializeÃandG!(Θ, Xv)
  @assert all(Θ.Ã[:,1] .≈ T(1.0)) #always true if Ã is initialized
  #first map each coefficient to its appropriate

  factored = genabsμₖ.(Ã', LÃ', ws, Rtws, RLws)

  return factored
end

#updates A₁ for a single pass by running the appropriate regression
#returns the loss function
#WARNING- we do not update all of A here after running the regression
function updateA₁!(loss, Θ::VolumePartsIter,
  Xv::AbstractXYIter{TM, TV,T},
  RegressionType::Val) where {TM, TV,T}

  local λ
  #try
    λ = loss(Θ, Xv, RegressionType)
  #=catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage) > 1000 ? errmessage[1:1000] : errmessage
    WARNING- I believe the below code is obsolete, so it is commented out
    @warn "Regression failed with error $errmessage\nAttempting fallback method"
    λ = loss(Θ, Xv, Val{:fallbackreg}())
  #@assert typeof(λ) <: Float64
  #@info "update: λ=$λ"
end=#
  return λ
end

  #set the initial value for g from the growth type
  #called for the inital vector at the beginning of optimization
function initializeΘgvector(g::TV, ::Type{T}, GrowthType::Val, G::Union{TV,Nothing}=nothing;
  initializeto1::Bool=true) where {TV, T}
  _initializeΘg(g::AbstractVector, ::Val{:identity}) = g |> Vector |> deepcopy
  _initializeΘg(g::AbstractVector, ::Val{:log}) = g |> Vector |> deepcopy .|> log
  function _initializeΘg(g::AbstractVector, ::Val{:twosidedlog},
    lowerbound::T=PARAM[:optimlowerbound], twosidedloglower::T=PARAM[:itertwosidedloglower],
    upperbound::T=PARAM[:optimupperbound], twosidedlogupper::T=PARAM[:itertwosidedlogupper])
    gout = g |> Vector |> deepcopy
    Δbound = upperbound
    Δtwosidedlog = twosidedlogupper - twosidedloglower
    gout .= -(twosidedloglower)*Δbound/Δtwosidedlog #solved to set G=1
    gout
  end

  if initializeto1
    g .= T(1.0)
  end
  gidentity = g |> Vector |> deepcopy #use the pretransformation version for verification later
  gout = _initializeΘg(g, GrowthType)
  #g .= gout |> TV
  copyto!(g, gout)
  @assert all(conditiongrowth(gout, GrowthType) |> Vector .≈ gidentity)
  return gout
end

conditiongrowth(g, ::Val{:identity})= g
conditiongrowth(g, ::Val{:log}) = exp.(g)
conditiongrowth(g::CuVector, ::Val{:log}) = CUDA.exp.(g)
function conditiongrowth(g::AbstractVector{T}, ::Val{:twosidedlog},
  lowerbound::T=PARAM[:optimlowerbound], twosidedloglower::T=PARAM[:itertwosidedloglower],
  upperbound::T=PARAM[:optimupperbound], twosidedlogupper::T=PARAM[:itertwosidedlogupper]
  ) where T

  (upperbound + lowerbound ≈ 0.0) || throw("bounds of twosided log must be symmetric about 0")
  Δbound = upperbound
  Δtwosidedlog = twosidedlogupper - twosidedloglower
  function maptolog(x) #maps the input to the two sided log of g e.g. [-10,10]->±exp([-2.5,2.5])
    normalizedabsx::T = abs(x)/Δbound
    mappedabsx::T = exp(normalizedabsx*Δtwosidedlog + twosidedloglower)
    return (signsans0(x)) * mappedabsx
  end

  return maptolog.(g)
end

function updateA₁!(x::TV, loss, Θ::VolumePartsIter,
  Xv::AbstractXYIter{TM, TV,T},
  RegressionType::Val, ::Val{:identity}
  ) where {TM, TV,T}

  #no need to copy if they are equivelent
  (Θ.g===x) && return throw("I disallowed Θ.g===x due to race error potential
    (It shouldn't make a big difference anyway)")
  #@info "updateA₁!: growth type identity"
  copyto!(Θ.g,x)
  return updateA₁!(loss,Θ,Xv,RegressionType)
end

function updateA₁!(x::TV, loss, Θ::VolumePartsIter,
  Xv::AbstractXYIter{TM, TV,T},
  RegressionType::Val, GrowthType::Val;
  ) where {TM, TV,T}

  (Θ.g===x) && throw("x should be a different vector if solving for GrowthType $GrowthType")
  #no need to copy if they are equivelent
  copyto!(Θ.g, conditiongrowth(x, GrowthType))
  #@info "updateA₁!: growth type identity"
  return updateA₁!(loss,Θ,Xv,RegressionType)
end

#handles the CUDA case by converting the supplied vector to a CUDA vector
updateA₁!(x::Vector, loss, Θ::VolumePartsIter,
  Xv::AbstractXYIter{TM, TV,T},
  RegressionType::Val, GrowthType::Val;
  ) where {TM<:CuMatrix, TV<:CuVector,T} = updateA₁!(
  x |> TV, loss,Θ,Xv,RegressionType,GrowthType)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module TestImports

using Test

module Indirection
import PublicAPI
end

module StrictImport
using ..Indirection
Indirection.PublicAPI.@strict import PublicAPI
end

function test_strict_import()
    @test !isdefined(StrictImport, Symbol("@public"))
    @test !isdefined(StrictImport, Symbol("@strict"))
    @test !isdefined(StrictImport, :of)
    @test isdefined(StrictImport.PublicAPI, Symbol("@public"))
    @test isdefined(StrictImport.PublicAPI, Symbol("@strict"))
    @test isdefined(StrictImport.PublicAPI, :of)
    @test !isdefined(StrictImport.PublicAPI, :Internal)
end

end  # module


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPUmfpack

import ..MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    SubVector, StrideOneVector, SparseMatrixCSC, get_tril_to_full,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    AbstractOptions, AbstractLinearSolver, set_options!, UMFPACK,
    introduce, factorize!, solve!, mul!, improve!, is_inertia, inertia

const INPUT_MATRIX_TYPE = :csc

const umfpack_default_ctrl = copy(UMFPACK.umf_ctrl)
const umfpack_default_info = copy(UMFPACK.umf_info)

@kwdef mutable struct Options <: AbstractOptions
    umfpack_pivtol::Float64 = 1e-4
    umfpack_pivtolmax::Float64 = 1e-1
    umfpack_sym_pivtol::Float64 = 1e-3
    umfpack_block_size::Float64 = 16
    umfpack_strategy::Float64 = 2.
end

mutable struct Solver <: AbstractLinearSolver
    inner::UMFPACK.UmfpackLU
    tril::SparseMatrixCSC{Float64}
    full::SparseMatrixCSC{Float64}
    tril_to_full_view::SubVector{Float64}

    p::Vector{Float64}

    tmp::Vector{Ptr{Cvoid}}
    ctrl::Vector{Float64}
    info::Vector{Float64}

    opt::Options
    logger::Logger
end

umfpack_di_numeric(colptr::StrideOneVector{Int32},rowval::StrideOneVector{Int32},
                   nzval::StrideOneVector{Float64},symbolic::Ptr{Nothing},
                   tmp::Vector{Ptr{Nothing}},ctrl::Vector{Float64},
                   info::Vector{Float64}) = ccall(
                       (:umfpack_di_numeric,:libumfpack),
                       Int32,
                       (Ptr{Int32},Ptr{Int32},Ptr{Float64},Ptr{Cvoid},Ptr{Cvoid},
                        Ptr{Float64},Ptr{Float64}),
                       colptr,rowval,nzval,symbolic,tmp,ctrl,info)
umfpack_di_solve(typ,colptr,rowval,nzval,x,b,numeric,ctrl,info) = ccall(
    (:umfpack_di_solve,:libumfpack),
    Int32,
    (Int32, Ptr{Int32}, Ptr{Int32}, Ptr{Float64},Ptr{Float64},
     Ptr{Float64}, Ptr{Cvoid}, Ptr{Float64},Ptr{Float64}),
    typ,colptr,rowval,nzval,x,b,numeric,ctrl,info)

function Solver(csc::SparseMatrixCSC;
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs)

    p = Vector{Float64}(undef,csc.n)
    full,tril_to_full_view = get_tril_to_full(csc)

    full.colptr.-=1; full.rowval.-=1

    inner = UMFPACK.UmfpackLU(C_NULL,C_NULL,full.n,full.n,full.colptr,full.rowval,full.nzval,0)
    UMFPACK.finalizer(UMFPACK.umfpack_free_symbolic,inner)
    UMFPACK.umfpack_symbolic!(inner)
    ctrl = copy(umfpack_default_ctrl)
    info = copy(umfpack_default_info)
    ctrl[4]=opt.umfpack_pivtol
    ctrl[12]=opt.umfpack_sym_pivtol
    ctrl[5]=opt.umfpack_block_size
    ctrl[6]=opt.umfpack_strategy

    tmp = Vector{Ptr{Cvoid}}(undef, 1)

    return Solver(inner,csc,full,tril_to_full_view,p,tmp,ctrl,info,opt,logger)
end

function factorize!(M::Solver)
    UMFPACK.umfpack_free_numeric(M.inner)
    M.full.nzval.=M.tril_to_full_view
    status = umfpack_di_numeric(M.inner.colptr,M.inner.rowval,M.inner.nzval,M.inner.symbolic,M.tmp,M.ctrl,M.info)
    M.inner.numeric = M.tmp[]

    M.inner.status = status
    return M
end
function solve!(M::Solver,rhs::StrideOneVector{Float64})
    status = umfpack_di_solve(1,M.inner.colptr,M.inner.rowval,M.inner.nzval,M.p,rhs,M.inner.numeric,M.ctrl,M.info)
    rhs .= M.p
    return rhs
end
is_inertia(::Solver) = false
inertia(M::Solver) = throw(InertiaException())

function improve!(M::Solver)
    if M.ctrl[4] == M.opt.umfpack_pivtolmax
        @debug(M.logger,"improve quality failed.")
        return false
    end
    M.ctrl[4] = min(M.opt.umfpack_pivtolmax,M.ctrl[4]^.75)
    @debug(M.logger,"improved quality: pivtol = $(M.ctrl[4])")
    return true

    return false
end
introduce(::Solver)="umfpack"

end # module


text: muladd(x,muladd(x,muladd(x,muladd(x,5,4),3),2),1)

max_stars_repo_path: Module 2/Chapter04/muladd.jl

max_issues_repo_path: Module 2/Chapter04/muladd.jl

max_forks_repo_path: Module 2/Chapter04/muladd.jl


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil: @require
using ParallelStencil.FiniteDifferences3D
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @init_parallel_stencil($package, Float64, 3)
        @require @is_initialized()
        nx, ny, nz = 7, 5, 6
        A       =  @rand(nx  , ny  , nz  );
        Ax      =  @rand(nx+1, ny  , nz  );
        Ay      =  @rand(nx  , ny+1, nz  );
        Az      =  @rand(nx  , ny  , nz+1);
        Axy     =  @rand(nx+1, ny+1, nz  );
        Axz     =  @rand(nx+1, ny  , nz+1);
        Ayz     =  @rand(nx  , ny+1, nz+1);
        Axyz    =  @rand(nx+1, ny+1, nz+1);
        Axyzz   =  @rand(nx+1, ny+1, nz+2);
        Axyyz   =  @rand(nx+1, ny+2, nz+1);
        Axxyz   =  @rand(nx+2, ny+1, nz+1);
        Axx     =  @rand(nx+2, ny  , nz  );
        Ayy     =  @rand(nx  , ny+2, nz  );
        Azz     =  @rand(nx  , ny  , nz+2);
        Axxyy   =  @rand(nx+2, ny+2, nz  );
        Axxzz   =  @rand(nx+2, ny  , nz+2);
        Ayyzz   =  @rand(nx  , ny+2, nz+2);
        Axyyzz  =  @rand(nx+1, ny+2, nz+2);
        Axxyzz  =  @rand(nx+2, ny+1, nz+2);
        Axxyyz  =  @rand(nx+2, ny+2, nz+1);
        Axxyyzz =  @rand(nx+2, ny+2, nz+2);
        R       = @zeros(nx  , ny  , nz  );
        Rxxyyzz = @zeros(nx+2, ny+2, nz+2);
        @testset "1. compute macros" begin
            @testset "differences" begin
                @require @is_initialized()
                @parallel  d_xa!(R, Ax)      = (@all(R) = @d_xa(Ax); return)
                @parallel  d_ya!(R, Ay)      = (@all(R) = @d_ya(Ay); return)
                @parallel  d_za!(R, Az)      = (@all(R) = @d_za(Az); return)
                @parallel  d_xi!(R, Axyyzz)  = (@all(R) = @d_xi(Axyyzz); return)
                @parallel  d_yi!(R, Axxyzz)  = (@all(R) = @d_yi(Axxyzz); return)
                @parallel  d_zi!(R, Axxyyz)  = (@all(R) = @d_zi(Axxyyz); return)
                @parallel d2_xi!(R, Axxyyzz) = (@all(R) = @d2_xi(Axxyyzz); return)
                @parallel d2_yi!(R, Axxyyzz) = (@all(R) = @d2_yi(Axxyyzz); return)
                @parallel d2_zi!(R, Axxyyzz) = (@all(R) = @d2_zi(Axxyyzz); return)
                R.=0; @parallel  d_xa!(R, Ax);       @test all(R .== Ax[2:end,    :,    :].-Ax[1:end-1,      :,      :])
                R.=0; @parallel  d_ya!(R, Ay);       @test all(R .== Ay[    :,2:end,    :].-Ay[      :,1:end-1,      :])
                R.=0; @parallel  d_za!(R, Az);       @test all(R .== Az[    :,    :,2:end].-Az[      :,      :,1:end-1])
                R.=0; @parallel  d_xi!(R, Axyyzz);   @test all(R .== Axyyzz[2:end  ,2:end-1,2:end-1].-Axyyzz[1:end-1,2:end-1,2:end-1])
                R.=0; @parallel  d_yi!(R, Axxyzz);   @test all(R .== Axxyzz[2:end-1,2:end  ,2:end-1].-Axxyzz[2:end-1,1:end-1,2:end-1])
                R.=0; @parallel  d_zi!(R, Axxyyz);   @test all(R .== Axxyyz[2:end-1,2:end-1,2:end  ].-Axxyyz[2:end-1,2:end-1,1:end-1])
                R.=0; @parallel d2_xi!(R, Axxyyzz);  @test all(R .== (Axxyyzz[3:end,2:end-1,2:end-1].-Axxyyzz[2:end-1,2:end-1,2:end-1]).-(Axxyyzz[2:end-1,2:end-1,2:end-1].-Axxyyzz[1:end-2,2:end-1,2:end-1]))
                R.=0; @parallel d2_yi!(R, Axxyyzz);  @test all(R .== (Axxyyzz[2:end-1,3:end,2:end-1].-Axxyyzz[2:end-1,2:end-1,2:end-1]).-(Axxyyzz[2:end-1,2:end-1,2:end-1].-Axxyyzz[2:end-1,1:end-2,2:end-1]))
                R.=0; @parallel d2_zi!(R, Axxyyzz);  @test all(R .== (Axxyyzz[2:end-1,2:end-1,3:end].-Axxyyzz[2:end-1,2:end-1,2:end-1]).-(Axxyyzz[2:end-1,2:end-1,2:end-1].-Axxyyzz[2:end-1,2:end-1,1:end-2]))
            end;
            @testset "selection" begin
                @require @is_initialized()
                @parallel all!(R, A)        = (@all(R) = @all(A); return)
                @parallel inn!(R, Axxyyzz)  = (@all(R) = @inn(Axxyyzz); return)
                @parallel inn_x!(R, Axx)    = (@all(R) = @inn_x(Axx); return)
                @parallel inn_y!(R, Ayy)    = (@all(R) = @inn_y(Ayy); return)
                @parallel inn_z!(R, Azz)    = (@all(R) = @inn_z(Azz); return)
                @parallel inn_xy!(R, Axxyy) = (@all(R) = @inn_xy(Axxyy); return)
                @parallel inn_xz!(R, Axxzz) = (@all(R) = @inn_xz(Axxzz); return)
                @parallel inn_yz!(R, Ayyzz) = (@all(R) = @inn_yz(Ayyzz); return)
                R.=0; @parallel all!(R, A);         @test all(R .== A)
                R.=0; @parallel inn!(R, Axxyyzz);   @test all(R .== Axxyyzz[2:end-1,2:end-1,2:end-1])
                R.=0; @parallel inn_x!(R, Axx);     @test all(R .== Axx[2:end-1,      :,      :])
                R.=0; @parallel inn_y!(R, Ayy);     @test all(R .== Ayy[      :,2:end-1,      :])
                R.=0; @parallel inn_z!(R, Azz);     @test all(R .== Azz[      :,      :,2:end-1])
                R.=0; @parallel inn_xy!(R, Axxyy);  @test all(R .== Axxyy[2:end-1,2:end-1,      :])
                R.=0; @parallel inn_xz!(R, Axxzz);  @test all(R .== Axxzz[2:end-1,      :,2:end-1])
                R.=0; @parallel inn_yz!(R, Ayyzz);  @test all(R .== Ayyzz[      :,2:end-1,2:end-1])
            end;
            @testset "averages" begin
                @require @is_initialized()
                @parallel av!(R, Axyz)      = (@all(R) = @av(Axyz); return)
                @parallel av_xa!(R, Ax)     = (@all(R) = @av_xa(Ax); return)
                @parallel av_ya!(R, Ay)     = (@all(R) = @av_ya(Ay); return)
                @parallel av_za!(R, Az)     = (@all(R) = @av_za(Az); return)
                @parallel av_xi!(R, Axyyzz) = (@all(R) = @av_xi(Axyyzz); return)
                @parallel av_yi!(R, Axxyzz) = (@all(R) = @av_yi(Axxyzz); return)
                @parallel av_zi!(R, Axxyyz) = (@all(R) = @av_zi(Axxyyz); return)
                @parallel av_xya!(R, Axy)   = (@all(R) = @av_xya(Axy); return)
                @parallel av_xza!(R, Axz)   = (@all(R) = @av_xza(Axz); return)
                @parallel av_yza!(R, Ayz)   = (@all(R) = @av_yza(Ayz); return)
                @parallel av_xyi!(R, Axyzz) = (@all(R) = @av_xyi(Axyzz); return)
                @parallel av_xzi!(R, Axyyz) = (@all(R) = @av_xzi(Axyyz); return)
                @parallel av_yzi!(R, Axxyz) = (@all(R) = @av_yzi(Axxyz); return)
                R.=0; @parallel av!(R, Axyz);       @test all(R .== (Axyz[1:end-1,1:end-1,1:end-1].+Axyz[2:end,1:end-1,1:end-1].+Axyz[2:end,2:end,1:end-1].+Axyz[2:end,2:end,2:end].+Axyz[1:end-1,2:end,2:end].+Axyz[1:end-1,1:end-1,2:end].+Axyz[2:end,1:end-1,2:end].+Axyz[1:end-1,2:end,1:end-1])*0.125)
                R.=0; @parallel av_xa!(R, Ax);      @test all(R .== (Ax[2:end,    :,    :].+Ax[1:end-1,      :,    :]).*0.5)
                R.=0; @parallel av_ya!(R, Ay);      @test all(R .== (Ay[    :,2:end,    :].+Ay[      :,1:end-1,    :]).*0.5)
                R.=0; @parallel av_za!(R, Az);      @test all(R .== (Az[    :,    :,2:end].+Az[      :,    :,1:end-1]).*0.5)
                R.=0; @parallel av_xi!(R, Axyyzz);  @test all(R .== (Axyyzz[2:end  ,2:end-1,2:end-1].+Axyyzz[1:end-1,2:end-1,2:end-1]).*0.5)
                R.=0; @parallel av_yi!(R, Axxyzz);  @test all(R .== (Axxyzz[2:end-1,2:end  ,2:end-1].+Axxyzz[2:end-1,1:end-1,2:end-1]).*0.5)
                R.=0; @parallel av_zi!(R, Axxyyz);  @test all(R .== (Axxyyz[2:end-1,2:end-1,2:end  ].+Axxyyz[2:end-1,2:end-1,1:end-1]).*0.5)
                R.=0; @parallel av_xya!(R, Axy);    @test all(R .== (Axy[1:end-1,1:end-1,:].+Axy[2:end,1:end-1,:].+Axy[1:end-1,2:end,:].+Axy[2:end,2:end,:])*0.25)
                R.=0; @parallel av_xza!(R, Axz);    @test all(R .== (Axz[1:end-1,:,1:end-1].+Axz[2:end,:,1:end-1].+Axz[1:end-1,:,2:end].+Axz[2:end,:,2:end])*0.25)
                R.=0; @parallel av_yza!(R, Ayz);    @test all(R .== (Ayz[:,1:end-1,1:end-1].+Ayz[:,2:end,1:end-1].+Ayz[:,1:end-1,2:end].+Ayz[:,2:end,2:end])*0.25)
                R.=0; @parallel av_xyi!(R, Axyzz);  @test all(R .== (Axyzz[1:end-1,1:end-1,2:end-1].+Axyzz[2:end,1:end-1,2:end-1].+Axyzz[1:end-1,2:end,2:end-1].+Axyzz[2:end,2:end,2:end-1])*0.25)
                R.=0; @parallel av_xzi!(R, Axyyz);  @test all(R .== (Axyyz[1:end-1,2:end-1,1:end-1].+Axyyz[2:end,2:end-1,1:end-1].+Axyyz[1:end-1,2:end-1,2:end].+Axyyz[2:end,2:end-1,2:end])*0.25)
                R.=0; @parallel av_yzi!(R, Axxyz);  @test all(R .== (Axxyz[2:end-1,1:end-1,1:end-1].+Axxyz[2:end-1,2:end,1:end-1].+Axxyz[2:end-1,1:end-1,2:end].+Axxyz[2:end-1,2:end,2:end])*0.25)
            end;
            @testset "others" begin
                @require @is_initialized()
                @parallel maxloc!(R, Axxyyzz) = (@all(R) = @maxloc(Axxyyzz); return)
                R.=0; @parallel maxloc!(R, Axxyyzz); @test all(R .== max.(max.(max.(max.(max.(max.(Axxyyzz[1:end-2,2:end-1,2:end-1],Axxyyzz[3:end,2:end-1,2:end-1]),Axxyyzz[2:end-1,2:end-1,2:end-1]),Axxyyzz[2:end-1,1:end-2,2:end-1]),Axxyyzz[2:end-1,3:end,2:end-1]),Axxyyzz[2:end-1,2:end-1,1:end-2]),Axxyyzz[2:end-1,2:end-1,3:end]))
            end;
        end;
        @testset "2. apply masks" begin
            @testset "selection" begin
                @require @is_initialized()
                @parallel inn_all!(Rxxyyzz, A)       = (@inn(Rxxyyzz) = @all(A); return)
                @parallel inn_inn!(Rxxyyzz, Axxyyzz) = (@inn(Rxxyyzz) = @inn(Axxyyzz); return)
                Rxxyyzz.=0; @parallel inn_all!(Rxxyyzz, A);        @test all(Rxxyyzz[2:end-1,2:end-1,2:end-1] .== A)
                Rxxyyzz[2:end-1,2:end-1,2:end-1].=0;               @test all(Rxxyyzz .== 0)  # Test that boundary values remained zero.
                Rxxyyzz.=0; @parallel inn_inn!(Rxxyyzz, Axxyyzz);  @test all(Rxxyyzz[2:end-1,2:end-1,2:end-1] .== Axxyyzz[2:end-1,2:end-1,2:end-1])
                Rxxyyzz[2:end-1,2:end-1,2:end-1].=0;               @test all(Rxxyyzz .== 0)  # Test that boundary values remained zero.
            end;
            @testset "differences" begin
                @require @is_initialized()
                @parallel  inn_d_xa!(Rxxyyzz, Ax)       = (@inn(Rxxyyzz) = @d_xa(Ax); return)
                @parallel  inn_d_yi!(Rxxyyzz, Axxyzz)   = (@inn(Rxxyyzz) = @d_yi(Axxyzz); return)
                @parallel  inn_d2_yi!(Rxxyyzz, Axxyyzz) = (@inn(Rxxyyzz) = @d2_yi(Axxyyzz); return)
                Rxxyyzz.=0; @parallel  inn_d_xa!(Rxxyyzz, Ax);       @test all(Rxxyyzz[2:end-1,2:end-1,2:end-1] .== Ax[2:end,    :,    :].-Ax[1:end-1,      :,    :])
                Rxxyyzz[2:end-1,2:end-1,2:end-1].=0;                 @test all(Rxxyyzz .== 0)  # Test that boundary values remained zero.
                Rxxyyzz.=0; @parallel  inn_d_yi!(Rxxyyzz, Axxyzz);   @test all(Rxxyyzz[2:end-1,2:end-1,2:end-1] .== Axxyzz[2:end-1,2:end  ,2:end-1].-Axxyzz[2:end-1,1:end-1,2:end-1])
                Rxxyyzz[2:end-1,2:end-1,2:end-1].=0;                 @test all(Rxxyyzz .== 0)  # Test that boundary values remained zero.
                Rxxyyzz.=0; @parallel inn_d2_yi!(Rxxyyzz, Axxyyzz);  @test all(Rxxyyzz[2:end-1,2:end-1,2:end-1] .== (Axxyyzz[2:end-1,3:end,2:end-1].-Axxyyzz[2:end-1,2:end-1,2:end-1]).-(Axxyyzz[2:end-1,2:end-1,2:end-1].-Axxyyzz[2:end-1,1:end-2,2:end-1]))
                Rxxyyzz[2:end-1,2:end-1,2:end-1].=0;                 @test all(Rxxyyzz .== 0)  # Test that boundary values remained zero.
            end;
        end;
        @reset_parallel_stencil()
    end;
)) end == nothing || true;


text:
"""
    rkf45(f,a,b,y0,M,delta)
berisi program untuk mencari solusi persamaan
differensial `y' = f(t,y)` dengan `y(a) = y0` pada interval `[a, b]`. Program ini secara
default berisi 6 masukan, yaitu fungsi `f(t,y)`,
titik ujung interval penyelesaian `[a,b]`,
nilai awal `y0`, jumlah sub-interval `M` dan nilai toleransi error `delta`

# Example
```jl
julia> f(t,y) = 1+y^2;

julia> sol = rkf45(f,0,1.4,0,7,2e-5)
14×2 Array{Float64,2}:
 0.0     0.0
 0.2     0.20271
 0.6     0.684165
 0.8     1.02968
 ⋮
 1.35    4.45586
 1.375   5.04273
 1.3875  5.39534
 1.4     5.79895
```
return solusi masalah nilai awal `sol`.
"""
function rkf45(f,a,b,y0,M,delta)
  M = round(M)
  a2=1/4;b2=1/4;a3=3/8;b3=3/32;c3=9/32;
  a4=12/13;b4=1932/2197;c4=-7200/2197;d4=7296/2197;
  a5=1;b5=439/216;c5=-8;d5=3680/513;e5=-845/4104;
  a6=1/2;b6=-8/27;c6=2;d6=-3544/2565;e6=1859/4104;f6=-11/40;
  r1=1/360;r3=-128/4275;r4=-2197/75240;r5=1/50;r6=2/55;
  n1=25/216;n3=1408/2565;n4=2197/4104;n5=-1/5;
  big=1e15;
  h=(b-a)/M;
  hmin=h/64;
  hmax=h*64;
  maxi=200;
  j=1;
  Y = y0;
  T = a;
  br= b-0.001*abs(b);
  err=NaN
  while T[j]<b
    if (T[j]+h)>br;h=b-T[j];end
    #% Hitung koefisien
    k1=h*f(T[j],Y[j]);
    y2=Y[j]+b2*k1;
    k2=h*f(T[j]+a2*h,y2);
    y3=Y[j]+b3*k1+c3*k2;
    k3=h*f(T[j]+a3*h,y3);
    y4=Y[j]+b4*k1+c4*k2+d4*k3;
    k4=h*f(T[j]+a4*h,y4);
    y5=Y[j]+b5*k1+c5*k2+d5*k3+e5*k4;
    k5=h*f(T[j]+a5*h,y5);
    y6=Y[j]+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5;
    k6=h*f(T[j]+a6*h,y6);
    err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
    ynew=Y[j]+n1*k1+n3*k3+n4*k4+n5*k5;
    #% Perbarui ukuran langkah
    if (err<delta) || (h<2*hmin)
      Y = [Y; ynew];
      if (T[j]+h)>br
        T = [T; b];
      else
        T = [T; T[j]+h];
      end
      j=j+1;
    end
    if (err==0)
      s=0;
    else
      s=0.84*(delta*h/err)^(0.25);
    end
    if (s<0.75)&&(h>2*hmin);h = h/2;end
    if (s>1.50)&&(2*h<hmax);h = 2*h;end
    if abs(Y[j])>big || maxi==j;break;end
  end
  sol = [T Y];
  return sol
end


text:
using Compat
using Compat.REPL
using Compat.REPL.Terminals
using TerminalMenus
using ObjectFile.ELF
using HTTP
import PkgLicenses
using JSON
using MbedTLS

# It's Magic!
export run_wizard

include("wizard/state.jl")
include("wizard/utils.jl")
include("wizard/obtain_source.jl")
include("wizard/interactive_build.jl")
include("wizard/deploy.jl")

# This is here so that if the wizard crashes, we may have a shot at resuming.
last_wizard_state = WizardState()

function run_wizard(state::WizardState = WizardState())
    global last_wizard_state

    if last_wizard_state.step != :done && last_wizard_state.step != :step1
        terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
        choice = request(terminal,
            "Would you like to resume the previous incomplete wizard run?",
            RadioMenu([
                "Resume previous run",
                "Start from scratch",
            ]),
        )

        if choice == 1
            state = last_wizard_state
        end
    end

    if state.step == :step1
        print_wizard_logo(state.outs)

        println(state.outs,
            "Welcome to the BinaryBuilder wizard.\n"*
            "We'll get you set up in no time.\n")
    end

    try
        while state.step != :done
            if state.step == :step1
                step1(state)
                state.step = :step2
            elseif state.step == :step2
                step2(state)
                state.step = :step3
            elseif state.step == :step3
                step34(state)
            elseif state.step == :step3_retry
                step3_retry(state)
            elseif state.step == :step5a
                step5a(state)
            elseif state.step == :step5b
                step5b(state)
            elseif state.step == :step5c
                step5c(state)
                state.step = :step6
            elseif state.step == :step6
                step6(state)
            elseif state.step == :step7
                step7(state)
                state.step = :done
            end
        end
    catch err
        last_wizard_state = state
        if isa(err, InterruptException)
            msg = "\n\nWizard stopped, use run_wizard() to resume.\n\n"
            printstyled(state.outs, msg, bold=true, color=:red)
        else
            bt = catch_backtrace()
            Base.showerror(stderr, err, bt)
            println(state.outs, "\n")
        end
        return state
    end

    println(state.outs, "\nWizard Complete. Press any key to exit...")
    read(state.ins, Char)

    state
end


text:
"""
    (data, lon, lat, time) = read_netcdf2d(file, varname="")

Read variable `varname` (optional) from netcdf `file`.
Returns data and coordinate and time axes.
"""
function  read_netcdf2d(file, varname="")

nc = NetCDF.open(file, mode=NC_NOWRITE)
if( haskey(nc.vars,"lon") )
   lat=ncread(file,"lat");
   lon=ncread(file,"lon");
elseif (haskey(nc.vars,"longitude"))
   lat=ncread(file,"latitude");
   lon=ncread(file,"longitude");
elseif (haskey(nc.vars,"x"))
   lat=ncread(file,"y");
   lon=ncread(file,"x");
elseif (haskey(nc.vars,"X"))
   lat=ncread(file,"Y");
   lon=ncread(file,"X");
else
   println("Input file does not contain lon or longitude dimensional variables")
   quit(1)
end

if(varname=="")
   for key in  keys(nc.vars)
      if(key!="lon" && key!="lat" && key!="time" && key!="longitude" && key!="latitude" && key!="time_bnds" && key!="x" && key!="y" && key!="X" && key!="Y" )
            varname=key
      end
   end
end

var=ncread(file,varname);

var=float(var);
if(haskey(nc.vars[varname].atts,"missing_value"))
    xmiss=nc.vars[varname].atts["missing_value"]
    xmiss=convert(typeof(var[1,1,1]),xmiss)
    var[var .== xmiss] .= NaN;
end
NetCDF.close(nc)

return var,lon,lat,varname
end

#function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
#    m, n = length(vy), length(vx)
#    vx = reshape(vx, 1, n)
#    vy = reshape(vy, m, 1)
#    (repmat(vx, m, 1), repmat(vy, 1, n))
#end


text:
const tfn = pso.TestFunctions

@testset "TestFunctions" begin

    @testset "sphere function" begin
        @test tfn.sphere([0]) == 0
        @test tfn.sphere([0.0]) == 0
        @test tfn.sphere(fill(0, 10)) == 0
        @test tfn.sphere([3.0, 4.0]) == 25
    end

    @testset "euclidean_distance function" begin
        @test tfn.euclidean_distance([0]) == 0
        @test tfn.euclidean_distance([0.0]) == 0
        @test tfn.euclidean_distance(fill(0, 10)) == 0
        @test tfn.euclidean_distance([3.0, 4.0]) == 5
    end

    @testset "rastrigin function" begin
        @test tfn.rastrigin([0]) == 0
        @test tfn.rastrigin([0.0]) == 0
        @test tfn.rastrigin(fill(0, 10)) == 0
        @test tfn.rastrigin([0.0, 0.0], A=-2.3) == 0.0
    end

    @testset "ackley's function" begin
        @test tfn.ackley([0,0]) == 0
        @test tfn.ackley([0.0, 0.0]) == 0
    end

    @testset "rosenbrock's function" begin
        @test tfn.rosenbrock([1.0]) == 0
        @test tfn.rosenbrock([1, 1]) == 0
        @test tfn.rosenbrock(fill(1, 10)) == 0.0
    end

    @testset "eggholder function" begin
        @test tfn.eggholder(512, 404.2319) ≈ -959.6407 atol=0.0001
        @test tfn.eggholder([512, 404.2319]) ≈ -959.6407 atol=0.0001
    end

    @testset "beale function" begin
        @test tfn.beale(3, 0.5) == 0
        @test tfn.beale([3, 0.5]) == 0
    end

    @testset "square_movingpeak function" begin
        @test tfn.square_movingpeak([0, 0]) == 0
        @test tfn.square_movingpeak([3, 4]) == 25
        @test tfn.square_movingpeak([1, 1], optimum=[1, 1]) == 0
        @test tfn.square_movingpeak([2, 2], optimum=[1, 1]) ≈ 2
    end

    @testset "ellipse function" begin
        @test tfn.ellipse(-0.5, -1.5) == 0
        @test tfn.ellipse([-0.5, -1.5]) == 0
    end

    @testset "stybliski function" begin
        for i in 1:5
            @test tfn.stybliski(fill(-2.903534, i)) ≈ -39.16616*i atol=0.0001
        end
    end

    @testset "mccormick function" begin
        @test tfn.mccormick(0, 0) == 1
        @test tfn.mccormick([0, 0]) == 1
        @test tfn.mccormick(-0.54719, -1.54719) ≈ -1.9133 atol=0.0001
    end

    @testset "himmelblau function" begin
        @test tfn.himmelblau(-0.270845, -0.923039) ≈ 181.617 atol=0.001
        @test tfn.himmelblau(3, 2) ≈ 0
        @test tfn.himmelblau([3, 2]) ≈ 0
        @test tfn.himmelblau(-2.805118,  3.131312) ≈ 0 atol=0.0001
        @test tfn.himmelblau(-3.77931 , -3.283186) ≈ 0 atol=0.0001
        @test tfn.himmelblau( 3.584428, -1.848126) ≈ 0 atol=0.0001
    end

    @testset "griewank function" begin
        @test tfn.griewank([0, 0]) == 0
    end

    @testset "schafferf6 function" begin
        @test tfn.schafferf6(0, 0) == 0
        @test tfn.schafferf6([0, 0]) == 0
    end

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Base.LinAlg.LAPACK.chklapackerror, (Int32,))
    precompile(Base.string, (ASCIIString, Char,))
    precompile(Base.Broadcast.broadcast_shape, (Array{Float64, 1},))
    precompile(Base.throw_boundserror, (Array{Float64, 2}, Tuple{Base.UnitRange{Int32}, Base.UnitRange{Int32}},))
    precompile(Base.Serializer.serialize_type, (Base.SerializationState{Base.Pipe}, DataType,))
    precompile(Base.Broadcast.broadcast_shape, (Array{Float64, 2},))
    precompile(Base.Broadcast.check_broadcast_shape, (Tuple{Int32}, Array{Float64, 1},))
    precompile(Base.Broadcast.check_broadcast_shape, (Tuple{Int32, Int32}, Array{Float64, 2},))
    precompile(Base.LinAlg._generic_matmatmul!, (Array{Float64, 2}, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.throw_checksize_error, (Array{Float64, 2}, Int32, Base.UnitRange{Int32},))
    precompile(Base.Filesystem.pathsep, (ASCIIString,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, UInt8,))
    precompile(Base.Serializer.serialize_array_data, (Base.Pipe, Array{UInt8, 1},))
    precompile(Base.copy!, (Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32}, Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32},))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, Array{UInt8, 1},))
    precompile(Base.write, (Base.Pipe, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base.resize!, (Array{Base.Dict{Int32, Any}, 1}, Int32,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, ASCIIString,))
    precompile(Base.Filesystem.lstat, (ASCIIString,))
    precompile(Base.check_open, (Base.PipeEndpoint,))
    precompile(Base.stride, (Array{Float64, 1}, Int32,))
    precompile(Base.rehash!, (Base.Dict{Int32, Any}, Int32,))
    precompile(Base.open, (Base.CmdRedirect, ASCIIString, Base.TTY,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, QuoteNode,))
    precompile(Base.LinAlg.istriu, (Array{Float64, 2},))
    precompile(Base.Serializer.write_as_tag, (Base.Pipe, Int32,))
    precompile(Base.malloc_julia_pipe!, (Base.PipeEndpoint,))
    precompile(Base.Serializer.serialize_mod_names, (Base.SerializationState{Base.Pipe}, Module,))
    precompile(Base.Broadcast.gen_broadcast_body_iter, (Int32, Int32, Base.MulFun,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base._link_pipe, (Ptr{Void}, Ptr{Void},))
    precompile(Base.copy, (Array{Float64, 2},))
    precompile(Base.write, (Base.Pipe, Array{UInt8, 1},))
    precompile(Base.unsafe_write, (Base.Pipe, Base.RefValue{Int32}, Int32,))
    precompile(Base.resize!, (Array{Base.Dict{Int32, Base.Dict{Int32, Any}}, 1}, Int32,))
    precompile(Base.fill!, (Array{Float64, 1}, Int32,))
    precompile(Base.copy_transpose!, (Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32}, Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32},))
    precompile(Base.unsafe_copy!, (Array{Float64, 2}, Int32, Array{Float64, 2}, Int32, Int32,))
    precompile(Base.LinAlg.generic_matvecmul!, (Array{Float64, 1}, Char, Array{Float64, 2}, Array{Float64, 1},))
    precompile(Base._jl_spawn, (UTF8String, Array{Union{UTF8String, ASCIIString}, 1}, Ptr{Void}, Base.Process, Base.PipeEndpoint, Base.TTY, Base.TTY,))
    precompile(Base.unsafe_convert, (Type{Ptr{Ptr{UInt8}}}, Base.RefArray{Ptr{UInt8}, Array{Ptr{UInt8}, 1}, Any},))
    precompile(Base.LinAlg.triu!, (Array{Float64, 2}, Int32,))
    precompile(Base.LinAlg.BLAS.gemv!, (Char, Float64, Array{Float64, 2}, Array{Float64, 1}, Float64, Array{Float64, 1},))
    precompile(Base.Broadcast.gen_broadcast_body_iter, (Int32, Int32, Function,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Tuple{Int32},))
    precompile(Base.zero, (Type{Float64},))
    precompile(Base.transposeblock!, (Array{Float64, 2}, Array{Float64, 2}, Int32, Int32, Int32, Int32,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base.julia_cmd, (UTF8String,))
    precompile(Base.getindex, (Tuple{Symbol}, Base.UnitRange{Int32},))
    precompile(Base.call, (Array{Any, 1}, Type{Base.Cmd}, Base.Cmd,))
    precompile(Base.LinAlg.generic_lufact!, (Array{Float64, 2}, Type{Base.Val{true}},))
    precompile(Base.LinAlg.LAPACK.geevx!, (Char, Char, Char, Char, Array{Float64, 2},))
    precompile(Base.LinAlg.LAPACK.getri!, (Array{Float64, 2}, Array{Int32, 1},))
    precompile(Base.LinAlg.istril, (Array{Float64, 2},))
    precompile(Base.LinAlg.matmul3x3!, (Array{Float64, 2}, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.uv_write, (Base.PipeEndpoint, Ptr{UInt8}, UInt32,))
    precompile(Base.Broadcast.gen_broadcast_function, (Function, Int32, Int32, Base.MulFun,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, LineNumberNode,))
    precompile(Base.LinAlg.LAPACK.chkfinite, (Array{Float64, 2},))
    precompile(Base.Serializer.serialize_any, (Base.SerializationState{Base.Pipe}, LineNumberNode,))
    precompile(Base.Broadcast.gen_broadcast_function, (Function, Int32, Int32, Function,))
    precompile(Base.LinAlg.lapack_size, (Char, Array{Float64, 2},))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Symbol,))
    precompile(Base.call, (Type{Base.Process}, Base.Cmd, Ptr{Void}, Base.Pipe, Base.TTY, Base.TTY,))
    precompile(Base.LinAlg.matmul2x2!, (Array{Float64, 2}, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Array{UInt8, 1},))
    precompile(Base.sizeof, (Array{Float64, 2},))
    precompile(Base.Broadcast.bitcache_eq, (Array{Float64, 1}, Int32, Int32, Int32, Array{Bool, 1},))
    precompile(Base.ht_keyindex2, (Base.Dict{Any, Base.Dict{Int32, Base.Dict{Int32, Any}}}, Base.MulFun,))
    precompile(Base.create_expr_cache, (ASCIIString, ASCIIString,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Int32,))
    precompile(Base.julia_cmd, ())
    precompile(Base.LinAlg.BLAS.gemm!, (Char, Char, Float64, Array{Float64, 2}, Array{Float64, 2}, Float64, Array{Float64, 2},))
    precompile(Base.Cartesian._ncall, (Int32, Expr, Symbol,))
    precompile(Base.mapreduce_impl, (Base.IdFun, Base.MinFun, Array{Float64, 1}, Int32, Int32,))
    precompile(Base.rehash!, (Base.Dict{Int32, Base.Dict{Int32, Any}}, Int32,))
    precompile(Base.promote_type, (Type{UInt8}, Type{Union{}},))
    precompile(Base.setindex!, (Array{Float64, 2}, Array{Float64, 1}, Base.StepRange{Int32, Int32},))
    precompile(Base.call, (Type{Base.RefArray{Union{UTF8String, ASCIIString}, Array{Union{UTF8String, ASCIIString}, 1}, Any}}, Array{Union{UTF8String, ASCIIString}, 1}, Int32, Bool,))
    precompile(Base.uvfinalize, (Base.Process,))
    precompile(Base.throw_boundserror, (Array{Float64, 1}, Base.UnitRange{Int32},))
    precompile(Base.LinAlg.char_uplo, (Symbol,))
    precompile(Base.cmd_gen, (Tuple{Tuple{Base.Cmd}, Tuple{ASCIIString}, Tuple{ASCIIString}, Tuple{ASCIIString}, Tuple{ASCIIString}, Tuple{ASCIIString}, Tuple{ASCIIString, ASCIIString}, Tuple{ASCIIString}, Tuple{ASCIIString}},))
    precompile(Base.convert, (Type{Base.IO}, Base.TTY,))
    precompile(Base.setup_stdio, (Function, Tuple{Base.Pipe, Base.TTY, Base.TTY},))
    precompile(Base.throw_setindex_mismatch, (Array{Float64, 1}, Tuple{Int32},))
    precompile(Base.fill!, (Array{Float64, 1}, Float64,))
    precompile(Base.spawn, (Array{Any, 1}, Base.Cmd, Tuple{Base.Pipe, Base.TTY, Base.TTY}, Bool, Bool,))
    precompile(Base.rehash!, (Base.Dict{Any, Base.Dict{Int32, Base.Dict{Int32, Any}}}, Int32,))
    precompile(Base.unsafe_copy!, (Array{Float64, 1}, Int32, Array{Float64, 1}, Int32, Int32,))
    precompile(Base.LinAlg.LAPACK.trtrs!, (Char, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.complex, (Array{Float64, 1}, Array{Float64, 1},))
    precompile(Base.Filesystem.unlink, (ASCIIString,))
    precompile(Base.LinAlg.copy!, (Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32}, Char, Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32},))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Module,))
    precompile(Base.LinAlg.inv, (Base.LinAlg.LowerTriangular{Float64, Array{Float64, 2}},))
    precompile(Base.abs2, (Array{Float64, 1},))
    precompile(Base._setindex!, (Base.Dict{Any, Base.Dict{Int32, Base.Dict{Int32, Any}}}, Base.Dict{Int32, Base.Dict{Int32, Any}}, Function, Int32,))
    precompile(Base._setindex!, (Base.Dict{Any, Base.Dict{Int32, Base.Dict{Int32, Any}}}, Base.Dict{Int32, Base.Dict{Int32, Any}}, Base.MulFun, Int32,))
    precompile(Base.promote_typeof, (Float64, Int32,))
    precompile(Base.with_output_color, (Function, Symbol, Base.TTY, ASCIIString,))
    precompile(Base.sizeof, (Array{Float64, 1},))
    precompile(Base._setindex!, (Base.Dict{Int32, Base.Dict{Int32, Any}}, Base.Dict{Int32, Any}, Int32, Int32,))
    precompile(Base.LinAlg.issym, (Array{Float64, 2},))
    precompile(Base.LinAlg.chol!, (Array{Float64, 2}, Type{Base.LinAlg.LowerTriangular},))
    precompile(Base.call, (Type{Base.Dict{Int32, Base.Dict{Int32, Any}}},))
    precompile(Base.compilecache, (ASCIIString,))
    precompile(Base.pipeline, (Array{Any, 1}, Base.Cmd,))
    precompile(Base.LinAlg.LAPACK.getrf!, (Array{Float64, 2},))
    precompile(Base.getindex, (Tuple{Expr, Expr}, Base.UnitRange{Int32},))
    precompile(Base.LinAlg.chol!, (Array{Float64, 2}, Type{Base.LinAlg.UpperTriangular},))
    precompile(Base.Cartesian._ncall, (Int32, Symbol, Expr,))
    precompile(Base.getindex, (Array{Float64, 1}, Int32,))
    precompile(Base.Cartesian._nref, (Int32, Symbol, Symbol,))
    precompile(Base.LinAlg.LAPACK.chktrans, (Char,))
    precompile(Base.vcat, (Array{Float64, 1}, Array{Float64, 1},))
    precompile(Base.all, (Base.BitArray{1},))
    precompile(Base.arg_gen, (ASCIIString, ASCIIString,))
    precompile(Base.transpose, (Array{Float64, 1},))
    precompile(Base.call, (Type{Base.RefArray{Ptr{UInt8}, Array{Ptr{UInt8}, 1}, Any}}, Array{Ptr{UInt8}, 1}, Int32, Array{Any, 1},))
    precompile(Base._unsafe_getindex, (Base.LinearFast, Array{Float64, 2}, Base.Colon, Base.UnitRange{Int32},))
    precompile(Base.repl_color, (ASCIIString, Symbol,))
    precompile(Base.hvcat_fill, (Array{Float64, 1}, Tuple{Float64, Float64, Float64, Int32, Int32, Int32},))
    precompile(Base.fill!, (Array{Float64, 2}, Float64,))
    precompile(Base.flush, (Base.PipeEndpoint,))
    precompile(Base.LinAlg.diagind, (Int32, Int32, Int32,))
    precompile(Base.cmd_gen, (Tuple{Tuple{UTF8String}, Tuple{ASCIIString, ASCIIString}, Tuple{ASCIIString, ASCIIString}},))
    precompile(Base.fill!, (Array{Float64, 2}, Int32,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, SimpleVector,))
    precompile(Base.Broadcast.longer_size, (Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.call, (Type{Base.Dict{Int32, Any}},))
    precompile(Base._setindex!, (Base.Dict{Int32, Any}, Function, Int32, Int32,))
    precompile(Base.push!, (Array{Float64, 1}, Float64,))
    precompile(Base.promote_shape, (Tuple{Int32, Int32}, Tuple{Int32, Int32},))
    precompile(Base.getindex, (Type{Union{UTF8String, ASCIIString}}, UTF8String,))
    precompile(Base.call, (Type{Array}, Type{Float64}, Tuple{Int32, Int32},))
    precompile(Base.throw_boundserror, (Array{Float64, 2}, Tuple{Base.Colon, Base.UnitRange{Int32}},))
    precompile(Base.Serializer.serialize_any, (Base.SerializationState{Base.Pipe}, QuoteNode,))
    precompile(Base.info, (Array{Any, 1}, Base.TTY, ASCIIString,))
    precompile(Base.vcat, (Float64, Float64, Float64,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Bool,))
    precompile(Base.Serializer.serialize_any, (Base.SerializationState{Base.Pipe}, UInt8,))
    precompile(Base.unsafe_convert, (Type{Ptr{Ptr{UInt8}}}, Ptr{Void},))
    precompile(Base.LinAlg.generic_matmatmul!, (Array{Float64, 2}, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.length, (Array{Float64, 1},))
    precompile(Base._uv_hook_close, (Base.Process,))
    precompile(Base.Broadcast.broadcast!, (Base.MulFun, Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.unsafe_convert, (Type{Ptr{UInt8}}, Ptr{Void},))
    precompile(Base.setindex!, (Array{Float64, 1}, Float64, Int32, Int32,))
    precompile(Base.promote_typeof, (Float64, Float64,))
    precompile(Base.LinAlg.LAPACK.chkdiag, (Char,))
    precompile(Base.Broadcast.broadcast!, (Function, Array{Float64, 1}, Array{Float64, 1}, Array{Float64, 1},))
    precompile(Base.append!, (Array{Union{UTF8String, ASCIIString}, 1}, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base.call, (Type{Array}, Type{Float64}, Int32, Int32,))
    precompile(Base.unsafe_convert, (Type{Ptr{UInt8}}, ASCIIString,))
    precompile(Base.LinAlg.LAPACK.syevr!, (Char, Char, Char, Array{Float64, 2}, Float64, Float64, Int32, Int32, Float64,))
    precompile(Base.containsnul, (UTF8String,))
    precompile(Base.Serializer.serialize_any, (Base.SerializationState{Base.Pipe}, ASCIIString,))
    precompile(Base.println_with_color, (Symbol, Base.TTY, ASCIIString,))
    precompile(Base.Cartesian._nloops, (Int32, Symbol, Symbol, Expr,))
    precompile(Base.LinAlg._cholfact!, (Array{Float64, 2}, Type{Base.Val{false}}, Symbol,))
    precompile(Base.eye, (Type{Float64}, Int32, Int32,))
    precompile(Base.Broadcast.longer_size, (Array{Float64, 1}, Array{Float64, 1},))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, Expr,))
    precompile(Base.uvfinalize, (Base.PipeEndpoint,))
    precompile(Base.Broadcast.longer_size, (Array{Float64, 2},))
    precompile(Base.LinAlg.gemm_wrapper!, (Array{Float64, 2}, Char, Char, Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.convert, (Type{Float64}, Int32,))
    precompile(Base.Serializer.serialize_cycle, (Base.SerializationState{Base.Pipe}, Expr,))
    precompile(Base.getindex, (Array{Union{UTF8String, ASCIIString}, 1}, Int32,))
    precompile(Base.ht_keyindex2, (Base.Dict{Any, Base.Dict{Int32, Base.Dict{Int32, Any}}}, Function,))
    precompile(Base.Filesystem.isdir, (ASCIIString,))
    precompile(Base.close, (Base.Pipe,))
    precompile(Base.unsafe_write, (Base.PipeEndpoint, Ptr{UInt8}, UInt32,))
    precompile(Base.setindex!, (Array{Float64, 1}, Int32, Int32,))
    precompile(Base.bytestring, (ASCIIString, ASCIIString,))
    precompile(Base.promote_op, (Base.AddFun, Type{Float64}, Type{Float64},))
    precompile(Base.cstr, (ASCIIString,))
    precompile(Base.Filesystem.rm, (Array{Any, 1}, ASCIIString,))
    precompile(Base.promote_op, (Base.SubFun, Type{Float64}, Type{Float64},))
    precompile(Base.spawn_opts_swallow, (Tuple{Base.Pipe, Base.TTY, Base.TTY},))
    precompile(Base.LinAlg.lufact!, (Array{Float64, 2}, Type{Base.Val{true}},))
    precompile(Base.transpose!, (Array{Float64, 2}, Array{Float64, 2},))
    precompile(Base.LinAlg.inv, (Base.LinAlg.LU{Float64, Array{Float64, 2}},))
    precompile(Base.promote_typeof, (Int32,))
    precompile(Base.cconvert, (Type{Ptr{Ptr{UInt8}}}, Ptr{Void},))
    precompile(Base.LinAlg.copy_transpose!, (Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32}, Char, Array{Float64, 2}, Base.UnitRange{Int32}, Base.UnitRange{Int32},))
    precompile(Base.copy, (Array{Float64, 1},))
    precompile(Base.promote_typeof, (Int32, Int32,))
    precompile(Base.Filesystem.joinpath, (ASCIIString, UTF8String,))
    precompile(Base.success, (Base.Process,))
    precompile(Base.setindex!, (Array{Float64, 1}, Int32, Int32, Int32,))
    precompile(Base.unsafe_convert, (Type{Ptr{UInt8}}, UTF8String,))
    precompile(Base.cconvert, (Type{Ptr{UInt8}}, UTF8String,))
    precompile(Base.LinAlg.inv, (Array{Float64, 2},))
    precompile(Base.indexed_next, (Tuple{Base.Pipe, Base.Process}, Int32, Int32,))
    precompile(Base.start, (Tuple{Base.Pipe, Base.Process},))
    precompile(Base.convert, (Type{Ptr{Void}}, Ptr{Void},))
    precompile(Base.link_pipe, (Base.PipeEndpoint, Bool, Base.PipeEndpoint, Bool,))
    precompile(Base.LinAlg.inv, (Base.LinAlg.UpperTriangular{Float64, Array{Float64, 2}},))
    precompile(Base.cconvert, (Type{Ptr{UInt8}}, Ptr{Void},))
    precompile(Base.setup_stdio, (Base.Pipe, Bool,))
    precompile(Base.promote_result, (Type{Bool}, Type{UInt8}, Type{UInt8}, Type{Union{}},))
    precompile(Base.LinAlg.parent, (Base.LinAlg.UpperTriangular{Float64, Array{Float64, 2}},))
    precompile(Base.Broadcast.longer_size, (Array{Float64, 1},))
    precompile(Base.ht_keyindex2, (Base.Dict{Int32, Base.Dict{Int32, Any}}, Int32,))
    precompile(Base.getindex, (Array{Function, 1}, Int32,))
    precompile(Base.call, (Type{Ref{Ptr{UInt8}}}, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base.eltype, (Array{Float64, 1},))
    precompile(Base.ht_keyindex2, (Base.Dict{Int32, Any}, Int32,))
    precompile(Base.Cartesian.lreplace!, (LineNumberNode, Base.Cartesian.LReplace{ASCIIString},))
    precompile(Base.promote_type, (Type{Bool}, Type{UInt8},))
    precompile(Base.promote_type, (Type{Float64}, Type{Float64},))
    precompile(Base.arg_gen, (Base.Cmd,))
    precompile(Base.LinAlg.eigvals!, (Array{Any, 1}, Array{Float64, 2},))
    precompile(Base.Serializer.serialize, (Base.Pipe, Expr,))
    precompile(Base.Serializer.serialize_type_data, (Base.SerializationState{Base.Pipe}, DataType,))
    precompile(Base.astcopy, (LineNumberNode,))
    precompile(Base.Serializer.serialize_array_data, (Base.Pipe, Array{Union{UTF8String, ASCIIString}, 1},))
    precompile(Base.LinAlg.gemv!, (Array{Float64, 1}, Char, Array{Float64, 2}, Array{Float64, 1},))
    precompile(Base.redir_err, (Base.Cmd, Base.TTY,))
    precompile(Base.ctranspose, (Array{Float64, 2},))
    precompile(Base.LinAlg.convert, (Type{Array{Float64, 2}}, Base.LinAlg.UpperTriangular{Float64, Array{Float64, 2}},))
    precompile(Base.promote_rule, (Type{UInt8}, Type{Bool},))
    precompile(Base._mapreduce, (Base.IdFun, Base.MinFun, Base.LinearFast, Array{Float64, 1},))
    precompile(Base.ctranspose, (Array{Float64, 1},))
    precompile(Base.LinAlg.eigvals, (Array{Any, 1}, Array{Float64, 2},))
    precompile(Base.unsafe_write, (Base.PipeEndpoint, Base.RefValue{UInt8}, Int32,))
    precompile(Base.convert, (Type{Base.IO}, Base.Pipe,))
    precompile(Base.minimum, (Array{Float64, 1},))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, DataType,))
    precompile(Base.unsafe_getindex, (Array{Float64, 1}, Base.UnitRange{Int32},))
    precompile(Base.init_pipe!, (Array{Any, 1}, Base.PipeEndpoint,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, ASCIIString,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, QuoteNode,))
    precompile(Base.cconvert, (Type{Ptr{UInt8}}, ASCIIString,))
    precompile(Base.spawn, (Array{Any, 1}, Base.CmdRedirect, Tuple{Base.Pipe, Base.TTY, Base.TTY}, Bool, Bool,))
    precompile(Base.arg_gen, (UTF8String,))
    precompile(Base.copy!, (Array{Float64, 2}, Int32, Array{Float64, 2}, Int32, Int32,))
    precompile(Base.Serializer.serialize, (Base.SerializationState{Base.Pipe}, LineNumberNode,))
    precompile(Base.getindex, (Array{Float64, 1}, Base.UnitRange{Int32},))
end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey


text:
include("main_tests.jl")
include("hits_tests.jl")
include("tools_tests.jl")
include("io_tests.jl")
include("math_tests.jl")
include("mc_tests.jl")

println("All tests complete")


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(particle)
q1[3] = 0.25
v1 = zeros(particle.nq)
v1[1] = 10.0

# ## Time
T = 100
h = 0.01

# ## Simulator
s = Simulator(particle, T, h=h)

# ## Simulate
simulate!(s, q1, v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s,
    fixed_camera=false)


text:
using LatBo.Indices.neighbor_index
facts("Checking caching") do
    lattice = LatBo.LB.D2Q9
    context("Input throws if incorrect") do
        @fact_throws ErrorException LatBo.Indices.cached_cartesian([5, 4, 2], lattice)
        @fact_throws ErrorException LatBo.Indices.cached_periodic([5, 4, 2], lattice)
    end

    context("Cached Cartesian") do
        cartesian = LatBo.Indices.Cartesian([5, 4])
        cached = LatBo.Indices.Cached(cartesian, lattice)

        @fact size(cached.cache) --> (length(lattice), length(cartesian))
        @fact typeof(cached.indexing) --> typeof(cartesian)
        @fact size(cached) --> size(cartesian)
        @fact length(cached) --> length(cartesian)
        for i = 1:length(cartesian), d = 1:length(lattice)
            const coords = gridcoords(cartesian, i)
            if all(coords .> 1) && all(coords .< size(cartesian))
                const expected = neighbor_index(cartesian, i, lattice, d)
                @fact neighbor_index(cached, i, lattice, d) --> expected
            else
                @fact neighbor_index(cached, i, lattice, d) --> 0
            end
        end
    end

    context("Cached Periodic") do
        periodic = LatBo.Indices.Periodic([5, 4])
        cached = LatBo.Indices.Cached(periodic, lattice)
        @fact size(cached.cache) --> (length(lattice), length(periodic))
        @fact typeof(cached.indexing) --> typeof(periodic)
        @fact size(cached) --> size(periodic)
        @fact length(cached) --> length(periodic)
        for i = 1:length(periodic), d = 1:length(lattice)
            @fact neighbor_index(cached, i, lattice, d) --> neighbor_index(periodic, i, lattice, d)
        end
    end
end


text:
using Empirikos
using Test
using InteractiveUtils
using FiniteDifferences
using ForwardDiff
using Hypatia
using QuadGK
using Documenter

# Doctests

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

# Other tests

@testset "Test set default" begin
    include("test_set_defaults.jl")
end

@testset "Test ci tools" begin
    include("test_ci_tools.jl")
end

@testset "Dict function" begin
    include("test_dict_function.jl")
end

@testset "EBayes Intervals" begin
    include("test_intervals.jl")
end

@testset "Sample ordering" begin
    include("test_ordering.jl")
end

@testset "Binomial samples" begin
    include("test_binomial.jl")
end

@testset "Normal samples" begin
    include("test_normal.jl")
end

@testset "Test targets" begin
    include("test_targets.jl")
end

@testset "Chi-squared F localization" begin
    include("test_lord_cressie.jl")
end

@testset "Further F localization interval tests" begin
    include("test_flocalization_intervals.jl")
end

@testset "Kernel density tests" begin
    include("test_kernel_density.jl")
end

@testset "Bernoulli samples" begin
    include("test_bernoulli.jl")
end

@testset "Two kinds of modulus problems" begin
    include("test_modulus.jl")
end

@testset "Compound and DKW F localization" begin
    include("test_compound.jl")
end


text:
using AbstractPlotting
import AbstractPlotting: plot, plot!

const DEFAULT_PYPLOT_OPTS = (color=:mediumblue, linewidth=2)
function plot!(s::Scene, c::Cell; kwargs...)
    cam3d!(s)
    for (i,poly) in enumerate(c)
        lines!(s, push!(getindex.(poly, 1), poly[1][1]),
                  push!(getindex.(poly, 2), poly[1][2]),
                  push!(getindex.(poly, 3), poly[1][3]);
                  DEFAULT_PYPLOT_OPTS..., # default options
                  kwargs...               # possible keyword overrides
              )
    end
    return s
end
plot(c::Cell; kwargs...) = plot!(Scene(), c; kwargs...)


text:
using Dispersal, Test

@testset "allele frequency on LC50" begin
    init =  [0.1 0.2 0.3;
             0.4 0.5 0.6;
             0.7 0.8 0.9]

    output = ArrayOutput(init; tspan=1:3)
    rule = Ruleset(DeltaAlleleFrequencySurv(exposure = 10, LC50=10, hillcoefficient=2.0, deviationPhenotype = 10.0, dominanceDegree = -1))
    sim!(output, rule)
    @test output[1] == init
    @test output[2] ≈ [0.1828  0.3504  0.5016;
                       0.6352  0.75    0.8448;
                       0.9184  0.9696  0.9972] atol=1e-4
    @test output[3] ≈ [0.322707  0.571209  0.751677;
                       0.873187  0.946875  0.984955;
                       0.999613  1.00184   1.00027] atol=1e-4

    #! <1.0 is not always satisfied. Is it because of @fastmath ??
    @test (output[1] .<= 1.01) == trues(3,3)
    @test (0.0 .<= output[1]) == trues(3,3)
    @test (output[2] .<= 1.01) == trues(3,3)
    @test (0.0 .<= output[2]) == trues(3,3)
    @test (output[3] .<= 1.01) == trues(3,3)
    @test (0.0 .<= output[3]) == trues(3,3)
end

@testset "allele frequency on LC50 Map" begin
    pFreqinit =  [0.1 0.2 0.3;
                  0.4 0.5 0.6;
                  0.7 0.8 0.9]

    exposure = [50.0 50.0 50.0;
                20.0 20.0 20.0;
                0.0 0.0 0.0]

    output = ArrayOutput(pFreqinit; tspan=1:5, aux=(exposure=exposure,))
    rule = Ruleset(DeltaAlleleFrequencySurvMap(layerkey=Val(:exposure),
                    LC50=10, hillcoefficient=2.0, deviationPhenotype = 10.0, dominanceDegree = -1.0))
    sim!(output, rule)

    @test output[1] == pFreqinit
    @test output[2] ≈ [ 0.259231  0.489231  0.687692;
                        0.77632   0.9       0.99168;
                        0.7       0.8       0.9] atol=1e-4
    @test output[3] ≈ [ 0.610737  0.968742  1.11622;
                        1.06951   1.05552   1.00618;
                        0.7       0.8       0.9] atol=1e-4

    #! <1.0 is not always satisfied. Is it because of @fastmath ??
    @test (output[1] .<= 1.1) == trues(3,3)
    @test (0.0 .<= output[1]) == trues(3,3)
    @test (output[2] .<= 1.1) == trues(3,3)
    @test (0.0 .<= output[2]) == trues(3,3)
    # @test (output[3] .<= 1.1) == trues(3,3)
    # @test (0.0 .<= output[3]) == trues(3,3)
    # @test (output[4] .<= 1.1) == trues(3,3)
    # @test (0.0 .<= output[4]) == trues(3,3)
    # @test (output[4] .<= 1.1) == trues(3,3)
    # @test (0.0 .<= output[4]) == trues(3,3)

    output = ArrayOutput(pFreqinit; tspan=1:5, aux=(exposure=exposure,))
    rule = Ruleset(DeltaAlleleFrequencySurvMap_noFastmath(layerkey=Val(:exposure),
                    LC50=10, hillcoefficient=2.0, deviationPhenotype = 10.0, dominanceDegree = 1))
    sim!(output, rule)

end


text:
using DiffEqFlux, Flux, Zygote

Zygote.@nograd diagm
include("eig.jl")

""" LDS condition for CTRNN0
        Enforce Ω(I - τA) + (I - τA')Ω ≻ 0,
        parametrize sqrt(Ω) """
struct ISSCTRNN0{P,Q1,Q2,Q3,F} <: Modifier
    name::String
    param_size::Int
    state_dim::Int
    initial_params::P
    get_Ω::Q1
    get_sqrtΩ::Q2
    get_sqrtΩinv::Q3
    ϵ::F
    lb::F
    logτ_idx::UnitRange{Int}
    A_idx::UnitRange{Int}
    enable::Bool
    function ISSCTRNN0(
        state_dim;
        enable = true,
        is_fixed = false,
        initb = Flux.zeros,
        ϵ = 1e-3,
        scale = 0.5,
        shift = -3.0,
    )
        # initial parameters
        initial_params() = vec(initb(state_dim))

        # function to fetch parameters: if fixed, ignore parameters
        if is_fixed
            Ω_init = ones(state_dim)
            get_Ω = (p) -> Ω_init
            get_sqrtΩ = (p) -> Ω_init
            get_sqrtΩinv = (p) -> Ω_init
        else
            # get_sqrtΩ = (p) -> 0.1 .+ log.(1.0 .+ p .^ 2) #Flux.huber_loss(zeros(state_dim), p, agg=identity)
            get_sqrtΩ = (p) -> softmax(p)
            get_sqrtΩinv = (p) -> 1.0 ./ get_sqrtΩ(p)

            # get_sqrtΩ = (p) -> scale*(1.0.+exp.(-p.+shift))
            # get_sqrtΩinv = (p) -> sigmoid.(p.-shift)/scale
            get_Ω = (p) -> get_sqrtΩ(p) .^ 2
        end

        # parameter indices: assume computed same way as by CTRNN0
        i = 1
        logτ_idx = i:i
        i += 1
        A_idx = i:i+state_dim*state_dim-1

        # logτ_idx = subsystem.index_dict["logτ"]
        # A_idx = subsystem.index_dict["A"]

        new{
            typeof(initial_params),
            typeof(get_Ω),
            typeof(get_sqrtΩ),
            typeof(get_sqrtΩinv),
            typeof(ϵ),
        }(
            "ISSCTRNN0",
            state_dim,
            state_dim,
            initial_params,
            get_Ω,
            get_sqrtΩ,
            get_sqrtΩinv,
            ϵ,
            shift,
            logτ_idx,
            A_idx,
            enable,
        )
    end
end

""" Modify p_sub to guarantee stability """
function (f::ISSCTRNN0)(p_sub, p_link, p; T = Float64, rhs = T(f.ϵ))
    if f.enable
        # Ω = f.get_Ω(p)
        sqrtΩ = f.get_sqrtΩ(p)
        invsqrtΩ = f.get_sqrtΩinv(p)

        τ = exp.(@view p_sub[f.logτ_idx])
        A = reshape(view(p_sub, f.A_idx), f.state_dim, f.state_dim)

        mat = T(0.5) * τ .* ((sqrtΩ .* A)' .* invsqrtΩ)'
        factor = one(T) / (one(T) + relu(maxeig(mat + mat') - one(T) + T(f.ϵ)))

        _p_sub = vcat(p_sub[1], factor * vec(A), p_sub[f.A_idx[end]+1:end])

        return _p_sub, p_link
    else
        return p_sub, p_link
    end
end

"""  Debug: Computes the largest eigenvalue to ensure Ω(τA - I) + (τA' - I)Ω ≺ 0 """
function debug(f::ISSCTRNN0, p_sub, p_link, p)
    Ω = f.get_Ω(p)

    τ = exp.(@view p_sub[f.logτ_idx])
    A = reshape(view(p_sub, f.A_idx), f.state_dim, f.state_dim)

    mat = τ .* Ω .* A - diagm(Ω)
    return maxeig(mat + mat')
end

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

""" LDS condition for CTRNN0
        Enforce Ω(I - τA) + (I - τA')Ω ≻ 0,
        parametrize sqrt(Ω) """
struct ISSCTRNN0Alt{P,Q1,Q2,Q3,F} <: Modifier
    name::String
    param_size::Int
    state_dim::Int
    initial_params::P
    get_Ω::Q1
    get_sqrtΩ::Q2
    get_sqrtΩinv::Q3
    ϵ::F
    lb::F
    logτ_idx::UnitRange{Int}
    A_idx::UnitRange{Int}
    enable::Bool
    function ISSCTRNN0Alt(
        state_dim;
        enable = true,
        is_fixed = false,
        initb = Flux.zeros,
        ϵ = 1e-3,
        lb = -3.0,
    )
        # initial parameters
        initial_params() = vec(initb(state_dim))

        # function to fetch parameters: if fixed, ignore parameters
        if is_fixed
            Ω_init = ones(state_dim)
            get_Ω = (p) -> Ω_init
            get_sqrtΩ = (p) -> Ω_init
            get_sqrtΩinv = (p) -> Ω_init
        else
            get_sqrtΩ = (p) -> 0.5 .+ exp.(p) #log.(1.0.+p.^2)
            get_sqrtΩinv = (p) -> 1.0 ./ get_sqrtΩ(p)
            get_Ω = (p) -> get_sqrtΩ(p) .^ 2
        end

        # parameter indices: assume computed same way as by CTRNN0
        i = 1
        logτ_idx = i:i
        i += 1
        A_idx = i:i+state_dim*state_dim-1

        # logτ_idx = subsystem.index_dict["logτ"]
        # A_idx = subsystem.index_dict["A"]

        new{
            typeof(initial_params),
            typeof(get_Ω),
            typeof(get_sqrtΩ),
            typeof(get_sqrtΩinv),
            typeof(ϵ),
        }(
            "ISSCTRNN0Alt",
            state_dim,
            state_dim,
            initial_params,
            get_Ω,
            get_sqrtΩ,
            get_sqrtΩinv,
            ϵ,
            lb,
            logτ_idx,
            A_idx,
            enable,
        )
    end
end

""" Modify p_sub to guarantee stability """
function (f::ISSCTRNN0Alt)(p_sub, p_link, p; T = Float64, rhs = T(f.ϵ))
    if f.enable
        sqrtΩ = f.get_sqrtΩ(p)
        Ω = f.get_Ω(p)

        τ = exp.(@view p_sub[f.logτ_idx])
        A = reshape(view(p_sub, f.A_idx), f.state_dim, f.state_dim)

        # mat = T(0.5)*τ.*((sqrtΩ.*A)'.*sqrtΩ)'
        mat = T(0.5) * τ .* sqrtΩ .* A
        factor = one(T) / (one(T) + relu(maxeig(mat + mat') - one(T) + T(f.ϵ)))

        _p_sub = vcat(p_sub[1], factor * vec((A' .* sqrtΩ)'), p_sub[f.A_idx[end]+1:end])

        return _p_sub, p_link
    else
        return p_sub, p_link
    end
end

"""  Debug: Computes the largest eigenvalue to ensure Ω(τA - I) + (τA' - I)Ω ≺ 0 """
function debug(f::ISSCTRNN0Alt, p_sub, p_link, p)
    Ω = f.get_Ω(p)

    τ = exp.(@view p_sub[f.logτ_idx])
    A = reshape(view(p_sub, f.A_idx), f.state_dim, f.state_dim)

    mat = τ .* Ω .* A - diagm(Ω)
    return maxeig(mat + mat')
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
InputFolder = './Images/NucleiNoisy/';
OutputFolder = './Results/Images/NucleiNoisy3/';

@iA = '*.tif';		% Image filter

@fxg_gDenoiseNLM [iA] > [D, An];
params.PatchSizeHalf = 5;
params.WindowSizeHalf = 7;
params.Sigma = 0.3;
params.TopOpen = 35;
params.AddNoiseVar = 0;
/endf

/show iA >;
/show D >;
/keep D > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
function test_simple_net_gradient(backend)
    println("-- Testing gradients on simple network (example for gradient checking code)")

    srand(12345678)
    ############################################################
    # Prepare Random Data
    ############################################################
    N = 5    # works with arbitrary minibatch size as long as
             # N == batch_size in MemoryDataLayer so it cycles through
             # and gets the same data during forward()
    M = 10
    P = 4

    X = rand(M, N)
    W = rand(M, P)
    B = rand(P, 1)

    Y = (W'*X .+ B)
    Y = Y + 0.01*randn(size(Y))

    ############################################################
    # Define network
    ############################################################
    data_layer = MemoryDataLayer(batch_size=N, data=Array[X,Y])

    w1 = InnerProductLayer(neuron=Neurons.Sigmoid(), name="ip1",output_dim=20, tops=[:a], bottoms=[:data])
    w2 = InnerProductLayer(neuron=Neurons.Identity(), name="ip2",output_dim=4, tops=[:b], bottoms=[:a])
    loss_layer = SquareLossLayer(name="loss", bottoms=[:b, :label] )

    net = Net("TEST", backend, [w1,w2, loss_layer, data_layer])

    # epsilon:     milage may vary 1e-4 - 1e-8
    # digit:       compare this many digits to check for 'identity'
    # visualize:   prints out correct and failed positions of gradients
    return test_gradients(net, epsilon=1e-8, digit=6, visual=false )
end

if test_cpu
    test_simple_net_gradient(backend_cpu)
end


text:
include("structs.jl")
include("cfd.jl")

function dump_total_liquid_flux(sys::System, p::Parameters,
        t::Integer)
    flux_init  = get_flux(sys.ρ[:, :, 2],
                          sys.s[:, :, 2],
                          sys.v[:, :, 2], 1)
    flux_final = get_flux(sys.ρ[:, :, 2],
                          sys.s[:, :, 2],
                          sys.v[:, :, 2], p.ny - 1)

    total_flux_init  = sum(flux_init)
    total_flux_final = sum(flux_final)

    if t == 1
        io_init  = open("../dump/total_flux_init.txt",  "w")
        io_final = open("../dump/total_flux_final.txt", "w")
    else
        io_init  = open("../dump/total_flux_init.txt",  "a+")
        io_final = open("../dump/total_flux_final.txt", "a+")
    end
    writedlm(io_init,  total_flux_init,  '\n')
    writedlm(io_final, total_flux_final, '\n')
end

function dump(sys::System, p::Parameters)
    # x = ideal_gas
    # pressure1 = @. R * x.T * sys.ρ[:, :, 1] / p.M[1] / sys.s
    # x = tait_C₅H₁₂
    # ρ̂₂ = @. sys.ρ[:, :, 2] / (1 - sys.s)
    # pressure2 = @. (x.B + x.P₀) * 10^((ρ̂₂ - x.ρ₀) / x.C /
    #                                   ρ̂₂) - x.B

    @views ρ̂ = sys.ρ ./ sys.s

    writedlm("../dump/pressure.txt", sys.P[:, :], ' ')
    # writedlm("../dump/pressure1.txt", pressure1 .- sys.P, ' ')
    # writedlm("../dump/pressure2.txt", pressure2 .- sys.P, ' ')
    writedlm("../dump/v1.txt", sys.v[:, :, 1], ' ')
    writedlm("../dump/v2.txt", sys.v[:, :, 2], ' ')
    writedlm("../dump/density1.txt", ρ̂[:, :, 1], ' ')
    writedlm("../dump/density2.txt", ρ̂[:, :, 2], ' ')
    writedlm("../dump/flow1.txt", sys.ρ[:, :, 1] .\
             sys.s[:, :, 1] .* sys.v[:, :, 1], ' ')
    writedlm("../dump/flow2.txt", sys.ρ[:, :, 2] .\
             sys.s[:, :, 2] .* sys.v[:, :, 2], ' ')
    writedlm("../dump/saturation.txt", sys.s[:, :, 1], ' ')
end

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration