text:
using Combinatorics
using SymPy

simplify(x::Real) = x
simplify(x::SymPy.Sym) = x
atan2(y::SymPy.Sym, x::SymPy.Sym) = atan(y, x)

abstract type AnalyticEquilibrium <: MagneticEquilibrium end
abstract type AnalyticPerturbation <: MagneticEquilibrium end

X(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("X() not implemented for ", ET)
Y(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("Y() not implemented for ", ET)
Z(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("Z() not implemented for ", ET)
R(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("R() not implemented for ", ET)
r(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("r() not implemented for ", ET)
θ(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("θ() not implemented for ", ET)
ϕ(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("ϕ() not implemented for ", ET)

J(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = error("J() not implemented for ", ET)

periodicity(x::AbstractArray{T,1}, equ::ET) where {T, ET <: MagneticEquilibrium} = zero(x)

A₁(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
A₂(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
A₃(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)

φ(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)

g₁₁(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₁₂(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₁₃(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₂₁(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₂₂(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₂₃(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₃₁(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₃₂(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)
g₃₃(x::AbstractArray{T,1}, equ::MagneticEquilibrium) where {T} = zero(T)

function A(x, equ)
    [A₁(x,equ), A₂(x,equ), A₃(x,equ)]
end

function g(x, equ)
    [g₁₁(x,equ)  g₁₂(x,equ)  g₁₃(x,equ);
     g₂₁(x,equ)  g₂₂(x,equ)  g₂₃(x,equ);
     g₃₁(x,equ)  g₃₂(x,equ)  g₃₃(x,equ);]
end

struct ZeroPerturbation <: AnalyticPerturbation
    name::String
    ZeroPerturbation() = new("ZeroPerturbation")
end

A₁(x::AbstractArray{T,1}, pert::ZeroPerturbation) where {T} = zero(T)
A₂(x::AbstractArray{T,1}, pert::ZeroPerturbation) where {T} = zero(T)
A₃(x::AbstractArray{T,1}, pert::ZeroPerturbation) where {T} = zero(T)

"Returns the i-th component of the vector corresponding to the one-form α"
function get_vector_component(α, g̅, i)
    simplify(factor(g̅[i,1] * α[1] + g̅[i,2] * α[2] + g̅[i,3] * α[3]))
end

"Returns the m-th component of the one-form corresponding to the two-form β"
function hodge²¹(β, g̅, J, m)
    α = Sym(0)

    for i in 1:3
        for j in 1:3
            for k in 1:3
                for l in 1:3
                    α += β[i,j] * g̅[i,k] * g̅[j,l] * levicivita([k,l,m])
                end
            end
        end
    end

    return simplify(factor(J*α))
end

"""
Generate functions for evaluating analytic equilibria.
"""
function generate_equilibrium_functions(equ::AnalyticEquilibrium, pert::AnalyticPerturbation; output=0)
    # define symbols for coordinates x = (x₁, x₂, x₃),
    # positive=true is set so that sqrt(x^2) does not become |x^2|
    x₁, x₂, x₃ = symbols("x₁, x₂, x₃", real=true, positive=true)
    x = [x₁, x₂, x₃]

    # check for compatible metric
    if typeof(pert) != ZeroPerturbation
        @assert J(x,equ) == J(x,pert)
        @assert g(x,equ) == g(x,pert)
    end

    # obtain metric
    gmat = g(x, equ)
    symprint("g", gmat, output, 2)

    # invert metric
    ginv = inv(gmat)
    for i in 1:3
        for j in 1:3
            ginv[i,j] = simplify(factor(ginv[i,j]))
        end
    end
    symprint("g⁻¹", ginv, output, 2)

    # compute Jacobian determinant
    # Jdet² = factor(det(gmat))
    # symprint("J²", Jdet², output, 2)

    # Jdet = sqrt(Jdet²)
    Jdet = J(x, equ)
    symprint("J", Jdet, output, 2)

    # obtain vector potential
    A¹ = A(x, equ) .+ A(x, pert)
    symprint("A¹", A¹, output, 2)

    # compute vector potential in contravariant coordinates
    Avec = [get_vector_component(A¹, ginv, i) for i in 1:3]
    symprint("Avec", Avec, output, 2)

    # compute Jacobian of vector potential A
    DA = [simplify(diff(A¹[i], x[j])) for i in 1:3, j in 1:3]
    symprint("DA", DA, output, 3)

    # compute second derivative of vector potential A
    DDA1 = [simplify(factor(diff(diff(A¹[1], x[i]), x[j]))) for i in 1:3, j in 1:3]
    symprint("DDA1", DDA1, output, 3)

    DDA2 = [simplify(factor(diff(diff(A¹[2], x[i]), x[j]))) for i in 1:3, j in 1:3]
    symprint("DDA2", DDA2, output, 3)

    DDA3 = [simplify(factor(diff(diff(A¹[3], x[i]), x[j]))) for i in 1:3, j in 1:3]
    symprint("DDA3", DDA3, output, 3)

    # compute components of magnetic field B
    Bᶜ = [simplify(factor(diff(A¹[3], x[2]) - diff(A¹[2], x[3]))),
          simplify(factor(diff(A¹[1], x[3]) - diff(A¹[3], x[1]))),
          simplify(factor(diff(A¹[2], x[1]) - diff(A¹[1], x[2])))]
    symprint("Bᶜ", Bᶜ, output, 2)

    # compute magnetic field two-form B²
    B² = [ 0      +Bᶜ[3]    -Bᶜ[2];
          -Bᶜ[3]   0        +Bᶜ[1];
          +Bᶜ[2]   -Bᶜ[1]    0    ] .* Rational(1,2)
    symprint("B²", B², output, 2)

    # compute magnetic field one-form B¹ = ⋆B²
    B¹ = [hodge²¹(B², ginv, Jdet, i) for i in 1:3]
    symprint("B¹", B¹, output, 2)

    # compute magnetic field in contravariant coordinates
    Bvec = [get_vector_component(B¹, ginv, i) for i in 1:3]
    symprint("Bvec", Bvec, output, 2)

    # compute absolute value |B| of B
    Babs = simplify(sqrt(factor(B¹[1] * Bvec[1] + B¹[2] * Bvec[2] + B¹[3] * Bvec[3])))
    symprint("|B|", Babs, output, 2)

    # compute magnetic unit one-form
    b¹ = [simplify(factor( B¹[i] / Babs )) for i in 1:3]
    symprint("b¹", b¹, output, 2)

    # compute unit magnetic field in contravariant coordinates
    bvec = [get_vector_component(b¹, ginv, i) for i in 1:3]
    symprint("bvec", bvec, output, 2)

    # compute Jacobian of magnetic field B
    DB = [diff(B¹[i], x[j]) for i in 1:3, j in 1:3]
    symprint("DB", DB, output, 3)

    # compute Jacobian of magnetic unit vector b
    Db = [diff(b¹[i], x[j]) for i in 1:3, j in 1:3]
    symprint("Db", Db, output, 3)

    # compute second derivative of magnetic unit vector b
    DDb1 = [diff(diff(b¹[1], x[i]), x[j]) for i in 1:3, j in 1:3]
    symprint("DDb1", DDb1, output, 3)

    DDb2 = [diff(diff(b¹[2], x[i]), x[j]) for i in 1:3, j in 1:3]
    symprint("DDb2", DDb2, output, 3)

    DDb3 = [diff(diff(b¹[3], x[i]), x[j]) for i in 1:3, j in 1:3]
    symprint("DDb3", DDb3, output, 3)

    # compute first derivatives of absolute value of magnetic field
    DBabs = [diff(Babs, x[j]) for j in 1:3]
    symprint("D|B|", DBabs, output, 3)

    # compute second derivatives of absolute value of magnetic field
    DDBabs = [diff(diff(Babs, x[i]), x[j]) for i in 1:3, j in 1:3]
    symprint("DD|B|", DDBabs, output, 3)

    # obtain scalar potential
    φ⁰ = φ(x, equ) .+ φ(x, pert)

    # compute components of electric field E
    E¹ = [simplify(factor(diff(-φ⁰, x[1]))),
          simplify(factor(diff(-φ⁰, x[2]))),
          simplify(factor(diff(-φ⁰, x[3])))]
    symprint("E¹", E¹, output, 2)

    # compute electric field in contravariant coordinates
    Evec = [get_vector_component(E¹, ginv, i) for i in 1:3]
    symprint("Evec", Evec, output, 2)

    # collect all functions to generate code for
    functions = Dict{String,Any}()
    indices   = ["₁", "₂", "₃"]
    indicesup = ["¹", "²", "³"]

    functions["X"] = X(x, equ)
    functions["Y"] = Y(x, equ)
    functions["Z"] = Z(x, equ)

    try functions["R"] = R(x, equ) catch end
    try functions["r"] = r(x, equ) catch end
    try functions["θ"] = θ(x, equ) catch end
    try functions["ϕ"] = ϕ(x, equ) catch end

    functions["J"] = Jdet
    functions["B"] = Babs
    functions["φ"] = φ⁰

    for i in 1:3
        functions["A" * indices[i]] = A¹[i]
        functions["B" * indices[i]] = B¹[i]
        functions["b" * indices[i]] = b¹[i]
        functions["E" * indices[i]] = E¹[i]

        functions["A" * indicesup[i]] = Avec[i]
        functions["B" * indicesup[i]] = Bvec[i]
        functions["b" * indicesup[i]] = bvec[i]
        functions["E" * indicesup[i]] = Evec[i]

        functions["dBdx" * indices[i]] = DBabs[i]
    end

    for i in 1:3
        for j in 1:3
            functions["g"  * indices[i]   * indices[j]]   = gmat[i,j]
            functions["g"  * indicesup[i] * indicesup[j]] = ginv[i,j]

            functions["B"  * indices[i]   * indices[j]]   = B²[i,j]

            functions["dA" * indices[i] * "dx" * indices[j]] = DA[i,j]
            functions["dB" * indices[i] * "dx" * indices[j]] = DB[i,j]
            functions["db" * indices[i] * "dx" * indices[j]] = Db[i,j]

            functions["d²A₁" * "dx" * indices[i] * "dx" * indices[j]] = DDA1[i,j]
            functions["d²A₂" * "dx" * indices[i] * "dx" * indices[j]] = DDA2[i,j]
            functions["d²A₃" * "dx" * indices[i] * "dx" * indices[j]] = DDA3[i,j]
            functions["d²b₁" * "dx" * indices[i] * "dx" * indices[j]] = DDb1[i,j]
            functions["d²b₂" * "dx" * indices[i] * "dx" * indices[j]] = DDb2[i,j]
            functions["d²b₃" * "dx" * indices[i] * "dx" * indices[j]] = DDb3[i,j]
            functions["d²B"  * "dx" * indices[i] * "dx" * indices[j]] = DDBabs[i,j]
        end
    end

    functions
end

"""
Generate code for evaluating analytic equilibria.
"""
function generate_equilibrium_code(equ, pert=ZeroPerturbation(); output=0)

    if output ≥ 1
        println("Generating code for ")
        println(equ)
        if typeof(pert) != ZeroPerturbation
            println("   and ")
            println(pert)
        end
        println()
    end

    parameters = fieldnames(typeof(equ))
    functions = generate_equilibrium_functions(equ, pert; output=output)

    equ_code = :(  )

    # generate Julia code and export parameters
    for param in parameters
        if param != :name
            value = getfield(equ, param)

            p_code = quote
                export $param
                $param = $value
            end

            # append p_code to equ_code
            push!(equ_code.args, p_code)
        end
    end

    # generate Julia code and export functions
    for (key, value) in functions
        f_symb = Symbol(key)
        f_expr = value

        output ≥ 1 ? println("Generating function ", key) : nothing

        # patch for removing Warnings with Julia v0.6 and SymPy v1.0
        f_str  = sympy.julia_code(f_expr)
        f_str  = replace(f_str, ".+" => " .+ ")
        f_str  = replace(f_str, ".-" => " .- ")
        f_str  = replace(f_str, ".*" => " .* ")
        f_str  = replace(f_str, "./" => " ./ ")
        f_str  = replace(f_str, ".^" => " .^ ")
        f_body = Meta.parse(f_str)

        output ≥ 2 ? println("   ", f_body) : nothing

        f_code = quote
            export $f_symb
            @inline function $f_symb(x₁, x₂, x₃)
                $f_body
            end
            @inline function $f_symb(t::Number, x::AbstractArray{T,1}) where {T <: Number}
                $f_symb(x[1],x[2],x[3])
            end
        end

        # append f_code to equ_code
        push!(equ_code.args, f_code)
    end

    output ≥ 1 ? println() : nothing

    return equ_code
end

"""
Evaluate functions for evaluating analytic equilibria.
"""
function load_equilibrium(equ, pert=ZeroPerturbation(); target_module=Main, output=0)
    equ_code = generate_equilibrium_code(equ, pert; output=output)
    Core.eval(target_module, equ_code)
end

function symprint(name, symexpr, output=1, detail_level=0)
    if output ≥ detail_level
        println(name, " = ", symexpr, "\n")
    end
end


text:
import IndexedTables: flatten

function mapslices(f, x::DNDSparse, dims; name=nothing)
    iterdims = setdiff([1:ndims(x);], map(d->keyindex(x, d), dims))
    if iterdims != [1:length(iterdims);]
        throw(ArgumentError("$dims must be the trailing dimensions of the table. You can use `permutedims` first to permute the dimensions."))
    end

    # Note: the key doesn't need to be put in a tuple, this is
    # also bad for sortperm, but is required since DArrays aren't
    # parameterized by the container type Columns
    vals = isempty(dims) ?  values(x) : (keys(x, (dims...)), values(x))
    tmp = ndsparse((keys(x, (iterdims...)),), vals,
                   allowoverlap=false, closed=true)

    cs = delayedmap(tmp.chunks) do c
        ks = isempty(dims) ? columns(columns(keys(c))[1]) : IndexedTables.concat_cols(columns(keys(c))[1], columns(values(c))[1])
        vs = isempty(dims) ? columns(values(c)) : columns(values(c))[2]
        y = ndsparse(ks, vs)
        mapslices(f, y, dims; name=name)
    end
    fromchunks(cs)
  # cache_thunks(mapchunks(y -> mapslices(f, y, dims, name=name),
  #                        t, keeplengths=false))
end

mapslices(f, x::DNDSparse, dims::Symbol; name=nothing) =
    mapslices(f, x, (dims,); name=name)

function flatten(x::DNextTable, col)
    fromchunks(delayedmap(t -> flatten(t, col), x.chunks))
end

ext: jl


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test
boundary_condition = BoundaryConditionDirichlet(initial_condition)

source_terms = source_terms_convergence_test

solver = DGSEM(polydeg=4, surface_flux=flux_hll)

# Gnomonic projection similar to how the coupled cubed sphere is created.
# Note that half the elements will have left-handed coordinates.
function gnomonic_projection(xi, eta, inner_radius, thickness, direction)
  alpha = xi * pi/4

  # Equiangular projection
  x = tan(alpha)

  # Coordinates on unit square per direction
  square_coordinates = [SVector(-1, x),
                        SVector( 1, x),
                        SVector(x, -1),
                        SVector(x,  1)]

  # Radius on square surface
  r = sqrt(1 + x^2)

  # Radius of the sphere
  R = inner_radius + thickness * (0.5 * (eta + 1))

  # Projection onto the sphere
  R / r * square_coordinates[direction]
end

function ring_mapping(inner_radius, thickness, direction)
  mapping(xi, eta) = gnomonic_projection(xi, eta, inner_radius, thickness, direction)
end

mapping_as_string(direction) = """
function gnomonic_projection(xi, eta, inner_radius, thickness, direction)
  alpha = xi * pi/4

  x = tan(alpha)

  r = sqrt(1 + x^2)
  R = inner_radius + thickness * (0.5 * (eta + 1))

  # Cube coordinates per direction
  cube_coordinates = [SVector(-1, x),
                      SVector( 1, x),
                      SVector(x, -1),
                      SVector(x,  1)]

  R / r * cube_coordinates[direction]
end; function ring_mapping(inner_radius, thickness, direction)
  mapping(xi, eta) = gnomonic_projection(xi, eta, inner_radius, thickness, direction)
end; mapping = ring_mapping(1, 1, $direction)
"""

cells_per_dimension = (8, 4)

mesh1 = StructuredMesh(cells_per_dimension, ring_mapping(1, 1, 1),
                   periodicity=false, mapping_as_string=mapping_as_string(1))

semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition, solver,
  source_terms=source_terms, boundary_conditions=(
    x_neg=BoundaryConditionCoupled(3, (:begin, :i_forward), Float64),
    x_pos=BoundaryConditionCoupled(4, (:begin, :i_forward), Float64),
    y_neg=boundary_condition,
    y_pos=boundary_condition,
  ))

mesh2 = StructuredMesh(cells_per_dimension, ring_mapping(1, 1, 2),
                   periodicity=false, mapping_as_string=mapping_as_string(2))

semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition, solver,
  source_terms=source_terms, boundary_conditions=(
    x_neg=BoundaryConditionCoupled(3, (:end, :i_forward), Float64),
    x_pos=BoundaryConditionCoupled(4, (:end, :i_forward), Float64),
    y_neg=boundary_condition,
    y_pos=boundary_condition,
  ))

mesh3 = StructuredMesh(cells_per_dimension, ring_mapping(1, 1, 3),
                   periodicity=false, mapping_as_string=mapping_as_string(3))

semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition, solver,
  source_terms=source_terms, boundary_conditions=(
    x_neg=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64),
    x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64),
    y_neg=boundary_condition,
    y_pos=boundary_condition,
  ))

mesh4 = StructuredMesh(cells_per_dimension, ring_mapping(1, 1, 4),
                   periodicity=false, mapping_as_string=mapping_as_string(4))

semi4 = SemidiscretizationHyperbolic(mesh4, equations, initial_condition, solver,
  source_terms=source_terms, boundary_conditions=(
    x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64),
    x_pos=BoundaryConditionCoupled(2, (:end, :i_forward), Float64),
    y_neg=boundary_condition,
    y_pos=boundary_condition,
  ))

semi = SemidiscretizationCoupled((semi1, semi2, semi3, semi4))

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)
###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
const SELECTION         = Dict{Symbol, Input{Matrix{Vector2{Int}}}}()
const SELECTION_QUERIES = Dict{Symbol, Rectangle{Int}}()

function insert_selectionquery!(name::Symbol, value::Rectangle)
  SELECTION_QUERIES[name] = value
  SELECTION[name]         = Input(Vector2{Int}[]')
  SELECTION[name]
end
function insert_selectionquery!(name::Symbol, value::Signal{Rectangle{Int}})
  lift(value) do v
    SELECTION_QUERIES[name] = v
  end
  SELECTION[name]         = Input(Array(Vector2{Int}, value.value.w, value.value.h))
  SELECTION[name]
end
function delete_selectionquery!(name::Symbol)
  delete!(SELECTION_QUERIES, name)
  delete!(SELECTION, name)
end

windowhints = [
  (GLFW.SAMPLES, 0),
  (GLFW.DEPTH_BITS, 0),
  (GLFW.ALPHA_BITS, 0),
  (GLFW.STENCIL_BITS, 0),
  (GLFW.AUX_BUFFERS, 0)
]

const ROOT_SCREEN = createwindow("Romeo", 1920, 1280, windowhints=windowhints, debugging=false)
insert_selectionquery!(:mouse_hover, lift(ROOT_SCREEN.inputs[:mouseposition]) do mpos
  Rectangle(int(mpos[1]), int(mpos[2]), 1,1)
end)

parameters = [
  (GL_TEXTURE_WRAP_S,  GL_CLAMP_TO_EDGE),
  (GL_TEXTURE_WRAP_T,  GL_CLAMP_TO_EDGE ),

  (GL_TEXTURE_MIN_FILTER, GL_NEAREST),
  (GL_TEXTURE_MAG_FILTER, GL_NEAREST)
]

global const RENDER_FRAMEBUFFER = glGenFramebuffers()
glBindFramebuffer(GL_FRAMEBUFFER, RENDER_FRAMEBUFFER)

framebuffsize = [ROOT_SCREEN.inputs[:framebuffer_size].value]

color   = Texture(RGBA{Ufixed8},     framebuffsize, parameters=parameters)
stencil = Texture(Vector2{GLushort}, framebuffsize, parameters=parameters)

glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, color.id, 0)
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, stencil.id, 0)

rboDepthStencil = GLuint[0]

glGenRenderbuffers(1, rboDepthStencil)
glBindRenderbuffer(GL_RENDERBUFFER, rboDepthStencil[1])
glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT32, framebuffsize...)
glFramebufferRenderbuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, rboDepthStencil[1])

lift(ROOT_SCREEN.inputs[:framebuffer_size]) do window_size
  if all(x->x>0, window_size)
    resize!(color, window_size)
    resize!(stencil, window_size)
    glBindRenderbuffer(GL_RENDERBUFFER, rboDepthStencil[1])
    glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT32, (window_size)...)
  end
end

function renderloop(ROOT_SCREEN)
  global ROOT_SCREEN
  yield()
  glBindFramebuffer(GL_FRAMEBUFFER, RENDER_FRAMEBUFFER)
  glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

  render(ROOT_SCREEN)

  #Read all the selection queries
  if !isempty(SELECTION_QUERIES)
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    for (key, value) in SELECTION_QUERIES
      if value.w < 1 || value.w > 5000
        println(value.w)
      end
      if value.h < 1 || value.h > 5000
        println(value.h)
      end
      const data = Array(Vector2{Uint16}, value.w, value.h)
      glReadPixels(value.x, value.y, value.w, value.h, stencil.format, stencil.pixeltype, data)
      push!(SELECTION[key], convert(Matrix{Vector2{Int}}, data))
    end
  end

  glReadBuffer(GL_COLOR_ATTACHMENT0)
  glBindFramebuffer(GL_READ_FRAMEBUFFER, RENDER_FRAMEBUFFER)
  glBindFramebuffer(GL_DRAW_FRAMEBUFFER, 0)
  glClear(GL_COLOR_BUFFER_BIT)
  glDisable(GL_SCISSOR_TEST)
  ROOT_SCREEN_size = ROOT_SCREEN.inputs[:framebuffer_size].value
  glBlitFramebuffer(0,0, ROOT_SCREEN_size..., 0,0, ROOT_SCREEN_size..., GL_COLOR_BUFFER_BIT, GL_LINEAR)
  GLFW.SwapBuffers(ROOT_SCREEN.nativewindow)
  GLFW.PollEvents()
end

glClearColor(0, 0, 0, 0)


text:
@testset "mountain_car" begin

    env = MountainCarEnv(; rng = StableRNG(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    env = ContinuousMountainCarEnv(; rng = StableRNG(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
using QuantumWalk
using LightGraphs
using Test
using LinearAlgebra,SparseArrays

macro test_no_error(ex)
    quote
        try
            $(esc(ex))
            true
        catch
            false
        end
    end
end

include("ctqw.jl")
include("szegedy.jl")
include("qwsearch.jl")
include("qwevolution.jl")
include("type_hierarchy.jl")
include("maximizing_function.jl")


text:
using WAVEFORM, LinearAlgebra

# Partitioning
n = 96
partition_size = 48
overlap = 0

y = partition(n, partition_size, overlap)

@test size(y) == (partition_size, Int(n / partition_size))

# Grid interpolator
n_fine = (100, 100)
n_coarse = (50, 50)
T = Float64
(f2c,c2f) = fine2coarse(n_fine,n_coarse,T)

x = randn(100, 100)
y = randn(50, 50)

@test dot(vec(y), f2c*vec(x)) - dot(vec(x), f2c'*vec(y)) < 1e-9
@test dot(vec(y), c2f'*vec(x)) - dot(vec(x), c2f*vec(y)) < 1e-9


text:
using AWS

config_file = isempty(ARGS) ? "config.jl" : ARGS[1]
include(config_file)

include("tests3.jl")
include("testsqs.jl")
include("testautoscale.jl")
include("ec2_utils.jl")
#include("testec2.jl")

const dbg = get(config, "dbg", false)
const key = get(config, "key", AWS.AWS_SECKEY)
const id = get(config, "key", AWS.AWS_ID)
const region = get(config, "region", AWS.AWS_REGION)

env = AWSEnv(timeout=60.0, id=id, key=key, region=region, dbg=dbg)

TestS3.runtests(env, config)
TestSQS.runtests(env, config)
TestEC2Utils.runtests(env, config)
TestAutoscale.runtests(env, config)


text:
using SafeTestsets

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

    @testset "run_program" begin
        using AdventOfCode.Computer
        @test run_program!([1002,4,3,4,33]).mem.vect == [1002,4,3,4,99]

        p = Program([3,0,4,0,99])
        run_program!(p, [42])
        @test p.output == 42

        function equals8(i)
            p = Program([3,9,8,9,10,9,4,9,99,-1,8])
            run_program!(p, i)
            p.output==1
        end
        @test !equals8(7)
        @test  equals8(8)
        @test !equals8(9)

        function lessthan8(i)
            p = Program([3,9,7,9,10,9,4,9,99,-1,8])
            run_program!(p, i)
            p.output == 1
        end
        @test  lessthan8(7)
        @test !lessthan8(8)

        function equals8_(i)
            p = Program([3,3,1108,-1,8,3,4,3,99])
            run_program!(p, i)
            p.output == 1
        end
        @test !equals8_(7)
        @test  equals8_(8)
        @test !equals8_(9)

        function lessthan8_(i)
            p = Program([3,3,1107,-1,8,3,4,3,99])
            run_program!(p, i)
            p.output == 1
        end
        @test  lessthan8_(7)
        @test !lessthan8_(8)

        function non_zero(i)
            p = Program([3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9])
            run_program!(p, i)
            p.output
        end
        @test non_zero(0) == 0
        @test non_zero(1) == 1
        @test non_zero(2) == 1

        function non_zero_(i)
            p = Program([3,3,1105,-1,9,1101,0,0,12,4,12,99,1])
            run_program!(p, i)
            p.output
        end
        @test non_zero_(0) == 0
        @test non_zero_(1) == 1
        @test non_zero_(2) == 1

        function cmp8(i)
            p = Program([3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
                         1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
                         999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99])
            run_program!(p,i)
            p.output
        end
        @test cmp8(7) ==  999
        @test cmp8(8) == 1000
        @test cmp8(9) == 1001
    end

    btest("5.1") do
        part1() == 13787043
    end

    btest("5.2") do
        part2() == 3892695
    end
end


text:
using HDF5
using Base.Test

include("plain.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
include("external.jl")
include("swmr.jl")


text:
"""
    phase_teleportation(zxd)
Reducing T-count of `zxd` with the algorithms in [arXiv:1903.10477](https://arxiv.org/abs/1903.10477).
"""
function phase_teleportation(cir::ZXDiagram{T, P}) where {T, P}
    zxg = ZXGraph(cir)
    ncir = zxg.master

    simplify!(Rule{:lc}(), zxg)
    simplify!(Rule{:p1}(), zxg)
    simplify!(Rule{:p2}(), zxg)
    simplify!(Rule{:p3}(), zxg)
    simplify!(Rule{:p1}(), zxg)
    match_id = match(Rule{:id}(), zxg)
    match_gf = match(Rule{:gf}(), zxg)
    while length(match_id) + length(match_gf) > 0
        rewrite!(Rule{:id}(), zxg, match_id)
        rewrite!(Rule{:gf}(), zxg, match_gf)
        simplify!(Rule{:lc}(), zxg)
        simplify!(Rule{:p1}(), zxg)
        simplify!(Rule{:p2}(), zxg)
        simplify!(Rule{:p3}(), zxg)
        simplify!(Rule{:p1}(), zxg)
        match_id = match(Rule{:id}(), zxg)
        match_gf = match(Rule{:gf}(), zxg)
    end

    simplify!(Rule{:i1}(), ncir)
    simplify!(Rule{:i2}(), ncir)
    return ncir
end

function phase_teleportation(bir::BlockIR)
    zxd = convert_to_zxd(bir)
    nzxd = phase_teleportation(zxd)
    chain = convert_to_chain(nzxd)
    return BlockIR(bir.parent, bir.nqubits, chain)
end


text:
#! /usr/bin/julia
#
#	text_delete.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_text = ARGS[1]
key_in = ARGS[2]
println(fname_text)
println(key_in)
#
dict_aa = text_to_dict_proc(fname_text)
#
if ( haskey(dict_aa,key_in))
	delete!(dict_aa, key_in)
	dict_to_text_proc(dict_aa,fname_text)
end
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using GeoPhyInv
using Base.Test
using BenchmarkTools
using Calculus

const JD=GeoPhyInv.Data

fields=[:P]
# testing misfit and gradient calculation for TD objects
tgrid1=range(0., stop=1., length=5);
tgrid2=range(0., stop=1., length=10);

ageom=GeoPhyInv.Acquisition.AGeom_fixed(10,10,10,10,10,10,10,10)

x=GeoPhyInv.Data.TD_ones(fields,tgrid1,ageom)

yvec=randn(length(x))
xvec=similar(yvec)

@testset "testing copy!" begin
	@btime copy!(x, yvec);

	@btime copy!(xvec, x);

	@test xvec == yvec
end

@testset "simple LS error: x and y same time grid" begin
	y=GeoPhyInv.Data.TD_ones(fields,tgrid1,ageom)

	randn!(x)
	randn!(y)

	pa=GeoPhyInv.Data.P_misfit(x,y);
	@time GeoPhyInv.Data.func_grad!(pa,:dJx);
	gg1=vec(pa.dJx)

	xvec=vec(pa.x)
	gg2=Calculus.gradient(x -> (copy!(pa.x, x); return GeoPhyInv.Data.func_grad!(pa)), xvec)

	# check gradient with Finite Differencing
	@test gg1 ≈ gg2
end

#=

rrrrrr

# loop over same time grid and different time grid (interp_flag on/off)
for y in [GeoPhyInv.Data.TD_ones(fields,tgrid2,ageom), GeoPhyInv.Data.TD_ones(fields,tgrid1,ageom)]
	println("#########################################")

	# generate some random data
	randn!(y.d[1,1])
	randn!(x.d[1,1])

	for func_attrib in [:cls]
		coup=GeoPhyInv.Coupling.TD_delta(y.tgrid, [0.1,0.1], 0.0,  x.fields, x.ageom)
		randn!(coup.ssf[1,1])
		pa=GeoPhyInv.Data.P_misfit(x,y, func_attrib=func_attrib, coup=coup);

		xvec=vec(pa.x.d[1,1])
		gg2=Calculus.gradient(x -> err(x), xvec)

		@time GeoPhyInv.Data.func_grad!(pa,:dJx);
		gg1=vec(pa.dJx.d[1,1])
		# check gradient with Finite Differencing
		@test gg1 ≈ gg2

		function errw(xvec)

			for i in eachindex(pa.coup.ssf[1,1])
				pa.coup.ssf[1,1][i]=xvec[i]
			end

			return GeoPhyInv.Data.func_grad!(pa)
		end

		xvec=vec(pa.coup.ssf[1,1])
		gg2=Calculus.gradient(x -> errw(x), xvec)

		@time GeoPhyInv.Data.func_grad!(pa,:dJssf);
		gg1=vec(pa.dJssf[1,1])
		# check gradient with Finite Differencing
		@test gg1 ≈ gg2

	end
end

=#


text:
import Base: map, merge, filter

if isdefined(Base, :foreach)
    import Base.foreach
end

export map,
       probe,
       filter,
       filterwhen,
       foldp,
       sampleon,
       merge,
       previous,
       delay,
       droprepeats,
       flatten,
       bind!,
       bindmap!,
       unbind!,
       bound_srcs,
       bound_dests

"""
    map(f, s::Signal...) -> signal

Transform signal `s` by applying `f` to each element. For multiple signal arguments, apply `f` elementwise.
"""
function map(
        f, input::Signal, inputsrest::Signal...;
        init = f(map(value, (input, inputsrest...))...),
        typ = typeof(init), name = auto_name!("map", input, inputsrest...)
    )
    n = Signal(typ, init, (input, inputsrest...); name = name)
    connect_map(f, n, input, inputsrest...)
    n
end

function connect_map(f, output, inputs...)
    add_action!(output) do
        set_value!(output, f(map(value, inputs)...))
    end
end

probe(node, name, io=stderr) =
    preserve(map(x -> (println(io, name, " >! ", x); x), node))

"""
    foreach(f, inputs...)

Same as `map`, but will be prevented from gc until all the inputs have gone out of scope. Should be used in cases where `f` does a side-effect.
"""
foreach(f, in1::Signal, inputs::Signal...; kwargs...) = preserve(map(f, in1, inputs...; kwargs...))

"""
    filter(f, default, signal)

remove updates from the `signal` where `f` returns `false`. The filter will hold
the value default until f(value(signal)) returns true, when it will be updated
to value(signal).
"""
function filter(f::Function, default, input::Signal{T}; name=auto_name!("filter", input)) where T
    n = Signal(T, f(value(input)) ? value(input) : default, (input,); name=name)
    connect_filter(f, default, n, input)
    n
end

"""
    filter(f, signal)

remove updates from the `signal` where `f` returns `false`. The filter will hold
the current value of the signal until `f(value(signal))` returns true.
"""
filter(f::Function, input::Signal{T}; kwargs...) where {T} = filter(f, value(input), input; kwargs...)

function connect_filter(f, default, output, input)
    add_action!(output) do
        val = value(input)
        if f(val)
            set_value!(output, val)
        else
            deactivate!(output)
        end
    end
end

"""
    filterwhen(switch::Signal{Bool}, default, input)

Keep updates to `input` only when `switch` is true.

If switch is false initially, the specified default value is used.
"""
function filterwhen(predicate::Signal{Bool}, default, input::Signal{T};
                     name=auto_name!("filterwhen", predicate, input)) where T
    n = Signal(T, value(predicate) ? value(input) : default, (input,); name=name)
    connect_filterwhen(n, predicate, input)
    n
end

function connect_filterwhen(output, predicate, input)
    add_action!(output) do
        if value(predicate)
            set_value!(output, value(input))
        else
            deactivate!(output)
        end
    end
end

"""
    foldp(f, init, inputs...)

[Fold](http://en.wikipedia.org/wiki/Fold_(higher-order_function)) over past values.

Accumulate a value as the input signals change. `init` is the initial value of the accumulator.
`f` should take `1 + length(inputs)` arguments: the first is the current accumulated value and the rest are the current input signal values. `f` will be called when one or more of the `inputs` updates. It should return the next accumulated value.
"""
function foldp(f::Function, v0, input::Signal, inputsrest::Signal...;
        typ=typeof(v0), name=auto_name!("foldp", input, inputsrest...))
    n = Signal(typ, v0, (input, inputsrest...); name=name)
    connect_foldp(f, v0, n, (input, inputsrest...))
    n
end

function connect_foldp(f, v0, output, inputs)
    add_action!(output) do
        vals = map(value, inputs)
        set_value!(output, f(output.value, vals...))
    end
end

"""
    sampleon(a, b)

Sample the value of `b` whenever `a` updates.
"""
function sampleon(sample_trigger, input::Signal{T}; name=auto_name!("sampleon", input)) where T
    n = Signal(T, value(input), (sample_trigger,); name=name)
    connect_sampleon(n, input)
    n
end

function connect_sampleon(output, input)
    # this will only get run when sampler updates, as sample_trigger is output's
    # only parent
    add_action!(output) do
        set_value!(output, input.value)
    end
end

"""
    merge(inputs...)

Merge many signals into one. Returns a signal which updates when
any of the inputs update. If many signals update at the same time,
the value of the *youngest* (most recently created) input signal is taken.
"""
function merge(in1::Signal, inputs::Signal...; name=auto_name!("merge", in1, inputs...))
    ins = (in1, inputs...)
    youngestid = maximum(map(x->x.id, ins))
    youngest_val = nodes[youngestid].value
    n = Signal(typejoin(map(eltype, ins)...), value(youngest_val), ins; name=name)
    connect_merge(n, in1, inputs...)
    n
end

function connect_merge(output, inputs...)
    function merge_action()
        lastactive = getlastactive(output)
        lastactive != nothing && set_value!(output, value(lastactive))
        # we don't deactivate! on lastactive == nothing, since I suppose the push
        # should propagate even if some of the nodes died just after updating.
    end
    add_action!(merge_action, output)
end

"""
`getlastactive(merge_node)`
Search backwards in nodes, and return the first active node that is one
of merge_node's parents
"""
function getlastactive(merge_node)
    i = merge_node.id - 1
    while i > 0
        node = nodes[i].value
        if isactive(node) && node in merge_node.parents
            return node
        end
        i -= 1
    end
    # If parent nodes have all been GC'd, but there is still a reference to the
    # merge in user code, then none of the parents should have been active,
    # so the merge action shouldn't run, so we shouldn't have got here. However,
    # in the rare case that the node got GC'd after it was found to be an active
    # parent of the merge but before we got here, then I guess the merge node
    # shouldn't change value.
    return nothing
end

"""
    previous(input, default=value(input))

Create a signal which holds the previous value of `input`.
You can optionally specify a different initial value.
"""
function previous(input::Signal{T}, default=value(input); name=auto_name!("previous", input)) where T
    n = Signal(T, default, (input,); name=name)
    connect_previous(n, input)
    n
end

function connect_previous(output, input)
    prev_value = value(input)
    add_action!(output) do
        set_value!(output, prev_value)
        prev_value = value(input)
    end
end

"""
    delay(input, default=value(input))

Schedule an update to happen after the current update propagates
throughout the signal graph.

Returns the delayed signal.
"""
function delay(input::Signal{T}, default=value(input); name=auto_name!("delay", input)) where T
    n = Signal(T, default, (input,); name=name)
    connect_delay(n, input)
    n
end

function connect_delay(output, input)
    function push_delayed(inpval)
        # only push when input is active (avoids it pushing to itself endlessly)
        push!(output, inpval)
        nothing
    end
    foreach(push_delayed, input; init=nothing)
end

"""
    droprepeats(input)

Drop updates to `input` whenever the new value is the same
as the previous value of the signal.
"""
function droprepeats(input::Signal{T}; name=auto_name!("droprepeats", input)) where T
    n = Signal(T, value(input), (input,); name=name)
    connect_droprepeats(n, input)
    n
end

function connect_droprepeats(output, input)
    prev_value = value(input)
    add_action!(output) do
        if prev_value != value(input)
            set_value!(output, value(input))
            prev_value = value(input)
        else
            deactivate!(output)
        end
    end
end

"""
    flatten(input::Signal{Signal}; typ=Any)

Flatten a signal of signals into a signal which holds the
value of the current signal. The `typ` keyword argument specifies
the type of the flattened signal. It is `Any` by default.
"""
function flatten(input::Signal; typ=Any, name=auto_name!("flatten", input))
    n = Signal(typ, input.value.value, (input,); name=name)
    connect_flatten(n, input)
    n
end

"""
`connect_flatten(output, input)`

`output` is the flatten node, `input` is the Signal{Signal} ("sigsig") node. The
flatten needs to update on changes to the input sigsig, or changes to the value
of the current sig (`current_node`). The former is achieved through a foreach `wire_flatten`
attached to the input sigsig. The latter is achieved through binding the flatten
to `current_node`.
"""
function connect_flatten(output, input)
    # input is a Signal{Signal} (aka sigsig), current_node is the signal/node
    # that is the input's current value. wire_flatten will run when the sigsig gets a new signal as its
    # value. This ensures that set_flatten_val will be run (and flatten output
    # node's value will update) when either the current_node updates, or when
    # the input sigsig updates.
    current_node = input.value
    wire_flatten() = begin
        # If the sigsig's value has changed update output's parents so it will
        # only update when the new current_node updates, and no longer
        # update when the previous signal updates.
        if current_node != input.value
            unbind!(output, current_node, false)
            current_node = input.value
            bind!(output, current_node, false)
            # the bind will have run downstream actions - avoid doubling up
            deactivate!(output)
        end
    end

    add_action!(wire_flatten, output)
    bind!(output, current_node, false)
end

const _bindings = Dict() # XXX GC Issue? can't use WeakKeyDict with Pairs...
const _active_binds = Dict()

"""
    `bind!(dest, src, twoway=true; initial=true)`

for every update to `src` also update `dest` with the same value and, if
`twoway` is true, vice-versa. If `initial` is false, `dest` will only be updated
to `src`'s value when `src` next updates, otherwise (if `initial` is true) both
`dest` and `src` will take `src`'s value immediately.
"""
function bind!(dest::Signal, src::Signal, twoway=true; initial=true)
    dest2src = twoway ? identity : nothing
    bindmap!(dest, identity, src, dest2src)
end

"""
    `bindmap!(dest::Signal, src2dest::Function, src::Signal, dest2src=nothing; initial=true)`

for every update to `src` also update `dest` with a modified value (using the
function `src2dest`) and, if `dest2src` is specified, a two-way update will hold.
If `initial` is false, `dest` will only be updated to `src`'s modified value
when `src` next updates, otherwise (if `initial` is true) both `dest` and `src`
will take their respective modified values immediately.
"""
function bindmap!(dest::Signal, src2dest::Function, src::Signal, dest2src = nothing; initial = true)
    twoway = dest2src ≠ nothing
    if haskey(_bindings, src=>dest)
        # subsequent bind!(dest, src) after initial should be a no-op
        # though we should allow a change in preference for twoway bind.
        if twoway
            bindmap!(src, dest2src, dest, initial = initial)
        end
        return
    end

    # We don't set src as a parent of dest, since a
    # two-way bind would technically introduce a cycle into the signal graph,
    # and I suppose we'd prefer not to have that. Instead we just set dest as
    # active when src updates, which will allow its downstream actions to run.

    ordered_pair = src.id < dest.id ? src=>dest : dest=>src # ordered by id
    twoway && (_active_binds[ordered_pair] = false)
    # the binder action comes after dest, so dest's downstream actions
    # won't run unless we arrange it.
    function bind_updater(srcval)
        if !haskey(_bindings, src=>dest)
            # will happen if has been unbound but node not gc'd
            return
        end
        is_twoway = haskey(_active_binds, ordered_pair)
        if is_twoway && _active_binds[ordered_pair]
            # The _active_binds flag stops the (infinite) cycle of src
            # updating dest updating src ... in the case of a two-way bind
            _active_binds[ordered_pair] = false
        else
            is_twoway && (_active_binds[ordered_pair] = true)
            # we "pause" the current push!, simulate a push! to dest with
            # run_push then resume processing the original push by reactivating
            # the previously active nodes.
            active_nodes = pause_push()
            # `true` below is for dont_remove_dead nodes - messes with active_nodes
            # TODO - check that - not sure it actually does, this may be a relic
            # of an earlier implementation which used the node's id's
            run_push(dest, src2dest(src.value), onerror_rethrow, true) # here is the only place where we implement the modifying function!
            foreach(activate!, active_nodes)
        end
        nothing
    end
    finalizer(src) do src
        unbind!(dest, src, twoway)
    end
    _bindings[src=>dest] = map(bind_updater, src; name="binder: $(src.name)=>$(dest.name)")
    initial && bind_updater(src.value) # init now that _bindings[src=>dest] is set

    if twoway
        bindmap!(src, dest2src, dest, initial = initial)
    end

end

"""
    `unbind!(dest, src, twoway=true)`

remove a link set up using `bind!`
"""
function unbind!(dest::Signal, src::Signal, twoway=true)
    if !haskey(_bindings, src=>dest)
        return
    end

    _bindings[src=>dest] != nothing && close(_bindings[src=>dest])
    delete!(_bindings, src=>dest)

    ordered_pair = src.id < dest.id ? src=>dest : dest=>src # ordered by id
    haskey(_active_binds, ordered_pair) && delete!(_active_binds, ordered_pair)

    if twoway
        unbind!(src, dest, false)
    end
end

"""
Pause a push by recording the active nodes and setting them to inactive.
The push can be resumed by reactivating the nodes.
"""
function pause_push()
    active_nodes = WeakRef[]
    for noderef in nodes
        node = noderef.value
        if isactive(node)
            push!(active_nodes, WeakRef(node))
            deactivate!(node)
        end
    end
    active_nodes
end

"""
`bound_dests(src::Signal)` returns a vector of all signals that will update when
`src` updates, that were bound using `bind!(dest, src)`
"""
bound_dests(s::Signal) = [dest for (src, dest) in keys(_bindings) if src == s]

"""
`bound_srcs(dest::Signal)` returns a vector of all signals that will cause
an update to `dest` when they update, that were bound using `bind!(dest, src)`
"""
bound_srcs(s::Signal) = [src for (src, dest) in keys(_bindings) if dest == s]


text:
using BinDeps
@BinDeps.setup

deps = [ libxgboostwrapper = library_dependency("xgboostwrapper", aliases = ["libxgboostwrapper.so"]) ]

prefix=joinpath(BinDeps.depsdir(libxgboostwrapper),"usr")
provides(BuildProcess,
           (@build_steps begin
               `rm -rf xgboost`
               `git clone https://github.com/dmlc/xgboost.git --recursive`
               CreateDirectory(prefix)
               CreateDirectory(joinpath(prefix, "lib"))
               @build_steps begin
                   ChangeDirectory("xgboost")
                   FileRule(joinpath(prefix,"lib","libxgboostwrapper.so"), @build_steps begin
                       `git checkout 0.47` # v0.40
                       `bash build.sh`
                       `cp wrapper/libxgboostwrapper.so $prefix/lib`
                   end)
               end
            end),
         libxgboostwrapper)

@BinDeps.install Dict(:xgboostwrapper => :_xgboost)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# fold(l|r) & mapfold(l|r)
@test foldl(-, 1:5) == -13
@test foldl(-, 10, 1:5) == -5

@test Base.mapfoldl(abs2, -, 2:5) == -46
@test Base.mapfoldl(abs2, -, 10, 2:5) == -44

@test Base.mapfoldl(abs2, /, 2:5) ≈ 1/900
@test Base.mapfoldl(abs2, /, 10, 2:5) ≈ 1/1440

@test Base.mapfoldl((x)-> x ⊻ true, &, true, [true false true false false]) == false
@test Base.mapfoldl((x)-> x ⊻ true, &, [true false true false false]) == false

@test Base.mapfoldl((x)-> x ⊻ true, |, [true false true false false]) == true
@test Base.mapfoldl((x)-> x ⊻ true, |, false, [true false true false false]) == true

@test foldr(-, 1:5) == 3
@test foldr(-, 10, 1:5) == -7

@test Base.mapfoldr(abs2, -, 2:5) == -14
@test Base.mapfoldr(abs2, -, 10, 2:5) == -4

# reduce & mapreduce
@test reduce((x,y)->"($x+$y)", 9:11) == "((9+10)+11)"
@test reduce(max, [8 6 7 5 3 0 9]) == 9
@test reduce(+, 1000, 1:5) == (1000 + 1 + 2 + 3 + 4 + 5)
@test reduce(+,1) == 1

@test mapreduce(-, +, [-10 -9 -3]) == ((10 + 9) + 3)
@test mapreduce((x)->x[1:3], (x,y)->"($x+$y)", ["abcd", "efgh", "01234"]) == "((abc+efg)+012)"

# sum

@test sum(Int8[]) === Int32(0)
@test sum(Int[]) === Int(0)
@test sum(Float64[]) === 0.0

@test sum(Int8(3)) === Int8(3)
@test sum(3) === 3
@test sum(3.0) === 3.0

@test sum([Int8(3)]) === Int32(3)
@test sum([3]) === 3
@test sum([3.0]) === 3.0

z = reshape(1:16, (2,2,2,2))
fz = float(z)
@test sum(z) === 136
@test sum(fz) === 136.0

@test_throws ArgumentError sum(sin, Int[])
@test sum(sin, 3) == sin(3.0)
@test sum(sin, [3]) == sin(3.0)
a = sum(sin, z)
@test a ≈ sum(sin, fz)
@test a ≈ sum(sin.(fz))

z = [-4, -3, 2, 5]
fz = float(z)
a = randn(32) # need >16 elements to trigger BLAS code path
b = complex(randn(32), randn(32))
@test sumabs(Float64[]) === 0.0
@test sumabs([Int8(-2)]) === Int32(2)
@test sumabs(z) === 14
@test sumabs(fz) === 14.0
@test sumabs(a) ≈ sum(abs.(a))
@test sumabs(b) ≈ sum(abs.(b))

@test sumabs2(Float64[]) === 0.0
@test sumabs2([Int8(-2)]) === Int32(4)
@test sumabs2(z) === 54
@test sumabs2(fz) === 54.0
@test sumabs2(a) ≈ sum(abs2.(a))
@test sumabs2(b) ≈ sum(abs2.(b))

# check variants of summation for type-stability and other issues (#6069)
sum2(itr) = invoke(sum, Tuple{Any}, itr)
plus(x,y) = x + y
sum3(A) = reduce(plus, A)
sum4(itr) = invoke(reduce, Tuple{Function, Any}, plus, itr)
sum5(A) = reduce(plus, 0, A)
sum6(itr) = invoke(reduce, Tuple{Function, Int, Any}, plus, 0, itr)
sum7(A) = mapreduce(x->x, plus, A)
sum8(itr) = invoke(mapreduce, Tuple{Function, Function, Any}, x->x, plus, itr)
sum9(A) = mapreduce(x->x, plus, 0, A)
sum10(itr) = invoke(mapreduce, Tuple{Function, Function, Int, Any}, x->x,plus,0,itr)
for f in (sum2, sum5, sum6, sum9, sum10)
    @test sum(z) == f(z)
    @test sum(Int[]) == f(Int[]) == 0
    @test sum(Int[7]) == f(Int[7]) == 7
    @test typeof(f(Int8[])) == typeof(f(Int8[1])) == typeof(f(Int8[1 7]))
end
for f in (sum3, sum4, sum7, sum8)
    @test sum(z) == f(z)
    @test_throws ArgumentError f(Int[])
    @test sum(Int[7]) == f(Int[7]) == 7
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

@test sum_kbn([1,1e100,1,-1e100]) == 2
@test sum_kbn(Float64[]) == 0.0

# prod

@test prod(Int[]) === 1
@test prod(Int8[]) === Int32(1)
@test prod(Float64[]) === 1.0

@test prod([3]) === 3
@test prod([Int8(3)]) === Int32(3)
@test prod([3.0]) === 3.0

@test prod(z) === 120
@test prod(fz) === 120.0

@test prod(1:big(16)) == big(20922789888000)
@test prod(big(typemax(Int64)):big(typemax(Int64))+16) == parse(BigInt,"25300281663413827620486300433089141956148633919452440329174083959168114253708467653081909888307573358090001734956158476311046124934597861626299416732205795533726326734482449215730132757595422510465791525610410023802664753402501982524443370512346073948799084936298007821432734720004795146875180123558814648586972474376192000")

@test typeof(prod(Array(trues(10)))) == Bool

# check type-stability
prod2(itr) = invoke(prod, Tuple{Any}, itr)
@test prod(Int[]) === prod2(Int[]) === 1
@test prod(Int[7]) === prod2(Int[7]) === 7
@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1, 7])) == Int32
@test typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7])) == Int32

# maximum & minimum & extrema

@test_throws ArgumentError maximum(Int[])
@test_throws ArgumentError minimum(Int[])

@test maximum(5) == 5
@test minimum(5) == 5
@test extrema(5) == (5, 5)

@test maximum([4, 3, 5, 2]) == 5
@test minimum([4, 3, 5, 2]) == 2
@test extrema([4, 3, 5, 2]) == (2, 5)

@test isnan(maximum([NaN]))
@test isnan(minimum([NaN]))
@test isequal(extrema([NaN]), (NaN, NaN))

@test isnan(maximum([NaN, 2., 3.]))
@test isnan(minimum([NaN, 2., 3.]))
@test isequal(extrema([NaN, 2., 3.]), (NaN,NaN))

@test isnan(maximum([4., 3., NaN, 5., 2.]))
@test isnan(minimum([4., 3., NaN, 5., 2.]))
@test isequal(extrema([4., 3., NaN, 5., 2.]), (NaN,NaN))

@test maxabs(Int[]) == 0
@test_throws ArgumentError Base.minabs(Int[])

@test maxabs(-2) == 2
@test minabs(-2) == 2
@test maxabs([1, -2, 3, -4]) == 4
@test minabs([-1, 2, -3, 4]) == 1

@test maximum(abs2, 3:7) == 49
@test minimum(abs2, 3:7) == 9

@test maximum(Int16[1]) === Int16(1)
@test maximum(collect(Int16(1):Int16(100))) === Int16(100)
@test maximum(Int32[1,2]) === Int32(2)

@test extrema(reshape(1:24,2,3,4),1) == reshape([(1,2),(3,4),(5,6),(7,8),(9,10),(11,12),(13,14),(15,16),(17,18),(19,20),(21,22),(23,24)],1,3,4)
@test extrema(reshape(1:24,2,3,4),2) == reshape([(1,5),(2,6),(7,11),(8,12),(13,17),(14,18),(19,23),(20,24)],2,1,4)
@test extrema(reshape(1:24,2,3,4),3) == reshape([(1,19),(2,20),(3,21),(4,22),(5,23),(6,24)],2,3,1)

# any & all

@test any([]) == false
@test any(Bool[]) == false
@test any([true]) == true
@test any([false, false]) == false
@test any([false, true]) == true
@test any([true, false]) == true
@test any([true, true]) == true
@test any([true, true, true]) == true
@test any([true, false, true]) == true
@test any([false, false, false]) == false

@test all([]) == true
@test all(Bool[]) == true
@test all([true]) == true
@test all([false, false]) == false
@test all([false, true]) == false
@test all([true, false]) == false
@test all([true, true]) == true
@test all([true, true, true]) == true
@test all([true, false, true]) == false
@test all([false, false, false]) == false

@test any(x->x>0, []) == false
@test any(x->x>0, Int[]) == false
@test any(x->x>0, [-3]) == false
@test any(x->x>0, [4]) == true
@test any(x->x>0, [-3, 4, 5]) == true

@test all(x->x>0, []) == true
@test all(x->x>0, Int[]) == true
@test all(x->x>0, [-3]) == false
@test all(x->x>0, [4]) == true
@test all(x->x>0, [-3, 4, 5]) == false

@test reduce(|, fill(trues(5), 24))  == trues(5)
@test reduce(|, fill(falses(5), 24)) == falses(5)
@test reduce(&, fill(trues(5), 24))  == trues(5)
@test reduce(&, fill(falses(5), 24)) == falses(5)

@test_throws TypeError any(x->0, [false])
@test_throws TypeError all(x->0, [false])

# short-circuiting any and all

let c = [0, 0], A = 1:1000
    any(x->(c[1]=x; x==10), A)
    all(x->(c[2]=x; x!=10), A)

    @test c == [10,10]
end

# any and all with functors

immutable SomeFunctor end
(::SomeFunctor)(x) = true

@test any(SomeFunctor(), 1:10)
@test all(SomeFunctor(), 1:10)

# in

@test in(1, Int[]) == false
@test in(1, Int[1]) == true
@test in(1, Int[2]) == false
@test in(0, 1:3) == false
@test in(1, 1:3) == true
@test in(2, 1:3) == true

# contains

@test contains("quick fox", "fox") == true
@test contains("quick fox", "lazy dog") == false

# count & countnz

@test count(x->x>0, Int[]) == 0
@test count(x->x>0, -3:5) == 5

@test countnz(Int[]) == 0
@test countnz(Int[0]) == 0
@test countnz(Int[1]) == 1
@test countnz([1, 0, 2, 0, 3, 0, 4]) == 4

## cumsum, cummin, cummax

z = rand(10^6)
let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
    @test (es - sum(z)) < es * 1e-13
    cs = cumsum(z)
    @test (es - cs[end]) < es * 1e-13
    @test (es2 - cs[10^5]) < es2 * 1e-13
end

@test sum(collect(map(UInt8,0:255))) == 32640
@test sum(collect(map(UInt8,254:255))) == 509

A = reshape(map(UInt8, 101:109), (3,3))
@test @inferred(sum(A)) == 945
@test @inferred(sum(view(A, 1:3, 1:3))) == 945

A = reshape(map(UInt8, 1:100), (10,10))
@test @inferred(sum(A)) == 5050
@test @inferred(sum(view(A, 1:10, 1:10))) == 5050

# issue #11618
@test sum([-0.0]) === -0.0
@test sum([-0.0, -0.0]) === -0.0
@test prod([-0.0, -0.0]) === 0.0

#contains
let A = collect(1:10)
    @test A ∋ 5
    @test A ∌ 11
    @test contains(==,A,6)
end

# issue #18695
test18695(r) = sum( t^2 for t in r )
@test @inferred(test18695([1.0,2.0,3.0,4.0])) == 30.0
@test_throws ArgumentError test18695(Any[])

max_stars_repo_name: MichaelHatherly/julia

max_issues_repo_name: MichaelHatherly/julia

max_forks_repo_name: MichaelHatherly/julia


text:
"""
This example learns to classify hand-written digits from the MNIST
dataset.  There are 60000 training and 10000 test examples. Each input
x consists of 784 pixels representing a 28x28 image.  The pixel values
have been normalized to [0,1]. Each output y is a ten-dimensional
one-hot vector (a vector that has a single non-zero component)
indicating the correct class (0-9) for a given image.

To run the demo, simply `include("mnist.jl")` and run `MNIST.train()`.
The dataset will be automatically downloaded.  You can provide the
initial weights as an optional argument to `train`, which should have
the form [w0,b0,w1,b1,...] where wi (with size = output x input) is
the weight matrix and bi (with size = output) is the bias vector for
layer i.  The function `MNIST.weights(h...)` can be used to create
random starting weights for a neural network with hidden sizes (h...).
If not specified, default weights are created using `MNIST.weights()`
which correspond to a 0 hidden layer network, i.e. a softmax model.
`train` also accepts the following keyword arguments: `lr` specifies
the learning rate, `epochs` gives number of epochs.  The cross entropy
loss and accuracy for the train and test sets will be printed at every
epoch and optimized parameters will be returned.
"""
module MNIST
using AutoGrad
using GZip
using Main
using Compat

function predict(w, x)
    i = 1
    while i+2 < length(w)
        x = max(0, w[i]*x .+ w[i+1])
        i += 2
    end
    return w[i]*x .+ w[i+1]
end

function loss(w, x, ygold)
    ypred = predict(w, x)
    ynorm = ypred .- log(sum(exp(ypred),1))
    -sum(ygold .* ynorm) / size(ygold,2)
end

function accuracy(w, x, ygold)
    ypred = predict(w, x)
    sum((ypred .== maximum(ypred,1)) & (ygold .== maximum(ygold,1))) / size(ygold,2)
end

function train(w=weights(); lr=.1, epochs=20)
    isdefined(MNIST,:dtrn) || loaddata()
    println((0, :ltrn, loss(w,xtrn,ytrn), :ltst, loss(w,xtst,ytst), :atrn, accuracy(w,xtrn,ytrn), :atst, accuracy(w,xtst,ytst)))
    gradfun = grad(loss)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = gradfun(w, x, y)
            for i in 1:length(w)
                w[i] -= lr * g[i]
            end
        end
        println((epoch, :ltrn, loss(w,xtrn,ytrn), :ltst, loss(w,xtst,ytst), :atrn, accuracy(w,xtrn,ytrn), :atst, accuracy(w,xtst,ytst)))
    end
    return w
end

function weights(h...; seed=nothing)
    seed==nothing || srand(seed)
    w = Any[]
    x = 28*28
    for y in [h..., 10]
        push!(w, convert(Array{Float32}, 0.1*randn(y,x)))
        push!(w, zeros(Float32, y))
        x = y
    end
    return w
end

function loaddata()
    info("Loading data...")
    global xtrn, xtst, ytrn, ytst, dtrn
    xshape(a)=reshape(a./255f0,784,div(length(a),784))
    yshape(a)=(a[a.==0]=10; full(sparse(convert(Vector{Int},a),1:length(a),1f0)))
    xtrn = xshape(gzread("train-images-idx3-ubyte.gz")[17:end])
    xtst = xshape(gzread("t10k-images-idx3-ubyte.gz")[17:end])
    ytrn = yshape(gzread("train-labels-idx1-ubyte.gz")[9:end])
    ytst = yshape(gzread("t10k-labels-idx1-ubyte.gz")[9:end])
    dtrn = minibatch(xtrn, ytrn, 100)
    info("Loading done...")
end

function gzread(file; path=joinpath(AutoGrad.datapath,file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = @compat read(f)
    close(f)
    return(a)
end

function minibatch(x, y, batchsize)
    data = Any[]
    nx = size(x,2)
    for i=1:batchsize:nx
        j=min(i+batchsize-1,nx)
        push!(data, (x[:,i:j], y[:,i:j]))
    end
    return data
end

end # module


text:
module TestCPortGraphs
using Test

using Catlab.Graphs, Catlab.WiringDiagrams, Catlab.CategoricalAlgebra
import Catlab.CategoricalAlgebra: migrate!

g = cycle_graph(Graph, 4)
pg = CPortGraph()
migrate!(pg, g)

@test subpart(pg, :box) == 1:4
@test subpart(pg, :src) == 1:4
@test subpart(pg, :tgt) == [2,3,4,1]

bar = Graph()
add_vertices!(bar, 2)
add_edge!(bar, 1,2)
Bar = migrate!(OpenCPortGraph(), bar)
@test subpart(Bar, :box) == 1:2
@test subpart(Bar, :con) == 1:2

opg = OpenCPortGraph(pg)
add_parts!(opg, :OuterPort, 2; con=[1,2])

g = OpenCPortGraph()
add_parts!(g, :Box, 2)
add_parts!(g, :Port, 4; box=[1,1,2,2])
add_parts!(g, :Wire, 2; src=[1,2], tgt=[3,4])
@test subpart(g, :src) == 1:2
@test subpart(g, :tgt) == 3:4
@test subpart(g, :box) == [1,1,2,2]
@test subpart(g, :con) == Int64[]

h = ocompose(g, [Bar, Bar])
@test subpart(h, :src) == [1,3,1,2]
@test subpart(h, :tgt) == [2,4,3,4]

barbell(k::Int) = begin
  g = OpenCPortGraph()
  add_parts!(g, :Box, 2)
  add_parts!(g, :Port, 2k; box=[fill(1,k); fill(2,k)])
  add_parts!(g, :Wire, k; src=1:k, tgt=k+1:2k)
  return g
end

@testset "Hexagon" begin
  g = OpenCPortGraph()
  add_parts!(g, :Box, 3)
  add_parts!(g, :Port, 6; box=[1,1,2,2,3,3])
  add_parts!(g, :Wire, 3; src=[2,4,6], tgt=[3,5,1])
  hex = ocompose(g, [Bar, Bar, Bar])
  @test subpart(hex, :src) == [1,3,5,2,4,6]
  @test subpart(hex, :tgt) == [2,4,6,3,5,1]
  @test subpart(hex, :con) == Int[]
  add_parts!(hex, :OuterPort, 6, con=1:6)

  twohex = ocompose(barbell(6), [hex, hex])
  g′ = migrate!(Graph(), migrate!(CPortGraph(), twohex))
  @test subpart(g′, :src) == [1,3,5,2,4,6,7,9,11,8,10,12,1,2,3,4,5,6]
  @test subpart(g′, :tgt) == [2,4,6,3,5,1,8,10,12,9,11,7,7,8,9,10,11,12]
end

@testset "Bundling CPGs" begin

  EdgeGraph() = begin
    g = Graph()
    add_parts!(g, :V, 2)
    add_parts!(g, :E, 1, src=[1], tgt=[2])
    return g
  end

  e₁ = EdgeGraph()
  oe₁ = migrate!(OpenCPortGraph(), (migrate!(CPortGraph(), e₁)))
  @test subpart(oe₁, :con) == 1:2
  mboe₁ = migrate!(BundledCPG(), oe₁)
  @test subpart(mboe₁, :bun) == 1:2

  boe₁ = BundledCPG(oe₁)
  add_parts!(boe₁, :Bundle, 1)
  set_subpart!(boe₁, :bun, [1,1])
  boe₁

  f = OpenCPortGraph(migrate!(CPortGraph(), e₁))
  f = BundledCPG(f)
  fee = ocompose(f, [boe₁, boe₁])
  @test fee[:, :src] == [1,3,1,2]
  @test fee[:, :tgt] == [2,4,3,4]

  add_parts!(f, :Bundle, 2)
  add_parts!(f, :OuterPort, 2, con=[1,2], bun=[1,2])
  fee′ = ocompose(f, [boe₁, boe₁])
  @test fee′[:, :con] == 1:4
  @test fee′[:, :bun] == [1,1,2,2]
  h = ocompose(f, [fee′, fee′])
  @test nparts(h, :Wire) == 10

  set_subpart!(fee′, :bun, 1)
  h = ocompose(f, [fee′, fee′])
  @test nparts(h, :Wire) == 12
end

end


text:
int_rules_1_2_3_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.3.2*(d*x)^m*(a+b*x^n+c*x^(2*n))^p =#
    #= Antiderivative(((~d')*(~x))^(~m')*((~b')*(~x)^(~n)+(~c')*(~x)^(~n2'))^(~p'),~x) :=  d^m*Antiderivative(x^(m+n*p)*(b+c*x^n)^p,x) <-- FreeQ([b,c,d,m,n],x) &&  EqQ(n2,2*n) && IntegerQ(p) && (IntegerQ(m) || GtQ(d,0)) =#
    #= Antiderivative(((~d')*(~x))^(~m')*((~b')*(~x)^(~n)+(~c')*(~x)^(~n2'))^(~p'),~x) :=  1/d^(n*p)*Antiderivative((d*x)^(m+n*p)*(b+c*x^n)^p,x) <-- FreeQ([b,c,d,m],x) &&  EqQ(n2,2*n) && IntegerQ(p) && IntegerQ(n) =#
    #= Antiderivative(((~d')*(~x))^(~m')*((~b')*(~x)^(~n)+(~c')*(~x)^(~n2'))^(~p'),~x) :=  (d*x)^m/x^m*Antiderivative(x^(m+n*p)*(b+c*x^n)^p,x) <-- FreeQ([b,c,d,m,n],x) &&  EqQ(n2,2*n) && IntegerQ(p) && Not(IntegerQ(m) || GtQ(d,0)) =#
    #= Antiderivative(((~d')*(~x))^(~m')*((~b')*(~x)^(~n)+(~c')*(~x)^(~n2'))^(~p),~x) :=  (b*x^n+c*x^(2*n))^p/((d*x)^(n*p)*(b+c*x^n)^p)*Antiderivative((d*x)^(m+n*p)*(b+c* x^2)^p,x) <-- FreeQ([b,c,d,m,n,p],x) && EqQ(n2,2*n) &&  Not(IntegerQ(p)) =#
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ(Simplify((~m - ~n) + 1), 0))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (IGtQ(~p, 0) && Not(IntegerQ(Simplify((~m + 1) / ~n)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => Antiderivative((~x) ^ (~m + 2 * ~n * ~p) * (~c + ~b * (~x) ^ -(~n) + ~a * (~x) ^ (-2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (ILtQ(~p, 0) && NegQ(~n)))
    #= Antiderivative(((~d')*(~x))^(~m')*((~a)+(~b')*(~x)^(~n')+(~c')*(~x)^(~n2'))^(~p),~x) :=  1/c^p*Antiderivative((d*x)^m*(b/2+c*x^n)^(2*p),x) <-- FreeQ([a,b,c,d,m,n,p],x) &&  EqQ(n2,2*n) && EqQ(b^2-4*a*c,0) && IntegerQ(p) =#
    #= Antiderivative(((~d')*(~x))^(~m')*((~a)+(~b')*(~x)^(~n)+(~c')*(~x)^(~n2'))^(~p),~x) :=  (d*x)^(m+1)*(a+b*x^n+c*x^(2*n))^(p+1)/(2*a*d*n*(p+1)*(2*p+1)) - (d*x)^(m+1)*(2*a+b*x^n)*(a+b*x^n+c*x^(2*n))^p/(2*a*d*n*(2*p+1)) <--  FreeQ([a,b,c,d,m,n,p],x) && EqQ(n2,2*n) && EqQ(b^2-4*a*c,0) &&  Not(IntegerQ(p)) && EqQ(m+2*n*(p+1)+1,0) && NeQ(2*p+1,0) =#
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~c) ^ IntPart(~p) * (~b / 2 + ~c * (~x) ^ ~n) ^ (2 * FracPart(~p)))) * Antiderivative((~d * ~x) ^ ~m * (~b / 2 + ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p - 1 / 2)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) / (1 + (2 * ~c * (~x) ^ ~n) / ~b) ^ (2 * FracPart(~p))) * Antiderivative((~d * ~x) ^ ~m * (1 + (2 * ~c * (~x) ^ ~n) / ~b) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(2 * ~p))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(Simplify((~m + 1) / ~n))))
    @apply_utils Antiderivative((~d * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (((~d) ^ IntPart(~m) * (~d * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(Simplify((~m + 1) / ~n))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([k = GCD(~m + 1, ~n)], (1 / k) * Subst(Antiderivative((~x) ^ ((~m + 1) / k - 1) * (~a + ~b * (~x) ^ (~n / k) + ~c * (~x) ^ ((2 * ~n) / k)) ^ ~p, ~x), ~x, (~x) ^ k) <-- k != 1) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && IntegerQ(~m))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([k = Denominator(~m)], (k / ~d) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + (~b * (~x) ^ (k * ~n)) / (~d) ^ ~n + (~c * (~x) ^ (2 * k * ~n)) / (~d) ^ (2 * ~n)) ^ ~p, ~x), ~x, (~d * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (FractionQ(~m) && IntegerQ(~p)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d) ^ (~n - 1) * (~d * ~x) ^ ((~m - ~n) + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p * (~b * ~n * ~p + ~c * (~m + ~n * (2 * ~p - 1) + 1) * (~x) ^ ~n)) / (~c * (~m + 2 * ~n * ~p + 1) * (~m + ~n * (2 * ~p - 1) + 1)) - ((~n * ~p * (~d) ^ ~n) / (~c * (~m + 2 * ~n * ~p + 1) * (~m + ~n * (2 * ~p - 1) + 1))) * Antiderivative((~d * ~x) ^ (~m - ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p - 1) * Simp(~a * ~b * ((~m - ~n) + 1) - (2 * ~a * ~c * (~m + ~n * (2 * ~p - 1) + 1) - (~b) ^ 2 * (~m + ~n * (~p - 1) + 1)) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && (GtQ(~m, ~n - 1) && (NeQ(~m + 2 * ~n * ~p + 1, 0) && NeQ(~m + ~n * (2 * ~p - 1) + 1, 0)))))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / (~d * (~m + 1)) - ((~n * ~p) / ((~d) ^ ~n * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + ~n) * (~b + 2 * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && LtQ(~m, -1)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / (~d * (~m + 2 * ~n * ~p + 1)) + ((~n * ~p) / (~m + 2 * ~n * ~p + 1)) * Antiderivative((~d * ~x) ^ ~m * (2 * ~a + ~b * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (IGtQ(~p, 0) && NeQ(~m + 2 * ~n * ~p + 1, 0)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d) ^ (~n - 1) * (~d * ~x) ^ ((~m - ~n) + 1) * (~b + 2 * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) - ((~d) ^ ~n / (~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d * ~x) ^ (~m - ~n) * (~b * ((~m - ~n) + 1) + 2 * ~c * (~m + 2 * ~n * (~p + 1) + 1) * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (ILtQ(~p, -1) && (GtQ(~m, ~n - 1) && LeQ(~m, 2 * ~n - 1))))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (-((~d) ^ (2 * ~n - 1)) * (~d * ~x) ^ ((~m - 2 * ~n) + 1) * (2 * ~a + ~b * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + ((~d) ^ (2 * ~n) / (~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d * ~x) ^ (~m - 2 * ~n) * (2 * ~a * ((~m - 2 * ~n) + 1) + ~b * (~m + ~n * (2 * ~p + 1) + 1) * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (ILtQ(~p, -1) && GtQ(~m, 2 * ~n - 1)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (-((~d * ~x) ^ (~m + 1)) * (((~b) ^ 2 - 2 * ~a * ~c) + ~b * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~d * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1) * Simp(((~b) ^ 2 * (~m + ~n * (~p + 1) + 1) - 2 * ~a * ~c * (~m + 2 * ~n * (~p + 1) + 1)) + ~b * ~c * (~m + ~n * (2 * ~p + 3) + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && ILtQ(~p, -1))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d) ^ (2 * ~n - 1) * (~d * ~x) ^ ((~m - 2 * ~n) + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~c * (~m + 2 * ~n * ~p + 1)) - ((~d) ^ (2 * ~n) / (~c * (~m + 2 * ~n * ~p + 1))) * Antiderivative((~d * ~x) ^ (~m - 2 * ~n) * Simp(~a * ((~m - 2 * ~n) + 1) + ~b * (~m + ~n * (~p - 1) + 1) * (~x) ^ ~n, ~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (GtQ(~m, 2 * ~n - 1) && (NeQ(~m + 2 * ~n * ~p + 1, 0) && IntegerQ(~p))))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~d * (~m + 1)) - (1 / (~a * (~d) ^ ~n * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + ~n) * (~b * (~m + ~n * (~p + 1) + 1) + ~c * (~m + 2 * ~n * (~p + 1) + 1) * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (LtQ(~m, -1) && IntegerQ(~p)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => (~d * ~x) ^ (~m + 1) / (~a * ~d * (~m + 1)) - (1 / (~a * (~d) ^ ~n)) * Antiderivative(((~d * ~x) ^ (~m + ~n) * (~b + ~c * (~x) ^ ~n)) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && LtQ(~m, -1))))
    @apply_utils Antiderivative((~x) ^ ~m / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => Antiderivative(PolynomialDivide((~x) ^ ~m, ~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && IGtQ(~m, 3 * ~n - 1))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => ((~d) ^ (2 * ~n - 1) * (~d * ~x) ^ ((~m - 2 * ~n) + 1)) / (~c * ((~m - 2 * ~n) + 1)) - ((~d) ^ (2 * ~n) / ~c) * Antiderivative(((~d * ~x) ^ (~m - 2 * ~n) * (~a + ~b * (~x) ^ ~n)) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && GtQ(~m, 2 * ~n - 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt(~a / ~c, 2)], With([r = Rt(2q - ~b / ~c, 2)], (1 / (2 * ~c * r)) * Antiderivative(((~x) ^ (~m - 3 * (~n / 2)) * (q + r * (~x) ^ (~n / 2))) / (q + r * (~x) ^ (~n / 2) + (~x) ^ ~n), ~x) - (1 / (2 * ~c * r)) * Antiderivative(((~x) ^ (~m - 3 * (~n / 2)) * (q - r * (~x) ^ (~n / 2))) / ((q - r * (~x) ^ (~n / 2)) + (~x) ^ ~n), ~x))) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n / 2, 0) && (IGtQ(~m, 0) && (GeQ(~m, (3 * ~n) / 2) && (LtQ(~m, 2 * ~n) && NegQ((~b) ^ 2 - 4 * ~a * ~c)))))))
    @apply_utils Antiderivative((~x) ^ ~(m') / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt(~a / ~c, 2)], With([r = Rt(2q - ~b / ~c, 2)], (1 / (2 * ~c * r)) * Antiderivative((~x) ^ (~m - ~n / 2) / ((q - r * (~x) ^ (~n / 2)) + (~x) ^ ~n), ~x) - (1 / (2 * ~c * r)) * Antiderivative((~x) ^ (~m - ~n / 2) / (q + r * (~x) ^ (~n / 2) + (~x) ^ ~n), ~x))) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n / 2, 0) && (IGtQ(~m, 0) && (GeQ(~m, ~n / 2) && (LtQ(~m, (3 * ~n) / 2) && NegQ((~b) ^ 2 - 4 * ~a * ~c)))))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((~d) ^ ~n / 2) * (~b / q + 1) * Antiderivative((~d * ~x) ^ (~m - ~n) / (~b / 2 + q / 2 + ~c * (~x) ^ ~n), ~x) - ((~d) ^ ~n / 2) * (~b / q - 1) * Antiderivative((~d * ~x) ^ (~m - ~n) / ((~b / 2 - q / 2) + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && GeQ(~m, ~n))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~c / q) * Antiderivative((~d * ~x) ^ ~m / ((~b / 2 - q / 2) + ~c * (~x) ^ ~n), ~x) - (~c / q) * Antiderivative((~d * ~x) ^ ~m / (~b / 2 + q / 2 + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -(Subst(Antiderivative((~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && IntegerQ(~m))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([k = Denominator(~m)], (-k / ~d) * Subst(Antiderivative((~a + ~b * (~d) ^ -(~n) * (~x) ^ (-k * ~n) + ~c * (~d) ^ (-2 * ~n) * (~x) ^ (-2 * k * ~n)) ^ ~p / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~d * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && FractionQ(~m))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -((~d) ^ IntPart(~m)) * (~d * ~x) ^ FracPart(~m) * ((~x) ^ -1) ^ FracPart(~m) * Subst(Antiderivative((~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && Not(RationalQ(~m)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * (~x) ^ (k * ~n) + ~c * (~x) ^ (2 * k * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && FractionQ(~n)))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~d) ^ IntPart(~m) * (~d * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && FractionQ(~n)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ Simplify(~n / (~m + 1)) + ~c * (~x) ^ Simplify((2 * ~n) / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~d) ^ IntPart(~m) * (~d * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative((~d * ~x) ^ ~m / ((~b - q) + 2 * ~c * (~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative((~d * ~x) ^ ~m / (~b + q + 2 * ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (-((~d * ~x) ^ (~m + 1)) * (((~b) ^ 2 - 2 * ~a * ~c) + ~b * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~d * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1) * Simp(((~b) ^ 2 * (~n * (~p + 1) + ~m + 1) - 2 * ~a * ~c * (~m + 2 * ~n * (~p + 1) + 1)) + ~b * ~c * (2 * ~n * ~p + 3 * ~n + ~m + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && ILtQ(~p + 1, 0)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) / ((1 + (2 * ~c * (~x) ^ ~n) / (~b + Rt((~b) ^ 2 - 4 * ~a * ~c, 2))) ^ FracPart(~p) * (1 + (2 * ~c * (~x) ^ ~n) / (~b - Rt((~b) ^ 2 - 4 * ~a * ~c, 2))) ^ FracPart(~p))) * Antiderivative((~d * ~x) ^ ~m * (1 + (2 * ~c * (~x) ^ ~n) / (~b + sqrt((~b) ^ 2 - 4 * ~a * ~c))) ^ ~p * (1 + (2 * ~c * (~x) ^ ~n) / (~b - sqrt((~b) ^ 2 - 4 * ~a * ~c))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && EqQ(~n2, 2 * ~n)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (~m - ~n * ~p) * (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n], ~x) && (EqQ(~mn, -(~n)) && (IntegerQ(~p) && PosQ(~n)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (((~x) ^ (~n * FracPart(~p)) * (~a + ~b / (~x) ^ ~n + ~c * (~x) ^ ~n) ^ FracPart(~p)) / (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~x) ^ (~m - ~n * ~p) * (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~mn, -(~n)) && (Not(IntegerQ(~p)) && PosQ(~n)))
    @apply_utils Antiderivative((~d * ~x) ^ ~(m') * (~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (((~d) ^ IntPart(~m) * (~d * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && EqQ(~mn, -(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * (~v) ^ ~n + ~(c') * (~v) ^ ~(n2')) ^ ~(p'), ~x) => (1 / Coefficient(~v, ~x, 1) ^ (~m + 1)) * Subst(Antiderivative(SimplifyIntegrand((~x - Coefficient(~v, ~x, 0)) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (LinearQ(~v, ~x) && (IntegerQ(~m) && NeQ(~v, ~x))))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * (~v) ^ ~n + ~(c') * (~v) ^ ~(n2')) ^ ~(p'), ~x) => ((~u) ^ ~m / (Coefficient(~v, ~x, 1) * (~v) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && LinearPairQ(~u, ~v, ~x))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.3 General/1.2.3.2 (d x)^m (a+b x^n+c x^(2 n))^p.jl


text:
"""
    fill_zero(keys_unique::NamedTuple; value_colnames)
    fill_zero(keys_unique::NamedTuple, df::DataFrame)
    fill_zero(df::DataFrame...)
    fill_zero(d::Dict...)
    fill_zero(keys_unique, d::Dict)

This function can be used to fill zeros in either a dictionary or DataFrame.
- Options for dictionary editing:
    - If only (a) dictionary/ies is/are input, the dictionaries will be edited such that
        they all contain all permutations of their key values. All dictionaries in a
        resultant list of dictionaries will be the same length.
    - If a dictionary is input with a list of keys, it will be edited to ensure that it
        includes all permutations.
    - If only a list of keys is input, a new dictionary will be created, containing all key
        permutations with values initialized to zero.
- Options for DataFrame editing:
    - If only (a) DataFrame(s) is/are input, the DataFrame(s) will be edited such that
        they all contain all permutations of their key values. All DataFrames in a
        resultant list of DataFrames will be the same length.
    - If a DataFrame is input with a NamedTuple, it will be edited to ensure that it
        includes all permutations of the NamedTuple's values.
    - If only a NamedTuple is input, a new DataFrame will be created, containing all key
        permutations with values initialized to zero.

# Arguments
- `keys_unique::Tuple`: A list of arrays whose permutations should be included in the
    resultant dictionary.
- `keys_unique::NamedTuple`: A list of arrays whose permutations should be included in the
    resultant dictionary. The NamedTuple's keys correspond to the DataFrame columns where
    they will be stored.
- `d::Dict...`: The dictionary/ies to edit.
- `df::DataFrame...`: The DataFrame(s) to edit.

# Keywords
- `value_colnames::Any = :value`: "value" column labels to add and set to zero when creating
    a new DataFrame. Default is `:value`.

# Returns
- `d::Dict...` if input included dictionaries and/or Tuples
- `df::DataFrame...` if input included DataFrames and/or NamedTuples
"""
function fill_zero(keys_fill::NamedTuple; value_colnames = :value)
    df_fill = DataFrame(permute(keys_fill))
    return edit_with(df_fill, Add.(convert_type.(Symbol, value_colnames), 0.))
end

function fill_zero(keys_fill::Any; permute_keys::Bool = true)
    permute_keys && (keys_fill = permute(keys_fill))
    return Dict(k => 0. for k in keys_fill)
end

function fill_zero(keys_fill::Vararg{Any}; permute_keys::Bool = true)
    permute_keys && (keys_fill = permute(keys_fill))
    return Dict(k => 0. for k in keys_fill)
end

# ----- EDIT EXISTING ----------------------------------------------------------------------

function fill_zero(
    df::Vararg{DataFrame};
    permute_keys::Bool = true,
    colnames = [],
    with::Union{Dict,NamedTuple} = Dict(),
)
    df = ensurearray(df)
    N = length(df)

    id = _generate_id(N)
    idx = findindex.(df)
    to = findvalue.(df)

    df = indexjoin(df)
    df = fill_zero(df; permute_keys=permute_keys, colnames=colnames, with=with)
    from = [propertynames_with(df, k) for k in id]

    df = [edit_with(df[:,[idx[ii];from[ii]]], Rename.(from[ii],to[ii])) for ii in 1:N]

    return Tuple(df)
end

function fill_zero(
    df::DataFrame;
    permute_keys::Bool = true,
    colnames = [],
    with::Union{Dict,NamedTuple} = Dict(),
)
    idx = findindex(df)

    !isempty(colnames) && intersect!(idx, ensurearray(colnames))

    if !isempty(with)
        df = indexjoin(df, _intersect_with(with, df))
    elseif permute_keys
        df = indexjoin(permute(df[:,idx]), df)
    end

    return df
end

function fill_zero(d::Vararg{Dict}; permute_keys::Bool = true)
    d = copy.(ensurearray(d))
    # Find all keys present in the input dictionary/ies and ensure all are present.
    keys_fill = unique([collect.(keys.(d))...;])
    d = [fill_zero(keys_fill, x; permute_keys = permute_keys) for x in d]
    return length(d) == 1 ? d[1] : Tuple(d)
end

# ----- GIVEN SET LIST ---------------------------------------------------------------------

function fill_zero(keys_fill::NamedTuple, df::DataFrame)
    @warn("Depreciated!")
    df = copy(df)
    df_fill = fill_zero(keys_fill)
    df = fill_zero(df, df_fill; permute_keys = false)[1]
    return df
end

function fill_zero(set::Dict, df::DataFrame)
    @warn("Depreciated!")
    idx = intersect(findindex(df), collect(keys(set)))
    val = [set[k] for k in idx]
    return fill_zero(NamedTuple{Tuple(idx,)}(val,), df)
end

function fill_zero(keys_fill::Any, d::Dict; permute_keys::Bool = true)
    @warn("Depreciated!")
    d = copy(d)
    # If permuting keys, find all possible permutations of keys that should be present
    # and determine which are missing. Then add missing keys to the dictionary and return.
    permute_keys && (keys_fill = permute(keys_fill))
    keys_missing = setdiff(keys_fill, collect(keys(d)))

    [push!(d, k => 0.) for k in keys_missing]
    return d
end

"""
Initialize a new DataFrame and fills it with the specified input value.
"""
function fill_with(keys_fill::NamedTuple, value::Any; value_colnames = :value)
    df = fill_zero(keys_fill; value_colnames = value_colnames)
    df = edit_with(df, Replace.(value_colnames, 0.0, value))
    return df
end

"""
    permute(x::Any)
    permute(x...)

# Arguments
- `x::Any`: input data to permute

# Returns
All possible permutations of the input values.
    - `x::DataFrame` or `x::NamedTuple`: Input type and key/column names will be preserved.
    - Given any other input `x`,
        - `x::Array{Tuple,1}` of the possible combinations of the input data.
            Each tuple will be ordered in the same way the input data was ordered.
        - `x::Any`: If `x` does not contain multiple sets to permute, `permute` will return
            `x`, unchanged.
"""
function permute(df::DataFrame)
    idx = propertynames(df)
    val = DataFrame(permute(unique.(eachcol(df))))
    return edit_with(val, Rename.(propertynames(val), idx))
end

function permute(x::NamedTuple)
    idx = keys(x)
    val = eachcol(DataFrame(ensuretuple.(permute(values(x)))))
    return NamedTuple{Tuple(idx, )}(val, )
end

permute(x::Any) = any(isarray.(x)) ? permute(ensurearray.(x)...) : x
permute(x::Vararg{Any}) = vec(collect(Iterators.product(x...)))

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
@inline function Config!(kwargs_dict::Dict)
    config_kwargs_dict = Dict{Symbol, Any}()
    params = [
        (:gh_hostname_for_api,   "REGISTRATOR_GH_HOSTNAME_FOR_API",   "https://api.github.com"),
        (:gh_hostname_for_clone, "REGISTRATOR_GH_HOSTNAME_FOR_CLONE", "github.com"),
        (:registry,              "REGISTRATOR_REGISTRY"),
        (:serverless_repo,       "REGISTRATOR_SERVERLESS_REPO"),
    ]
    for param in params
        key = param[1]
        env_var = param[2]
        if haskey(kwargs_dict, key)
            config_kwargs_dict[key] = kwargs_dict[:registry]
            delete!(kwargs_dict, :registry)
        else
            if length(param) == 3
                default_value = param[3]
                config_kwargs_dict[key] = get(ENV, env_var, default_value)
            else
                config_kwargs_dict[key] = ENV[env_var]
            end
        end
    end
    github_token = Credential()
    if haskey(kwargs_dict, :github_token)
        github_token[] = kwargs_dict[:github_token]
        delete!(kwargs_dict, :github_token)
    else
        github_token[] = ENV["REGISTRATOR_GITHUB_TOKEN"]
    end
    config = Config(;
        github_token,
        config_kwargs_dict...
    )
    return config

    if haskey(kwargs_dict, :registry)
        registry = kwargs_dict[:registry]
        delete!(kwargs_dict, :registry)
    else
        registry = ENV[""]
    end
    if haskey(kwargs_dict, :serverless_repo)
        serverless_repo = kwargs_dict[:serverless_repo]
        delete!(kwargs_dict, :serverless_repo)
    else
        serverless_repo = ENV["REGISTRATOR_SERVERLESS_REPO"]
    end

end


text:
const ρ = 1000.0 # For lack of better information (kg/m3)
const c_p = 4013.0 # From Onken et a. (2007)


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "src")
AUX_DIR = joinpath(SRC_DIR, "auxiliary")
OUT_DIR=joinpath(Base.source_dir(), "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
#include(joinpath(SRC_DIR, "fitzHughNagumo.jl"))
#include(joinpath(SRC_DIR, "fitzHughNagumo_conjugateUpdt.jl"))

include(joinpath(SRC_DIR, "types.jl"))
include(joinpath(SRC_DIR, "vern7.jl"))
include(joinpath(SRC_DIR, "tsit5.jl"))
include(joinpath(SRC_DIR, "rk4.jl"))
include(joinpath(SRC_DIR, "ralston3.jl"))
include(joinpath(SRC_DIR, "priors.jl"))
include(joinpath(SRC_DIR, "guid_prop_bridge.jl"))

include(joinpath(SRC_DIR, "bounded_diffusion_domain.jl"))
include(joinpath(SRC_DIR, "lorenz_system.jl"))

include(joinpath(SRC_DIR, "random_walk.jl"))
include(joinpath(SRC_DIR, "blocking_schedule.jl"))
include(joinpath(SRC_DIR, "starting_pt.jl"))
include(joinpath(SRC_DIR, "mcmc.jl"))
include(joinpath(SRC_DIR, "path_to_wiener.jl"))

using StaticArrays
using Distributions # to define priors
using Random        # to seed the random number generator

# Let's generate the data
# -----------------------
using Bridge
include(joinpath(AUX_DIR, "data_simulation_fns.jl"))
Random.seed!(4)
θ₀ = [10.0, 28.0, 8.0/3.0, 3.0, 3.0, 3.0]
Pˣ = Lorenz(θ₀...)

x0, dt, T = ℝ{3}(1.5, -1.5, 25.0), 1/5000, 4.0
tt = 0.0:dt:T
XX, _ = simulate_segment(ℝ{3}(0.0, 0.0, 0.0), x0, Pˣ, tt)
XX
using Plots

plt = plot3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
             title="Lorenz Attractor", marker = 2)
idxes = [i for i in 1:100:length(tt)]
@gif for i in idxes
    push!(plt, XX.yy[i][1], XX.yy[i][2], XX.yy[i][3])
end every 1


text:
#Functions to deal with unified backend

export getActiveBackend, setBackend, getAvailableBackends, getBackendCount, getBackend
export AF_BACKEND_CPU, AF_BACKEND_CUDA, AF_BACKEND_OPENCL

function getActiveBackend()
    backend = af_getActiveBackend()
    if backend == AF_BACKEND_CPU
        println("CPU Backend")
    elseif backend == AF_BACKEND_CUDA
        println("CUDA Backend")
    elseif backend == AF_BACKEND_OPENCL
        println("OPENCL Backend")
    end
end

setBackend(backend::Cxx.CppEnum{:af_backend}) = af_setBackend(backend)

function getAvailableBackends()
    n = af_getAvailableBackends()
    if n == 0
        println("None")
    elseif n == 1
        println("CPU")
    elseif n == 2
        println("CUDA")
    elseif n == 3
        println("CPU and CUDA")
    elseif n == 4
        println("OpenCL")
    elseif n == 5
        println("CPU and OpenCL")
    elseif n == 6
        println("CUDA and OpenCL")
    elseif n == 7
        println("CPU, CUDA and OpenCL")
    end
end

getBackendCount() = af_getBackendCount()
function getBackend(a::AFAbstractArray)
    b = af_getBackendId(a)
    if b == AF_BACKEND_CPU
        println("CPU Backend")
    elseif b == AF_BACKEND_OPENCL
        println("OpenCL Backend")
    elseif b == AF_BACKEND_CUDA
        println("CUDA Backend")
    end
end


text:
module Test_Prismatic_zAxis

using  Modia3D
import Modia3D.ModiaMath

@assembly FallingBall(;h=1.0) begin
   world  = Object3D()
   sphere = Object3D( Modia3D.Solid(Modia3D.SolidSphere(0.1), "Aluminium") )

   # Constrain sphere movement (initial placement at position [0,h,0])
   prismatic = Modia3D.Prismatic(world, sphere, axis=3, s_start=h)
end

gravField = Modia3D.UniformGravityField(n=[0,0,-1])
fallBall = FallingBall(h=3.5, sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.3))
simulationModel = Modia3D.SimulationModel( fallBall, stopTime=1.0 )
result          = ModiaMath.simulate!(simulationModel)
ModiaMath.plot(result, ("prismatic.s", "prismatic.v"))

println("... success of Test_Prismatic_zAxis.jl!")

end


text:
#=
Created November 2, 2013

Author: Spencer Lyon

This Julia script implements routines described in 'Numerically Stable
and Accurate Stochastic Simulation Approaches for Solving Dynamic
Economic Models' by Kenneth L. Judd, Lilia Maliar and Serguei Maliar,
(2011), Quantitative Economics 2/2, 173-210 (henceforth, JMM, 2011).'

As such, it is adapted from the file 'Num_Stab_Approx.m'.
=#
using MathProgBase  # load in linprog function

function normalize_data(X, Y, intercept=true)

    T, n = size(X)
    X1 = (X[:, 2:n] .- ones(T) * mean(X[:, 2:n], 1)) ./ (ones(T) * std(X[:, 2:n], 1))
    Y1 = (Y .- ones(T) * mean(Y)) ./ (ones(T) * std(Y))
    return X1, Y1
end

function de_normalize(X, Y, beta)

    # Allocate memory for new coefs
    T, n = size(X)
    B = Array(Float64, size(beta, 1) + 1, size(beta, 2))

    # Infer de-normalized slope coefficients
    B[2:end, :] = (1.0 ./ std(X[:, 2:n], 1)') * std(Y) .* beta

    # Infer intercept from others.
    B[1, :] .= mean(Y) .- mean(X[:, 2:n], 1) * B[2:end, :]

    return B
end

function OLS(X, Y, normalize=true)
    # Normal OLS.
    # Verified on 11-2-13
    T, n = size(X)
    if normalize
        X1, Y1 = normalize_data(X, Y)
        B1 = X1 \ Y1
        B = de_normalize(X, Y, B1)
    else
        B = X \ Y
    end

    return B
end

function LS_SVD(X, Y, normalize=true)
    # OLS using singular value decomposition
    # Verified on 11-2-13
    if normalize
        X1, Y1 = normalize_data(X, Y)
        U, S, V = svd(X1, thin=true)
        S_inv = diagm(0 => 1.0 ./ S)
        B1 = V*S_inv * U' * Y1
        B = de_normalize(X, Y, B1)
    else
        U, S, V = svd(X, thin=true)
        S_inv = diagm(1.0 ./ S)
        B = V * S_inv * U' * Y
    end
    return B
end

function LAD_PP(X, Y, normalize=true)
    T, n1 = size(X)
    N = size(Y, 2)

    if normalize
        n1 -= 1
        X1, Y1 = normalize_data(X, Y)
    else
        X1 = X
        Y1 = Y
    end

    #lower and upper bound
    LB = [zeros(n1)-100; zeros(2*T)]
    UB = [zeros(n1)+100; fill(Inf, 2T)]

    f = [zeros(n1); ones(2*T)]
    Aeq =  [X1 Matrix(I, T, T) -Matrix(I, T, T)]
    B1 = zeros(size(X1, 2), N)

    for j = 1:N
       beq = Y1[:, j]
       sol = linprog(f, Aeq, '=', beq, LB, UB)
       B1[:, j] = sol.sol[1:n1]
    end

    B = normalize ? de_normalize(X, Y, B1) : B1

    return B
end

function LAD_DP(X, Y, normalize=true)
    T, n1 = size(X)
    N = size(Y, 2)

    if normalize
        n1 -= 1
        X1, Y1 = normalize_data(X, Y)
    else
        X1 = X
        Y1 = Y
    end

    LB = - ones(T)
    UB = ones(T)
    Aeq = X1'
    beq = zeros(n1)
    B1 = zeros(size(X1, 2), N)

    for j=1:N
        f = - Y1[:, j]
        sol = linprog(f, Aeq, '=', beq, LB, UB)
        B1[:, j] .= -sol.attrs[:lambda][1:n1]
    end

    B = normalize ? de_normalize(X, Y, B1) : B1

    return B
end

function RLS_T(X, Y, penalty=-7)

    T, n = size(X)
    n1 = n - 1
    X1, Y1 = normalize_data(X, Y)
    # B1 = inv(X1' * X1 + T / n1 * eye(n1) * 10.0^ penalty) * X1' * Y1
    B1 = inv(X1' * X1 + T / n1 * I * 10.0^ penalty) * X1' * Y1
    B = de_normalize(X, Y, B1)
    return B
end

function RLS_TSVD(X, Y, penalty=7)
    T, n = size(X)
    n1 = n - 1
    X1, Y1 = normalize_data(X, Y)
    U, S, V = svd(X1; thin=true)
    r = sum((maximum(S)./ S) .<= 10.0^penalty)
    Sr_inv = zeros(Float64, n1, n1)
    Sr_inv[1:r, 1:r] = diagm(1./ S[1:r])
    B1 = V*Sr_inv*U'*Y1
    B = de_normalize(X, Y, B1)
    return B
end

function RLAD_PP(X, Y, penalty=7)
    # TODO: There is a bug here. linprog returns wrong answer, even when
    #       MATLAB gets it right (lame)

    T, n1 = size(X)
    N = size(Y, 2)

    n1 -= 1
    X1, Y1 = normalize_data(X, Y)

    LB = 0.0  # lower bound is 0
    UB =  Inf  # no upper bound

    f = [10.0^penalty*ones(n1*2)*T/n1; ones(2*T)]
    # Aeq =  [X1 -X1 eye(T, T) -eye(T,T)]
    Aeq =  [X1 -X1 I -I]

    B1 = zeros(size(X1, 2), N)

    for j=1:N
        beq = Y1[:, j]
        sol = linprog(f, Aeq, '=', beq, LB, UB)
        xlp = sol.sol
        B1[:, j] = xlp[1:n1] - xlp[n1+1:2*n1]
    end

    B = de_normalize(X, Y, B1)

    return B

end