text:
using Images, LinearAlgebra, Plots

function low_rank_approx(x::SVD{T}, k) where T <: Real
    #take a low rank (k) approx of the SVD of a matrix
    #via outer products
    return x.U[:, 1:k] * Diagonal(x.S[1:k]) * x.V[:, 1:k]'
end

image_path = joinpath(dirname(@__DIR__), "data/20210501_171236.jpg")

img_big = load(image_path)
img = Images.imresize(img_big, ratio = 1/4)

channels = channelview(img)
svd_r = svd(channels[1,:,:])
svd_g = svd(channels[2,:,:])
svd_b = svd(channels[3,:,:])

p1 = plot(svd_r.S, c = :red, label = "red singular values")
plot!(p1, svd_g.S, c = :green, label = "green singular values")
plot!(p1, svd_b.S, c = :blue, label = "blue singular values")

#recall trace norm identity

# ||A||f  = √sum(σ(A)^2)

#Frobenius norm => Pretend matrix is big vector and take l2 norm
# hints that truncated svd might contain a good representation on matrix

#in case you don't trust me...
norm(channels[1,:,:]) ≈ norm(svd_r.S)
norm(channels[2,:,:]) ≈ norm(svd_g.S)
norm(channels[3,:,:]) ≈ norm(svd_b.S)

approx_error = [norm(low_rank_approx(svd_r, i) - channels[1,:,:]) for i = 1 :10: length(svd_r.S)]

p2 = plot(approx_error./norm(channels[1,:,:]))

#what does it look like though...

reconstructed_images = []

for k = 1:10:60
    r_low = low_rank_approx(svd_r, k)
    g_low = low_rank_approx(svd_g, k)
    b_low = low_rank_approx(svd_b, k)
    push!(reconstructed_images, colorview(RGB, r_low, g_low, b_low))
end

output_plot = [plot(ri) for ri in reconstructed_images]

plot(output_plot...)

#So we can get a linear reduction in size by keeping the components of the
#outer product U[:, 1:k], S[1:k], V[:, 1:k] and reconstructing via matrix
#multiplication

#So...

@show sizeof(img_big)/1e6

#So to get approx 1% reconstruction error in Frobenius norm we need about half the singular values

percent_reduction_01 = findfirst(approx_error./norm(channels[1,:,:]) .< 0.01)/length(approx_error)

@show percent_reduction_01

#so the final "saved" image would be

@show (sizeof(img_big)/1e6) * percent_reduction_01

#not bad, but look at what jpeg gets...


text:
export MType, zerotype, montes_initial, phi_expansion, valuationf, extend_valuation

type MType
  phi
  f
  m
  F
  z
  NewtonOp
  ResOp
end

valuation{T <: Integer}(x::fmpz, y::T) = valuation(x, ZZ(y))

function phi_expansion(f::fmpz_poly, phi::fmpz_poly)
  (lead(f) != 1) && (lead(phi) != 1) && error("must be monic")
  q,r = pseudodivrem(f,phi)
  arr = [r]

  while degree(q) >= degree(phi)
    q,r = pseudodivrem(q,phi)
    push!(arr,r)
  end
  q,r = pseudodivrem(q,phi)
  push!(arr,r)
  return arr
end

function valuationf(p::Int)
  function v(x::fmpz)
    return valuation(x,p)[1]
  end
  return v
end

function extend_valuation(v::Function)
  function vext(f::fmpz_poly)
    a = Array(Int, length(f))
    for i in 1:length(f)
      a[i] = v(coeff(f,i-1))
    end
    return minimum(a)
  end
  return vext
end

function montes_initial(f::fmpz_poly, p::Int)
  kx,x = PolynomialRing(ResidueRing(ZZ, p), "x")
  fac = factor(kx(f))

  res = Array(MType, length(fac))

  for phi in keys(fac)

    F = NmodPolyRingResField(phi, :x)
    z = F(gen(parent(phi)))

    res[i] = MType(phi,degree(phi),degree(phi), F, z, 1, 1)
  end

  return res
end

function zerotype(f::nmod_poly)
end

ext: jl


text:
function get_data(  dataobject::HydroDataType,
                    vars::Array{Symbol,1},
                    units::Array{Symbol,1},
                    direction::Symbol,
                    center::Array{<:Real,1},
                    mask::MaskType,
                    ref_time::Real)

    boxlen = dataobject.boxlen
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)
    vars_dict = Dict()
    #vars = unique(vars)

    if direction == :z
        apos = :cx
        bpos = :cy
        cpos = :cz

        avel = :vx
        bvel = :vy
        cvel = :vz

    elseif direction == :y
        apos = :cz
        bpos = :cx
        cpos = :cy

        avel = :vz
        bvel = :vx
        cvel = :vy
    elseif direction == :x
        apos = :cz
        bpos = :cy
        cpos = :cx

        avel = :vz
        bvel = :vy
        cvel = :vx
    end

    column_names = propertynames(dataobject.data.columns)

    for i in vars

        # quantities that are in the datatable
        if in(i, column_names)

            selected_unit = getunit(dataobject, i, vars, units)
            if i == :cx
                if isamr
                    vars_dict[i] =  select(dataobject.data, apos) .- 2 .^getvar(dataobject, :level) .* center[1]
                else # if uniform grid
                    vars_dict[i] =  select(dataobject.data, apos) .- 2^lmax .* center[1]
                end
            elseif i == :cy
                if isamr
                    vars_dict[i] =  select(dataobject.data, bpos) .- 2 .^getvar(dataobject, :level) .* center[2]
                else # if uniform grid
                    vars_dict[i] =  select(dataobject.data, bpos) .- 2^lmax .* center[2]
                end
            elseif i == :cx
                if isamr
                    vars_dict[i] =  select(dataobject.data, cpos) .- 2 .^getvar(dataobject, :level) .* center[3]
                else # if uniform grid
                    vars_dict[i] =  select(dataobject.data, cpos) .- 2^lmax .* center[3]
                end
            else
                #if selected_unit != 1.
                    #println(i)
                    vars_dict[i] = select(dataobject.data, i) .* selected_unit
                #else
                    #vars_dict[i] = select(dataobject.data, i)
                #end
            end

        # quantities that are derived from the variables in the data table
        elseif i == :cellsize
            selected_unit = getunit(dataobject, :cellsize, vars, units)
            if isamr
                vars_dict[:cellsize] =  map(row-> dataobject.boxlen / 2^row.level * selected_unit , dataobject.data)
            else # if uniform grid
                vars_dict[:cellsize] =  map(row-> dataobject.boxlen / 2^lmax * selected_unit , dataobject.data)
            end
        elseif i == :volume
            selected_unit = getunit(dataobject, :volume, vars, units)
            vars_dict[:volume] =  convert(Array{Float64,1}, getvar(dataobject, :cellsize) .^3 .* selected_unit)

        elseif i == :jeanslength
            selected_unit = getunit(dataobject, :jeanslength, vars, units)
            vars_dict[:jeanslength] = getvar(dataobject, :cs, unit=:cm_s)  .*
                                        sqrt(3. * pi / (32. * dataobject.info.constants.G))  ./
                                        sqrt.( getvar(dataobject, :rho, unit=:g_cm3) ) ./ dataobject.info.scale.cm  .*  selected_unit
        elseif i == :jeansnumber
            selected_unit = getunit(dataobject, :jeansnumber, vars, units)
            vars_dict[:jeansnumber] = getvar(dataobject, :jeanslength) ./ getvar(dataobject, :cellsize) ./ selected_unit

        elseif i == :freefall_time
            selected_unit = getunit(dataobject, :freefall_time, vars, units)
            vars_dict[:freefall_time] = sqrt.( 3. * pi / (32. * dataobject.info.constants.G) ./ getvar(dataobject, :rho, unit=:g_cm3)  ) .* selected_unit

        elseif i == :mass
            selected_unit = getunit(dataobject, :mass, vars, units)
            vars_dict[:mass] =  getmass(dataobject) .* selected_unit

        elseif i == :cs
            selected_unit = getunit(dataobject, :cs, vars, units)
            vars_dict[:cs] =   sqrt.( dataobject.info.gamma .*
                                        select( dataobject.data, :p) ./
                                        select( dataobject.data, :rho) ) .* selected_unit

        elseif i == :vx2
            selected_unit = getunit(dataobject, :vx2, vars, units)
            vars_dict[:vx2] =  select(dataobject.data, :vx).^2  .* selected_unit.^2
        elseif i == :vy2
            selected_unit = getunit(dataobject, :vy2, vars, units)
            vars_dict[:vy2] =  select(dataobject.data, :vy).^2  .* selected_unit.^2
        elseif i == :vz2
            selected_unit = getunit(dataobject, :vz2, vars, units)
            vars_dict[:vz2] =  select(dataobject.data, :vz).^2  .* selected_unit.^2

        elseif i == :v
            selected_unit = getunit(dataobject, :v, vars, units)
            vars_dict[:v] =  sqrt.(select(dataobject.data, :vx).^2 .+
                                   select(dataobject.data, :vy).^2 .+
                                   select(dataobject.data, :vz).^2 ) .* selected_unit
        elseif i == :v2
           selected_unit = getunit(dataobject, :v2, vars, units)
           vars_dict[:v2] =      (select(dataobject.data, :vx).^2 .+
                                  select(dataobject.data, :vy).^2 .+
                                  select(dataobject.data, :vz).^2 ) .* selected_unit .^2

        elseif i == :vϕ_cylinder

            radius = getvar(dataobject, :r_cylinder, center=center)
            x = getvar(dataobject, :x, center=center)
            y = getvar(dataobject, :y, center=center)
            vx = getvar(dataobject, :vx)
            vy = getvar(dataobject, :vy)

            # selected_unit = getunit(dataobject, :vϕ, vars, units)
            # vϕ = (x .* vy .- y .* vx) ./ radius .* selected_unit
            # vϕ[isnan.(vϕ)] .= 0 # overwrite NaN due to radius = 0
            # vars_dict[:vϕ] = vϕ

            # vϕ = omega x radius
            selected_unit = getunit(dataobject, :vϕ_cylinder, vars, units)
            a = (-1 .* y) .^2 + x .^2
            b = ( x .* vy .- y .* vx) .^2
            vϕ_cylinder =  sqrt.( a .* b  ) ./ radius .^2 .* selected_unit
            #(y .* (y .* vx .- x .* vy) ).^2 .- ( x .* (y .* vx .- x .* vy) ) .^2

            #(x .* vy .- y .* vx) ./ radius .* selected_unit
            vϕ_cylinder[isnan.(vϕ_cylinder)] .= 0. # overwrite NaN due to radius = 0
            vars_dict[:vϕ_cylinder] = vϕ_cylinder

        elseif i == :vϕ_cylinder2
            #radius = getvar(dataobject, :r_cylinder, center=center)
            #x = getvar(dataobject, :x, center=center)
            #y = getvar(dataobject, :y, center=center)
            #vx = getvar(dataobject, :vx)
            #vy = getvar(dataobject, :vy)

            selected_unit = getunit(dataobject, :vϕ_cylinder2, vars, units)
            #vϕ2 = ((x .* vy .- y .* vx) ./ radius .* selected_unit ).^2
            #vϕ2[isnan.(vϕ2)] .= 0 # overwrite NaN due to radius = 0
            #vars_dict[:vϕ2] = vϕ2
            #vϕ_cylinder2 = ( sqrt.( (y .^2 .* (y .* vx .- x .* vy) .^2 ) .- ( x .^2 .* (y .* vx .- x .* vy) .^2 )  ) ./ radius .^2 .* selected_unit ) .^2
            #selected_unit = getunit(dataobject, :vϕ_cylinder2, vars, units)

            #vϕ_cylinder2 = ((x .* vy .- y .* vx) ./ radius .* selected_unit ).^2
            #vϕ_cylinder2[isnan.(vϕ_cylinder2)] .= 0 # overwrite NaN due to radius = 0
            vars_dict[:vϕ_cylinder2] = (getvar(dataobject, :vϕ_cylinder, center=center) .* selected_unit).^2

        elseif i == :vz2

            vz = getvar(dataobject, :vz)
            selected_unit = getunit(dataobject, :vz2, vars, units)
            vars_dict[:vz2] =  (vz .* selected_unit ).^2

        elseif i == :vr_cylinder

            radius = getvar(dataobject, :r_cylinder, center=center )

            x = getvar(dataobject, :x, center=center)
            y = getvar(dataobject, :y, center=center)
            vx = getvar(dataobject, :vx)
            vy = getvar(dataobject, :vy)

            selected_unit = getunit(dataobject, :vr_cylinder, vars, units)
            vr = (x .* vx .+ y .* vy) ./ radius .* selected_unit
            vr[isnan.(vr)] .= 0 # overwrite NaN due to radius = 0
            vars_dict[:vr_cylinder] =  vr

        elseif i == :vr_cylinder2

            selected_unit = getunit(dataobject, :vr_cylinder2, vars, units)
            vars_dict[:vr_cylinder2] = (getvar(dataobject, :vr_cylinder, center=center) .* selected_unit).^2

        elseif i == :x
            selected_unit = getunit(dataobject, :x, vars, units)
            if isamr
                vars_dict[:x] =  (getvar(dataobject, apos) .* boxlen ./ 2 .^getvar(dataobject, :level) .-  boxlen * center[1] )  .* selected_unit
            else # if uniform grid
                vars_dict[:x] =  (getvar(dataobject, apos) .* boxlen ./ 2^lmax .-  boxlen * center[1] )  .* selected_unit
            end
        elseif i == :y
            selected_unit = getunit(dataobject, :y, vars, units)
            if isamr
                vars_dict[:y] =  (getvar(dataobject, bpos) .* boxlen ./ 2 .^getvar(dataobject, :level) .- boxlen * center[2] )  .* selected_unit
            else # if uniform grid
                vars_dict[:y] =  (getvar(dataobject, bpos) .* boxlen ./ 2^lmax .- boxlen * center[2] )  .* selected_unit
            end
        elseif i == :z
            selected_unit = getunit(dataobject, :z, vars, units)
            if isamr
                vars_dict[:z] =  (getvar(dataobject, cpos) .* boxlen ./ 2 .^getvar(dataobject, :level) .- boxlen * center[3] )  .* selected_unit
            else # if uniform grid
                vars_dict[:z] =  (getvar(dataobject, cpos) .* boxlen ./ 2^lmax .- boxlen * center[3] )  .* selected_unit
            end

        elseif i == :mach #no unit needed
            vars_dict[:mach] = getvar(dataobject, :v) ./ getvar(dataobject, :cs)

        elseif i == :ekin
            selected_unit = getunit(dataobject, :ekin, vars, units)
            vars_dict[:ekin] =   0.5 .* getmass(dataobject)  .* getvar(dataobject, :v).^2 .* selected_unit

        elseif i == :r_cylinder
            selected_unit = getunit(dataobject, :r_cylinder, vars, units)
            if isamr
                vars_dict[:r_cylinder] = convert(Array{Float64,1}, select( dataobject.data, (apos, bpos, :level)=>p->
                                                selected_unit * sqrt( (p[apos] * boxlen / 2^p.level - boxlen * center[1] )^2 +
                                                                   (p[bpos] * boxlen / 2^p.level - boxlen * center[2] )^2 ) ) )
            else # if uniform grid

                vars_dict[:r_cylinder] = convert(Array{Float64,1}, select( dataobject.data, (apos, bpos)=>p->
                                                selected_unit * sqrt( (p[apos] * boxlen / 2^lmax - boxlen * center[1] )^2 +
                                                                   (p[bpos] * boxlen / 2^lmax - boxlen * center[2] )^2 ) ) )
            end
        elseif i == :r_sphere
            selected_unit = getunit(dataobject, :r_sphere, vars, units)
            if isamr
                vars_dict[:r_sphere] = select( dataobject.data, (apos, bpos, cpos, :level)=>p->
                                        selected_unit * sqrt( (p[apos] * boxlen / 2^p.level -  boxlen * center[1]  )^2 +
                                                               (p[bpos] * boxlen / 2^p.level -  boxlen * center[2] )^2  +
                                                               (p[cpos] * boxlen / 2^p.level -  boxlen * center[3] )^2 ) )
            else # if uniform grid
                vars_dict[:r_sphere] = select( dataobject.data, (apos, bpos, cpos)=>p->
                                        selected_unit * sqrt( (p[apos] * boxlen / 2^lmax -  boxlen * center[1]  )^2 +
                                                               (p[bpos] * boxlen / 2^lmax -  boxlen * center[2] )^2  +
                                                               (p[cpos] * boxlen / 2^lmax -  boxlen * center[3] )^2 ) )
            end

        end

    end

    if length(mask) > 1
        for i in keys(vars_dict)
            vars_dict[i]=vars_dict[i][mask]
        end
    end

    if length(vars)==1
            return vars_dict[vars[1]]
    else
            return vars_dict
    end

end


text:
using CrystalNets
using Test
using Random
using PeriodicGraphs
using StaticArrays
using Base.Threads

function stress_test(g::Union{AbstractString,PeriodicGraph3D}, num=0,
                     refgraph=topological_genome(PeriodicGraph3D(g)), skiplong=false)
    ref = string(refgraph)
    control = 50
    m::Int = if iszero(num)
        clamp(round(Int, begin
            _counter = 0
            timebefore = time_ns()
            while time_ns() - timebefore < 0.2e9
                _ref = string(topological_genome(refgraph))
                @assert _ref == ref
                _counter += 1
            end
            timeafter = time_ns()
            elapsed = (timeafter - timebefore)*1.0e-9 / _counter
            if elapsed < 0.5
                20 + 30 / elapsed
            elseif elapsed < 5
                control = 10
                10 + 19 / elapsed
            elseif elapsed < 30
                control = 3
                3 + 59 / elapsed
            else
                skiplong && throw("(skipped: over 30s)")
                control = 1
                120 / elapsed
            end
        end), 2, 200)
    else
        num
    end
    println("(m = ", m, ')')
    n = nv(refgraph)
    # graph = PeriodicGraph3D(refgraph)
    # counter = 0
    graphs = [PeriodicGraph3D(refgraph) for _ in 1:nthreads()]
    # timings = Float64[]
    counter = Atomic{Int}(1)
    #Juno.progress() do progressid
        @threads for i in 1:m
        # for i in 1:m
            id = threadid()
            # r = collect(1:n)
            # r[end], r[1], r[10] = r[10], r[end], r[1]
            r = randperm(n)
            # offsets = [zero(SVector{3,Int}) for _ in 1:n]
            offsets = [SVector{3,Int}([rand(-3:3) for _ in 1:3]) for _ in 1:n]
            graph = swap_axes!(offset_representatives!(graphs[id][r], offsets), randperm(3))
            graphs[id] = graph
            # tic = time_ns()
            res = string(topological_genome(graph))
            # tac = time_ns()
            if res != ref
                #@info "stress_test main loop" progress="done" _id=progressid
                throw(string("TEST FAILED: different representations for res and ref:\n",
                        "ref = \"", ref, "\"\nres = \"", res, '\"', '\n'))
            end
            iterations = atomic_add!(counter, 1)
            # counter += 1
            #if iterations % (1 + m÷50) == 0
            #    @info "stress_test main loop" progress=iterations/m _id=progressid
            #end
            if iterations % control == 0
            # if counter % control == 0
                GC.gc(true)
                # @info "stress_test main loop" progress=counter/m _id=progressid
            end
        end
        #@info "stress_test main loop" progress="done" _id=progressid
    #end
    # foreach(println, graphs)
    nothing
end
function stress_test(net::CrystalNet, m=nothing)
    stress_test(net.graph, m)
end

function firsttests(arc)
    g1 = PeriodicGraph3D("3 1 2 2 0 0 1 2 0 2 0 1 2 0 0 2 1 2 0 0 0")
    g2 = PeriodicGraph3D("3 1 2 0 0 0 1 2 0 0 1 1 2 0 1 0 1 2 1 0 0")
    @assert topological_genome(g1) == topological_genome(g2)
    nothing
end

function test_arc(path, newarc=nothing, onlyone=true, skipslows=false)
    flag, arc = CrystalNets.parse_arc(path)
    @assert flag
    firsttests(arc)
    timings = Dict{String,Float64}()
    timing = time_ns()
    for (key, id) in arc
        key[1] != '3' && continue
        if skipslows
            id ∈ ("tep", "fav", "moo-a", "tei", "gng", "mov", "nul", "lcz",
                      "soa", "mtn-e-a", "nun", "txt", "tcc-x", "unr", "fiv",
                      "sig") && continue
        end
        print("Testing ", id, ' ')
        try
            graph = try
                tic = time_ns()
                x = topological_genome(PeriodicGraph3D(key))
                tac = time_ns()
                timing = tac - tic
                x
            catch e
                if e == "UNTESTED"
                    println("(untested)")
                    continue
                end
                rethrow()
            end
            if !(isnothing(newarc))
                open(newarc, "a") do f
                    println(f, "key\t\t", string(graph))
                    println(f, "id\t\t", id)
                    println(f)
                end
            end
            timings[id] = timing * 1e-9
            if timing > 5e9
                GC.gc(true)
                ccall(:malloc_trim, Cvoid, (Cint,), 0)
            end
            if onlyone
                println("(m = 1)")
            else
                # append!(timings[id],
                    stress_test(key, 0, graph, false)
                # )
            end
        catch e
            if e isa String
                println("\"\"\"", e, "\"\"\"")
            elseif e isa InterruptException || e isa TaskFailedException && e.task.exception isa InterruptException
                println("INTERRUPTION CAUGHT")
                rethrow()
            elseif e isa TaskFailedException
                println("TaskFailedException: ", e.task.exception, '\n')
            elseif e isa AssertionError
                println("AssertionError: ", e, '\n')
            else
                println("\nUNKNOWN ERROR happened:")
                showerror(stdout, e)
                println('\n')
            end
        end
    end
    if !(isnothing(newarc))
        open(newarc, "a") do f
            println(f, '\n')
        end
    end
    return timings
end

test_arc(joinpath(@__DIR__, "..", "src", "default.arc"))


text:
using DemeNet: Signer, DemeSpec, Deme, save, Consensus, DemeID
using PeaceVote.BraidChains.Braiders: BraiderConfig, Braider, Mixer, braid!
using PeaceCypher

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:PeaceVote)
save(demespec) ### Necessary to connect with Mixer

deme = Deme(demespec)
uuid = demespec.uuid

mixer = Signer(deme,"mixer")
mixerserver = Mixer(1999,deme,mixer)

server = Signer(deme,"server")

MIXER_ID = mixer.id
SERVER_ID = server.id

config = BraiderConfig(1998,1999,UInt8(3),UInt8(64),SERVER_ID,DemeID(uuid,MIXER_ID))

braider = Braider(config,deme,server)

for i in 1:3
    account = "account$i"
    member = Signer(deme,account * "/member")
    push!(braider.voters,member.id)
end

# ### Users do:

@sync for i in 1:3
    @async begin
        account = "account$i"
        member = Signer(deme,account * "/member")
        voter = Signer(deme,account * "/voters/$(member.id)")
        braid!(config,deme,voter,member)
    end
end

### After that gatekeeper gets ballot

contractbraid = take!(braider)
@show consentbraid = Consensus(contractbraid,deme.notary)


text:
################################################################################
#
#  Promote rules
#
################################################################################

AbstractAlgebra.promote_rule(::Type{S}, ::Type{fmpz}) where S <: NumFieldElem = S

AbstractAlgebra.promote_rule(::Type{fmpz}, ::Type{S}) where S <: NumFieldElem = S

AbstractAlgebra.promote_rule(::Type{S}, ::Type{fmpq}) where S <: NumFieldElem = S

AbstractAlgebra.promote_rule(::Type{fmpq}, ::Type{S}) where S <: NumFieldElem = S

AbstractAlgebra.promote_rule(::Type{T}, ::Type{S}) where {S <: NumFieldElem, T <: Integer} = S

AbstractAlgebra.promote_rule(::Type{S}, ::Type{T}) where {S <: NumFieldElem, T <: Integer} = S

AbstractAlgebra.promote_rule(::Type{NfRelElem{nf_elem}}, ::Type{nf_elem}) = NfRelElem{nf_elem}

AbstractAlgebra.promote_rule(::Type{NfRelNSElem{nf_elem}}, ::Type{nf_elem}) = NfRelNSElem{nf_elem}

AbstractAlgebra.promote_rule(::Type{NfRelElem{NfAbsNSElem}}, ::Type{NfAbsNSElem}) = NfRelElem{NfAbsNSElem}

AbstractAlgebra.promote_rule(::Type{NfRelNSElem{NfAbsNSElem}}, ::Type{NfAbsNSElem}) = NfRelNSElem{NfAbsNSElem}

AbstractAlgebra.promote_rule(::Type{nf_elem}, ::Type{NfRelElem{nf_elem}}) = NfRelElem{nf_elem}

AbstractAlgebra.promote_rule(::Type{nf_elem}, ::Type{NfRelNSElem{nf_elem}}) = NfRelNSElem{nf_elem}

AbstractAlgebra.promote_rule(::Type{NfAbsNSElem}, ::Type{NfRelElem{NfAbsNSElem}}) = NfRelElem{NfAbsNSElem}

AbstractAlgebra.promote_rule(::Type{NfAbsNSElem}, ::Type{NfRelNSElem{NfAbsNSElem}}) = NfRelNSElem{NfAbsNSElem}

function AbstractAlgebra.promote_rule(::Type{NfRelNSElem{T}}, ::Type{NfAbsNSElem}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, NfAbsNSElem)
    return NfRelNSElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{NfAbsNSElem}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, NfAbsNSElem)
    return NfRelElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfAbsNSElem}, ::Type{NfRelNSElem{T}}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, NfAbsNSElem)
    return NfRelNSElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfAbsNSElem}, ::Type{NfRelElem{T}}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, NfAbsNSElem)
    return NfRelElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelNSElem{T}}, ::Type{nf_elem}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, nf_elem)
    return NfRelNSElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{nf_elem}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, nf_elem)
    return NfRelElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{nf_elem}, ::Type{NfRelNSElem{T}}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, nf_elem)
    return NfRelNSElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{nf_elem}, ::Type{NfRelElem{T}}) where T <: NumFieldElem
  if T === AbstractAlgebra.promote_rule(T, nf_elem)
    return NfRelElem{T}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelElem{S}}, ::Type{NfRelElem{S}}) where S <: NumFieldElem
  return NfRelElem{S}
end

function AbstractAlgebra.promote_rule(::Type{NfRelNSElem{S}}, ::Type{NfRelNSElem{S}}) where S <: NumFieldElem
  return NfRelNSElem{S}
end

function AbstractAlgebra.promote_rule(::Type{NfRelElem{S}}, ::Type{NfRelElem{T}}) where {S <: NumFieldElem, T <: NumFieldElem}
  if S === T
    return NfRelElem{S}
  end
  U = AbstractAlgebra.promote_rule(S, T)
  if U === T
    return NfRelElem{T}
  elseif U === S
    return NfRelElem{S}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelNSElem{S}}, ::Type{NfRelNSElem{T}}) where {S <: NumFieldElem, T <:NumFieldElem}
  if S === T
    return NfRelNSElem{S}
  end
  U = AbstractAlgebra.promote_rule(S, T)
  if U === T
    return NfRelNSElem{T}
  elseif U === S
    return NfRelNSElem{S}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelNSElem{S}}, ::Type{NfRelElem{T}}) where {T <: NumFieldElem, S <: NumFieldElem}
  if T === NfRelNSElem{S}
    return NfRelElem{T}
  elseif S === NfRelElem{T}
    return NfRelNSElem{S}
  elseif AbstractAlgebra.promote_rule(S, NfRelElem{T}) === NfRelElem{T}
    return NfRelElem{T}
  elseif AbstractAlgebra.promote_rule(T, NfRelNSElem{S}) == NfRelNSElem{S}
    return NfRelNSElem{S}
  else
    return Union{}
  end
end

function AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{NfRelNSElem{S}}) where {T <: NumFieldElem, S <: NumFieldElem}
  if T === NfRelNSElem{S}
    return NfRelElem{T}
  elseif S === NfRelElem{T}
    return NfRelNSElem{S}
  elseif AbstractAlgebra.promote_rule(S, NfRelElem{T}) === NfRelElem{T}
    return NfRelElem{T}
  elseif AbstractAlgebra.promote_rule(T, NfRelNSElem{S}) == NfRelNSElem{S}
    return NfRelNSElem{S}
  else
    return Union{}
  end
end

################################################################################
#
#  Coercion generic
#
################################################################################

(K::NfRel{T})(x::NfRelElem{T}) where {T <: NumFieldElem} = K(x.data)
(K::NfRelNS{T})(x::NfRelNSElem{T}) where {T <: NumFieldElem} = K(x.data)

function (K::NfRel{S})(a::T) where {S <: NumFieldElem, T <: NumFieldElem}
  if S === T && parent(a) == base_field(K)
    return K(parent(K.pol)(a))
  end
  ET = elem_type(K)
  if AbstractAlgebra.promote_rule(T, ET) !== ET
    return force_coerce_throwing(K, a)
  end
  el = base_field(K)(a)
  return K(el)
end

function (K::NfRelNS{S})(a::T) where {S <: NumFieldElem, T <: NumFieldElem}
  if S === T
    return K(parent(K.pol[1])(a))
  end
  ET = elem_type(K)
  if AbstractAlgebra.promote_rule(T, ET) !== ET
    error("Element not coercible!")
  end
  el = base_field(K)(a)
  return K(el)
end


text:
# # Grid operations

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
## Basic differential operations
There are a variety of (purely) grid-based operators that are useful for carrying
out calculations in immersed layer problems. We will demonstrate a few of them
here.
We will start by generating the cache, just as we did in [Immersed layer caches](@ref)
=#

using ImmersedLayers
using CartesianGrids
#!jl using Plots

#=
### Set up a grid and cache
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)

#=
We still generate a cache for these operations, but
now, we only supply the grid. There are no immersed surfaces
for this demonstration.
=#
cache = SurfaceScalarCache(g,scaling=GridScaling)

#=
To demonstrate, let's generate a Gaussian
=#
p = zeros_grid(cache)
xg, yg = x_grid(cache), y_grid(cache)
p .= exp.(-(xg∘xg)-(yg∘yg))

#=
Now, let's generate the gradient of these data
=#
v = zeros_gridgrad(cache)
grad!(v,p,cache)
#!jl plot(v,cache)

#=
We can then compute the derivative of this data
=#
divv = zeros_grid(cache)
divergence!(divv,v,cache)
#!jl plot(divv,cache)

#=
## Convective derivatives
And finally, let's compute convective derivatives. First, we will compute

$$\mathbf{v}\cdot\nabla p$$

For this operation, we create a special additional cache using [`ConvectiveDerivativeCache`](@ref).
This extra cache holds additional memory for making the calculation of the convective derivative faster
if we compute it often.
=#
cdcache = ConvectiveDerivativeCache(cache)
vdp = zeros_grid(cache)
convective_derivative!(vdp,v,p,cache,cdcache) #hide
@time convective_derivative!(vdp,v,p,cache,cdcache)
nothing #hide

# Plot it
#!jl plot(vdp,cache)

#=
Now, let's compute

$$\mathbf{v}\cdot\nabla\mathbf{v}$$

For this, we create a cache for `VectorGridData`, and a new instance of
`ConvectiveDerivativeCache` to go along with it.
=#
vcache = SurfaceVectorCache(g,scaling=GridScaling)
vdv = zeros_grid(vcache)
cdvcache = ConvectiveDerivativeCache(vcache)
convective_derivative!(vdv,v,vcache,cdvcache) #hide
@time convective_derivative!(vdv,v,vcache,cdvcache)
nothing #hide

# Plot it
#!jl plot(vdv,vcache)

#md # ## Grid operator functions
#md # ```@docs
#md # divergence!
#md # grad!
#md # curl!
#md # convective_derivative!
#md # convective_derivative
#md # ConvectiveDerivativeCache
#md # laplacian!
#md # inverse_laplacian!
#md # ```


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 9720d140-8bf9-11eb-3391-7da0608d03d5
using DelimitedFiles

# ╔═╡ 5e2a99c0-8c4e-11eb-2261-3596c3b51071
using CSV,DataFrames

# ╔═╡ f7b7af70-8b35-11eb-0f84-0f9081e07bd0
begin
		rxn1 = [-1 0 -1 1 0 0 0 0 0 -1 1 1 0 0 0 0 0 0]'
		rxn2 = [0 0 0 -1 1 1 0 0 0 0 0 0 0 0 0 0 0 0]'
		rxn3 = [0 0 0 0 0 -1 1 1 0 0 0 0 -1 0 0 0 0 0]'
		rxn4 = [0 -1 1 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0]'
		rxn5f = [0 0 2 0 0 -2 0 0 0 0 0 0 4 -4 -3 -3 2 3]'
		rxn5r = rxn5f.*-1
		trans1 = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]'
		trans2 = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]'
		trans3 = [0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0]'
		trans4 = [0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0]'
		trans5 = [0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0]'
		trans6 = [0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0]'
		trans7 = [0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0]'
		trans8 = [0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0]'
		trans9 = [0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0]'
		trans10 = trans9.*-1
		trans11 = [0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]'
		trans12 = trans11.*-1
		trans13 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0]'
		trans14 = trans13.*-1
		trans15 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0]'
		trans16 = trans15.*-1
		trans17 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]'
		trans18 = trans17.*-1
		trans19 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]'
		trans20 = trans19.*-1

	s = hcat(rxn1,rxn2,rxn3,rxn4,rxn5f,rxn5r,trans1,trans2,trans3,trans4,trans5,trans6,trans7,trans8,trans9,trans10,trans11,trans12,trans13,trans14,trans15,trans16,trans17,trans18,trans19,trans20)
	print(s)
end

# ╔═╡ 84d03cd0-8c1a-11eb-0038-7717a9b98516
begin
	abc = open("notes.txt","w")
	writedlm(abc,s)
	close(abc)
	# the file f is automatically closed after this block finishes
end

# ╔═╡ ea9aa420-8b5b-11eb-399a-11706b03ddb8
# Checking elemental balances
s_rxn = s[:,1:6]

# ╔═╡ 04f6d650-8c50-11eb-2209-1b64eefbb63b
# reading the csv file
	d = CSV.read("C:\\Users\\bidst\\Google Drive\\Dropbox\\Documents\\School\\PhD Classes\\Chem_Eng 7770\\HW\\PS 3\\element_matrix_construction.csv",DataFrame)

# ╔═╡ 8ad9457e-8cbb-11eb-3748-bfdfeab03d9a
d_rxn = convert(Matrix, d[1:end-1,2:end])

# ╔═╡ 01f07c10-8cbc-11eb-30d7-4d4a4f9b5f38
elementcheckmatrix = d_rxn*s_rxn

# ╔═╡ Cell order:
# ╠═f7b7af70-8b35-11eb-0f84-0f9081e07bd0
# ╠═9720d140-8bf9-11eb-3391-7da0608d03d5
# ╠═84d03cd0-8c1a-11eb-0038-7717a9b98516
# ╠═ea9aa420-8b5b-11eb-399a-11706b03ddb8
# ╠═5e2a99c0-8c4e-11eb-2261-3596c3b51071
# ╠═04f6d650-8c50-11eb-2209-1b64eefbb63b
# ╠═8ad9457e-8cbb-11eb-3748-bfdfeab03d9a
# ╠═01f07c10-8cbc-11eb-30d7-4d4a4f9b5f38


text:
using GridArrays , LinearAlgebra, Test

@testset "gauss" begin
    x1,w1 = GridArrays.gausschebyshev(10)
    x2,w2 = GridArrays.gausschebyshev(BigFloat,10)
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    x1,w1 = GridArrays.gausschebyshevu(10)
    x2,w2 = GridArrays.gausschebyshevu(BigFloat,10)
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    x1,w1 = GridArrays.gausslegendre(10)
    x2,w2 = GridArrays.gausslegendre(BigFloat,10)
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    a = rand()
    x1,w1 = GridArrays.gausslaguerre(10, a)
    x2,w2 = GridArrays.gausslaguerre(10, big(a))
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    x1,w1 = GridArrays.gausshermite(10)
    x2,w2 = GridArrays.gausshermite(10)
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    a= rand();b=rand()
    x1,w1 = GridArrays.gaussjacobi(10,a,b)
    x2,w2 = GridArrays.gaussjacobi(10,big(a),big(b))
    @test norm(x1-x2)+norm(w1-w2) < 1e-12

    a= -.5
    x1,w1 = GridArrays.gaussjacobi(10,a,a)
    x2,w2 = GridArrays.gaussjacobi(10,big(a),big(a))
    x,w = GridArrays.gausschebyshev(10)
    @test w≈w1
    @test x≈x1
    x,w = GridArrays.gausschebyshev(BigFloat,10)
    @test x≈x2
    @test w≈w2

    a= .5
    x1,w1 = GridArrays.gaussjacobi(10,a,a)
    x2,w2 = GridArrays.gaussjacobi(10,big(a),big(a))
    x,w = GridArrays.gausschebyshevu(10)
    @test w≈w1
    @test x≈x1
    x,w = GridArrays.gausschebyshevu(BigFloat,10)
    @test x≈x2
    @test w≈w2
end


text:
@testset "Metropolis Hastings general" begin

    @testset "1D" begin
        Random.seed!(123456)
        prop(mu) = MvNormal(mu, 3)

        target(x) = log(pdf(MvNormal([2], 1), x))

        samps, acc = metropolis_hastings(target, prop, [2], 2000, 200)

        @test mean(samps) ≈ 2 atol = 0.1
        @test std(samps) ≈ 1 atol = 0.1
    end

    @testset "3D" begin
        Random.seed!(123456)
        prop(mu) = MvNormal(mu, [3 0 0; 0 3 0; 0 0 3])

        target(x) = log(pdf(MvNormal([-2, 2, 0], [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1]), x))

        samps, acc = metropolis_hastings(target, prop, [-2, 2, 0], 2000, 200)

        μ = mean(samps, dims = 1)
        σ = std(samps, dims = 1)
        corrs = cor(samps)
        @test vec(μ) ≈ [-2; 2; 0] atol = 0.2
        @test vec(σ) ≈ [1; 1; 1] atol = 0.1
        @test corrs ≈ [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1] atol = 0.2
    end

    @testset "No logged target" begin
        Random.seed!(123456)
        prop(mu) = MvNormal(mu, [3 0 0; 0 3 0; 0 0 3])

        target(x) = pdf(MvNormal([-2, 2, 0], [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1]), x)

        samps, acc = metropolis_hastings(target, prop, [-2, 2, 0], 2000, 200, islogged = false)

        μ = mean(samps, dims = 1)
        σ = std(samps, dims = 1)
        corrs = cor(samps)
        @test vec(μ) ≈ [-2; 2; 0] atol = 0.2
        @test vec(σ) ≈ [1; 1; 1] atol = 0.1
        @test corrs ≈ [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1] atol = 0.2
    end
end

@testset "Metropolis Hastings simple" begin

    @testset "1D" begin
        Random.seed!(123456)
        proprnd(mu) = rand(MvNormal(mu, 3))

        target(x) = log(pdf(MvNormal([2], 1), x))

        samps, acc = metropolis_hastings_simple(target, proprnd, [2], 2000, 200)

        @test mean(samps) ≈ 2 atol = 0.1
        @test std(samps) ≈ 1 atol = 0.1

    end

    @testset "3D" begin
        Random.seed!(123456)
        proprnd(mu) = rand(MvNormal(mu, [3 0 0; 0 3 0; 0 0 3]))

        target(x) = log(pdf(MvNormal([-2, 2, 0], [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1]), x))

        samps, acc = metropolis_hastings_simple(target, proprnd, [-2,2, 0], 2000, 200)

        μ = mean(samps, dims = 1)
        σ = std(samps, dims = 1)
        corrs = cor(samps)
        @test vec(μ) ≈ [-2; 2; 0] atol = 0.2
        @test vec(σ) ≈ [1; 1; 1] atol = 0.1
        @test corrs ≈ [1 0.5 0.5; 0.5 1 0.5; 0.5 0.5 1] atol = 0.2

    end
    Random.seed!()
end


text:
# # Эксперимент №3 Зашумленные

# ## Набор параметров
α       = 0.004;        # Параметр регуляризации
w       = 0.0001;       # Эмпирический параметр регуляризации
S       = 3000;        # Количество итераций
β       = 0.001;        # Шаг минимизации
#
x_tp    = 0.05;         # Стартовое местоположение фронта
T_end   = 0.47;         # Регулируем конечное местоположение фронта
ε       = 10^(-1.5);    # Крутизна фронта
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;         # Число интервалов по ``T``
δ       = 0.00008;      # Модуль величины помех

using Dates;
timestamp           = String(Dates.format(now(), "yyyy-mm-dd_HH-MM"));
wdir                = String("Expirement3___" * timestamp);
prefix              = String("Exp3_" * timestamp * "_");
n_direct_draft      = prefix * "direct_draft.png";
n_direct_gif        = prefix * "direct.gif";
n_noised_data       = prefix * "noised_data.png";
n_noised_velocity   = prefix * "noised_velocity.png";
n_initial_guess     = prefix * "initial_guess.png";
n_Js_min_data_guess = prefix * "Js_withguess.jld";
n_Qs_min_data_guess = prefix * "Qs_withguess.jld";
n_Js_min_data       = prefix * "Js_noguess.jld";
n_Qs_min_data       = prefix * "Qs_noguess.jld";
n_guess_draft       = prefix * "minimization_draft_withguess.png"
n_noguess_draft     = prefix * "minimization_draft_noguess.png"
n_minim_gif_guess   = prefix * "minimization_withguess.gif";
n_minim_gif_noguess = prefix * "minimization_noguess.gif";
n_minim_avi_guess   = prefix * "minimization_withguess.avi";
n_minim_avi_noguess = prefix * "minimization_noguess.avi";
n_minim_mp4_guess   = prefix * "minimization_withguess.mp4";
n_minim_mp4_noguess = prefix * "minimization_noguess.mp4";

using Markdown
info = """Вы зпустили эксперимент №3. Сценарий скрипта:
  - Создадим каталог `$(wdir)`.
  - Перейдем в него и все результаты поместим там.
  - Все файлы будем сохранять с префиксом `$(prefix)`
  - Сгенерируем экспериментальные данные на сетке `$(Nx+1) × $(Mt+1)`
    * Сохраним эскиз прямого решения `$(n_direct_draft)`
    * Сохраним анимацию решения прямой задачи, с помощью которой генерировали
        экспериментальные даные `$(n_direct_gif)`.
  - Зашумим экспериментальные данные f₁, f₂ равномерным распределение с модулем `$(δ)`
    * Сохраним график зашумленных экспериментальных данных `$(n_noised_data)`
    * Сохраним график скорости фронта из зашумленных данных `$(n_noised_velocity)`
  - Найдем начальное приближение с параметром регуляризации `$(α)`
    * Сохраним начальное приближение `$(n_initial_guess)`
  - Начнем процесс минимизации используя найденное приближение
    * Сериализуем данные и сохраним `$(n_Js_min_data_guess)`, `$(n_Qs_min_data_guess)`
    * Сохраним Черновик процесса минимизации `$(n_guess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_guess)`
  - Начнем процесс минимизации с нулевого начального приближения
    * Сериализуем данные и сохраним `$(n_Js_min_data)`, `$(n_Qs_min_data)`
    * Сохраним Черновик процесса минимизации `$(n_noguess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_noguess)`
"""
display(Markdown.parse(info))

@info "Переходим в $(wdir)"
run(`mkdir $(wdir)`)
run(`cd $(wdir)`)

# -----------------------------------------------------------------------------
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity_map;
using Serialization;
using Plots; pyplot();
using Dierckx;
# -----------------------------------------------------------------------------

# ## Решение на точных данных

# ### Решение прямой задачи для генерирования экспериментальной информации ----------------------------------------------------------------------------- a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
# -----------------------------------------------------------------------------
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения")
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# В следующей строке, мы записываем все параметры в latex формате,
# не забываем экранировать все спецсимволы.
# Дальше, мы отобразим все эти параметры на графике.
using Printf;
params = "\$w = $(w), \\varepsilon = $(@sprintf("%.4f", ε)), "*
"\\beta = $(β), f_1 \\in $(@sprintf("[%.2f, %.2f]", extrema(f1_data)...)), "*
"\\alpha = $(α)\$"
nothing #hide
# -----------------------------------------------------------------------------

# ## Слабо зашумленные данные
# -----------------------------------------------------------------------------
using Random
using NonLinearReactionAdvectionDiffusionWithFrontData: front_velocity;
rng = MersenneTwister(1234)

# Создадим мелкий шум.
f1_data_noised = f1_data + randn(rng, length(f1_data)) * δ;
f2_data_noised = f2_data + randn(rng, length(f2_data)) * δ;
h = Xₙ[2] - Xₙ[1];
δₙ = sqrt(sum( [ (f1_data[n] - f1_data_noised[n])^2 * h for n in 1:N+1]));
@info δₙ

plot(f1_data_noised[1:5:end], Tₘ[1:5:end], ylims = (0, T_end * 1.1), label = "", title="Зашумленные \$f_1, δₙ = $(δₙ)\$", xlim = (0,1))
plot!(f1_data[1:5:end], Tₘ[1:5:end], label = "")

#plot!(Tₘ, f2_data_noised, label="Зашумленные \$f_2\$")
savefig(n_noised_data);

# Убедимся, что численное дифференцирование зашумленных функций --- это плохо.
v_f1 = front_velocity(f1_data_noised, Tₘ, M);
plot(Tₘ, v_f1[5:end-5], label="Скорость зашумленного фронта \$\\dfrac{df_1}{dt}\$")
savefig(n_noised_velocity);

# -----------------------------------------------------------------------------
# ------------------ Старт с нулевого приближения -----------------------------
# -----------------------------------------------------------------------------
q₀ = zero(Xₙ)
# -----------------------------------------------------------------------------
@time qs, noised_noguess_Js, noised_noguess_Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize(n_Js_min_data, noised_noguess_Js);
serialize(n_Qs_min_data, noised_noguess_Qs);
# -----------------------------------------------------------------------------

get_k() = f1_data[end] - f1_data[1];
l2diff(V::Vector) = sqrt( sum( [ (qₙ[n] - V[n])^2 * h for n in 1:N+1]) )

plot(xlabel = "k", xlims = (0,1), ylims = (0, 0.7), ylabel = "|| q^inv - q||_l2");
plot!(ks[1, :], map(l2diff, [s3, m3, l3]), label = "x_tp = 0.35")
plot!(ks[2,:], map(l2diff, [s2, m2, l2]), label = "x_tp = 0.25")
plot!(ks[3,:], map(l2diff, [s1, m1, l1]), label = "x_tp = 0.05")
savefig("9p.eps")
savefig("9p.png")

r_9p = [
          s3 m3 l3
          s2 m2 l2
          s1 m1 l1
         ];
r_ks = [
          ks3 km3 kl3
          ks2 km2 kl2
          ks1 km1 kl1
         ]

s3 = p9[1:N, 1]
m3 = p9[1:N, 2]
l3 = p9[1:N, 3]
s2 = p9[N+1:N+N+1, 1]
m2 = p9[N+1:N+N+1, 2]
l2 = p9[N+1:N+N+1, 3]
s1 = p9[2N+1:end, 1]
m1 = p9[2N+1:end, 2]
l1 = p9[2N+1:end, 3]


text:
function all_values(abaco, ts, ne, snap)
    ovals = origins_vals(abaco, ne, ts)
    merge(snap.vals, ovals)
end

function deep_count(abaco, ne, domains)
    #@debug "deep_count: $ne, domains: $domains"
    count = 0
    # TODO: filter on tag value
    if length(domains) === 1
        return length(filter(el->el.tag == domains[1], abaco.origins[ne]))
    else
        for node in abaco.origins[ne]
            count += deep_count(abaco, node.ne, domains[2:end])
        end
    end
    count
end

function origins_vals(abaco::Context, ne, ts)
    vals = Dict()
    index = abaco.interval == -1 ? 1 : mvindex(ts, abaco.interval, abaco.ages)

    if haskey(abaco.origins, ne)
        for origin_elem in abaco.origins[ne]
            ropts = span(ts, abaco.interval)
            snap = origin_elem.snap[index]
            # @debug "$ropts --> origin_elem [$(origin_elem.tag).$(origin_elem.ne)]: $snap"

            for (var, val) in all_values(abaco, ts, origin_elem.ne, snap)
                newvar = "$(origin_elem.tag).$var"
                #@debug "[$ne] getting var $newvar"
                if !haskey(vals, newvar)
                    nodes = split(newvar, ".")[1:end-1]
                    value = LValue(deep_count(abaco, ne, nodes), [])
                    vals[newvar] = value
                end
                #@debug "$ne var $newvar: appending $(val.value) (was $(vals[newvar].value))"
                # if interval == -1 considers all values collected at different times
                if (abaco.interval == -1 || ropts === snap.ts)
                    append!(vals[newvar].value, val.value)
                end
            end
        end
    end
    vals
end

function propagate(abaco::Context, snap, ne, formula_name)
    ts = snap.ts
    (etype, target) = abaco.target[ne]
    var = "$etype.$formula_name"
    #@debug "[$ne]: ts:$ts - propagating $var to node [$target]"
    trigger_formulas(abaco, getsnap(abaco, ts, target), target, var)
end

function trigger_formulas(abaco, snap, ne, var::String)
    poll_formulas(abaco, snap, ne, var)

    if haskey(abaco.target, ne)
        propagate(abaco, snap, ne, var)
    end
end

function trigger_formulas(abaco, snap, ne, vars)
    poll_formulas(abaco, snap, ne, vars)

    for var in vars
        # propagate input variables to target
        if haskey(abaco.target, ne)
            propagate(abaco, snap, ne, var)
        end
    end
end

update_inputs(abaco, snap, ne, formula_name, result::Real) = snap_add(snap, ne, formula_name, result)

function update_inputs(abaco, snap, ne, formula_name, result::PValue)
    snap_add(snap, ne, formula_name, result)
    trigger_formulas(abaco, snap, ne, formula_name)
end

#
#     poll_formulas(abaco, snap, ne, vars)
#
# Check the completion status of all formulas that depends on `vars`.
#
function poll_formulas(abaco, snap, ne, vars)
    if snap === nothing
        return
    end
    tag = etype(abaco, ne)
    cfg = snapsetting(abaco, tag)
    formulas = dependents(abaco, tag, vars)

    @debug "[$ne] formulas that depends on [$vars]: $formulas"
    for formula_name in formulas
        fstate = get(snap.outputs, formula_name, FormulaState(false, formula_name))

        # pull up the origins variables
        ovals = origins_vals(abaco, ne, snap.ts)
        allvals = merge(snap.vals, ovals)
        result = poll(cfg, fstate, allvals)

        #@debug "[$ne] allvals: $allvals"
        @debug "[$ne] poll($formula_name) = $result"
        if result !== nothing

            #snap_add(snap, ne, formula_name, result)
            update_inputs(abaco, snap, ne, formula_name, result)

            # propagate to target
            if haskey(abaco.target, ne)
                #@info "propagate $ne --> $formula_name"
                propagate(abaco, snap, ne, formula_name)
            end

            if cfg.oncomplete !== nothing
                if cfg.handle === nothing
                    cfg.oncomplete(snap.ts,
                                   ne,
                                   formula_name,
                                   result,
                                   allvals)
                else
                    cfg.oncomplete(cfg.handle,
                                   snap.ts,
                                   ne,
                                   formula_name,
                                   result,
                                   allvals)
                end
            end
        end
    end
end

function trigger_formulas(abaco, ::Nothing, ne, var, value)
end

#
#     poll(abaco, formula_state::FormulaState, vals::Dict)
#
# Returns the formula value if the formula is computable.
#
# Returns `NaN` when some input variables are missing because
# the formula cannot run to completion.
#
function poll(setting::SnapsSetting, formula_state::FormulaState, vals::Dict)
    formula = setting.formula[formula_state.f]

    poll_simple(setting, formula, formula_state, vals)
end

function poll_simple(setting::SnapsSetting, formula::Formula, formula_state::FormulaState, vals::Dict)
    # Step 1: decide if the formula may be evaluated: all variables are collected
    for v in formula.inputs
        #if !(haskey(vals, v) && isready(vals[v]))
        if !(haskey(vals, v))
            return nothing
        end
    end

    # Step 2: evaluate the formula
    if !formula_state.done
        if setting.emitone
            formula_state.done = true
        end
        result = eval(formula, vals)
        return result
    end
    return nothing
end


text:
# This file is a part of BitOperations.jl, licensed under the MIT License (MIT).

using BitOperations
using Test

@testset "bit operations" begin
    @testset "bsizeof" begin
        @test bsizeof(one(Int64)) == 64
        @test bsizeof(Int16) == 16
    end

    @testset "multi-bit bmask" begin
        @test bmask(UInt64, 0:0) == 0x1
        @test bmask(UInt64, 0:3) == 0xf
        @test bmask(UInt64, 0:63) == 0xffffffffffffffff
        @test bmask(UInt64, 60:63) == 0xf000000000000000
        @test bmask(UInt64, 60:63) == 0xf000000000000000
        @test bmask(UInt64, 63:63) == 0x8000000000000000
    end

    @testset "single-bit lsbmask and msbmask" begin
        @test lsbmask(UInt64) == 0x1
        @test msbmask(UInt64) == 0x8000000000000000

        @test lsbmask(Int64) == 1
        @test msbmask(Int64) == signed(0x8000000000000000)
    end

    @testset "multi-bit lsbmask and msbmask" begin
        @test lsbmask(UInt64, 0) == 0x0
        @test lsbmask(UInt64, 64) == 0xffffffffffffffff
        @test lsbmask(UInt64, 6) == 0x000000000000003f

        @test msbmask(UInt64, 0) == 0x0
        @test msbmask(UInt64, 64) == 0xffffffffffffffff
        @test msbmask(UInt64, 6) == 0xfc00000000000000
    end

    @testset "single-bit bget and bset" begin
        @test bget(8, 0) == false
        @test bget(8, 3) == true
        @test bget(0x8000000000000000, 63) == true
        @test bget(0x8000000000000000, 64) == false

        @test bset(0x04, 7, true) == 0x84
        @test bset(0xff, 7, false) == 0x7f
    end

    @testset "multi-bit bget and bset" begin
        @test bget(0x1234, 4:11) == 0x23
        @test bset(0x15a4, 4:11, 0x23) == 0x1234
    end

    @testset "single-bit bflip" begin
        @test bflip(0b1001000110100, 5) == 0b1001000010100
        @test bflip(0b1001000110100, 9) == 0b1000000110100
    end

    @testset "multi-bit bflip" begin
        @test bflip(0xa5a5, 4:11) == 0xaa55
    end

    @testset "single-bit lsbget and msbget" begin
        @test lsbget(0) == false
        @test lsbget(1) == true

        @test msbget(0x7fffffffffffffff) == false
        @test msbget(0x8000000000000000) == true
    end

    @testset "multi-bit lsbget and msbget" begin
        x = 0xabffffffffffffd5

        @test lsbget(x, 0) == 0
        @test lsbget(x, 6) == 0x15
        @test lsbget(x, 64) == x

        @test msbget(x, 0) == 0
        @test msbget(x, 6) == 0x2a
        @test msbget(x, 64) == x
    end
end


text:
# include("../preamble.jl")
# include((@__DIR__)*"/../../preamble.jl")

T_reflecting = [-1.1 1.1; 1.0 -1.0]
c_reflecting = [1.0; -1.0]
b_reflecting = 10.0

P_upr_reflecting = [0.0 1.0]
P_lwr_reflecting = [1.0 0.0] # both reflecting boundaries

model = BoundedFluidQueue(T_reflecting,c_reflecting,P_lwr_reflecting,P_upr_reflecting,b_reflecting)


text:
#=
    borrow from https://github.com/Roger-luo/BatchedRoutines.jl
=#

import LinearAlgebra

#batched cpu gemm by BatchedRoutines.jl
for (gemm, elty) in
    ((:dgemm_,:Float64),
     (:sgemm_,:Float32),)
    @eval begin
        function batched_gemm!(transA::AbstractChar,
                               transB::AbstractChar,
                               alpha::($elty),
                               A::AbstractArray{$elty, 3},
                               B::AbstractArray{$elty, 3},
                               beta::($elty),
                               C::AbstractArray{$elty, 3})
            @assert !Base.has_offset_axes(A, B, C)
            @assert size(A, 3) == size(B, 3) == size(C, 3) "batch size mismatch"
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            LinearAlgebra.BLAS.chkstride1(A)
            LinearAlgebra.BLAS.chkstride1(B)
            LinearAlgebra.BLAS.chkstride1(C)

            ptrA = Base.unsafe_convert(Ptr{$elty}, A)
            ptrB = Base.unsafe_convert(Ptr{$elty}, B)
            ptrC = Base.unsafe_convert(Ptr{$elty}, C)

            for k in 1:size(A, 3)
                ccall((LinearAlgebra.BLAS.@blasfunc($gemm), LinearAlgebra.BLAS.libblas), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{LinearAlgebra.BLAS.BlasInt}, Ref{LinearAlgebra.BLAS.BlasInt},
                     Ref{LinearAlgebra.BLAS.BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{LinearAlgebra.BLAS.BlasInt},
                     Ptr{$elty}, Ref{LinearAlgebra.BLAS.BlasInt}, Ref{$elty}, Ptr{$elty},
                     Ref{LinearAlgebra.BLAS.BlasInt}),
                     transA, transB, m, n,
                     ka, alpha, ptrA, max(1,stride(A,2)),
                     ptrB, max(1,stride(B,2)), beta, ptrC,
                     max(1,stride(C,2)))

                ptrA += size(A, 1) * size(A, 2) * sizeof($elty)
                ptrB += size(B, 1) * size(B, 2) * sizeof($elty)
                ptrC += size(C, 1) * size(C, 2) * sizeof($elty)
            end

            C
        end
    end
end


text:
function mergesort(v)
    lv = length(v)
    if lv == 1
        return v
    else
        return fusion(fusionsort(@view v[1:lv÷2]), fusionsort(@view v[lv÷2+1:end]))
    end
end

function merge(v1, v2)
    if length(v1) == 0
        return v2
    elseif length(v2) == 0
        return v1
    elseif v2[1] < v1[1]
        return vcat(v2[1], merge(v1, @view v2[2:end]))
    else
        return vcat(v1[1], merge((@view v1[2:end]), v2))
    end
end

v = [144, 89, 4, 9, 95, 12, 86, 25]
println("unsorted array = ",v)

sorted_v = mergesort(v)

println("merge sorted array = ", sorted_v)


text:
"""
Main module of Modia.

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

"""
module Modia

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.8.2"
const Date = "2022-03-08"
const modelsPath = joinpath(Modia.path, "models")

print(" \n\nWelcome to ")
print("Mod")
printstyled("ia", bold=true, color=:red)
print(" - ")
printstyled("Dynamic ", color=:light_black)
print("Mod")
printstyled("eling and Simulation with Jul", color=:light_black)
printstyled("ia", bold=true, color=:red)

println()
println("Version $Version ($Date)")

# Defalut switch settings
logStatistics = false
logExecution = false
logCalculations = false

useNewCodeGeneration = false
experimentalTranslation = false

using Reexport

@reexport using Unitful                 # export Unitful symbols
@reexport using DifferentialEquations   # export DifferentialEquations symbols

export ModiaBase
export CVODE_BDF, IDA
export instantiateModel, @instantiateModel, assert, stringifyDefinition
export stripUnit

export simulate!, linearize!, get_result
export @usingModiaPlot, usePlotPackage, usePreviousPlotPackage, currentPlotPackage
export resultInfo, printResultInfo, rawSignal, getPlotSignal, defaultHeading
export signalNames, timeSignalName, hasOneTimeSignal, hasSignal

export SimulationModel, measurementToString, get_lastValue
export positive, negative, previous, edge, after, reinit, pre
export initial, terminal, isInitial, isTerminal, initLinearEquationsIteration!
export get_xNames
export registerExtraSimulateKeywordArguments
export get_extraSimulateKeywordArgumentsDict

export modelToJSON, JSONToModel, writeModel, readModel

import Sundials
const  CVODE_BDF = Sundials.CVODE_BDF
const  IDA = Sundials.IDA

using Base.Meta: isexpr
using OrderedCollections: OrderedDict

using ModiaBase
using ModiaBase.Symbolic
using ModiaBase.Simplify
using ModiaBase.BLTandPantelidesUtilities
using ModiaBase.BLTandPantelides
using ModiaBase.Differentiate

import ModiaResult
import ModiaResult: usePlotPackage, usePreviousPlotPackage, currentPlotPackage
import ModiaResult: resultInfo, printResultInfo, rawSignal, getPlotSignal, defaultHeading
import ModiaResult: signalNames, timeSignalName, hasOneTimeSignal, hasSignal

import StaticArrays   # Make StaticArrays available for the tests

using  Measurements
import MonteCarloMeasurements
using JSON
#using Profile
import Test
using  TimerOutputs
using  InteractiveUtils

global to = TimerOutput()

Unitful.unit(      v::MonteCarloMeasurements.AbstractParticles{T,N}) where {T,N} = unit(T)
Unitful.upreferred(v::MonteCarloMeasurements.AbstractParticles{T,N}) where {T,N} = uconvert(upreferred(unit(v)), v)

# append! as needed in EquationAndStateInfo.jl and in CodeGeneration.jl
appendVariable!(v1::Vector{FloatType}, s::FloatType) where {FloatType} = push!(v1,s)
appendVariable!(v1::Vector{FloatType}, v2)           where {FloatType} = append!(v1,v2)

"""
    stripUnit(v)

Convert the unit of `v` to the preferred units (default are the SI units),
and then strip the unit. For details see `upreferred` and `preferunits` in
[Unitful](https://painterqubits.github.io/Unitful.jl/stable/conversion/)

The function is defined as: `stripUnit(v) = ustrip.(upreferred.(v))`.
"""
stripUnit(v) = ustrip.(upreferred.(v))

"""
     str = unitAsString( unitOfQuantity::Unitful.FreeUnits )

Return a string representation of the unit of a quantity that can be used in a unit string macro
(see also Unitful [issue 412](https://github.com/PainterQubits/Unitful.jl/issues/412)).

# Example
```
v1 = 2.0u"m/s"
v1_unit = unitAsString( unit(v1) )   # = "m*s^-1"
v2_withoutUnit = 2.0
code = :( \$v2_withoutUnit@u_str(\$v1_unit) )  # = 2.0u"m*s^-1"
v2 = eval(code)
@show v1
@show v1_unit
@show v2
```

# Notes
Transforms unit to string representation that is parseable again
(see also Unitful [issue 412](https://github.com/PainterQubits/Unitful.jl/issues/412)).
This implementation is a hack and only works in common cases.
Implementation is performed in the following way:

1. Transform to string and display exponents on units not as Unicode superscripts (= default on macOS).
2. Replace " " by "*", since Unitful removes "*" when converting to string.
"""
unitAsString(unitOfQuantity::Unitful.FreeUnits) = replace(repr(unitOfQuantity,context = Pair(:fancy_exponent,false)), " " => "*")

"""
    quantityType = quantity(numberType, numberUnit::Unitful.FreeUnits)

Return the quantity type given the numberType and the numberUnit.

# Example
```julia
mutable struct Data{FloatType <: AbstractFloat}
    velocity::quantity(FloatType, u"m/s")
end

v = Data{Float64}(2.0u"mm/s")
@show v                         # v =
```
"""
quantity(numberType, numberUnit::Unitful.FreeUnits) = Quantity{numberType, dimension(numberUnit), typeof(numberUnit)}

quantityTypes(::Type{Unitful.Quantity{T,D,U}}) where {T,D,U} = (T,D,U)

include("EquationAndStateInfo.jl")
include("StateSelection.jl")

include("ModelCollections.jl")
include("EvaluateParameters.jl")
include("EventHandler.jl")
include("CodeGeneration.jl")
# include("GenerateGetDerivatives.jl")
include("Synchronous.jl")
include("SimulateAndPlot.jl")
include("ReverseDiffInterface.jl")
include("PathPlanning.jl")
include("JSONModel.jl")
using .JSONModel

# include("IncidencePlot.jl")
# using .IncidencePlot
# Base.Experimental.@optlevel 0

const drawIncidence = false

include("Symbolic.jl")
include("ModiaLang.jl")

end


text:
##############################################################################
##
## FixedEffect
##
##############################################################################

type FixedEffect{R <: Integer, W <: AbstractVector{Float64}, I <: AbstractVector{Float64}}
    refs::Vector{R}         # refs of the original PooledDataVector
    sqrtw::W                # weights
    scale::Vector{Float64}  # 1/(∑ sqrt(w) * interaction) within each group
    interaction::I          # the continuous interaction
    factorname::Symbol      # Name of factor variable
    interactionname::Symbol # Name of continuous variable in the original dataframe
    id::Symbol              # Name of new variable if save = true
end

##############################################################################
##
## Constructor
##
##############################################################################

# Constructor
function FixedEffect{R <: Integer}(
    refs::Vector{R}, l::Integer, sqrtw::AbstractVector{Float64},
    interaction::AbstractVector{Float64}, factorname::Symbol,
    interactionname::Symbol, id::Symbol)
    scale = zeros(Float64, l)
    @inbounds @simd for i in 1:length(refs)
         scale[refs[i]] += abs2(interaction[i] * sqrtw[i])
    end
    @inbounds @simd for i in 1:l
           scale[i] = scale[i] > 0 ? (1.0 / sqrt(scale[i])) : 0.
       end
    FixedEffect(refs, sqrtw, scale, interaction, factorname, interactionname, id)
end

# Constructors from dataframe + terms
function FixedEffect(df::AbstractDataFrame, terms::Terms, sqrtw::AbstractVector{Float64})
    out = FixedEffect[]
    for term in terms.terms
        result = FixedEffect(df, term, sqrtw)
        if isa(result, FixedEffect)
            push!(out, result)
        elseif isa(result, Vector{FixedEffect})
            append!(out, result)
        end
    end
    return out
end

# Constructors from dataframe + symbol
function FixedEffect(df::AbstractDataFrame, a::Symbol, sqrtw::AbstractVector{Float64})
    v = df[a]
    if isa(v, PooledDataVector)
        return FixedEffect(v.refs, length(v.pool), sqrtw, Ones{Float64}(length(v)), a, :none, a)
    else
        # x from x*id -> x + id + x&id
        return nothing
    end
end

# Constructors from dataframe + expression
function FixedEffect(df::AbstractDataFrame, a::Expr, sqrtw::AbstractVector{Float64})
    _check(a) || throw("Expression $a shouyld only contain & and variable names")
    factorvars, interactionvars = _split(df, allvars(a))
    if isempty(factorvars)
        # x1&x2 from (x1&x2)*id
        return nothing
    end
    z = group(df, factorvars)
    interaction = _multiply(df, interactionvars)
    factorname = _name(factorvars)
    interactionname = _name(interactionvars)
    id = _name(allvars(a))
    l = length(z.pool)
    return FixedEffect(z.refs, l, sqrtw, interaction, factorname, interactionname, id)
end

function _check(a::Expr)
    a.args[1] == :& && check(a.args[2]) && check(a.args[3])
end
check(a::Symbol) = true

function _name(s::Vector{Symbol})
    if isempty(s)
        out = :none
    else
        out = convert(Symbol, reduce((x1, x2) -> string(x1)*"x"*string(x2), s))
    end
    return out
end

function _split(df::AbstractDataFrame, ss::Vector{Symbol})
    catvars, contvars = Symbol[], Symbol[]
    for s in ss
        isa(df[s], PooledDataVector) ? push!(catvars, s) : push!(contvars, s)
    end
    return catvars, contvars
end

function _multiply(df, ss::Vector{Symbol})
    if isempty(ss)
        out = Ones(size(df, 1))
    else
        if isa(df[ss[1]], Vector{Float64})
            out = deepcopy(df[ss[1]])
        else
            out = convert(Vector{Float64}, df[ss[1]])
        end
        for i in 2:length(ss)
            broadcast!(*, out, out, df[ss[i]])
        end
    end
    return out
end


text:
export Payload

struct Payload <: Rocket
    mass::Float64
end

mass(payload::Payload) = payload.mass
thrust(r::Payload) = 0
update!(r::Payload, t::Number, Δt::Number) = nothing


text:
# update_Cov.jl

function getW(κ::T, D::Matrix{T}, corParams::CorParams) where T <: Real
    κCor = covMat(D, κ, corParams) # not the inteded use of the function, but this works.
    W = PDMat(κCor + I)
    Winv1 = W \ ones(T, size(W,1))
    sumWinv = sum(Winv1)
    return (W, sumWinv)
end

function getSufficients_W(y::Vector{T}, W::PDMat, sumWinv::T,
    mixcomp_type::MixComponentNormal=MixComponentNormal()) where T <: Real

    yhat = sum(W \ y) / sumWinv
    d = y .- yhat
    Winvd = W \ d
    s2 = d'Winvd

    return (yhat, s2)
end

function lmarg_W(W::PDMat, sumWinv::T, yhat::T, s2::T, σ2::T,
    priorMixcomponent::PriorMixcomponent_Normal) where T <: Real

    a1 = logdet(W)
    a2 = s2 / σ2
    a3 = log(sumWinv)
    a40 = σ2 / sumWinv + priorMixcomponent.μ.σ2
    a4 = log(a40)
    a5 = (yhat - priorMixcomponent.μ.μ)^2 / a40

    return -0.5*( a1 + a2 + a3 + a4 + a5 )
end

function ldens_proptoIG_Jacobian(lx::T, shape::T, scale::T) where T <: Real
    ## Jacobian included for log scale
    x = exp(lx)
    return -shape * lx - scale / x
end

function lprior_κCor(lκ::T, lcorParam::T, κHypers::SNR_Hyper_ScInvChiSq,
    corParamHypers::MaternHyper_ScInvChiSq) where T <: Real

    la = ldens_proptoIG_Jacobian(lκ, 0.5*κHypers.κ_ν, 0.5*κHypers.κ_ν*κHypers.κ0)
    lb = ldens_proptoIG_Jacobian(lcorParam, 0.5*corParamHypers.lenscale_ν,
        0.5*corParamHypers.lenscale_ν*corParamHypers.lenscale0)
    return la + lb
end

function update_Cov!(mixcomp::MixComponentNormal, y::Vector{T}, # assumes y is full
    priorMixcomponent::PriorMixcomponent_Normal,
    κHypers_type::SNR_Hyper_ScInvChiSq, corParamHypers_type::MaternHyper_ScInvChiSq) where T <: Real

    nζon = length(mixcomp.ζon_indx)

    lκ_old = log(mixcomp.κ)
    lcorParam_old = log(mixcomp.corParams.lenscale)

    if nζon > 0

        W_old, sumWinv_old = getW(mixcomp.κ, mixcomp.D[mixcomp.ζon_indx, mixcomp.ζon_indx], mixcomp.corParams)
        yhat_old, s2_old =  getSufficients_W(y[mixcomp.ζon_indx], W_old, sumWinv_old)

        llik_old = lmarg_W(W_old, sumWinv_old, yhat_old, s2_old, mixcomp.σ2, priorMixcomponent)
        lpri_old = lprior_κCor(lκ_old, lcorParam_old, mixcomp.κ_hypers, mixcomp.corHypers)

        cand = [ lκ_old, lcorParam_old ] + rand(mixcomp.rng, MvNormal(mixcomp.cSig.mat))
        κ_cand = exp(cand[1])
        corParams_cand = deepcopy(mixcomp.corParams)
        corParams_cand.lenscale = exp(cand[2])

        W_cand, sumWinv_cand = getW(κ_cand, mixcomp.D[mixcomp.ζon_indx, mixcomp.ζon_indx], corParams_cand)
        yhat_cand, s2_cand =  getSufficients_W(y[mixcomp.ζon_indx], W_cand, sumWinv_cand)

        llik_cand = lmarg_W(W_cand, sumWinv_cand, yhat_cand, s2_cand, mixcomp.σ2, priorMixcomponent)
        lpri_cand = lprior_κCor(cand[1], cand[2], mixcomp.κ_hypers, mixcomp.corHypers)

        lar = llik_cand + lpri_cand - llik_old - lpri_old

        if log(rand(mixcomp.rng)) < lar
            mixcomp.κ = deepcopy(κ_cand)
            mixcomp.corParams = deepcopy(corParams_cand)
            mixcomp.Cor = corrMat(mixcomp.D, mixcomp.corParams)
            mixcomp.accpt += 1

            lparams_out = deepcopy(cand)
            out = Dict(:sumWinv => sumWinv_cand, :yhat => yhat_cand, :s2 => s2_cand)
        else
            lparams_out = vcat(deepcopy(lκ_old), deepcopy(lcorParam_old))
            out = Dict(:sumWinv => sumWinv_old, :yhat => yhat_old, :s2 => s2_old)
        end

    else # if no observations assigned to the updated lag

        lpri_old = lprior_κCor(lκ_old, lcorParam_old, mixcomp.κ_hypers, mixcomp.corHypers)

        cand = [ lκ_old, lcorParam_old ] + rand(mixcomp.rng, MvNormal(mixcomp.cSig.mat))
        κ_cand = exp(cand[1])
        corParams_cand = deepcopy(mixcomp.corParams)
        corParams_cand.lenscale = exp(cand[2])

        lpri_cand = lprior_κCor(cand[1], cand[2], mixcomp.κ_hypers, mixcomp.corHypers)

        lar = lpri_cand - lpri_old

        if log(rand(mixcomp.rng)) < lar
            mixcomp.κ = deepcopy(κ_cand)
            mixcomp.corParams = deepcopy(corParams_cand)
            mixcomp.Cor = corrMat(mixcomp.D, mixcomp.corParams)
            mixcomp.accpt += 1

            lparams_out = deepcopy(cand)
            out = Dict()
        else
            lparams_out = vcat(deepcopy(lκ_old), deepcopy(lcorParam_old))
            out = Dict()
        end
    end

    if mixcomp.adapt
        mixcomp.adapt_iter += 1
        mixcomp.runningsum_Met += lparams_out
        runningmean = mixcomp.runningsum_Met / float(mixcomp.adapt_iter)
        runningdev = ( lparams_out - runningmean )
        mixcomp.runningSS_Met = runningdev * runningdev' # the mean is changing, but this approx. is fine.
    end

    return out
end


text:
using PyCall
pygui(:qt)
using PyPlot
matplotlib[:rcdefaults]()

using InterpNd

x = logspace(-1, 0, 30)
y = sin(4pi*x)

xi = linspace(x[1],x[end])
yi = interp1d(y)(x,xi)

figure(1)
clf()
plot(xi,yi)
plot(x,y,".")
margins(0,0.1)

savefig("example.png", transparent=true)


text:
using DataValues
using Test

@testset "DataValueArray: DataVector" begin

    A = [1:10...]
    B = [1:5...]

    X = DataValueArray(A)
    Y = DataValueArray(B)

# --- Base.push!

# Base.push!{T, V}(X::DataValueVector{T}, v::V)
    Z = DataValueVector{Int}()
    push!(Z, 5)
    @test isequal(Z[1], DataValue(5))

# Base.push!{T, V}(X::DataValueVector{T}, v::DataValue{V})
    push!(Z, DataValue())
    @test isequal(Z, DataValueArray([5, 1], [false, true]))
    push!(Z, DataValue(5))
    @test isequal(Z, DataValueArray([5, 1, 5], [false, true, false]))

    Z2 = copy(Z)
    push!(Z2, DataValue{Int}())
    @test isequal(Z2, DataValueArray([5, NA, 5, NA]))

    ZZ = DataValueVector{Int}()
    push!(ZZ, 5, DataValue(), DataValue(5))
    @test isequal(ZZ, Z)

# --- Base.pop!

# Base.pop!(X::DataValueArray)
    @test isequal(pop!(Z), DataValue(5))
    @test isequal(pop!(Z), DataValue{Int}())
    @test isequal(pop!(Z), DataValue(5))

# --- Base.pushfirst!

# Base.pushfirst!(X::DataValueVector, v)
    @test isequal(pushfirst!(X, 3), DataValueArray(vcat(3, [1:10...])))

# Base.pushfirst!(X::DataValueVector, v::DataValue)
    @test isequal(pushfirst!(X, DataValue(2)),
                DataValueArray(vcat([2, 3], [1:10...])))
    @test isequal(pushfirst!(X, DataValue{Int}()),
                DataValueArray(vcat([1, 2, 3], [1:10...]),
                            vcat(true, fill(false, 12))))

    @test isequal(pushfirst!(Y, 1, DataValue(), DataValue(3)),
                DataValueArray([1, 2, 3, 1, 2, 3, 4, 5],
                            [false, true, false, false,
                            false, false, false, false]))

# --- Base.popfirst!

# Base.popfirst!{T}(X::DataValueVector{T})
    Z = DataValueArray([1:10...])

    @test isequal(popfirst!(Z), DataValue(1))
    @test isequal(Z, DataValueArray([2:10...]))

    pushfirst!(Z, DataValue{Int}())

    @test isequal(popfirst!(Z), DataValue{Int}())

# --- test Base.splice!

# Base.splice!(X::DataValueVector, i::Integer, ins=_default_splice)
    A = [1:10...]
    B = [1:10...]
    X = DataValueArray(B)

    i = rand(1:10)
    @test isequal(splice!(X, i), DataValue(splice!(A, i)))
    @test isequal(X, DataValueArray(A))

    i = rand(1:9)
    j = rand(1:9)
    @test isequal(splice!(X, i, j), DataValue(splice!(A, i, j)))
    @test isequal(X, DataValueArray(A))

    A = [1:10...]
    B = [1:10...]
    X = DataValueArray(B)
    i = rand(1:5)
    n = rand(3:5)
    @test isequal(splice!(X, i, [1:n...]), DataValue(splice!(A, i, [1:n...])))
    @test isequal(X, DataValueArray(A))

# Base.splice!{T<:Integer}(X::DataValueVector,
#                          rng::UnitRange{T},
#                          ins=_default_splice)

# test with length(rng) > length(ins)
    A = [1:20...]
    B = [1:20...]
    X = DataValueArray(B)
    f = rand(1:7)
    d = rand(3:5)
    l = f + d
    ins = [1, 2]
    @test isequal(splice!(X, f:l, ins),
                DataValueArray(splice!(A, f:l, ins)))
    @test isequal(X, DataValueArray(A))

    i = rand(1:length(X))
    @test isequal(splice!(X, 1:i), DataValueArray(splice!(A, 1:i)))

    A = [1:20...]
    B = [1:20...]
    X = DataValueArray(B)
    f = rand(10:15)
    d = rand(3:5)
    l = f + d
    ins = [1, 2]
    @test isequal(splice!(X, f:l, ins),
                DataValueArray(splice!(A, f:l, ins)))
    @test isequal(X, DataValueArray(A))

# test with length(rng) < length(ins)
    A = [1:20...]
    B = [1:20...]
    X = DataValueArray(B)
    f = rand(1:7)
    d = rand(3:5)
    l = f + d
    ins = [1, 2, 3, 4, 5, 6, 7]
    @test isequal(splice!(X, f:l, ins),
                DataValueArray(splice!(A, f:l, ins)))
    @test isequal(X, DataValueArray(A))

    A = [1:20...]
    B = [1:20...]
    X = DataValueArray(B)
    f = rand(10:15)
    d = rand(3:5)
    l = f + d
    ins = [1, 2, 3, 4, 5, 6, 7]
    @test isequal(splice!(X, f:l, ins),
                DataValueArray(splice!(A, f:l, ins)))
    @test isequal(X, DataValueArray(A))

# --- test Base.deleteat!

# Base.deleteat!(X::DataValueArray, inds)
    X = DataValueArray(1:10)
    @test isequal(deleteat!(X, 1), DataValueArray(2:10))

# --- test Base.append!

# Base.append!(X::DataValueVector, items::AbstractVector)
    @test isequal(append!(X, [11, 12]),
                DataValueArray(2:12))
    @test isequal(append!(X, [DataValue(13), DataValue(14)]),
                DataValueArray(2:14))
    @test isequal(append!(X, [DataValue(15), DataValue{Int}()]),
                DataValueArray(DataValue{Int}[2:15..., DataValue()]))

# --- test Base.prepend!

# Base.prepend!(X::DataValueVector, items::AbstractVector)

    X = DataValueArray(3:12)
    @test isequal(prepend!(X, [1, 2]),
                DataValueArray(1:12))
    @test isequal(prepend!(X, [DataValue(-1), DataValue(0)]),
                DataValueArray(-1:12))
    @test isequal(prepend!(X, [DataValue{Int}(), DataValue(-2)]),
                DataValueArray(DataValue{Int}[DataValue(), -2:12...]))

    X2 = DataValueArray([1,NA,4])
    prepend!(X2, X2)
    @test isequal(X2, DataValueArray([1,NA,4,1,NA,4]))

# --- test Base.sizehint!

# Base.sizehint!(X::DataValueVector, newsz::Integer)
    sizehint!(X, 20)

# --- test padna!

# padna!{T}(X::DataValueVector{T}, front::Integer, back::Integer)
    X = DataValueArray(1:5)
    padna!(X, 2, 3)
    @test length(X.values) == 10
    @test X.isna == vcat(true, true, fill(false, 5), true, true, true)

# padna(X::DataValueVector, front::Integer, back::Integer)
    X = DataValueArray(1:5)
    Y = padna(X, 2, 3)
    @test length(Y.values) == 10
    @test Y.isna == vcat(true, true, fill(false, 5), true, true, true)

# --- test Base.reverse!/Base.reverse

    y = DataValueArray([NA, 2, 3, 4, NA, 6])
    @assert isequal(reverse(y),
                DataValueArray([6, NA, 4, 3, 2, NA]))

# check case where only nothing occurs in final position
    @assert isequal(unique(DataValueArray([1, 2, 1, NA])),
                DataValueArray([1, 2, NA]))

# Base.reverse!(X::DataValueVector, s=1, n=length(X))
# check for case where isbitstype(eltype(X)) = false
    Z = DataValueArray(Array{Int,1}[[1, 2], [3, 4], [5, 6]])
    @test isequal(reverse!(Z),
                DataValueArray(Array{Int,1}[[5, 6], [3, 4], [1, 2]]))

# Base.reverse!(X::DataValueVector, s=1, n=length(X))
# check for case where isbitstype(eltype(X)) = false & any(isna, X) = true
    A = fill([1,2], 20)
    Z = DataValueArray{Array{Int,1}}(20)
    i = rand(2:7)
    for i in [i - 1, i, i + 1, 20 - (i + 2), 20 - (i - 1)]
        Z[i] = [1, 2]
    end
    vals = Z.values
    nulls = Z.isna
    @test isequal(reverse!(Z), DataValueArray(A, reverse!(Z.isna)))

## empty!(X::DataValueVector)
    n = rand(1:1_000)
    A, M = rand(n), rand(Bool, n)
    X = DataValueArray(A, M)
    empty!(X)
    @test isempty(X)

end


text:
module Threshold

using CUDA

include("poisson.jl")

end


text:
"""
    check_duplicate_reaction(rxn::Reaction, rxns::Dict{String, Reaction}; only_metabolites=true)

Check if `rxn` already exists in `rxns` but has another `id`.
If `only_metabolites` is `true` then only the metabolite `id`s are checked.
Otherwise, compares metabolite `id`s and the absolute value of their stoichiometric coefficients to those of `rxn`.
If `rxn` has the same reaction equation as another reaction in `rxns`, the return the `id`.
Otherwise return `nothing`.

See also: [`reaction_mass_balanced`](@ref)
"""
function check_duplicate_reaction(
    crxn::Reaction,
    rxns::OrderedDict{String,Reaction};
    only_metabolites = true,
)
    for (k, rxn) in rxns
        if rxn.id != crxn.id # skip if same ID
            if only_metabolites # only check if metabolites are the same
                if issetequal(keys(crxn.metabolites), keys(rxn.metabolites))
                    return k
                end
            else # also check the stoichiometric coefficients
                reaction_checker = true
                for (kk, vv) in rxn.metabolites # get reaction stoich
                    if abs(get(crxn.metabolites, kk, 0)) != abs(vv) # if at least one stoich doesn't match
                        reaction_checker = false
                        break
                    end
                end
                if reaction_checker &&
                   issetequal(keys(crxn.metabolites), keys(rxn.metabolites))
                    return k
                end
            end
        end
    end
    return nothing
end

"""
    is_boundary(rxn_dict::Dict{String, Float64})

Return true if the reaction denoted by `rxn_dict` is a boundary reaction, otherwise return false.
Checks if on boundary by inspecting the number of metabolites in `rxn_dict`.
Boundary reactions have only one metabolite, e.g. an exchange reaction, or a sink/demand reaction.
"""
function is_boundary(rxn_dict::Dict{String,Float64})::Bool
    length(keys(rxn_dict)) == 1
end

is_boundary(model::MetabolicModel, rxn_id::String) =
    is_boundary(reaction_stoichiometry(model, rxn_id))

is_boundary(rxn::Reaction) = is_boundary(rxn.metabolites)

is_boundary(model::StandardModel, rxn::Reaction) = is_boundary(rxn) # for consistency with functions below

"""
    reaction_atom_balance(model::StandardModel, rxn)

Returns a dictionary mapping the stoichiometry of atoms through a single reaction. Uses the
metabolite information in `model` to determine the mass balance. Accepts a reaction
dictionary, a reaction string id or a `Reaction` as an argument for `rxn`.

See also: [`reaction_mass_balanced`](@ref)
"""
function reaction_atom_balance(model::StandardModel, reaction_dict::Dict{String,Float64})
    atom_balances = Dict{String,Float64}()
    for (met, stoich_rxn) in reaction_dict
        adict = metabolite_formula(model, met)
        isnothing(adict) &&
            throw(ErrorException("Metabolite $met does not have a formula assigned to it."))
        for (atom, stoich_molecule) in adict
            atom_balances[atom] =
                get(atom_balances, atom, 0.0) + stoich_rxn * stoich_molecule
        end
    end
    return atom_balances
end

function reaction_atom_balance(model::StandardModel, rxn_id::String)
    reaction_atom_balance(model, model.reactions[rxn_id].metabolites)
end

reaction_atom_balance(model::StandardModel, rxn::Reaction) =
    reaction_atom_balance(model, rxn.id)

"""
    reaction_mass_balanced(model::StandardModel, rxn)

Checks if `rxn` is atom balanced. Returns a boolean for whether the reaction is balanced,
and the associated balance of atoms for convenience (useful if not balanced). Calls
`reaction_atom_balance` internally.

See also: [`check_duplicate_reaction`](@ref), [`reaction_atom_balance`](@ref)
"""
reaction_mass_balanced(model::StandardModel, rxn_id::String) =
    all(values(reaction_atom_balance(model, rxn_id)) .== 0)

reaction_mass_balanced(model::StandardModel, rxn::Reaction) =
    reaction_mass_balanced(model, rxn.id)

reaction_mass_balanced(model::StandardModel, reaction_dict::Dict{String,Float64}) =
    all(values(reaction_atom_balance(model, reaction_dict)) .== 0)

"""
    stoichiometry_string(rxn_dict::Dict{String, Float64}; format_id = x -> x)

Return the reaction equation as a string. The metabolite strings can be manipulated by
setting `format_id`.

# Example
```
julia> req = Dict("coa_c" => -1, "for_c" => 1, "accoa_c" => 1, "pyr_c" => -1)
julia> stoichiometry_string(req)
"coa_c + pyr_c = for_c + accoa_c"

julia> stoichiometry_string(req; format_id = x -> x[1:end-2])
"coa + pyr = for + accoa"
```
"""
function stoichiometry_string(req; format_id = x -> x)
    count_prefix(n) = abs(n) == 1 ? "" : string(abs(n), " ")
    substrates =
        join((string(count_prefix(n), format_id(met)) for (met, n) in req if n < 0), " + ")
    products =
        join((string(count_prefix(n), format_id(met)) for (met, n) in req if n >= 0), " + ")
    return substrates * " = " * products
end

"""
    stoichiometry_string(rxn::Reaction; kwargs)

Alternative of [`stoichiometry_string`](@ref) take takes a `Reaction` as an argument.
"""
stoichiometry_string(rxn::Reaction; kwargs...) =
    stoichiometry_string(rxn.metabolites; kwargs...)


text:
"A wrapper for one value in a PostgreSQL result."
struct PQValue{OID}
    "PostgreSQL result"
    jl_result::Result

    "Row index of the result (0-indexed)"
    row::Cint

    "Column index of the result (0-indexed)"
    col::Cint

    function PQValue{OID}(jl_result::Result, row::Integer, col::Integer) where OID
        return new{OID}(jl_result, row - 1, col - 1)
    end
end

"""
    PQValue(jl_result::Result, row::Integer, col::Integer) -> PQValue
    PQValue{OID}(jl_result::Result, row::Integer, col::Integer) -> PQValue{OID}

Construct a `PQValue` wrapping one value in a PostgreSQL result.
Row and column positions are provided 1-indexed.
If the `OID` type parameter is not provided, the Oid of the field will be retrieved from
the result.
"""
function PQValue(jl_result::Result, row::Integer, col::Integer)
    oid = libpq_c.PQftype(jl_result.result, col - 1)

    return PQValue{oid}(jl_result, row, col)
end

"""
    isnull(jl_result::Result, row::Integer, col::Integer) -> Bool

Return whether the result value at the specified row and column (1-indexed) is `NULL`.
"""
function isnull(jl_result::Result, row::Integer, col::Integer)
    return libpq_c.PQgetisnull(jl_result.result, row - 1, col - 1) == 1
end

"""
    num_bytes(pqv::PQValue) -> Cint

The length in bytes of the `PQValue`'s corresponding data.
LibPQ.jl currently always uses text format, so this is equivalent to C's `strlen`.

See also: [`data_pointer`](@ref)
"""
num_bytes(pqv::PQValue) = libpq_c.PQgetlength(pqv.jl_result.result, pqv.row, pqv.col)

"""
    data_pointer(pqv::PQValue) -> Ptr{UInt8}

Get a raw pointer to the data for one value in a PostgreSQL result.
This data will be freed by libpq when the result is cleared, and should only be used
temporarily.
"""
data_pointer(pqv::PQValue) = libpq_c.PQgetvalue(pqv.jl_result.result, pqv.row, pqv.col)

"""
    unsafe_string(pqv::PQValue) -> String

Construct a `String` from a `PQValue` by copying the data.
"""
function Base.unsafe_string(pqv::PQValue)
    return unsafe_string(data_pointer(pqv), num_bytes(pqv))
end

"""
    string_view(pqv::PQValue) -> String

Wrap a `PQValue`'s underlying data in a `String`.
This function uses [`data_pointer`](@ref) and [`num_bytes`](@ref) and does not copy.

!!! note

    The underlying data will be freed by libpq when the result is cleared, and should only
    be used temporarily.

See also: [`bytes_view`](@ref)
"""
function string_view(pqv::PQValue)
    return String(unsafe_wrap(Vector{UInt8}, data_pointer(pqv), num_bytes(pqv)))
end

"""
    bytes_view(pqv::PQValue) -> Vector{UInt8}

Wrap a `PQValue`'s underlying data in a vector of bytes.
This function uses [`data_pointer`](@ref) and [`num_bytes`](@ref) and does not copy.

This function differs from [`string_view`](@ref) as it keeps the `\0` byte at the end.
`PQValue` parsing functions should use `bytes_view` when the data returned by PostgreSQL
is not in UTF-8.

!!! note

    The underlying data will be freed by libpq when the result is cleared, and should only
    be used temporarily.
"""
bytes_view(pqv::PQValue) = unsafe_wrap(Vector{UInt8}, data_pointer(pqv), num_bytes(pqv) + 1)

Base.String(pqv::PQValue) = unsafe_string(pqv)
Base.convert(::Type{String}, pqv::PQValue) = String(pqv)
Base.length(pqv::PQValue) = length(string_view(pqv))
Base.lastindex(pqv::PQValue) = lastindex(string_view(pqv))

# Fallback, because Base requires string iteration state to be indices into the string.
# In an ideal world, PQValue would be an AbstractString and this particular method would
# not be necessary.
"""
    parse(::Type{T}, pqv::PQValue) -> T

Parse a value of type `T` from a `PQValue`.
By default, this uses any existing `parse` method for parsing a value of type `T` from a
`String`.

You can implement default PostgreSQL-specific parsing for a given type by overriding
`pqparse`.
"""
Base.parse(::Type{T}, pqv::PQValue) where {T} = pqparse(T, string_view(pqv))

"""
    LibPQ.pqparse(::Type{T}, str::AbstractString) -> T

Parse a value of type `T` from any `AbstractString`.
This is used to parse PostgreSQL's output format.
"""
function pqparse end

# Fallback method
pqparse(::Type{T}, str::AbstractString) where {T} = parse(T, str)

# allow parsing as a Symbol anything which works as a String
pqparse(::Type{Symbol}, str::AbstractString) = Symbol(str)

## integers
_DEFAULT_TYPE_MAP[:int2] = Int16
_DEFAULT_TYPE_MAP[:int4] = Int32
_DEFAULT_TYPE_MAP[:int8] = Int64

## floating point
_DEFAULT_TYPE_MAP[:float4] = Float32
_DEFAULT_TYPE_MAP[:float8] = Float64

## oid
_DEFAULT_TYPE_MAP[:oid] = Oid

## numeric
_DEFAULT_TYPE_MAP[:numeric] = Decimal

# no default for monetary; needs lconv and lc_monetary from result/connection

## character
# bpchar is char(n)
function pqparse(::Type{String}, str::AbstractString)
    return String(rstrip(str, ' '))
end
# char is "char"
_DEFAULT_TYPE_MAP[:char] = PQChar
pqparse(::Type{PQChar}, str::AbstractString) = PQChar(first(str))
pqparse(::Type{Char}, str::AbstractString) = Char(pqparse(PQChar, str))
# varchar, text, and name are all String

## binary data

_DEFAULT_TYPE_MAP[:bytea] = Vector{UInt8}

# Needs it's own `parse` method as it uses bytes_view instead of string_view
function Base.parse(::Type{Vector{UInt8}}, pqv::PQValue{PQ_SYSTEM_TYPES[:bytea]})
    pqparse(Vector{UInt8}, bytes_view(pqv))
end

function pqparse(::Type{Vector{UInt8}}, bytes::Array{UInt8,1})
    byte_length = Ref{Csize_t}(0)

    unescaped_ptr = libpq_c.PQunescapeBytea(bytes, byte_length)

    if unescaped_ptr == C_NULL
        error("Could not unescape byte sequence $(String(bytes))")
    end

    unescaped_vec = copy(unsafe_wrap(Vector{UInt8}, unescaped_ptr, byte_length[]))

    libpq_c.PQfreemem(unescaped_ptr)

    return unescaped_vec
end

## bool
# TODO: check whether we ever need this or if PostgreSQL always gives t or f
_DEFAULT_TYPE_MAP[:bool] = Bool
const BOOL_TRUE = r"^\s*(t|true|y|yes|on|1)\s*$"i
const BOOL_FALSE = r"^\s*(f|false|n|no|off|0)\s*$"i
function pqparse(::Type{Bool}, str::AbstractString)
    if occursin(BOOL_TRUE, str)
        return true
    elseif occursin(BOOL_FALSE, str)
        return false
    else
        error("\"$str\" is not a valid boolean")
    end
end

## dates and times
# ISO, YMD

# Cut off digits after the third after the decimal point,
# since DateTime in Julia currently handles only milliseconds
# see https://github.com/invenia/LibPQ.jl/issues/33
_trunc_seconds(str) = replace(str, r"(\.[\d]{3})\d+" => s"\g<1>")

_DEFAULT_TYPE_MAP[:timestamp] = DateTime
const TIMESTAMP_FORMAT = dateformat"y-m-d HH:MM:SS.s"  # .s is optional here
function pqparse(::Type{DateTime}, str::AbstractString)
    if str == "infinity"
        return typemax(DateTime)
    elseif str == "-infinity"
        return typemin(DateTime)
    end

    # Cut off digits after the third after the decimal point,
    # since DateTime in Julia currently handles only milliseconds, see Issue #33
    str = replace(str, r"(\.[\d]{3})\d+" => s"\g<1>")
    return parse(DateTime, str, TIMESTAMP_FORMAT)
end

# ISO, YMD
_DEFAULT_TYPE_MAP[:timestamptz] = ZonedDateTime
const TIMESTAMPTZ_FORMATS = (
    dateformat"y-m-d HH:MM:SSz",
    dateformat"y-m-d HH:MM:SS.sz",
    dateformat"y-m-d HH:MM:SS.ssz",
    dateformat"y-m-d HH:MM:SS.sssz",
)
function pqparse(::Type{ZonedDateTime}, str::AbstractString)
    if str == "infinity"
        return ZonedDateTime(typemax(DateTime), tz"UTC")
    elseif str == "-infinity"
        return ZonedDateTime(typemin(DateTime), tz"UTC")
    end

    for fmt in TIMESTAMPTZ_FORMATS[1:end-1]
        parsed = tryparse(ZonedDateTime, str, fmt)
        parsed !== nothing && return parsed
    end

    return parse(ZonedDateTime, _trunc_seconds(str), TIMESTAMPTZ_FORMATS[end])
end

_DEFAULT_TYPE_MAP[:date] = Date
function pqparse(::Type{Date}, str::AbstractString)
    if str == "infinity"
        return typemax(Date)
    elseif str == "-infinity"
        return typemin(Date)
    end

    return parse(Date, str)
end

_DEFAULT_TYPE_MAP[:time] = Time
function pqparse(::Type{Time}, str::AbstractString)
    try
        return parse(Time, str)
    catch err
        if !(err isa InexactError)
            rethrow(err)
        end
    end

    return parse(Time, _trunc_seconds(str))
end

# UNIX timestamps
function Base.parse(::Type{DateTime}, pqv::PQValue{PQ_SYSTEM_TYPES[:int8]})
    unix2datetime(parse(Int64, pqv))
end

function Base.parse(::Type{ZonedDateTime}, pqv::PQValue{PQ_SYSTEM_TYPES[:int8]})
    TimeZones.unix2zdt(parse(Int64, pqv))
end

## intervals
# iso_8601
_DEFAULT_TYPE_MAP[:interval] = Dates.CompoundPeriod
const INTERVAL_REGEX = Ref{Regex}()  # set at __init__

function _interval_regex()
    function _field_match(period_type, number_match="-?\\d+")
        name = nameof(period_type)
        letter = first(String(name))
        return "(?:(?<$name>$number_match)$letter)?"
    end

    io = IOBuffer()
    print(io, "^P")
    for long_type in (Year, Month, Day)
        print(io, _field_match(long_type))
    end
    print(io, "(?:T")
    for long_type in (Hour, Minute)
        print(io, _field_match(long_type))
    end
    print(io,
        _field_match(Second, "(?<whole_seconds>-?\\d+)(?:\\.(?<frac_seconds>\\d{1,9}))?"),
    )
    print(io, ")?\$")

    return Regex(String(take!(io)))
end

# parse the iso_8601 interval output format
# https://www.postgresql.org/docs/10/datatype-datetime.html#DATATYPE-INTERVAL-OUTPUT
function pqparse(::Type{Dates.CompoundPeriod}, str::AbstractString)
    interval_regex = INTERVAL_REGEX[]
    matched = match(interval_regex, str)

    if matched === nothing
        error("Couldn't parse $str as interval using regex $interval_regex")
    end

    periods = Period[]
    sizehint!(periods, 7)
    for period_type in (Year, Month, Day, Hour, Minute)
        period_str = matched[nameof(period_type)]
        if period_str !== nothing
            push!(periods, period_type(parse(Int, period_str)))
        end
    end

    if matched["Second"] !== nothing
        whole_seconds_str = matched["whole_seconds"]
        whole_seconds = parse(Int, whole_seconds_str)
        if whole_seconds != 0
            push!(periods, Second(whole_seconds))
        end

        #=
        We need to parse the fractional seconds as a period.
        Here we try to keep to the largest period type possible for representing the
        fractional seconds.

        For example, 1 is 100 Milliseconds, but 0001 is 100 Microseconds
        =#
        frac_seconds_str = matched["frac_seconds"]
        if frac_seconds_str !== nothing
            len = length(frac_seconds_str)
            frac_periods = [Millisecond, Microsecond, Nanosecond]
            period_coeff = fld1(len, 3)
            period_type = frac_periods[period_coeff]  # field regex prevents BoundsError

            frac_seconds = parse(Int, frac_seconds_str) * 10 ^ (3 * period_coeff - len)
            if frac_seconds != 0
                push!(periods, period_type(frac_seconds))
            end
        end
    end

    return Dates.CompoundPeriod(periods)
end

## ranges
_DEFAULT_TYPE_MAP[:int4range] = Interval{Int32}
_DEFAULT_TYPE_MAP[:int8range] = Interval{Int64}
_DEFAULT_TYPE_MAP[:numrange] = Interval{Decimal}
_DEFAULT_TYPE_MAP[:tsrange] = Interval{DateTime}
_DEFAULT_TYPE_MAP[:tstzrange] = Interval{ZonedDateTime}
_DEFAULT_TYPE_MAP[:daterange] = Interval{Date}

# Matches anything but the start or end of an interval or a comma
const RANGE_ITEM = "[^\\[\\(\\]\\),]+"
# Makes sure the string starts and ends with a bracket or parentheses, has two items in the
# interval, and a comma to separate the two items.
const RANGE_REGEX = Regex("^([\\[\\(])($RANGE_ITEM),($RANGE_ITEM)([\\]\\)])\$")
get_inclusivity(ch) = ch in ("[", "]") ? true : false

function pqparse(::Type{Interval{T}}, str::AbstractString) where {T}
    if str == "empty"
        return Interval{T}()
    end

    matched = match(RANGE_REGEX, str)
    if matched === nothing
        error("Couldn't parse $str as range using regex $RANGE_REGEX")
    else
        matched = matched.captures
    end

    start_inclusivity = get_inclusivity(matched[1])
    end_inclusivity = get_inclusivity(matched[4])

    # Datetime formats have quotes around them so we strip those out
    start = pqparse(T, strip(matched[2], ['"']))
    endpoint = pqparse(T, strip(matched[3], ['"']))

    return Interval(start, endpoint, start_inclusivity, end_inclusivity)
end

## arrays
# numeric arrays never have double quotes and always use ',' as a separator
parse_numeric_element(::Type{T}, str) where T = parse(T, str)

parse_numeric_element(::Type{Union{T, Missing}}, str) where T =
    str == "NULL" ? missing : parse(T, str)

function parse_numeric_array(eltype::Type{T}, str::AbstractString) where T
    eq_ind = findfirst(isequal('='), str)

    if eq_ind !== nothing
        offset_str = str[1:eq_ind-1]
        range_strs = split(str[1:eq_ind-1], ['[',']']; keepempty=false)

        ranges = map(range_strs) do range_str
            lower, upper = split(range_str, ':'; limit=2)
            return parse(Int, lower):parse(Int, upper)
        end

        arr = OffsetArray{T}(undef, ranges...)
        el_iter = eachmatch(r"[^\}\{,]+", str[eq_ind+1:end])
    else
        arr = Array{T}(undef, array_size(str)...)
        el_iter = eachmatch(r"[^\}\{,]+", str)
    end

    idx_iter = imap(reverse, product(reverse(axes(arr))...))
    for (idx, num_match) in zip(idx_iter, el_iter)
        arr[idx...] = parse_numeric_element(T, num_match.match)
    end

    return arr
end

function array_size(str)
    ndims = something(findfirst(c -> c != '{', str), 0) - 1
    dims = zeros(Int, ndims)

    curr_dim = ndims
    curr_pos = ndims
    open_braces = ndims
    last_ind = lastindex(str)
    el_count = 0
    while curr_dim > 0 && curr_pos < last_ind
        curr_pos = nextind(str, curr_pos)

        if str[curr_pos] == '}'
            open_braces -= 1

            if open_braces < curr_dim
                dims[curr_dim] = el_count
                curr_dim -= 1
                el_count = 1
            end
        elseif str[curr_pos] == '{'
            open_braces += 1
        elseif str[curr_pos] == ','
            if open_braces == curr_dim
                el_count += 1
            end
        else
            if open_braces == curr_dim && el_count == 0
                el_count = 1
            end
        end
    end

    return dims
end

for pq_eltype in ("int2", "int4", "int8", "float4", "float8", "oid", "numeric")
    array_oid = PQ_SYSTEM_TYPES[Symbol("_$pq_eltype")]
    jl_type = _DEFAULT_TYPE_MAP[Symbol(pq_eltype)]
    jl_missingtype = Union{jl_type, Missing}

    # could be an OffsetArray or Array of any dimensionality
    _DEFAULT_TYPE_MAP[array_oid] = AbstractArray{jl_missingtype}

    for jl_eltype in (jl_type, jl_missingtype)
        @eval function pqparse(
            ::Type{A}, str::AbstractString
        ) where A <: AbstractArray{$jl_eltype}
            parse_numeric_array($jl_eltype, str)::A
        end
    end
end

struct FallbackConversion <: AbstractDict{Tuple{Oid, Type}, Base.Callable}
end

function Base.getindex(cmap::FallbackConversion, oid_typ::Tuple{Integer, Type})
    _, typ = oid_typ

    return function parse_type(pqv::PQValue)
        parse(typ, pqv)
    end
end

Base.haskey(cmap::FallbackConversion, oid_typ::Tuple{Integer, Type}) = true

"""
A fallback conversion mapping (like [`PQConversions`](@ref) which holds a single function
for converting PostgreSQL data of a given Oid to a given Julia type, using the [`parse`](@ref)
function.
"""
const _FALLBACK_CONVERSION = FallbackConversion()