text:
relpath_safe(path, startpath) = isempty(startpath) ? path : relpath(path, startpath)

function Base.relpath(filename, pkgdata::PkgData)
    if isabspath(filename)
        # `Base.locate_package`, which is how `pkgdata` gets initialized, might strip pieces of the path.
        # For example, on Travis macOS the paths returned by `abspath`
        # can be preceded by "/private" which is not present in the value returned by `Base.locate_package`.
        idx = findfirst(basedir(pkgdata), filename)
        if idx !== nothing
            idx = first(idx)
            if idx > 1
                filename = filename[idx:end]
            end
            filename = relpath_safe(filename, basedir(pkgdata))
        end
    elseif startswith(filename, "compiler")
        # Core.Compiler's pkgid includes "compiler/" in the path
        filename = relpath(filename, "compiler")
    end
    return filename
end

function iswritable(file::AbstractString)  # note this trashes the Base definition, but we don't need it
    return uperm(stat(file)) & 0x02 != 0x00
end

function unique_dirs(iter)
    udirs = Set{String}()
    for file in iter
        dir, basename = splitdir(file)
        push!(udirs, dir)
    end
    return udirs
end

function file_exists(filename)
    filename = normpath(filename)
    isfile(filename) && return true
    alt = get(cache_file_key, filename, nothing)
    alt === nothing && return false
    return isfile(alt)
end

function use_compiled_modules()
    return Base.JLOptions().use_compiled_modules != 0
end

function firstline(ex::Expr)
    for a in ex.args
        isa(a, LineNumberNode) && return a
        if isa(a, Expr)
            line = firstline(a)
            isa(line, LineNumberNode) && return line
        end
    end
    return nothing
end
firstline(rex::RelocatableExpr) = firstline(rex.ex)

newloc(methloc::LineNumberNode, ln, lno) = fixpath(ln)

location_string(file::AbstractString, line) = abspath(file)*':'*string(line)
location_string(file::Symbol, line) = location_string(string(file), line)

function linediff(la::LineNumberNode, lb::LineNumberNode)
    (isa(la.file, Symbol) && isa(lb.file, Symbol) && (la.file::Symbol === lb.file::Symbol)) || return typemax(Int)
    return abs(la.line - lb.line)
end

# Return the only non-trivial expression in ex, or ex itself
function unwrap(ex::Expr)
    if ex.head === :block || ex.head === :toplevel
        for (i, a) in enumerate(ex.args)
            if isa(a, Expr)
                for j = i+1:length(ex.args)
                    istrivial(ex.args[j]) || return ex
                end
                return unwrap(a)
            elseif !istrivial(a)
                return ex
            end
        end
        return nothing
    end
    return ex
end
unwrap(rex::RelocatableExpr) = RelocatableExpr(unwrap(rex.ex))

istrivial(a) = a === nothing || isa(a, LineNumberNode)

isgoto(stmt) = isa(stmt, Core.GotoNode) | isexpr(stmt, :gotoifnot)

function pushex!(exsigs::ExprsSigs, ex::Expr)
    uex = unwrap(ex)
    if is_doc_expr(uex)
        body = uex.args[4]
        if isa(body, Expr) && body.head !== :call   # don't trigger for docexprs like `"docstr" f(x::Int)`
            exsigs[RelocatableExpr(body)] = nothing
        end
        if length(uex.args) < 5
            push!(uex.args, false)
        else
            uex.args[5] = false
        end
    end
    exsigs[RelocatableExpr(ex)] = nothing
    return exsigs
end

## WatchList utilities
function systime()
    # It's important to use the same clock used by the filesystem
    tv = Libc.TimeVal()
    tv.sec + tv.usec/10^6
end
function updatetime!(wl::WatchList)
    wl.timestamp = systime()
end
Base.push!(wl::WatchList, filenameid::Pair{<:AbstractString,PkgId}) =
    push!(wl.trackedfiles, filenameid)
Base.push!(wl::WatchList, filenameid::Pair{<:AbstractString,PkgFiles}) =
    push!(wl, filenameid.first=>filenameid.second.id)
Base.push!(wl::WatchList, filenameid::Pair{<:AbstractString,PkgData}) =
    push!(wl, filenameid.first=>filenameid.second.info)
WatchList() = WatchList(systime(), Dict{String,PkgId}())
Base.in(file, wl::WatchList) = haskey(wl.trackedfiles, file)

@static if Sys.isapple()
    # HFS+ rounds time to seconds, see #22
    # https://developer.apple.com/library/archive/technotes/tn/tn1150.html#HFSPlusDates
    newer(mtime, timestamp) = ceil(mtime) >= floor(timestamp)
else
    newer(mtime, timestamp) = mtime >= timestamp
end

"""
    success = throwto_repl(e::Exception)

Try throwing `e` from the REPL's backend task. Returns `true` if the necessary conditions
were met and the throw can be expected to succeed. The throw is generated from another
task, so a `yield` will need to occur before it happens.
"""
function throwto_repl(e::Exception)
    if isdefined(Base, :active_repl_backend) &&
            Base.active_repl_backend.backend_task.state === :runnable &&
            isempty(Base.Workqueue) &&
            Base.active_repl_backend.in_eval
        @async Base.throwto(Base.active_repl_backend.backend_task, e)
        return true
    end
    return false
end

function printf_maxsize(f::Function, io::IO, args...; maxchars::Integer=500, maxlines::Integer=20)
    # This is dumb but certain to work
    iotmp = IOBuffer()
    for a in args
        print(iotmp, a)
    end
    print(iotmp, '\n')
    seek(iotmp, 0)
    str = read(iotmp, String)
    if length(str) > maxchars
        str = first(str, (maxchars+1)÷2) * "…" * last(str, maxchars - (maxchars+1)÷2)
    end
    lines = split(str, '\n')
    if length(lines) <= maxlines
        for line in lines
            f(io, line)
        end
        return
    end
    half = (maxlines+1) ÷ 2
    for i = 1:half
        f(io, lines[i])
    end
    maxlines > 1 && f(io, ⋮)
    for i = length(lines) - (maxlines-half) + 1:length(lines)
        f(io, lines[i])
    end
end
println_maxsize(args...; kwargs...) = println_maxsize(stdout, args...; kwargs...)
println_maxsize(io::IO, args...; kwargs...) = printf_maxsize(println, io, args...; kwargs...)

"""
    trim_toplevel!(bt)

Truncate a list of instruction pointers, as obtained from `backtrace()` or `catch_backtrace()`,
at the first "top-level" call (e.g., as executed from the REPL prompt) or the
first entry corresponding to a method in Revise or its dependencies.

This is used to make stacktraces obtained with Revise more similar to those obtained
without Revise, while retaining one entry to reveal Revise's involvement.
"""
function trim_toplevel!(bt)
    # return bt       # uncomment this line if you're debugging Revise itself
    n = itoplevel = length(bt)
    for (i, t) in enumerate(bt)
        sfs = StackTraces.lookup(t)
        for sf in sfs
            if sf.func === Symbol("top-level scope") || (isa(sf.linfo, Core.MethodInstance) && isa(sf.linfo.def, Method) && ((sf.linfo::Core.MethodInstance).def::Method).module ∈ (JuliaInterpreter, LoweredCodeUtils, Revise))
                itoplevel = i
                break
            end
        end
        itoplevel < n && break
    end
    deleteat!(bt, itoplevel+1:length(bt))
    return bt
end


text:
mutable struct GraphicsPipeline
    mVkDevice::vk.VkDevice
    mHandleRef::Ref{vk.VkPipeline}
    mCreateInfo::GraphicsPipelineCreateInfo

    function GraphicsPipeline(createInfo::GraphicsPipelineCreateInfo)
        createInfoRef = handleRef(createInfo)
        @assert createInfoRef[].layout != vk.VK_NULL_HANDLE ("GraphicsPipeline requires a PipelineLayout " *
                                     "to be set in its CreateInfo.")
        this = new()
        this.mHandleRef = Ref{vk.VkPipeline}()
        this.mCreateInfo = createInfo
        this.mVkDevice = getLogicalDevice(createInfo.mLayout)
        if (vk.vkCreateGraphicsPipelines(this.mVkDevice, C_NULL, 1, createInfoRef, C_NULL, this.mHandleRef) != vk.VK_SUCCESS)
            error("Failed to create graphics pipeline!")
        end
        return this
    end
end

@class GraphicsPipeline

function handleRef(this::GraphicsPipeline)::Ref{vk.VkPipeline}
    this.mHandleRef
end

function destroy(this::GraphicsPipeline)
    vk.vkDestroyPipeline(this.mVkDevice, this.mHandleRef[], C_NULL)
end

function getLayout(this::GraphicsPipeline)::PipelineLayout
    return this.mCreateInfo.mLayout
end

ext: jl


text:
module LogisticRegression
    using Pluto

    include("logistic_regression.jl")

    run_notebook_logistic_regression() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/logistic_regression.jl"))
    run_notebook_candy_contain_chocolate() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/candy_contain_chocolate.jl"))
end


text:
@service S3

BUCKET_NAME = "aws-jl-test-issues---" * _now_formatted()

try
    S3.create_bucket(BUCKET_NAME)

    @testset "issue 223" begin
        # https://github.com/JuliaCloud/AWS.jl/issues/223
        body = "Hello World!"
        file_name = "contains spaces"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))
            resp = S3.get_object(BUCKET_NAME, file_name)

            @test String(resp) == body
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 227" begin
        @testset "s3 public bucket" begin
            # https://github.com/JuliaCloud/AWS.jl/issues/227
            config = AWSConfig(; creds=nothing)
            resp = S3.get_object("www.invenia.ca", "index.html"; aws_config=config)

            @test !isempty(resp)
        end

        @testset "s3 private bucket" begin
            bucket_name = "aws-jl-test-issues---" * _now_formatted()
            file_name = "hello_world"

            try
                S3.create_bucket(bucket_name)
                S3.put_object(bucket_name, file_name)

                @test_throws AWSException S3.get_object(
                    bucket_name, file_name; aws_config=AWSConfig(; creds=nothing)
                )
            finally
                S3.delete_object(bucket_name, file_name)
                S3.delete_bucket(bucket_name)
            end
        end

        @testset "lambda" begin
            @service Lambda

            @test_throws NoCredentials Lambda.list_functions(;
                aws_config=AWSConfig(; creds=nothing)
            )
        end
    end

    @testset "issue 324" begin
        body = "Hello World!"
        file_name = "streaming.bin"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))
            resp = S3.get_object(BUCKET_NAME, file_name)
            @test String(resp) == body

            # ERROR: MethodError: no method matching iterate(::Base.BufferStream)
            #   => BUG: header `response_stream` is pushed into the query...
            io = Base.BufferStream()
            S3.get_object(
                BUCKET_NAME,
                file_name,
                Dict("response_stream" => io, "return_stream" => true),
            )
            if bytesavailable(io) > 0
                @test String(readavailable(io)) == body
            else
                @test "no body data was available" == body
            end

        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 466" begin
        file_name = "hang.txt"

        try
            S3.put_object(BUCKET_NAME, file_name)

            # The tests below validate the current behavior of how streams are handled.
            # Note: Avoid using `eof` for these tests can hang when using an unclosed `Base.BufferStream`

            stream = S3.get_object(BUCKET_NAME, file_name, Dict("return_stream" => true))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                @test isopen(stream)
            end

            stream = Base.BufferStream()
            S3.get_object(BUCKET_NAME, file_name, Dict("response_stream" => stream))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                # See: https://github.com/JuliaCloud/AWS.jl/issues/471
                @test_broken isopen(stream)
            end

            stream = Base.BufferStream()
            S3.get_object(
                BUCKET_NAME,
                file_name,
                Dict("response_stream" => stream, "return_stream" => true),
            )
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test !isopen(stream)
            else
                @test isopen(stream)
            end
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

    @testset "issue 474" begin
        body = "foo\0bar"
        expected = Vector{UInt8}(body)
        file_name = "null.txt"

        try
            S3.put_object(BUCKET_NAME, file_name, Dict("body" => body))

            raw = S3.get_object(BUCKET_NAME, file_name, Dict("return_raw" => true))
            @test raw == expected

            stream = S3.get_object(BUCKET_NAME, file_name, Dict("return_stream" => true))
            if AWS.DEFAULT_BACKEND[] isa AWS.HTTPBackend
                @test stream isa Base.BufferStream
                @test !isopen(stream)

                if !isopen(stream)
                    @test read(stream) == expected
                end
            else
                @test stream isa IOBuffer
                @test isopen(stream)

                seekstart(stream)
                @test read(stream) == expected
            end
        finally
            S3.delete_object(BUCKET_NAME, file_name)
        end
    end

finally
    S3.delete_bucket(BUCKET_NAME)
end


text:
class DevonshireCream {
    serveOn() {
        return "Scones";
    }
}

print DevonshireCream; // Expect: "DevonshireCream"

class Bagel {}
var bagel = Bagel();
print bagel; // Expect: "Bagel instance"

class C {}
var c = C();
c.field = "Hello, World!";
print c.field;

class Bacon { eat() { print "Crunch crunch crunch!"; } }
Bacon().eat(); // Prints "Crunch crunch crunch!".

class Cake {
    taste() {
        var adj = "declicous";
        print "The " + this.flavor + " cake is " + adj + "!";
    }
}

var cake = Cake();
cake.flavor = "German chocolate";
cake.taste(); // Prints "The German chocolate cake is declicious!".

class Thing {
    getCallback() {
        fun localFunction() {
            print this;
        }
        return localFunction;
    }
}

var callback = Thing().getCallback();
callback();

class Foo {
    init() {
        print this;
    }
}

var foo = Foo();

class Base {}
class Derived < Base {}
var derived = Derived();

class Doughnut {
    cook() {
        print "Fry until golden brown.";
    }
}
class BostonCream < Doughnut {}
BostonCream().cook();

class A {
    method() {
        print "A method";
    }
}

class B < A {
    method() {
        print "B method";
    }

    test() {
        super.method();
    }
}

class C < B {}

C().test();


text:
function test_integral_form()

    using EffectiveWaves

    # physical parameters
    θin = 0.0
    k=1.;
    ho = 2;

    medium = Medium(1.0,1.0+0.0im)
    ω = real(k*medium.c)
    # specie = Specie(ρ=0.1,r=1.0, c=0.5, volfrac=0.2)
    specie = Specie(ρ=0.1,r=1.0, c=0.5, volfrac=0.1)
    # specie = Specie(ρ=0.1,r=0.1, c=0.5, volfrac=0.1)
    specie2 = Specie(ρ=2.0, r=1.0, c=0.1, volfrac=0.15)

    # From effective wave theory
    k_eff0 = wavenumber_low_volfrac(ω, medium, [specie]; hankel_order = ho)
    k_effs = wavenumbers(ω, medium, [specie]; mesh_points = 10, tol = 1e-8, hankel_order = ho)
    k_effs = sort(k_effs, by=imag)

    eff_medium = effective_medium(medium, [specie])
    ω/eff_medium.c
    k_effs2 = wavenumbers(ω, medium, [specie2]; tol = 1e-8, hankel_order = ho)
    k_effs2 = sort(k_effs2, by=imag)

    using OffsetArrays, ApproxFun, IterTools
    include("src/integral_form/integral_form.jl")

    # (x, (MM_quad,b_mat)) = average_wave_system(ω, medium, specie; θin = θin, mesh_points = 501, hankel_order=ho);
    X = 0.0:0.005:30.0
    (MM_quad,b_mat) = average_wave_system(ω, X, medium, specie;θin = θin, hankel_order=ho);

    # discretization parameters
    J = length(collect(x)) - 1

    len = (J + 1) * (2ho + 1)
    MM_mat = reshape(MM_quad, (len, len));
    b = reshape(b_mat, (len));

    As = MM_mat\b;
    As_mat = reshape(As, (J+1, 2ho+1));

    amps_eff1 = scattering_amplitudes_average(ω, x, medium, [specie];
            k_eff = k_effs[1], hankel_order = ho, θin=θin, tol=1e-8)
    amps_eff2 = scattering_amplitudes_average(ω, x, medium, [specie];
            k_eff = k_effs[2], hankel_order = ho, θin=θin, tol=1e-8)
    amps_eff3 = scattering_amplitudes_average(ω, x, medium, [specie];
            k_eff = k_effs[end], hankel_order = ho, θin=θin, tol=1e-8)

    amps0_eff = scattering_amplitudes_average(ω, x, medium, [specie];
            k_eff = k_eff0, hankel_order = ho, θin=θin)
    amps2_eff = scattering_amplitudes_average(ω, x, medium, [specie2];
            k_eff = k_effs2[2], hankel_order = ho, θin=θin, tol=1e-8)
    # sanity check: the abs of reflection coefficients should always be smaller than one.
    reflection_coefficient(ω, medium, specie2; amps = amps2_eff, θin = θin)
    reflection_coefficient(ω, medium, specie; amps = amps_eff2, θin = θin)
    reflection_coefficient(ω, medium, specie; amps = amps_eff3, θin = θin)
    reflection_coefficient(ω, medium, specie; amps = amps_eff1, θin = θin)

    error0_eff = reshape( abs.(MM_mat*amps0_eff.amplitudes[:] .- b), (J+1, 2ho+1))
    error_eff1 = reshape( abs.(MM_mat*amps_eff1.amplitudes[:] .- b), (J+1, 2ho+1))
    error_eff2 = reshape( abs.(MM_mat*amps_eff2.amplitudes[:] .- b), (J+1, 2ho+1))
    error_eff3 = reshape( abs.(MM_mat*amps_eff3.amplitudes[:] .- b), (J+1, 2ho+1))

    error2_eff = reshape( abs.(MM_mat*amps2_eff.amplitudes[:] .- b), (J+1, 2ho+1))
    # With a completely wrong wavenumber I can still get a small error by just scalling the amplitudes
    amps2_eff = scattering_amplitudes_average(ω, x, medium, [specie];
            k_eff = k_effs2[2], hankel_order = ho, θin=θin)
    error2_eff = reshape( abs.(MM_mat*amps2_eff.amplitudes[:] .- b), (J+1, 2ho+1))
    # scale_amplitudes_effective(ω, k_effs[1], amps2_eff.amplitudes, medium, [specie]; θin = θin)

    error = reshape( abs.((MM_mat*As)./b .- 1.0+0.0im), (J+1, 2ho+1))

    using Plots; pyplot(linewidth=2) # unicodeplots()
    plot(xlabel = "depth (1 wavelength = 2π )", ylabel = "error %", ylims=(-0.1,0.5), title="Transmitted wave errors")
    plot!(x,error_eff1[:,ho+1], label = "Eff. error")
    plot!(x,error0_eff[:,ho+1], linestyle=:dash, label = "Eff. low φ error")
    plot!(x,error2_eff[:,ho+1], linestyle=:dot, label = "Eff. wrong k_eff error")
    plot!(x,error[:,ho+1], linestyle=:dashdot, label = "Integral method error")

    plot(x, [real.(As_mat[:,ho+1]),imag.(As_mat[:,ho+1])], labels = ["real sol." "imag sol."])
    plot!(x, [real.(amps_eff1.amplitudes[:,ho+1]),imag.(amps_eff1.amplitudes[:,ho+1])],
        labels = ["real eff. 1" "imag eff. 1"], linestyle=:dash)
    plot!(x, [real.(amps_eff2.amplitudes[:,ho+1]),imag.(amps_eff2.amplitudes[:,ho+1])],
        labels = ["real eff. 2" "imag eff. 2"], linestyle=:dash)

    plot(x, [abs.(As_mat[:,ho+1] .- amps_eff1.amplitudes[:,ho+1])],
            labels = ["diff abs eff. " "diff imag eff."])

    plot(x, [real.(As_mat[:,ho+2]),imag.(As_mat[:,ho+2])], labels = ["real sol." "imag sol."])
    plot!(x, [real.(amps_eff1.amplitudes[:,ho+2]),imag.(amps_eff1.amplitudes[:,ho+2])],
        labels = ["real eff. 1" "imag eff. 1"], linestyle=:dash)

    # plot!(x, [real.(amps_eff3.amplitudes[:,ho+1]),imag.(amps_eff3.amplitudes[:,ho+1])],
    #     labels = ["real eff. 3" "imag eff. 3"], linestyle=:dash)
    # plot!(x, [real.(amps0_eff.amplitudes[:,ho+1]),imag.(amps0_eff.amplitudes[:,ho+1])],
    #     labels = ["real φ eff." "imag φ eff."], linestyle=:dot)

    is = 250:(length(collect(x))-1)
    plot(x[is], log.(abs.(As_mat[is,ho+1])), labels = "abs sol.")
    plot!(x[is], log.(abs.(As_eff_mat[is,ho+1])), labels = "abs eff.")
    # b_mat ≈ [ sum(MM_quad[l,m,j,n]*A_mat[j,n] for j=1:(J+1), n=1:(2M+1)) for l=1:(J+1), m=1:(2M+1)]
end

# ints = [ check_integration(1.0+1.0im; h = 1./n) for n=10:30:310]

function test_check_integration()
#from Mathematics, had several errors when running
    math = [
            [65.0544 - 146.891im, 24.1498 - 8.37225im, -5.10012 + 3.74221im
                , 0.940428 - 3.4201im, 7.25312 - 0.0534466im, 22.9058 + 38.434im, -49.5574 + 92.4377im]
            ,[-32.7224 + 1.52201im, -12.033 - 16.2499im, 5.93884 - 10.8663im
                , 8.51838 + 3.68509im, -4.94285 + 9.32462im, -12.0994 - 5.48456im, 3.6755 - 13.9906im]
            ,[-9.82076 - 18.6418im, 10.2188 - 12.2782im, 11.8699 + 5.34705im
                , -2.698 + 12.0431im, -13.1158 - 0.218682im, -2.92199 - 13.8296im, 12.9993 - 6.4169im]
         ];
    # result of ints = check_integration(h = 1/400., max_x = 8.0)
    julia = [
            [64.9882-146.903im, 24.1561-8.37069im, -5.10105+3.74215im
                , 0.940588-3.4198im, 7.25286-0.0526824im, 22.9026+38.4402im, -49.5882+92.5011im]
            ,[-32.8359+1.66127im, -12.0424-16.2309im, 5.93782-10.8642im
                , 8.51816+3.68541im, -4.94155+9.32627im, -12.1151-5.49728im, 3.77515-13.8555im]
            ,[-9.88661-18.6492im, 10.2111-12.2772im, 11.8693+5.34716im
                , -2.69806+12.0429im, -13.1162-0.21808im, -2.92135-13.8369im, 12.9773-6.359im]
        ];

    using Plots; pyplot()

    scatter(-3:3, abs.(julia[1]./math[1] .- 1), xlab = "hankel_order", ylab = "rel. error", lab ="", title="Mathematica vs Julia for x = 0.0")
    scatter(-3:3, abs.(julia[2]./math[2] .- 1), xlab = "hankel_order", ylab = "rel. error", lab ="", title="Mathematica vs Julia for x = 1.0")
    scatter(-3:3, abs.(julia[3]./math[3] .- 1), xlab = "hankel_order", ylab = "rel. error", lab ="", title="Mathematica vs Julia for x = 2.0")

#from julia
# using JLD
# ints =  first(values(load("integrated_As.jld")))
# data = [ [ints[i][j][k] for i in eachindex(ints)] for j=1:3, k=1:7];
# # i=2;j=3;k=5; ints[i][j][k] == data[j,k][i]
# # x = 0. => j = 1
# M = 3;
#
# hs =  [ 1./n for n=10:30:310]
# j = 3
#  plot()
#  for k=1:7
#      plot!(hs[4:end], abs.(data[j,k][4:end]), label = "m = $(k - M - 1)", xlab = "mesh element size")
#      scatter!([hs[end]], [abs(math[j][k])])
#  end
#   gui()
#
#  plot()
#  for k=1:7
#      plot!(hs[4:end], 1 - abs.(data[j,k][4:end])./abs(data[j,k][end]), label = "m = $(k - M - 1)")
#  end
#   gui()
end

# tests the integration scheme
function check_integration(k_eff::Complex{Float64} = 1.0+1.0im; k=1.,a=1., h = a*k/55., max_x = 6.0)

    A(n,x) = im^Float64(n)*exp(-im*n*θin)*exp(im*x*k_eff)
    k=1.; a=1.;
    # physical parameters
    θin = 0.3

    # discretization parameters
    M = 3;
    J = Int(round(max_x/h)) # choose an even number for Integration schemes
    x = OffsetArray{Float64}(undef, 0:J)
    x[0:J] = (0:J)*h

    # trapezoidal
    σ =  OffsetArray(trap_scheme(collect(x); xn=max_x), 0:J)

    PQ_quad = intergrand_kernel(x, a*k; θin = θin, M = M)

    # Apply integration scheme
    for j=0:J
        PQ_quad[:,:,j+1,:] *= σ[j]
    end
    # PQ_quad = [σ[j]*PQ_quad[l+1,m+M+1,j+1,n+M+1] for l=0:J, m=-M:M, j=0:J, n=-M:M]

    len = (J + 1) * (2M + 1)
    PQ_mat = reshape(PQ_quad, (len, len))

    A_mat = [ A(m,x[l]) for l = 0:J, m = -M:M]
    As = reshape(A_mat, (len))
    As_integrated = reshape(PQ_mat*As, (J+1,2M+1))

    i1 = Int(round(1/h)) # for x1 = 1
    i2 = Int(round(2/h)) # for x1 = 2
    [As_integrated[1,:], As_integrated[1+i1,:], As_integrated[1+i2,:]]
end


text:
def foo# 3

deftype Foo (A -> B)

;; here is a comment
; here is another comment

deftype (Foo a b)
  (a -> b)

(qat :: A -> b)
def tdeftype (lambda a b c -> (a + b))

def (add a b)
  (a + b)

def wat [[1,2,3], [4,5,6]]

def (catstrs strs)
  (foldr f
         (head strs)
         (tail strs))

def strs ["aa", "bb"]

def (mymap f xs)
  if ((length xs) == 0)
      then
xs
  else
    ((f (head xs))
     : (mymap f (tail xs)))

def empty []

def getFile
  (readFile "./parse.js")

;;def fileLines
;;  (getFile >>=
;;   ((mapM_ putStrLn) . lines))

def (testUnary n)
  ((-n) + n)

def (foo bar)
  let {
    lol =  [1,
   (lambda qwerty blah ->
      [qerty, blah,
       (lambda whatever -> whatever)])]
  }
  if bar
    then [lol,(- 1.2),"lulz",lol]
  else if something
    then [,]
  else
    somethingelse

def (splitHelp acc xs ys)
  if (null xs)
    then ((reverse acc), ys)
  else if (null (tail xs))
    then ((reverse acc), ys)
  else
    (splitHelp ((head ys) : acc)
               (tail (tail xs))
               (tail ys))

def (splitxs xs)
  (splitHelp [] xs xs)

def r def
  {
    a = 4
  }
    a

def main
  let {
    (f a) -> a
    unary = (print (testUnary 6))
    splitted = def {
                xs = (fst (splitxs [12,3,4,56]))
                } (xs ++ [0,9])
      }
      if False
        then superduper
      else
        (unary +
         fileLines +
         (print splitted))

def blah (3 / 4)


text:
function subplotlayout(sz::@compat(Tuple{Int,Int}))
  GridLayout(sz...)
end

function subplotlayout(rowcounts::AVec{Int})
  FlexLayout(sum(rowcounts), rowcounts)
end

function subplotlayout(numplts::Int, nr::Int, nc::Int)

  # figure out how many rows/columns we need
  if nr == -1
    if nc == -1
      nr = round(Int, sqrt(numplts))
      nc = ceil(Int, numplts / nr)
    else
      nr = ceil(Int, numplts / nc)
    end
  else
    nc = ceil(Int, numplts / nr)
  end

  # if it's a perfect rectangle, just create a grid
  if numplts == nr * nc
    return GridLayout(nr, nc)
  end

  # create the rowcounts vector
  i = 0
  rowcounts = Int[]
  for r in 1:nr
    cnt = min(nc, numplts - i)
    push!(rowcounts, cnt)
    i += cnt
  end

  FlexLayout(numplts, rowcounts)
end

Base.length(layout::FlexLayout) = layout.numplts
Base.start(layout::FlexLayout) = 1
Base.done(layout::FlexLayout, state) = state > length(layout)
function Base.next(layout::FlexLayout, state)
  r = 1
  c = 0
  for i = 1:state
    c += 1
    if c > layout.rowcounts[r]
      r += 1
      c = 1
    end
  end
  (r,c), state + 1
end

nrows(layout::FlexLayout) = length(layout.rowcounts)
ncols(layout::FlexLayout, row::Int) = row < 1 ? 0 : (row > nrows(layout) ? 0 : layout.rowcounts[row])

# get the plot index given row and column
Base.getindex(layout::FlexLayout, r::Int, c::Int) = sum(layout.rowcounts[1:r-1]) + c

Base.length(layout::GridLayout) = layout.nr * layout.nc
Base.start(layout::GridLayout) = 1
Base.done(layout::GridLayout, state) = state > length(layout)
function Base.next(layout::GridLayout, state)
  r = div(state-1, layout.nc) + 1
  c = mod1(state, layout.nc)
  (r,c), state + 1
end

nrows(layout::GridLayout) = layout.nr
ncols(layout::GridLayout) = layout.nc
ncols(layout::GridLayout, row::Int) = layout.nc

# get the plot index given row and column
Base.getindex(layout::GridLayout, r::Int, c::Int) = (r-1) * layout.nc + c

Base.getindex(subplt::Subplot, args...) = subplt.plts[subplt.layout[args...]]

# handle "linking" the subplot axes together
# each backend should implement the _remove_axis and _expand_limits methods
function link_axis(subplt::Subplot, isx::Bool)

  # collect the list of plots and the expanded limits for those plots that should be linked on this axis
  includedPlots = Any[]
  # lims = [Inf, -Inf]
  lims = Dict{Int,Any}()  # maps column to xlim
  for (i,(r,c)) in enumerate(subplt.layout)

    # shouldlink will be a bool or nothing.  if nothing, then use linkx/y (which is true if we get to this code)
    shouldlink = subplt.linkfunc(r,c)[isx ? 1 : 2]
    if shouldlink == nothing || shouldlink
      plt = subplt.plts[i]

      # if we don't have this
      k = isx ? c : r
      if (firstone = !haskey(lims, k))
        lims[k] = [Inf, -Inf]
      end

      isinner = (isx && r < nrows(subplt.layout)) || (!isx && !firstone)
      push!(includedPlots, (plt, isinner, k))

      _expand_limits(lims[k], plt, isx)
    end

  end

  # do the axis adjustments
  for (plt, isinner, k) in includedPlots
    if isinner
      _remove_axis(plt, isx)
    end
    (isx ? xlims! : ylims!)(plt, lims[k]...)
  end
end

# ------------------------------------------------------------

Base.string(subplt::Subplot) = "Subplot{$(subplt.backend) p=$(subplt.p) n=$(subplt.n)}"
Base.print(io::IO, subplt::Subplot) = print(io, string(subplt))
Base.show(io::IO, subplt::Subplot) = print(io, string(subplt))

getplot(subplt::Subplot, idx::Int = subplt.n) = subplt.plts[mod1(idx, subplt.p)]
getplotargs(subplt::Subplot, idx::Int) = getplot(subplt, idx).plotargs
convertSeriesIndex(subplt::Subplot, n::Int) = ceil(Int, n / subplt.p)

# ------------------------------------------------------------

function validateSubplotSupported()
  if !subplotSupported()
    error(CURRENT_BACKEND.sym, " does not support the subplot/subplot! commands at this time.  Try one of: ", join(filter(pkg->subplotSupported(backendInstance(pkg)), backends()),", "))
  end
end

"""
Create a series of plots:
```
  y = rand(100,3)
  subplot(y; n = 3)             # create an automatic grid, and let it figure out the nr/nc... will put plots 1 and 2 on the first row, and plot 3 by itself on the 2nd row
  subplot(y; n = 3, nr = 1)     # create an automatic grid, but fix the number of rows to 1 (so there are n columns)
  subplot(y; n = 3, nc = 1)     # create an automatic grid, but fix the number of columns to 1 (so there are n rows)
  subplot(y; layout = [1, 2])   # explicit layout by row... plot #1 goes by itself in the first row, plots 2 and 3 split the 2nd row (note the n kw is unnecessary)
  subplot(plts, n; nr = -1, nc = -1)  # build a layout from existing plots
  subplot(plts, layout)               # build a layout from existing plots
```
"""
function subplot(args...; kw...)
  validateSubplotSupported()
  d = Dict(kw)
  preprocessArgs!(d)

  # for plotting recipes, swap out the args and update the parameter dictionary
  args = _apply_recipe(d, args...; kw..., issubplot=true)

  # figure out the layout
  layoutarg = get(d, :layout, nothing)
  if layoutarg != nothing
    layout = subplotlayout(layoutarg)
  else
    n = get(d, :n, -1)
    if n < 0
      error("You must specify either layout or n when creating a subplot: ", d)
    end
    layout = subplotlayout(n, get(d, :nr, -1), get(d, :nc, -1))
  end

  # initialize the individual plots
  pkg = backend()
  plts = Plot{typeof(pkg)}[]
  for i in 1:length(layout)
    di = getPlotArgs(pkg, d, i)
    di[:subplot] = true
    dumpdict(di, "Plot args (subplot $i)")
    push!(plts, _create_plot(pkg; di...))
  end

  # create the object and do the plotting
  subplt = Subplot(nothing, plts, pkg, length(layout), 0, layout, d, false, false, false, (r,c) -> (nothing,nothing))
  subplot!(subplt, args...; d...)

  subplt
end

# ------------------------------------------------------------------------------------------------

# NOTE: for the subplot calls building from existing plots, we need the first plot to be separate to ensure dispatch calls this instead of the more general subplot(args...; kw...)

# grid layout
function subplot{P}(plt1::Plot{P}, plts::Plot{P}...; kw...)
  d = Dict(kw)
  layout = subplotlayout(length(plts)+1, get(d, :nr, -1), get(d, :nc, -1))
  subplot(vcat(plt1, plts...), layout, d)
end

# explicit layout
function subplot{P,I<:Integer}(pltsPerRow::AVec{I}, plt1::Plot{P}, plts::Plot{P}...; kw...)
  layout = subplotlayout(pltsPerRow)
  subplot(vcat(plt1, plts...), layout, Dict(kw))
end

# this will be called internally
function subplot{P<:PlottingPackage}(plts::AVec{Plot{P}}, layout::SubplotLayout, d::Dict)
  validateSubplotSupported()
  p = length(layout)
  n = sum([plt.n for plt in plts])
  subplt = Subplot(nothing, collect(plts), P(), p, n, layout, Dict(), false, false, false, (r,c) -> (nothing,nothing))

  _preprocess_subplot(subplt, d)
  _postprocess_subplot(subplt, d)

  subplt
end

# TODO: hcat/vcat subplots and plots together arbitrarily

# ------------------------------------------------------------------------------------------------

function _preprocess_subplot(subplt::Subplot, d::Dict, args = ())
  validateSubplotSupported()
  preprocessArgs!(d)

  # for plotting recipes, swap out the args and update the parameter dictionary
  args = _apply_recipe(d, args...; d..., issubplot=true)

  dumpdict(d, "After subplot! preprocessing")

  # get the full plotargs, overriding any new settings
  # TODO: subplt.plotargs should probably be merged sooner and actually used
  #       for color selection, etc.  (i.e. if we overwrite the subplot palettes to [:heat :rainbow])
  #       then we need to overwrite plt[1].plotargs[:color_palette] to :heat before it's actually used
  #       for color selection!

  # first merge the new args into the subplot's plotargs. then process the plot args and merge
  # those into the plot's plotargs.  (example... `palette = [:blues :reds]` goes into subplt.plotargs,
  # then the ColorGradient for :blues/:reds is merged into plot 1/2 plotargs, which is then used for color selection)
  for i in 1:length(subplt.layout)
    subplt.plts[i].plotargs = getPlotArgs(backend(), merge(subplt.plts[i].plotargs, d), i)
  end
  merge!(subplt.plotargs, d)

  # process links.  TODO: extract to separate function
  for s in (:linkx, :linky, :linkfunc)
    if haskey(d, s)
      setfield!(subplt, s, d[s])
      delete!(d, s)
    end
  end

  args
end

function _postprocess_subplot(subplt::Subplot, d::Dict)
  # init (after plot creation)
  if !subplt.initialized
    subplt.initialized = _create_subplot(subplt, false)
  end

  # add title, axis labels, ticks, etc
  for (i,plt) in enumerate(subplt.plts)
    di = plt.plotargs
    dumpdict(di, "Updating sp $i")
    _update_plot(plt, di)
  end

  _update_plot_pos_size(subplt, d)

  # handle links
  subplt.linkx && link_axis(subplt, true)
  subplt.linky && link_axis(subplt, false)

  # set this to be current
  current(subplt)
end

# ------------------------------------------------------------------------------------------------

"""
Adds to a subplot.
"""

# current subplot
function subplot!(args...; kw...)
  validateSubplotSupported()
  subplot!(current(), args...; kw...)
end

# not allowed:
function subplot!(plt::Plot, args...; kw...)
  error("Can't call subplot! on a Plot!")
end

# # this adds to a specific subplot... most plot commands will flow through here
function subplot!(subplt::Subplot, args...; kw...)
  # validateSubplotSupported()

  d = Dict(kw)
  args = _preprocess_subplot(subplt, d, args)

  # create the underlying object (each backend will do this differently)
  # note: we call it once before doing the individual plots, and once after
  #       this is because some backends need to set up the subplots and then plot,
  #       and others need to do it the other way around
  if !subplt.initialized
    subplt.initialized = _create_subplot(subplt, true)
  end

  # handle grouping
  group = get(d, :group, nothing)
  if group == nothing
    groupargs = []
  else
    groupargs = [extractGroupArgs(d[:group], args...)]
    delete!(d, :group)
  end

  kwList, xmeta, ymeta = createKWargsList(subplt, groupargs..., args...; d...)

  # TODO: something useful with meta info?

  for (i,di) in enumerate(kwList)

    subplt.n += 1
    plt = getplot(subplt)
    plt.n += 1

    # cleanup the dictionary that we pass into the plot! command
    di[:show] = false
    di[:subplot] = true
    for k in (:title, :xlabel, :xticks, :xlims, :xscale, :xflip,
                      :ylabel, :yticks, :ylims, :yscale, :yflip)
      delete!(di, k)
    end
    dumpdict(di, "subplot! kwList $i")
    dumpdict(plt.plotargs, "plt.plotargs before plotting")

    _add_series_subplot(plt; di...)
  end

  _postprocess_subplot(subplt, d)

  # show it automatically?
  if haskey(d, :show) && d[:show]
    gui()
  end

  subplt
end

function _add_series_subplot(plt::Plot, args...; kw...)
  d = Dict(kw)

  setTicksFromStringVector(d, d, :x, :xticks)
  setTicksFromStringVector(d, d, :y, :yticks)

  _add_series(plt.backend, plt; d...)

  _add_annotations(plt, d)
  warnOnUnsupportedScales(plt.backend, d)
end

# --------------------------------------------------------------------

function Base.copy(subplt::Subplot)
  subplot(subplt.plts, subplt.layout, subplt.plotargs)
end


text:
export free_device, free_devices

function gethandle(dev::Int)
    ref = Ref{Ptr{Cvoid}}()
    @nvml :nvmlDeviceGetHandleByIndex (Cuint,Ptr{Ptr{Cvoid}}) dev ref
    ref[]
end

function getcount()
    ref = Ref{Cuint}()
    @nvml :nvmlDeviceGetCount (Ptr{Cuint},) ref
    Int(ref[])
end

"""
    free_devices(maxcount::Int)

Returns GPU device ids with no running processes
"""
function free_devices(maxcount::Int)
    devs = Int[]
    for i = 0:ndevices()-1
        length(devs) >= maxcount && break
        h = gethandle(i)
        res = @nvml_nocheck(:nvmlDeviceGetComputeRunningProcesses,
            (Ptr{Cvoid},Ptr{Cuint},Ptr{Cvoid}),
            h, Ref{Cuint}(0), C_NULL)
        (res == 0 || res == 7) && push!(devs,i)
    end
    devs
end
free_devices() = free_devices(ndevices())
function free_device()
    devs = free_devices(1)
    isempty(devs) && throw("No free CUDA device.")
    devs[1]
end


text:
function print(io::IO, exprs::Vector{Expr})
    for ex in exprs
        ex = ExprPrettify.prettify(ex)
        if ex.head === :block
            for a in ex.args
                if !isnothing(a)
                    println(io, a)
                end
            end
        else
            println(io, ex)
        end
        println(io)
    end
    return
end


text:
mutable struct Iterator
    val::Any
    st::Any
end
Iterator(x) = Iterator(x, nothing)
Base.length(x::Iterator) = length(x.val)

const pyjlitertype = pynew()

function pyjliter_next(self::Iterator)
    val = self.val
    st = self.st
    if st === nothing
        z = iterate(val)
    else
        z = iterate(val, something(st))
    end
    if z === nothing
        errset(pybuiltins.StopIteration)
        pynew()
    else
        r, newst = z
        self.st = Some(newst)
        Py(r)
    end
end

function init_jlwrap_iter()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class IteratorValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __iter__(self):
            return self
        def __next__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjliter_next)))
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlitertype, jl.IteratorValue)
end

pyjl(v::Iterator) = pyjl(pyjlitertype, v)


text:
########################################################
#   EXPERIMENTO SINTÉTICO DE NUESTRO algoritmo
# Visualiza para ciertos tamaños de muestra el error obtenido
########################################################
using Random
using Plots
using TOML
FICHERO_CONFIGURACION = "Experimentos/.config.toml"
config = TOML.parsefile(FICHERO_CONFIGURACION)["visualizacion-inicializacion-pesos-R"]
img_path = config["DIRECTORIO_IMAGENES"]

Random.seed!(1)
include("../../OptimizedNeuralNetwork.jl/src/OptimizedNeuralNetwork.jl")
using .OptimizedNeuralNetwork

M = 1
K_range = 3
f_regression(x)=(x<1) ? exp(-x)-4 : log(x)
for (data_set_size,n) in zip([3,4,5, 8,15,23,51,73,100, 103],[2,3,5,7,10,20,51,72,90, 100])

    println("EXPERIMENTO SINTÉTICO")
    println("n=$n y tamaño conjunto $data_set_size")
    # Partición del conjunto de muestra
    X_train= Vector(LinRange(-K_range, K_range, n))
    Y_train = map(f_regression, X_train)
    # Cálculo de la red neuronal con pesos inicializados
    h = nn_from_data(X_train, Y_train, n, M)
    # Función de evaluación por forward propagation
    evaluate(x)=forward_propagation(h,
        RampFunction,x)
    # Visualización
    interval = [-K_range,K_range]
    file_name = "f_ideal_y_rn_con_$(n)_neuronas"
    plot(x->evaluate([x])[1],
        -K_range,K_range,
        label="red neuronal n=$n"
    )
    plot!(f_regression,
        label="f ideal",
        title="Comparativa función ideal y red neuronal n=$n"
    )
    png(img_path*file_name)
end


text:
function load_trigger_values(experimental_paradigm::String)

    if experimental_paradigm == "regsoi"
        # Keys and time of soi (s)
        triggers = Dict(
        "1" => 1.0,
        "2" => 0.3,
        "3" => 0.4,
        "4" => 0.5,
        "5" => 0.6,
        "6" => 0.7,
        "7" => 0.8,
        "8" => 0.9,
        "9" => 1.1,
        "10" => 1.2,
        "11" => 1.3,
        "12" => 1.4,
        "13" => 1.5,
        "14" => 1.6,
        "15" => 1.7,
        "16" => 1.8,
        "17" => 1.9,
        "18" => 2.0,
        "19" => 1.001, # a second round of 1s to see an effect of fatigue
        "20" => "Noise",
    )

    elseif experimental_paradigm == "soistream"
        # Keys (in ms) and time of soi (s)
        triggers = Dict(
        "250" => 0.25,
        "500" => 0.5,
        "750" => 0.75,
        "1000" => 1.0,
        "1500" => 1.5,
        "2000" => 2.0,
        "3000" => 3.0,
        "4000" => 4.0,
        "5000" => 5.0,
        "7000" => 7.0,
    )

    else
    @error "Experimental paradigm not found, make yours and add it into src/triggerkeys.jl"
    end
    return triggers

end

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
"""
    The function that will modulate tha activation value of
    each non-input neuron in the network. It maps the z vale (∑wx + b)
    of the neurons to ℝ. By default its taken the sigmoid function
    one(z)/(one(z) + exp(-z)). To change it, just redefine it!!!
    But remember change the act_function_prime function too!!!!!!
"""
function act_function(z::T)::T where T<:Real
    return sigmoid(z);
end

"""
    The first derivative of the act_function.
"""
function act_function_prime(z::T) where T<:Real
    return sigmoid_prime(z);
end

function sigmoid(z::T)::T where T<:Real
    return one(z)/(one(z) + exp(-z));
end
function sigmoid_prime(z::T) where T<:Real
    return exp(-z)/((one(z) + exp(-z))^2);
end

ReLU_ng = 0.1;

function ReLU(z::T)::T where T<:Real
    return max(ReLU_ng * z, z);
end

function ReLU_prime(z::T)::T where T<:Real
    if z < 0
        return ReLU_ng;
    else
        return one(z);
    end
end


text:
# Linear model equations, for ocean model implicit vertical diffusion
# convective adjustment step.
#
using ClimateMachine.DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder, NumericalFluxSecondOrder, NumericalFluxGradient

using ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient, CentralNumericalFluxSecondOrder

using LinearAlgebra: I, dot, Diagonal

using ClimateMachine.BalanceLaws:
   BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux

using ClimateMachine.VariableTemplates

import ClimateMachine.DGMethods:
     init_state_auxiliary!, update_auxiliary_state!, update_auxiliary_state_gradient!, vars_state, VerticalDirection, boundary_state!, compute_gradient_flux!, init_state_prognostic!, flux_first_order!, flux_second_order!, source!, wavespeed, compute_gradient_argument!

using StaticArrays
using Random

"""
 IVDCModel{M} <: BalanceLaw

 This code defines DG `BalanceLaw` terms for an operator, L, that is evaluated from iterative
 implicit solver to solve an equation of the form

 (L + 1/Δt) ϕ^{n+1} = ϕ^{n}/Δt

  where L is a vertical diffusion operator with a spatially varying diffusion
  coefficient.

 # Usage

"""

# Create a new linear model instance
abstract type AbstractIVDCModel <: BalanceLaw end
struct IVDCModel{FT} <: AbstractIVDCModel
 dt::FT
 κᶻ::FT
 κᶜ::FT
 cʰ::FT
 cᶻ::FT
 gconf
 function IVDCModel{FT}(
  ;
  dt = FT(5400),
  κᶻ = FT(1e-4),
  κᶜ = FT(0.1),
  cʰ = FT(0),
  cᶻ = FT(0),
  gconf=()
 ) where {FT<: AbstractFloat}
    return new{FT}(
      dt,
      κᶻ,
      κᶜ,
      cʰ,
      cᶻ,
      gconf
    )
 end
end

function IVDCDGModel(
   bl::IVDCModel,
   grid,
   nfnondiff,
   nfdiff,
   gnf;
   kwargs...,
   )

   modeldata=(dt=bl.dt,κᶻ=bl.κᶻ,κᶜ=bl.κᶜ,cʰ=bl.cʰ,cᶻ=bl.cᶻ,gconf=bl.gconf)

   return DGModel(bl,grid,nfnondiff,nfdiff,gnf;kwargs...,modeldata=modeldata,)
end

"""
 Set model state variables and operators
"""

# State variable and initial value, just one for now, θ
##
vars_state(m::IVDCModel, ::Prognostic, FT) = @vars(θ::FT)

init_state_prognostic!( args...) = ( init_state_cp!( args... ) )
function init_state_cp!( m::IVDCModel, Q::Vars, A::Vars, coords, t,)
  @inbounds begin
    x = coords[1]
    y = coords[2]
    z = coords[3]
    # Wavy with stable vertical stratification
    # Q.θ = (5 + 4 * cos(y * π / p.Lʸ)) * (1 + z / H)

    # Constant everywhere
    # Q.θ = 9

    # Big number
    # Q.θ = 50

    # Unstable column
    #
    # Define function with inversion near surface and stable stratification below
    #
    # $ \theta(z) = Ae^{\frac{-(z+L)}{L}} - Ce^{\frac{-D(z+L)}{L}} + B + Ez$
    #
    H₀=m.gconf.H/4
    H₁=m.gconf.H*2
    ΔH=H₁-H₀
    sx=x/m.gconf.Lˣ
    sy=y/m.gconf.Lʸ
    D₀=2.
    D₁=3.
    ΔD=D₁-D₀
    # rz=rand(length(z)).-0.5
    rz=rand()-0.5
    H=H₀+sx*ΔH
    D=D₀+sy*ΔD
    L=H/10.;A=20.;C=50.;B=8.;E=5.e-4;
    ft(xx,L)=exp(-xx/L);
    th1(zz)=A*ft.(-zz .+ L, L);
    th2(zz)=C*ft(D*(-zz .+ L), L);
    phi1=th1.(z)
    phi2=th2.(z)
    val=phi1 .- phi2 .+ B .+ E .* z;
    # Q.θ = val .+ 0.005.*( val.*rz )
    Q.θ = val + 0.005*val*rz
  end
  return nothing
end
##

##
vars_state(m::IVDCModel, ::Auxiliary, FT) = @vars(θ_init::FT)

function init_state_auxiliary!(m::IVDCModel,A::Vars, _...)
  @inbounds begin
    A.θ_init = -0
  end
  return nothing
end
##
#

# Variables and operations used in differentiating first derivatives
##
vars_state(m::IVDCModel, ::Gradient, FT) = @vars(∇θ::FT, ∇θ_init::FT,)
@inline function compute_gradient_argument!(
    m::IVDCModel,
    G::Vars,
    Q::Vars,
    A,
    t,
)
    G.∇θ = Q.θ
    G.∇θ_init = A.θ_init

    return nothing
end
##
#

# Variables and operations used in differentiating second derivatives
##
vars_state(m::IVDCModel, ::GradientFlux, FT) = @vars(κ∇θ::SVector{3, FT})
@inline function compute_gradient_flux!(
    m::IVDCModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    κ = diffusivity_tensor(m, G.∇θ_init[3])
    D.κ∇θ = -κ * G.∇θ
    return nothing
end
##

##
## Set vertical diffusivity profile based on vertical hydrography profile
@inline function diffusivity_tensor(m::IVDCModel, ∂θ∂z)
    κᶻ = m.κᶻ
    κᶜ = m.κᶜ
    ∂θ∂z < 0 ? κ = (@SVector [0, 0, κᶜ]) : κ = (@SVector [0, 0, κᶻ])
    return Diagonal(-κ)
end
##
#

# Function to apply I to state variable

##
@inline function source!(
    m::IVDCModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t,
    direction,
)
    ivdc_dt = m.dt
    @inbounds begin
     S.θ = Q.θ/ivdc_dt
    end

    return nothing
end
##
#

# Numerical fluxes and boundaries

##
function flux_first_order!(::IVDCModel, _...) end

function flux_second_order!(
    ::IVDCModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    H::Vars,
    A::Vars,
    t,
)
    F.θ += D.κ∇θ
end

function wavespeed(m::IVDCModel, n⁻, _...)
    C = abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)
    return C
end
##

function boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient, CentralNumericalFluxGradient},
    m::IVDCModel,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    bctype,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    return nothing
end

###    From -  function numerical_boundary_flux_gradient! , DGMethods/NumericalFluxes.jl
###    boundary_state!(
###        numerical_flux,
###        balance_law,
###        state_conservative⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        aux1⁻,
###    )

function boundary_state!(
    nf::Union{NumericalFluxSecondOrder,CentralNumericalFluxSecondOrder},
    m::IVDCModel,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    bctype,
    t,
    _...,
)
    Q⁺.θ = Q⁻.θ
    D⁺.κ∇θ = n⁻ * -0
    # D⁺.κ∇θ = Q/(rho*cp)   ( at top boundary for heat flux )
    return nothing
end

###    boundary_state!(
###        numerical_flux,
###        balance_law,
###        state_conservative⁺,
###        state_gradient_flux⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_gradient_flux⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        diff1⁻,
###        aux1⁻,
###    )

###    boundary_flux_second_order!(
###        numerical_flux,
###        balance_law,
###        Grad{S}(flux),
###        state_conservative⁺,
###        state_gradient_flux⁺,
###        state_hyperdiffusive⁺,
###        state_auxiliary⁺,
###        normal_vector,
###        state_conservative⁻,
###        state_gradient_flux⁻,
###        state_hyperdiffusive⁻,
###        state_auxiliary⁻,
###        bctype,
###        t,
###        state1⁻,
###        diff1⁻,
###        aux1⁻,
###    )

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
rand(rng::AbstractRNG, ::Random.SamplerType{E}) where {T<:AbstractFloat, E<:RealLog{T}} =
	exp(E, -randexp(rng, T))
rand(rng::AbstractRNG, ::Random.SamplerType{E}) where {E<:RealLog} =
	exp(E, -randexp(rng))

# todo: sample CLogarithmic
# (note that rand(Complex) samples uniformly from the unit square)


text:
###############################################################################
#
#   AbsSeries.jl : Power series over rings, capped relative precision
#
###############################################################################

export O, valuation, exp, precision, max_precision, set_prec!

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

doc"""
    O{T <: RingElement}(a::Nemo.AbsSeriesElem{T})
> Returns $0 + O(x^\mbox{deg}(a))$. Usually this function is called with $x^n$
> as parameter. Then the function returns the power series $0 + O(x^n)$, which
> can be used to set the precision of a power series when constructing it.
"""
function O(a::Nemo.AbsSeriesElem{T}) where T <: RingElement
   if iszero(a)
      return deepcopy(a)    # 0 + O(x^n)
   end
   prec = length(a) - 1
   prec < 0 && throw(DomainError())
   return parent(a)(Array{T}(0), 0, prec)
end

parent_type(::Type{AbsSeries{T}}) where {T <: RingElement} = AbsSeriesRing{T}

elem_type(::Type{AbsSeriesRing{T}}) where {T <: RingElement} = AbsSeries{T}

###############################################################################
#
#   Basic manipulation
#
###############################################################################

length(x::Nemo.AbsSeriesElem) = x.length

precision(x::Nemo.AbsSeriesElem) = x.prec

doc"""
    max_precision(R::SeriesRing)
> Return the maximum absolute precision of power series in the given power
> series ring.
"""
max_precision(R::AbsSeriesRing) = R.prec_max

function normalise(a::AbsSeries, len::Int)
   while len > 0 && iszero(a.coeffs[len])
      len -= 1
   end
   return len
end

function coeff(a::AbsSeries, n::Int)
   n < 0  && throw(DomainError())
   return n >= length(a) ? zero(base_ring(a)) : a.coeffs[n + 1]
end

doc"""
    gen{T <: RingElement}(R::AbsSeriesRing{T})
> Return the generator of the power series ring, i.e. $x + O(x^n)$ where
> $n$ is the precision of the power series ring $R$.
"""
function gen(R::AbsSeriesRing{T}) where T <: RingElement
   S = base_ring(R)
   return R([S(0), S(1)], 2, max_precision(R))
end

doc"""
    iszero(a::SeriesElem)
> Return `true` if the given power series is arithmetically equal to zero to
> its current precision, otherwise return `false`.
"""
iszero(a::SeriesElem) = length(a) == 0

doc"""
    isone(a::AbsSeries)
> Return `true` if the given power series is arithmetically equal to one to
> its current precision, otherwise return `false`.
"""
function isone(a::AbsSeries)
   return (length(a) == 1 && isone(coeff(a, 0))) || precision(a) == 0
end

doc"""
    isgen(a::AbsSeries)
> Return `true` if the given power series is arithmetically equal to the
> generator of its power series ring to its current precision, otherwise return
> `false`.
"""
function isgen(a::AbsSeries)
   return (valuation(a) == 1 && length(a) == 2 && isone(coeff(a, 1))) ||
           precision(a) == 0
end

doc"""
    isunit(a::Nemo.AbsSeriesElem)
> Return `true` if the given power series is arithmetically equal to a unit,
> i.e. is invertible, otherwise return `false`.
"""
isunit(a::Nemo.AbsSeriesElem) = valuation(a) == 0 && isunit(coeff(a, 0))

doc"""
    valuation(a::Nemo.AbsSeriesElem)
> Return the valuation of the given power series, i.e. the degree of the first
> nonzero term (or the precision if it is arithmetically zero).
"""
function valuation(a::Nemo.AbsSeriesElem)
   for i = 1:length(a)
      if !iszero(coeff(a, i - 1))
         return i - 1
      end
   end
   return precision(a)
end

function deepcopy_internal(a::AbsSeries{T}, dict::ObjectIdDict) where {T <: RingElement}
   coeffs = Array{T}(length(a))
   for i = 1:length(a)
      coeffs[i] = deepcopy(coeff(a, i - 1))
   end
   return parent(a)(coeffs, length(a), precision(a))
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::Nemo.AbsSeriesElem)
   len = length(x)

   if len == 0
      print(io, zero(base_ring(x)))
   else
      coeff_printed = false
      for i = 0:len - 1
         c = coeff(x, i)
         if !iszero(c)
            if coeff_printed
               print(io, "+")
            end
            if i != 0
               if !isone(c)
                  print(io, "(")
                  print(io, c)
                  print(io, ")")
                  if i != 0
                     print(io, "*")
                  end
               end
               print(io, string(var(parent(x))))
               if i != 1
                  print(io, "^")
                  print(io, i)
               end
            else
               print(io, c)
            end
            coeff_printed = true
         end
      end
   end
   print(io, "+O(", string(var(parent(x))), "^", precision(x), ")")
end

###############################################################################
#
#   Unary operators
#
###############################################################################

doc"""
    -(a::Nemo.AbsSeriesElem)
> Return $-a$.
"""
function -(a::Nemo.AbsSeriesElem)
   len = length(a)
   z = parent(a)()
   set_prec!(z, precision(a))
   fit!(z, len)
   for i = 1:len
      z = setcoeff!(z, i - 1, -coeff(a, i - 1))
   end
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

doc"""
    +{T <: RingElement}(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T})
> Return $a + b$.
"""
function +(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   prec = min(precision(a), precision(b))
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = max(lena, lenb)
   z = parent(a)()
   fit!(z, lenz)
   set_prec!(z, prec)
   i = 1
   while i <= min(lena, lenb)
      z = setcoeff!(z, i - 1, coeff(a, i - 1) + coeff(b, i - 1))
      i += 1
   end
   while i <= lena
      z = setcoeff!(z, i - 1, coeff(a, i - 1))
      i += 1
   end
   while i <= lenb
      z = setcoeff!(z, i - 1, coeff(b, i - 1))
      i += 1
   end
   set_length!(z, normalise(z, i - 1))
   return z
end

doc"""
    -{T <: RingElement}(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T})
> Return $a - b$.
"""
function -(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   prec = min(precision(a), precision(b))
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = max(lena, lenb)
   z = parent(a)()
   fit!(z, lenz)
   set_prec!(z, prec)
   i = 1
   while i <= min(lena, lenb)
      z = setcoeff!(z, i - 1, coeff(a, i - 1) - coeff(b, i - 1))
      i += 1
   end
   while i <= lena
      z = setcoeff!(z, i - 1, coeff(a, i - 1))
      i += 1
   end
   while i <= lenb
      z = setcoeff!(z, i - 1, -coeff(b, i - 1))
      i += 1
   end
   set_length!(z, normalise(z, i - 1))
   return z
end

doc"""
    *{T <: RingElement}(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T})
> Return $a\times b$.
"""
function *(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)

   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(precision(a) + bval, precision(b) + aval)
   prec = min(prec, max_precision(parent(a)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   if lena == 0 || lenb == 0
      return parent(a)(Array{T}(0), 0, prec)
   end
   t = base_ring(a)()
   lenz = min(lena + lenb - 1, prec)
   d = Array{T}(lenz)
   for i = 1:min(lena, lenz)
      d[i] = coeff(a, i - 1)*coeff(b, 0)
   end
   if lenz > lena
      for j = 2:min(lenb, lenz - lena + 1)
          d[lena + j - 1] = coeff(a, lena - 1)*coeff(b, j - 1)
      end
   end
   for i = 1:lena - 1
      if lenz > i
         for j = 2:min(lenb, lenz - i + 1)
            t = mul!(t, coeff(a, i - 1), coeff(b, j - 1))
            d[i + j - 1] = addeq!(d[i + j - 1], t)
         end
      end
   end
   z = parent(a)(d, lenz, prec)
   set_length!(z, normalise(z, lenz))
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

doc"""
    *{T <: RingElem}(a::T, b::Nemo.AbsSeriesElem{T})
> Return $a\times b$.
"""
function *(a::T, b::Nemo.AbsSeriesElem{T}) where {T <: RingElem}
   len = length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*coeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   return z
end

doc"""
    *(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.AbsSeriesElem)
> Return $a\times b$.
"""
function *(a::Union{Integer, Rational, AbstractFloat}, b::Nemo.AbsSeriesElem)
   len = length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*coeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   return z
end

doc"""
    *{T <: RingElem}(a::Nemo.AbsSeriesElem{T}, b::T)
> Return $a\times b$.
"""
*(a::Nemo.AbsSeriesElem{T}, b::T) where {T <: RingElem} = b*a

doc"""
    *(a::Nemo.AbsSeriesElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a\times b$.
"""
*(a::Nemo.AbsSeriesElem, b::Union{Integer, Rational, AbstractFloat}) = b*a

###############################################################################
#
#   Shifting
#
###############################################################################

doc"""
    shift_left(x::Nemo.AbsSeriesElem, n::Int)
> Return the power series $f$ shifted left by $n$ terms, i.e. multiplied by
> $x^n$.
"""
function shift_left(x::Nemo.AbsSeriesElem{T}, len::Int) where {T <: RingElement}
   len < 0 && throw(DomainError())
   xlen = length(x)
   prec = precision(x) + len
   prec = min(prec, max_precision(parent(x)))
   if xlen == 0
      z = zero(parent(x))
      set_prec!(z, prec)
      return z
   end
   zlen = min(prec, xlen + len)
   z = parent(x)()
   fit!(z, zlen)
   set_prec!(z, prec)
   for i = 1:len
      z = setcoeff!(z, i - 1, zero(base_ring(x)))
   end
   for i = 1:xlen
      z = setcoeff!(z, i + len - 1, coeff(x, i - 1))
   end
   set_length!(z, normalise(z, zlen))
   return z
end

doc"""
    shift_right(f::Nemo.AbsSeriesElem, n::Int)
> Return the power series $f$ shifted right by $n$ terms, i.e. divided by
> $x^n$.
"""
function shift_right(x::Nemo.AbsSeriesElem{T}, len::Int) where {T <: RingElement}
   len < 0 && throw(DomainError())
   xlen = length(x)
   if len >= xlen
      z = zero(parent(x))
      set_prec!(z, max(0, precision(x) - len))
      return z
   end
   z = parent(x)()
   fit!(z, xlen - len)
   set_prec!(z, precision(x) - len)
   for i = 1:xlen - len
      z = setcoeff!(z, i - 1, coeff(x, i + len - 1))
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

doc"""
    truncate(a::Nemo.AbsSeriesElem, n::Int)
> Return $a$ truncated to $n$ terms.
"""
function truncate(a::Nemo.AbsSeriesElem{T}, prec::Int) where {T <: RingElement}
   prec < 0 && throw(DomainError())
   len = length(a)
   if precision(a) <= prec
      return a
   end
   z = parent(a)()
   fit!(z, prec)
   set_prec!(z, prec)
   for i = 1:min(prec, len)
      z = setcoeff!(z, i - 1, coeff(a, i - 1))
   end
   for i = len + 1:prec
      z = setcoeff!(z, i - 1, zero(base_ring(a)))
   end
   set_length!(z, normalise(z, prec))
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

doc"""
    ^{T <: RingElement}(a::Nemo.AbsSeriesElem{T}, b::Int)
> Return $a^b$. We require $b \geq 0$.
"""
function ^(a::Nemo.AbsSeriesElem{T}, b::Int) where {T <: RingElement}
   b < 0 && throw(DomainError())
   # special case powers of x for constructing power series efficiently
   if b == 0
      z = one(parent(a))
      set_prec!(z, precision(a))
      return z
   elseif precision(a) > 0 && isgen(a) && b > 0
      return shift_left(a, b - 1)
   elseif length(a) == 1
      z = parent(a)(coeff(a, 0)^b)
      set_prec!(z, precision(a))
      return z
   elseif b == 1
      return deepcopy(a)
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit !=0
         z = z*z
         if (UInt(bit) & b) != 0
            z *= a
         end
         bit >>= 1
      end
      return z
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

doc"""
    =={T <: RingElement}(x::Nemo.AbsSeriesElem{T}, y::Nemo.AbsSeriesElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`. Recall
> that power series to different precisions may still be arithmetically
> equal to the minimum of the two precisions.
"""
function ==(x::Nemo.AbsSeriesElem{T}, y::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   check_parent(x, y)
   prec = min(precision(x), precision(y))
   m1 = min(length(x), length(y))
   m2 = max(length(x), length(y))
   m1 = min(m1, prec)
   m2 = min(m2, prec)
   if length(x) >= m2
      for i = m1 + 1: m2
         if !iszero(coeff(x, i - 1))
            return false
          end
      end
   else
      for i = m1 + 1: m2
         if !iszero(coeff(y, i - 1))
            return false
          end
      end
   end
   for i = 1:m1
      if coeff(x, i - 1) != coeff(y, i - 1)
         return false
      end
   end
   return true
end

doc"""
    isequal{T <: RingElement}(x::Nemo.AbsSeriesElem{T}, y::Nemo.AbsSeriesElem{T})
> Return `true` if $x == y$ exactly, otherwise return `false`. Only if the
> power series are precisely the same, to the same precision, are they declared
> equal by this function.
"""
function isequal(x::Nemo.AbsSeriesElem{T}, y::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   if parent(x) != parent(y)
      return false
   end
   if precision(x) != precision(y) || length(x) != length(y)
      return false
   end
   for i = 1:length(x)
      if !isequal(coeff(x, i - 1), coeff(y, i - 1))
         return false
      end
   end
   return true
end

###############################################################################
#
#   Approximation
#
###############################################################################

function Base.isapprox(f::Nemo.AbsSeriesElem, g::Nemo.AbsSeriesElem; atol::Real=sqrt(eps()))
   check_parent(f, g)
   nmin = min(precision(f), precision(g))
   i = 1
   while i <= nmin
      if !isapprox(coeff(f, i - 1), coeff(g, i - 1); atol=atol)
         return false
      end
      i += 1
   end
   return true
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

doc"""
    =={T <: RingElem}(x::Nemo.AbsSeriesElem{T}, y::T)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::Nemo.AbsSeriesElem{T}, y::T) where {T <: RingElem} = precision(x) == 0 ||
      ((length(x) == 0 && iszero(y)) || (length(x) == 1 && coeff(x, 0) == y))

doc"""
    =={T <: RingElem}(x::T, y::Nemo.AbsSeriesElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::T, y::Nemo.AbsSeriesElem{T}) where {T <: RingElem} = y == x

doc"""
    ==(x::Nemo.AbsSeriesElem, y::Union{Integer, Rational, AbstractFloat})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::Nemo.AbsSeriesElem, y::Union{Integer, Rational, AbstractFloat}) = precision(x) == 0 || ((length(x) == 0 && iszero(y))
                                       || (length(x) == 1 && coeff(x, 0) == y))

doc"""
    ==(x::Union{Integer, Rational, AbstractFloat}, y::Nemo.AbsSeriesElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::Union{Integer, Rational, AbstractFloat}, y::Nemo.AbsSeriesElem) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

doc"""
    divexact{T <: RingElement}(a::Nemo.AbsSeriesElem{T}, b::Nemo.AbsSeriesElem{T})
> Return $a/b$. Requires $b$ to be invertible.
"""
function divexact(x::Nemo.AbsSeriesElem{T}, y::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   v2 = valuation(y)
   if v2 != 0
      v1 = valuation(x)
      if v1 >= v2
         x = shift_right(x, v2)
         y = shift_right(y, v2)
      end
   end
   y = truncate(y, precision(x))
   return x*inv(y)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

doc"""
    divexact(a::Nemo.AbsSeriesElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::Nemo.AbsSeriesElem, y::Union{Integer, Rational, AbstractFloat})
   y == 0 && throw(DivideError())
   lenx = length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(coeff(x, i - 1), y))
   end
   return z
end

doc"""
    divexact{T <: RingElem}(a::Nemo.AbsSeriesElem{T}, b::T)
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::Nemo.AbsSeriesElem{T}, y::T) where {T <: RingElem}
   iszero(y) && throw(DivideError())
   lenx = length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(coeff(x, i - 1), y))
   end
   return z
end

###############################################################################
#
#   Inversion
#
###############################################################################

doc"""
   inv(a::Nemo.AbsSeriesElem)
> Return the inverse of the power series $a$, i.e. $1/a$.
"""
function inv(a::Nemo.AbsSeriesElem)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   a1 = coeff(a, 0)
   ainv = parent(a)()
   fit!(ainv, precision(a))
   set_prec!(ainv, precision(a))
   if precision(a) != 0
      ainv = setcoeff!(ainv, 0, divexact(one(base_ring(a)), a1))
   end
   a1 = -a1
   for n = 2:precision(a)
      s = coeff(a, 1)*coeff(ainv, n - 2)
      for i = 2:min(n, length(a)) - 1
         s += coeff(a, i)*coeff(ainv, n - i - 1)
      end
      ainv = setcoeff!(ainv, n - 1, divexact(s, a1))
   end
   set_length!(ainv, normalise(ainv, precision(a)))
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

doc"""
    exp(a::Nemo.AbsSeriesElem)
> Return the exponential of the power series $a$.
"""
function Base.exp(a::Nemo.AbsSeriesElem)
   if iszero(a)
      z = one(parent(a))
      set_prec!(z, precision(a))
      return z
   end
   z = parent(a)()
   fit!(z, precision(a))
   set_prec!(z, precision(a))
   z = setcoeff!(z, 0, Nemo.exp(coeff(a, 0)))
   len = length(a)
   for k = 1 : precision(a) - 1
      s = zero(base_ring(a))
      for j = 1 : min(k + 1, len) - 1
         s += j * coeff(a, j) * coeff(z, k - j)
      end
      !isunit(base_ring(a)(k)) && error("Unable to divide in exp")
      z = setcoeff!(z, k, divexact(s, k))
   end
   set_length!(z, normalise(z, precision(a)))
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function fit!(c::AbsSeries{T}, n::Int) where {T <: RingElement}
   if length(c.coeffs) < n
      t = c.coeffs
      c.coeffs = Array{T}(n)
      for i = 1:c.length
         c.coeffs[i] = t[i]
      end
      for i = length(c) + 1:n
         c.coeffs[i] = zero(base_ring(c))
      end
   end
   return nothing
end

function setcoeff!(c::AbsSeries{T}, n::Int, a::T) where {T <: RingElement}
   if (!iszero(a) && precision(c) > n) || n + 1 <= c.length
      fit!(c, n + 1)
      c.coeffs[n + 1] = a
      c.length = max(length(c), n + 1)
      # don't normalise
   end
   return c
end

function mul!(c::AbsSeries{T}, a::AbsSeries{T}, b::AbsSeries{T}) where {T <: RingElement}
   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(precision(a) + bval, precision(b) + aval)
   prec = min(prec, max_precision(parent(c)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   if lena == 0 || lenb == 0
      c.length = 0
   else
      t = base_ring(a)()

      lenc = min(lena + lenb - 1, prec)
      fit!(c, lenc)

      for i = 1:min(lena, lenc)
         c.coeffs[i] = mul!(c.coeffs[i], coeff(a, i - 1), coeff(b, 0))
      end

      if lenc > lena
         for i = 2:min(lenb, lenc - lena + 1)
            c.coeffs[lena + i - 1] = mul!(c.coeffs[lena + i - 1], coeff(a, lena - 1), coeff(b, i - 1))
         end
      end

      for i = 1:lena - 1
         if lenc > i
            for j = 2:min(lenb, lenc - i + 1)
               t = mul!(t, coeff(a, i - 1), coeff(b, j - 1))
               c.coeffs[i + j - 1] = addeq!(c.coeffs[i + j - 1], t)
            end
         end
      end

      c.length = normalise(c, lenc)
   end
   c.prec = prec
   return c
end

function addeq!(c::AbsSeries{T}, a::AbsSeries{T}) where {T <: RingElement}
   lenc = length(c)
   lena = length(a)

   prec = min(precision(a), precision(c))

   lena = min(lena, prec)
   lenc = min(lenc, prec)

   len = max(lenc, lena)
   fit!(c, len)
   for i = 1:lena
      c.coeffs[i] = addeq!(c.coeffs[i], coeff(a, i - 1))
   end
   c.length = normalise(c, len)
   c.prec = prec
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{AbsSeries{T}}, ::Type{AbsSeries{T}}) where T <: RingElement = AbsSeries{T}

function promote_rule(::Type{AbsSeries{T}}, ::Type{U}) where {T <: RingElement, U <: RingElement}
   promote_rule(T, U) == T ? AbsSeries{T} : Union{}
end

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::AbsSeriesRing{T} where {T <: RingElement})(b::RingElement)
   return a(base_ring(a)(b))
end

function (a::AbsSeriesRing{T})() where {T <: RingElement}
   z = AbsSeries{T}(Array{T}(0), 0, a.prec_max)
   z.parent = a
   return z
end

function (a::AbsSeriesRing{T})(b::Union{Integer, Rational, AbstractFloat}) where {T <: RingElement}
   if b == 0
      z = AbsSeries{T}(Array{T}(0), 0, a.prec_max)
   else
      z = AbsSeries{T}([base_ring(a)(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::AbsSeriesRing{T})(b::T) where {T <: RingElem}
   parent(b) != base_ring(a) && error("Unable to coerce to power series")
   if b == 0
      z = AbsSeries{T}(Array{T}(0), 0, a.prec_max)
   else
      z = AbsSeries{T}([b], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::AbsSeriesRing{T})(b::Nemo.AbsSeriesElem{T}) where {T <: RingElement}
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::AbsSeriesRing{T})(b::Array{T, 1}, len::Int, prec::Int) where {T <: RingElement}
   if length(b) > 0
      parent(b[1]) != base_ring(a) && error("Unable to coerce to power series")
   end
   z = AbsSeries{T}(b, len, prec)
   z.parent = a
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

# see RelSeries.jl


text:
# @system LeafAppearance(Stage, Emergence, LeafInitiation) begin
#     LTARa_max: maximum_phyllochron_asymptote => 0.4421 ~ preserve(u"d^-1", parameter)
#     _SDm => 117.7523 ~ preserve(u"d", parameter)
#     _k => 0.0256 ~ preserve(u"d^-1", parameter)

#     LTAR_max0(LTARa_max, SD, _SDm, _k): initial_maximum_phyllochron => begin
#         LTARa_max / (1 + exp(-k * (SD - SDm)))
#     end ~ preserve(u"d^-1", parameter)

#     LTAR_max(LTAR_max0, LTARa_max, n=leaves_appeared, ng=leaves_generic): maximum_phyllochron => begin
#         n0 = 0
#         r0 = LTAR_max0
#         r1 = LTARa_max / 2
#         r0 + (r1 - r0) * (clamp(n, n0, ng) - n0) / (ng - n0)
#     end ~ track(u"d^-1")

### Jim modification of Leaf Appearance

@system LeafAppearance(Stage, Emergence, LeafInitiation) begin
    LTARa_max: maximum_phyllochron_asymptote => 0.4421 ~ preserve(u"d^-1", parameter)

    LTAR_max(LTARa_max, n=leaves_appeared, ng=leaves_generic): maximum_phyllochron => begin
        n0 = 0
        r0 = LTARa_max
        r1 = LTARa_max / 2
        r0 + (r1 - r0) * (clamp(n, n0, ng) - n0) / (ng - n0)
    end ~ track(u"d^-1")

    LTA(r=LTAR_max, β=BF.ΔT): leaf_tip_appearance => r*β ~ accumulate(when=leaf_appearing)

    leaf_appearable(emerged) ~ flag
    leaf_appeared(leaves_appeared, leaves_initiated) => begin
        #HACK ensure leaves are initiated
        leaves_appeared >= leaves_initiated > 0
    end ~ flag
    leaf_appearing(leaf_appearable & !leaf_appeared) ~ flag

    #HACK: take account of sprout leaf appeared in emergence stage
    ILTA: initial_leaf_tip_appearance => 1 ~ track::int(when=emerged)
    leaves_appeared(ILTA, LTA) => (ILTA + LTA) ~ track::int(round=:floor)
end


text:
module githubactions4julia

    # Write your package code here.
    function plus(x,y)
        return x+y
    end

    export plus
end


text:
function tr2!(r::Robot) #главная функция
    for side in (Nord, West, Sud, Ost) # - перебор всех возможных направлений
        move_corner(r, side)# перемещаем робота в угол
        put_a_marker(r, side)#ставим маркер
    end
    show!(r)
end
function put_a_marker(r::Robot, side::HorizonSide)#ставим маркер
    for side in (Ost, Sud, West, Nord)
        while isborder(r, side) == false
            putmarker!(r)
            move!(r, side)
        end
    end
end
function move_corner(r::Robot, side::HorizonSide) # перемещаем робота в угол
    while isborder(r, Nord) == false
        move!(r, Nord)
    end
end