text:
using Test
using Spandex

@testset "setindex! 1" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    a[5] = 'e'
    a[3] = 'c'

    @test 3 == a.nnz

    a[2] = 'b'
    a[4] = 'd'

    @test 5 == a.nnz
end

@testset "setindex! 2" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    @test 1 == a.nnz

    a[1] = 'a'
    @test 1 == a.nnz

    a[1] = 'b'
    @test 1 == a.nnz
end

@testset "contains 1" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    a[5] = 'e'
    a[3] = 'c'

    @test contains(a, 1)
    @test contains(a, 3)
    @test contains(a, 5)
    @test !contains(a, 2)
    @test !contains(a, 4)
    @test !contains(a, 6)

    a[2] = 'b'
    a[4] = 'd'

    @test contains(a, 1)
    @test contains(a, 2)
    @test contains(a, 3)
    @test contains(a, 4)
    @test contains(a, 5)
    @test !contains(a, 6)
end

@testset "contains 2" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    @test contains(a, 1)
    @test !contains(a, 9)

    a[1] = 'a'
    @test contains(a, 1)
    @test !contains(a, 9)

    a[1] = 'b'
    @test contains(a, 1)
    @test !contains(a, 9)
end

@testset "getindex 1" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    a[5] = 'e'
    a[3] = 'c'

    @test 'a' == a[1]
    @test 'c' == a[3]
    @test 'e' == a[5]

    a[2] = 'b'
    a[4] = 'd'

    @test 'a' == a[1]
    @test 'b' == a[2]
    @test 'c' == a[3]
    @test 'd' == a[4]
    @test 'e' == a[5]
end

@testset "getindex 2" begin
    local a = SparseArray{Char}(10)

    a[1] = 'a'
    @test 'a' == a[1]

    a[1] = 'a'
    @test 'a' == a[1]

    a[1] = 'b'
    @test 'b' == a[1]
end

@testset "equals 1" begin
    local a = SparseArray{Char}(10)
    a[1] = 'a'
    a[2] = 'b'
    a[3] = 'c'
    a[4] = 'd'
    a[5] = 'e'

    local b = SparseArray{Char}(10)
    b[1] = 'a'
    b[5] = 'e'
    b[3] = 'c'
    b[2] = 'b'
    b[4] = 'd'

    @test equals(a, b)
end

@testset "equals 2" begin
    local a = SparseArray{Char}(20)
    a[1] = 'a'
    a[2] = 'b'
    a[3] = 'c'
    a[4] = 'd'
    a[5] = 'e'

    local b = SparseArray{Char}(10)
    b[1] = 'a'
    b[5] = 'e'
    b[3] = 'c'
    b[2] = 'b'
    b[4] = 'd'

    @test !equals(b, a)
end

@testset "equals 3" begin
    local a = SparseArray{Char}(20)
    a[1] = 'a'

    local b = SparseArray{Char}(10)
    b[1] = 'b'

    @test !equals(a, b)
end


text:
function commandSave(ed::Editor, args::String)
    editorSave(ed, args)
end

addCommand(:save, commandSave,
           help="save: save the file")

ext: jl


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   This file contains helpers related to space indices.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

for sym in [:F10, :F10obs, :F10adj, :F10M, :F10Mobs, :F10Madj, :Kp, :Ap,
            :Kp_vect, :Ap_vect, :S10, :S81a, :M10, :M81a, :Y10, :Y81a, :DstΔTc]

    qsym = Meta.quot(sym)

    @eval begin
        export $sym
        $sym() = Val($qsym)
    end
end


text:
"""
**Checks that the queries for occurrences searches are well formatted**

This is used internally.

Everything this function does is derived from the GBIF API documentation,
including (and especially) the values for enum types. This modifies the
queryset. Filters that are not allowed are removed, and filters that have
incorrect values are dropped too.

This feels like the most conservative option -- the user can always filter the
results when they are returned.
"""
function check_records_parameters!(q::Dict)

  # List of fields from GBIF
  allowed_fields = "q", "basisOfRecord", "catalogNumber", "collectionCode",
    "continent", "country", "datasetKey", "decimalLatitude", "decimalLongitude",
    "depth", "elevation", "eventDate", "geometry", "hasCoordinate",
    "hasGeospatialIssue", "institutionCode", "issue", "lastInterpreted",
    "mediaType", "month", "occurrenceId", "organismId", "protocol", "license",
    "publishingCountry", "publishingOrg", "crawlId", "recordedBy", "recordNumber",
    "scientificName", "locality", "stateProvince", "waterBody", "taxonKey",
    "kingdomKey", "phylumKey", "classKey", "orderKey", "familyKey", "genusKey",
    "subGenusKey", "speciesKey", "year", "establishmentMeans", "repatriated",
    "typeStatus", "facet", "facetMincount", "facetMultiselect", "limit", "offset"

  for (k, v) in q
    if !(k ∈ allowed_fields)
      @warn "$(k) is not a supported field -- will be dropped from the queryset"
      delete!(q, k)
    end
  end

  # Country must be a two-letters country code
  if "country" ∈ keys(q)
    if length(q["country"]) != 2
      @warn "$(q["country"]) is not a two letter country code -- will be dropped from the queryset"
      delete!(q, "country")
    end
  end

  # Latitude and longitudes
  # TODO lat -90/90 lon -180/180, can be "min,max"

  # ENUMs
  for (k, v) in q
    if k ∈ keys(gbifenums)
      okvals = filter(x -> x ∈ gbifenums[k], v)
      if length(okvals) != length(v)
        @warn "Some values in $(k) were invalid -- will be dropped from the queryset"
      end
      q[k] = okvals
    end
  end

end


text:
"""
    EarthEnv{Union{HabitatHeterogeneity,LandCover}} <: RasterDataSource

Data from the `EarthEnv` including `HabitatHeterogeneity` and `LandCover`

See: [www.earthenv.org](http://www.earthenv.org/)
"""
struct EarthEnv{X} <: RasterDataSource end

rasterpath(::Type{EarthEnv}) = joinpath(rasterpath(), "EarthEnv")

rasterurl(::Type{EarthEnv}) = URI(scheme="https", host="data.earthenv.org")

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
function parse_all(doc, server)
    ps = CSTParser.ParseState(doc._content)
    if endswith(doc._uri, ".jmd")
        doc.code.cst, ps = parse_jmd(ps, doc._content)
    else
        doc.code.cst, ps = CSTParser.parse(ps, true)
    end
    # update_includes(doc, server)
    empty!(doc.diagnostics)
    ls_diags = []

    for err in ps.errors
        if err.description == "Expected end."
            rng2 = max(0, first(err.loc)-1):last(err.loc)
            stack, offsets = StaticLint.get_stack(doc.code.cst, first(rng2))
            for i = length(stack):-1:1
                if stack[i] isa CSTParser.EXPR{T} where T <: Union{CSTParser.Begin,CSTParser.Quote,CSTParser.ModuleH,CSTParser.Function,CSTParser.Macro,CSTParser.For,CSTParser.While,CSTParser.If} && last(stack[i].args) isa CSTParser.EXPR{CSTParser.ErrorToken} && stack[i].args[end].args[1] isa CSTParser.KEYWORD
                    rng1 = offsets[i] .+ stack[i].args[1].span
                    push!(ls_diags, Diagnostic(Range(doc, rng1), 1, "ERROR", "ERROR", "Closing end is missing.", nothing))
                end
            end
            push!(ls_diags, Diagnostic(Range(doc, rng2), 1, "ERROR", "ERROR", err.description, nothing))
        else
            rng = max(0, first(err.loc)-1):last(err.loc)
            push!(ls_diags, Diagnostic(Range(doc, rng), 1, "ERROR", "ERROR", err.description, nothing))
        end
    end

    if server.runlinter
        if doc._runlinter
            StaticLint.pass(doc.code)
            bindings = StaticLint.build_bindings(server, find_root(doc, server).code);
            empty!(doc.code.rref)
            empty!(doc.code.uref)
            StaticLint.resolve_refs(doc.code.state.refs, bindings, doc.code.rref, doc.code.uref);

            for (i, r) in enumerate(doc.code.uref)
                r isa StaticLint.Reference{CSTParser.BinarySyntaxOpCall} && continue
                rng = 0:last(r.val.span)
                push!(ls_diags ,Diagnostic(Range(doc, r.loc.offset .+ rng), 2, "missing variable", "missing variable", "$i missing variable: $(string(Expr(r.val))) at $(r.loc)", nothing))
            end
        end
    end
    send(JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, ls_diags)), server)
end

StaticLint.getfile(server::LanguageServerInstance, path) = server.documents[URI2(filepath2uri(path))].code
StaticLint.setfile(server::LanguageServerInstance, path, x) = server.documents[URI2(filepath2uri(path))] = x
StaticLint.is_loaded(server::LanguageServerInstance, path) = haskey(server.documents, URI2(filepath2uri(path)))
function StaticLint.load_file(server::LanguageServerInstance, path::String, index, nb, parent)
    code = read(path, String)
    uri = filepath2uri(path)
    doc = Document(uri, code, true, server, index, nb, parent)
    StaticLint.setfile(server, path, doc)
    return doc.code
end

function convert_diagnostic(h::LSDiagnostic{T}, doc::Document) where {T}
    rng = Range(doc, h.loc)
    code = 1
    Diagnostic(rng, code, string(T), string(typeof(h).name), isempty(h.message) ? string(T) : h.message, nothing)
end

function publish_diagnostics(doc::Document, server)
    ls_diags = map(doc.diagnostics) do diag
        convert_diagnostic(diag, doc)
    end
    publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, ls_diags)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, publishDiagnosticsParams)
    send(response, server)
end

function clear_diagnostics(uri::URI2, server)
    doc = server.documents[uri]
    empty!(doc.diagnostics)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, Diagnostic[]))
    send(response, server)
end

function clear_diagnostics(server)
    for (uri, doc) in server.documents
        clear_diagnostics(uri, server)
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((d+e*x)^4*(a+b*arctanh(c*x)), x) == :(1//5*e^-1*(d+e*x)^5*(a+b*arctanh(c*x))+1//20*b*c^-1*e^4*x^4+-1//10*b*c^-5*e^-1*(-1*e+c*d)^5*log(1+c*x)+1//3*b*d*c^-1*e^3*x^3+1//10*b*c^-5*e^-1*(e+c*d)^5*log(1+-1*c*x)+1//10*b*c^-3*e^2*x^2*(e^2+10*c^2*d^2)+b*d*e*x*c^-3*(e^2+2*c^2*d^2))
@test integrate((d+e*x)^3*(a+b*arctanh(c*x)), x) == :(1//4*e^-1*(d+e*x)^4*(a+b*arctanh(c*x))+1//12*b*c^-1*e^3*x^3+(1/2)*b*d*c^-1*e^2*x^2+-1//8*b*c^-4*e^-1*(-1*e+c*d)^4*log(1+c*x)+1//4*b*e*x*c^-3*(e^2+6*c^2*d^2)+1//8*b*c^-4*e^-1*(e+c*d)^4*log(1+-1*c*x))
@test integrate((d+e*x)^2*(a+b*arctanh(c*x)), x) == :(1//3*e^-1*(d+e*x)^3*(a+b*arctanh(c*x))+1//6*b*c^-1*e^2*x^2+b*d*e*x*c^-1+-1//6*b*c^-3*e^-1*(-1*e+c*d)^3*log(1+c*x)+1//6*b*c^-3*e^-1*(e+c*d)^3*log(1+-1*c*x))
@test integrate((a+b*arctanh(c*x))*(d+e*x), x) == :((1/2)*e^-1*(d+e*x)^2*(a+b*arctanh(c*x))+(1/2)*b*e*x*c^-1+-1//4*b*c^-2*e^-1*(-1*e+c*d)^2*log(1+c*x)+1//4*b*c^-2*e^-1*(e+c*d)^2*log(1+-1*c*x))
@test integrate((d+e*x)^-1*(a+b*arctanh(c*x)), x) == :(e^-1*(a+b*arctanh(c*x))*log(2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+(1/2)*b*e^-1*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+-1*e^-1*(a+b*arctanh(c*x))*log(2*(1+c*x)^-1)+-1//2*b*e^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x)))
@test integrate((d+e*x)^-2*(a+b*arctanh(c*x)), x) == :(-1*e^-1*(d+e*x)^-1*(a+b*arctanh(c*x))+-1*b*c*(-1*e^2+c^2*d^2)^-1*log(d+e*x)+b*c*e^-1*(-2*e+2*c*d)^-1*log(1+c*x)+-1//2*b*c*e^-1*(e+c*d)^-1*log(1+-1*c*x))
@test integrate((d+e*x)^-3*(a+b*arctanh(c*x)), x) == :(-1//2*e^-1*(d+e*x)^-2*(a+b*arctanh(c*x))+b*c*(d+e*x)^-1*(-2*e^2+2*c^2*d^2)^-1+-1*b*d*c^3*(-1*e^2+c^2*d^2)^-2*log(d+e*x)+-1//4*b*c^2*e^-1*(e+c*d)^-2*log(1+-1*c*x)+1//4*b*c^2*e^-1*(-1*e+c*d)^-2*log(1+c*x))
@test integrate((a+b*arctanh(c*x))^2*(d+e*x)^3, x) == :(1//4*e^-1*(a+b*arctanh(c*x))^2*(d+e*x)^4+d*c^-3*(a+b*arctanh(c*x))^2*(e^2+c^2*d^2)+-1//4*c^-4*e^-1*(a+b*arctanh(c*x))^2*(e^4+c^4*d^4+6*c^2*d^2*e^2)+1//12*b^2*c^-4*e^3*log(1+-1*c^2*x^2)+1//12*b^2*c^-2*e^3*x^2+d*x*b^2*c^-2*e^2+-1*d*b^2*c^-3*e^2*arctanh(c*x)+-1*d*b^2*c^-3*(e^2+c^2*d^2)*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+1//4*e*b^2*c^-4*(e^2+6*c^2*d^2)*log(1+-1*c^2*x^2)+1//6*b*c^-1*e^3*x^3*(a+b*arctanh(c*x))+b*d*c^-1*e^2*x^2*(a+b*arctanh(c*x))+(1/2)*a*b*e*x*c^-3*(e^2+6*c^2*d^2)+(1/2)*e*x*b^2*c^-3*(e^2+6*c^2*d^2)*arctanh(c*x)+-2*b*d*c^-3*(a+b*arctanh(c*x))*(e^2+c^2*d^2)*log(2*(1+-1*c*x)^-1))
@test integrate((a+b*arctanh(c*x))^2*(d+e*x)^2, x) == :(1//3*c^-3*(a+b*arctanh(c*x))^2*(e^2+3*c^2*d^2)+1//3*e^-1*(a+b*arctanh(c*x))^2*(d+e*x)^3+-1//3*d*e^-1*(a+b*arctanh(c*x))^2*(d^2+3*c^-2*e^2)+-1//3*b^2*c^-3*e^2*arctanh(c*x)+-1//3*b^2*c^-3*(e^2+3*c^2*d^2)*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+1//3*x*b^2*c^-2*e^2+d*e*b^2*c^-2*log(1+-1*c^2*x^2)+-2//3*b*c^-3*(a+b*arctanh(c*x))*(e^2+3*c^2*d^2)*log(2*(1+-1*c*x)^-1)+1//3*b*c^-1*e^2*x^2*(a+b*arctanh(c*x))+2*a*b*d*e*x*c^-1+2*d*e*x*b^2*c^-1*arctanh(c*x))
@test integrate((a+b*arctanh(c*x))^2*(d+e*x), x) == :(d*c^-1*(a+b*arctanh(c*x))^2+(1/2)*e^-1*(a+b*arctanh(c*x))^2*(d+e*x)^2+-1//2*e^-1*(a+b*arctanh(c*x))^2*(d^2+c^-2*e^2)+(1/2)*e*b^2*c^-2*log(1+-1*c^2*x^2)+-1*d*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+a*b*e*x*c^-1+e*x*b^2*c^-1*arctanh(c*x)+-2*b*d*c^-1*(a+b*arctanh(c*x))*log(2*(1+-1*c*x)^-1))
@test integrate((a+b*arctanh(c*x))^2*(d+e*x)^-1, x) == :(e^-1*(a+b*arctanh(c*x))^2*log(2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+(1/2)*b^2*e^-1*Polylogarithms.polylog(3,1+-2*(1+c*x)^-1)+-1*e^-1*(a+b*arctanh(c*x))^2*log(2*(1+c*x)^-1)+-1//2*b^2*e^-1*Polylogarithms.polylog(3,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+b*e^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+-1*b*e^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x)))
@test integrate((a+b*arctanh(c*x))^2*(d+e*x)^-2, x) == :(-1*e^-1*(a+b*arctanh(c*x))^2*(d+e*x)^-1+c*b^2*(-1*e^2+c^2*d^2)^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+-1*c*b^2*(-1*e^2+c^2*d^2)^-1*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+c*b^2*e^-1*(-2*e+2*c*d)^-1*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+(1/2)*c*b^2*e^-1*(e+c*d)^-1*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+-2*b*c*(-1*e^2+c^2*d^2)^-1*(a+b*arctanh(c*x))*log(2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+2*b*c*(-1*e^2+c^2*d^2)^-1*(a+b*arctanh(c*x))*log(2*(1+c*x)^-1)+b*c*e^-1*(e+c*d)^-1*(a+b*arctanh(c*x))*log(2*(1+-1*c*x)^-1)+-1*b*c*e^-1*(-1*e+c*d)^-1*(a+b*arctanh(c*x))*log(2*(1+c*x)^-1))
@test integrate((a+b*arctanh(c*x))^3*(d+e*x)^3, x) == :(1//4*e^-1*(a+b*arctanh(c*x))^3*(d+e*x)^4+d*c^-3*(a+b*arctanh(c*x))^3*(e^2+c^2*d^2)+-1//4*b^3*c^-4*e^3*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+-1//4*b^3*c^-4*e^3*arctanh(c*x)+-1//4*c^-4*e^-1*(a+b*arctanh(c*x))^3*(e^4+c^4*d^4+6*c^2*d^2*e^2)+1//4*b*c^-4*e^3*(a+b*arctanh(c*x))^2+1//4*x*b^3*c^-3*e^3+-3//2*b*d*c^-3*e^2*(a+b*arctanh(c*x))^2+-3//4*e*b^3*c^-4*(e^2+6*c^2*d^2)*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+-1//2*b^2*c^-4*e^3*(a+b*arctanh(c*x))*log(2*(1+-1*c*x)^-1)+1//4*b*c^-1*e^3*x^3*(a+b*arctanh(c*x))^2+1//4*b^2*c^-2*e^3*x^2*(a+b*arctanh(c*x))+3//2*d*b^3*c^-3*e^2*log(1+-1*c^2*x^2)+3//2*d*b^3*c^-3*(e^2+c^2*d^2)*Polylogarithms.polylog(3,1+-2*(1+-1*c*x)^-1)+3//4*b*e*c^-4*(a+b*arctanh(c*x))^2*(e^2+6*c^2*d^2)+-3*b*d*c^-3*(a+b*arctanh(c*x))^2*(e^2+c^2*d^2)*log(2*(1+-1*c*x)^-1)+-3*d*b^2*c^-3*(a+b*arctanh(c*x))*(e^2+c^2*d^2)*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+3*a*d*x*b^2*c^-2*e^2+3*d*x*b^3*c^-2*e^2*arctanh(c*x)+-3//2*e*b^2*c^-4*(a+b*arctanh(c*x))*(e^2+6*c^2*d^2)*log(2*(1+-1*c*x)^-1)+3//2*b*d*c^-1*e^2*x^2*(a+b*arctanh(c*x))^2+3//4*b*e*x*c^-3*(a+b*arctanh(c*x))^2*(e^2+6*c^2*d^2))
@test integrate((a+b*arctanh(c*x))^3*(d+e*x)^2, x) == :(1//3*c^-3*(a+b*arctanh(c*x))^3*(e^2+3*c^2*d^2)+1//3*e^-1*(a+b*arctanh(c*x))^3*(d+e*x)^3+(1/2)*b^3*c^-3*e^2*log(1+-1*c^2*x^2)+(1/2)*b^3*c^-3*(e^2+3*c^2*d^2)*Polylogarithms.polylog(3,1+-2*(1+-1*c*x)^-1)+-1//2*b*c^-3*e^2*(a+b*arctanh(c*x))^2+-1//3*d*e^-1*(a+b*arctanh(c*x))^3*(d^2+3*c^-2*e^2)+a*x*b^2*c^-2*e^2+x*b^3*c^-2*e^2*arctanh(c*x)+(1/2)*b*c^-1*e^2*x^2*(a+b*arctanh(c*x))^2+-1*b*c^-3*(a+b*arctanh(c*x))^2*(e^2+3*c^2*d^2)*log(2*(1+-1*c*x)^-1)+-1*b^2*c^-3*(a+b*arctanh(c*x))*(e^2+3*c^2*d^2)*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+-3*d*e*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+3*b*d*e*c^-2*(a+b*arctanh(c*x))^2+-6*d*e*b^2*c^-2*(a+b*arctanh(c*x))*log(2*(1+-1*c*x)^-1)+3*b*d*e*x*c^-1*(a+b*arctanh(c*x))^2)
@test integrate((a+b*arctanh(c*x))^3*(d+e*x), x) == :(d*c^-1*(a+b*arctanh(c*x))^3+(1/2)*e^-1*(a+b*arctanh(c*x))^3*(d+e*x)^2+-1//2*e^-1*(a+b*arctanh(c*x))^3*(d^2+c^-2*e^2)+-3//2*e*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+3//2*b*e*c^-2*(a+b*arctanh(c*x))^2+3//2*d*b^3*c^-1*Polylogarithms.polylog(3,1+-2*(1+-1*c*x)^-1)+-3*b*d*c^-1*(a+b*arctanh(c*x))^2*log(2*(1+-1*c*x)^-1)+-3*d*b^2*c^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1)+-3*e*b^2*c^-2*(a+b*arctanh(c*x))*log(2*(1+-1*c*x)^-1)+3//2*b*e*x*c^-1*(a+b*arctanh(c*x))^2)
@test integrate((a+b*arctanh(c*x))^3*(d+e*x)^-1, x) == :(e^-1*(a+b*arctanh(c*x))^3*log(2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+-1*e^-1*(a+b*arctanh(c*x))^3*log(2*(1+c*x)^-1)+-3//4*b^3*e^-1*Polylogarithms.polylog(4,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+3//4*b^3*e^-1*Polylogarithms.polylog(4,1+-2*(1+c*x)^-1)+-3//2*b*e^-1*(a+b*arctanh(c*x))^2*Polylogarithms.polylog(2,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+-3//2*b^2*e^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(3,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+3//2*b*e^-1*(a+b*arctanh(c*x))^2*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+3//2*b^2*e^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(3,1+-2*(1+c*x)^-1))
@test integrate((a+b*arctanh(c*x))^3*(d+e*x)^-2, x) == :(-1*e^-1*(a+b*arctanh(c*x))^3*(d+e*x)^-1+-3*c*b^3*(-2*e^2+2*c^2*d^2)^-1*Polylogarithms.polylog(3,1+-2*(1+c*x)^-1)+3*c*b^3*(-2*e^2+2*c^2*d^2)^-1*Polylogarithms.polylog(3,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+-3*b*c*(a+b*arctanh(c*x))^2*(-1*e^2+c^2*d^2)^-1*log(2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+-3*c*b^2*(-1*e^2+c^2*d^2)^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+3*b*c*(a+b*arctanh(c*x))^2*(-1*e^2+c^2*d^2)^-1*log(2*(1+c*x)^-1)+3*c*b^2*(-1*e^2+c^2*d^2)^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*c*(1+c*x)^-1*(e+c*d)^-1*(d+e*x))+3*c*b^3*e^-1*(-4*e+4*c*d)^-1*Polylogarithms.polylog(3,1+-2*(1+c*x)^-1)+-3//4*c*b^3*e^-1*(e+c*d)^-1*Polylogarithms.polylog(3,1+-2*(1+-1*c*x)^-1)+-3*b*c*e^-1*(a+b*arctanh(c*x))^2*(-2*e+2*c*d)^-1*log(2*(1+c*x)^-1)+3*c*b^2*e^-1*(-2*e+2*c*d)^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*(1+c*x)^-1)+3//2*b*c*e^-1*(a+b*arctanh(c*x))^2*(e+c*d)^-1*log(2*(1+-1*c*x)^-1)+3//2*c*b^2*e^-1*(e+c*d)^-1*(a+b*arctanh(c*x))*Polylogarithms.polylog(2,1+-2*(1+-1*c*x)^-1))
@test integrate((d+e*x)^3*(a+b*arctanh(c*x^2)), x) == :(1//4*e^-1*(d+e*x)^4*(a+b*arctanh(c*x^2))+1//4*b*c^-1*e^3*x^2+b*d*c^-3//2*(-1*e^2+c*d^2)*arctan(x*c^(1/2))+-1*b*d*c^-3//2*(e^2+c*d^2)*arctanh(x*c^(1/2))+2*b*d*x*c^-1*e^2+-1//8*b*c^-2*e^-1*(e^4+c^2*d^4+-6*c*d^2*e^2)*log(1+c*x^2)+1//8*b*c^-2*e^-1*(e^4+c^2*d^4+6*c*d^2*e^2)*log(1+-1*c*x^2))
@test integrate((d+e*x)^2*(a+b*arctanh(c*x^2)), x) == :(1//3*e^-1*(d+e*x)^3*(a+b*arctanh(c*x^2))+-1//3*b*c^-3//2*(e^2+3*c*d^2)*arctanh(x*c^(1/2))+1//3*b*c^-3//2*(-1*e^2+3*c*d^2)*arctan(x*c^(1/2))+2//3*b*x*c^-1*e^2+-1//6*b*d*c^-1*e^-1*(-3*e^2+c*d^2)*log(1+c*x^2)+1//6*b*d*c^-1*e^-1*(3*e^2+c*d^2)*log(1+-1*c*x^2))
@test integrate((a+b*arctanh(c*x^2))*(d+e*x), x) == :((1/2)*e^-1*(d+e*x)^2*(a+b*arctanh(c*x^2))+b*d*c^-1//2*arctan(x*c^(1/2))+-1*b*d*c^-1//2*arctanh(x*c^(1/2))+-1//4*b*c^-1*e^-1*(-1*e^2+c*d^2)*log(1+c*x^2)+1//4*b*c^-1*e^-1*(e^2+c*d^2)*log(1+-1*c*x^2))
@test integrate((d+e*x)^-1*(a+b*arctanh(c*x^2)), x) == :(e^-1*(a+b*arctanh(c*x^2))*log(d+e*x)+(1/2)*b*e^-1*Polylogarithms.polylog(2,c^(1/2)*(e+d*c^(1/2))^-1*(d+e*x))+(1/2)*b*e^-1*Polylogarithms.polylog(2,c^(1/2)*(-1*e+d*c^(1/2))^-1*(d+e*x))+-1//2*b*e^-1*Polylogarithms.polylog(2,(-1c)^(1/2)*(e+d*(-1c)^(1/2))^-1*(d+e*x))+-1//2*b*e^-1*Polylogarithms.polylog(2,(-1c)^(1/2)*(-1*e+d*(-1c)^(1/2))^-1*(d+e*x))+(1/2)*b*e^-1*log(e*(e+d*c^(1/2))^-1*(1+-1*x*c^(1/2)))*log(d+e*x)+(1/2)*b*e^-1*log(-1*e*(-1*e+d*c^(1/2))^-1*(1+x*c^(1/2)))*log(d+e*x)+-1//2*b*e^-1*log(e*(e+d*(-1c)^(1/2))^-1*(1+-1*x*(-1c)^(1/2)))*log(d+e*x)+-1//2*b*e^-1*log(-1*e*(-1*e+d*(-1c)^(1/2))^-1*(1+x*(-1c)^(1/2)))*log(d+e*x))
@test integrate((d+e*x)^-2*(a+b*arctanh(c*x^2)), x) == :(-1*e^-1*(d+e*x)^-1*(a+b*arctanh(c*x^2))+b*c^(1/2)*(e^2+c*d^2)^-1*arctan(x*c^(1/2))+-1*b*c^(1/2)*(-1*e^2+c*d^2)^-1*arctanh(x*c^(1/2))+(1/2)*b*c*d*e^-1*(e^2+c*d^2)^-1*log(1+c*x^2)+2*b*c*d*e*(-1*e^4+c^2*d^4)^-1*log(d+e*x)+-1//2*b*c*d*e^-1*(-1*e^2+c*d^2)^-1*log(1+-1*c*x^2))
@test integrate((d+e*x)^-3*(a+b*arctanh(c*x^2)), x) == :(-1//2*e^-1*(d+e*x)^-2*(a+b*arctanh(c*x^2))+b*d*c^3//2*(e^2+c*d^2)^-2*arctan(x*c^(1/2))+-1*b*d*c^3//2*(-1*e^2+c*d^2)^-2*arctanh(x*c^(1/2))+b*c*e*(-1*e^4+c^2*d^4)^-2*(e^4+3*c^2*d^4)*log(d+e*x)+-1*b*c*d*e*(d+e*x)^-1*(-1*e^4+c^2*d^4)^-1+-1//4*b*c*e^-1*(-1*e^2+c*d^2)^-2*(e^2+c*d^2)*log(1+-1*c*x^2)+1//4*b*c*e^-1*(e^2+c*d^2)^-2*(-1*e^2+c*d^2)*log(1+c*x^2))
@test integrate((a+b*arctanh(c*x^2))^2*(d+e*x), x) == :(d*x*a^2+(1/2)*e*c^-1*(a+b*arctanh(c*x^2))^2+(1/2)*e*x^2*(a+b*arctanh(c*x^2))^2+d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+x*c^(1/2))^-1)+d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*c^(1/2))^-1)+-1*d*b^2*c^-1//2*arctanh(x*c^(1/2))^2+-1//2*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*c^(1/2)*(1+x*c^(1/2))^-1*(c^(1/2)+(-1c)^(1/2))^-1*(1+x*(-1c)^(1/2)))+-1//2*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+2*c^(1/2)*(1+x*c^(1/2))^-1*((-1c)^(1/2)+-1*c^(1/2))^-1*(1+-1*x*(-1c)^(1/2)))+-1//2*e*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+-1*c*x^2)^-1)+1//4*d*x*b^2*log(1+c*x^2)^2+1//4*d*x*b^2*log(1+-1*c*x^2)^2+im*d*b^2*c^-1//2*arctan(x*c^(1/2))^2+im*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+im*x*c^(1/2))^-1)+im*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+-1*im*x*c^(1/2))^-1)+a*b*d*x*log(1+c*x^2)+d*b^2*c^-1//2*arctan(x*c^(1/2))*log((1+-1*im*x*c^(1/2))^-1*(1+im)*(1+-1*x*c^(1/2)))+d*b^2*c^-1//2*arctan(x*c^(1/2))*log((1+-1*im*x*c^(1/2))^-1*(1+-1im)*(1+x*c^(1/2)))+d*b^2*c^-1//2*arctan(x*c^(1/2))*log(1+c*x^2)+d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(-2*c^(1/2)*(1+x*c^(1/2))^-1*((-1c)^(1/2)+-1*c^(1/2))^-1*(1+-1*x*(-1c)^(1/2)))+d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(2*c^(1/2)*(1+x*c^(1/2))^-1*(c^(1/2)+(-1c)^(1/2))^-1*(1+x*(-1c)^(1/2)))+d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(1+-1*c*x^2)+-1*a*b*d*x*log(1+-1*c*x^2)+-1*b*e*c^-1*(a+b*arctanh(c*x^2))*log(2*(1+-1*c*x^2)^-1)+-1*d*b^2*c^-1//2*arctan(x*c^(1/2))*log(1+-1*c*x^2)+-1*d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(1+c*x^2)+-2*a*b*d*c^-1//2*arctanh(x*c^(1/2))+-2*d*b^2*c^-1//2*arctan(x*c^(1/2))*log(2*(1+-1*im*x*c^(1/2))^-1)+-2*d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(2*(1+x*c^(1/2))^-1)+2*a*b*d*c^-1//2*arctan(x*c^(1/2))+2*d*b^2*c^-1//2*arctan(x*c^(1/2))*log(2*(1+im*x*c^(1/2))^-1)+2*d*b^2*c^-1//2*arctanh(x*c^(1/2))*log(2*(1+-1*x*c^(1/2))^-1)+-1//2*im*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-1*(1+-1*im*x*c^(1/2))^-1*(1+im)*(1+-1*x*c^(1/2)))+-1//2*im*d*b^2*c^-1//2*Polylogarithms.polylog(2,1+-1*(1+-1*im*x*c^(1/2))^-1*(1+-1im)*(1+x*c^(1/2)))+-1//2*d*x*b^2*log(1+c*x^2)*log(1+-1*c*x^2))
@test integrate((d+e*x)^2*(a+b*arctanh(c*x^3)), x) == :(1//3*e^-1*(d+e*x)^3*(a+b*arctanh(c*x^3))+(1/2)*b*c^-1//3*d^2*log(1+-1*c^2//3*x^2)+-1//4*b*c^-1//3*d^2*log(1+c^2//3*x^2+c^4//3*x^4)+(1/2)*b*3^(1/2)*c^-1//3*d^2*arctan(1//3*3^(1/2)*(1+2*c^2//3*x^2))+-1*b*d*e*c^-2//3*arctanh(x*c^1//3)+-1//4*b*d*e*c^-2//3*log(1+x*c^1//3+c^2//3*x^2)+-1//6*b*c^-1*e^-1*(-1*e^3+c*d^3)*log(1+c*x^3)+1//4*b*d*e*c^-2//3*log(1+c^2//3*x^2+-1*x*c^1//3)+1//6*b*c^-1*e^-1*(e^3+c*d^3)*log(1+-1*c*x^3)+(1/2)*b*d*e*3^(1/2)*c^-2//3*arctan(1//3*3^(1/2)+2//3*x*3^(1/2)*c^1//3)+-1//2*b*d*e*3^(1/2)*c^-2//3*arctan(1//3*3^(1/2)+-2//3*x*3^(1/2)*c^1//3))
@test integrate((a+b*arctanh(c*x^3))*(d+e*x), x) == :((1/2)*e^-1*(d+e*x)^2*(a+b*arctanh(c*x^3))+(1/2)*b*d*c^-1//3*log(1+-1*c^2//3*x^2)+-1//2*b*e*c^-2//3*arctanh(x*c^1//3)+-1//2*b*d^2*e^-1*arctanh(c*x^3)+-1//4*b*d*c^-1//3*log(1+c^2//3*x^2+c^4//3*x^4)+-1//8*b*e*c^-2//3*log(1+x*c^1//3+c^2//3*x^2)+1//8*b*e*c^-2//3*log(1+c^2//3*x^2+-1*x*c^1//3)+(1/2)*b*d*3^(1/2)*c^-1//3*arctan(1//3*3^(1/2)*(1+2*c^2//3*x^2))+-1//4*b*e*3^(1/2)*c^-2//3*arctan(1//3*3^(1/2)+-2//3*x*3^(1/2)*c^1//3)+1//4*b*e*3^(1/2)*c^-2//3*arctan(1//3*3^(1/2)+2//3*x*3^(1/2)*c^1//3))
@test integrate((d+e*x)^-1*(a+b*arctanh(c*x^3)), x) == :(e^-1*(a+b*arctanh(c*x^3))*log(d+e*x)+(1/2)*b*e^-1*Polylogarithms.polylog(2,c^1//3*(e+d*c^1//3)^-1*(d+e*x))+(1/2)*b*e^-1*Polylogarithms.polylog(2,c^1//3*(d*c^1//3+e*(-1)^2//3)^-1*(d+e*x))+(1/2)*b*e^-1*Polylogarithms.polylog(2,c^1//3*(d*c^1//3+-1*e*(-1)^1//3)^-1*(d+e*x))+-1//2*b*e^-1*Polylogarithms.polylog(2,c^1//3*(-1*e+d*c^1//3)^-1*(d+e*x))+-1//2*b*e^-1*Polylogarithms.polylog(2,c^1//3*(d*c^1//3+e*(-1)^1//3)^-1*(d+e*x))+-1//2*b*e^-1*Polylogarithms.polylog(2,c^1//3*(d*c^1//3+-1*e*(-1)^2//3)^-1*(d+e*x))+(1/2)*b*e^-1*log(e*(e+d*c^1//3)^-1*(1+-1*x*c^1//3))*log(d+e*x)+(1/2)*b*e^-1*log(-1*e*(d*c^1//3+-1*e*(-1)^1//3)^-1*((-1)^1//3+x*c^1//3))*log(d+e*x)+(1/2)*b*e^-1*log(e*(-1)^2//3*(d*c^1//3+e*(-1)^2//3)^-1*(1+x*(-1)^1//3*c^1//3))*log(d+e*x)+-1//2*b*e^-1*log(-1*e*(-1*e+d*c^1//3)^-1*(1+x*c^1//3))*log(d+e*x)+-1//2*b*e^-1*log(-1*e*(d*c^1//3+-1*e*(-1)^2//3)^-1*((-1)^2//3+x*c^1//3))*log(d+e*x)+-1//2*b*e^-1*log(e*(-1)^1//3*(d*c^1//3+e*(-1)^1//3)^-1*(1+x*(-1)^2//3*c^1//3))*log(d+e*x))
@test integrate((d+e*x)^-2*(a+b*arctanh(c*x^3)), x) == :(-1*e^-1*(d+e*x)^-1*(a+b*arctanh(c*x^3))+b*c^1//3*(-2*e^3+2*c*d^3)^-1*(e+d*c^1//3)*log(1+x*c^1//3)+b*c^1//3*(2*e^3+2*c*d^3)^-1*(-1*e+d*c^1//3)*log(1+-1*x*c^1//3)+-1*b*3^(1/2)*c^1//3*(2*e^2+2*c^2//3*d^2+2*d*e*c^1//3)^-1*arctan(1//3*3^(1/2)*(1+-2*x*c^1//3))+-1*b*c^1//3*(-4*e^3+4*c*d^3)^-1*(e+d*c^1//3)*log(1+c^2//3*x^2+-1*x*c^1//3)+-1*b*c^1//3*(4*e^3+4*c*d^3)^-1*(-1*e+d*c^1//3)*log(1+x*c^1//3+c^2//3*x^2)+(1/2)*b*c*d^2*e^-1*(-1*e^3+c*d^3)^-1*log(1+c*x^3)+-1*b*3^(1/2)*c^1//3*(2*e^3+2*c*d^3)^-1*(e+d*c^1//3)*arctan(1//3*3^(1/2)*(1+2*x*c^1//3))+-3*b*c*d^2*e^2*(-1*e^6+c^2*d^6)^-1*log(d+e*x)+-1//2*b*c*d^2*e^-1*(e^3+c*d^3)^-1*log(1+-1*c*x^3))
@test integrate(x^3*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b*c^-8*Polylogarithms.polylog(2,1+-2*(1+-1*c*x^(1/2))^-1)+c^-8*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+-1*c*x^(1/2))^-1)+-1*x*c^-6*(a+b*arctanh(c*x^(1/2)))+-1*b^-1*c^-8*(a+b*arctanh(c*x^(1/2)))^2+-11//6*b*c^-7*x^(1/2)+-5//18*b*c^-5*x^3//2+-1//2*c^-4*x^2*(a+b*arctanh(c*x^(1/2)))+-1//3*c^-2*x^3*(a+b*arctanh(c*x^(1/2)))+-1//15*b*c^-3*x^5//2+11//6*b*c^-8*arctanh(c*x^(1/2)))
@test integrate(x^2*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b*c^-6*Polylogarithms.polylog(2,1+-2*(1+-1*c*x^(1/2))^-1)+c^-6*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+-1*c*x^(1/2))^-1)+-1*x*c^-4*(a+b*arctanh(c*x^(1/2)))+-1*b^-1*c^-6*(a+b*arctanh(c*x^(1/2)))^2+-3//2*b*c^-5*x^(1/2)+-1//2*c^-2*x^2*(a+b*arctanh(c*x^(1/2)))+-1//6*b*c^-3*x^3//2+3//2*b*c^-6*arctanh(c*x^(1/2)))
@test integrate(x*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b*c^-4*Polylogarithms.polylog(2,1+-2*(1+-1*c*x^(1/2))^-1)+b*c^-4*arctanh(c*x^(1/2))+c^-4*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+-1*c*x^(1/2))^-1)+-1*b*c^-3*x^(1/2)+-1*x*c^-2*(a+b*arctanh(c*x^(1/2)))+-1*b^-1*c^-4*(a+b*arctanh(c*x^(1/2)))^2)
@test integrate((1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b*c^-2*Polylogarithms.polylog(2,1+-2*(1+-1*c*x^(1/2))^-1)+c^-2*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+-1*c*x^(1/2))^-1)+-1*b^-1*c^-2*(a+b*arctanh(c*x^(1/2)))^2)
@test integrate(x^-1*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b^-1*(a+b*arctanh(c*x^(1/2)))^2+(2a+2*b*arctanh(c*x^(1/2)))*log(2+-2*(1+c*x^(1/2))^-1)+-1*b*Polylogarithms.polylog(2,-1+2*(1+c*x^(1/2))^-1))
@test integrate(x^-2*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(-1*x^-1*(a+b*arctanh(c*x^(1/2)))+b*c^2*arctanh(c*x^(1/2))+b^-1*c^2*(a+b*arctanh(c*x^(1/2)))^2+-1*b*c*x^-1//2+-1*b*c^2*Polylogarithms.polylog(2,-1+2*(1+c*x^(1/2))^-1)+2*c^2*(a+b*arctanh(c*x^(1/2)))*log(2+-2*(1+c*x^(1/2))^-1))
@test integrate(x^-3*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(-1//2*x^-2*(a+b*arctanh(c*x^(1/2)))+b^-1*c^4*(a+b*arctanh(c*x^(1/2)))^2+-1*b*c^4*Polylogarithms.polylog(2,-1+2*(1+c*x^(1/2))^-1)+-1*c^2*x^-1*(a+b*arctanh(c*x^(1/2)))+2*c^4*(a+b*arctanh(c*x^(1/2)))*log(2+-2*(1+c*x^(1/2))^-1)+-3//2*b*c^3*x^-1//2+-1//6*b*c*x^-3//2+3//2*b*c^4*arctanh(c*x^(1/2)))
@test integrate(x^-4*(1+-1*x*c^2)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(-1//3*x^-3*(a+b*arctanh(c*x^(1/2)))+b^-1*c^6*(a+b*arctanh(c*x^(1/2)))^2+-1*b*c^6*Polylogarithms.polylog(2,-1+2*(1+c*x^(1/2))^-1)+-1*c^4*x^-1*(a+b*arctanh(c*x^(1/2)))+2*c^6*(a+b*arctanh(c*x^(1/2)))*log(2+-2*(1+c*x^(1/2))^-1)+-11//6*b*c^5*x^-1//2+-5//18*b*c^3*x^-3//2+-1//2*c^2*x^-2*(a+b*arctanh(c*x^(1/2)))+-1//15*b*c*x^-5//2+11//6*b*c^6*arctanh(c*x^(1/2)))
@test integrate(x^2*(d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :((1/2)*e^-1*x^2*(a+b*arctanh(c*x^(1/2)))+b*d^2*e^-3*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+d^2*e^-3*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+d^2*e^-3*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+(1/2)*b*c^-3*e^-1*x^(1/2)+-1*d*x*e^-2*(a+b*arctanh(c*x^(1/2)))+-2*d^2*e^-3*(a+b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1)+-1//2*b*c^-4*e^-1*arctanh(c*x^(1/2))+-1//2*b*d^2*e^-3*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1//2*b*d^2*e^-3*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+1//6*b*c^-1*e^-1*x^3//2+b*d*c^-2*e^-2*arctanh(c*x^(1/2))+-1*b*d*c^-1*e^-2*x^(1/2))
@test integrate(x*(d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(x*e^-1*(a+b*arctanh(c*x^(1/2)))+b*c^-1*e^-1*x^(1/2)+(1/2)*b*d*e^-2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+(1/2)*b*d*e^-2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+-1*b*d*e^-2*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+-1*b*c^-2*e^-1*arctanh(c*x^(1/2))+-1*d*e^-2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1*d*e^-2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+2*d*e^-2*(a+b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1))
@test integrate((d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(b*e^-1*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+e^-1*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+e^-1*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+-1*e^-1*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1)+-1//2*b*e^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1//2*b*e^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2))))
@test integrate(x^-1*(d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(a*d^-1*log(x)+b*d^-1*Polylogarithms.polylog(2,c*x^(1/2))+d^-1*(2a+2*b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1)+(1/2)*b*d^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+(1/2)*b*d^-1*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+-1*b*d^-1*Polylogarithms.polylog(2,-1*c*x^(1/2))+-1*b*d^-1*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+-1*d^-1*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1*d^-1*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2))))
@test integrate(x^-2*(d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(-1*d^-1*x^-1*(a+b*arctanh(c*x^(1/2)))+b*e*d^-2*Polylogarithms.polylog(2,-1*c*x^(1/2))+b*e*d^-2*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+b*c^2*d^-1*arctanh(c*x^(1/2))+e*d^-2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+e*d^-2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+-1*a*e*d^-2*log(x)+-1*b*c*d^-1*x^-1//2+-1*b*e*d^-2*Polylogarithms.polylog(2,c*x^(1/2))+-2*e*d^-2*(a+b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1)+-1//2*b*e*d^-2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1//2*b*e*d^-2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2))))
@test integrate(x^-3*(d+e*x)^-1*(a+b*arctanh(c*x^(1/2))), x) == :(-1//2*d^-1*x^-2*(a+b*arctanh(c*x^(1/2)))+a*d^-3*e^2*log(x)+b*d^-3*e^2*Polylogarithms.polylog(2,c*x^(1/2))+e*d^-2*x^-1*(a+b*arctanh(c*x^(1/2)))+(1/2)*b*c^4*d^-1*arctanh(c*x^(1/2))+(1/2)*b*d^-3*e^2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+(1/2)*b*d^-3*e^2*Polylogarithms.polylog(2,1+-2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+-1*b*d^-3*e^2*Polylogarithms.polylog(2,-1*c*x^(1/2))+-1*b*d^-3*e^2*Polylogarithms.polylog(2,1+-2*(1+c*x^(1/2))^-1)+-1*d^-3*e^2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+e^(1/2)*x^(1/2)))+-1*d^-3*e^2*(a+b*arctanh(c*x^(1/2)))*log(2*c*(1+c*x^(1/2))^-1*(-1*e^(1/2)+c*(-1d)^(1/2))^-1*((-1d)^(1/2)+-1*e^(1/2)*x^(1/2)))+2*d^-3*e^2*(a+b*arctanh(c*x^(1/2)))*log(2*(1+c*x^(1/2))^-1)+-1//2*b*c^3*d^-1*x^-1//2+-1//6*b*c*d^-1*x^-3//2+b*c*e*d^-2*x^-1//2+-1*b*e*c^2*d^-2*arctanh(c*x^(1/2)))

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.3 (d+e x)^m (a+b arctanh(c x^n))^p.jl


text:
"""
For shifting the modelled trace to be non-negative.
Shift the modelled trace to ensure no negative values.

Arguments:  modelled = the modelled trace.
Returns:    Array{Float64, 1}, zscored and non-negative fluorescence trace
"""

function shiftFluorescence(modelled::Array{Int})
  z_modelled = zscore(modelled)
  adjustment = abs(minimum([0.0, minimum(z_modelled)]))
  return z_modelled + adjustment
end


text:
# Code to asses the effect of moment-compensation on the diffusion signal
# MOST OF THE SYNTAX IN THIS FILE HAS CHANGED, RUN WITH CAUTION
## Read
#Waveform
using Koma, PlotlyJS
b = [0, 50, 200, 400, 600, 790]
E = Array{Any, 2}(undef,6,4)
FID = false
path_file = "/home/ccp/Documents/Koma.jl/"
#Timings
f = readlines(path_file*"qte_timings.txt")
timings = tryparse.(Float64, split(f[2]," "))
δ1, rf180, δ2 = timings[timings .≠ nothing]*1e-3
#Normal DTI
G, n1, n2 = Koma.read_diff_fwf(path_file*"QTE_Waveforms/MX_MC2_b790.txt")
Gmax = 21e-3
g1, g2 = G[1:20:n1,1], G[1:20:n2,4]
seq_dti = Sequence(Gmax*[Grad(g1,δ1) Delay(rf180) -Grad(g2,δ2)])
#ACQ
if FID
    seq_ACQ = Sequence([Grad(0,0)], [RF(0,0)], ADC(1,0))
else
    FOV = 4e-2 #10cm
    Nxy = 20
    ΔT = 20e-6 #100um
    Gmax_acq = 60e-3
    seq_ACQ,_,_,_ = PulseDesigner.EPI_base(FOV, Nxy, ΔT, Gmax_acq)
end
#Phantom
Ns = 10_000; Nsx = floor(Int64, sqrt(Ns) ); Ns = Nsx*Nsx
ADC = 2e-9 #Normal diffusion coefficient
v = 1e-2 # 2mm/s Capilar flow
Δx = 1e-2 #3mm
x = range(-Δx/2,Δx/2,length=Nsx) #spin coordinates
y = range(-Δx/2,Δx/2,length=Nsx)  #spin coordinates
xx, yy = x .+ y'*0, x*0 .+ y' #grid points
T = dur(seq_dti)+dur(seq_ACQ)/2
phantom  = Phantom(x=xx[:], y=yy[:], Dλ1=ADC*ones(Ns), Dλ2=ADC*ones(Ns))
phantom2 = Phantom(x=xx[:], y=yy[:], ux=(x,y,z,t)-> (x .> Δx*.4).*(-v*t) .+ (x .< -Δx*.4).*(v*t), Dλ1=ADC*ones(Ns), Dλ2=ADC*ones(Ns))
#Iter
for idx = [2]
#Seq
G, n1, n2 = Koma.read_diff_fwf(path_file*"QTE_Waveforms/MX_MC2_b$(b[idx+1]).txt")
g1, g2 = G[1:20:n1,1], G[1:20:n2,4]
seq = Sequence(Gmax*[Grad(g1,δ1) Delay(rf180) -Grad(g2,δ2)])
seq += seq_ACQ
#Simulate
f = sqrt(b[idx+1]/790)
simParams = Dict(:Δt => 10e-6, :Nblocks=> 1)
if FID
    recParams = Dict(:skip_rec => true)
else
    recParams = Dict(:epi => true, :recon => :fft, :Nx => Nxy)
end
S_MM = simulate(phantom, seq; simParams) #Eventualmente se agregará el objeto hardware
S_DTI = simulate(phantom, f*seq_dti+seq_ACQ; simParams)
S_MM2 = simulate(phantom2, seq; simParams)
S_DTI2 = simulate(phantom2, f*seq_dti+seq_ACQ; simParams)
S0 = simulate(phantom, 0*seq; simParams)
S02 = simulate(phantom2, 0*seq; simParams)
#Saving results
E[idx+1,1] = S_MM
E[idx+1,2] = S_DTI
E[idx+1,3] = S_MM2
E[idx+1,4] = S_DTI2
E[1,1] = S0
E[1,2] = S0
E[1,3] = S02
E[1,4] = S02
end

## Plot
i = 3
zmax = 15
plot(heatmap(z=abs.(E[i,1]), transpose=false, zmin=0, zmax=zmax), Layout(;title="No mov. MM-DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,2]), transpose=false, zmin=0, zmax=zmax), Layout(;title="No mov. DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,3]), transpose=false, zmin=0, zmax=zmax), Layout(;title="Mov. MM-DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,4]), transpose=false, zmin=0, zmax=zmax), Layout(;title="Mov. DTI b$(b[i])"))

p = Koma.plot_grads_moments(seq_dti+seq_ACQ)


text:
module Deployer

export execute

function execute(args::Dict{Symbol,Any})
  println("Deploying...")
end

end # module


text:
"""
    Maybe{X}

Type shortcut for "`X` or nothing" or "nullable `X`" in javaspeak. Name
got inspired by our functional friends.
"""
const Maybe{X} = Union{Nothing,X}

"""
    VPtr

A convenience wrapper for "any" (C `void`) pointer.
"""
const VPtr = Ptr{Cvoid}


text:
using StatsFuns
using Base.Test
using Compat

# xlogx & xlogy

println("\ttesting xlogx & xlogy ...")

@test xlogx(0) === 0.0
@test_approx_eq xlogx(2) 2.0 * log(2.0)

@test xlogy(0, 1) === 0.0
@test_approx_eq xlogy(2, 3) 2.0 * log(3.0)

# logistic & logit

println("\ttesting logistic & logit ...")

@test_approx_eq logistic(2) 1.0 / (1.0 + exp(-2.0))
@test_approx_eq logit(0.5) 0.0
@test_approx_eq logit(logistic(2)) 2.0

# log1psq

println("\ttesting log1psq ...")

@test_approx_eq log1psq(0.0) 0.0
@test_approx_eq log1psq(1.0) log1p(1.0)
@test_approx_eq log1psq(2.0) log1p(4.0)

# log1pexp, log1mexp, log2mexp & logexpm1

println("\ttesting log1pexp ...")

@test_approx_eq log1pexp(2.0) log(1.0 + exp(2.0))
@test_approx_eq log1pexp(-2.0) log(1.0 + exp(-2.0))
@test_approx_eq log1pexp(10000) 10000.0
@test_approx_eq log1pexp(-10000) 0.0

println("\ttesting log1mexp ...")

@test_approx_eq log1mexp(-1.0) log1p(- exp(-1.0))
@test_approx_eq log1mexp(-10.0) log1p(- exp(-10.0))

println("\ttesting log2mexp ...")

@test_approx_eq log2mexp(0.0) 0.0
@test_approx_eq log2mexp(-1.0) log(2.0 - exp(-1.0))

println("\ttesting logexpm1 ...")

@test_approx_eq logexpm1(2.0) log(exp(2.0) - 1.0)
@test_approx_eq logexpm1(log1pexp(2.0)) 2.0
@test_approx_eq logexpm1(log1pexp(-2.0)) -2.0

# log1pmx

println("\ttesting log1pmx ...")

@test_approx_eq log1pmx(0.0) 0.0
@test_approx_eq log1pmx(1.0) log(2.0) - 1.0
@test_approx_eq log1pmx(2.0) log(3.0) - 2.0

println("\ttesting logmxp1 ...")

@test_approx_eq logmxp1(1.0) 0.0
@test_approx_eq logmxp1(2.0) log(2.0) - 1.0
@test_approx_eq logmxp1(3.0) log(3.0) - 2.0

# logsumexp

println("\ttesting logsumexp ...")

@test_approx_eq logsumexp(2.0, 3.0) log(exp(2.0) + exp(3.0))
@test_approx_eq logsumexp(10002, 10003) 10000 + logsumexp(2.0, 3.0)

@test_approx_eq logsumexp([1.0, 2.0, 3.0]) 3.40760596444438
@test_approx_eq logsumexp([1.0, 2.0, 3.0] .+ 1000.) 1003.40760596444438

# softmax

println("\ttesting softmax ...")

x = [1.0, 2.0, 3.0]
r = @compat exp.(x) ./ sum(exp.(x))
@test_approx_eq softmax([1.0, 2.0, 3.0]) r
softmax!(x)
@test_approx_eq x r


text:
import NonlinearEquations
import PyPlot
import SparseArrays

#this block corresponds to the "Parameter values" table
D = 1e-9#m^2/s
M = 4.40098e-2#kg/mol
q = 1.5491e-10#m/s
V_ϕ = 34.8395e-6#m^3/mol
κ = 1e-13#m^2
μ = 6.7549e-5#kg/m/s
ρ_s = 1072.4#kg/m^3
ρ_w = 1066.5#kg/m^3
θ = 0.0844#unitless
g = 9.8067#m/s
Λ = κ / (μ * θ)
Γ = -g * (M - ρ_w * V_ϕ)

#this block corresponds to the "Initial conditions" table
#nx = 1000
#nz = 100
#nx = 50
#nz = 5
#nx = 100
#nz = 10
nx = 100
nz = 10
mins = [0, 0]
maxs = [1000, 100]
deltax = (maxs[1] - mins[1]) / nx
deltaz = (maxs[2] - mins[2]) / nz
tfinal = 3.1536e10#seconds
c0 = [j == nz && i * deltax >= 100 && i * deltax <= 500 ? 750.6 : 4.4457e-2 for j = 1:nz, i = 1:nx]#this is similar to the initial conditions in the document, but matches the dirichlet boundary condition

function plot(args...; filename=false)
	numplots = div(length(args), 2)
	fig, axs = PyPlot.subplots(1, numplots, figsize=(4 * numplots + 1, 4))
	for i = 1:numplots
		x = args[(i - 1) * 2 + 1]
		title = args[(i - 1) * 2 + 2]
		if startswith(title, "log")
			img = axs[i].imshow(x, origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], 10 * maxs[2]], vmin=-3)
		else
			img = axs[i].imshow(x, origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], 10 * maxs[2]])
		end
		axs[i].set(xlabel="x [meters]", ylabel="z [decimeters]")
		axs[i].set_title(title)
		fig.colorbar(img, ax=axs[i])
	end
	fig.tight_layout()
	display(fig)
	println()
	if  filename != false
		fig.savefig(filename)
	end
	PyPlot.close(fig)
end

#maybe the key is just that it is diffusing in for their example. change the transport boundary conditions so that there’s no CO2 in the reservoir initially, and it only diffuses in
#still, something weird is going on with the pressure -- gravity seems to be absent from the pressure

z(j) = j * deltaz
ij2k(i, j) = j + (i - 1) * nz
@NonlinearEquations.equations function pressure(p, c)
	NonlinearEquations.setnumequations(length(p))
	#do the internal nodes
	for i = 2:nx - 1, j = 1:nz
		k = ij2k(i, j)
		ju = j < nz ? j + 1 : j
		jd = j > 1 ? j - 1 : j
		il = i > 1 ? i - 1 : i
		ir = i < nx ? i + 1 : i
		kl = ij2k(il, j)
		kr = ij2k(ir, j)
		kd = ij2k(i, jd)
		ku = ij2k(i, ju)
		NonlinearEquations.addterm(k, (p[kr] - 2 * p[k] + p[kl]) / deltax ^ 2)
		NonlinearEquations.addterm(k, (p[ku] - 2 * p[k] + p[kd]) / deltaz ^ 2)
		NonlinearEquations.addterm(k, ((ρ_w * g - Γ * c[ku])* (ju - j) - (ρ_w * g - Γ * c[k]) * (j - jd)) * deltaz / deltaz ^ 2)#gravity term
	end
	#do the bottom nodes on the east/west boundary
	k_east = ij2k(nx, nz)
	k_west = ij2k(1, nz)
	#NonlinearEquations.addterm(k_east, p[k_east] - 1.99998975026e7)
	NonlinearEquations.addterm(k_east, p[k_east] - 2e7)
	NonlinearEquations.addterm(k_west, p[k_west] - 2e7)
	#do the east/west boundaries
	for j = 1:nz - 1
		k_east = ij2k(nx, j)
		ku_east = ij2k(nx, j + 1)
		k_west = ij2k(1, j)
		ku_west = ij2k(1, j + 1)
		#NonlinearEquations.addterm(k_east, p[k_east] - p[ku_east] - g * (ρ_w * (1 - c[k_east] * V_ϕ) + M * c[k_east]) * deltaz)
		#NonlinearEquations.addterm(k_west, p[k_west] - p[ku_west] - g * (ρ_w * (1 - c[k_west] * V_ϕ) + M * c[k_west]) * deltaz)
		#NonlinearEquations.addterm(k_west, p[k_west] - g * ρ_w * (nz - 1 + j) * deltaz)
		#NonlinearEquations.addterm(k_east, p[k_east] + g * ρ_w * (nz - 1 + j) * deltaz)
		NonlinearEquations.addterm(k_west, p[k_west] - 2e7 - g * ρ_w * (nz - 1 - j) * deltaz)
		NonlinearEquations.addterm(k_east, p[k_east] - 2e7 - g * ρ_w * (nz - 1 - j) * deltaz)
	end
end
@NonlinearEquations.equations function concentration(p, c)
	NonlinearEquations.setnumequations(length(c))
	cell_volume = deltax * deltaz * 1#assume it is 1m in the y direction
	areax = deltaz * 1
	areaz = deltax * 1
	for i = 1:nx, j = 1:nz
		k = ij2k(i, j)
		if j == nz && i * deltax >= 100 && i * deltax <= 500
			NonlinearEquations.addterm(k, c[k] - 750.6)
		else
			kl = i > 1 ? ij2k(i - 1, j) : k#node to the left
			kr = i < nx ? ij2k(i + 1, j) : k#node to the right
			kd = j > 1 ? ij2k(i, j - 1) : k#node below
			ku = j < nz ? ij2k(i, j + 1) : k#node above
			for ko in [kl, kr]
				k_upstream = p[ko] - p[k] > 0 ? ko : k
				NonlinearEquations.addterm(k, -Λ * (p[ko] - p[k]) / deltax * c[k_upstream] * areax / cell_volume)
			end
			for (ko, gravity_sign) in zip([kd, ku], [1, -1])
				k_viscous_upstream = p[ko] - p[k] > 0 ? ko : k
				NonlinearEquations.addterm(k, -Λ * (p[ko] - p[k]) / deltaz * c[k_viscous_upstream] * areaz / cell_volume)
				if k != ko#don't let gravity move stuff across the boundary
					k_gravity_upstream = gravity_sign == 1 ? k : ko
					NonlinearEquations.addterm(k, gravity_sign * Λ * (g * ρ_w - Γ * c[k_gravity_upstream]) / deltaz * c[k_gravity_upstream] * areaz / cell_volume)
				end
			end
		end
	end
end

#do explicit time stepping
dt = 1e4 / nz^2 * 10^2
c = c0# .+ 1e-3 * rand(size(c0))
p = zeros(size(c0))
t = 0.0
i = 0
frame_counter = 1
#=
	c = zeros(nz, nx)
	A = pressure_p(zeros(nz, nx), c)
	b = -pressure_residuals(zeros(nz, nx), c)
	p = reshape(A \ b, nz, nx)
	@show hash(p)
	plot(log10.(c), "log10(concentration)", p, "pressure")
	error("blah")
	=#
A = pressure_p(zeros(nz, nx), c0)
b = -pressure_residuals(zeros(nz, nx), c0)
p0 = reshape(A \ b, nz, nx)
p = p0
while t < 365 * 10
	global c, t, p, i, frame_counter
	A = pressure_p(zeros(nz, nx), c)
	b = -pressure_residuals(zeros(nz, nx), c)
	lastp = p
	p = reshape(A \ b, nz, nx)
	dc = reshape(-concentration_residuals(p, c), nz, nx)
	lastc = c
	c += dt * dc
	t += dt / 3600 / 24
	if mod(i, 100) == 0
		tstring = round(t; digits=1)
		plot(c, "log10(concentration) at t=$tstring", p, "pressure", p - p0, "p-p0", c - c0, "c-c0", p - lastp, "dp", c - lastc, "dc")
		frame_counter += 1
	end
	i += 1
end
tstring = round(t; digits=1)
plot(log10.(c), "log10(concentration) at t=$tstring days", p, "pressure at t=$tstring days", c - c0, "c-c0")
#plot(log10.(c), "log10(concentration) at t=$tstring days", p, "pressure at t=$tstring days"; filename="figs/frame_$(lpad(frame_counter, 6, '0')).png")
#run(`ffmpeg -y -r 25 -f image2 -s 1920x1080 -i figs/frame_%06d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p movie.mp4`)


text:
"""
    Sim(m::MDP, p::Policy[, initialstate]; kwargs...)
    Sim(m::POMDP, p::Policy[, updater[, initial_belief[, initialstate]]]; kwargs...)

Create a Sim object that contains everything needed to run and record a single simulation, including model, initial conditions, and metadata.

A vector of `Sim` objects can be executed with [`run`](@ref) or [`run_parallel`](@ref).

## Keyword Arguments
- `rng::AbstractRNG=Random.GLOBAL_RNG`
- `max_steps::Int=typemax(Int)`
- `simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps)`
- `metadata::NamedTuple a named tuple (or dictionary) of metadata for the sim that will be recorded, e.g. `(solver_iterations=500,)`.
"""
abstract type Sim end

struct POMDPSim <: Sim
    simulator::Simulator
    pomdp::POMDP
    policy::Policy
    updater::Updater
    initial_belief::Any
    initialstate::Any
    metadata::NamedTuple
end

problem(sim::POMDPSim) = sim.pomdp

struct MDPSim <: Sim
    simulator::Simulator
    mdp::MDP
    policy::Policy
    initialstate::Any
    metadata::NamedTuple
end

problem(sim::MDPSim) = sim.mdp

"""
    Sim(m::POMDP, p::Policy, metadata=(note="a note",))
    Sim(m::POMDP, p::Policy[, updater[, initial_belief[, initialstate]]]; kwargs...)

Create a `Sim` object that represents a POMDP simulation.
"""
function Sim(pomdp::POMDP,
             policy::Policy,
             up=updater(policy),
             initial_belief=initialstate_distribution(pomdp),
             initialstate=nothing;
             rng::AbstractRNG=Random.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata = NamedTuple()
            )

    if initialstate == nothing && statetype(pomdp) != Nothing
        is = rand(rng, initial_belief)
    else
        is = initialstate
    end
    return POMDPSim(simulator, pomdp, policy, up, initial_belief, is, merge(NamedTuple(), metadata))
end

"""
    Sim(m::MDP, p::Policy, metadata=(note="a note",))
    Sim(m::MDP, p::Policy[, initialstate]; kwargs...)

Create a `Sim` object that represents a MDP simulation.
"""
function Sim(mdp::MDP,
             policy::Policy,
             initialstate=nothing;
             rng::AbstractRNG=Random.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata = NamedTuple()
            )

    if initialstate == nothing && statetype(mdp) != Nothing
        is = POMDPs.initialstate(mdp, rng)
    else
        is = initialstate
    end
    return MDPSim(simulator, mdp, policy, is, merge(NamedTuple(), metadata))
end

POMDPs.simulate(s::POMDPSim) = simulate(s.simulator, s.pomdp, s.policy, s.updater, s.initial_belief, s.initialstate)
POMDPs.simulate(s::MDPSim) = simulate(s.simulator, s.mdp, s.policy, s.initialstate)

default_process(s::Sim, r::Real) = (reward=r,)
default_process(s::Sim, hist::SimHistory) = default_process(s, discounted_reward(hist))

run_parallel(queue::AbstractVector; kwargs...) = run_parallel(default_process, queue; kwargs...)

"""
    run_parallel(queue::Vector{Sim})
    run_parallel(f::Function, queue::Vector{Sim})

Run `Sim` objects in `queue` in parallel and return results as a `DataFrame`.

By default, the `DataFrame` will contain the reward for each simulation and the metadata provided to the sim.

# Arguments
- `queue`: List of `Sim` objects to be executed
- `f`: Function to process the results of each simulation
This function should take two arguments, (1) the `Sim` that was executed and (2) the result of the simulation, by default a `SimHistory`. It should return a named tuple that will appear in the dataframe. See Examples below.

## Keyword Arguments
- `show_progress::Bool`: whether or not to show a progress meter
- `progress::ProgressMeter.Progress`: determines how the progress meter is displayed

# Examples

```julia
run_parallel(queue) do sim, hist
    return (n_steps=n_steps(hist), reward=discounted_reward(hist))
end
```
will return a dataframe with with the number of steps and the reward in it.
"""
function run_parallel(process::Function, queue::AbstractVector, pool::AbstractWorkerPool=default_worker_pool();
                      progress=Progress(length(queue), desc="Simulating..."),
                      proc_warn::Bool=true, show_progress::Bool=true)

    if nworkers(pool) == 1 && proc_warn
        @warn("""
             run_parallel(...) was started with only 1 worker in the pool, so simulations will be run in serial.

             To supress this warning, use run_parallel(..., proc_warn=false).

             To use multiple processes, use addprocs() or the -p option (e.g. julia -p 4) and make sure the correct worker pool is assigned to argument `pool` in the call to run_parallel.
             """)
    end

    if progress in (nothing, false)
        progstr = (progress == nothing) ? "nothing" : "false"
        @warn("run_parallel(..., progress=$progstr) is deprecated. Use run_parallel(..., show_progress=false) instead.")
        show_progress = false
    end

    map_function(args...) = (show_progress ?
                             progress_pmap(args..., progress=progress) : pmap(args...))

    frame_lines = map_function(pool, queue) do sim
        result = simulate(sim)
        output = process(sim, result)
        return merge(sim.metadata, output)
    end

    return create_dataframe(frame_lines)
end

Base.run(queue::AbstractVector) = run(default_process, queue)

"""
    run(queue::Vector{Sim})
    run(f::Function, queue::Vector{Sim})

Run the `Sim` objects in `queue` on a single process and return the results as a dataframe.

See `run_parallel` for more information.
"""
function Base.run(process::Function, queue::AbstractVector; show_progress=true)
    lines = []
    if show_progress
        @showprogress for sim in queue
            result = simulate(sim)
            output = process(sim, result)
            line = merge(sim.metadata, output)
            push!(lines, line)
        end
    else
        for sim in queue
            result = simulate(sim)
            output = process(sim, result)
            line = merge(sim.metadata, output)
            push!(lines, line)
        end
    end
    return create_dataframe(lines)
end

function create_dataframe(lines::Vector)
    master = Dict{Symbol, AbstractVector}()
    for line in lines
        push_line!(master, line)
    end
    return DataFrame(master)
end

function push_line!(d::Dict{Symbol, AbstractVector}, line::NamedTuple)
    if isempty(d)
        len = 0
    else
        len = length(first(values(d)))
    end
    for (key, val) in pairs(line)
        T = typeof(val)
        if !haskey(d, key)
            d[key] = Vector{Union{T,Missing}}(missing, len)
        end
        data = d[key]
        if !isa(val,eltype(data))
            d[key] = convert(Array{promote_type(typeof(val), eltype(data)),1}, data)
        end
        push!(d[key], val)
    end
    for da in values(d)
        if length(da) < len + 1
            push!(da, missing)
        end
    end
    return d
end


text:
module ToyStencils

using ArgCheck: @check
using Preferences: @delete_preferences!, @load_preference, @set_preferences!

if isdefined(Base.Experimental, :Tapir)
    using Base.Experimental.Tapir: @sync, @spawn
else
    using Base.Threads: @sync, @spawn
end

include("debug.jl")
include("trapezoid.jl")

end


text:
using Base.Threads: @spawn
using Base.Iterators: drop
using ProgressMeter: @showprogress, Progress, next!
using NamedTupleTools: namedtuple

read_parquet(path, cols::Vector{Symbol}; kwargs...) = read_parquet(path, String.(cols); kwargs...)

read_parquet(path; kwargs...) = read_parquet(path, String[]; kwargs...)

function read_parquet(path, cols::Vector{String}; multithreaded=true, verbose = false)
	"""function for reading parquet"""

    par = ParFile(path)
	nc = ncols(par)

    colnames = [sch.name for sch in  drop(par.schema.schema, 1)]

    close(par)

	if length(cols) == 0
		colnums = collect(1:nc)
	else
		colnums = [findfirst(==(c), colnames) for c in cols]
	end

	results = Vector{Any}(undef, length(colnums))

    filemetadata = metadata(path)

    symbol_col_names = collect(Symbol(col) for col in colnames[colnums])

    p = Progress(length(colnums))
    if multithreaded
		for (i, j) in enumerate(colnums)
            results[i] = @spawn begin
                # next!(p)
                res = read_column(path, filemetadata, j)
                res
            end
        end
        results = fetch.(results)
    else

		for (i, j) in enumerate(colnums)
            results[i] = read_column(path, filemetadata, j)
            next!(p)
		end
	end

    return namedtuple(symbol_col_names, results)
end


text:
module TestFirstOrder2

using ModiaLang
using DifferentialEquations
@usingModiaPlot
using Test

# using RuntimeGeneratedFunctions
# RuntimeGeneratedFunctions.init(@__MODULE__)

inputSignal(t) = sin(t)

FirstOrder1 = Model(
    T = 0.2,
    x = Var(init=0.3),
    equations = :[u = inputSignal(time/u"s"),
                  T * der(x) + x = u,
                  y = 2*x]
)

FirstOrder2 = FirstOrder1 | Map(T = 0.3, x = Var(init=0.6))

firstOrder = @instantiateModel(FirstOrder2, logCode=false)

simulate!(firstOrder, Tsit5(), stopTime = 10, merge = Map(T = 0.4, x = 0.9),
          log=false, logParameters=true, logStates=true,
          requiredFinalStates = [-0.17964872595554535])

# Test get_result(instantiatedModel)
println()
result1 = get_result(firstOrder)
@show(result1[1:10,:])
println()
@show(result1[1:10, ["time", "u", "y"]])

println()
result2 = get_result(firstOrder, onlyStates=true, extraNames=["y"])
@show(result2[1:10,:])

println()
result3 = get_result(firstOrder, extraNames=["y"])
@show(result3[1:10,:])

# Linearize
println("\n... Linearize at stopTime = 0 and 10:")
(A_0 , x_0)  = linearize!(firstOrder, analytic = true)
(A_10, x_10) = linearize!(firstOrder, stopTime=10, analytic = true)
(A_10_numeric, x_10_numeric) = linearize!(firstOrder, stopTime=10, analytic=false)
xNames = get_xNames(firstOrder)
@show xNames
@show A_0 , x_0
@show A_10, x_10
@show A_10_numeric, x_10_numeric
@test isapprox(A_0,[-1/0.4])
@test isapprox(A_0, A_10)

plot(result1, [("u", "x"), "der(x)", "y"])

FirstOrder3 = Model(
    T = 2u"hr",
    x = Var(init=1.0),
    equations = :[u = if after(1.5u"hr"); 1.0 else 0.0 end,
                  T * der(x) + x = u]
)
firstOrder3 = @instantiateModel(FirstOrder3, logCode=false)
simulate!(firstOrder3, Tsit5(), stopTime = 10u"hr")
plot(firstOrder3, [("u", "x"), "der(x)"], figure=2)

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2, false),),
    ((2, true),),
    ((4, false),),
    ((4, true),),
    ]
insts["fast"] = [
    ((50, false),),
    ((50, true),),
    ((1000, false),),
    ((250, true),),
    ]
insts["various"] = [
    ((1000, false),),
    ((500, true),),
    ((4000, false),),
    ((1000, true),),
    ((8000, false),),
    ((2000, true),),
    ((10000, false),),
    ]
return (DiscreteMaxLikelihood, insts)


text:
module Variables
export initialize_decision_variables

using JuMP
function initialize_decision_variables(model_dict)
    """
    Decision variables are undefined at model initialization, the purpose of the
    optimisation is to find the values that the decision variables should take,
    based on a combination of those decision variables being minimized/maximized.

    The decision variables of interest are:

    energy_cap : Capacity (in kW) of a technology at a given location,
        limiting produced energy
    carrier_prod : Produced energy at every timestep (in kWh) from a technology
        at a given location
    carrier_con : Consumed energy at every timestep (in kWh) from a technology
        at a given location
    cost : Cost associated with a technology (sum of investment and operational)
    resource_area : available area (m2) over which to capture an energy resource
    storage_cap : Stored volume capacity (in kWh) of a technology at a given
        location, limiting stored energy.
    storage : Stored volume at every timestep (in kWh) for a technology at a
        given location
    resource_cap : Resource capacity (in kW or kW/m2 is resource_area used) of
        a technology at a given location, limiting available energy to consume
    resource : Resource consumed at every timestep (in kWh) from a technology
        at a given location
    carrier_export : Energy exported out of the system at every timestep (i.e.
        energy is produced but not consumed via another technology's energy_con,
        instead it is sent outside the system)
    cost_om : Operational and maintenance costs associated with a technology at
        a given location, in every timestep
    cost_investment : Fixed investment costs associated with a technology at a
        given location
    purchased : Binary variable stating whether or not a technology has a
        non-zero capacity in the system at a given location
    units : Integer variable stating how many of a technology (with fixed
        capacity) exist in the system at a given location
    operating_units : Integer variable stating how many integer units of a
        technology are committed in a given timestep, at a given location

    """
    backend_model = model_dict["backend_model"]
    sets = model_dict["sets"]
    variable_dict = Dict()
    #
    ## Variables which are always assigned
    #

    # Capacity
    variable_dict["energy_cap"] = @variable(
        backend_model,
        energy_cap[loc_techs=sets["loc_techs"]] >= 0
    );

    # Dispatch
    variable_dict["carrier_prod"] = @variable(
        backend_model,
        carrier_prod[loc_tech_carriers=sets["loc_tech_carriers_prod"],
                     timesteps=sets["timesteps"]] >= 0
    );
    variable_dict["carrier_con"] = @variable(
        backend_model,
        carrier_con[loc_tech_carriers=sets["loc_tech_carriers_con"],
                    timesteps=sets["timesteps"]] <= 0
    );

    # Costs
    variable_dict["cost"] = @variable(
        backend_model,
        cost[costs=sets["costs"], loc_techs=sets["loc_techs_cost"]]
    );

    #
    ## Conditionally assigned variables
    #

    if haskey(sets, "loc_techs_area")
        # Capacity
        variable_dict["resource_area"] = @variable(
            backend_model,
            resource_area[loc_techs=sets["loc_techs_area"]] >= 0
        );
    end

    if haskey(sets, "loc_techs_store")
        # Capacity
        variable_dict["storage_cap"] = @variable(
            backend_model,
            storage_cap[loc_techs=sets["loc_techs_store"]] >= 0
        );
        # Dispatch
        variable_dict["storage"] = @variable(
            backend_model,
            storage[loc_techs=sets["loc_techs_store"],
                    timesteps=sets["timesteps"]] >= 0
        );
    end

    if haskey(sets, "loc_techs_finite_resource_supply_plus")
        # Capacity
        variable_dict["resource_cap"] = @variable(
            backend_model,
            resource_cap[loc_techs=sets["loc_techs_finite_resource_supply_plus"]] >= 0
        );
    end

    if haskey(sets, "loc_techs_finite_resource")
        # Dispatch
        variable_dict["resource_con"] = @variable(
            backend_model,
            resource_con[loc_techs=sets["loc_techs_finite_resource"],
                         timesteps=sets["timesteps"]]
        );
    end

    if haskey(sets, "loc_techs_export")
        # Dispatch
        variable_dict["carrier_export"] = @variable(
            backend_model,
            carrier_export[loc_tech_carriers=sets["loc_tech_carriers_export"],
                           timesteps=sets["timesteps"]] >= 0
        );
    end

    if haskey(sets, "loc_techs_om_cost")
        # Costs
        variable_dict["cost_var"] = @variable(
            backend_model,
            cost_var[costs=sets["costs"], loc_techs=sets["loc_techs_om_cost"],
                     timesteps=sets["timesteps"]]
        );
    end

    if haskey(sets, "loc_techs_investment_cost")
        # Costs
        variable_dict["cost_investment"] = @variable(
            backend_model,
            cost_investment[costs=sets["costs"],
                            loc_techs=sets["loc_techs_investment_cost"]]
        );
    end

    # Binary/Integer variables
    if haskey(sets, "loc_techs_purchase")
        # Capacity
        variable_dict["purchased"] = @variable(
            backend_model,
            purchased[loc_techs=sets["loc_techs_purchase"]],
            Bin
        );
    end

    if haskey(sets, "loc_techs_milp")
        # Capacity
        variable_dict["units"] = @variable(
            backend_model,
            units[loc_techs=sets["loc_techs_milp"]] >= 0,
            Int
        );
        # Dispatch
        variable_dict["operating_units"] = @variable(
            backend_model,
            operating_units[loc_techs=sets["loc_techs_milp"],
                            timesteps=sets["timesteps"]] >= 0,
            Int
        );
    end

    return variable_dict
end

end

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


text:
function sem_summary(sem_fit::SemFit; show_fitmeasures = false, color = :light_cyan, digits = 2)
    print("\n")
    println("Fitted Structural Equation Model")
    print("\n")
    printstyled("--------------------------------- Properties --------------------------------- \n"; color = color)
    print("\n")
    println("Optimization algorithm:      $(optimizer(sem_fit))")
    println("Convergence:                 $(convergence(sem_fit))")
    println("No. iterations/evaluations:  $(n_iterations(sem_fit))")
    print("\n")
    println("Number of parameters:        $(n_par(sem_fit))")
    println("Number of observations:      $(n_obs(sem_fit))")
    print("\n")
    printstyled("----------------------------------- Model ----------------------------------- \n"; color = color)
    print("\n")
    print(sem_fit.model)
    print("\n")
    if show_fitmeasures
        printstyled("--------------------------------- Fitmeasures --------------------------------- \n"; color = color)
        print("\n")
        a = fit_measures(sem_fit)
        goal_length = maximum(length.(string.(keys(a)))) + 4
        for k in keys(a)
            key_length = length(string(k))
            print(k); print(repeat(" ", goal_length - key_length)); print(round(a[k]; digits = 2)); print("\n")
        end
    end
    print("\n")
end

function sem_summary(partable::ParameterTable; color = :light_cyan, secondary_color = :light_yellow, digits = 2, show_variables = true)

    if show_variables
        print("\n")
        printstyled("--------------------------------- Variables --------------------------------- \n"; color = color)
        print("\n")
        printstyled("Latent variables:    "; color = color); for var in partable.variables[:latent_vars] print("$var ") end; print("\n")
        printstyled("Observed variables:  "; color = color); for var in partable.variables[:observed_vars] print("$var ") end; print("\n")
        if haskey(partable.variables, :sorted_vars) && (length(partable.variables[:sorted_vars]) > 0)
            printstyled("Sorted variables:    "; color = color); for var in partable.variables[:sorted_vars] print("$var ") end; print("\n")
        end
    end

    print("\n")
    printstyled("---------------------------- Parameter Estimates ----------------------------- \n"; color = color)
    print("\n")

    columns = keys(partable.columns)

    printstyled("Loadings: \n"; color = color)
    print("\n")

    sorted_columns = [:to, :estimate, :identifier, :value_fixed, :start]
    loading_columns = sort_partially(sorted_columns, columns)

    for var in partable.variables[:latent_vars]
        indicator_indices =
            findall(
                (partable.columns[:from] .== var) .&
                (partable.columns[:parameter_type] .== :→) .&
                (partable.columns[:to] .∈ [partable.variables[:observed_vars]])
        )
        loading_array = reduce(hcat, check_round(partable.columns[c][indicator_indices]; digits = digits) for c in loading_columns)

        printstyled(var; color = secondary_color); print("\n")
        print("\n")
        pretty_table(loading_array; header = loading_columns, tf = PrettyTables.tf_borderless, alignment = :l)
        print("\n")

    end

    printstyled("Directed Effects: \n"; color = color)

    regression_indices =
            findall(
                (partable.columns[:parameter_type] .== :→) .&
                (
                    (
                        (partable.columns[:to] .∈ [partable.variables[:observed_vars]]) .&
                        (partable.columns[:from] .∈ [partable.variables[:observed_vars]])
                    ) .|
                    (
                        (partable.columns[:to] .∈ [partable.variables[:latent_vars]]) .&
                        (partable.columns[:from] .∈ [partable.variables[:observed_vars]])
                    ) .|
                    (
                        (partable.columns[:to] .∈ [partable.variables[:latent_vars]]) .&
                        (partable.columns[:from] .∈ [partable.variables[:latent_vars]])
                    )
                )
            )

    sorted_columns = [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start]
    regression_columns = sort_partially(sorted_columns, columns)

    regression_array = reduce(hcat, check_round(partable.columns[c][regression_indices]; digits = digits) for c in regression_columns)
    regression_columns[2] = Symbol("")

    print("\n")
    pretty_table(regression_array; header = regression_columns, tf = PrettyTables.tf_borderless, alignment = :l)
    print("\n")

    printstyled("Variances: \n"; color = color)

    variance_indices =
            findall(
                (partable.columns[:parameter_type] .== :↔) .&
                (partable.columns[:to] .== partable.columns[:from])
            )

    sorted_columns = [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start]
    variance_columns = sort_partially(sorted_columns, columns)

    variance_array = reduce(hcat, check_round(partable.columns[c][variance_indices]; digits = digits) for c in variance_columns)
    variance_columns[2] = Symbol("")

    print("\n")
    pretty_table(variance_array; header = variance_columns, tf = PrettyTables.tf_borderless, alignment = :l)
    print("\n")

    printstyled("Covariances: \n"; color = color)

    variance_indices =
            findall(
                (partable.columns[:parameter_type] .== :↔) .&
                (partable.columns[:to] .!= partable.columns[:from])
            )

    sorted_columns = [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start]
    variance_columns = sort_partially(sorted_columns, columns)

    variance_array = reduce(hcat, check_round(partable.columns[c][variance_indices]; digits = digits) for c in variance_columns)
    variance_columns[2] = Symbol("")

    print("\n")
    pretty_table(variance_array; header = variance_columns, tf = PrettyTables.tf_borderless, alignment = :l)
    print("\n")

    mean_indices =
        findall(
            (partable.columns[:parameter_type] .== :→) .&
            (partable.columns[:from] .== Symbol("1"))
        )

    if length(mean_indices) > 0

        printstyled("Means: \n"; color = color)

        sorted_columns = [:from, :parameter_type, :to, :estimate, :identifier, :value_fixed, :start]
        variance_columns = sort_partially(sorted_columns, columns)

        variance_array = reduce(hcat, check_round(partable.columns[c][mean_indices]; digits = digits) for c in variance_columns)
        variance_columns[2] = Symbol("")

        print("\n")
        pretty_table(variance_array; header = variance_columns, tf = PrettyTables.tf_borderless, alignment = :l)
        print("\n")
    end

    #printstyled("""No need to copy and paste results, you can use CSV.write(DataFrame(my_partable), "myfile.csv")"""; hidden = true)

end

function sem_summary(partable::EnsembleParameterTable; color = :light_cyan, secondary_color = :light_yellow, digits = 2, show_variables = true)

    if show_variables
        print("\n")
        printstyled("--------------------------------- Variables --------------------------------- \n"; color = color)
        print("\n")
        let partable = partable.tables[[keys(partable.tables)...][1]]
            printstyled("Latent variables:    "; color = color); for var in partable.variables[:latent_vars] print("$var ") end; print("\n")
            printstyled("Observed variables:  "; color = color); for var in partable.variables[:observed_vars] print("$var ") end; print("\n")
            if haskey(partable.variables, :sorted_vars) && (length(partable.variables[:sorted_vars]) > 0)
                printstyled("Sorted variables:    "; color = color); for var in partable.variables[:sorted_vars] print("$var ") end; print("\n")
            end
        end
    end
    print("\n")

    for k in keys(partable.tables)
        print("\n")
        printstyled(rpad(" Group: $k", 78), reverse = true)
        print("\n")
        sem_summary(partable.tables[k]; color = color, secondary_color = secondary_color, digits = digits, show_variables = false)
    end

    # printstyled("""No need to copy and paste results, you can use CSV.write(DataFrame(my_partable), "myfile.csv")"""; hidden = true)

end

function check_round(vec; digits)
    if eltype(vec) <: Number
        return round.(vec; digits = digits)
    end
    return vec
end

function sort_partially(sorted, to_sort)
    out = Vector{eltype(to_sort)}()
    for el in sorted
        if el ∈ to_sort
            push!(out, el)
        end
    end
    remaining = setdiff(to_sort, sorted)
    append!(out, remaining)
    return out
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module QAliases

using ..QT
using ..Empty
using ..SK
using ..EA
using ..PercLinear
using ..PercStep
using ..CommStep

export GraphQ0T, GraphQSKT, GraphQSKNormalT, GraphQEAT, GraphQPercLinearT, GraphQPercStepT, GraphQCommStepT

const GraphQ0T{fourK} = GraphQuant{fourK,GraphEmpty}

"""
    GraphQ0T(N::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph

Shortcut for GraphQuant(N, M, Γ, β, GraphEmpty, N).

See [`GraphQuant`](@ref).

Intended for testing/debugging purposes.
"""
GraphQ0T(Nk::Integer, M::Integer, Γ::Float64, β::Float64) = GraphQuant(Nk, M, Γ, β, GraphEmpty, Nk)

const GraphQSKT{fourK} = GraphQuant{fourK,GraphSK}

"""
    GraphQSKT(N::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph

Shortcut for GraphQuant(N, M, Γ, β, GraphSK, N), slightly optimized.

See [`GraphQuant`](@ref).
"""
GraphQSKT(Nk::Integer, M::Integer, Γ::Float64, β::Float64) = GraphQuant(Nk, M, Γ, β, GraphSK, SK.gen_J(Nk))

const GraphQSKNormalT{fourK} = GraphQuant{fourK,GraphSKNormal}
GraphQSKNormalT(Nk::Integer, M::Integer, Γ::Float64, β::Float64) = GraphQuant(Nk, M, Γ, β, GraphSKNormal, SK.gen_J_gauss(Nk))

const GraphQEAT{fourK,twoD} = GraphQuant{fourK,GraphEANormal{twoD}}

# """
#     GraphQEAT(L::Integer, D::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphQEAT(L::Integer, D::Integer, M::Integer, Γ::Float64, β::Float64)
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    A = EA.gen_EA(L, D)
    N = length(A)
    J = EA.gen_J(Float64, N, A) do
        4*rand() - 2
    end

    GraphQuant(N, M, Γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphQEAT(fname::AbstractString, M::Integer, Γ::Float64, β::Float64)
    L, D, A, J = EA.gen_AJ(fname)
    N = length(A)
    @assert N == L^D
    GraphQuant(N, M, Γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphQEAT(X::GraphEANormal{twoD}, M::Integer, Γ::Float64, β::Float64) where {twoD}
    @assert iseven(twoD)
    D = twoD ÷ 2
    N = X.N
    L = round(Int, N^(1/D))
    @assert L^D == N
    GraphQuant(N, M, Γ, β, GraphEANormal{twoD}, L, X.A, X.J)
end

const GraphQPercLinearT{fourK} = GraphQuant{fourK,GraphPercLinear}

# """
#     GraphQPercLinearT(N::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphQPercLinearT(N::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64)
    ξ, ξv = PercLinear.gen_ξ(N, P)
    GraphQuant(N, M, Γ, β, GraphPercLinear, ξ, ξv)
end

function GraphQPercLinearT(X::GraphPercLinear, M::Integer, Γ::Float64, β::Float64)
    GraphQuant(X.N, M, Γ, β, GraphPercLinear, X.ξ, X.ξv)
end

const GraphQPercStepT{fourK} = GraphQuant{fourK,GraphPercStep}

# """
#     GraphQPercStepT(N::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphQPercStepT(N::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64)
    ξ, ξv = PercStep.gen_ξ(N, P)
    GraphQuant(N, M, Γ, β, GraphPercStep, ξ, ξv)
end

function GraphQPercStepT(X::GraphPercStep, M::Integer, Γ::Float64, β::Float64)
    GraphQuant(X.N, M, Γ, β, GraphPercStep, X.ξ, X.ξv)
end

const GraphQCommStepT{fourK} = GraphQuant{fourK,GraphCommStep}

# """
#     GraphQCommStepT(K1::Integer, K2::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphQCommStepT(K1::Integer, K2::Integer, P::Integer, M::Integer, Γ::Float64, β::Float64)
    ξ, ξv = CommStep.gen_ξ(K1, P)
    N = K1 * K2
    GraphQuant(N, M, Γ, β, GraphCommStep, K2, ξ, ξv)
end

function GraphQCommStepT(X::GraphCommStep, M::Integer, Γ::Float64, β::Float64)
    GraphQuant(X.N, M, Γ, β, GraphCommStep, X.K2, X.ξ, X.ξv)
end

end # module


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import SDPNAL
# With the default tolerance (`tol=1e-6`), the termination status is often
# `NUMERICAL_ERROR` when there are inequality constraints.
# See https://github.com/jump-dev/SDPNAL.jl/pull/1
const OPTIMIZER = SDPNAL.Optimizer(tol=1e-4)
MOI.set(OPTIMIZER, MOI.Silent(), true)
const CONFIG = MOIT.TestConfig(atol=1e-3, rtol=1e-3)

# TODO move to MOI.Test
@testset "Upper bound for $set" for set in [MOI.Nonnegatives(1), MOI.PositiveSemidefiniteConeTriangle(1)]
    x, cx = MOI.add_constrained_variables(OPTIMIZER, set)
    fx = MOI.SingleVariable(x[1])
    c = MOI.add_constraint(OPTIMIZER, fx, MOI.LessThan(1.0))
    MOI.set(OPTIMIZER, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    obj = 1.0fx
    MOI.set(OPTIMIZER, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.optimize!(OPTIMIZER)
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), c) ≈ 1.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), c) ≈ -1.0 atol = CONFIG.atol
end
@testset "Lower bound for $set" for set in [MOI.Nonnegatives(1), MOI.PositiveSemidefiniteConeTriangle(1)]
    x, cx = MOI.add_constrained_variables(OPTIMIZER, set)
    fx = MOI.SingleVariable(x[1])
    c = MOI.add_constraint(OPTIMIZER, fx, MOI.GreaterThan(1.0))
    MOI.set(OPTIMIZER, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = 1.0fx
    MOI.set(OPTIMIZER, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.optimize!(OPTIMIZER)
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), c) ≈ 1.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), c) ≈ 1.0 atol = CONFIG.atol
end
@testset "Off-diagonal bounds for $set" for set in [MOI.Interval(-2.0, 2.0), MOI.Interval(-2.0, 3.0), MOI.Interval(-3.0, 2.0)]
    x, cx = MOI.add_constrained_variables(OPTIMIZER, MOI.PositiveSemidefiniteConeTriangle(2))
    a, b, c = x
    fa = MOI.SingleVariable(a)
    fb = MOI.SingleVariable(b)
    fc = MOI.SingleVariable(c)
    con_a = MOI.add_constraint(OPTIMIZER, fa, set)
    con_b = MOI.add_constraint(OPTIMIZER, fb, set)
    con_c = MOI.add_constraint(OPTIMIZER, fc, set)
    MOI.set(OPTIMIZER, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    obj = 1.0fa + 1.0fc
    MOI.set(OPTIMIZER, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.optimize!(OPTIMIZER)
    @test MOI.get(OPTIMIZER, MOI.ObjectiveValue()) ≈ 2set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), a) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), b) ≈ 0.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), c) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), cx) ≈ zeros(3) atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_a) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_a) ≈ -1.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_b) ≈ 0.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_b) ≈ 0.0 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_c) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_c) ≈ -1.0 atol = CONFIG.atol
end
@testset "Off-diagonal bounds for $set" for set in [MOI.Interval(-2.0, 2.0), MOI.Interval(-2.0, 3.0), MOI.Interval(-3.0, 2.0)]
    x, cx = MOI.add_constrained_variables(OPTIMIZER, MOI.PositiveSemidefiniteConeTriangle(2))
    a, b, c = x
    fa = MOI.SingleVariable(a)
    fb = MOI.SingleVariable(b)
    fc = MOI.SingleVariable(c)
    con_a = MOI.add_constraint(OPTIMIZER, fa, set)
    con_b = MOI.add_constraint(OPTIMIZER, fb, set)
    con_c = MOI.add_constraint(OPTIMIZER, fc, set)
    MOI.set(OPTIMIZER, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    obj = 1.0fb
    MOI.set(OPTIMIZER, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.optimize!(OPTIMIZER)
    @test MOI.get(OPTIMIZER, MOI.ObjectiveValue()) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), a) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), b) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.VariablePrimal(), c) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), cx) ≈ [1/4, -1/4, 1/4] atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_a) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_a) ≈ -1/4 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_b) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_b) ≈ -1/4 atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintPrimal(), con_c) ≈ set.upper atol = CONFIG.atol
    @test MOI.get(OPTIMIZER, MOI.ConstraintDual(), con_c) ≈ -1/4 atol = CONFIG.atol
end

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "SDPNAL"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(OPTIMIZER, false)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, true)
end

# UniversalFallback is needed for starting values, even if they are ignored by SDPNAL
const CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
MOI.empty!(OPTIMIZER)
const CACHED = MOIU.CachingOptimizer(CACHE, OPTIMIZER)
const BRIDGED = MOIB.full_bridge_optimizer(CACHED, Float64)
MOIB.remove_bridge(BRIDGED, MOIB.Constraint.ScalarSlackBridge{Float64})

@testset "Options" begin
    optimizer = SDPNAL.Optimizer(printlevel = 1)
    @test MOI.get(optimizer, MOI.RawParameter("printlevel")) == 1

    param = MOI.RawParameter("bad_option")
    err = MOI.UnsupportedAttribute(param)
    @test_throws err SDPNAL.Optimizer(bad_option = 1)
end

@testset "Unit" begin
    MOIT.unittest(BRIDGED, CONFIG, [
        # Unbounded problem not supported
        "solve_unbounded_model",
        # `TimeLimitSec` not supported.
        "time_limit_sec",
        # `NumberOfThreads` not supported.
        "number_threads",
        # Quadratic functions are not supported
        "solve_qcp_edge_cases", "solve_qp_edge_cases",
        # Integer and ZeroOne sets are not supported
        "solve_integer_edge_cases", "solve_objbound_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end
@testset "Continuous Linear" begin
    # See explanation in `MOI/test/Bridges/lazy_bridge_OPTIMIZER.jl`.
    # This is to avoid `Variable.VectorizeBridge` which does not support
    # `ConstraintSet` modification.
    #MOIB.remove_bridge(BRIDGED, MOIB.Constraint.ScalarSlackBridge{Float64}) # Already removed above
    MOIT.contlineartest(BRIDGED, CONFIG, String[
        # Error using mexMatvec
        # mexMatvec: dimension of 4TH input is incorrect
        # Error in matvecSSN
        # Error in psqmrAB
        # Error in SSNfun
        # Error in sdpnalplus_main_default
        # Error in sdpnalplus
        # >>
        # Intel MKL ERROR: Parameter 6 was incorrect on entry to DGEMV .
        # Error using save
        # Variable 'jx_sdpnalplus_arg_out_1' not found.
        "linear9",
        # Infeasible problems not supported
        "linear8a", "linear8b", "linear8c", "linear12",
        "partial_start"
    ])
end
@testset "Continuous Conic" begin
    MOIT.contconictest(BRIDGED, CONFIG, [
        # Infeasible problems not supported
        "lin3", "norminf2", "normone2", "rotatedsoc2", "soc3",
        # Error using mexMatvec
        # mexMatvec: dimension of 4TH input is incorrect
        # Error in matvecSSN
        # Error in psqmrAB
        # Error in SSNfun
        # Error in sdpnalplus_main_default
        # Error in sdpnalplus
        # >>
        # Intel MKL ERROR: Parameter 6 was incorrect on entry to DGEMV .
        # Error using save
        # Variable 'jx_sdpnalplus_arg_out_1' not found.
        "lin4",
        # >> s2      [iter=1753,(b) reset ADM.maxiter to 200]
        # s2      [iter=1966,(b) reset ADM.maxiter to 300]
        #       [iter=2281,(b) reset ADM.maxiter to 400]
        # s2s2s2s2s2s2s2s2s2s2s2s2s2s2s2s2s2      [iter=2694,(b) reset ADM.maxiter to 500]
        #       [iter=3210,(b) reset ADM.maxiter to 600]
        # s2      [iter=3824,(b) reset ADM.maxiter to 700]
        # s2s2s2s2s2s2s2s2s2      [iter=4537,(b) reset ADM.maxiter to 800]
        # Test Failed at /home/blegat/.julia/packages/MathOptInterface/izyVf/src/Test/contconic.jl:1262
        #   Expression: MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        #    Evaluated: MathOptInterface.NUMERICAL_ERROR == MathOptInterface.OPTIMAL
        "geomean1f", "geomean1v",
        # `ExponentialCone` and `PowerCone` not supported.
        "exp", "dualexp", "pow", "dualpow", "logdet",
        # `RootDetConeSquare` -> `RootDetConeTriangle` bridge missing.
        "rootdets"
    ])
end


text:
#
# Evaluates definite integral
# ∫ erf(b-y)/(sqrt(2)*ℓ), y∈[c,d]
#
# Wolfram code: Integrate[erf(b-y)/(sqrt(2)*ℓ), {y, c, d}]
#

function sub_1(; b = b, c = c, d = d, ℓ = ℓ)

    @assert(c < d)

    sqrt(2/π) * ℓ * ( exp(-(b - c)^2/(2 * ℓ^2)) - exp(-(b - d)^2/(2 * ℓ^2))) +
        (b - c) * erf((b - c)/(sqrt(2) * ℓ)) +
        (d - b) * erf((b - d)/(sqrt(2) * ℓ))

end

function verify_sub_1(; b = b, c = c, d = d, ℓ = ℓ)

    # integrand
    f(y) = erf((b-y)/(sqrt(2)*ℓ))

    # own numerical verification on 1d regular grid
    dx = 1e-6
    X  = c:dx:d
    num = 0.0
    for x in X
        num += f(x)*dx
    end

    # numerical verification via Cuba.jl library
    h = Cuba.vegas((x, out) -> out[1] = (d-c) * f(c + (d-c)*x[1]), 1, 1, maxevals=1_000_000)

    display(h)

    # numerical verification via HCubature.jl library
    hc = hquadrature(f, c, d)

    # return exact and numerical results
    sub_1(; b = b, c = c, d = d, ℓ = ℓ),
    h.integral[1],
    hc[1],
    num

end

#
# Evaluates definite integral
# ∫∫ exp(- (x-y)² / (2ℓ²) ) dx dy, x∈[a,b], y∈[c,d]
#
# Wolfram code: Integrate[E^(-(x - y)^2/(2 (ℓ^2))), {x, a, b}]
#

function integralfunction_1(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    @assert(lowerx < upperx)
    @assert(lowery < uppery)

    # easier names
    a = lowerx
    b = upperx
    c = lowery
    d = uppery

    # Result of ∫ exp(- (x-y)² / (2ℓ²) ) dx  x∈[a,b] is:
    # sqrt(π/2) r (erf((b - y)/(sqrt(2) ℓ)) - erf((a - y)/(sqrt(2) ℓ)))
    # We need to carry out the second integral over y.
    # The result of this second integral over y is:

    sqrt(π/2) * ℓ * (sub_1(b=b, c=c, d=d, ℓ=ℓ) - sub_1(b=a, c=c, d=d, ℓ=ℓ))

end

#
# Verification: numerical evaluation for definite integral
# ∫∫ exp(- (x-y)² / (2ℓ²) ) dx dy, x∈[a, b], y∈[c, d]
#

function verify_1(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    @assert(lowerx < upperx)
    @assert(lowery < uppery)

    # intergrand
    f(x,y) =  exp(-(x-y)^2 / (2*(ℓ^2)))

    # numerical verification on regular grid
    Δx = 1e-3
    X  = lowerx:Δx:upperx
    Y  = lowery:Δx:uppery

    num = 0.0
    for i in 1:length(X)
        for j in 1:length(Y)
            @inbounds num += f(X[i], Y[j])
        end
    end
    num *= Δx^2

    # numerical verification via Cuba.jl libary
    h = divonne((x, out) -> out[1] = (uppery-lowery)*(upperx-lowerx) * f(lowerx + (upperx-lowerx)*x[1], lowery + (uppery-lowery)*x[2]), 2, 1)

    display(h)

    # return exact and numerical results to verify
    integralfunction_1(lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ), num, h.integral[1]

end

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

#
# Evaluates definite integral
# ∫ exp(- (a-y)² / (2ℓ²) ) dy,  y ∈ [c, d]
#
function sub_2a(; a = a, c = c, d = d, ℓ = ℓ)

    # Wolfram code: Integrate[E^(-(a - y)^2/(2 ℓ^2)), {y, c, d}]
    # Solution returned in plain text:
    # integral_c^d e^(-(a - y)^2/(2 ℓ^2)) dy = sqrt(π/2) r (erf((a - c)/(sqrt(2) ℓ)) - erf((a - d)/(sqrt(2) ℓ)))

    sqrt(π/2) * ℓ * (erf((a - c)/(sqrt(2) * ℓ)) - erf((a - d)/(sqrt(2) * ℓ)))

end

#
# Verifies definite integral
# ∫ exp(- (a-y)² / (2ℓ²) ) dy,  y ∈ [c, d]
#
function verify_sub_2a(; a = a, c = c, d = d, ℓ = ℓ)

    f(y) = exp(-(a-y)^2/(2*ℓ^2))

    # own numerical verification on 1d regular grid
    dx = 1e-6
    X  = c:dx:d
    num = 0.0
    for x in X
        num += f(x)*dx
    end

    # numerical verification via Cuba.jl libary
    h = Cuba.vegas((x, out) -> out[1] = (d-c) * f(c + (d-c)*x[1]), 1, 1, maxevals=3_000_000)

    display(h)

    # return exact and numerical results
    sub_2a(; a = a, c = c, d = d, ℓ = ℓ),
    h.integral[1],
    num

end

#
# Evaluates definite integral
# ∫ y erf((a-y)/(sqrt(2)ℓ) dy,  y ∈ [c, d]
#
function sub_2b(; a = a, c = c, d = d, ℓ = ℓ)

    # Wolfram code: Integrate[y Erf[(a - y)/(Sqrt[2] ℓ)], {y, c, d}]
    # Solution returned in plain text:
    # integral_c^d y erf((a - y)/(sqrt(2) ℓ)) dy = 1/2 (a^2 - c^2 + ℓ^2) erf((a - c)/(sqrt(2) ℓ)) - 1/2 (a^2 - d^2 + ℓ^2) erf((a - d)/(sqrt(2) ℓ)) + (ℓ (a + c) e^(-(a - c)^2/(2 ℓ^2)))/sqrt(2 π) - (ℓ (a + d) e^(-(a - d)^2/(2 ℓ^2)))/sqrt(2 π)

    1/2 * (a^2 - c^2 + ℓ^2) * erf((a - c)/(sqrt(2) * ℓ)) -
    1/2 * (a^2 - d^2 + ℓ^2) * erf((a - d)/(sqrt(2) * ℓ)) +
    (ℓ * (a + c) * exp(-(a - c)^2/(2 * ℓ^2)))/sqrt(2*π) - (ℓ*(a + d) * exp(-(a - d)^2/(2*ℓ^2)))/sqrt(2π)

end

#
# Verifies definite integral
# ∫ y erf((a-y)/(sqrt(2)ℓ) dy,  y ∈ [c, d]
#
function verify_sub_2b(; a = a, c = c, d = d, ℓ = ℓ)

    f(y) = y * erf((a-y)/(sqrt(2)*ℓ))

    # own numerical verification on 1d regular grid
    dx = 1e-6
    X  = c:dx:d
    num = 0.0
    for x in X
        num += f(x)*dx
    end

    # numerical verification via Cuba.jl libary
    h = Cuba.vegas((x, out) -> out[1] = (d-c) * f(c + (d-c)*x[1]), 1, 1, maxevals=1_000_000)

    display(h)

    # return exact and numerical results
    sub_2b(; a = a, c = c, d = d, ℓ = ℓ),
    h.integral[1],
    num

end

#
# Evaluates definite integral
# ∫∫ x exp(- (x-y)² / (2ℓ²) ) dx dy, x∈[a, b], y∈[c, d]
#
function integralfunction_2(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    @assert(lowerx < upperx)
    @assert(lowery < uppery)

    # easier names
    a = lowerx
    b = upperx
    c = lowery
    d = uppery

    # Intermediate step:
    # Solution of integral ∫ᵇₐ x exp(-(x-y)^2 / (2ℓ²)) dx is
    # integral_a^b x e^(-(x - y)^2/(2 ℓ^2)) dx = 1/2 ℓ (2 ℓ (e^(-(a - y)^2/(2 ℓ^2)) - e^(-(b - y)^2/(2 ℓ^2))) - sqrt(2 π) y erf((a - y)/(sqrt(2) ℓ)) + sqrt(2 π) y erf((b - y)/(sqrt(2) ℓ)))
    # (Wolfram code: Integrate[x * E^(-((x - y)^2)/(2*ℓ^2)), {x, a, b}]
    #

    #1/2 * ℓ *(2*ℓ * (exp(-(a - y)^2/(2*ℓ^2)) - exp(-(b - y)^2/(2*ℓ^2))) -
    #          sqrt(2π) * y * erf((a - y)/(sqrt(2)*ℓ)) + sqrt(2π) * y * erf((b - y) /(sqrt(2)*ℓ)))

    1/2 * ℓ *(2*ℓ * (sub_2a(a=a, c=c, d=d, ℓ=ℓ) - sub_2a(a=b, c=c, d=d, ℓ=ℓ)) -
              sqrt(2π) * (sub_2b(a = a, c = c, d = d, ℓ = ℓ) -  sub_2b(a = b, c = c, d = d, ℓ = ℓ)))

end

function verify_2(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    f(x,y) =  x * exp(- (x-y)^2 / (2*ℓ^2) )

    # numerical verification on regular grid
    Δx = 1e-3
    X  = lowerx:Δx:upperx
    Y  = lowery:Δx:uppery

    num = 0.0
    for i in 1:length(X)
        for j in 1:length(Y)
            @inbounds num += f(X[i], Y[j])
        end
    end
    num *= Δx^2

    # numerical verification via Cuba.jl libary
    h = divonne((x, out) -> out[1] = (uppery-lowery)*(upperx-lowerx) * f(lowerx + (upperx-lowerx)*x[1], lowery + (uppery-lowery)*x[2]), 2, 1)

    display(h)

    # return exact and numerical results to verify
    integralfunction_2(lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ), num, h.integral[1]

end

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

#
# Evaluate integral:
# ∫ y² erf( (a-y) / (sqrt(2) ℓ) ) dy
#

function sub_3a(; a = a, c = c, d = d, ℓ = ℓ)

    # Wolram code: Integrate[y^2 Erf[(a - y)/(Sqrt[2] ℓ)], {y, c, d}]
    #
    # solution:
    # integral_c^d y y erf((a - y)/(sqrt(2) ℓ)) dy = 1/3 (-(a^3 + 3 a ℓ^2) erf((c - a)/(sqrt(2) ℓ)) + (a^3 + 3 a ℓ^2) erf((d - a)/(sqrt(2) ℓ)) + sqrt(2/π) ℓ e^(-(a - c)^2/(2 ℓ^2)) (a^2 + a c + c^2 + 2 ℓ^2) - sqrt(2/π) ℓ e^(-(a - d)^2/(2 ℓ^2)) (a^2 + a d + d^2 + 2 ℓ^2) + c^3 (-erf((a - c)/(sqrt(2) ℓ))) + d^3 erf((a - d)/(sqrt(2) ℓ)))

    1/3 * (-(a^3 + 3*a*ℓ^2) * (erf((c - a)/(sqrt(2)*ℓ)) - erf((d - a)/(sqrt(2)*ℓ))) +
              sqrt(2/π) * ℓ * (exp(-(a - c)^2/(2 * ℓ^2)) * (a^2 + a * c + c^2 + 2 * ℓ^2) -
                               exp(-(a - d)^2/(2 * ℓ^2)) * (a^2 + a * d + d^2 + 2 * ℓ^2)) +
          c^3 * (-erf((a - c)/(sqrt(2) * ℓ))) + d^3 * erf((a - d)/(sqrt(2) * ℓ)))

end

#
# Verify integral:
# ∫ y² erf( (a-y) / (sqrt(2) ℓ) ) dy
#

function verify_sub_3a(; a = a, c = c, d = d, ℓ = ℓ)

    f(y) = y^2 * erf( (a-y) / (sqrt(2) * ℓ) )

    # own numerical verification on 1d regular grid
    dx = 1e-6
    X  = c:dx:d
    num = 0.0
    for x in X
        num += f(x)*dx
    end

    # numerical verification via Cuba.jl libary
    h = Cuba.vegas((x, out) -> out[1] = (d-c) * f(c + (d-c)*x[1]), 1, 1, maxevals=10_000_000)

    display(h)

    # return exact and numerical results
    sub_3a(; a = a, c = c, d = d, ℓ = ℓ),
    h.integral[1],
    num

end

#
# Evaluate integral:
# ∫ y exp( -(a-y)² / (2 ℓ²) ) dy
#

function sub_3b(; a = a, c = c, d = d, ℓ = ℓ)

    # Wolram code: Integrate[y/E^((a - y)^2/(2 ℓ^2)), {y, c, d}]
    #
    # solution:
    # integral_c^d y e^(-(a - y)^2/(2 ℓ^2)) dy = 1/2 ℓ (2 ℓ (e^(-(a - c)^2/(2 ℓ^2)) - e^(-(a - d)^2/(2 ℓ^2))) + sqrt(2 π) a erf((a - c)/(sqrt(2) ℓ)) - sqrt(2 π) a erf((a - d)/(sqrt(2) ℓ)))

    1/2 * ℓ * (2 * ℓ * (exp(-(a - c)^2/(2*ℓ^2)) - exp(-(a - d)^2/(2*ℓ^2))) +
            sqrt(2*π) * a * erf((a - c)/(sqrt(2)*ℓ)) - sqrt(2*π) * a * erf((a - d)/(sqrt(2)*ℓ)))

end

#
# Verify integral:
# ∫ y exp( -(a-y)² / (2 ℓ²) ) dy
#

function verify_sub_3b(; a = a, c = c, d = d, ℓ = ℓ)

    f(y) = y * exp( -(a-y)^2 / (2 * ℓ^2) )

    # own numerical verification on 1d regular grid
    dx = 1e-6
    X  = c:dx:d
    num = 0.0
    for x in X
        num += f(x)*dx
    end

    # numerical verification via Cuba.jl libary
    h = Cuba.vegas((x, out) -> out[1] = (d-c) * f(c + (d-c)*x[1]), 1, 1, maxevals=10_000_000)

    display(h)

    # return exact and numerical results
    sub_3b(; a = a, c = c, d = d, ℓ = ℓ),
    h.integral[1],
    num

end

#
# Evaluate integral:
# ∫∫ y x exp(- (x-y)² / (2ℓ²) ) dx dy, x∈[a, b], y∈[c, d]
#
function integralfunction_3(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    @assert(lowerx < upperx)
    @assert(lowery < uppery)

    # easier names
    a = lowerx
    b = upperx
    c = lowery
    d = uppery

    # Intermediate step: ∫ y x exp(- (x-y)² / (2ℓ²) ) dx dy, x∈[a, b]
    # Wolfram code: Integrate[(x y)/E^((x - y)^2/(2 ℓ^2)), {x, a, b}]
    # solution:
    # integral_a^b y x e^(-(x - y)^2/(2 ℓ^2)) dx = 1/2 ℓ y (2 ℓ (e^(-(a - y)^2/(2 ℓ^2)) - e^(-(b - y)^2/(2 ℓ^2))) - sqrt(2 π) y erf((a - y)/(sqrt(2) ℓ)) + sqrt(2 π) y erf((b - y)/(sqrt(2) ℓ)))
    #
    # We need to integrate the above solution over y y∈[c, d]

    1/2 * ℓ * (2 * ℓ * (sub_3b(; a = a, c = c, d = d, ℓ = ℓ) - sub_3b(; a = b, c = c, d = d, ℓ = ℓ)) -
                    sqrt(2*π) * sub_3a(; a = a, c = c, d = d, ℓ = ℓ) +
                    sqrt(2*π) * sub_3a(; a = b, c = c, d = d, ℓ = ℓ))

end

function verify_3(; lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ)

    # integrand
    f(x,y) =  y * x * exp(- (x-y)^2 / (2ℓ^2) )

    # numerical verification on regular grid
    Δx = 1e-3
    X  = lowerx:Δx:upperx
    Y  = lowery:Δx:uppery

    num = 0.0
    for i in 1:length(X)
        for j in 1:length(Y)
            @inbounds num += f(X[i], Y[j])
        end
    end
    num *= Δx^2

    # numerical verification via Cuba.jl libary
    h = divonne((x, out) -> out[1] = (uppery-lowery)*(upperx-lowerx) * f(lowerx + (upperx-lowerx)*x[1], lowery + (uppery-lowery)*x[2]), 2, 1)

    display(h)

    # return exact and numerical results to verify
    integralfunction_3(lowerx = lowerx, upperx = upperx, lowery = lowery, uppery = uppery, ℓ = ℓ), num, h.integral[1]

end


text:
println("查看某类型最小值，最大值")

println("Max Int =" , typemax(Int))
println("Min Float32 = ", typemin(Float32))

println("查看某个值二进制")
println("bitstring(1.0_01e12) = ", bitstring(1.0_01e12))

println("机器精度 给出的是1.0和下一个Julia可表示的相邻浮点数之差")
println("eps(Float64) = ", eps(Float64))
println("eps(0.0) = ", eps(0.0))

println("nextfloat(0.0) = ", nextfloat(0.0))
println("prevfloat(0.0) = ", prevfloat(0.0))

println("BigFloat可用setprecision和setrounding函数改变默认精度和舍入模式，通常是全局的")
println("但如果放在do-block中，则可改变特定代码精度和舍入模式")

setrounding(BigFloat, RoundUp) do
    x = BigFloat(1) + parse(BigFloat, "0.1")
println("set RoundUp, x = ", x)
end

setrounding(BigFloat, RoundDown) do
    x = BigFloat(1) + parse(BigFloat, "0.1")
println("set RoundDown, x = ", x)
end

setprecision(40) do
    x = BigFloat(1) + parse(BigFloat, "0.1")
println("set precision 40, x = ", x)
end

println("数值后可以不含空格直接接变量，暗指乘法，例如3x=3 * x")
x = 2
println("x=2, 3x^2 + x = 3 * (2^2) + 2 = ", 3x^2 + x)

println("逻辑右移 -2 >>> 1 = ", -2 >>> 1)
println("算术右移 -2 >> 1 = ", -2 >> 1)

println("dot向量化操作，例如 [1, 2, 3] .* 2 = [1 * 2, 2 * 2, 3 * 2]")
println("[1, 2, 3] .* 2 = ", [1, 2, 3] .* 2)

println("链式比较语法糖 1 < 2 <= 2 < 3 == 3 > 2 >= 1 == 1 < 3 != 5 : ", 1 < 2 <= 2 < 3 == 3 > 2 >= 1 == 1 < 3 != 5)

println("Float64(1//3) =", Float64(1 // 3))

max_stars_repo_name: SenZhangAI/example

max_issues_repo_name: SenZhangAI/example

max_forks_repo_name: SenZhangAI/example


text:
using Pkg
using PackageCompiler
using ICI_simulateur

cd("ICI_simulateur/src")

create_sysimage(
    # Liste des paquets à inclure dans l'image système
    [:ICI_simulateur],

    # Chemin de l'image générée
    #sysimage_path = "ICI_simulateur.so",

    # Script de précompilation
    precompile_execution_file = "run_for_precompile.jl",
    #joinpath(@__DIR__, "test", "runtests.jl"),

    # Ligne à ajouter pour générer une image système "portable" pour
    # toutes les architectures x86-64
    #cpu_target = "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)"
)