text:
# Robot and scene models

abstract type DrawObject <: Function end

# Modified ROV model from GrabCAD
# http://grabcad.com/library/rov-7
mutable struct DrawROV <: DrawObject
  data
  visid::Int
  symbol::Symbol
  offset::AffineMap
end

mutable struct DrawScene <: DrawObject
  data
  symbol::Symbol
  offset::AffineMap
end

function (dmdl::DrawROV)(vc,
                         am::AbstractAffineMap  )
  #
  settransform!(vc[:models][dmdl.symbol], am ∘ dmdl.offset)
  nothing
end
function (dmdl::DrawROV)(vc,
                         t::Translation,
                         R::Rotation  )
  #
  dmdl(vc, Translation ∘ LinearMap(R))
end
function (dmdl::DrawROV)(vc)
  setgeometry!(vc[:models][dmdl.symbol], dmdl.data)
  dmdl(vc, Translation(0.,0,0) ∘ LinearMap(Rotations.Quat(1.0,0,0,0)) )
  nothing
end

function (dmdl::DrawScene)(vc,
                           am::AbstractAffineMap  )
  #
  settransform!(vc[:env][dmdl.symbol], am ∘ dmdl.offset)
  nothing
end
function (dmdl::DrawScene)(vc,
                           t::Translation,
                           R::Rotation  )
  #
  dmdl(vc, Translation ∘ LinearMap(R))
end
function (dmdl::DrawScene)(vc)
  setgeometry!(vc[:env][dmdl.symbol], dmdl.data)
  dmdl(vc, Translation(0.,0,0) ∘ LinearMap(Rotations.Quat(1.0,0,0,0)) )
  nothing
end

ext: jl


text:
####################################################################################################
####################################### Colors #####################################################
####################################################################################################

const LIBRARIES = [:Octavian, :MKL, :OpenBLAS, :blis, :Tullio, :Gaius, :LoopVectorization, :Generic, :RecursiveFactorization, :TriangularSolve];
"""
Defines the mapping between libraries and colors
"""# #0071c5 == Intel Blue
# make sure colors are distinguishable against white background by adding white to the seed list,
# then deleting it from the resultant palette
palette = distinguishable_colors(length(LIBRARIES) + 2, [colorant"white", colorant"black", colorant"#66023C", colorant"#0071c5"])
deleteat!(palette, 1); deleteat!(palette, 1)
const COLOR_MAP = Dict(zip(LIBRARIES, palette))
getcolor(l::Symbol) = COLOR_MAP[l]
for (alias,ref) ∈ [(:BLIS,:blis),(:generic,:Generic),(:GENERIC,:Generic)]
    COLOR_MAP[alias] = COLOR_MAP[ref]
end

const JULIA_LIBS = Set(["Octavian", "Tullio", "Gaius", "Generic", "GENERIC", "generic", "RecursiveFactorization", "TriangularSolve"])
isjulialib(x) = x ∈ JULIA_LIBS

####################################################################################################
####################################### Plots ######################################################
####################################################################################################

function pick_suffix(desc = "")
    suffix = if Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f)))
        "AVX512"
    elseif Bool(VectorizationBase.has_feature(Val(:x86_64_avx2)))
        "AVX2"
    elseif Bool(VectorizationBase.has_feature(Val(:x86_64_avx)))
        "AVX"
    else
        "REGSIZE$(Int(VectorizationBase.register_size()))"
    end
    if desc != ""
        suffix *= '_' * desc
    end
    "$(Sys.CPU_NAME)_$suffix"
end

function _pkgdir()
    return dirname(dirname(@__FILE__))
end

"""
    default_plot_directory()
"""
function default_plot_directory()
    return joinpath(_pkgdir(), "docs", "src", "assets")
end

"""
    default_plot_filename(br::BenchmarkResult;
                          desc,
                          logscale)
"""
function default_plot_filename(br::BenchmarkResult{T};
                               desc::AbstractString,
                               logscale::Bool) where {T}
    l, u = extrema(br.sizes)
    if logscale
        desc *= "_logscale"
    end
    desc = (br.threaded ? "_multithreaded" : "_singlethreaded") * desc
    suffix = pick_suffix(desc)
    return "gemm_$(string(T))_$(l)_$(u)_$(suffix)"
end

"""
    plot(br::BenchmarkResult;
         desc = "",
         logscale = false,
         width = 1200,
         height = 600,
         measure = :minimum,
         plot_directory = default_plot_directory(),
         plot_filename = default_plot_filename(br; desc = desc, logscale = logscale),
         file_extensions = ["svg", "png"],
         displayplot = true)

`measure` refers to the BenchmarkTools summary on times. Valid options are:
`:minimum`, `:medain`, `:mean`, `:maximum`, and `:hmean`.

 -  `:minimum` would yield the maximum `GFLOPS`, and would be the usual estimate used in Julia.
 - `:hmean`, the harmonic mean of the times, is usful if you want an average GFLOPS, instead of a GFLOPS computed with the average times.
"""
function Gadfly.plot(br::BenchmarkResult{T}; kwargs...) where {T}
    _plot(br; kwargs...)
end
roundint(x) = round(Int,x)
# `_plot` is just like `plot`, except _plot returns the filenames
function _plot(
    br::BenchmarkResult{T};
    desc::AbstractString = "",
    logscale::Bool = false,
    width = 12inch,
    height = 8inch,
    measure = :minimum,
    plot_directory::AbstractString = default_plot_directory(),
    plot_filename::AbstractString = default_plot_filename(br; desc = desc, logscale = logscale),
    file_extensions = ["svg", "png"],
    displayplot = true
) where {T}
    j = get_measure_index(measure) # throw early if `measure` invalid
    colors = getcolor.(br.libraries);
    libraries = string.(br.libraries)
    xscale = logscale ? Scale.x_log10(labels=string ∘ roundint ∘ exp10) : Scale.x_continuous
    plt = plot(
        Gadfly.Guide.manual_color_key("Libraries", libraries, colors),
        Guide.xlabel("Size"), Guide.ylabel("GFLOPS"), xscale#, xmin = minsz, xmax = maxsz
    )
    for i ∈ eachindex(libraries)
        linestyle = isjulialib(libraries[i]) ? :solid : :dash
        l = layer(
            x = br.sizes, y = br.gflops[:,i,j],
            Geom.line, Theme(default_color = colors[i], line_style = [linestyle])
        )
        push!(plt, l)
    end
    minsz, maxsz = extrema(br.sizes)
    if logscale
      l10min = log10(minsz); l10max = log10(maxsz);
      push!(plt, Stat.xticks(ticks = range(l10min, l10max, length=round(Int,(1+2*(l10max-l10min))))))
    end
    displayplot && display(plt)
    mkpath(plot_directory)
    _filenames = String[]
    extension_dict = Dict("svg" => SVG, "png" => PNG, "pdf" => PDF, "ps" => PS)
    for ext in file_extensions
        _filename = joinpath(plot_directory, "$(plot_filename).$(ext)")
        draw(extension_dict[ext](_filename, width, height), plt)
        push!(_filenames, _filename)
    end
    return _filenames
end


text:
# Autogenerated wrapper script for Leptonica_jll for aarch64-linux-gnu
export convertfilestopdf, convertfilestops, convertformat, convertsegfilestopdf, convertsegfilestops, converttopdf, converttops, fileinfo, imagetops, liblept, xtractprotos

using Giflib_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using libwebp_jll
using Zlib_jll
using OpenJpeg_jll
JLLWrappers.@generate_wrapper_header("Leptonica")
JLLWrappers.@declare_executable_product(convertfilestopdf)
JLLWrappers.@declare_executable_product(convertfilestops)
JLLWrappers.@declare_executable_product(convertformat)
JLLWrappers.@declare_executable_product(convertsegfilestopdf)
JLLWrappers.@declare_executable_product(convertsegfilestops)
JLLWrappers.@declare_executable_product(converttopdf)
JLLWrappers.@declare_executable_product(converttops)
JLLWrappers.@declare_executable_product(fileinfo)
JLLWrappers.@declare_executable_product(imagetops)
JLLWrappers.@declare_library_product(liblept, "liblept.so.5")
JLLWrappers.@declare_executable_product(xtractprotos)
function __init__()
    JLLWrappers.@generate_init_header(Giflib_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, libwebp_jll, Zlib_jll, OpenJpeg_jll)
    JLLWrappers.@init_executable_product(
        convertfilestopdf,
        "bin/convertfilestopdf",
    )

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

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

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

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

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

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

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

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

    JLLWrappers.@init_library_product(
        liblept,
        "lib/liblept.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, OptimalTransmissionRouting

Documenter.makedocs(
    modules = OptimalTransmissionRouting,
    format = Documenter.HTML(),
    sitename = "OptimalTransmissionRouting",
    authors = "Hakan Ergun",
    pages = [
        "Home" => "index.md"
    ]
)

Documenter.deploydocs(
    target = "build",
    repo = "github.com/Electa-Git/OptimalTransmissionRouting.jl.git",
    branch = "gh-pages",
    devbranch = "main",
    versions = ["stable" => "v^", "v#.#"],
    push_preview = false
)


text:
int_rules_9_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 9.4*Miscellaneous*integration*rules =#
    @apply_utils Antiderivative((~(a') + ~(b') * (~F)((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~(n') / (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2), ~x) => ((2 * ~e * ~g) / (~C * (~e * ~f - ~d * ~g))) * Subst(Antiderivative((~a + ~b * (~F)(~c * ~x)) ^ ~n / ~x, ~x), ~x, sqrt(~d + ~e * ~x) / sqrt(~f + ~g * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~C, ~F], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && (EqQ(~B * ~e * ~g - ~C * (~e * ~f + ~d * ~g), 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * (~F)((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~(n') / (~(A') + ~(C') * (~x) ^ 2), ~x) => ((2 * ~e * ~g) / (~C * (~e * ~f - ~d * ~g))) * Subst(Antiderivative((~a + ~b * (~F)(~c * ~x)) ^ ~n / ~x, ~x), ~x, sqrt(~d + ~e * ~x) / sqrt(~f + ~g * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~C, ~F], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && (EqQ(~e * ~f + ~d * ~g, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * (~F)((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~n / (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2), ~x) => Unintegrable((~a + ~b * (~F)((~c * sqrt(~d + ~e * ~x)) / sqrt(~f + ~g * ~x))) ^ ~n / (~A + ~B * ~x + ~C * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~C, ~F, ~n], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && (EqQ(~B * ~e * ~g - ~C * (~e * ~f + ~d * ~g), 0) && Not(IGtQ(~n, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * (~F)((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~n / (~A + ~(C') * (~x) ^ 2), ~x) => Unintegrable((~a + ~b * (~F)((~c * sqrt(~d + ~e * ~x)) / sqrt(~f + ~g * ~x))) ^ ~n / (~A + ~C * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~C, ~F, ~n], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && (EqQ(~e * ~f + ~d * ~g, 0) && Not(IGtQ(~n, 0))))
    @apply_utils Antiderivative(~u / ~y, ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * log(RemoveContent(~y, ~x)) <-- Not(FalseQ(q)))
    @apply_utils Antiderivative(~u / (~y * ~w), ~x) => With([q = DerivativeDivides(~y * ~w, ~u, ~x)], q * log(RemoveContent(~y * ~w, ~x)) <-- Not(FalseQ(q)))
    @apply_utils Antiderivative(~u * (~y) ^ ~(m'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], (q * (~y) ^ (~m + 1)) / (~m + 1) <-- Not(FalseQ(q))) <-- FreeQ(~m, ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(~u * (~y) ^ ~(m') * (~z) ^ ~(n'), ~x) => With([q = DerivativeDivides(~y * ~z, ~u * (~z) ^ (~n - ~m), ~x)], (q * (~y) ^ (~m + 1) * (~z) ^ (~m + 1)) / (~m + 1) <-- Not(FalseQ(q))) <-- FreeQ([~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(~u, ~x) => With([v = SimplifyIntegrand(~u, ~x)], Antiderivative(v, ~x) <-- SimplerIntegrandQ(v, ~u, ~x))
    @apply_utils Antiderivative(~(u') * (~(e') * sqrt(~(a') + ~(b') * (~x) ^ ~(n')) + ~(f') * sqrt(~(c') + ~(d') * (~x) ^ ~(n'))) ^ ~m, ~x) => (~a * (~e) ^ 2 - ~c * (~f) ^ 2) ^ ~m * Antiderivative(ExpandIntegrand(~u * (~e * sqrt(~a + ~b * (~x) ^ ~n) - ~f * sqrt(~c + ~d * (~x) ^ ~n)) ^ -(~m), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (ILtQ(~m, 0) && EqQ(~b * (~e) ^ 2 - ~d * (~f) ^ 2, 0))
    @apply_utils Antiderivative(~(u') * (~(e') * sqrt(~(a') + ~(b') * (~x) ^ ~(n')) + ~(f') * sqrt(~(c') + ~(d') * (~x) ^ ~(n'))) ^ ~m, ~x) => (~b * (~e) ^ 2 - ~d * (~f) ^ 2) ^ ~m * Antiderivative(ExpandIntegrand(~u * (~x) ^ (~m * ~n) * (~e * sqrt(~a + ~b * (~x) ^ ~n) - ~f * sqrt(~c + ~d * (~x) ^ ~n)) ^ -(~m), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (ILtQ(~m, 0) && EqQ(~a * (~e) ^ 2 - ~c * (~f) ^ 2, 0))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') * (~u) ^ ~n + ~v) ^ ~(p') * ~w, ~x) => Antiderivative((~u) ^ (~m + ~n * ~p) * (~a + (~u) ^ -(~n) * ~v) ^ ~p * ~w, ~x) <-- FreeQ([~a, ~m, ~n], ~x) && (IntegerQ(~p) && (Not(GtQ(~n, 0)) && Not(FreeQ(~v, ~x))))
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * ~y) ^ ~(m') * (~(c') + ~(d') * ~v) ^ ~(n'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && EqQ(~v, ~y)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * ~y) ^ ~(m') * (~(c') + ~(d') * ~v) ^ ~(n') * (~(e') + ~(f') * ~w) ^ ~(p'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (EqQ(~v, ~y) && EqQ(~w, ~y))
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * ~y) ^ ~(m') * (~(c') + ~(d') * ~v) ^ ~(n') * (~(e') + ~(f') * ~w) ^ ~(p') * (~(g') + ~(h') * ~z) ^ ~(q'), ~x) => With([r = DerivativeDivides(~y, ~u, ~x)], r * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x), ~x, ~y) <-- Not(FalseQ(r))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x) && (EqQ(~v, ~y) && (EqQ(~w, ~y) && EqQ(~z, ~y)))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~y) ^ ~n), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], ~a * Antiderivative(~u, ~x) + ~b * q * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~y) ^ ~n) ^ ~p, ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~n, ~p], ~x)
    @apply_utils Antiderivative(~(u') * (~v) ^ ~(m') * (~(a') + ~(b') * (~y) ^ ~n) ^ ~(p'), ~x) => Module([q, r], q * r * Subst(Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(r = Divides((~y) ^ ~m, (~v) ^ ~m, ~x))) && Not(FalseQ(q = DerivativeDivides(~y, ~u, ~x)))) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~y) ^ ~n + ~(c') * (~v) ^ ~(n2')) ^ ~p, ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ(~v, ~y))
    @apply_utils Antiderivative(~(u') * (~A + ~(B') * (~y) ^ ~n) * (~(a') + ~(b') * (~v) ^ ~n + ~(c') * (~w) ^ ~(n2')) ^ ~(p'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~A + ~B * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~v, ~y) && EqQ(~w, ~y)))
    @apply_utils Antiderivative(~(u') * (~A + ~(B') * (~y) ^ ~n) * (~(a') + ~(c') * (~w) ^ ~(n2')) ^ ~(p'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~A + ~B * (~x) ^ ~n) * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~c, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ(~w, ~y))
    @apply_utils Antiderivative(~(u') * (~v) ^ ~(m') * (~(a') + ~(b') * (~y) ^ ~n + ~(c') * (~w) ^ ~(n2')) ^ ~(p'), ~x) => Module([q, r], q * r * Subst(Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(r = Divides((~y) ^ ~m, (~v) ^ ~m, ~x))) && Not(FalseQ(q = DerivativeDivides(~y, ~u, ~x)))) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ(~w, ~y))
    @apply_utils Antiderivative(~(u') * (~z) ^ ~(m') * (~A + ~(B') * (~y) ^ ~n) * (~(a') + ~(b') * (~v) ^ ~n + ~(c') * (~w) ^ ~(n2')) ^ ~(p'), ~x) => Module([q, r], q * r * Subst(Antiderivative((~x) ^ ~m * (~A + ~B * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(r = Divides((~y) ^ ~m, (~z) ^ ~m, ~x))) && Not(FalseQ(q = DerivativeDivides(~y, ~u, ~x)))) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~v, ~y) && EqQ(~w, ~y)))
    @apply_utils Antiderivative(~(u') * (~z) ^ ~(m') * (~A + ~(B') * (~y) ^ ~n) * (~(a') + ~(c') * (~w) ^ ~(n2')) ^ ~(p'), ~x) => Module([q, r], q * r * Subst(Antiderivative((~x) ^ ~m * (~A + ~B * (~x) ^ ~n) * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(r = Divides((~y) ^ ~m, (~z) ^ ~m, ~x))) && Not(FalseQ(q = DerivativeDivides(~y, ~u, ~x)))) <-- FreeQ([~a, ~c, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && EqQ(~w, ~y))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~y) ^ ~n) ^ ~(m') * (~(c') + ~(d') * (~v) ^ ~n) ^ ~(p'), ~x) => With([q = DerivativeDivides(~y, ~u, ~x)], q * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~m * (~c + ~d * (~x) ^ ~n) ^ ~p, ~x), ~x, ~y) <-- Not(FalseQ(q))) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && EqQ(~v, ~y)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~y) ^ ~n) ^ ~(m') * (~(c') + ~(d') * (~v) ^ ~n) ^ ~(p') * (~(e') + ~(f') * (~w) ^ ~n) ^ ~(q'), ~x) => With([r = DerivativeDivides(~y, ~u, ~x)], r * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~m * (~c + ~d * (~x) ^ ~n) ^ ~p * (~e + ~f * (~x) ^ ~n) ^ ~q, ~x), ~x, ~y) <-- Not(FalseQ(r))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p, ~q], ~x) && (EqQ(~v, ~y) && EqQ(~w, ~y))
    @apply_utils Antiderivative(~u * (~F) ^ ~v, ~x) => With([q = DerivativeDivides(~v, ~u, ~x)], (q * (~F) ^ ~v) / log(~F) <-- Not(FalseQ(q))) <-- FreeQ(~F, ~x)
    @apply_utils Antiderivative(~u * (~w) ^ ~(m') * (~F) ^ ~v, ~x) => With([q = DerivativeDivides(~v, ~u, ~x)], q * Subst(Antiderivative((~x) ^ ~m * (~F) ^ ~x, ~x), ~x, ~v) <-- Not(FalseQ(q))) <-- FreeQ([~F, ~m], ~x) && EqQ(~w, ~v)
    @apply_utils Antiderivative(~u * (~a + ~(b') * (~v) ^ ~(p') * (~w) ^ ~(p')) ^ ~(m'), ~x) => With([c = Simplify(~u / (~w * D(~v, ~x) + ~v * D(~w, ~x)))], c * Subst(Antiderivative((~a + ~b * (~x) ^ ~p) ^ ~m, ~x), ~x, ~v * ~w) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p], ~x) && IntegerQ(~p)
    @apply_utils Antiderivative(~u * (~a + ~(b') * (~v) ^ ~(p') * (~w) ^ ~(q')) ^ ~(m') * (~v) ^ ~(r'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) + ~q * ~v * D(~w, ~x)))], ((c * ~p) / (~r + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ (~p / (~r + 1))) ^ ~m, ~x), ~x, (~v) ^ (~r + 1) * ~w) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q, ~r], ~x) && (EqQ(~p, ~q * (~r + 1)) && (NeQ(~r, -1) && IntegerQ(~p / (~r + 1))))
    @apply_utils Antiderivative(~u * (~a + ~(b') * (~v) ^ ~(p') * (~w) ^ ~(q')) ^ ~(m') * (~v) ^ ~(r') * (~w) ^ ~(s'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) + ~q * ~v * D(~w, ~x)))], ((c * ~p) / (~r + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ (~p / (~r + 1))) ^ ~m, ~x), ~x, (~v) ^ (~r + 1) * (~w) ^ (~s + 1)) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q, ~r, ~s], ~x) && (EqQ(~p * (~s + 1), ~q * (~r + 1)) && (NeQ(~r, -1) && IntegerQ(~p / (~r + 1))))
    @apply_utils Antiderivative(~u * (~(a') * (~v) ^ ~(p') + ~(b') * (~w) ^ ~(q')) ^ ~(m'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) - ~q * ~v * D(~w, ~x)))], c * ~p * Subst(Antiderivative((~b + ~a * (~x) ^ ~p) ^ ~m, ~x), ~x, ~v * (~w) ^ (~m * ~q + 1)) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q], ~x) && (EqQ(~p + ~q * (~m * ~p + 1), 0) && (IntegerQ(~p) && IntegerQ(~m)))
    #= Antiderivative((~u)*((~a')*(~v)^(~p')+(~b')*(~w)^(~q'))^(~m'),~x) :=  With([c=Simplify(u/(p*w*D(v,x)-q*v*D(w,x)))], -c*q*Subst(Antiderivative((a+b*x^q)^m,x),x,v^(m*p+1)*w) <-- FreeQ(c,x)) <--  FreeQ([a,b,m,p,q],x) && EqQ(p+q*(m*p+1),0) && IntegerQ(q) &&  IntegerQ(m) =#
    @apply_utils Antiderivative(~u * (~(a') * (~v) ^ ~(p') + ~(b') * (~w) ^ ~(q')) ^ ~(m') * (~v) ^ ~(r'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) - ~q * ~v * D(~w, ~x)))], -c * ~q * Subst(Antiderivative((~a + ~b * (~x) ^ ~q) ^ ~m, ~x), ~x, (~v) ^ (~m * ~p + ~r + 1) * ~w) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q, ~r], ~x) && (EqQ(~p + ~q * (~m * ~p + ~r + 1), 0) && (IntegerQ(~q) && IntegerQ(~m)))
    @apply_utils Antiderivative(~u * (~(a') * (~v) ^ ~(p') + ~(b') * (~w) ^ ~(q')) ^ ~(m') * (~w) ^ ~(s'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) - ~q * ~v * D(~w, ~x)))], ((-c * ~q) / (~s + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ (~q / (~s + 1))) ^ ~m, ~x), ~x, (~v) ^ (~m * ~p + 1) * (~w) ^ (~s + 1)) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q, ~s], ~x) && (EqQ(~p * (~s + 1) + ~q * (~m * ~p + 1), 0) && (NeQ(~s, -1) && (IntegerQ(~q / (~s + 1)) && IntegerQ(~m))))
    @apply_utils Antiderivative(~u * (~(a') * (~v) ^ ~(p') + ~(b') * (~w) ^ ~(q')) ^ ~(m') * (~v) ^ ~(r') * (~w) ^ ~(s'), ~x) => With([c = Simplify(~u / (~p * ~w * D(~v, ~x) - ~q * ~v * D(~w, ~x)))], ((-c * ~q) / (~s + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ (~q / (~s + 1))) ^ ~m, ~x), ~x, (~v) ^ (~m * ~p + ~r + 1) * (~w) ^ (~s + 1)) <-- FreeQ(c, ~x)) <-- FreeQ([~a, ~b, ~m, ~p, ~q, ~r, ~s], ~x) && (EqQ(~p * (~s + 1) + ~q * (~m * ~p + ~r + 1), 0) && (NeQ(~s, -1) && (IntegerQ(~q / (~s + 1)) && IntegerQ(~m))))
    @apply_utils Antiderivative(~u * (~x) ^ ~(m'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative(SubstFor((~x) ^ (~m + 1), ~u, ~x), ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ(~m, ~x) && (NeQ(~m, -1) && FunctionOfQ((~x) ^ (~m + 1), ~u, ~x))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfLinear(u, x)], ShowStep("", "Antiderivative(F((a+b*x)^(1/n),x),x)", "n/b*Subst(Antiderivative(x^(n-1)*F(x,-a/b+x^n/b),x),x,(a+b*x)^(1/n))", Hold( lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))))) <-- Not(FalseQ(lst)) && SubstForFractionalPowerQ(u, lst((3)), x)) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfLinear(u, x)], lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))) <-- Not(FalseQ(lst)) && SubstForFractionalPowerQ(u, lst((3)), x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfQuotientOfLinears(u, x)], ShowStep("", "Antiderivative(F(((a+b*x)/(c+d*x))^(1/n),x),x)", "n*(b*c-a*d)*Subst(Antiderivative(x^(n-1)*F(x,(-a+c*x^n)/(b-d*x^n))/(b-d*x^ n)^2,x),x,((a+b*x)/(c+d*x))^(1/n))", Hold( lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))))) <-- Not(FalseQ(lst))) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfQuotientOfLinears(u, x)], lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))) <-- Not(FalseQ(lst))))
    @apply_utils Antiderivative(~(u') * (~(a') * (~v) ^ ~(m') * (~w) ^ ~(n') * (~z) ^ ~(q')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a * (~v) ^ ~m * (~w) ^ ~n * (~z) ^ ~q) ^ FracPart(~p)) / ((~v) ^ (~m * FracPart(~p)) * (~w) ^ (~n * FracPart(~p)) * (~z) ^ (~q * FracPart(~p)))) * Antiderivative(~u * (~v) ^ (~m * ~p) * (~w) ^ (~n * ~p) * (~z) ^ (~p * ~q), ~x) <-- FreeQ([~a, ~m, ~n, ~p, ~q], ~x) && (Not(IntegerQ(~p)) && (Not(FreeQ(~v, ~x)) && (Not(FreeQ(~w, ~x)) && Not(FreeQ(~z, ~x)))))
    @apply_utils Antiderivative(~(u') * (~(a') * (~v) ^ ~(m') * (~w) ^ ~(n')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a * (~v) ^ ~m * (~w) ^ ~n) ^ FracPart(~p)) / ((~v) ^ (~m * FracPart(~p)) * (~w) ^ (~n * FracPart(~p)))) * Antiderivative(~u * (~v) ^ (~m * ~p) * (~w) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(FreeQ(~v, ~x)) && Not(FreeQ(~w, ~x))))
    @apply_utils Antiderivative(~(u') * (~(a') * (~v) ^ ~(m')) ^ ~p, ~x) => (((~a) ^ IntPart(~p) * (~a * (~v) ^ ~m) ^ FracPart(~p)) / (~v) ^ (~m * FracPart(~p))) * Antiderivative(~u * (~v) ^ (~m * ~p), ~x) <-- FreeQ([~a, ~m, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(FreeQ(~v, ~x)) && (Not(EqQ(~a, 1) && EqQ(~m, 1)) && Not(EqQ(~v, ~x) && EqQ(~m, 1)))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~b) ^ IntPart(~p) * (~a + ~b * (~x) ^ ~n) ^ FracPart(~p)) / ((~x) ^ (~n * FracPart(~p)) * (1 + (~a * (~x) ^ -(~n)) / ~b) ^ FracPart(~p))) * Antiderivative(~u * (~x) ^ (~n * ~p) * (1 + (~a * (~x) ^ -(~n)) / ~b) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (Not(IntegerQ(~p)) && (ILtQ(~n, 0) && (Not(RationalFunctionQ(~u, ~x)) && IntegerQ(~p + 1 / 2))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~v) ^ ~n) ^ ~p, ~x) => ((~a + ~b * (~v) ^ ~n) ^ FracPart(~p) / ((~v) ^ (~n * FracPart(~p)) * (~b + ~a * (~v) ^ -(~n)) ^ FracPart(~p))) * Antiderivative(~u * (~v) ^ (~n * ~p) * (~b + ~a * (~v) ^ -(~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (Not(IntegerQ(~p)) && (ILtQ(~n, 0) && (BinomialQ(~v, ~x) && Not(LinearQ(~v, ~x)))))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~(m') * (~v) ^ ~n) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~m * (~v) ^ ~n) ^ FracPart(~p) / ((~v) ^ (~n * FracPart(~p)) * (~b * (~x) ^ ~m + ~a * (~v) ^ -(~n)) ^ FracPart(~p))) * Antiderivative(~u * (~v) ^ (~n * ~p) * (~b * (~x) ^ ~m + ~a * (~v) ^ -(~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~m, ~p], ~x) && (Not(IntegerQ(~p)) && (ILtQ(~n, 0) && BinomialQ(~v, ~x)))
    @apply_utils Antiderivative(~(u') * (~(a') * (~x) ^ ~(r') + ~(b') * (~x) ^ ~(s')) ^ ~m, ~x) => With([v = (~a * (~x) ^ ~r + ~b * (~x) ^ ~s) ^ FracPart(~m) / ((~x) ^ (~r * FracPart(~m)) * (~a + ~b * (~x) ^ (~s - ~r)) ^ FracPart(~m))], v * Antiderivative(~u * (~x) ^ (~m * ~r) * (~a + ~b * (~x) ^ (~s - ~r)) ^ ~m, ~x) <-- NeQ(Simplify(v), 1)) <-- FreeQ([~a, ~b, ~m, ~r, ~s], ~x) && (Not(IntegerQ(~m)) && PosQ(~s - ~r))
    @apply_utils Antiderivative(~u / (~a + ~(b') * (~x) ^ ~n), ~x) => With([v = RationalFunctionExpand(~u / (~a + ~b * (~x) ^ ~n), ~x)], Antiderivative(v, ~x) <-- SumQ(v)) <-- FreeQ([~a, ~b], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(~u * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(~p) && Not(AlgebraicFunctionQ(~u, ~x)))))
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(~u * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (Not(IntegerQ(~p)) && Not(AlgebraicFunctionQ(~u, ~x)))))
    @apply_utils Antiderivative(~u / (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')), ~x) => With([v = RationalFunctionExpand(~u / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x)], Antiderivative(v, ~x) <-- SumQ(v)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && IGtQ(~n, 0))
    @apply_utils Antiderivative(~(u') / (~(a') * (~x) ^ ~(m') + ~(b') * sqrt(~(c') * (~x) ^ ~n)), ~x) => Antiderivative((~u * (~a * (~x) ^ ~m - ~b * sqrt(~c * (~x) ^ ~n))) / ((~a) ^ 2 * (~x) ^ (2 * ~m) - (~b) ^ 2 * ~c * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n], ~x)
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = FunctionOfLinear(u, x)], ShowStep("", "Antiderivative(F(a+b*x),x)", "Subst(Antiderivative(F(x),x),x,a+b*x)/b", Hold( Dist(1/lst((3)), Subst(Antiderivative(lst((1)), x), x, lst((2)) + lst((3))*x), x))) <-- Not(FalseQ(lst))) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([lst = FunctionOfLinear(u, x)], Dist(1/lst((3)), Subst(Antiderivative(lst((1)), x), x, lst((2)) + lst((3))*x), x) <-- Not(FalseQ(lst))))
    If(TrueQ($LoadShowSteps), Antiderivative((~u)/(~x), ~x) := With([lst = PowerVariableExpn(u, 0, x)], ShowStep("", "Antiderivative(F((c*x)^n)/x,x)", "Subst(Antiderivative(F(x)/x,x),x,(c*x)^n)/n", Hold( 1/lst((2))* Subst(Antiderivative(NormalizeIntegrand(Simplify(lst((1))/x), x), x), x, (lst((3))*x)^lst((2))))) <-- Not(FalseQ(lst)) && NeQ(lst((2)), 0)) <-- SimplifyFlag && NonsumQ(u) && Not(RationalFunctionQ(u, x)), Antiderivative((~u)/(~x), ~x) := With([lst = PowerVariableExpn(u, 0, x)], 1/lst((2))* Subst(Antiderivative(NormalizeIntegrand(Simplify(lst((1))/x), x), x), x, (lst((3))*x)^lst((2))) <-- Not(FalseQ(lst)) && NeQ(lst((2)), 0)) <-- NonsumQ(u) && Not(RationalFunctionQ(u, x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u)*(~x)^(~m'), ~x) := With([lst = PowerVariableExpn(u, m + 1, x)], ShowStep("If*g=GCD(m+1,n)>1,", "Antiderivative(x^m*F(x^n),x)", "1/g*Subst(Antiderivative(x^((m+1)/g-1)*F(x^(n/g)),x),x,x^g)", Hold( 1/lst((2))* Subst(Antiderivative(NormalizeIntegrand(Simplify(lst((1))/x), x), x), x, (lst((3))*x)^lst((2))))) <-- Not(FalseQ(lst)) && NeQ(lst((2)), m + 1)) <-- SimplifyFlag && IntegerQ(m) && NeQ(m, -1) && NonsumQ(u) && (GtQ(m, 0) || Not(AlgebraicFunctionQ(u, x))), Antiderivative((~u)*(~x)^(~m'), ~x) := With([lst = PowerVariableExpn(u, m + 1, x)], 1/lst((2))* Subst(Antiderivative(NormalizeIntegrand(Simplify(lst((1))/x), x), x), x, (lst((3))*x)^lst((2))) <-- Not(FalseQ(lst)) && NeQ(lst((2)), m + 1)) <-- IntegerQ(m) && NeQ(m, -1) && NonsumQ(u) && (GtQ(m, 0) || Not(AlgebraicFunctionQ(u, x))))
    @apply_utils Antiderivative((~x) ^ ~m * ~u, ~x) => With([k = Denominator(~m)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * ReplaceAll(~u, (~x->begin
                                            #= none:1 =#
                                            #= none:3 =#
                                            (~x) ^ k
                                        end)), ~x), ~x, (~x) ^ (1 / k))) <-- FractionQ(~m)
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], ShowStep("", "Antiderivative(F(sqrt(a+b*x+c*x^2),x),x)", "2*Subst(Antiderivative(F((c*sqrt(a)-b*x+sqrt(a)*x^2)/(c-x^2),(-b+2* sqrt(a)*x)/(c-x^2))* (c*sqrt(a)-b*x+sqrt(a)*x^2)/(c-x^2)^2,x),x,(-sqrt(a)+sqrt(a+ b*x+c*x^2))/x)", Hold(2*Subst(Antiderivative(lst((1)), x), x, lst((2))))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 1)) <-- SimplifyFlag && EulerIntegrandQ(u, x), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], 2*Subst(Antiderivative(lst((1)), x), x, lst((2))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 1)) <-- EulerIntegrandQ(u, x))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], ShowStep("", "Antiderivative(F(sqrt(a+b*x+c*x^2),x),x)", "2*Subst(Antiderivative(F((a*sqrt(c)+b*x+sqrt(c)*x^2)/(b+2*sqrt(c)*x),(- a+x^2)/(b+2*sqrt(c)*x))* (a*sqrt(c)+b*x+sqrt(c)*x^2)/(b+2*sqrt(c)*x)^2,x),x,sqrt(c)* x+sqrt(a+b*x+c*x^2))", Hold(2*Subst(Antiderivative(lst((1)), x), x, lst((2))))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 2)) <-- SimplifyFlag && EulerIntegrandQ(u, x), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], 2*Subst(Antiderivative(lst((1)), x), x, lst((2))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 2)) <-- EulerIntegrandQ(u, x))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], ShowStep("", "Antiderivative(F(sqrt(a+b*x+c*x^2),x),x)", "-2*sqrt(b^2-4*a*c)*Subst(Antiderivative(F(-sqrt(b^2-4*a*c)*x/(c-x^2), (b*c+c*sqrt(b^2-4*a*c)+(-b+sqrt(b^2-4*a*c))*x^2)/(-2*c*(c- x^2)))*x/(c-x^2)^2,x), x,2*c*sqrt(a+b*x+c*x^2)/(b-sqrt(b^2-4*a*c)+2*c*x))", Hold(2*Subst(Antiderivative(lst((1)), x), x, lst((2))))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 3)) <-- SimplifyFlag && EulerIntegrandQ(u, x), Antiderivative((~u), ~x) := With([lst = FunctionOfSquareRootOfQuadratic(u, x)], 2*Subst(Antiderivative(lst((1)), x), x, lst((2))) <-- Not(FalseQ(lst)) && EqQ(lst((3)), 3)) <-- EulerIntegrandQ(u, x))
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~v) ^ 2), ~x) => (1 / (2 * ~a)) * Antiderivative(Together(1 / (1 - ~v / Rt(-(~a) / ~b, 2))), ~x) + (1 / (2 * ~a)) * Antiderivative(Together(1 / (1 + ~v / Rt(-(~a) / ~b, 2))), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~v) ^ ~n), ~x) => Dist(2 / (~a * ~n), Sum(Antiderivative(Together(1 / (1 - (~v) ^ 2 / ((-1) ^ ((4k) / ~n) * Rt(-(~a) / ~b, ~n / 2)))), ~x), [k, 1, ~n / 2]), ~x) <-- FreeQ([~a, ~b], ~x) && IGtQ(~n / 2, 1)
    @apply_utils Antiderivative(1 / (~a + ~(b') * (~v) ^ ~n), ~x) => Dist(1 / (~a * ~n), Sum(Antiderivative(Together(1 / (1 - ~v / ((-1) ^ ((2k) / ~n) * Rt(-(~a) / ~b, ~n)))), ~x), [k, 1, ~n]), ~x) <-- FreeQ([~a, ~b], ~x) && IGtQ((~n - 1) / 2, 0)
    @apply_utils Antiderivative(~v / (~a + ~(b') * (~u) ^ ~(n')), ~x) => Antiderivative(ReplaceAll(ExpandIntegrand(PolynomialInSubst(~v, ~u, ~x) / (~a + ~b * (~x) ^ ~n), ~x), (~x->begin
                            #= none:1 =#
                            #= none:3 =#
                            ~u
                        end)), ~x) <-- FreeQ([~a, ~b], ~x) && (IGtQ(~n, 0) && PolynomialInQ(~v, ~u, ~x))
    @apply_utils Antiderivative(~u, ~x) => With([v = NormalizeIntegrand(~u, ~x)], Antiderivative(v, ~x) <-- v !== ~u)
    @apply_utils Antiderivative(~u, ~x) => With([v = ExpandIntegrand(~u, ~x)], Antiderivative(v, ~x) <-- SumQ(v))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~(m')) ^ ~(p') * (~(c') + ~(d') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (((~a + ~b * (~x) ^ ~m) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q) / (~x) ^ (~m * ~p)) * Antiderivative(~u * (~x) ^ (~m * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p, ~q], ~x) && (EqQ(~a + ~d, 0) && (EqQ(~b + ~c, 0) && (EqQ(~m + ~n, 0) && EqQ(~p + ~q, 0))))
    @apply_utils Antiderivative(~u * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (sqrt(~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) / ((4 * ~c) ^ (~p - 1 / 2) * (~b + 2 * ~c * (~x) ^ ~n))) * Antiderivative(~u * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p - 1 / 2)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfLinear(u, x)], ShowStep("", "Antiderivative(F((a+b*x)^(1/n),x),x)", "n/b*Subst(Antiderivative(x^(n-1)*F(x,-a/b+x^n/b),x),x,(a+b*x)^(1/n))", Hold( lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))))) <-- Not(FalseQ(lst))) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([lst = SubstForFractionalPowerOfLinear(u, x)], lst((2))*lst((4))* Subst(Antiderivative(lst((1)), x), x, lst((3))^(1/lst((2)))) <-- Not(FalseQ(lst))))
    Antiderivative((~u), (~x)) := CannotIntegrate(u, x)
end

max_stars_repo_path: src/rules/9 Miscellaneous/9.4 Miscellaneous integration rules.jl


text:
using Distributions: Gamma

"type-parameterized implementation of dirichlet multinomial models"
type DirichletMultinomial{T} # T specifies the outcome value type
    n              :: Int # number of outcomes
    params         :: Dict{T, Float64} # hyperparameters
    initial_params :: Dict{T, Float64}
end
params(dm::DirichletMultinomial) = dm.params # accessor function

function DirichletMultinomial(n, params) # constructor method
    DirichletMultinomial(n, params, params)
end

"""
Return the log of the probability of the outcome d.
The outcome d is modeled by a muliset which is implemented
as a dictionary, mapping outcome values to their counts, or
any iterable object of outcome values.
"""
function logscore{T}(dm::DirichletMultinomial{T}, d::Dict{T, Int})
    s = sum(values(dm.params))
    LogProb((factorial(dm.n)*gamma(s) / gamma(dm.n+s)) *
    prod(gamma(d[k]+dm.params[k])/(factorial(d[k])*gamma(dm.params[k])) for k in keys(d)))
end

function logscore{T}(dm::DirichletMultinomial{T}, iter)
    d = Dict{T, Int}()
    for x in iter
        if haskey(d, x)
            d[x] += 1
        else
            d[x] = 1
        end
    end
    logscore(dm, d)
end

"update the parameters of dm with params"
function set_params!{T}(dm::DirichletMultinomial{T}, params)
    dm.params = dm.initial_params + params
end

function sample(dm::DirichletMultinomial)
    multinomial_weights = Dict(k => rand(Gamma(dm.params[k])) for k in keys(dm.params))
    samples = Dict(k=>0 for k in keys(dm.params))
    for i in 1:dm.n
        rand_number = rand() * sum(values(multinomial_weights))
        bound = 0.0
        for k in keys(dm.params)
            bound += multinomial_weights[k]
            if rand_number <= bound
                samples[k] += 1
                break
            end
        end
    end
    samples
end

add_obs!(dm::DirichletMultinomial, obs) = dm.params[obs] += 1

"""
Dirichlet categorical models are Dirichlet multinomial models
where the number of outcomes is fixed to one.
"""
type DirichletCategorical{T}
    dm :: DirichletMultinomial{T}
end

function DirichletCategorical{T}(params::Dict{T, Float64})
    DirichletCategorical(DirichletMultinomial(1, params))
end

function logscore{T}(dc::DirichletCategorical{T}, x::T)
    logscore(dc.dm, (x,))
end

function set_params!{T}(dc::DirichletCategorical{T}, params)
    set_params!(dc.dm, params)
end

function sample(dc::DirichletCategorical)
    d = sample(dc.dm)
    for key in keys(d)
        if d[key] == 1
            return key
        end
    end
end

add_obs!(dc::DirichletCategorical, obs) = add_obs!(dc.dm, obs)

# dirichlet categorical distribution with integer hyperparameters
type UrnModel{T}
    counts     :: Dict{T,Int}
    sum_counts :: Int
end
counts(um::UrnModel) = um.counts
sum_counts(um::UrnModel) = um.sum_counts

UrnModel{T}(dict::Dict{T,Int}) = UrnModel(dict, sum(values(dict)))
UrnModel(iter) = UrnModel(Dict{eltype(iter),Int}(x=>1 for x in iter))

logscore{T}(um::UrnModel{T}, x::T) = LogProb(counts(um)[x] / sum_counts(um))

function add_obs!(um::UrnModel, obs)
    counts(um)[obs] += 1
    um.sum_counts += 1
    obs
end

function remove_obs!(um::UrnModel, obs)
    @assert counts(um)[obs] > 1
    counts(um)[obs] -= 1
    um.sum_counts -= 1
    nothing
end

function sample(um::UrnModel)
    vals = keys(counts(um))
    weights = values(counts(um))
    rand_num = rand(1:sum_counts(um))
    bound = 0
    for (v,w) in zip(vals,weights)
        bound += w
        if rand_num <= bound
            return v
        end
    end
end

type SimpleCond{C, D, T} # C...context, D...distribution
    dists :: Dict{C, D}
    support :: Vector{T}
end

SimpleCond(dists) = SimpleCond(dists, dists |> keys |> collect)

logscore(cond::SimpleCond, x, context) = logscore(cond.dists[context], x)
sample(cond::SimpleCond, context) = sample(cond.dists[context])
function add_obs!{C,D}(cond::SimpleCond{C,D}, obs, context)
    if !haskey(cond.dists, context)
        cond.dists[context] = D(cond.support)
    end
    add_obs!(cond.dists[context], obs)
end

# urn model conditional
type UrnModelCond{T,C} # C ... context type
    support :: Vector{T}
    dists   :: Dict{C,UrnModel{T}}
end
UrnModelCond(support, C::DataType) =
    UrnModelCond(collect(support), Dict{C,UrnModel{eltype(support)}}())

function logscore(cond::UrnModelCond, x, context)
    if !haskey(cond.dists, context)
        cond.dists[context] = UrnModel(cond.support)
    end
    logscore(cond.dists[context], x)
end
function sample(cond::UrnModelCond, context)
    if !haskey(cond.dists, context)
        cond.dists[context] = UrnModel(cond.support)
    end
    sample(cond.dists[context])
end
function add_obs!(cond::UrnModelCond, obs, context)
    if !haskey(cond.dists, context)
        cond.dists[context] = UrnModel(cond.support)
    end
    add_obs!(cond.dists[context], obs)
end
function remove_obs!(cond::UrnModelCond, obs, context)
    remove_obs!(cond.dists[context], obs)
end


text:
ng = NEAT.NodeGene(1, :INPUT)
@test typeof(ng.id) == Int64
@test typeof(ng.ntype) == Symbol
@test typeof(ng.bias) == Float64
@test typeof(ng.response) == Float64
@test typeof(ng.activation) == Symbol
@test typeof(ng.timeConstant) == Float64

ng = NEAT.NodeGene(1, :INPUT, 0.,4.9, :sigm, 1.5)
@test ng.ntype == :INPUT
@test ng.ntype != :HIDDEN
@test ng.ntype != :OUTPUT

ng = NEAT.NodeGene(1, :HIDDEN, 0.,4.9, :sigm, 1.5)
@test ng.ntype != :INPUT
@test ng.ntype == :HIDDEN
@test ng.ntype != :OUTPUT

ng = NEAT.NodeGene(1, :OUTPUT, 0.,4.9, :sigm, 1.5)
@test ng.ntype != :INPUT
@test ng.ntype != :HIDDEN
@test ng.ntype == :OUTPUT

child = NEAT.get_child(ng,ng)
@test child.id           == ng.id
@test child.ntype        == ng.ntype
@test child.bias         == ng.bias
@test child.response     == ng.response
@test child.activation   == ng.activation
@test child.timeConstant == ng.timeConstant

ng = NEAT.NodeGene(1, :HIDDEN, 0., 1., :sigm)
NEAT.mutate_bias!(ng, config)
@test ng.bias != 0

ng = NEAT.NodeGene(1, :HIDDEN, 0., 1., :sigm)
NEAT.mutate_response!(ng,config)
@test ng.response != 1

ng = NEAT.NodeGene(1, :HIDDEN, 0., 1., :sigm)
NEAT.mutate!(ng, config) # make sure it runs :()

g = NEAT.Global(config)
cg =  NEAT.ConnectionGene(g, 1, 2, .5, true)
cg2 =  NEAT.ConnectionGene(g, 1, 2, .5, true)

# test that two connection with identical in and outs DO have same key and innvotion numbers
@test cg.innovNumber == cg2.innovNumber
@test cg.key == cg2.key

# test that two connection with diffent outs do NOT have same key or innvotion numbers
cg3 =  NEAT.ConnectionGene(g, 1, 3, .5, true)
@test cg3.innovNumber != cg2.innovNumber
@test cg3.key != cg2.key

cg =  NEAT.ConnectionGene(NEAT.Global(config), 1, 2, .5, true)
w = cg.weight
NEAT.mutate_weight!(cg,config)
@test cg.weight != w

cg =  NEAT.ConnectionGene(NEAT.Global(config), 1, 2, .5, true); w = cg.weight
NEAT.weight_replaced!(cg,config)
@test cg.weight != w

# test that mutate someitme actual does mutate (better test to come)
# @test sum(map(x->(NEAT.mutate!(cg, config)== true) ? 1 : 0,[1:100])) > 0

g = NEAT.Global(config)
cg =  NEAT.ConnectionGene(g, 1, 2, .5, true)
@test length(g.innovations) == 1
cg2, cg3 = NEAT.split(g,cg,3)
@test length(g.innovations) == 3
@test cg2.inId == cg.inId
@test cg2.outId != cg.outId
@test cg3.inId != cg.inId
@test cg3.outId == cg.outId

@test NEAT.is_same_innov(cg, cg) == true
@test NEAT.is_same_innov(cg, cg2) == false
@test NEAT.is_same_innov(cg, cg3) == false
@test NEAT.is_same_innov(cg2, cg3) == false

# test for choice between one or other
cg3 = NEAT.get_child(cg,cg2)
@test (cg3==cg || cg3 == cg2) == true
@test (cg3==cg && cg3 == cg2) == false

# different object but with same values
cg2 = NEAT.deepcopy(cg)
@test cg2 != cg
@test cg2.inId == cg.inId
@test cg2.outId == cg.outId
@test cg2.weight == cg.weight
@test cg2.enable == cg.enable
@test cg2.innovNumber == cg.innovNumber


text:
# module Bukdu

include("renderers/text.jl")
include("renderers/json.jl")
include("renderers/html.jl")
include("renderers/markdown.jl")
include("renderers/mustache.jl")

module ViewFilter
filters = Dict()
end # module Bukdu.ViewFilter

for func in [plugins, before, after]
    name = Base.function_name(func)
    function add_view_filter(block::Function, render_func::Function, typ_name)
        params = tuple(methods(block).mt.defs.func.sig.parameters[2:end]...)
        key = (render_func,typ_name,params)
        ViewFilter.filters[(func,key)] = block
    end
    @eval begin
        function $name{AL<:ApplicationLayout}(block::Function, render_func::Function, D::LayoutDivision{AL})
            typ_name = viewlayout_symbol(D)
            $add_view_filter(block,render_func,typ_name)
        end

        function $name(block::Function, render_func::Function, modul::Module)
            typ_name = Val{Base.module_name(modul)}
            $add_view_filter(block,render_func,typ_name)
        end

        function $name(block::Function, render_func::Function, T::Type)
            typ_name = T.name.name
            $add_view_filter(block,render_func,typ_name)
        end
    end
end

function filtering(render_block::Function, render_func::Function, T::Type, args...)::Conn
    typ_name = (:Val == T.name.name) ? T : T.name.name
    params = map(x->Any, args)
    key = (render_func,typ_name,params)
    if haskey(ViewFilter.filters, (plugins,key))
        ViewFilter.filters[(plugins,key)](args...)
    end
    if haskey(ViewFilter.filters, (before,key))
        f = ViewFilter.filters[(before,key)]
        ViewFilter.filters[(before,key)](args...)
    end
    conn = render_block()
    if haskey(ViewFilter.filters, (after,key))
        ViewFilter.filters[(after,key)](args...)
    end
    conn
end

function render(modul::Module, args...; kw...)::Conn
    V = Val{Base.module_name(modul)}
    render(V, args...; kw...)
end

function render{AL<:ApplicationLayout}(D::LayoutDivision{AL}, args...; kw...)::Conn
    V = isa(D.dividend, Module) ? Val{Base.module_name(D.dividend)} : D.dividend
    L = D.divisor
    params = map(x->Any,args)
    key = (render,viewlayout_symbol(D),params)
    if haskey(ViewFilter.filters, (plugins,key))
        ViewFilter.filters[(plugins,key)](args...)
    end
    if haskey(ViewFilter.filters, (before,key))
        ViewFilter.filters[(before,key)](args...)
    end
    conn::Conn = render(V, args...; kw...)
    conn_body = conn.resp_body
    if isempty(kw)
        bodies = tuple(conn_body, args[2:end]...)
    else
        argstuple = isempty(args) ? tuple() : map(typeof, args)
        if method_exists(layout, tuple(L, Any, argstuple..., Dict))
            bodies = tuple(conn_body, args..., Dict(kw))
        else
            bodies = tuple(conn_body, args...)
        end
    end
    conn.resp_body = layout(L(), bodies...)
    if haskey(ViewFilter.filters, (after,key))
        ViewFilter.filters[(after,key)](args...)
    end
    conn
end


text:
# too many allocations for time type, probably we need another strategy for this
function buff_parser(res, lbuff, cc, nd, current_line, df, ::Type{T}) where  T <: TimeType
    flag = 0
    if cc > nd
        val = nothing
    else
        val = Dates.tryparsenext_internal(T, lbuff, cc, nd, df, false)
    end

    if val === nothing
        @simd for i in cc:nd
            @inbounds if (lbuff.data[i] != 0x20 && lbuff.data[i] != 0x2e)
                flag = 1

            end
        end
        res[current_line[]] = missing
    else
        res[current_line[]] = T(val[1]...)
    end
    flag
end

function buff_parser(res, lbuff, cc, nd, current_line, ::Type{T}; base = 10) where T <: Integer
    val = Base.tryparse_internal(T, lbuff, cc, nd, base, false)
    flag = 0
    # hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
    # (Ptr{UInt8},Csize_t,Csize_t), lbuff, cc-1, nd - cc +1)
    # hasvalue ? res[current_line[]] = T(val) : res[current_line[]] = missing
    if val === nothing
        @simd for i in cc:nd
            @inbounds if (lbuff.data[i] != 0x20 && lbuff.data[i] != 0x2e)
                flag = 1

            end
        end
        res[current_line[]] = missing
    else
        res[current_line[]] = val
    end
    flag
    # (x, code, startpos, value_len, total_len) = Parsers.xparse(T, lbuff, cc, nd)
    # code == 33 ? res[current_line[]] = x : x = missing
end
function buff_parser(res, lbuff, cc, nd, current_line, ::Type{T}) where T <: Real
    hasvalue, val = ccall(:jl_try_substrtod, Tuple{Bool, Float64},
    (Ptr{UInt8},Csize_t,Csize_t), lbuff, cc-1, nd - cc +1)
    flag = 0
    if hasvalue
        res[current_line[]] = val
    else
        @simd for i in cc:nd
            @inbounds if (lbuff[i] != 0x20 && lbuff[i] != 0x2e)
                flag = 1

            end
        end
        res[current_line[]] = missing
    end
    flag
    # RES = Parsers.xparse(T, lbuff, cc, nd)
    # RES.code == 33 ? res[current_line[]] = RES.val : res[current_line[]] = missing
    # return 0
end
function buff_parser(res, lbuff, cc, nd, current_line, ::Type{Float32})
    hasvalue, val = ccall(:jl_try_substrtof, Tuple{Bool, Float32},
    (Ptr{UInt8},Csize_t,Csize_t), lbuff, cc-1, nd - cc +1)
    flag = 0
    if hasvalue
        res[current_line[]] = val
    else
        @simd for i in cc:nd
            @inbounds if (lbuff[i] != 0x20 && lbuff[i] != 0x2e)
                flag = 1

            end
        end
        res[current_line[]] = missing
    end
    flag
    # (x, code, startpos, value_len, total_len) = Parsers.xparse(Float32, lbuff, cc, nd)
    # code == 33 ? res[current_line[]] = x : x = missing
end

function buff_parser(res, lbuff, cc, nd, current_line, ::Type{BigFloat})
    newlo = cc
    newhi = nd
    for i in cc:nd
        lbuff[i] == 0x20 ? newlo += 1 : break
    end
    for i in nd:-1:newlo
        lbuff[i] == 0x20 ? newhi -= 1 : break
    end
    if newhi-newlo+1 == 0
        res[current_line[]] = missing
        return 0
    end
    z = BigFloat(precision=Base.MPFR.DEFAULT_PRECISION[])
    err = ccall((:mpfr_set_str, :libmpfr), Int32, (Ref{BigFloat}, Cstring, Int32, Base.MPFR.MPFRRoundingMode), z, SubString(lbuff, newlo, newhi), 0, Base.MPFR.ROUNDING_MODE[])
    if err == 0
        res[current_line[]] = z
        return 0
    else
        res[current_line[]] = missing
        return 1
    end
end

function buff_parser(res, lbuff, cc, nd, current_line, ::Type{String})
    # (x, code, startpos, value_len, total_len) = Parsers.xparse(String, lbuff, cc, nd, Parsers.Options(ignoreemptylines=true))
    # code == 33 ? res[current_line[]] = x : x = missing
    l = nd - cc + 1
    cnt = 0
    for i in cc:nd
        cnt += lbuff[i] === 0x20
    end
    l == 0 || l == cnt ? res[current_line[]] = missing : res[current_line[]] = unsafe_string(pointer(lbuff, cc), l)
    return 0

end

function buff_parser(res, lbuff, cc, nd, current_line, trim, ::Type{String})
    # (x, code, startpos, value_len, total_len) = Parsers.xparse(String, lbuff, cc, nd, Parsers.Options(ignoreemptylines=true))
    # code == 33 ? res[current_line[]] = x : x = missing
    newlo = cc
    newhi = nd
    for i in nd:-1:cc
        lbuff[i] == 0x20 ? newhi -= 1 : break
    end
    l = newhi - newlo + 1
    l == 0 ? res[current_line[]] = missing : res[current_line[]] = unsafe_string(pointer(lbuff, cc), l)
    return 0
end

function buff_parser(res, lbuff, cc, nd, current_line, ::Type{Bool})
    # (x, code, startpos, value_len, total_len) = Parsers.xparse(String, lbuff, cc, nd, Parsers.Options(ignoreemptylines=true))
    # code == 33 ? res[current_line[]] = x : x = missing
    flag = 0
    newlo = cc
    newhi = nd
    for i in cc:nd
        lbuff[i] == 0x20 ? newlo += 1 : break
    end
    for i in nd:-1:newlo
        lbuff[i] == 0x20 ? newhi -= 1 : break
    end
    if length(newlo:newhi) == 1
        if lbuff[newlo] == 0x30
            res[current_line[]] = false
        elseif lbuff[newlo] == 0x31
            res[current_line[]] = true
        else
            res[current_line[]] = missing
            flag = 1
        end
    elseif length(newlo:newhi) == 0
        res[current_line[]] = missing
    else
        res[current_line[]] = missing
        flag = 1
    end

    return flag

end
# function buff_parser(res, lbuff, cc, nd, current_line, ::Type{T}) where T <: InlineString
#     (x, code, startpos, value_len, total_len) = Parsers.xparse(T, lbuff, cc, nd, Parsers.Options())
#     code == 33 ? res[current_line[]] = x : x = missing
# end
#
# function (::Type{T})(buf::Vector{UInt8}, pos, len) where {T <: InlineString}
#    if T === InlineString1
#        len == 1 || InlineStrings.stringtoolong(T, sizeof(x))
#        return Base.bitcast(InlineString1, buf[pos])
#    else
#        length(buf) < len && InlineStrings.buftoosmall()
#        len < sizeof(T) || InlineStrings.stringtoolong(T, len)
#        y = GC.@preserve buf unsafe_load(convert(Ptr{T}, pointer(buf, pos)))
#        sz = 8 * (sizeof(T) - len)
#        return Base.or_int(Base.shl_int(Base.lshr_int(InlineStrings._bswap(y), sz), sz), Base.zext_int(T, UInt8(len)))
#    end
# end

function buff_parser(res, lbuff, cc, nd, current_line, ::Type{T}) where T <: InlineString
    if cc>nd
        res[current_line[]] = missing
    else
        l = nd - cc + 1
        cnt = 0
        for i in cc:nd
            cnt += lbuff[i] === 0x20
        end
        l == cnt ? res[current_line[]] = missing : res[current_line[]] = T(lbuff, cc, nd-cc+1)
    end
    return 0
end
function buff_parser(res, lbuff, cc, nd, current_line, char_buff, ::Type{T}) where T <: Characters
    if cc>nd
        res[current_line[]] = missing
    else
        l = nd - cc + 1
        cnt = 0
        for i in cc:nd
            cnt += lbuff[i] === 0x20
        end
        l == cnt ? res[current_line[]] = missing : res[current_line[]] = T(lbuff, char_buff, cc, nd)
    end
    return 0
end

function buff_parser(res, lbuff, cc, nd, current_line, char_buff, ::Type{T}) where T <: DT
    if cc>nd
        nothing
    else
        copyto!(res, crrent_line[], lbuff, cc, min(strlength(T), nd - cc +1))
    end
    return 0
end


text:
using JuMP
using PyPlot
using DataFrames
using StructArrays
using CSV
using LinearAlgebra
using KNITRO
using Gurobi
using Statistics

#####TESTS FOR SCR-NCP QUANTITY, TESTED VIA DIAGONALIZATION COMPLEX#######
include("EPECs/EPEC_quantity.jl")
include("diagonalization.jl")

###############DECLARE PARTICIPATING AGENTS###########
###############DECLARE PARTICIPATING AGENTS###########
rel_path = "DATA"
filepath = joinpath(@__DIR__, rel_path)
cd(filepath)
generators_dat = CSV.read("generators_data.csv", header=true, delim= ";")
MC_g = generators_dat[1,2:8]
MC_g_str = generators_dat[1,9]
#CAPACITIES OF GENERATORS
capacity_g = generators_dat[2,2:9]
capacity_g = values(capacity_g) ./ 100
capacity_g_nonstr = values(capacity_g[1:7])
capacity_g_str = values(capacity_g[8])

#RAMPING LIMIT OF THE STRATEGIC GENERATOR:
Rᴳ = 1.0
res_dat = CSV.read("wind_PV.csv", header=true, delim= ";")
eight_res  = res_dat[1:3:end,:] #selecting every 3rd of the res_dat
renewable_nonstr = ["PV"]
renewable_str = ["WIND"]
capacity_re_nonstr = eight_res.PV/100
capacity_re_str =  (eight_res.Wind/100) * 1.5

#LOADING THE DEMAND
load_dat = CSV.read("load.csv",header=true, delim= ";")
load = load_dat.load[1:3:end,:]/100

eff_ch = eff_dch = generators_dat.ESS[end]
energy_max = (generators_dat.ESS[end-2]/100)*1.5
power_max = (generators_dat.ESS[end-4]/100)*1.5

MC_re_str = MC_re_nonstr = 0 #Zero marginal cost

str_stor = ["ESS_STR"]
str_gens = ["GEN_STR", "WIND_STR"]
str_agents = vcat(str_stor, str_gens)

generators = names(MC_g)
renewables = [:PV]

non_str_gens = vcat(generators)
non_str_res = vcat(renewables)
non_str_agents = vcat(generators, renewables)

#NOTE:ATM There is no non-str ESS involved in the pool
#MARKET PRICE LIMITS
p_cap = 100
p_floor = 0

function t_interval(lenght, resolution)
    t_steps = [i for i in range(1, stop = lenght, step = resolution)]
    return t_steps
end
t_steps = t_interval(8,1)

###########INIT BID, TECH LIMIT DFs#########
function return_data()
    #DATA OF THE STR AGENTS, DIFFERENTIATE BETWEEN BIDS AND TRUE VALUES
    dat_str = DataFrame()
    for i in range(1, stop = length(str_gens)), l in range(1, stop = length(str_stor))
        #This needs to be revisited....
        if i ==  1 #if str agent is GEN
            dat_str[i] =  StructArray(
            P_max = ones(length(t_steps)) * capacity_g_str,
            G = zeros(length(t_steps)),
            c_g = zeros(length(t_steps)),
            MC_g = ones(length(t_steps)) * MC_g_str)
        elseif i == 2 #if str agent is RE
            dat_str[i] =  StructArray(
            P_max = capacity_re_str,
            G = zeros(length(t_steps)),
            c_g = zeros(length(t_steps)),
            MC_g = zeros(length(t_steps)))
        end

        dat_str[i + l] =  StructArray(
        E_max = ones(length(t_steps)) * energy_max,
        P_max = ones(length(t_steps)) * power_max,
        DCH = zeros(length(t_steps)),
        MC_dch = ones(length(t_steps)) * p_cap,
        c_dch = zeros(length(t_steps)),
        CH = zeros(length(t_steps)),
        c_ch = zeros(length(t_steps)),
        MC_ch = ones(length(t_steps)) * p_floor)
    end

    #NOTE: For non-str agents set the quantity bid equal to the maximum capacity
    #NON STR AGENTS THE BIDS ARE EQUAL TO THE TRUE VALUES
    dat_nonstr = DataFrame()
    for i in range(1, stop = length(non_str_gens)), k in range(1, stop = length(non_str_res))
        dat_nonstr[i] =  StructArray(
        G = ones(length(t_steps)) * capacity_g_nonstr[i],
        c_g = ones(length(t_steps)) * MC_g[1+i-1])
        dat_nonstr[i + k] = StructArray(
        G = capacity_re_nonstr,
        c_g = ones(length(t_steps)) * MC_re_nonstr #as zero marginal cost
        )
    end

    #NOW RENAME THEM
    names!(dat_str, vcat(Symbol.(str_gens), Symbol.(str_stor)))
    names!(dat_nonstr, Symbol.(non_str_agents))

    return dat_str, dat_nonstr
end
dat_str, dat_nonstr = return_data()

function compute_profits(EPEC)
    lambda = value.(EPEC[:λ][:])
    #NON-STRATEGIC GENERATORS
    gen_1 = value.(EPEC[:g]["GEN1",:])
    gen_2 = value.(EPEC[:g]["GEN2",:])
    gen_3 = value.(EPEC[:g]["GEN3",:])
    gen_5 = value.(EPEC[:g]["GEN5",:])
    gen_6 = value.(EPEC[:g]["GEN6",:])
    gen_7 = value.(EPEC[:g]["GEN7",:])
    gen_8 = value.(EPEC[:g]["GEN8",:])

    #DEMAND
    d = value.(EPEC[:d])
    #NON-STRATEGIC RES
    res_PV = value.(EPEC[:w][:PV,:])
    #DISPATCHED QUANTITIES FROMS STR AGENTS
    dch = value.(EPEC[:dch_st])
    ch = value.(EPEC[:ch_st])
    gen = value.(EPEC[:g_st])
    res = value.(EPEC[:w_st])
    #CALCULATING PROFITS
    earned_ESS = sum(dch[t] * lambda[t] for t in t_steps)
    paid_ESS = sum(ch[t] * lambda[t] for t in t_steps)
    profit_ESS =  earned_ESS - paid_ESS
    # println("PROFIT_ESS: ", profit_ESS)

    earned_GEN = sum(gen[t] * lambda[t] for t in t_steps)
    paid_GEN = sum(dat_str.GEN_STR.MC_g[t] * gen[t] for t in t_steps)
    profit_GEN =  earned_GEN - paid_GEN
    # println("PROFIT_GEN: ", profit_GEN)

    earned_RES = sum(res[t] * lambda[t] for t in t_steps)
    paid_RES = sum(dat_str.WIND_STR.MC_g[t] * res[t] for t in t_steps)
    profit_RES =  earned_RES - paid_RES
    # println("PROFIT_RES: ", profit_RES)

    return profit_ESS, profit_GEN, profit_RES
end

#Possible obj_types = ["empty", "competitive", "collusive", "ESS_fav"]
obj_types = ["empty"]
#The initial parameters for starting the relaxation, NOTE: The lower the value the strickter the relaxation is.
param_adj = [1,10]
flag_sequanctial_red = true
flag_statistics = true

optimizer = with_optimizer(KNITRO.Optimizer,
                # algorithm = 1,
                hessian_no_f=1,
                datacheck=0,
                outlev=0,
                # ms_enable= 1,
                # ms_maxsolves = 50
                )

#INIT DataFrame for the results
simu_EPEC_quantity = DataFrame(time_build = [], time = [], bidding = String[], SC_obj = String[], ESS = Float64[], GEN = Float64[], WIND_STR = Float64[],
                                            SW= Float64[], PS= Float64[], Avg_MP = Float64[], ϵ1 = Float64[], ϵ2 = Float64[],
                                            confirmedNE = Bool[], solver = [], ESSᴰ = Float64[], GENᴰ = Float64[], WIND_STRᴰ = Float64[])

for o in obj_types
    for n in param_adj
        global feasibility = true

        SP_obj = o
        bidding = "quantity"
        println("....Building EPEC model with objective: $o, and epsilon: $n ...")
        t_1_build = time()
        EPEC = build_EPEC_quantity_mod(t_steps, non_str_gens, non_str_res, str_agents, optimizer, SP_obj)
        t_2_build = time()
        t_build_EPEC = t_2_build - t_1_build
        println("....Solving EPEC model with objective: $o, and epsilon: $n ...")
        if flag_sequanctial_red == true
            global iter_limit = 40
            global iterations = 0
            global tol = 0.0001
            global epsilon1_new = n
            global epsilon2_new = n

            JuMP.fix(EPEC[:ϵ1], epsilon1_new; force = true)
            JuMP.fix(EPEC[:ϵ2], epsilon2_new; force = true)
            optimize!(EPEC)
            objective_value(EPEC)
            profit_ESS, profit_GEN, profit_RES = compute_profits(EPEC)
            println("PROFITS:
                            PROFIT_ESS: $profit_ESS,
                            PROFIT_GEN: $profit_GEN,
                            PROFIT_RES: $profit_RES")

            ############THE SMALL SCRIPT SNIPPET FOR UPDATING THE REGULARIZATION TERM
            t1 = time()
            while epsilon1_new >= tol || epsilon2_new >= tol
                #After reaching iteration limit exit:
                if iterations >= iter_limit
                    break
                end
                #########SETTING THE WARM START POINT TO THE PREVIOUS SOLUTION#########
                #if terminition status is local optimal, otherwise exit
                term_stat = termination_status(EPEC)
                if term_stat == MOI.OPTIMAL || term_stat == MOI.LOCALLY_SOLVED
                    set_start_value.(all_variables(EPEC), value.(all_variables(EPEC)))
                else
                    feasibility = false
                    break
                end
                ###################FIXING EPRSILONS TO THE REDUCED PARAMS##############
                global epsilon1_new = value.(EPEC[:ϵ1])*0.5
                if epsilon2_new >= tol
                    global epsilon2_new = value.(EPEC[:ϵ2])*0.5
                end
                # else
                #     set_start_value.(all_variables(EPEC), prev_vals)
                #     global epsilon_new = value.(EPEC[:ϵ])*0.99
                # end
                fix(EPEC[:ϵ1], epsilon1_new; force = true)
                fix(EPEC[:ϵ2], epsilon2_new; force = true)

                println("
                #########
                THE NEW RELAXATION TERMS:
                            # ϵ1: $epsilon1_new and
                            ϵ2: $epsilon2_new and
                            #ITERATIONS: $iterations
                            #TERMINATION: $term_stat
                #########")
                ################OPTIMIZING############
                optimize!(EPEC)
                ##############CACLCULATING PROFIT############
                profit_ESS, profit_GEN, profit_RES = compute_profits(EPEC)
                println("PROFITS:
                                PROFIT_ESS: $profit_ESS,
                                PROFIT_GEN: $profit_GEN,
                                PROFIT_RES: $profit_RES")
                ############PROGRESSING IN TIME##############
                global iterations += 1
                global prev_vals = value.(all_variables(EPEC))
            end
            t2 = time()
            t_EPEC = t2 - t1
        end

        if feasibility == true
            println("....Model was solved to a feasible point in $t_EPEC sec :) ...")
            lambda = value.(EPEC[:λ][:])
            #NON-STRATEGIC GENERATORS
            gen_1 = value.(EPEC[:g]["GEN1",:])
            gen_2 = value.(EPEC[:g]["GEN2",:])
            gen_3 = value.(EPEC[:g]["GEN3",:])
            gen_5 = value.(EPEC[:g]["GEN5",:])
            gen_6 = value.(EPEC[:g]["GEN6",:])
            gen_7 = value.(EPEC[:g]["GEN7",:])
            gen_8 = value.(EPEC[:g]["GEN8",:])

            #DEMAND
            d = value.(EPEC[:d])
            #NON-STRATEGIC RES
            res_PV = value.(EPEC[:w][:PV,:])
            #DISPATCHED QUANTITIES FROMS STR AGENTS
            dch = value.(EPEC[:dch_st])
            ch = value.(EPEC[:ch_st])
            gen = value.(EPEC[:g_st])
            res = value.(EPEC[:w_st])

            flag_statistics = true
            if flag_statistics == true
                using Statistics
                avg_mp = mean(lambda)
                max_mp = maximum(lambda)
                min_mp = minimum(lambda)
                var_mp = var(lambda)

                #CALCULATING SOCIAL WELFARE
                SW =
                    (sum(
                    d[t] .* p_cap
                    for t in t_steps)
                        +
                    sum( # minus the producers surplus
                        # -
                        # dch[t] * p_dch_bid_strat[t]
                        # +
                        # ch[t] * p_ch_bid_strat[t]
                        -
                        res[t] * dat_str.WIND_STR.MC_g[t]
                        -
                        res_PV[t] * dat_nonstr.PV.c_g[t]
                        -
                        gen_1[t] * dat_nonstr["GEN1"].c_g[t]
                        -
                        gen_2[t] * dat_nonstr["GEN2"].c_g[t]
                        -
                        gen_3[t] * dat_nonstr["GEN3"].c_g[t]
                        -
                        gen_5[t] * dat_nonstr["GEN5"].c_g[t]
                        -
                        gen_6[t] * dat_nonstr["GEN6"].c_g[t]
                        -
                        gen_7[t] * dat_nonstr["GEN7"].c_g[t]
                        -
                        gen_8[t] * dat_nonstr["GEN8"].c_g[t]
                        -
                        gen[t] * dat_str.GEN_STR.MC_g[t]
                    for t in t_steps)
                )

                println("STATISTICS:
                ....Avg market price: $avg_mp
                ....Maximum: $max_mp
                ....Minimum: $min_mp
                ....Social Welfare: $SW")
            end

            global bidding = "complex"
            iter_diag, profit_ESS_MILP_ESS, profit_GEN_MILP_GEN, profit_RES_MILP_WIND = run_diagonaluization(EPEC)

            if iter_diag == 0
                confirmed_by_diag = true
            else
                confirmed_by_diag = false
            end

            push!(simu_EPEC_quantity, (t_build_EPEC, t_EPEC, "SCR-NC", "$SP_obj", profit_ESS, profit_GEN, profit_RES,
                    SW, profit_ESS + profit_GEN + profit_RES, avg_mp, n, n, confirmed_by_diag, "KNITRO",
                    profit_ESS_MILP_ESS, profit_GEN_MILP_GEN, profit_RES_MILP_WIND))

        else
            println("....Model converged to an infeasible point, after $iterations iterations :( ...")
        end

    end
end

@show mean(simu_EPEC_quantity.SW)
@show mean(simu_EPEC_quantity.PS)
@show mean(simu_EPEC_quantity.ESS)
@show mean(simu_EPEC_quantity.GEN)
@show mean(simu_EPEC_quantity.WIND_STR)


text:
using DimensionalData
using Statistics

function slidingcarpet(x, y, Z; resolution=(800, 800), kwargs...)
    # x is a vector of times labelling columns of Z, , y is a vector labelling rows of Z and Z is the matrix of time series data
    fig = Figure(;resolution=resolution)
    # sl_a = Makie.Slider(fig[1, 1], range = 0:0.01:10, startvalue = 3)
    # sl_t = Makie.Slider(fig[2, 1], range = 0:0.01:10, startvalue = 3)
    sla = labelslider!(fig, "Window", (2:length(x)))
    fig[1, 2] = sla.layout
    slb = labelslider!(fig, "Time     ", 0.01:0.005:1)
    fig[2, 2] = slb.layout

    u = lift(sla.slider.value, slb.slider.value) do sl...
        r = floor(length(x)*sl[1]./length(x)) # The number of time points to show
        s = ceil((length(x) - r)*sl[2]) # The starting position

        extrema([s, s+r-1])
    end
    limits = lift(x -> (to_value(x)[1], to_value(x)[2], nothing, nothing), u)
    colsize!(fig.layout, 1, Relative(1/6))
    ax = Axis(fig[3, 2], limits=limits)
    heatmap!(ax, Z; kwargs...)
    xx = round.(LinRange(1, round(length(x), sigdigits=1), 5), sigdigits=1)
    xxx = round.(x[Int.(xx)], sigdigits=2)
    xxx[ismissing.(xxx)] .= NaN
    #ax.xticks = (xx, string.(xxx))
    y[ismissing.(y)] .= ""
    ax.yticks = (1:size(Z, 2), y)
    return fig
end
export slidingcarpet

function slidingcarpet(X::DimArray; resolution=(800, 400), kwargs...)
    x = dims(X, Ti).val
    y = dims(X, :channel).val
    Z = Array(X)
    slidingcarpet(x, y, Z; resolution=resolution, kwargs...)
end

function neuroslidingcarpet(X::DimArray; resolution=(800, 400), kwargs...)
    time = dims(X, Ti).val
    channels = AN.getstructureacronyms(Meta.parse.(string.(dims(X, :channel).val)))
    fig = slidingcarpet(time, channels, X.-mean(Array(X), dims=1); resolution=resolution, kwargs...)
end


text:
include("lib_for_robot.jl")
function task9(r)
    i=0
    side=West
    fl=true
    while fl
        for tmp in 0:i
            move!(r,side)
            if ismarker(r)==true
                fl=false
                break
            end
        end
        side=reverseside(side)
        if fl
            for tmp in 0:i
                move!(r,side)
                if ismarker(r)==true
                    fl=false
                    break
                end
            end
        end
        side=reverseside(side)
        i+=1
    end
end


text:
"""
LARGEST_COMPONENT
-----------------
    return the largest connected component of A.
    Acc = largest_component(A) returns the largest connected component
    of the graph A.  If A is directed, this returns the largest
    strongly connected component.

    Acc = largest_component(A,true) returns the largest connected piece of
    a directed graph where connectivity is undirected.  Algorithmically,
    this takes A, drops the directions, then components the largest component
    and returns just this piece of the original _directed_ network.  So the
    output Acc is directed in this case.

    (Acc,p) = largest_component(A) also returns a logical vector
    indicating which vertices in A were chosen.

Functions
---------
- (Acc,p) = largest_component{T}(A::SparseMatrixCSC{T,Int64})
- (Acc,p) = largest_component{T}(A::SparseMatrixCSC{T,Int64},sym::Bool)

Example
-------
~~~
A = load_matrix_network("dfs_example")
(Acc,p) = largest_component(A)
~~~
"""
function largest_component(A::SparseMatrixCSC{T,Int64}) where T
    return largest_component(A,false)
end

function largest_component(A::SparseMatrixCSC{T,Int64},sym::Bool) where T
    if sym
        # As = A|A' until Julia implements this on sparse matrices, use this:
        As = max.(A,A')
        cc = scomponents(As)
    else
        cc = scomponents(A)
    end
    cind = argmax(cc.sizes)
    p = cc.map .== cind
    # no logical indexing so:
    idx = findall(p)
    Acc = A[idx,idx]
    return (Acc,p)
end


text:
greet = "Hello"
whom = "world"
greet * ", " * whom * "."

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
immutable Client
	url::String
	auth::Union{String,Tuple{String,String}}

	Client(url,auth) = new(rstrip(url,'/'),auth)

end

function Client(url::String, user::String, password::String)
	return Client(url,(user,password))
end

function doRequest(c::Client,method::Function, path::String; kwargs...)
	kwargsDict = Dict{String,Any}(string(k) => v for (k,v) in kwargs)

	url = c.url * "/api/v1" * path

	if isa(c.auth,String) && !haskey(get(kwargsDict,"headers",Dict{String,String}()),"Authorization")
		get!(kwargsDict,"headers",Dict{String,String}())["Authorization"] = "token $(c.auth)"
	else
		replace(url,"://","://$(c.auth[1]):$(c.auth[2])@",1)
	end

	return method(url; Dict(Symbol(k) => v for (k,v) in kwargsDict)...)
end

function getResponse(c::Client, method::Function, path::String; kwargs...)
	resp = doRequest(c,method,path;kwargs...)

	status = Requests.statuscode(resp)
	status == 401 && throw(GiteaError(status,HttpCommon.STATUS_CODES[status]))
	status == 403 && throw(GiteaError(status,HttpCommon.STATUS_CODES[status]))
	status == 404 && throw(GiteaError(status,HttpCommon.STATUS_CODES[status]))

	if status÷100 != 2
		errMap = Requests.json(resp)
		DEBUG && JSON.print(errMap,2)
		throw(GiteaError(status,errMap))
	end

	return resp
end

function getParsedResponse{T}(::Type{T}, c::Client, method::Function, path::String;kwargs...)
	resp = getResponse(c,method,path;kwargs...)
	DEBUG && JSON.print(Requests.json(resp),2)
	return unmarshalJSON(T,Requests.json(resp))
end

function getStatusCode(c::Client, method::Function, path::String; kwargs...)
	return Requests.status(doRequest(c,method,path;kwargs...))
end


text:
module Op
using LinearAlgebra
using ..Interp:getD
using ..Interp:getCpts
using ..Interp:getCpts1
using ..Interp:baryW1
using ..Interp:barymat
using ..Interp:forward
using ..Interp:rev
using ..EdgeDetect:find

function setup(; N = 20, bc1 = [:a,0,0.0], bc2 = [:b,0,0.0], ranges=zeros(0,2), contin = Bool[], A0=0.0, A1=0.0, A2=0.0, B=1.0, a = -1.0, b = 1.0, dosplit=true)

    a = Float64(a)
    b = Float64(b)

    if typeof(A0) <: Number
        function A0f(x); return A0; end
        ranges0 = zeros(0,2)
        contin0 = Bool[]
    else
        A0f = A0
        if dosplit && length(ranges) == 0
            ranges0, contin0 = find(A0f,a=a,b=b)
        end
    end

    if typeof(A1) <: Number
        function A1f(x); return A1; end
        ranges1 = zeros(0,2)
        contin1 = Bool[]
    else
        A1f = A1
        if dosplit && length(ranges) == 0
            ranges1, contin1 = find(A1f,a=a,b=b)
        end
    end

    secondorder=true
    if typeof(A2) <: Number
        if abs(A2) < 1e-15
            secondorder=false
        end
        function A2f(x); return A2; end
        ranges2 = zeros(0,2)
        contin2 = Bool[]
    else
        A2f = A2
        if dosplit && length(ranges) == 0
            ranges2, contin2 = find(A2f,a=a,b=b)
        end
    end

    if typeof(B) <: Number
        function Bf(x); return B; end
        rangesB = zeros(0,2)
        continB = Bool[]
    else
        Bf = B
        if dosplit && length(ranges) == 0
            rangesB, continB = find(Bf,a=a,b=b)
        end
    end

    if length(ranges) == 0 && dosplit

        if dosplit && length(ranges) == 0
            #merge ranges
            ranges = [ranges0; ranges1;ranges2;rangesB]
            contin = [contin0; contin1; contin2; continB]
            s = sortperm(ranges[:,1])
            ranges = ranges[s,:]
            contin = contin[s,:]

#            println("ranges A")
#            println(ranges)

            if size(ranges)[1] > 0
                keep = [1]
                for i in 2:(size(ranges)[1])
                    if sum(abs.(ranges[i-1,:] - ranges[i,:])) > 1e-4
                        push!(keep, i)
                    end
                end
                ranges = ranges[keep,:]
                contin = contin[keep,:]
                if length(ranges) > 0
                    println("detected edges")
                    println(ranges)
#                    println("contin")
#                    println(contin)
#                    println()
                end
            end
        end
    end

#    println("asdf")
#    println("typeof(A2f) ", typeof(A2f))

    if secondorder
        return setup_fn(; N=N, bc1=bc1, bc2=bc2, ranges=ranges, contin = contin, A0=A0f, A1=A1f, A2=A2f, B=Bf, a = a, b = b)
    else
        return setup_fn_firstorder(; N=N, bc1=bc1, ranges=ranges, contin = contin, A0=A0f, A1=A1f, B=Bf, a = a, b = b)
    end

end

function getmats(N; a=-1.0, b = 1.0)

    D = getD(N);
    pts = getCpts(N);

    Dab = D/ ((b-a)/2.0)
    pts_f = rev(pts,a,b)

    return pts_f, Dab

end

function setup_fn(; N = 20, bc1 = [:a,0,0.0], bc2 = [:b,0,0.0], ranges=zeros(0,2), contin = Bool[], A0::Function, A1::Function, A2::Function, B::Function, a = -1.0, b = 1.0)

#    println("b")

#    ranges = ranges[contin,:]

    if typeof(N) <: Number
        N = ones(Int64, 1+size(ranges)[1]) * N
    end

    DL = []
    VL = []

    a1 = a
    b1 = a
    for c in 1:size(ranges)[1]
        b1 = ranges[c,1]
        pts, D = getmats(N[c]; a = a1+1e-20, b = b1-1e-20)
        push!(DL, D)
        push!(VL, pts)
        a1 = ranges[c,2]

    end

    b1 = b
    Vpts, Dab = getmats( N[end]; a = a1+1e-20, b = b1-1e-20)
    println("Vpts ", Vpts[1:4])
    push!(DL, Dab)
    push!(VL, Vpts)

    num_parts = length(DL)
    hsize = 0
    for (c,(d,v)) in enumerate(zip(DL,VL))
        hsize += size(d)[1]
    end

    H = zeros(hsize, hsize)
    S = zeros(hsize, hsize)

    rhs = zeros(hsize)

    Bound = zeros(0,hsize)

    function apply_bc(bc)
        if ismissing(bc)
            return hcat(1.0, zeros(1,hsize-1)), 0.0
        end
        if bc[1] == :a && bc[2] == 0
            return hcat(1.0, zeros(1,hsize-1)), bc[3]
        elseif bc[1] == :b && bc[2] == 0
            return hcat(zeros(1,hsize-1), 1.0), bc[3]
        elseif bc[1] == :a && bc[2] == 1
            return hcat(DL[1][1,:]', zeros(1,hsize-size(DL[1])[1])), bc[3]
        elseif bc[1] == :b && bc[2] == 1
#            println("s ", size(zeros(1,hsize-size(DL[end])[1])))
#            println("s2 ",  size(DL[end][end,:]))
            return hcat( zeros(1,hsize-size(DL[end])[1]), DL[end][end,:]'),bc[3]
        else
            println("WARNING, bc not recognized $bc")
            return hcat(1.0, zeros(1,hsize-1)), bc[3]
        end
    end

    println("bc1 $bc1 bc2 $bc2")

    Bound = zeros(0,hsize)
    BC1, v1 = apply_bc(bc1)

#    println("hsize $hsize")
#    println("size BC1 ", size(BC1))
#    println("size BC2 ", size(BC2))

    Bound = [Bound; BC1]

    if !ismissing(bc2)
        BC2, v2 = apply_bc(bc2)
        Bound = [Bound; BC2]
    end

    if ismissing(bc2)
        rhs[end] = v1
    else
        rhs[end-1] = v1
        rhs[end] = v2
    end

#    Bound = hcat(1.0, zeros(1,hsize-1))  #boundary conditions
#    Bound = [Bound; hcat(zeros(1,hsize-1), 1.0)]

    n1r = 0
    n2r = 0
    n1c = 0
    n2c = 0

    for (c,(d,pts)) in enumerate(zip(DL,VL))

        #        ptsmat = repeat(pts', size(d)[1], 1)

        ptsmat = repeat(pts,1, size(d)[1])

#        println("A2 ", A2.(pts))
        D2 = (A2.(ptsmat)) .*(d*d)
        D1 = (A1.(ptsmat)) .*(d)

        #D2 = A2(1.0)*(d*d)
        #D1 = A1(1.0)*(d)

        println("A0 ", A0.(pts[1:3]))

        v = diagm(A0.(pts))
        #v = A0.(ptsmat)

        N = size(d)[1] - 1

        x = getCpts(N)
        w = baryW1(N)
        y = getCpts1(N-2)
        P = barymat(y, x, w)

        H0 = v + D2 + D1

        n2r += size(v)[1] - 2
        n2c += size(v)[1]
        H[(n1r+1):n2r,(n1c+1):n2c] += P*H0
        #        S[(n1r+1):n2r,(n1c+1):n2c] += P*I(N+1)
        S[(n1r+1):n2r,(n1c+1):n2c] += P*diagm(B.(pts))

        rhs[(n1r+1):n2r] = P*B.(pts)

        if c < num_parts #add continuity of psi, psi'
            Bound = [hcat(zeros(1,n2c-1), 1, -1, zeros( 1,hsize-(n2c+1))); Bound] #make interface psi continuous
            st = size(DL[2])[2]
            Bound = [vcat(zeros(n1c,1), DL[c][end,:], -DL[c+1][1,:], zeros(hsize-st-n2c ,1))'; Bound] #make deriviates contin
        end

        n2r = n2r
        n1r = n2r
        n2c = n2c
        n1c = n2c

    end

    sb = size(Bound)[1]
    H[end-sb+1:end, :] = Bound
    S[end-sb+1:end, :] .= 0.0

    return H, S, rhs
end

function setup_fn_firstorder(; N = 20, bc1 = [:a,0,0.0], ranges=zeros(0,2), contin = Bool[], A0::Function, A1::Function, B::Function, a = -1.0, b = 1.0)

    println("solve first order")

#    ranges = ranges[contin,:]

    if typeof(N) <: Number
        N = ones(Int64, 1+size(ranges)[1]) * N
    end

    DL = []
    VL = []

    a1 = a
    b1 = a
    for c in 1:size(ranges)[1]
        b1 = ranges[c,1]
        pts, D = getmats(N[c]; a = a1+1e-15, b = b1-1e-15)
        push!(DL, D)
        push!(VL, pts)
        a1 = ranges[c,2]

    end
    b1 = b
    Vpts, Dab = getmats( N[end]; a = a1+1e-15, b = b1-1e-15)
    push!(DL, Dab)
    push!(VL, Vpts)

    num_parts = length(DL)
    hsize = 0
    for (c,(d,v)) in enumerate(zip(DL,VL))
        hsize += size(d)[1]
    end

    H = zeros(hsize, hsize)
    S = zeros(hsize, hsize)

    rhs = zeros(hsize)

    Bound = zeros(0,hsize)

    function apply_bc(bc)
        if ismissing(bc)
            return hcat(1.0, zeros(1,hsize-1)), 0.0
        end
        if bc[1] == :a && bc[2] == 0
            return hcat(1.0, zeros(1,hsize-1)), bc[3]
        elseif bc[1] == :b && bc[2] == 0
            return hcat(zeros(1,hsize-1), 1.0), bc[3]
        elseif bc[1] == :a && bc[2] == 1
            return hcat(DL[1][1,:]', zeros(1,hsize-size(DL[1])[1])), bc[3]
        elseif bc[1] == :b && bc[2] == 1
#            println("s ", size(zeros(1,hsize-size(DL[end])[1])))
#            println("s2 ",  size(DL[end][end,:]))
            return hcat( zeros(1,hsize-size(DL[end])[1]), DL[end][end,:]'),bc[3]
        else
            println("WARNING, bc not recognized")
            return hcat(1.0, zeros(1,hsize-1)), bc[3]
        end
    end

    println("bc1 $bc1 ")

    Bound = zeros(0,hsize)
    BC1, v1 = apply_bc(bc1)

#    println("hsize $hsize")
#    println("size BC1 ", size(BC1))
#    println("size BC2 ", size(BC2))

    Bound = [Bound; BC1]

    rhs[end] = v1

    n1r = 0
    n2r = 0
    n1c = 0
    n2c = 0

    for (c,(d,pts)) in enumerate(zip(DL,VL))

        ptsmat = repeat(pts,1, size(d)[1])

#        println("A2 ", A2.(pts))
        D1 = (A1.(ptsmat)) .*(d)

        #D2 = A2(1.0)*(d*d)
        #D1 = A1(1.0)*(d)

        v = diagm(A0.(pts))
        #v = A0.(ptsmat)

        N = size(d)[1] - 1

        x = getCpts(N)
        w = baryW1(N)
        y = getCpts1(N-1)
        P = barymat(y, x, w)

        H0 = v + D1

        n2r += size(v)[1] - 1
        n2c += size(v)[1]
        H[(n1r+1):n2r,(n1c+1):n2c] += P*H0
        S[(n1r+1):n2r,(n1c+1):n2c] += P*I(N+1)

        rhs[(n1r+1):n2r] = P*B.(pts)

        if c < num_parts #add continuity of psi, psi'
            Bound = [hcat(zeros(1,n2c-1), 1, -1, zeros( 1,hsize-(n2c+1))); Bound] #make interface psi continuous
            st = size(DL[2])[2]
            Bound = [vcat(zeros(n1c,1), DL[c][end,:], -DL[c+1][1,:], zeros(hsize-st-n2c ,1))'; Bound] #make deriviates contin
        end

        n2r = n2r
        n1r = n2r
        n2c = n2c
        n1c = n2c

    end

    sb = size(Bound)[1]
    H[end-sb+1:end, :] = Bound
    S[end-sb+1:end, :] .= 0.0

    return H, S, rhs
end

end #end module


text:
"""
    StudentTLikelihood(ν::T, σ::Real=one(T))

## Arguments
- `ν::Real` : degrees of freedom of the student-T
- `σ::Real` : standard deviation of the local scale

[Student-t likelihood](https://en.wikipedia.org/wiki/Student%27s_t-distribution) for regression:
```math
    p(y|f,ν,σ) = Γ(0.5(ν+1))/(sqrt(νπ) σ Γ(0.5ν)) * (1+(y-f)^2/(σ^2ν))^(-0.5(ν+1))
```
`ν` is the number of degrees of freedom and `σ` is the standard deviation for local scale of the data.

---

For the analytical solution, it is augmented via:
```math
    p(y|f,ω) = N(y|f,σ^2 ω)
```
Where `ω ~ IG(0.5ν,,0.5ν)` where `IG` is the inverse gamma distribution
See paper [Robust Gaussian Process Regression with a Student-t Likelihood](http://www.jmlr.org/papers/volume12/jylanki11a/jylanki11a.pdf)
"""
mutable struct StudentTLikelihood{T<:Real,A<:AbstractVector{T}} <: RegressionLikelihood{T}
    ν::T
    α::T
    σ::T
    c::A
    θ::A
    function StudentTLikelihood{T}(ν::T, σ::T=one(T)) where {T<:Real}
        ν > 0.5 || error("ν should be greater than 0.5")
        return new{T,Vector{T}}(ν, (ν + one(T)) / 2.0, σ)
    end
    function StudentTLikelihood{T}(
        ν::T, σ::T, c::A, θ::A
    ) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(ν, (ν + one(T)) / 2.0, σ, c, θ)
    end
end

function StudentTLikelihood(ν::T, σ::T=one(T)) where {T<:Real}
    return StudentTLikelihood{T}(ν, σ)
end

function implemented(
    ::StudentTLikelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function init_likelihood(
    likelihood::StudentTLikelihood{T},
    inference::AbstractInference{T},
    ::Int,
    nSamplesUsed::Int,
) where {T}
    if inference isa AnalyticVI || inference isa GibbsSampling
        StudentTLikelihood{T}(
            likelihood.ν, likelihood.σ, rand(T, nSamplesUsed), zeros(T, nSamplesUsed)
        )
    else
        StudentTLikelihood{T}(likelihood.ν, likelihood.σ)
    end
end

function (l::StudentTLikelihood)(y::Real, f::Real)
    return gamma(l.α) / (sqrt(l.ν * π) * gamma(l.ν / 2)) * (1 + abs2((y - f) / l.σ))^(-l.α)
    # tdistpdf(l.ν, (y - f) / l.σ) uses R so not differentiable
end

function Distributions.loglikelihood(l::StudentTLikelihood, y::Real, f::Real)
    return log(l(y, f))
    # tdistlogpdf(l.ν, (y - f) / l.σ) uses R so not differentiable
end

function Base.show(io::IO, model::StudentTLikelihood{T}) where {T}
    return print(io, "Student-t likelihood")
end

function compute_proba(
    l::StudentTLikelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    return μ, max.(σ², zero(σ²)) .+ 0.5 * l.ν * l.σ^2 / (0.5 * l.ν - 1)
end

## Local Updates ##

function local_updates!(
    l::StudentTLikelihood{T}, y::AbstractVector, μ::AbstractVector, diag_cov::AbstractVector
) where {T}
    l.c .= 0.5 * (diag_cov + abs2.(μ - y) .+ l.σ^2 * l.ν)
    return l.θ = l.α ./ l.c
end

function sample_local!(
    l::StudentTLikelihood{T}, y::AbstractVector, f::AbstractVector
) where {T}
    l.c .= rand.(InverseGamma.(l.α, 0.5 * (abs2.(f - y) .+ l.σ^2 * l.ν)))
    set_ω!(l, inv.(l.c))
    return nothing
end

## Global Gradients ##

@inline ∇E_μ(l::StudentTLikelihood, ::AOptimizer, y::AbstractVector) = (l.θ .* y,)
@inline ∇E_Σ(l::StudentTLikelihood, ::AOptimizer, ::AbstractVector) = (0.5 .* l.θ,)

## ELBO Section ##

function expec_loglikelihood(
    l::StudentTLikelihood{T},
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
) where {T}
    tot = -0.5 * length(y) * (log(twoπ * l.σ^2))
    tot += -sum(log.(l.c) .- digamma(l.α))
    tot +=
        -0.5 * (
            dot(l.θ, diag_cov) + dot(l.θ, abs2.(μ)) - 2.0 * dot(l.θ, μ .* y) +
            dot(l.θ, abs2.(y))
        )
    return tot
end

AugmentedKL(l::StudentTLikelihood, ::AbstractVector) = InverseGammaKL(l)

function InverseGammaKL(l::StudentTLikelihood{T}) where {T}
    α_p = l.ν / 2
    β_p = α_p * l.σ^2
    return InverseGammaKL(l.α, l.c, α_p, β_p)
end

## PDF and Log PDF Gradients ## (verified gradients)

function ∇loglikelihood(l::StudentTLikelihood{T}, y::Real, f::Real) where {T<:Real}
    return (one(T) + l.ν) * (y - f) / ((f - y)^2 + l.σ^2 * l.ν)
end

function hessloglikelihood(l::StudentTLikelihood{T}, y::Real, f::Real) where {T<:Real}
    v = l.ν * l.σ^2
    Δ² = (f - y)^2
    return (one(T) + l.ν) * (-v + Δ²) / (v + Δ²)^2
end


text:
module problem130

description = """
Composites with prime repunit property
A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k; for example, R(6) = 111111.
Given that n is a positive integer and GCD(n, 10) = 1, it can be shown that there always exists a value, k, for which R(k) is divisible by n,
  and let A(n) be the least such value of k; for example, A(7) = 6 and A(41) = 5.
You are given that for all primes, p > 5, that p − 1 is divisible by A(p). For example, when p = 41, A(41) = 5, and 40 is divisible by 5.
However, there are rare composite values for which this is also true; the first five examples being 91, 259, 451, 481, and 703.

Find the sum of the first twenty-five composite values of n for which GCD(n, 10) = 1 and n − 1 is divisible by A(n).
"""

using problem129

function search(x :: Int)
  results = Int[]
  for n = 2:Inf
    if gcd(n, 10) == 1 && !isprime(n) && (n-1) % problem129.A(n) == 0
      push!(results, n)
      if length(results) == x
        return sum(results)
      end
    end
  end
end

using Base.Test

@test search(5) == sum([91, 259, 451, 481, 703])

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions