text:
# OP1
"""
`faceswith(mesh::Mesh, idx)`

Return a vector `I` of the indices of faces containing vertex with given `idx`.
"""
faceswith(mesh::Mesh, idx) =
  findall(face -> idx in face, mesh.faces)

function faceswith(mesh::HalfEdgeMesh, idx)
  faces = UInt32[]
  itersameorigin(mesh, idx) do edge
    push!(faces, edge.face)
  end
  faces
end

# OP2
"""
`vneighbors(mesh::Mesh, idx; depth=1)`

Return a vector `I` containing the indices of vertices surrounding vertex with given `idx`.
"""
function vneighbors(mesh::Mesh, idx; depth=1)
  neighbors = Set{UInt32}()
  distance = zerobased(fill(depth + 1, length(mesh.vertices)))
  usedfaces = fill(false, length(mesh.faces))
  distance[idx] = 0

  for depth in 1:depth, (i, face) in enumerate(mesh.faces)
    usedfaces[i] && continue
    dist = minimum(v -> distance[v], face)
    if dist < depth
      usedfaces[i] = true
      for v in face
        distance[v] = min(distance[v], dist + 1) # works only for triangles
        push!(neighbors, v)
      end
    end
  end
  delete!(neighbors, idx)

  collect(neighbors)
end

function vneighbors(mesh::HalfEdgeMesh, idx; depth=1)
  neighbors = Set{UInt32}()
  distance = zerobased(fill(depth + 1, length(mesh.vertices)))
  distance[idx] = 0

  queue = Queue{UInt32}()
  enqueue!(queue, idx)
  while !isempty(queue)
    idxv = dequeue!(queue)
    itersameorigin(mesh, idxv) do edge
      for vert in (mesh.halfedges[edge.next].origin, mesh.halfedges[edge.prev].origin)
        if distance[idxv] + 1 < distance[vert]
          distance[vert] = distance[idxv] + 1
          push!(neighbors, vert)
          if distance[vert] != depth
            enqueue!(queue, vert)
          end
        end
      end
    end
  end
  delete!(neighbors, idx)

  collect(neighbors)
end

# OP3
"""
`fneighbors(mesh::Mesh, idx; depth=1)`

Return a vector `I` containing the indices of faces surrounding face with given `idx`.
"""
function fneighbors(mesh::TriangleMesh, idx; depth=1)
  neighbors = Set{UInt32}()
  distance = zerobased(fill(depth + 1, length(mesh.vertices)))
  usedfaces = fill(false, length(mesh.faces))

  distance[face(mesh, idx)] .= 0
  usedfaces[idx] = true

  for depth in 1:depth, (i, face) in enumerate(mesh.faces)
    usedfaces[i] && continue
    dist = sort(map(v -> distance[v], face))
    if dist[2] < depth
      usedfaces[i] = true
      push!(neighbors, i)
      for v in face
        distance[v] = min(distance[v], dist[2] + 1) # works only for triangles
      end
    end
  end
  delete!(neighbors, idx)

  collect(neighbors)
end

function fneighbors(mesh::HalfEdgeMesh, idx; depth=1)
  neighbors = Set{UInt32}()
  distance = fill(depth + 1, length(mesh.faces))
  distance[idx] = 0

  queue = Queue{UInt32}()
  enqueue!(queue, idx)
  while !isempty(queue)
    facev = dequeue!(queue)
    itersameface(mesh, facev) do edge
      edge.twin != 0 || return # next edge
      twin = mesh.halfedges[edge.twin].face
      if distance[facev] + 1 < distance[twin]
        distance[twin] = distance[facev] + 1
        push!(neighbors, twin)
        if distance[twin] != depth
          enqueue!(queue, twin)
        end
      end
    end
  end
  delete!(neighbors, idx)

  collect(neighbors)
end

# OP4
"""
`traverse(mesh::Mesh{3,Point})`
"""
traverse(mesh::Mesh{3,Point}, faceidx, target::Point) where Point =
Channel{Int32}() do channel
  # compute surface norms
  usedfaces = fill(false, length(mesh.faces))

  @inbounds while true
    @label start
    usedfaces[faceidx] && break
    usedfaces[faceidx] = true
    facev = face(mesh, faceidx)
    put!(channel, faceidx) # yield

    # find all edges
    edgepoints = [(vertex(mesh, S), vertex(mesh, T)) for (S, T) in edges(facev)]

    # compute edge norms
    edgenorms = [normalize((T .- S) × facenorm(mesh, faceidx)) for (S, T) in edgepoints]
    # compute directional vecotrs
    dirs = [normalize(target .- midpoint(S, T)) for (S, T) in edgepoints]

    # compute dot products
    dotprods = [e ⋅ d for (e, d) in zip(edgenorms, dirs)] # e,d normalized
    # find smallest angle
    ord = sortperm(dotprods; rev=true) # cosθ₁ < cosθ₂ => θ₁ > θ₂

    dotprods[ord[1]] <= 0 && break # point is inside the current face
    for v in ord
      S, T = facev[[v, v % 3 + 1]]
      # check if edge has neighbouring face
      for newfaceidx in eachindex(mesh.faces)
        (newfaceidx == faceidx || usedfaces[newfaceidx]) && continue
        newface = mesh.faces[newfaceidx]
        if S in newface && T in newface
          faceidx = newfaceidx
          @goto start
        end
      end
    end
  end
end

# OP5
"""
`flipdiag!(mesh::AbstractMesh{N,V,Triangle}, faceidx₁, faceidx₂)`
"""
function flipdiag!(mesh::TriangleMesh, Aidx, Bidx)
  A = face(mesh, Aidx)
  B = face(mesh, Bidx)

  olddiag = A ∩ B
  @assert length(olddiag) == 2 "Face $Aidx should have exacly one edge in common with face $Bidx."

  newdiag = (setdiff(A, olddiag)[1], setdiff(B, olddiag)[1])

  mesh.faces[Aidx] = Triangle(newdiag[1], newdiag[2], olddiag[1])
  mesh.faces[Bidx] = Triangle(olddiag[2], newdiag[2], newdiag[1])
  newdiag
end

# OP4
traverse(mesh::HalfEdgeMesh{3,Point}, faceidx, target::Point) where Point =
Channel{Int32}() do channel
  usedfaces = fill(false, length(mesh.faces))

  @inbounds while true
    usedfaces[faceidx] && break
    usedfaces[faceidx] = true
    halfedges, facev = face(mesh, faceidx; keephalfedges=true)
    put!(channel, faceidx) # yield

    # find all edges
    edgepoints = [(vertex(mesh, S), vertex(mesh, T)) for (S, T) in edges(facev)]

    # compute edge norms
    edgenorms = [normalize((T .- S) × facenorm(mesh, faceidx)) for (S, T) in edgepoints]
    # compute directional vecotrs
    dirs = [normalize(target .- midpoint(S, T)) for (S, T) in edgepoints]

    # compute dot products
    dotprods = [e ⋅ d for (e, d) in zip(edgenorms, dirs)] # e,d normalized
    # find smallest angle
    ord = sortperm(dotprods; rev=true) # cosθ₁ < cosθ₂ => θ₁ > θ₂

    dotprods[ord[1]] <= 0 && break # point is inside the current face
    for v in ord
      edge = halfedges[v]
      # check if edge has neighbouring face
      twin = mesh.halfedges[edge].twin
      twin == 0 && continue
      newfaceidx = mesh.halfedges[twin].face
      (newfaceidx == faceidx || usedfaces[newfaceidx]) && continue
      faceidx = newfaceidx
      break
    end
  end
end

# OP5
function flipdiag!(mesh::HalfEdgeMesh{N,V,Triangle}, Aidx, Bidx) where {N,V}
  # find common edge
  start = ait = mesh.faces[Aidx]
  found = false
  while ait != start
    edge = mesh.halfedges[ait]
    if edge.twin != 0 && mesh.halfedges[edge.twin].face == Bidx
      found = true
      break
    end
    ait = edge.next
  end

  @assert found "Face $Aidx should have exacly one edge in common with face $Bidx."

  # swap
  @inbounds begin
    bit = mesh.halfedges[ait].twin

    Aprev = mesh.halfedges[ait].prev
    Bprev = mesh.halfedges[bit].prev

    Anext = mesh.halfedges[ait].next
    Bnext = mesh.halfedges[bit].next

    mesh.vertexhalfedge[meshes.halfedges[ait].origin] = Bnext
    mesh.vertexhalfedge[meshes.halfedges[bit].origin] = Anext

    mesh.halfedges[ait].origin = mesh.halfedges[Aprev].origin
    mesh.halfedges[bit].origin = mesh.halfedges[Bprev].origin

    A = (Anext, ait, Bprev, Anext, ait)
    B = (Bnext, bit, Aprev, Bnext, bit)

    for i in 2:4
      mesh.halfedges[A[i]].next = A[i+1]
      mesh.halfedges[A[i]].prev = A[i-1]
      mesh.halfedges[B[i]].next = B[i+1]
      mesh.halfedges[B[i]].prev = B[i-1]
    end

    mesh.halfedges[Bprev].face = Aidx
    mesh.halfedges[Aprev].face = Bidx

    mesh.faces[Aidx] = ait
    mesh.faces[Bidx] = bit

    (mesh.halfedges[ait].origin, mesh.halfedges[bit].origin)
  end
end


text:
mutable struct RotationVar{T,N} <: AbstractExtDataAdv
    decfl::Any
    decint::Any
    function RotationVar(adv::Advection{T,N}) where {T,N}
        return new{T,N}(missing,missing)
    end
end

"""
    initcoef!(pv::RotationVar{T, Nsp, Nv}, self::AdvectionData{T, Nsp, Nv, Nsum})

Implementation of the interface function that is called at the begining of each advection
    This is implementation for Vlasov-Poisson equation

"""

function initcoef!(
    pv::RotationVar{T,N},
    self::AdvectionData{T,N, timeopt, NoTimeAlg},
) where {T,N, timeopt}
    st_cur, st_other = getst(self).perm
    mesh_cur = self.adv.t_mesh[st_cur]
    mesh_other = self.adv.t_mesh[st_other]
    sign = ( st_cur == 1 ) ? -1 : 1
    pv.decfl = sign * getcur_t(self)/step(mesh_cur) * mesh_other.points

    #  if isvelocitystate(self)
    #     #        println("sp trace init moins")
    #     pv.bufcur = (getcur_t(self) / step(mesh_v)) * mesh_sp.points
    # else
    #     mesh_sp = self.adv.t_mesh_sp[state_dim]
    #     #        println("sp trace init moins")
    #     pv.bufcur = (-getcur_t(self) / step(mesh_sp)) * mesh_v.points
    # end
end
function initcoef!(
    pv::RotationVar{T,2},
    self::Union{AdvectionData{T,2,timeopt,ABTimeAlg_ip},AdvectionData{T,2,timeopt,ABTimeAlg_new}},
) where {T, timeopt}
    adv = self.adv
    sz = sizeall(adv)
    buf1 = -getcur_t(self)/step(adv.t_mesh[1]) * adv.t_mesh[2].points
    buf2 = getcur_t(self)/step(adv.t_mesh[2]) * adv.t_mesh[1].points
    if ismissing(self.bufcur)
        self.bufcur = zeros(OpTuple{2,T}, sz)
    end
    for ind in CartesianIndices(sizeall(adv))
        self.bufcur[ind] = OpTuple((buf1[ind.I[2]], buf2[ind.I[1]]))
    end

 #   @show self.bufcur
end

    #  if isvelocitystate(self)
    #     #        println("sp trace init moins")
    #     pv.bufcur = (getcur_t(self) / step(mesh_v)) * mesh_sp.points
    # else
    #     mesh_sp = self.adv.t_mesh_sp[state_dim]
    #     #        println("sp trace init moins")
    #     pv.bufcur = (-getcur_t(self) / step(mesh_sp)) * mesh_v.points
    # end

function getrotationvar(adv::Advection)
    return RotationVar(adv)
end

"""
    getalpha(pv::RotationVar, self::AdvectionData, ind)

Implementation of the interface function that is called before each interpolation in advection

"""
getalpha(pv::RotationVar, self::AdvectionData, ind) = (pv.decfl[ind],)

ext: jl


text:
convert_string(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject{String} = AetherlangObject(string(arg))
convert_str2num(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject{String})::AetherlangObject{Float64} = AetherlangObject(parse(Float64, arg.dataref))
convert_str2sym(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject{String})::AetherlangObject{String} = AetherlangObject(arg.dataref, true)
convert_num2int(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject{F}where F<:AbstractFloat)::AetherlangObject{I}where I<:Integer = AetherlangObject(Int(trunc(arg.dataref)))
function convert_char2int(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject{String})::AetherlangObject
    if length(arg.dataref) > 1
        throw(AetherlangError("char2num can only operate on strings of length 1"))
    end
    AetherlangObject(Int(arg.dataref[1]))
end
function convert_int2char(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject{I}where I<:Integer)::AetherlangObject
    AetherlangObject(Char(arg.dataref))
end

convert_namespace_modify = Namespace(
    "num2int" => AetherlangObject(convert_num2int),
    "string" => AetherlangObject(convert_string),
    "str2num" => AetherlangObject(convert_str2num),
    "str2sym" => AetherlangObject(convert_str2sym),
    "char2int" => AetherlangObject(convert_char2int),
    "int2char" => AetherlangObject(convert_int2char)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
#=
Model: MNA_5 (10913 variables, 9 inputs)
=#
using Reachability, MAT, SparseArrays

mna5(o::Pair{Symbol, <:Any}...) = mna5(Options(Dict{Symbol, Any}(o)))

function mna5(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    # dynamics
    file = matopen(@relpath "mna5.mat")
    A = sparse(read(file, "A"))
    n = size(A, 1)
    b = sparsevec(19:27, [fill(-0.1, 5); fill(-0.2, 4)], n)

    # initial set
    X0 = Hyperrectangle([fill(0.000225, 10); zeros(10903)],
                        [fill(0.000025, 10); zeros(10903)])

    # instantiate continuous LTI system
    S = IVP(CACS(A, b, nothing), X0)

    # property: x1 < 0.2 && x2 < 0.15
    property = SafeStatesProperty(HPolyhedron([
        HalfSpace(sparsevec([1], [1.0], n), 0.2),
        HalfSpace(sparsevec([2], [1.0], n), 0.15)]))

    # =======================
    # Problem default options
    # =======================
    # 2D blocks except for the last block, which is 1D
    partition = vcat([(2*i-1:2*i) for i in 1:5456], [n:n])

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [1],
                                  :partition => partition,
                                  :plot_vars => [0, 1],
                                  :assume_sparse => true,
                                  :lazy_expm => true)

    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [1, 2],
                                  :partition => partition,
                                  :property => property,
                                  :assume_sparse => true,
                                  :lazy_expm => true)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using HorizonSideRobots

function obhod(r)
    c=0::Int
    count=0::Int
    while !isborder(r,Ost)
        move!(r,Ost)
        if ismarker(r)
            c += temperature(r)
            count += 1
        end
    end
    if !isborder(r,Nord)
        move!(r,Nord)
        if ismarker(r)
            c += temperature(r)
            count += 1
        end
    end
    while !isborder(r,West)
        move!(r,West)
        if ismarker(r)
            c += temperature(r)
            count += 1
        end
    end
    return c,count
end

function main(r)
    x=0
    y=0
    i=0
    j=0
    while !isborder(r,Nord)
        x, y = obhod(r)
        i += x
        j += y
    end
    return print(i / j)
end

main(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
using StatsMonitor
using Test
using Sockets

opts = StatsMonitor.Config(interval=1)

@async StatsMonitor.tcp_server(opts)

s1() = """{ "myId": { "m": $(randn()), "tags": ["t1", "t2"] } }"""
s2() = """{ "myId2": { "c": 1, "tags": ["t2", "t3"] } }"""
s3() = """{ "myId3": { "s": $(rand(1:10)), "tags": ["t2", "t3"] } }"""

function test(n = 1000)
    sock = UDPSocket()

    for _ in 1:n
        send(sock, ip"127.0.0.1", 8125, s1())
        send(sock, ip"127.0.0.1", 8125, s2())
        send(sock, ip"127.0.0.1", 8125, s3())
    end

    close(sock)
end
test()

# function test2(n=10_000)
#     s1() = """{ "myId": { "mean": $(randn()), "tags": ["t1", "t2"] } }"""
#     s2() = """{ "myId2": { "omean": $(randn()), "tags": ["t2", "t3"] } }"""

#     sock = UDPSocket()
#     f1() = for i in 1:n; send(sock, ip"127.0.0.1", opts.port, s1()); end
#     f2() = for i in 1:n; send(sock, ip"127.0.0.1", opts.port, s2()); end

#     @time f1()
#     @time f1()
#     @time f2()
#     @time f2()
#     close(sock)
# end
# test2()


text:
"""
    rps10()

RPS Serial Chains for 10 positions on a circular hyperboloid

## References

Hai-Jun Su and J. Micheal McCarthy:
"Kinematics Synthesis of RPS Serial Chains",
In Proceedings of the ASME Design Engineering Technical Conferences
(CDROM).  Paper DETC03/DAC-48813.  Chicago, IL, Sept.02-06, 2003.
"""
function rps10()
    @polyvar g1 g2 g3 p1 p2 p3 q0 q1 q2 q3
    equations = [
    -0.1279703687075118*g1^2 - 0.48596123125526264*g1*g2 + 0.30699556370717496*g2^2 + 0.3778977698527674*g1*g3 - 0.23404544076569642*g2*g3 + 0.01563626178508072*g3^2 + 0.327228678790004*g1^2*p1 + 0.8426829275672494*g1*g2*p1 + 0.6075645757034159*g2^2*p1 - 1.1371405598667543*g1*g3*p1 + 0.229293271620915*g2*g3*p1 - 0.21948911177437957*g3^2*p1 - 0.2075154964282774*g1^2*p1^2 - 0.37702968479068544*g1*g2*p1^2 - 0.16688906819159421*g2^2*p1^2 + 0.7986954318323025*g1*g3*p1^2 + 0.866826144775651*g2*g3*p1^2 + 0.37440456461987165*g3^2*p1^2 + 1.5614616440131446*g1^2*p2 - 1.7388380675822595*g1*g2*p2 + 0.06790915713070725*g2^2*p2 - 0.4309121044684771*g1*g3*p2 + 0.9086272006283425*g2*g3*p2 - 0.2764931751394603*g3^2*p2 - 1.8163349832174116*g1^2*p1*p2 - 0.9167144057621401*g1*g2*p1*p2 + 1.0203368504488892*g2^2*p1*p2 - 0.23194646823111892*g1*g3*p1*p2 + 0.539670777307627*g2*g3*p1*p2 + 0.7959981327685224*g3^2*p1*p2 + 0.08717268867521591*g1^2*p2^2 + 0.9504154644263471*g1*g2*p2^2 - 0.48206756571420756*g2^2*p2^2 - 1.065062423127697*g1*g3*p2^2 + 0.1209952909274163*g2*g3*p2^2 + 0.3948948770389917*g3^2*p2^2 + 0.289766299873838*g1^2*p3 - 1.2778927965251532*g1*g2*p3 + 0.9087896778886251*g2^2*p3 - 0.5812612591154215*g1*g3*p3 - 0.7595904624983555*g2*g3*p3 + 0.5084892760496751*g3^2*p3 - 0.3268802641947883*g1^2*p1*p3 + 0.657630238424344*g1*g2*p1*p3 + 1.1093919363972093*g2^2*p1*p3 + 0.4551393419480071*g1*g3*p1*p3 + 1.8553852513069364*g2*g3*p1*p3 - 0.7825116722024211*g3^2*p1*p3 + 0.5810469298461638*g1^2*p2*p3 - 1.1557382363783264*g1*g2*p2*p3 - 0.11367961187637783*g2^2*p2*p3 + 1.7077140933509898*g1*g3*p2*p3 - 0.36547942767108677*g2*g3*p2*p3 - 0.4673673179697859*g3^2*p2*p3 + 0.12034280775306151*g1^2*p3^2 - 0.5733857796356615*g1*g2*p3^2 + 0.6489566339058018*g2^2*p3^2 + 0.2663669912953945*g1*g3*p3^2 - 0.9878214357030672*g2*g3*p3^2 - 0.7692994416588633*g3^2*p3^2 - 0.19466145678474384*q0 - 0.7153041427190404*p1*q0 - 1.3528776260043915*p2*q0 - 1.7070452538121381*p3*q0 - 1.0516635822669562*q1 + 1.2244185478631853*p1*q1 - 0.05844567698552443*p2*q1 - 0.37706149953585283*p3*q1 + 0.580102254517945*q2 + 1.2898860704586343*p1*q2 - 0.6655948497180294*p2*q2 + 0.697758704890495*p3*q2 - 0.042921436747585445*q3 + 0.5172073855756967*p1*q3 + 0.6917094054122289*p2*q3 - 1.4579672250860476*p3*q3,

    0.16011034303688113*g1^2 - 0.9005468824403076*g1*g2 - 0.3519015838689263*g2^2 + 0.5202586158306898*g1*g3 + 0.908682123022068*g2*g3 - 0.4464562170645777*g3^2 - 0.13844524415679324*g1^2*p1 + 1.5568085644333742*g1*g2*p1 + 1.6863862382239232*g2^2*p1 - 1.7409458121154344*g1*g3*p1 - 0.13872356093602894*g2*g3*p1 - 0.5159047084859331*g3^2*p1 - 0.2741643484200128*g1^2*p1^2 - 0.34212012775550327*g1*g2*p1^2 - 0.07542436599114127*g2^2*p1^2 + 0.37458987278720324*g1*g3*p1^2 + 0.4782561996467687*g2*g3*p1^2 + 0.3495887144111541*g3^2*p1^2 + 0.41377445473869573*g1^2*p2 - 1.5789383736211624*g1*g2*p2 + 1.268319517294935*g2^2*p2 + 0.6163793667190677*g1*g3*p2 - 0.43374574206406646*g2*g3*p2 - 0.2061458017243186*g3^2*p2 + 0.14555549639831628*g1^2*p1*p2 - 1.1674745895517964*g1*g2*p1*p2 - 0.9428064489876502*g2^2*p1*p2 + 0.0024916775818734295*g1*g3*p1*p2 + 0.5291621555283466*g2*g3*p1*p2 + 0.7972509525893339*g3^2*p1*p2 + 0.1807885464109201*g1^2*p2^2 + 0.9404541869824675*g1*g2*p2^2 - 0.5780030515551372*g2^2*p2^2 - 1.0257418447585547*g1*g3*p2^2 + 0.09251778173989315*g2*g3*p2^2 + 0.39721450514421713*g3^2*p2^2 + 0.40272988912109214*g1^2*p3 - 0.8272484673958682*g1*g2*p3 + 1.057139636924469*g2^2*p3 - 0.12353226665002319*g1*g3*p3 - 2.5741855761862396*g2*g3*p3 + 1.560474007685759*g3^2*p3 - 0.6150996832616941*g1^2*p1*p3 + 0.09937192239106099*g1*g2*p1*p3 + 0.8226042775491553*g2^2*p1*p3 + 0.4732438203631739*g1*g3*p1*p3 + 1.6946050580334677*g2*g3*p1*p3 - 0.20750459428746135*g3^2*p1*p3 + 0.7556396990592089*g1^2*p2*p3 - 1.412614951501404*g1*g2*p2*p3 - 0.09676545515565128*g2^2*p2*p3 - 0.9781576342585658*g1*g3*p2*p3 + 2.5006617995144724*g2*g3*p2*p3 - 0.6588742439035575*g3^2*p2*p3 + 0.09337580200909272*g1^2*p3^2 - 0.5983340592269643*g1*g2*p3^2 + 0.6534274175462785*g2^2*p3^2 + 0.6511519719713513*g1*g3*p3^2 - 0.5707739813866619*g2*g3*p3^2 - 0.7468032195553712*g3^2*p3^2 + 0.6382474578966228*q0 - 1.032036285581197*p1*q0 - 1.4759481703093122*p2*q0 - 3.0203435337313205*p3*q0 - 0.22812438675350769*q1 - 0.2157590670168509*p1*q1 - 0.1270558344695696*p2*q1 - 0.5148593639524484*p3*q1 + 1.3667793800860086*q2 - 0.06171123442924746*p1*q2 - 0.7314954155886625*p2*q2 + 0.7189348075213543*p3*q2 - 0.902118536026858*q3 + 0.43214823742186254*p1*q3 + 0.6677624868260497*p2*q3 + 0.5162571144422815*p3*q3,

    0.20816475809219404*g1^2 - 0.44624795696445435*g1*g2 + 0.1573457781818856*g2^2 + 0.2432511536576595*g1*g3 - 0.4587424991969163*g2*g3 - 0.048748564896809544*g3^2 - 0.542424299098038*g1^2*p1 + 0.412606879197033*g1*g2*p1 + 0.026228493490255755*g2^2*p1 - 0.7135454436169615*g1*g3*p1 + 0.38566725023570736*g2*g3*p1 + 0.06641935708182738*g3^2*p1 + 0.1946156279601214*g1^2*p1^2 + 1.0726514255671113*g1*g2*p1^2 - 0.29746151974577967*g2^2*p1^2 + 0.3474150051655493*g1*g3*p1^2 + 1.5662794253637933*g2*g3*p1^2 + 0.10284589178565828*g3^2*p1^2 + 0.8655281158446179*g1^2*p2 - 1.4227007533612923*g1*g2*p2  - 0.3561608986253729*g2^2*p2 - 0.7002053827479838*g1*g3*p2 + 0.1451263721376322*g2*g3*p2 - 0.4446427929457582*g3^2*p2  - 0.35089234105147404*g1^2*p1*p2 - 1.8637325747105546*g1*g2*p1*p2 - 0.7643910878410862*g2^2*p1*p2 + 0.7421389633104346*g1*g3*p1*p2 + 0.8043890896223826*g2*g3*p1*p2 + 1.1152834288925604*g3^2*p1*p2 - 0.05058692105297476*g1^2*p2^2 - 0.8545531093164939*g1*g2*p2^2 - 0.25045809562785276*g2^2*p2^2 - 1.482438556873845*g1*g3*p2^2 - 0.2760311985894717*g2*g3*p2^2 + 0.30104501668082756*g3^2*p2^2 + 0.41615153726461007*g1^2*p3 - 1.6031132124173149*g1*g2*p3 + 1.1652768530802575*g2^2*p3 + 0.1236694347662175*g1*g3*p3  - 0.033510271732486586*g2*g3*p3 + 0.6625023868605743*g3^2*p3 - 0.06941899872446193*g1^2*p1*p3 - 0.5612725019588681*g1*g2*p1*p3 + 1.4835363108262836*g2^2*p1*p3 - 0.8310204341509994*g1*g3*p1*p3 + 1.3650887611787323*g2*g3*p1*p3 - 1.4141173121018216*g3^2*p1*p3 - 0.2915853970368523*g1^2*p2*p3 - 1.2521117933146961*g1*g2*p2*p3 + 0.38706376702247*g2^2*p2*p3 + 1.2309129178715645*g1*g3*p2*p3 + 2.001338697637118*g2*g3*p2*p3 - 0.09547836998561768*g3^2*p2*p3 - 0.14402870690714664*g1^2*p3^2 - 0.2180983162506176*g1*g2*p3^2 + 0.5479196153736324*g2^2*p3^2 + 1.1350235517082958*g1*g3*p3^2 - 1.2902482267743214*g2*g3*p3^2 - 0.40389090846648584*g3^2*p3^2 - 0.31676197137727014*q0 + 0.44977644852595483*p1*q0 - 0.06472442427348668*p2*q0 - 2.2439307772054415*p3*q0 - 0.5547165223690258*q1 + 0.23831878651082344*p1*q1 + 0.031977776730485255*p2*q1 + 0.16687455406564522*p3*q1 + 0.9423377906275198*q2 + 1.376589178886685*p1*q2 + 0.5306523901876015*p2*q2 + 0.4754891181933043*p3*q2 + 0.09673230093655334*q3 + 0.0892904130224598*p1*q3 + 0.943551163213123*p2*q3 - 1.2527250130712726*p3*q3,
     -0.04095049824628835*g1^2 + 0.043116025511842154*g1*g2 + 0.003940499198786224*g2^2 + 0.7629770334036455*g1*g3 - 0.8492350760146794*g2*g3 - 0.2784174783424625*g3^2 + 0.4281444292173086*g1^2*p1 - 0.8828960936117035*g1*g2*p1 + 0.5676682886279524*g2^2*p1 - 1.3924527881029736*g1*g3*p1 - 0.08671339002537767*g2*g3*p1 - 0.7256104095984146*g3^2*p1 - 0.5422313181564682*g1^2*p1^2 + 0.4871946471731439*g1*g2*p1^2 + 0.33193585698170985*g2^2*p1^2 - 0.10273772900088107*g1*g3*p1^2 + 0.6624874115365778*g2*g3*p1^2 + 0.21029546117475836*g3^2*p1^2 + 1.1356911567255628*g1^2*p2 - 1.222100685178249*g1*g2*p2 - 0.09334002143332033*g2^2*p2 - 1.6524959396527132*g1*g3*p2 + 1.5569725124184146*g2*g3*p2 + 0.7669386068453008*g3^2*p2 - 1.5587675264538823*g1^2*p1*p2 - 1.1317581527003464*g1*g2*p1*p2 + 0.01829080736739283*g2^2*p1*p2 + 1.0020439814840232*g1*g3*p1*p2 + 0.7984049760283556*g2*g3*p1*p2 + 1.5404767190864894*g3^2*p1*p2 + 0.4331904808414006*g1^2*p2^2 - 0.1389285705830233*g1*g2*p2^2 - 0.012183710127155694*g2^2*p2^2 - 0.23769385136664706*g1*g3*p2^2 + 0.24382837068218804*g2*g3*p2^2 - 0.4210067707142449*g3^2*p2^2 + 0.21198555697618326*g1^2*p3 - 0.01647709154129892*g1*g2*p3 - 0.18945218115272705*g2^2*p3 + 0.77528301920843*g1*g3*p3 - 2.0514046696465*g2*g3*p3 + 0.10432028494512646*g3^2*p3 - 0.1676311045685404*g1^2*p1*p3 - 1.325546381572095*g1*g2*p1*p3 + 0.8078478514339609*g2^2*p1*p3 - 0.9669176278885212*g1*g3*p1*p3 - 0.36872926793739896*g2*g3*p1*p3 - 0.6402167468654205*g3^2*p1*p3 + 0.6297207100844667*g1^2*p2*p3 - 2.067771321161895*g1*g2*p2*p3  - 0.17862819697751522*g2^2*p2*p3 + 0.020391323549034297*g1*g3*p2*p3  + 2.777563669744398*g2*g3*p2*p3 - 0.45109251310695153*g3^2*p2*p3 + 0.10904083731506761*g1^2*p3^2 - 0.3482660765901206*g1*g2*p3^2 - 0.31975214685455416*g2^2*p3^2 + 0.34043158036752813*g1*g3*p3^2 - 0.9063157822187657*g2*g3*p3^2 + 0.21071130953948652*g3^2*p3^2 + 0.31542747738996463*q0 - 0.27020230824684643*p1*q0 - 1.8092897421375431*p2*q0 - 0.12685366076858268*p3*q0 - 0.9721054060313574*q1 + 1.4332583965298273*p1*q1 - 0.3658292969614953*p2*q1 - 0.458292808629767*p3*q1 + 1.1057480001700448*q2 + 0.649216154064302*p1*q2 + 0.1435470147844548*p2*q2 + 1.8049686045262234*p3*q2 + 0.3619641675513017*q3 + 1.0386298649000567*p1*q3 - 0.2739870731830222*p2*q3 - 0.38992289294835114*p3*q3,
     -0.41615764608945516*g1^2 - 1.2331171001793817*g1*g2 + 0.10423594498637195*g2^2 + 0.4451741240918564*g1*g3 - 0.0807794759847403*g2*g3 + 0.015584822151867354*g3^2 + 0.5169791211840113*g1^2*p1 - 0.3281633186673521*g1*g2*p1 + 0.10768082059655043*g2^2*p1 - 1.78387184821123*g1*g3*p1 + 0.1962385955438586*g2*g3*p1 - 0.0932755727182936*g3^2*p1 - 0.12307658314371513*g1^2*p1^2 + 1.3649915585405705*g1*g2*p1^2 - 0.01641144275933561*g2^2*p1^2 + 1.282522294958988*g1*g3*p1^2 + 0.13118389677242223*g2*g3*p1^2 + 0.13948802590305073*g3^2*p1^2 - 0.4784260518169776*g1^2*p2 - 2.092134198423298*g1*g2*p2 - 0.2652478875380973*g2^2*p2 - 0.9322070346912057*g1*g3*p2 + 0.28229645793462466*g2*g3*p2 - 0.07438003692790207*g3^2*p2 + 0.1376969425780227*g1^2*p1*p2 - 1.3462021315216954*g1*g2*p1*p2 - 0.35831006800801096*g2^2*p1*p2 - 0.3236078908735904*g1*g3*p1*p2 - 0.23312973570099904*g2*g3*p1*p2 + 0.2206131254299883*g3^2*p1*p2 + 0.03212463900566726*g1^2*p2^2 - 0.543299069419884*g1*g2*p2^2 - 0.1081313237618179*g2^2*p2^2 - 1.8348908578280814*g1*g3*p2^2 - 0.3379484782876818*g2*g3*p2^2 + 0.07600668475615065*g3^2*p2^2 + 0.6002572265406737*g1^2*p3 + 0.1802055521069689*g1*g2*p3 + 1.5644331744196656*g2^2*p3 - 0.5267676534246675*g1*g3*p3 - 1.676157338774887*g2*g3*p3 + 0.4008307913692461*g3^2*p3 - 0.5875157312206984*g1^2*p1*p3 - 0.26794777918572443*g1*g2*p1*p3 + 1.8161766879761405*g2^2*p1*p3 + 0.08217873131336825*g1*g3*p1*p3 + 0.41625228622759664*g2*g3*p1*p3 - 1.228660956755442*g3^2*p1*p3 + 0.801736157469905*g1^2*p2*p3 - 0.06119905069237429*g1*g2*p2*p3 + 0.5320140756032581*g2^2*p2*p3 - 0.5120321678148483*g1*g3*p2*p3 + 1.3994959245799465*g2*g3*p2*p3 - 1.3337502330731632*g3^2*p2*p3 + 0.09095194413804787*g1^2*p3^2 - 0.8216924891206866*g1*g2*p3^2 + 0.12454276652115351*g2^2*p3^2 + 0.5523685628690934*g1*g3*p3^2 + 0.20676458151525956*g2*g3*p3^2 - 0.2154947106592014*g3^2*p3^2 + 0.29633687895121585*q0 - 0.5313843690622682*p1*q0 + 0.8180539762829769*p2*q0 - 2.5655211923295855*p3*q0 - 1.2602089223582702*q1 + 1.1609512634985952*p1*q1 - 1.5289512288758575*p2*q1 - 0.17424561505966216*p3*q1 + 0.10364901876603111*q2 - 0.013973408764994696*p1*q2 - 0.30080374272031296*p2*q2 + 1.5473286276087392*p3*q2 - 0.04276145686338927*q3 + 0.12497482938060817*p1*q3 + 0.06337974500071619*p2*q3 - 1.1223229062458282*p3*q3,

    0.0995239208560676*g1^2 - 0.4069835646127759*g1*g2 + 0.2822661867004282*g2^2 + 0.8451107786883508*g1*g3 - 0.5316708978452792*g2*g3 - 0.8798108963881374*g3^2 - 0.08841286657967666*g1^2*p1 - 0.5437440767827946*g1*g2*p1 + 0.8310187547943032*g2^2*p1 - 2.300828301805621*g1*g3*p1 + 1.7576212044508612*g2*g3*p1 + 0.8123567734335738*g3^2*p1 - 0.06820768665581543*g1^2*p1^2 + 1.3043060864108098*g1*g2*p1^2 - 0.009054065062096384*g2^2*p1^2 + 1.3077574118576432*g1*g3*p1^2 + 0.06768474224643237*g2*g3*p1^2 + 0.07726175171791182*g3^2*p1^2 + 0.21964796928975452*g1^2*p2 - 0.7213869971382741*g1*g2*p2 + 0.06428880062023586*g2^2*p2 + 0.6769693907722261*g1*g3*p2 + 1.0938889186830445*g2*g3*p2 + 1.033778176369894*g3^2*p2 + 0.04579106142198317*g1^2*p1*p2 - 1.6136177673887118*g1*g2*p1*p2 - 0.08313705215573315*g2^2*p1*p2 - 0.3180870028357637*g1*g3*p1*p2 - 0.1323253404673372*g2*g3*p1*p2 + 0.03734599073374998*g3^2*p1*p2 + 0.06270407236710761*g1^2*p2^2 - 0.19162792586107563*g1*g2*p2^2 - 0.009068232422623026*g2^2*p2^2 - 1.9648502282295741*g1*g3*p2^2 - 0.09696983966017489*g2*g3*p2^2 - 0.053635839944484585*g3^2*p2^2 + 0.09136896607569936*g1^2*p3 - 1.5831271499881143*g1*g2*p3  + 1.1586203238882382*g2^2*p3 + 1.561944037084699*g1*g3*p3 + 0.5373765622781439*g2*g3*p3 - 0.7510078085771675*g3^2*p3 + 0.03313601431114031*g1^2*p1*p3 + 0.999575794109116*g1*g2*p1*p3 + 0.0005673645755714211*g2^2*p1*p3 - 1.8768450822200746*g1*g3*p1*p3 + 0.15574381848783894*g2*g3*p1*p3 - 0.03370337888671173*g3^2*p1*p3 + 0.040671941002624216*g1^2*p2*p3 - 2.5271478321380254*g1*g2*p2*p3 - 0.11221443152128503*g2^2*p2*p3 + 0.5426249873354888*g1*g3*p2*p3 + 0.01844664972668835*g2*g3*p2*p3 + 0.07154249051866082*g3^2*p2*p3 + 0.005503614288707821*g1^2*p3^2 - 1.1126781605497342*g1*g2*p3^2 + 0.01812229748471941*g2^2*p3^2 + 0.657092816371931*g1*g3*p3^2 + 0.029285097413742515*g2*g3*p3^2 - 0.02362591177342723*g3^2*p3^2 + 0.49802078883164164*q0 - 1.5549626616482004*p1*q0 - 1.3177149462798845*p2*q0 - 0.49898148138677*p3*q0 - 0.12338824573791841*q1 - 0.06111510706623446*p1*q1 - 0.04157014428073322*p2*q1 - 0.007996832327295606*p3*q1 + 0.3953545453627241*q2 + 1.168253359907005*p1*q2 + 0.1289788318383644*p2*q2 + 1.6166074402033737*p3*q2 - 1.1634376524391858*q3 + 1.1753697910765462*p1*q3 + 1.303746265623441*p2*q3 - 0.9550529463247988*p3*q3,
     -0.43946259392041137*g1^2 - 0.2794194400312886*g1*g2 + 0.28633513817241923*g2^2 + 0.4837842932289669*g1*g3 - 0.2789447622513483*g2*g3 + 0.00840391192868436*g3^2 - 1.5744202005495247*g1^2*p1 - 0.2861663506179279*g1*g2*p1 + 0.7386941869051702*g2^2*p1 + 0.5331769017373268*g1*g3*p1 + 0.3754313238644043*g2*g3*p1 - 0.49557272342979514*g3^2*p1 - 0.1697121445535322*g1^2*p1^2 + 0.17649121044520624*g1*g2*p1^2 + 0.230772384999563*g2^2*p1^2 + 1.5487222154552227*g1*g3*p1^2 + 1.1549190062177301*g2*g3*p1^2 - 0.061060240446030825*g3^2*p1^2 + 0.3428887153466959*g1^2*p2 - 2.0520172954900926*g1*g2*p2 + 0.2548434377395737*g2^2*p2 - 1.1661156986017325*g1*g3*p2 + 1.1418016311188108*g2*g3*p2 - 0.03717326023724154*g3^2*p2 - 1.1298456541690676*g1^2*p1*p2 - 2.0395244570443807*g1*g2*p1*p2 - 0.46158802071478294*g2^2*p1*p2 - 0.7937227851031279*g1*g3*p1*p2 + 1.3925863980743391*g2*g3*p1*p2 + 1.5914336748838505*g3^2*p1*p2 + 0.5461657368466324*g1^2*p2^2 - 0.40076456254394*g1*g2*p2^2 - 0.5789448963236655*g2^2*p2^2 - 1.2606913391501393*g1*g3*p2^2 - 0.8839870841626388*g2*g3*p2^2 + 0.03277915947703309*g3^2*p2^2 - 0.8186969313136523*g1^2*p3 - 0.02925076747851844*g1*g2*p3 + 0.7392135200738452*g2^2*p3 + 0.12848993732517067*g1*g3*p3 - 0.9246721961161062*g2*g3*p3 + 0.34711806557393554*g3^2*p3 - 1.3114980542204153*g1^2*p1*p3 + 0.7458693987541652*g1*g2*p1*p3  + 1.309443056027262*g2^2*p1*p3 - 0.5609050703954811*g1*g3*p1*p3 - 0.560863866588139*g2*g3*p1*p3 + 0.0020549981931533984*g3^2*p1*p3 + 0.21040044489655793*g1^2*p2*p3 - 2.094674294769259*g1*g2*p2*p3 + 0.9014935747766074*g2^2*p2*p3 - 0.7741828342315165*g1*g3*p2*p3 + 1.0179705141740856*g2*g3*p2*p3 - 1.1118940196731655*g3^2*p2*p3 - 0.37645359229310016*g1^2*p3^2 + 0.22427335209873375*g1*g2*p3^2 + 0.34817251132410243*g2^2*p3^2 - 0.2880308763050832*g1*g3*p3^2 - 0.27093192205509137*g2*g3*p3^2 + 0.028281080968997733*g3^2*p3^2 + 0.14472354381930777*q0 + 1.3312987370741496*p1*q0 - 0.560558892849028*p2*q0 - 0.26763465433412853*p3*q0 + 0.34474968158768826*q1 + 1.1967487089081013*p1*q1 - 0.5416788729536276*p2*q1 + 0.35732161373948407*p3*q1 + 0.4086834090152563*q2 + 0.9057045158335053*p1*q2 + 0.790785041402257*p2*q2 + 0.25330644272284614*p3*q2 - 0.3854162586341408*q3 - 0.04738889206888738*p1*q3 + 1.2360880142692436*p2*q3 - 0.03146357181747231*p3*q3,
     -0.3364985136694329*g1^2 - 0.6220713377668979*g1*g2  + 0.30253663989085705*g2^2 + 0.09691241616903694*g1*g3 - 0.23475287081412524*g2*g3 + 0.03759204880218851*g3^2 - 1.0655428995245861*g1^2*p1 + 1.2015525036653503*g1*g2*p1 + 0.7601103825805371*g2^2*p1 - 0.19738187400246698*g1*g3*p1 + 0.6151569823342218*g2*g3*p1 - 0.19896444557080828*g3^2*p1 + 0.12867336914719074*g1^2*p1^2 + 1.2422164684570598*g1*g2*p1^2 - 0.32698199621342966*g2^2*p1^2 + 0.988903111633841*g1*g3*p1^2 + 0.015998020944011523*g2*g3*p1^2 + 0.19830862706623895*g3^2*p1^2 - 0.36506543203523656*g1^2*p2 - 1.8191248166962*g1*g2*p2 - 0.04702353062492582*g2^2*p2 + 0.22650992191881464*g1*g3*p2 + 0.6757394159618821*g2*g3*p2 - 0.20784006514510822*g3^2*p2 - 1.3823450011564475*g1^2*p1*p2 - 0.10921823662847378*g1*g2*p1*p2 + 0.9411204746791397*g2^2*p1*p2 + 0.6416176172245012*g1*g3*p1*p2 + 0.6050178988297725*g2*g3*p1*p2 + 0.4412245264773078*g3^2*p1*p2 + 0.06735301456610396*g1^2*p2^2 - 0.8585408362609629*g1*g2*p2^2 - 0.3090762258515019*g2^2*p2^2 - 0.4001418483728859*g1*g3*p2^2 + 0.4779675484647562*g2*g3*p2^2 + 0.24172321128539795*g3^2*p2^2 - 0.5137007261636957*g1^2*p3 - 0.9820327320175907*g1*g2*p3 + 0.926095293449879*g2^2*p3 - 0.7003487284409724*g1*g3*p3 - 0.9430374895060032*g2*g3*p3 + 0.08868217302851986*g3^2*p3 - 0.8028419556001375*g1^2*p1*p3 + 1.4979362188623169*g1*g2*p1*p3 + 0.6798639272820932*g2^2*p1*p3 - 1.0625669164194511*g1*g3*p1*p3  + 1.6499737530190743*g2*g3*p1*p3 + 0.12297802831804432*g3^2*p1*p3 - 0.2829951083769471*g1^2*p2*p3 - 1.67180387717421*g1*g2*p2*p3 + 0.22862840829725906*g2^2*p2*p3 - 0.45573353517584536*g1*g3*p2*p3 - 0.030399560249614144*g2*g3*p2*p3 + 0.05436670007968804*g3^2*p2*p3 - 0.1960263837132947*g1^2*p3^2 - 0.3836756321960968*g1*g2*p3^2 + 0.6360582220649316*g2^2*p3^2 - 0.5887612632609552*g1*g3*p3^2 - 0.49396556940876774*g2*g3*p3^2 - 0.4400318383516369*g3^2*p3^2 - 0.0036301750236126602*q0 + 0.5043969625148573*p1*q0 + 0.6199290278052706*p2*q0 - 0.5010767403147031*p3*q0 + 0.27845755063810307*q1 + 0.9141858220671707*p1*q1 - 0.04474642022970836*p2*q1 + 0.21510107355697375*p3*q1 + 0.6261152163591358*q2 - 0.23441023070199246*p1*q2  + 0.5860551333365724*p2*q2 + 0.6514610848268766*p3*q2 - 0.1414398428437503*q3 + 0.18837597417330545*p1*q3 + 0.23529746996939166*p2*q3 + 0.34502969365002717*p3*q3,
     -0.4997153826800627*g1^2 + 1.1379854833548109*g1*g2 - 0.6474309248194395*g2^2 + 0.5441177673182162*g1*g3  - 0.6356209362624222*g2*g3 - 0.00249167421336129*g3^2 - 1.5544041627556424*g1^2*p1 + 2.3739950333676623*g1*g2*p1  - 0.626498203472872*g2^2*p1 - 1.251167893581656*g1*g3*p1 + 1.014528145879352*g2*g3*p1 - 0.07360314669929015*g3^2*p1 - 0.11700211538167525*g1^2*p1^2 + 1.4292168664824174*g1*g2*p1^2 - 0.08561887506760898*g2^2*p1^2 + 0.11756149787203274*g1*g3*p1^2 + 1.2753331275027957*g2*g3*p1^2 + 0.20262099044928422*g3^2*p1^2 - 0.5322121916934887*g1^2*p2 - 0.5705836762208184*g1*g2*p2 + 1.338552997747994*g2^2*p2 + 0.1626686039427204*g1*g3*p2 + 0.8255938648528275*g2*g3*p2 - 0.2877403739032339*g3^2*p2 - 1.7484788356484942*g1^2*p1*p2 - 0.7462669895235802*g1*g2*p1*p2 + 0.9024930142580471*g2^2*p1*p2 - 0.836428661723292*g1*g3*p1*p2 + 0.10315211869079574*g2*g3*p1*p2 + 0.8459858213904471*g3^2*p1*p2 + 0.0524174914208927*g1^2*p2^2 - 1.0349080883999526*g1*g2*p2^2 - 0.4456138068634661*g2^2*p2^2 - 0.8702119648626987*g1*g3*p2^2 + 0.1347633798985455*g2*g3*p2^2 + 0.3931963154425734*g3^2*p2^2 + 0.7263228121576023*g1^2*p3 - 2.2210658961930227*g1*g2*p3 + 1.6108892968930464*g2^2*p3 - 1.0504867050871953*g1*g3*p3 + 1.0659099359741329*g2*g3*p3 + 0.5672833993063003*g3^2*p3 - 0.06811413086649659*g1^2*p1*p3 - 0.31001429706461026*g1*g2*p1*p3 + 1.4153449120254278*g2^2*p1*p3 - 1.8526453839326555*g1*g3*p1*p3 + 1.3321822393355465*g2*g3*p1*p3 - 1.3472307811589312*g3^2*p1*p3 + 0.8918357940981352*g1^2*p2*p3 - 1.500420306224534*g1*g2*p2*p3 - 0.43633002867072346*g2^2*p2*p3 - 1.087845994748532*g1*g3*p2*p3 - 0.15102803926606495*g2*g3*p2*p3 - 0.4555057654274117*g3^2*p2*p3 + 0.06458462396078254*g1^2*p3^2 - 0.39430877808246473*g1*g2*p3^2 + 0.5312326819310751*g2^2*p3^2 + 0.752650466990666*g1*g3*p3^2 - 1.4100965074013412*g2*g3*p3^2 - 0.5958173058918577*g3^2*p3^2 + 1.1496379817128635*q0 + 2.2545055129278047*p1*q0 - 0.5186004321512713*p2*q0 - 2.904495508356949*p3*q0 + 0.3811371457709*q1 + 1.1561384774363503*p1*q1 - 0.03459051968276095*p2*q1 - 0.5860104174434722*p3*q1 - 0.42262305475993184*q2 - 0.06960283700014054*p1*q2 + 0.6975228406191866*p2*q2 + 1.1761445768437886*p3*q2 - 0.4132504788764938*q3 + 1.0477095469315294*p1*q3 + 0.5581768684334514*p2*q3 + 0.4360576481532085*p3*q3,

    -1. + 0.9336143308746049*g1 + 1.1781580271766483*g2 + 0.551650235747964*g3
    ]
    #

    return TestSystem(equations, multi_bezout_number=(18432, [[g1, g2, g3], [p1, p2, p3, q0, q1, q2, q3]]), mixed_volume=1024)
end


text:
# # Time-Varying Bayesian Linear Regression

using AbstractGPs
using ColorTypes
using FixedPointNumbers
using Plots
using Random
using Stheno

# ## Define and inspect our model

#=
g1 and g2 are time-varying basis functions. In a real application, these might be specified
by data.

w1 and w2 are the time-varying basis functions for these bases, specified by slowly
varying GPs.

f is the prediction of the regressor.

y is the addition of f and rough temporally-correlated "noise".
=#
f = @gppp let
    g1 = x->x / 4
    g2 = cos
    w1 = stretch(GP(SEKernel()), 0.2)
    w2 = stretch(GP(SEKernel()), 1)
    f = g1 * w1 + g2 * w2
    y = f + 0.3 * GP(Matern12Kernel())
end;

# Sample from the prior from plotting and for conditioning.
rng, N, Nplot, S = MersenneTwister(123456), 250, 500, 100;
x = GPPPInput(:y, sort(rand(rng, N) * 10));
ŷ = rand(rng, f(x));

# Compute the posterior.
f′ = posterior(f(x), ŷ);

# Sample from the posterior and write to file.
xp_ = range(-2.5, stop=12.5, length=Nplot);
xp_w1 = GPPPInput(:w1, xp_);
xp_w2 = GPPPInput(:w2, xp_);
xp_y = GPPPInput(:y, xp_);

# Generate joint posterior samples.
xp = BlockData(xp_w1, xp_w2, xp_y);
f′_xp = rand(rng, f′(xp, 1e-9), S);
w1′s, w2′s, y′s = split(xp, f′_xp);

# ## Plot results

gr();
posterior_plot = plot(
    legend=:topleft,
    legendfont=Plots.Font(
        "sans-serif",
        10,
        :hcenter,
        :vcenter,
        0.0,
        RGB{Normed{UInt8, 8}}(0.0,0.0,0.0)
    ),
    background_color_legend=RGBA(1, 1, 1, 0),
    foreground_color_legend=RGBA(1, 1, 1, 0),
);

# Plot posterior over w1.
plot!(posterior_plot, xp_, f′(xp_w1); color=:green, label="w1");
plot!(posterior_plot, xp_, w1′s; color=:green, label="", linewidth=1, alpha=0.2);

# Plot posterior over w2.
plot!(posterior_plot, xp_, f′(xp_w2); color=:magenta, label="w2");
plot!(posterior_plot, xp_, w2′s; color=:magenta, label="", linewidth=1, alpha=0.2);

# Plot x1 and x2
plot!(posterior_plot, xp_, (x->x / 4).(xp_);
    linecolor=:black,
    linewidth=1.0,
    label="x / 4",
);
plot!(posterior_plot, xp_, cos.(xp_);
    linecolor=:black,
    linewidth=1.0,
    linestyle=:dash,
    label="cos",
);

# Plot samples against which we're regressing.
scatter!(posterior_plot, x.x, ŷ;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

# Plot posterior over `y`.
plot!(posterior_plot, xp_, f′(xp_y); color=:blue, label="y");
plot!(posterior_plot, xp_, y′s; color=:blue, label="", linewidth=1, alpha=0.2);

posterior_plot


text:
wider(::Type{Int64}) = Int128
wider(::Type{Int128}) = BigInt

# include a non-inlined version in case of widening (otherwise, all widened cases would fully inline)
@noinline _typeparser(::Type{T}, source, pos, len, b, code, options::Options{ignorerepeated, Q, debug, S, D, DF}, ::Type{IntType}) where {T <: AbstractFloat, ignorerepeated, Q, debug, S, D, DF, IntType} =
    typeparser(T, source, pos, len, b, code, options, IntType)

@inline function typeparser(::Type{T}, source, pos, len, b, code, options::Options{ignorerepeated, Q, debug, S, D, DF}, ::Type{IntType}=Int64) where {T <: AbstractFloat, ignorerepeated, Q, debug, S, D, DF, IntType}
    startpos = pos
    origb = b
    x = zero(T)
    digits = zero(IntType)
    if debug
        println("float parsing")
    end
    neg = b == UInt8('-')
    if neg || b == UInt8('+')
        pos += 1
        incr!(source)
    end
    if eof(source, pos, len)
        code |= INVALID | EOF
        @goto done
    end
    b = peekbyte(source, pos)
    if b == options.decimal
        @goto decimalcheck
    end
    b -= UInt8('0')
    if debug
        println("float 1) $(b + UInt8('0'))")
    end
    if b > 0x09
        # character isn't a digit, check for special values, otherwise INVALID
        b += UInt8('0')
        if b == UInt8('n') || b == UInt8('N')
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto invalid
            end
            b = peekbyte(source, pos)
            if b == UInt8('a') || b == UInt8('A')
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto invalid
                end
                b = peekbyte(source, pos)
                if b == UInt8('n') || b == UInt8('N')
                    x = T(NaN)
                    pos += 1
                    incr!(source)
                    if eof(source, pos, len)
                        code |= EOF
                    end
                    code |= OK
                    @goto done
                end
            end
        elseif b == UInt8('i') || b == UInt8('I')
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto invalid
            end
            b = peekbyte(source, pos)
            if b == UInt8('n') || b == UInt8('N')
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto invalid
                end
                b = peekbyte(source, pos)
                if b == UInt8('f') || b == UInt8('F')
                    x = ifelse(neg, T(-Inf), T(Inf))
                    code |= OK
                    pos += 1
                    incr!(source)
                    if eof(source, pos, len)
                        code |= EOF
                        @goto done
                    end
                    b = peekbyte(source, pos)
                    if b == UInt8('i') || b == UInt8('I')
                        pos += 1
                        incr!(source)
                        if eof(source, pos, len)
                            code |= EOF
                            @goto done
                        end
                        b = peekbyte(source, pos)
                        if b == UInt8('n') || b == UInt8('N')
                            pos += 1
                            incr!(source)
                            if eof(source, pos, len)
                                code |= EOF
                                @goto done
                            end
                            b = peekbyte(source, pos)
                            if b == UInt8('i') || b == UInt8('I')
                                pos += 1
                                incr!(source)
                                if eof(source, pos, len)
                                    code |= EOF
                                    @goto done
                                end
                                b = peekbyte(source, pos)
                                if b == UInt8('t') || b == UInt8('T')
                                    pos += 1
                                    incr!(source)
                                    if eof(source, pos, len)
                                        code |= EOF
                                        @goto done
                                    end
                                    b = peekbyte(source, pos)
                                    if b == UInt8('y') || b == UInt8('Y')
                                        pos += 1
                                        incr!(source)
                                        if eof(source, pos, len)
                                            code |= EOF
                                        end
                                        @goto done
                                    end
                                end
                            end
                        end
                    end
                    @goto done
                end
            end
        else
        end
@label invalid
        fastseek!(source, startpos)
        code |= INVALID
        @goto done
    end
    while true
        digits = IntType(10) * digits + b
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            x = T(ifelse(neg, -digits, digits))
            code |= OK | EOF
            @goto done
        end
        b = peekbyte(source, pos) - UInt8('0')
        if debug
            println("float 2) $(b + UInt8('0'))")
        end
        b > 0x09 && break
        if overflows(IntType) && digits > overflowval(IntType)
            fastseek!(source, startpos)
            return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
        end
    end
    b += UInt8('0')
    if debug
        println("float 3) $(Char(b))")
    end
@label decimalcheck
    if b == options.decimal
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            x = T(ifelse(neg, -digits, digits))
            code |= OK | EOF
            @goto done
        end
        b = peekbyte(source, pos)
        if b - UInt8('0') > 0x09 && !(b == UInt8('e') || b == UInt8('E') || b == UInt8('f') || b == UInt8('F'))
            x = T(ifelse(neg, -digits, digits))
            code |= OK
            @goto done
        end
    end
    frac = 0
    if debug
        println("float 4) $(Char(b))")
    end
    b -= UInt8('0')
    if b < 0x0a
        while true
            digits = IntType(10) * digits + b
            pos += 1
            incr!(source)
            frac += 1
            if eof(source, pos, len)
                x = scale(T, digits, -frac, neg)
                code |= OK | EOF
                @goto done
            end
            b = peekbyte(source, pos) - UInt8('0')
            if debug
                println("float 5) $b")
            end
            b > 0x09 && break
            if overflows(IntType) && digits > overflowval(IntType)
                fastseek!(source, startpos)
                return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
            end
        end
    end
    b += UInt8('0')
    if debug
        println("float 6) $(Char(b))")
    end
    # check for exponent notation
    if b == UInt8('e') || b == UInt8('E') || b == UInt8('f') || b == UInt8('F')
        pos += 1
        incr!(source)
        # error to have a "dangling" 'e'
        if eof(source, pos, len)
            code |= INVALID | EOF
            @goto done
        end
        b = peekbyte(source, pos)
        if debug
            println("float 7) $(Char(b))")
        end
        exp = zero(IntType)
        negexp = b == UInt8('-')
        if negexp || b == UInt8('+')
            pos += 1
            incr!(source)
        end
        b = peekbyte(source, pos) - UInt8('0')
        if debug
            println("float 8) $b")
        end
        if b > 0x09
            # invalid to have a "dangling" 'e'
            code |= INVALID
            @goto done
        end
        while true
            exp = IntType(10) * exp + b
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                x = scale(T, digits, ifelse(negexp, -exp, exp) - frac, neg)
                code |= OK | EOF
                @goto done
            end
            b = peekbyte(source, pos) - UInt8('0')
            if debug
                println("float 9) $b")
            end
            if b > 0x09
                x = scale(T, digits, ifelse(negexp, -exp, exp) - frac, neg)
                code |= OK
                @goto done
            end
            if overflows(IntType) && exp > overflowval(IntType)
                fastseek!(source, startpos)
                return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
            end
        end
    else
        x = scale(T, digits, -frac, neg)
        code |= OK
    end

@label done
    return x, code, pos
end

using Base.GMP, Base.GMP.MPZ

const ONES = [BigInt(1)]
const NUMS = [BigInt()]
const QUOS = [BigInt()]
const REMS = [BigInt()]
const SCLS = [BigInt()]

const BIG_E = UInt8('E')
const LITTLE_E = UInt8('e')

const bipows5 = [big(5)^x for x = 0:325]

function roundQuotient(num, den)
    @inbounds quo, rem = MPZ.tdiv_qr!(QUOS[Threads.threadid()], REMS[Threads.threadid()], num, den)
    q = Int64(quo)
    cmpflg = cmp(MPZ.mul_2exp!(rem, 1), den)
    return ((q & 1) == 0 ? 1 == cmpflg : -1 < cmpflg) ? q + 1 : q
end

maxsig(::Type{Float16}) = 2048
maxsig(::Type{Float32}) = 16777216
maxsig(::Type{Float64}) = 9007199254740992

ceillog5(::Type{Float16}) = 5
ceillog5(::Type{Float32}) = 11
ceillog5(::Type{Float64}) = 23

const F16_SHORT_POWERS = [exp10(Float16(x)) for x = 0:2ceillog5(Float16)-1]
const F32_SHORT_POWERS = [exp10(Float32(x)) for x = 0:2ceillog5(Float32)-1]
const F64_SHORT_POWERS = [exp10(Float64(x)) for x = 0:2ceillog5(Float64)-1]

pow10(::Type{Float16}, e) = (@inbounds v = F16_SHORT_POWERS[e+1]; return v)
pow10(::Type{Float32}, e) = (@inbounds v = F32_SHORT_POWERS[e+1]; return v)
pow10(::Type{Float64}, e) = (@inbounds v = F64_SHORT_POWERS[e+1]; return v)

significantbits(::Type{Float16}) = 11
significantbits(::Type{Float32}) = 24
significantbits(::Type{Float64}) = 53

bitlength(this) = GMP.MPZ.sizeinbase(this, 2)
bits(::Type{T}) where {T <: Union{Float16, Float32, Float64}} = 8sizeof(T)

BigInt!(y::BigInt, x::BigInt) = x
BigInt!(y::BigInt, x::Union{Clong,Int32}) = MPZ.set_si!(y, x)
# copied from base/gmp.jl:285
function BigInt!(y::BigInt, x::Integer)
    x == 0 && return y
    nd = ndigits(x, base=2)
    z = GMP.MPZ.realloc2!(y, nd)
    s = sign(x)
    s == -1 && (x = -x)
    x = unsigned(x)
    size = 0
    limbnbits = sizeof(GMP.Limb) << 3
    while nd > 0
        size += 1
        unsafe_store!(z.d, x % GMP.Limb, size)
        x >>>= limbnbits
        nd -= limbnbits
    end
    z.size = s*size
    z
end

function scale(::Type{T}, v, exp) where {T <: Union{Float16, Float32, Float64}}
    ms = maxsig(T)
    cl = ceillog5(T)
    if v < ms
        # fastest path
        if 0 <= exp < cl
            return T(v) * pow10(T, exp)
        elseif -cl < exp < 0
            return T(v) / pow10(T, -exp)
        end
    end
    v == 0 && return zero(T)
    @inbounds mant = BigInt!(NUMS[Threads.threadid()], v)
    if 0 <= exp < 327
        num = MPZ.mul!(mant, bipows5[exp+1])
        bex = bitlength(num) - significantbits(T)
        bex <= 0 && return ldexp(T(num), exp)
        @inbounds one = MPZ.mul_2exp!(MPZ.set_si!(ONES[Threads.threadid()], 1), bex)
        quo = roundQuotient(num, one)
        return ldexp(T(quo), bex + exp)
    elseif -327 < exp < 0
        maxpow = length(bipows5) - 1
        @inbounds scl = SCLS[Threads.threadid()]
        if -exp <= maxpow
            MPZ.set!(scl, bipows5[-exp+1])
        else
            # this branch isn't tested
            MPZ.set!(scl, bipows5[maxpow+1])
            MPZ.mul!(scl, bipows5[-exp-maxpow+1])
        end
        bex = bitlength(mant) - bitlength(scl) - significantbits(T)
        bex = min(bex, 0)
        num = MPZ.mul_2exp!(mant, -bex)
        quo = roundQuotient(num, scl)
        if (bits(T) - leading_zeros(quo) > significantbits(T)) || exp == -324
            bex += 1
            quo = roundQuotient(num, MPZ.mul_2exp!(scl, 1))
        end
        if exp <= -324
            return T(ldexp(BigFloat(quo), bex + exp))
        else
            return ldexp(T(quo), bex + exp)
        end
    else
        return exp > 0 ? T(Inf) : T(0)
    end
end

@inline function scale(::Type{T}, lmant, exp, neg) where {T <: Union{Float16, Float32, Float64}}
    result = scale(T, lmant, exp)
    return ifelse(neg, -result, result)
end


text:
module DDL
import SearchLight: query
import SearchLight.Migrations: create_table, column, columns, primary_key, add_index, drop_table, add_indices
import TimeZones
using SearchLight, TimeZones
export  create_versions_trigger, create_validity_intervals_constraints, up, down

createTsrTriggerFun = """
  CREATE OR REPLACE FUNCTION f_bitempranges ()
  RETURNS trigger AS
  \$\$
       DECLARE

       BEGIN
            RAISE NOTICE 'NEW: %', NEW;
            NEW.tsrworld := tstzrange(NEW.tsworld_validfrom,NEW.tsworld_invalidfrom,'[)');
            NEW.tsrdb := tstzrange(NEW.tsdb_validfrom,NEW.tsdb_invalidfrom,'[)');
            RETURN NEW;
       END;
  \$\$ LANGUAGE 'plpgsql';
  """
  createBitempTrigger = """
  CREATE TRIGGER versions_trig
  BEFORE INSERT OR UPDATE ON validityIntervals
  FOR EACH ROW EXECUTE PROCEDURE f_bitempranges();
  """

"""
create_versions_trigger()
  propagate skalars ref_validfrom and ref_invalidfrom to range ref_valid
  skalars are needed enable foreign key constraints, which in POSTGRES cannot be declared
  onto bounds of intervals directly.
"""
function create_versions_trigger()
  q = """
  CREATE OR REPLACE FUNCTION f_versionrange ()
  RETURNS trigger AS
  \$\$
       DECLARE

       BEGIN
            RAISE NOTICE 'NEW: %', NEW;
            NEW.ref_valid := int8range(NEW.ref_validfrom,NEW.ref_invalidfrom,'[)');
            RETURN NEW;
       END;
  \$\$ LANGUAGE 'plpgsql';
  """
  SearchLight.query(q)
end

"""
create_validity_intervals_constraints()
  Install the GIST constraint that prevents creation of overlapping intervals, i.e.
  that at most one version is valid for a given 2 dimensional point in transaction and reference time
"""
function create_validity_intervals_constraints()
  q="""
  ALTER TABLE validityIntervals
  ADD CONSTRAINT bitemp EXCLUDE USING GIST (ref_version WITH =, is_committed WITH =, tsrworld WITH &&, tsrdb WITH &&)
  """
  SearchLight.query(q)
end

"""
up()
callback from SearchLight.Migrations.up
creating the Database Schema for bitemporal Transactions
including test dummy tables
"""
function up()
  createGistExtension = "CREATE EXTENSION IF NOT EXISTS btree_gist;"

  create_table(:histories) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:dummy, :integer)
    ]
  end

  create_table(:versions) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:ref_history, :bigint,"REFERENCES histories(id) ON DELETE CASCADE")
    ]
  end

  create_table(:workflows) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:ref_history, :bigint,"REFERENCES histories(id) ON DELETE CASCADE")
      column(:ref_version, :bigint,"REFERENCES versions(id) ON DELETE CASCADE")
      column(:tsw_validfrom, :timestamptz,"(3)")
      column(:tsdb_validfrom, :timestamptz,"(3)")
      column(:is_committed, :bigint)
    ]
  end

  create_table(:testtstzranges) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:i, :tstzrange)
    ]
  end

  create_table(:validityIntervals) do
      [
        column(:id,:bigserial,"PRIMARY KEY")
        column(:ref_history, :bigint,"REFERENCES histories(id) ON DELETE CASCADE")
        column(:ref_version, :bigint,"REFERENCES versions(id) ON DELETE CASCADE")
        column(:tsworld_validfrom, :timestamptz,"(3)")
        column(:tsworld_invalidfrom, :timestamptz,"(3)")
        column(:tsdb_validfrom, :timestamptz,"(3)")
        column(:tsdb_invalidfrom, :timestamptz,"(3)")
        column(:tsrworld, :tstzrange)
        column(:tsrdb, :tstzrange)
        column(:is_committed, :bigint)
      ]
  end

  create_table(:testdummyComponents) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:ref_history, :bigint,"REFERENCES histories(id) ON DELETE CASCADE")
      column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
    ]
    end

    create_table(:testdummyComponentRevisions) do
      [
        column(:id,:bigserial,"PRIMARY KEY")
        column(:ref_component, :bigint, "REFERENCES testdummyComponents(id)")
        column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
        column(:ref_invalidfrom, :bigint, "DEFAULT 9223372036854775807 REFERENCES versions(id) ON DELETE SET DEFAULT")
        column(:ref_valid, :int8range)
        column(:description, :string)
      ]
    end

  create_table(:testdummySubComponents) do
    [
      column(:id,:bigserial,"PRIMARY KEY")
      column(:ref_history, :bigint,"REFERENCES testdummyComponents(id) ON DELETE CASCADE")
      column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
      column(:ref_super, :bigint,"REFERENCES histories(id) ON DELETE CASCADE")
    ]
    end

    create_table(:testdummySubComponentRevisions) do
      [
        column(:id,:bigserial,"PRIMARY KEY")
        column(:ref_component, :bigint, "REFERENCES testdummySubComponents(id)")
        column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
        column(:ref_invalidfrom, :bigint, "DEFAULT 9223372036854775807 REFERENCES versions(id) ON DELETE SET DEFAULT")
        column(:ref_valid, :int8range)
        column(:description, :string)
      ]
    end

    createTestdummyComponentRevisionsTrigger = """
    CREATE TRIGGER tr_versions_trig
    BEFORE INSERT OR UPDATE ON testdummyComponentRevisions
    FOR EACH ROW EXECUTE PROCEDURE f_versionrange();
    """

    createTestdummyComponentRevisionsConstraints = """
    ALTER TABLE testdummycomponentrevisions
    ADD CONSTRAINT testdummycomponentrevisionsversionrange EXCLUDE USING GIST (ref_component WITH =, ref_valid WITH &&)
    """

    createTestdummySubComponentRevisionsTrigger = """
    CREATE TRIGGER tr_versions_trig
    BEFORE INSERT OR UPDATE ON testdummySubComponentRevisions
    FOR EACH ROW EXECUTE PROCEDURE f_versionrange();
    """

    createTestdummySubComponentRevisionsConstraints = """
    ALTER TABLE testdummySubComponentrevisions
    ADD CONSTRAINT testdummysubcomponentrevisionsversionrange EXCLUDE USING GIST (ref_component WITH =, ref_valid WITH &&)
    """

SearchLight.query(createGistExtension)
create_validity_intervals_constraints()
SearchLight.query(createTsrTriggerFun)
SearchLight.query(createBitempTrigger)
create_versions_trigger()
SearchLight.query(createTestdummyComponentRevisionsTrigger)
SearchLight.query(createTestdummyComponentRevisionsConstraints)
SearchLight.query(createTestdummySubComponentRevisionsTrigger)
SearchLight.query(createTestdummySubComponentRevisionsConstraints)
MaxDate =  ZonedDateTime(DateTime(2038, 1, 19,14,7), tz"UTC")
MaxDateSQL = SQLInput(MaxDate)
InfinityKey = 9223372036854775807 ::Integer

SearchLight.query("""
INSERT INTO histories VALUES($InfinityKey,0)
"""
)
SearchLight.query("""
INSERT INTO versions VALUES($InfinityKey)
"""
)

end

"""
down()
callback from SearchLight.Migrations.down
tearing down the Database Schema for bitemproal Transactions
including test dummy tables
"""
function down()
  drop_table(:testdummySubComponentRevisions)
  drop_table(:testdummySubComponents)
  drop_table(:testdummyComponentRevisions)
  drop_table(:testdummyComponents)
  drop_table(:validityIntervals)
  drop_table(:workflows)
  drop_table(:versions)
  drop_table(:histories)
end
end #module


text:
module Unroll

copy_and_substitute_tree(e, varname, newtext) = e

copy_and_substitute_tree(e::Symbol, varname, newtext) =
    e == varname? newtext : e

function copy_and_substitute_tree(e::Expr, varname, newtext)
    e2 = Expr(e.head)
    for subexp in e.args
        push!(e2.args, copy_and_substitute_tree(subexp, varname, newtext))
    end
    if e.head == :if
        newe = e2
        try
            u = eval(current_module(), e2.args[1])
            if u == true
                newe = e2.args[2]
            elseif u == false
                if length(e2.args) == 3
                    newe = e2.args[3]
                else
                    newe = :nothing
                end
            end
        end
        e2 = newe
    end
    e2
end

macro unroll(expr)
    if expr.head != :for || length(expr.args) != 2 ||
        expr.args[1].head != :(=) ||
        typeof(expr.args[1].args[1]) != Symbol ||
        expr.args[2].head != :block
        error("Expression following unroll macro must be a for-loop as described in the documentation")
    end
    varname = expr.args[1].args[1]
    ret = Expr(:block)
    for k in eval(current_module(), expr.args[1].args[2])
        e2 = copy_and_substitute_tree(expr.args[2], varname, k)
        push!(ret.args, e2)
    end
    esc(ret)
end

macro tuplegen(expr)
    if expr.head != :comprehension || length(expr.args) != 2 ||
        expr.args[2].head != :(=) ||
        typeof(expr.args[2].args[1]) != Symbol
        error("Expression following tuplegen macro must be a comprehension as described in the documentation")
    end
    varname = expr.args[2].args[1]
    ret = Expr(:tuple)
    for k in eval(current_module(), expr.args[2].args[2])
        e2 = copy_and_substitute_tree(expr.args[1], varname, k)
        push!(ret.args,e2)
    end
    esc(ret)
end

export @unroll
export @tuplegen

end


text:
using Test

import MutableArithmetics
const MA = MutableArithmetics

using MathOptInterface
const MOI = MathOptInterface

_zero(::Type{MOI.SingleVariable}, T::Type) = zero(MOI.ScalarAffineFunction{T})
_zero(x, ::Type) = zero(x)
function promote_operation_test(op::Function, T, x::Type, y::Type)
    f() = MA.promote_operation(op, x, y)
    @test typeof(op(_zero(x, T), _zero(y, T))) == f()
    @test 0 == @allocated f()
end

@testset "promote_operation allocation with $T" for T in [Int, Float64, Float32]
    AffType = MOI.ScalarAffineFunction{T}
    QuadType = MOI.ScalarQuadraticFunction{T}
    for op in [+, -, *]
        promote_operation_test(op, T, T, MOI.SingleVariable)
        promote_operation_test(op, T, MOI.SingleVariable, T)
        promote_operation_test(op, T, T, AffType)
        promote_operation_test(op, T, AffType, T)
        promote_operation_test(op, T, T, QuadType)
        promote_operation_test(op, T, QuadType, T)
        promote_operation_test(op, T, MOI.SingleVariable, AffType)
        promote_operation_test(op, T, AffType, MOI.SingleVariable)
        if op != *
            promote_operation_test(op, T, MOI.SingleVariable, QuadType)
            promote_operation_test(op, T, QuadType, MOI.SingleVariable)
            promote_operation_test(op, T, AffType, QuadType)
            promote_operation_test(op, T, QuadType, AffType)
        end
    end
end

@testset "promote_operation with $T" for T in [Int, Float64, Float32]
    @test MA.promote_operation(*, MOI.SingleVariable, T) == MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, T, MOI.SingleVariable) == MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, MOI.ScalarAffineFunction{T}, T) == MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, T, MOI.ScalarAffineFunction{T}) == MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, MOI.ScalarQuadraticFunction{T}, T) == MOI.ScalarQuadraticFunction{T}
    @test MA.promote_operation(*, T, MOI.ScalarQuadraticFunction{T}) == MOI.ScalarQuadraticFunction{T}
end

@testset "scaling with $T" for T in [Float64, Float32]
    x = MOI.VariableIndex(1)
    fx = MOI.SingleVariable(x)
    @test T(3) == MA.scaling(T(0)fx + T(3))
    f = T(2)fx + T(3)
    err = InexactError(:convert, T, f)
    @test_throws err MA.scaling(f)
end

@testset "Unary `-` with $T" for T in [Float64, Float32]
    x = MOI.VariableIndex(1)
    fx = MOI.SingleVariable(x)
    for f in [T(2)fx + T(3), T(4)*fx*fx + T(2)fx + T(3)]
        g = -f
        @test g ≈ MA.operate!(-, f)
        @test g ≈ f
        @test -g ≈ MA.operate!(-, f)
        @test -g ≈ f
    end
end

function all_tests(T::Type, a, b, c, d, e, f, g)
    exclude = ["scalar_uniform_scaling", "cube"]
    @testset "Scalar" begin
        MA.Test.scalar_test(a, exclude = exclude)
        MA.Test.scalar_test(b, exclude = exclude)
        MA.Test.scalar_test(c, exclude = exclude)
        MA.Test.scalar_test(d, exclude = exclude)
        MA.Test.scalar_test(e, exclude = exclude)
        MA.Test.scalar_test(f, exclude = exclude)
        MA.Test.scalar_test(g, exclude = exclude)
    end
    if !(a isa MOI.ScalarQuadraticFunction)
        exclude = String[]
        if T <: Integer
            push!(exclude, "quadratic_division")
        end
        @testset "Quadratic" begin
            MA.Test.quadratic_test(a, b, c, d, exclude = exclude)
            MA.Test.quadratic_test(b, c, d, e, exclude = exclude)
            MA.Test.quadratic_test(c, d, e, f, exclude = exclude)
            MA.Test.quadratic_test(d, e, f, g, exclude = exclude)
            MA.Test.quadratic_test(e, f, g, a, exclude = exclude)
            MA.Test.quadratic_test(f, g, a, b, exclude = exclude)
            MA.Test.quadratic_test(g, a, b, c, exclude = exclude)
        end
    end
    exclude = String[]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "sparse_quadratic")
    end
    @testset "Sparse" begin
        MA.Test.sparse_test(a, b, [a b c; b c a; a b a], exclude = exclude)
    end
    exclude = String[]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "matrix_vector")
        push!(exclude, "dot")
    end
    @testset "Vector" begin
        MA.Test.array_test([a, b, c], exclude = exclude)
        MA.Test.array_test([b, c, d], exclude = exclude)
        MA.Test.array_test([c, d, e], exclude = exclude)
        MA.Test.array_test([d, e, a], exclude = exclude)
        MA.Test.array_test([e, a, b], exclude = exclude)
    end
    exclude = ["broadcast_division", "matrix_uniform_scaling", "symmetric_matrix_uniform_scaling", "matrix_vector_division"]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "dot")
        push!(exclude, "broadcast_multiplication")
        push!(exclude, "sum_multiplication")
    end
    @testset "Matrix" begin
        MA.Test.array_test([a b; c d], exclude = exclude)
        MA.Test.array_test([c e; e d], exclude = exclude)
        MA.Test.array_test([a b c; b c a; a b a], exclude = exclude)
        MA.Test.array_test([d b c; d c e; e b a], exclude = exclude)
    end
end

x = MOI.VariableIndex(1)
y = MOI.VariableIndex(2)
fx = MOI.SingleVariable(x)
fy = MOI.SingleVariable(y)

@testset "Affine in $T" for T in [Int]
    @testset "Int" begin
        MA.Test.int_test(MOI.ScalarAffineFunction{T}, exclude = ["int_mul", "int_add_mul"])
    end
    a = T(2) * fx + T(1)
    b = T(4) * fy + T(2)
    c = T(3) * fx - T(2) * fy - T(3)
    d = T(3) * fx - T(2) * fy + T(4) * fx
    e = T(5) * fx - T(5)
    f = T(1) * fy - T(2) * fx + T(2)
    g = T(2) * fx + T(3) * fx
    all_tests(T, a, b, c, d, e, f, g)
end

@testset "MutableArithmetics with polynomials in $T" for T in [Int]
    @testset "Int" begin
        MA.Test.int_test(MOI.ScalarQuadraticFunction{T}, exclude = ["int_mul", "int_add_mul"])
    end
    a = T(2) * fx + T(1) + T(4) * fx * fy
    b = T(4) * fy + T(3) * fy * fy - T(3) * fy * fx + T(2)
    c = T(2) * fx * fx + T(3) * fx - T(2) * fy - T(3)
    d = T(2) * fx * fx + T(3) * fx - T(2) * fy + T(4) * fx - T(1) * fy * fy
    e = T(5) * fx - T(5) - T(4) * fx * fy
    f = T(1) * fy + T(2) * fy * fy - T(2) * fx + T(2)
    g = T(2) * fx + T(3) * fx + T(3) * fx * fy
    all_tests(T, a, b, c, d, e, f, g)
end


text:
using SummarizedExperiments
using Documenter
import DataStructures
import DataFrames

makedocs(
    sitename="SummarizedExperiments.jl",
    modules = [SummarizedExperiments],
    pages=[
        "Home" => "index.md"
    ]
)

deploydocs(;
    repo="github.com/LTLA/SummarizedExperiments.jl",
)


text:
###############################################################################
#
# Calculate soil water content from pressure
#
###############################################################################
"""
    soil_erwc(sh::AbstractSoilVC{FT}, p_25::FT) where {FT<:AbstractFloat}

Returns the Effective relative water content of the soil
    ``\\frac{Θ - Θr}{Θs - Θr}``, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `p_25` Matrix water potential equivalent to 25 degree C, with surface tension
correction
"""
function soil_erwc(
            sh::BrooksCorey{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    if p_25 < 0
        @unpack b, ϕs = sh;

        return (-ϕs/p_25) ^ (1/b)
    else
        return FT(1)
    end
end

function soil_erwc(
            sh::VanGenuchten{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    if p_25 < 0
        @unpack m, n, α = sh;

        return ( 1 / ( 1 + (-p_25 * α) ^ n ) ) ^ m
    else
        return FT(1)
    end
end

"""
    soil_rwc(sh::AbstractSoilVC{FT}, p_25::FT) where {FT<:AbstractFloat}

Returns the relative soil water content, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `p_25` Matrix water potential equivalent to 25 degree C, with surface tension
correction
"""
function soil_rwc(
            sh::BrooksCorey{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    @unpack Θs = sh;

    return soil_swc(sh, p_25) / Θs
end

function soil_rwc(
            sh::VanGenuchten{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    @unpack Θs = sh;

    return soil_swc(sh, p_25) / Θs
end

"""
    soil_swc(sh::AbstractSoilVC{FT}, p_25::FT) where {FT<:AbstractFloat}

Returns the soil water content, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `p_25` Matrix water potential equivalent to 25 degree C, with surface tension
correction
"""
function soil_swc(
            sh::BrooksCorey{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    return soil_erwc(sh, p_25) * (Θs - Θr) + Θr
end

function soil_swc(
            sh::VanGenuchten{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    return soil_erwc(sh, p_25) * (Θs - Θr) + Θr
end

###############################################################################
#
# Calculate soil k from water content
#
###############################################################################
"""
    soil_k_ratio_erwc(
                sh::AbstractSoilVC{FT},
                erwc::FT
    ) where {FT<:AbstractFloat}

Return the soil k relative to maximal k, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `erwc` Effective relative soil water content (``\\frac{Θs - Θr}{Θs - Θ}``)
"""
function soil_k_ratio_erwc(
            sh::BrooksCorey{FT},
            erwc::FT
) where {FT<:AbstractFloat}
    k_ratio = erwc ^ (2*sh.b+3);

    return max(k_ratio, FT(1e-20))
end

function soil_k_ratio_erwc(
            sh::VanGenuchten{FT},
            erwc::FT
) where {FT<:AbstractFloat}
    @unpack m = sh;
    k_ratio = sqrt(erwc) * (1 - (1 - erwc^(1/m)) ^ m)^2;

    return max(k_ratio, FT(1e-20))
end

"""
    soil_k_ratio_rwc(sh::AbstractSoilVC{FT}, rwc::FT) where {FT<:AbstractFloat}

Return the soil k relative to maximal k, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `rwc` Relative soil water content
"""
function soil_k_ratio_rwc(
            sh::BrooksCorey{FT},
            rwc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (rwc * Θs - Θr) / (Θs - Θr)));

    return soil_k_ratio_erwc(sh, erwc)
end

function soil_k_ratio_rwc(
            sh::VanGenuchten{FT},
            rwc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (rwc * Θs - Θr) / (Θs - Θr)));

    return soil_k_ratio_erwc(sh, erwc)
end

"""
    soil_k_ratio_swc(sh::AbstractSoilVC{FT}, swc::FT) where {FT<:AbstractFloat}

Return the soil k relative to maximal k, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `swc` Relative soil water content
"""
function soil_k_ratio_swc(
            sh::BrooksCorey{FT},
            swc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (swc - Θr) / (Θs - Θr)));

    return soil_k_ratio_erwc(sh, erwc)
end

function soil_k_ratio_swc(
            sh::VanGenuchten{FT},
            swc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (swc - Θr) / (Θs - Θr)));

    return soil_k_ratio_erwc(sh, erwc)
end

###############################################################################
#
# Calculate soil k from pressure
#
###############################################################################
"""
    soil_k_ratio_p25(
                sh::AbstractSoilVC{FT},
                p_25::FT
    ) where {FT<:AbstractFloat}

Return the soil k relative to maximal k, given
- `sh` [`AbstractSoilVC`](@ref) type soil hydraulics
- `p_25` Matrix water potential equivalent to 25 degree C, with surface tension
"""
function soil_k_ratio_p25(
            sh::AbstractSoilVC{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    erwc = soil_erwc(sh, p_25);

    return soil_k_ratio_erwc(sh, erwc)
end

###############################################################################
#
# Calculate soil pressure from water content
#
###############################################################################
"""
    soil_p_25_erwc(sh::AbstractSoilVC{FT}, erwc::FT) where {FT<:AbstractFloat}

Returns the Relative water content of the soil, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `erwc` Effectibe relative soil water content (``\\frac{Θs - Θr}{Θs - Θ}``)

Note that this function returns the matrix potential of water, not water
potential. Also, the potential is that at 25 degree C, not yet been corrected
for temperature effect on surface tension.
"""
function soil_p_25_erwc(
            sh::BrooksCorey{FT},
            erwc::FT
) where {FT<:AbstractFloat}
    if erwc < 1
        @unpack b, ϕs = sh;

        return -ϕs / (erwc ^ b)
    else
        return FT(0)
    end
end

function soil_p_25_erwc(
            sh::VanGenuchten{FT},
            erwc::FT
) where {FT<:AbstractFloat}
    if erwc < 1
        @unpack m, n, α = sh;

        return -1 * (erwc ^ (-1/m) - 1) ^ (1/n) / α
    else
        return FT(0)
    end
end

"""
    soil_p_25_rwc(sh::AbstractSoilVC{FT}, rwc::FT) where {FT<:AbstractFloat}

Returns the Relative water content of the soil, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `rwc` Relative soil water content

Note that this function returns the matrix potential of water, not water
potential. Also, the potential is that at 25 degree C, not yet been corrected
for temperature effect on surface tension.
"""
function soil_p_25_rwc(
            sh::BrooksCorey{FT},
            rwc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (rwc * Θs - Θr) / (Θs - Θr)));

    return soil_p_25_erwc(sh, erwc)
end

function soil_p_25_rwc(
            sh::VanGenuchten{FT},
            rwc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (rwc * Θs - Θr) / (Θs - Θr)));

    return soil_p_25_erwc(sh, erwc)
end

"""
    soil_p_25_swc(sh::AbstractSoilVC{FT}, rwc::FT) where {FT<:AbstractFloat}

Returns the Relative water content of the soil, given
- `sh` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type soil hydraulics
- `swc` Soil water content

Note that this function returns the matrix potential of water, not water
potential. Also, the potential is that at 25 degree C, not yet been corrected
for temperature effect on surface tension.
"""
function soil_p_25_swc(
            sh::BrooksCorey{FT},
            swc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (swc - Θr) / (Θs - Θr)));

    return soil_p_25_erwc(sh, erwc)
end

function soil_p_25_swc(
            sh::VanGenuchten{FT},
            swc::FT
) where {FT<:AbstractFloat}
    @unpack Θr, Θs = sh;

    erwc = min(1, max(0, (swc - Θr) / (Θs - Θr)));

    return soil_p_25_erwc(sh, erwc)
end

###############################################################################
#
# Fit Soil VC
#
###############################################################################
"""
    fit_soil_VC!(
                vc_vG::VanGenuchten{FT},
                vc_BC::BrooksCorey{FT}
    ) where {FT<:AbstractFloat}

Update BrooksCorey setup from known VanGenuchten parameters, given
- `vc_vG` [`VanGenuchten`](@ref) type soil VC
- `vc_BC` [`BrooksCorey`](@ref) type soil VC
"""
function fit_soil_VC!(
            vc_vG::VanGenuchten{FT},
            vc_BC::BrooksCorey{FT}
) where {FT<:AbstractFloat}
    # generate data to fit
    _Θ    = range(vc_vG.Θr+FT(1e-2); stop=vc_vG.Θs-FT(1e-2), length=30);
    _P_vG = -1 .* soil_p_25_swc.([vc_vG], _Θ);
    _P_BC = -1 .* soil_p_25_swc.([vc_BC], _Θ);

    f(x) = (
        vc_BC.b  = x[1];
        vc_BC.ϕs = x[2];
        _P_BC   .= -1 .* soil_p_25_swc.([vc_BC], _Θ);
        _diff    = sum( (log.(_P_BC) .- log.(_P_vG)) .^ 2 );
        return -_diff
    );

    st = SolutionToleranceND{FT}([1e-3, 1e-6], 30);
    ms = ReduceStepMethodND{FT}(x_mins = FT[1e-3, 1e-6],
                                x_maxs = FT[ 100, 1000],
                                x_inis = [vc_BC.b, vc_BC.ϕs],
                                Δ_inis = FT[0.1, 1e-3]);
    bc = find_peak(f, ms, st);
    vc_BC.b  = bc[1];
    vc_BC.ϕs = bc[2];

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using BenchmarkTools;

# is it an issue?

function L0_tup(arg1; kwarg1, kwarg2, kwarg3, kwarg4)
  arg1 + kwarg1 + kwarg2 + kwarg3 + kwarg4
end

function L0_arr(arg1, kwarg_arr)
  kwarg1, kwarg2, kwarg3, kwarg4 = kwarg_arr
  arg1 + kwarg1 + kwarg2 + kwarg3 + kwarg4
end

function L0_tup_wrapped(arg1, kwarg_arr)
  symbs = (:kwarg1, :kwarg2, :kwarg3, :kwarg4)
  return L0_tup(arg1; NamedTuple{symbs}(kwarg_arr)...)
end

# test array
vars = [5, 4, 3, 2]

# @btime L0_arr(5, $vars) # 1.808 ns
# @btime L0_tup_wrapped(5, $vars) # 468 ns

# the typ overhead adds ~ 400 ns overhead. Is this signifcant relative to Mittag-Leffler function?

using MittagLeffler
t = -Float64.(0.0:0.001:1.0) # 1000 time points, representative data set
# @btime mittleff.(0.75, $t); # 37.430 ms, ~100000x more time required than for creating the tuple
# @btime mittleff.(0.25, $t); # 39.961 ms, ~100000x more time required than for creating the tuple

# the time for sending array -> namedtuple is negligible compared to mittleff.

# a neatened up PoC of the original
-----------------------------------
struct dummy_type{T<:Function}
  free_params::Tuple
  fixed_params::NamedTuple
  afunc::T
end

function freeze_params(orig::dummy_type, tofix::NamedTuple)
  newfree =  Tuple([i for i in orig.free_params if !(i in keys(tofix))])
  newfixed = merge(orig.fixed_params, tofix)
  return dummy_type(newfree, newfixed, orig.afunc)
end

function get_frozen_func(orig::dummy_type)
  f = let func = orig.afunc, fixed = orig.fixed_params
    (t; kwargs...) -> func(t; kwargs..., fixed...)
  end
  return f
end

function func1(t; a, b, c)
  t + a + b + c
end

t1 = dummy_type((:a, :b, :c), NamedTuple{}(), func1)
t2 = freeze_params(t1, (a = 1,))
func2 = get_frozen_func(t2)

full_kwargs = (a = 1, b = 2, c = 3)
part_kwargs = (b = 2, c = 3)

@btime func1.(t; $full_kwargs...); # 1.085 us
@btime func2.(t; $part_kwargs...); # 1.085 us, exactly the same

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
# revisit after issue #14 is closed
# ""
# function constraint_gen_contingency_ptdf_thermal_limit_from_soft(pm::_PM.AbstractDCPModel, n::Int, i::Int, cut_map, rate, gen_set, gen_alpha, gen_bus)
# end

""
function constraint_voltage_magnitude_link(pm::_PM.AbstractDCPModel, n_1::Int, n_2::Int, i::Int)
    # do nothing voltages are assumed to be the same
end

function constraint_ohms_yt_from_goc(pm::_PM.AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr  = var(pm, n,  :p, f_idx)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    JuMP.@constraint(pm.model, p_fr == -b*(va_fr - va_to))
    # omit reactive constraint
end


text:
"""
    M_ABG(r, rs, rhos, alpha, beta, gamma)

Compute the enclosed mass profile for an αβγ model.

* `r`: radii in kpc
* `rs`: scale radius in kpc
* `rhos`: scale density in Msun / kpc3
* `alpha`: transition sharpness
* `beta`: negative outer log slope
* `gamma`: negative inner log slope
"""
function M_ABG(r, rs, rhos, alpha, beta, gamma)
    x = r / rs
    y = 4π * rhos * rs^3
    omega = 3.0 - gamma
    hyp = hyp2f1.(omega / alpha,
                 (beta - gamma) / alpha,
                 1 + omega / alpha,
                 -(x .^ alpha))
    return @. y * x ^ omega * hyp / omega
end

"""
    rho_ABG(r, rs, rhos, alpha, beta, gamma)

Compute the local volume density for an αβγ model.

* `r`: radii in kpc
* `rs`: scale radius in kpc
* `rhos`: scale density in Msun / kpc3
* `alpha`: transition sharpness
* `beta`: negative outer log slope
* `gamma`: negative inner log slope
"""
function rho_ABG(r, rs, rhos, alpha, beta, gamma)
    x = r / rs
    return @. rhos * (x ^ -gamma * (1 + x ^ alpha) ^ ((-beta - gamma) / alpha))
end

"""
    drhodr_ABG(r, rs, rhos, alpha, beta, gamma)

Compute the linear density derivative for αβγ model.

* `r`: radii in kpc
* `rs`: scale radius in kpc
* `rhos`: scale density in Msun / kpc3
* `alpha`: transition sharpness
* `beta`: negative outer log slope
* `gamma`: negative inner log slope
"""
function drhodr_ABG(r, rs, rhos, alpha, beta, gamma)
    x = r / rs
    α = alpha
    β = beta
    γ = gamma
    term1 = @. -γ * x ^ (-γ - 1) * (1 + x ^ α)^(-(β - γ) / α)
    term2 = @. -(β - γ) * x ^ (α - γ - 1) * (1 + x ^ α)^(-(β - γ) / α - 1)
    return @. rhos / rs * (term1 + term2)
end

"""
    ABGModel(rs, rhos, alpha, beta, gamma)

αβγ double-power law density model.

[Merritt et al. 2006](http://adsabs.harvard.edu/abs/2006AJ....132.2685M)

* `rs`: scale radius in kpc
* `rhos`: scale density in Msun / kpc3
* `alpha`: transition sharpness
* `beta`: negative outer log slope
* `gamma`: negative inner log slope
"""
struct ABGModel <: HaloModel
    rs::Float64
    rhos::Float64
    alpha::Float64
    beta::Float64
    gamma::Float64
end

# default to 1e12 Mvir / Msun NFW halo
ABGModel() = ABGModel(25.3, 3.7e6, 1, 3, 1)
density(halo::ABGModel, r) = rho_ABG(r, halo.rs, halo.rhos,
                                     halo.alpha, halo.beta, halo.gamma)
mass(halo::ABGModel, r) = M_ABG(r, halo.rs, halo.rhos,
                                halo.alpha, halo.beta, halo.gamma)
scale_radius(halo::ABGModel) = begin
    ((2 - halo.gamma) / (halo.beta - 2)) ^ (1 / halo.alpha) * halo.rs
end

"""
    ABG_from_virial(Mvir, cvir, alpha, beta, gamma; <keyword arguments>)

Construct a ABGModel halo from the virial mass, concentration, and shape parameters.

# Arguments
- `mdef::AbstractString = default_mdef`: halo mass definition (e.g., "200c", "vir")
- `cosmo::AbstractCosmology = default_cosmo`: cosmology under which to evaluate the overdensity
- `z::Real = 0.0`: redshift at which to evaluate the overdensity
"""
function ABG_from_virial(Mvir, cvir, alpha, beta, gamma;
                         mdef = default_mdef,
                         cosmo = default_cosmo,
                         z = 0.0)
    Rvir = Rvir_from_Mvir(Mvir; mdef = mdef, cosmo = cosmo, z = z)
    r2 = Rvir / cvir
    rs = ((2 - gamma) / (beta - 2))^(-1 / alpha) * r2
    rhos = Mvir / M_ABG(Rvir, rs, 1.0, alpha, beta, gamma)
    return ABGModel(rs, rhos, alpha, beta, gamma)
end

"""
    ABG_from_virial(Mvir, cvir, Mstar; <keyword arguments>)

Construct a ABGModel halo from the virial mass, concentration, and stellar mass, using the
scaling relations from DiCintio et al. 2014

# Arguments
- `mdef::AbstractString = default_mdef`: halo mass definition (e.g., "200c", "vir")
- `cosmo::AbstractCosmology = default_cosmo`: cosmology under which to evaluate the overdensity
- `z::Real = 0.0`: redshift at which to evaluate the overdensity
"""
function ABG_from_virial(Mvir, cvir, Mstar;
                         mdef = default_mdef,
                         cosmo = default_cosmo,
                         z = 0.0)
    Rvir = Rvir_from_Mvir(Mvir; mdef = mdef, cosmo = cosmo, z = z)
    alpha, beta, gamma = abg_from_logshm(log10(Mstar / Mvir))
    r2 = Rvir / cvir
    rs = ((2 - gamma) / (beta - 2))^(-1 / alpha) * r2
    rhos = Mvir / M_ABG(Rvir, rs, 1.0, alpha, beta, gamma)
    return ABGModel(rs, rhos, alpha, beta, gamma)
end


text:
using Propagators, Test
using Propagators: Cell, propagator, propagate!

let
  inc(x) = x+1

  x = Cell()
  y = Cell()

  propagator(inc)(x, y)

  x[] = 3
  propagate!()

  @test y[] == 4
end

include("../demos/sqrt.jl")

let
  x = Cell()
  y = Cell()
  sqrt_network(x, y)
  x[] = 2
  propagate!()
  @test y[] ≈ sqrt(2)
end

using Propagators: constraint, quadratic

let
  x = Cell()
  y = Cell()
  quadratic(x, y)
  x[] = 2
  propagate!()
  @test y[] == 4

  x = Cell()
  y = Cell()
  quadratic(x, y)
  y[] = 4
  propagate!()
  @test x[] == 2
end

let
  x = Cell()
  y = Cell()
  z = Cell()
  constraint(*)(x, y, z)
  x[] = 2
  y[] = 3
  propagate!()
  @test z[] == 6

  x = Cell()
  y = Cell()
  z = Cell()
  constraint(*)(x, y, z)
  x[] = 2
  z[] = 6
  propagate!()
  @test y[] == 3
end

let
  x = Cell()
  y = Cell()
  z = Cell()
  constraint(*)(x, y, z)
  x[] = 2..5
  y[] = 3..6
  propagate!()
  @test z[] == x[]*y[]

  z[] = interval(5, 10)
  length(Propagators.workqueue)
  propagate!()
  @test x[] == interval(2, 10/3)
  @test y[] == interval(3, 5)
end


text:
using Documenter
using GoogleDrive

makedocs(

         sitename="GoogleDrive",
         modules =[GoogleDrive],
          pages = [
        "Home" => "index.md",
        "Features" => "features.md",

    ],
    )

deploydocs(repo = "github.com/tejasvaidhyadev/GoogleDrive.jl.git",
           target = "build",
           branch = "gh-pages",
           latest = "mater",
           osname = "linux"
    push_preview = true,
)


text:
using VisualStringDistances
using Test
using VisualStringDistances: glyph!, Glyph, GlyphCoordinates, ConstantVector
using UnbalancedOptimalTransport: fdot, KL, sinkhorn_divergence!
using LinearAlgebra: dot

printglyph_dashes = (io, g) -> printglyph(io, g; symbols=("#", "-"))
@testset "VisualStringDistances.jl" begin

    @testset "Glyphs" begin
        g = glyph!(hex2bytes("0000000018242442427E424242420000"))
        @test sprint(printglyph_dashes, g) == """
             --------
             --------
             --------
             --------
             ---##---
             --#--#--
             --#--#--
             -#----#-
             -#----#-
             -######-
             -#----#-
             -#----#-
             -#----#-
             -#----#-
             --------
             --------
             """

        @test_throws ArgumentError glyph!(hex2bytes("0000000018242442427E424242420"))

        @test_throws ErrorException Glyph(Char(0x12480))

        g = glyph!(hex2bytes("00000000000003C0042004200840095008E01040100010002000200000000000"))
        @test sprint(printglyph_dashes, g) == """
            ----------------
            ----------------
            ----------------
            ------####------
            -----#----#-----
            -----#----#-----
            ----#----#------
            ----#--#-#-#----
            ----#---###-----
            ---#-----#------
            ---#------------
            ---#------------
            --#-------------
            --#-------------
            ----------------
            ----------------
            """
    end

    @testset "Printing abc many ways" begin

        abc_printed_rep = """
                          ------------------------
                          ------------------------
                          ------------------------
                          ---------#--------------
                          ---------#--------------
                          ---------#--------------
                          --####---#-###----####--
                          -#----#--##---#--#----#-
                          ------#--#----#--#------
                          --#####--#----#--#------
                          -#----#--#----#--#------
                          -#----#--#----#--#------
                          -#---##--##---#--#----#-
                          --###-#--#-###----####--
                          ------------------------
                          ------------------------
                          """

        @test sprint(printglyph_dashes, Glyph("abc")) == abc_printed_rep
        @test sprint(printglyph_dashes, "abc") == abc_printed_rep
        @test sprint(printglyph_dashes, hcat(Glyph("a"), Glyph("bc"))) == abc_printed_rep
        @test sprint(printglyph_dashes, GlyphCoordinates("abc")) == abc_printed_rep

        abc_substring = Glyph(SubString("abcd", 1:3))
        @test sprint(printglyph_dashes, abc_substring) == abc_printed_rep
    end

    @testset "More GlyphCoordinates" begin

        @test GlyphCoordinates('a') ==
              GlyphCoordinates("a") ==
              GlyphCoordinates{Float64}("a")
        @test length(GlyphCoordinates('a')) ≈ sum(!iszero, Glyph("a"))
        # test indexing
        @test GlyphCoordinates('a')[1] == collect(Tuple(findfirst(!iszero, Glyph("a"))))

        gc =  GlyphCoordinates('a')
        @test gc[1:length(gc)] == gc[:] == gc.v

        gcF32 =  GlyphCoordinates{Float32}('a')
        @test gcF32[:] ≈ gc[:]
    end

    @testset "ConstantVector" begin
        for constant in (5.0, 2.2 + im * 3.2, 1f0)
            T = typeof(constant)
            c = ConstantVector{constant,T}(10)
            @test collect(c) isa Vector{T}
            @test length(c) == 10
            @test c == fill(constant, 10) == collect(c)
            @test sum(c) ≈ sum(collect(c))
            f(x) = sin(x) + 10.0
            d = randn(10)
            @test fdot(f, c, d) ≈ fdot(f, collect(c), d)
            @test fdot(f, d, c) ≈ fdot(f, d, collect(c))

            @test dot(c, d) ≈ dot(collect(c), d)
            @test dot(d, c) ≈ dot(d, collect(c))
        end
    end

    @testset "`visual_distance`" begin
        for T in (Float32, Float64), ϵ in (T(0.1), T(0.2)), ρ in (T(1.0), T(5.0))
            v1 = visual_distance(T, "abc", "def"; D=KL(ρ), ϵ=ϵ)
            v2 = visual_distance(T, "def", "ghi"; D=KL(ρ), ϵ=ϵ)
            v3 = visual_distance(T, "abc", "ghi"; D=KL(ρ), ϵ=ϵ)

            @test v1 >= 0
            @test v1 ≈ visual_distance(T, "def", "abc"; D=KL(ρ), ϵ=ϵ) rtol = 1e-3

            # Note: triangle inequality doesn't necessary hold in general
            # (it's not proven, as far as I know)
            # However, it does in this case!
            @test v3 <= v1 + v2

            v4 = visual_distance(T, "abc", "abd"; D=KL(ρ), ϵ=ϵ)
            @test v4 <= v1
        end

        # Defaults
        v1 = visual_distance("abc", "def")
        v2 = visual_distance("def", "ghi")
        v3 = visual_distance("abc", "ghi")

        @test v1 >= 0
        @test v1 ≈ visual_distance("def", "abc") rtol = 1e-3

        @test v3 <= v1 + v2

        v4 = visual_distance("abc", "abd")
        @test v4 <= v1

        abc_measure = word_measure("abc")
        def_measure = word_measure("def")
        @test v1 ≈ sinkhorn_divergence!(KL(1.0), abc_measure, def_measure, 0.1)

        # Make sure we can use non-String types
        abc_substring = SubString("abcd", 1:3)
        def_substring = SubString("defh", 1:3)
        @test v1 ≈ visual_distance(abc_substring, def_substring)

        # Test normalization
        @test visual_distance("abc", "def", normalize=sqrt) ≈ v1 / sqrt(3)
        @test visual_distance("abc", "def", normalize=identity) ≈ v1 / 3
    end
end


text:
import LightGraphs: LevelSynchronousBFS, bfs_tree
import Base.Threads: Atomic, @threads

@testset "Parallel_BFS" begin

    @testset "Thread Queue " begin
        next = @inferred(LightGraphs.ThreadQueue(Int, 5)) # Initialize threadqueue
        @test isempty(next) == true
        push!(next, 1)
        @test next[1][] == 1
        @threads for i=2:5
            push!(next, i)
        end
        @test Set([i[] for i in next[1:5]]) == Set([1,2,3,4,5])
        first = shift!(next)
        @test first == 1
    end

    g5 = DiGraph(4)
    add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
    g6 = smallgraph(:house)

    for g in testdigraphs(g5)
      T = eltype(g)
      z = @inferred(bfs_tree(LevelSynchronousBFS(), g, T(1)))
      next = LightGraphs.ThreadQueue(T, nv(g)) # Initialize threadqueue
      parents = [Atomic{T}(0) for i=1:nv(g)] # Create parents array
      LightGraphs.bfs_tree!(LevelSynchronousBFS(), next, g, T(1), parents)
      t = [i[] for i in parents]
      @test t == [T(1),T(1),T(1),T(3)]
      @test nv(z) == T(4) && ne(z) == T(3) && !has_edge(z, 2, 3)
    end

    function istree{T<:Integer}(parents::Vector{Atomic{T}}, maxdepth, n::T)
        flag = true
        for i in one(T):n
            s = i
            depth = 0
            while parents[s][] > 0 && parents[s][] != s
                s = parents[s][]
                depth += 1
                if depth > maxdepth
                    return false
                end
            end
        end
        return flag
    end

    for g in testgraphs(g6)
        n = nv(g)
        T = eltype(g)
        next = LightGraphs.ThreadQueue(eltype(g), nv(g)) # Initialize threadqueue
        parents = [Atomic{T}(0) for i=1:nv(g)] # Create parents array
        @test length(next.data) == n
        @inferred(LightGraphs.bfs_tree!(LevelSynchronousBFS(), next, g, T(1), parents))
        @test istree(parents, n, n)
        p = [i[] for i in parents]
        t = tree(p)
        @test is_directed(t)
        @test typeof(t) <: AbstractGraph
        @test ne(t) < nv(t)
    end

end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# See https://github.com/JuliaLang/julia/pull/25102
if VERSION < v"0.7.0-DEV.3073"
  bindir = JULIA_HOME
else
  bindir = Sys.BINDIR
end

if VERSION >= v"0.7.0-DEV.3382"
  using Libdl
end

println("VERSION: $(VERSION)")

jl_share = abspath(joinpath(bindir, Base.DATAROOTDIR, "julia"))
println("JL_HOME: $(jl_share)")

libdir = abspath(dirname(Libdl.dlpath("libjulia")))
println("libdir: $(libdir)")


text:
using GViz
using Test

@testset "GViz.jl" begin
    # Write your tests here.
end


text:
# Tests for validated_integ

using TaylorModels
using LinearAlgebra: norm
using Test

setformat(:full)

# NOTE: IntervalArithmetic v0.16.0 includes this function; but
# IntervalRootFinding is bounded to use v0.15.x
interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)
interval_rand(X::IntervalBox) = interval_rand.(X)

@testset "Tests for `validated_integ`" begin
    @testset "Forward integration" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) =
            (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 0.0, 10.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, Val(2))

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ(falling_ball!, q0, δq0,
            tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end
    end

    @testset "Backward integration" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) = (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 10.0, 0.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, Val(2))

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ(falling_ball!, q0, δq0,
            tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end
    end
end


text:
using Checkpoints
using Test

setdepot(joinpath(@__DIR__, "checkpoints"))
clear()
f(x) = (sleep(x); return x^2)

@testset "Testing Checkpoint" begin

    sleeptime = 3
    file = "test1.checkpoint"
    path = joinpath(@__DIR__, "checkpoints", file)
    @test !ispath(path)
    runtime = @elapsed( result = checkpoint(() -> f(sleeptime), file) )

    @test runtime > sleeptime
    @test result == sleeptime^2
    @test ispath(path)

    # Run the function again, verify that it takes very little time and returns the same
    # result.
    runtime = @elapsed( result = checkpoint(() -> f(2 * sleeptime), file) )
    @test runtime < 1

    # Note - the results should be the same as the last invocation because we cached the
    # results. To recompute the function, we need to clean the file
    @test result == sleeptime^2

    # Test cleaning
    Checkpoints.clear(file)
    @test !ispath(path)

    runtime = @elapsed( result = checkpoint(() ->f(2 * sleeptime), file) )
    @test runtime > 2 * sleeptime
    @test result == f( 2 * sleeptime )
end

@testset "Testing Macro" begin
    sleeptime = 3
    file = "test2.checkpoint"
    path = joinpath(@__DIR__, "checkpoints", file)
    @test !ispath(path)

    runtime = @elapsed begin
        result = @checkpoint f(sleeptime) file
    end

    @test runtime > sleeptime
    @test result == f(sleeptime)
    @test ispath(path)

    runtime = @elapsed begin
        result = @checkpoint f(2 * sleeptime) file
    end

    @test runtime < 1
    @test result == f(sleeptime)

    # Test the "begin -> end" syntax
    runtime = @elapsed begin
        y = @checkpoint begin
            sleep(10)
            return 10
        end "test3.checkpoint"
    end

    @test runtime > 10
    @test y == 10

    # Use the same checkpoint directory.
    runtime = @elapsed begin
        x = @checkpoint begin
            sleep(10)
            return 20
       end "test3.checkpoint"
    end

    @test runtime < 1
    @test x == 10
end

@testset "Testing Keyword" begin
    # Test function with just keyword arguments
    f(;a = 1, b = 2, c = 3) = (a, b, c)

    x = @checkpoint f(;a = 2) "test4.checkpoint"
    @test x == (2, 2, 3)

    x = @checkpoint f(;a = 2) "test4.checkpoint"
    @test x == (2, 2, 3)

    y = @checkpoint f(;b = 4) "test5.checkpoint"
    @test y == (1,4,3)

    y = @checkpoint f(;b = 4) "test5.checkpoint"
    @test y == (1,4,3)

    z = @checkpoint f(;a = 2, c = 1) "test6.checkpoint"
    @test z == (2,2,1)

    z = @checkpoint f(;a = 2, c = 1) "test6.checkpoint"
    @test z == (2,2,1)
    clear()
end