text:
#  _____ _____ ____ _____ ____        ____        ____  ___   ___  ____
# |_   _| ____/ ___|_   _/ ___|      |  _ \      / ___|/ _ \ / _ \|  _ \
#   | | |  _| \___ \ | | \___ \ _____| |_) |____| |  _| | | | | | | | | |
#   | | | |___ ___) || |  ___) |_____|  _ <_____| |_| | |_| | |_| | |_| |
#   |_| |_____|____/ |_| |____/      |_| \_\     \____|\___/ \___/|____/

# base attribute reference types
@test attr("foo").name == "foo"
@test attr(:foo).name == "foo"

@test [e.name for e=attr("foo", "bar", "baz").attrs] == ["foo", "bar", "baz"]
@test [e.name for e=attr(:foo, :bar, :baz).attrs] == ["foo", "bar", "baz"]

@test attr(:foo)[1].attr.name == "foo"
@test attr(:foo)[1].idx == 1

#   ____                _ _ _   _                   _
#  / ___|___  _ __   __| (_) |_(_) ___  _ __   __ _| |
# | |   / _ \| '_ \ / _` | | __| |/ _ \| '_ \ / _` | |
# | |__| (_) | | | | (_| | | |_| | (_) | | | | (_| | |
#  \____\___/|_| |_|\__,_|_|\__|_|\___/|_| |_|\__,_|_|
#  _____                              _
# | ____|_  ___ __  _ __ ___  ___ ___(_) ___  _ __  ___
# |  _| \ \/ / '_ \| '__/ _ \/ __/ __| |/ _ \| '_ \/ __|
# | |___ >  <| |_) | | |  __/\__ \__ \ | (_) | | | \__ \
# |_____/_/\_\ .__/|_|  \___||___/___/_|\___/|_| |_|___/
#            |_|

function check_expression(expr, expected; attrs = Dict(), vals = Dict())
    refs = DynamoDB.refs_tracker()
    @test expected == DynamoDB.write_expression(refs, expr)

    @test refs.attrs == attrs
    @test refs.vals == vals
end

# simple expressions
check_expression(attr("foo"), "#1";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(attr("foo")[1], "#1[1]";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(attr("foo", "bar")[1], "#1.#2[1]";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(DynamoDB.value_or_literal(1), ":1";
                 attrs = Dict(), vals = Dict(":1" => Dict("N" => "1")))

check_expression(size(attr("foo")), "size(#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())

check_expression(begins_with(attr("foo"), "bar"), "begins_with(#1, :2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("S" => "bar")))
check_expression(begins_with(attr("foo"), attr("bar")), "begins_with(#1, #2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(contains(attr("foo"), "bar"), "contains(#1, :2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("S" => "bar")))
check_expression(contains(attr("foo"), attr("bar")), "contains(#1, #2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(is_string(attr("foo")), "attribute_type(#1, \"S\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_string_set(attr("foo")), "attribute_type(#1, \"SS\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_real(attr("foo")), "attribute_type(#1, \"N\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_real_set(attr("foo")), "attribute_type(#1, \"NS\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_binary(attr("foo")), "attribute_type(#1, \"B\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_binary_set(attr("foo")), "attribute_type(#1, \"BS\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_bool(attr("foo")), "attribute_type(#1, \"BOOL\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_null(attr("foo")), "attribute_type(#1, \"NULL\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_list(attr("foo")), "attribute_type(#1, \"L\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_map(attr("foo")), "attribute_type(#1, \"M\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(is_document(attr("foo")), "attribute_type(#1, \"M\")";
                 attrs = Dict("#1" => "foo"), vals = Dict())

check_expression(exists(attr("foo")), "attribute_exists(#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(not_exists(attr("foo")), "attribute_not_exists(#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())

check_expression(attr("foo") < 1, "(#1) < (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 < attr("foo"), "(:1) < (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") < attr("bar"), "(#1) < (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(attr("foo") > 1, "(#1) > (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 > attr("foo"), "(:1) > (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") > attr("bar"), "(#1) > (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(attr("foo") <= 1, "(#1) <= (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 <= attr("foo"), "(:1) <= (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") <= attr("bar"), "(#1) <= (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(attr("foo") >= 1, "(#1) >= (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 >= attr("foo"), "(:1) >= (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") >= attr("bar"), "(#1) >= (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(attr("foo") == 1, "(#1) = (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(attr("foo") == WeakRef(1), "(#1) = (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 == attr("foo"), "(:1) = (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(WeakRef(1) == attr("foo"), "(:1) = (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") == attr("bar"), "(#1) = (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(attr("foo") != 1, "(#1) <> (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(1 != attr("foo"), "(:1) <> (#2)";
                 attrs = Dict("#2" => "foo"), vals = Dict(":1" => Dict("N" => "1")))
check_expression(attr("foo") != attr("bar"), "(#1) <> (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(and(attr("foo"), attr("bar")), "(#1) AND (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_expression(or(attr("foo"), attr("bar")), "(#1) OR (#2)";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())
check_expression(or(attr("foo") > 7, attr("bar") < 7), "((#1) > (:2)) OR ((#3) < (:2))";
                 attrs = Dict("#1" => "foo", "#3" => "bar"), vals = Dict(":2" => Dict("N" => "7")))
@test DynamoDB.can_write_expression(or(attr("foo"), no_conditions())) == false
@test DynamoDB.can_write_expression(or(no_conditions(), attr("foo"))) == false
@test DynamoDB.can_write_expression(or(attr("foo") > 7, no_conditions())) == false
@test DynamoDB.can_write_expression(or(no_conditions(), attr("foo") > 7)) == false
@test DynamoDB.can_write_expression(or(no_conditions(), no_conditions())) == false

check_expression(!attr("foo"), "NOT (#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(!(attr("foo") > 7), "NOT ((#1) > (:2))";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "7")))
check_expression(not(attr("foo")), "NOT (#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())
check_expression(not(attr("foo") > 7), "NOT ((#1) > (:2))";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "7")))

check_expression(between(attr("foo"), 1, 2), "#1 BETWEEN :2 AND :3";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1"), ":3" => Dict("N" => "2")))
check_expression(between(attr("foo"), attr("bar"), 2), "#1 BETWEEN #2 AND :3";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict(":3" => Dict("N" => "2")))
check_expression(between(attr("foo"), 1, attr("bar")), "#1 BETWEEN :2 AND #3";
                 attrs = Dict("#1" => "foo", "#3" => "bar"), vals = Dict(":2" => Dict("N" => "1")))
check_expression(between(attr("foo"), attr("bar"), attr("baz")), "#1 BETWEEN #2 AND #3";
                 attrs = Dict("#1" => "foo", "#2" => "bar", "#3" => "baz"), vals = Dict())

# compound expressions
check_expression(and(attr("foo"), size(attr("bar")) > 4), "(#1) AND ((size(#2)) > (:3))";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict(":3" => Dict("N" => "4")))
check_expression(and(size(attr("bar")) > 4, attr("foo")), "((size(#1)) > (:2)) AND (#3)";
                 attrs = Dict("#3" => "foo", "#1" => "bar"), vals = Dict(":2" => Dict("N" => "4")))
check_expression(or(attr("foo"), size(attr("bar")) > 4), "(#1) OR ((size(#2)) > (:3))";
                 attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict(":3" => Dict("N" => "4")))
check_expression(or(size(attr("bar")) > 4, attr("foo")), "((size(#1)) > (:2)) OR (#3)";
                 attrs = Dict("#3" => "foo", "#1" => "bar"), vals = Dict(":2" => Dict("N" => "4")))

# duplicate references are combined into one
check_expression(attr("foo") > attr("foo"), "(#1) > (#1)";
                 attrs = Dict("#1" => "foo"), vals = Dict())

# duplicate values are combined into one
check_expression(and(attr("foo") >= 1, attr("foo") <= 1), "((#1) >= (:2)) AND ((#1) <= (:2))";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))

# non-conditions are removed by and
check_expression(and(no_conditions(), attr("foo") != 1), "(#1) <> (:2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))

@test DynamoDB.can_write_expression(nothing) == false

#  _   _           _       _
# | | | |_ __   __| | __ _| |_ ___
# | | | | '_ \ / _` |/ _` | __/ _ \
# | |_| | |_) | (_| | (_| | ||  __/
#  \___/| .__/ \__,_|\__,_|\__\___|
#       |_|
#  _____                              _
# | ____|_  ___ __  _ __ ___  ___ ___(_) ___  _ __  ___
# |  _| \ \/ / '_ \| '__/ _ \/ __/ __| |/ _ \| '_ \/ __|
# | |___ >  <| |_) | | |  __/\__ \__ \ | (_) | | | \__ \
# |_____/_/\_\ .__/|_|  \___||___/___/_|\___/|_| |_|___/
#            |_|

function check_updates(expr :: Array, expected;
                       attrs = Dict(), vals = Dict())
    refs = DynamoDB.refs_tracker()
    @test expected == DynamoDB.serialize_updates(expr, refs)

    @test refs.attrs == attrs
    @test refs.vals == vals
end
check_updates(x :: DynamoDB.DynamoUpdateExpression, expected; attrs = Dict(), vals = Dict()) =
    check_updates([x], expected; attrs=attrs, vals=vals)

# maybe should be above -- this function is only for update expressions though
check_expression(get_or_else(attr("foo"), 1), "if_not_exists(#1, :2)";
                 attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "1")))

# lone operations
check_updates(append_to_list(attr("foo"), 3), "SET #1 = list_append(#1, :2)";
              attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "3")))
check_updates(append_to_list(attr("foo"), attr("bar")), "SET #1 = list_append(#1, #2)",
              attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_updates(assign(attr("foo"), 3), "SET #1 = :2";
              attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "3")))
check_updates(assign(attr("foo"), attr("bar")), "SET #1 = #2";
              attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_updates(add_to_set(attr("foo"), 3), "ADD #1 :2";
              attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "3")))
check_updates(add_to_set(attr("foo"), attr("bar")), "ADD #1 #2";
              attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_updates(remove_from_set(attr("foo"), 3), "DELETE #1 :2";
              attrs = Dict("#1" => "foo"), vals = Dict(":2" => Dict("N" => "3")))
check_updates(remove_from_set(attr("foo"), attr("bar")), "DELETE #1 #2";
              attrs = Dict("#1" => "foo", "#2" => "bar"), vals = Dict())

check_updates(delete(attr("foo")), "REMOVE #1";
              attrs = Dict("#1" => "foo"), vals = Dict())

# combo time!
check_updates([delete(attr("foo")), append_to_list(attr("bar"), 7), delete(attr("baz")), assign(attr("boo"), 11)],
              "SET #2 = list_append(#2, :3), #5 = :6 REMOVE #1, #4";
              attrs = Dict("#1" => "foo", "#2" => "bar", "#4" => "baz", "#5" => "boo"),
              vals = Dict(":3" => Dict("N" => "7"), ":6" => Dict("N" => "11")))


text:
export polar,
       eigenHermitian,
       factorize

import LinearAlgebra.qr
function qr(A::ITensor,
            Linds...;
            kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,qr")
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  QT,RT = qr(tensor(A),Lpos,Rpos)
  Q,R = ITensor(QT),ITensor(RT)
  q = commonindex(Q,R)
  settags!(Q,tags,q)
  settags!(R,tags,q)
  q = settags(q,tags)
  return Q,R,q
end

# TODO: allow custom tags in internal indices?
function polar(A::ITensor,
               Linds...;
               kwargs...)
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,PT = polar(tensor(A),Lpos,Rpos)
  U,P = ITensor(UT),ITensor(PT)
  u = commoninds(U,P)
  p = uniqueinds(P,U)
  replaceinds!(U,u,p')
  replaceinds!(P,u,p')
  return U,P,commoninds(U,P)
end

import LinearAlgebra.svd
"""
    svd(A::ITensor,
        leftind1::Index,
        leftind2::Index,
        ...
        ;kwargs...)

Singular value decomposition (SVD) of an ITensor A, computed
by treating the "left indices" provided collectively
as a row index, and the remaining "right indices" as a
column index (matricization of a tensor).

Whether the SVD performs a trunction depends on the keyword
arguments provided. The following keyword arguments are recognized:
* `maxdim` [Int]
* `mindim` [Int]
* `cutoff` [Float64]
* `absoluteCutoff` [Bool] Default value: false.
* `doRelCutoff` [Bool] Default value: true.
* `utags` [String] Default value: "Link,u".
* `vtags` [String] Default value: "Link,v".
* `fastSVD` [Bool] Defaut value: false.
"""
function svd(A::ITensor,
             Linds...;
             kwargs...)
  utags::TagSet = get(kwargs,:utags,"Link,u")
  vtags::TagSet = get(kwargs,:vtags,"Link,v")
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,ST,VT = svd(tensor(A),Lpos,Rpos;kwargs...)
  U,S,V = ITensor(UT),ITensor(ST),ITensor(VT)
  u₀ = commonindex(U,S)
  v₀ = commonindex(S,V)

  u = settags(u₀,utags)
  v = settags(u₀,vtags)

  U *= δ(dag(u₀),u)
  S = δ(dag(u₀),u)*S*δ(dag(v₀),v)
  V *= δ(dag(v₀),v)

  return U,S,V,u,v
end

function _factorize_center(A::ITensor,
                           Linds...;
                           kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  v = commonindex(S,V)
  for ss = 1:dim(u)
    S[ss,ss] = sqrt(S[ss,ss])
  end
  FU = settags(U*S,tags,v)
  FV = settags(S*V,tags,u)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_left_svd(A::ITensor,
                                  Linds...;
                                  kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  FU = settags(U,tags,u)
  FV = settags(S*V,tags,u)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_right_svd(A::ITensor,
                                   Linds...;
                                   kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  v = commonindex(S,V)
  FU = settags(U*S,tags,v)
  FV = settags(V,tags,v)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_left_eigen(A::ITensor,
                                    Linds...;
                                    kwargs...)
  Lis = commoninds(inds(A),IndexSet(Linds...))
  A² = A*prime(dag(A),Lis)
  FU,D = eigenHermitian(A²,Lis,prime(Lis); ispossemidef=true,
                                           kwargs...)
  FV = dag(FU)*A
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_right_eigen(A::ITensor,
                                     Linds...;
                                     kwargs...)
  Ris = uniqueinds(inds(A),IndexSet(Linds...))
  A² = A*prime(dag(A),Ris)
  FV,D = eigenHermitian(A²,Ris,prime(Ris); ispossemidef=true,
                                           kwargs...)
  FU = A*dag(FV)
  return FU,FV,commonindex(FU,FV)
end

import LinearAlgebra.factorize
"""
factorize(A::ITensor, Linds...; dir = "center", which_factorization = "svd", cutoff = 0.0)

Do a low rank factorization of A either using an SVD or an eigendecomposition of A'A or AA'.
"""
function factorize(A::ITensor,
                   Linds...;
                   kwargs...)
  dir::String = get(kwargs,:dir,"center")
  if dir == "center"
    return _factorize_center(A,Linds...;kwargs...)
  end
  which_factorization::String = get(kwargs,:which_factorization,"svd")
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  use_eigen = false
  if which_factorization == "eigen" || (which_factorization == "automatic" && cutoff > 1e-12)
    use_eigen = true
  end
  if dir == "fromleft"
    if use_eigen
      return _factorize_from_left_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_left_svd(A,Linds...;kwargs...)
    end
  elseif dir == "fromright"
    if use_eigen
      return _factorize_from_right_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_right_svd(A,Linds...;kwargs...)
    end
  end
  throw(ArgumentError("In factorize, no dir = $dir supported. Use center, fromleft or fromright."))
end

function eigenHermitian(A::ITensor,
                        Linds=findinds(A,("",0)),
                        Rinds=prime(IndexSet(Linds));
                        kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,eigen")
  lefttags::TagSet = get(kwargs,:lefttags,tags)
  righttags::TagSet = get(kwargs,:righttags,prime(tags))
  Lis = commoninds(inds(A),IndexSet(Linds))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,DT = eigenHermitian(tensor(A),Lpos,Rpos;kwargs...)
  U,D = ITensor(UT),ITensor(DT)
  u = commonindex(U,D)
  settags!(U,lefttags,u)
  settags!(D,lefttags,u)
  u = settags(u,lefttags)
  v = uniqueindex(D,U)
  D *= δ(v,settags(u,righttags))
  return U,D,u,v
end

import LinearAlgebra.eigen
function eigen(A::ITensor,
               Linds=findinds(A,("",0)),
               Rinds=prime(IndexSet(Linds));
               kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,eigen")
  Lis = commoninds(inds(A),IndexSet(Linds))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,DT = eigen(tensor(A),Lpos,Rpos;kwargs...)
  U,D = ITensor(UT),ITensor(DT)
  u = commonindex(U,D)
  settags!(U,tags,u)
  settags!(D,tags,u)
  u = settags(u,tags)
  v = uniqueindex(D,U)
  D *= δ(v,u')
  return U,D,u
end

ext: jl


text:
import Base.Order: Ordering, Forward, ord, lt

# merge sorted vectors vl and vr into v
# from indices lo to hi in v
function mergesorted!(v::AbstractVector,
                      lo::Int, hi::Int,
                      vl::AbstractVector,
                      lol::Int, hil::Int,
                      vr::AbstractVector,
                      lor::Int, hir::Int,
                      order::Ordering)
    c = lol
    p = lor
    nl = hil
    nr = hir
    i = lo
    @inbounds while c <= nl && p <= nr && i <= hi
        if lt(order, vr[p], vl[c])
            v[i] = vr[p]
            p = p+1
            i = i+1
        else
            v[i] = vl[c]
            c = c+1
            i = i+1
        end
    end
    @inbounds while p <= nr && i <= hi
        v[i] = vr[p]
        i = i+1
        p = p+1
    end
    @inbounds while c <= nl && i <= hi
        v[i] = vl[c]
        i = i+1
        c = c+1
    end
    v
end

function mergesorted!(v::AbstractVector, vl::AbstractVector,
                      vr::AbstractVector, order::Ordering)
    inds = indices(v,1)
    indsl = indices(vl,1)
    indsr = indices(vr,1)
    mergesorted!(v,first(inds),last(inds),vl,first(indsl),last(indsl),
                 vr,first(indsr),last(indsr),order)
end

function mergesorted!(v::AbstractVector,
                      vl::AbstractVector,
                      vr::AbstractVector,
                      lt=isless,
                      by=identity,
                      rev::Bool=false,
                      order::Ordering=Forward)
    ordr = ord(lt,by,rev,order)
    mergesorted!(v, vl, vr, ordr)
end

function mergesorted(vl::AbstractVector,
                      vr::AbstractVector,
                      lt=isless,
                      by=identity,
                      rev::Bool=false,
                      order::Ordering=Forward)
    v = similar(promote_type(typeof(vl),typeof(vr)), length(vl)+length(vr))
    ordr = ord(lt,by,rev,order)
    mergesorted!(v, vl, vr, ordr)
end


text:
""" General-purpose tools for metaprogramming in Julia.
"""
module Meta
export Expr0, JuliaFunction, JuliaFunctionSig, parse_docstring, parse_function,
  parse_function_sig, generate_docstring, generate_function,
  append_expr!, concat_expr, replace_symbols, strip_lines

using Base.Meta: ParseError
using Compat
using AutoHashEquals
using MLStyle: @match

# Data types
############

const Expr0 = Union{Symbol,Expr}

@auto_hash_equals struct JuliaFunction
  call_expr::Expr
  return_type::Union{Expr0,Nothing}
  impl::Union{Expr,Nothing}
  doc::Union{String,Nothing}

  function JuliaFunction(call_expr::Expr, return_type=nothing,
                         impl=nothing, doc=nothing)
    new(call_expr, return_type, impl, doc)
  end
end

@auto_hash_equals struct JuliaFunctionSig
  name::Symbol
  types::Vector{Expr0}
end

# Parsing Julia functions
#########################

""" Parse Julia expression that is (possibly) annotated with docstring.
"""
function parse_docstring(expr::Expr)::Tuple{Union{String,Nothing},Expr}
  expr = strip_lines(expr)
  if expr.head == :macrocall && (
      # XXX: It seems that the @doc macro can show up in two forms.
      expr.args[1] == GlobalRef(Core, Symbol("@doc")) ||
      expr.args[1] == Expr(:core, Symbol("@doc")))
    (expr.args[2], expr.args[3])
  else
    (nothing, expr)
  end
end

""" Parse Julia function definition into standardized form.
"""
function parse_function(expr::Expr)::JuliaFunction
  doc, expr = parse_docstring(expr)
  fun_expr, impl = @match expr begin
    Expr(:(=), args...) => args
    Expr(:function, args...) => args
    _ => throw(ParseError("Ill-formed function definition $expr"))
  end
  @match fun_expr begin
    (Expr(:(::), Expr(:call, args...), return_type) =>
      JuliaFunction(Expr(:call, args...), return_type, impl, doc))
    (Expr(:call, args...) =>
      JuliaFunction(Expr(:call, args...), nothing, impl, doc))
    _ => throw(ParseError("Ill-formed function header $fun_expr"))
  end
end

""" Parse signature of Julia function.
"""
function parse_function_sig(call_expr::Expr)::JuliaFunctionSig
  name, args = @match call_expr begin
    Expr(:call, name::Symbol, Expr(:parameters, kw...), args...) => (name, args)
    Expr(:call, name::Symbol, args...) => (name, args)
    _ => throw(ParseError("Ill-formed function signature $call_expr"))
  end
  types = map(args) do expr
    @match expr begin
      Expr(:(::), _, type) => type
      Expr(:(::), type) => type
      _ => :Any
    end
  end
  JuliaFunctionSig(name, types)
end
parse_function_sig(fun::JuliaFunction) = parse_function_sig(fun.call_expr)

""" Wrap Julia expression with docstring.
"""
function generate_docstring(expr::Expr, doc::Union{String,Nothing})::Expr
  if isnothing(doc)
    expr
  else
    Expr(:macrocall, GlobalRef(Core, Symbol("@doc")),
         LineNumberNode(0), doc, expr)
  end
end

""" Generate Julia expression for function definition.
"""
function generate_function(fun::JuliaFunction)::Expr
  head = if isnothing(fun.return_type)
    fun.call_expr
  else
    Expr(:(::), fun.call_expr, fun.return_type)
  end
  body = if isnothing(fun.impl)
    Expr(:block)
  else
    # Wrap implementation inside block if not already.
    impl = fun.impl
    impl.head == :block ? impl : Expr(:block, impl)
  end
  expr = Expr(:function, head, body)
  generate_docstring(expr, fun.doc)
end

# Operations on Julia expressions
#################################

""" Append a Julia expression to a block expression.
"""
function append_expr!(block::Expr, expr)::Expr
  @assert block.head == :block
  @match expr begin
    Expr(:block, args...) => append!(block.args, args)
    _ => push!(block.args, expr)
  end
  block
end

""" Concatenate two Julia expressions into a block expression.
"""
function concat_expr(expr1::Expr, expr2::Expr)::Expr
  @match (expr1, expr2) begin
    (Expr(:block, a1...), Expr(:block, a2...)) => Expr(:block, a1..., a2...)
    (Expr(:block, a1...), _) => Expr(:block, a1..., expr2)
    (_, Expr(:block, a2...)) => Expr(:block, expr1, a2...)
    _ => Expr(:block, expr1, expr2)
  end
end

""" Replace symbols occurring anywhere in a Julia function (except the name).
"""
function replace_symbols(bindings::AbstractDict, f::JuliaFunction)::JuliaFunction
  JuliaFunction(
    Expr(f.call_expr.head, f.call_expr.args[1],
         (replace_symbols(bindings, a) for a in f.call_expr.args[2:end])...),
    isnothing(f.return_type) ? nothing : replace_symbols(bindings, f.return_type),
    isnothing(f.impl) ? nothing : replace_symbols(bindings, f.impl),
    f.doc
  )
end

""" Replace symbols occuring anywhere in a Julia expression.
"""
function replace_symbols(bindings::AbstractDict, expr)
  recurse(expr) = replace_symbols(bindings, expr)
  @match expr begin
    Expr(head, args...) => Expr(head, map(recurse,args)...)
    sym::Symbol => get(bindings, sym, sym)
    _ => expr
  end
end

""" Remove all LineNumberNodes from a Julia expression.
"""
function strip_lines(expr::Expr; recurse::Bool=false)::Expr
  args = [ x for x in expr.args if !isa(x, LineNumberNode) ]
  if recurse
    args = [ isa(x, Expr) ? strip_lines(x; recurse=true) : x for x in args ]
  end
  Expr(expr.head, args...)
end

end


text:
# using AmbulanceDeployment, Dates
#
# Pkg.add("Query")
# using Query
using AmbulanceDeployment

struct EMSEngine{T}
    eventlog::DataFrame
    eventqueue::PriorityQueue{T,Int,Base.Order.ForwardOrdering}
    shortqueue::PriorityQueue{T,Int,Base.Order.ForwardOrdering}
    guiArray::Array{Any, 1}
end

struct gui_event
    #call made / call responded / ambulance arrived / shortfall
    event_type::String
    #neighborhood where call is made from
    neighborhood_id::Int
    # station where call is made from
    deployment_id::Int
    # id of the event generated
    event_id::Int
    ambulance_id::Int
    arrival_time::Int
    remaining_amb::Int
    # current timestamp of event
    timestamp::DateTime
end

## generates the priority queue and instantiates the EMSEngine struct
function generateCalls(problem::DispatchProblem)
    ncalls = nrow(problem.emergency_calls)
    eventlog = DataFrame(
        id = 1:ncalls,
        dispatch_from = zeros(Int, ncalls),
        waittime = fill(0.0, ncalls),
        responsetime = fill(Inf, ncalls),
        scenetime = fill(Inf, ncalls),
        conveytime = fill(Inf, ncalls),
        returntime = fill(Inf, ncalls),
        return_to = zeros(Int, ncalls),
        return_type = fill(:station, ncalls),
        hospital = zeros(Int, ncalls),
        ambulance = zeros(Int, ncalls)
    )
    eventqueue = PriorityQueue{Tuple{Symbol,Int,Int,Int},Int}()
    shortqueue = PriorityQueue{Tuple{Symbol,Int,Int,Int},Int}()
    guiArray = Any[]
    for i in 1:nrow(problem.emergency_calls)
        t = problem.emergency_calls[i, :arrival_seconds]
        enqueue!(eventqueue, (:call, i, t, problem.emergency_calls[i, :neighborhood]), t)
    end
    EMSEngine{Tuple{Symbol,Int,Int,Int}}(eventlog, eventqueue, shortqueue, guiArray) ## returns the ems engine struct
end

function simulate_events!(problem::DispatchProblem, dispatch::ClosestDispatch)
    engine = generateCalls(problem)
    while(~isempty(engine.eventqueue))
        (event, id, t, value) = dequeue!(engine.eventqueue)

        if event == :call
            # if the event is reachable then change event to station2call
            # if not it is a shortfall
            println("calling event id: $id time: $t value: $value")
            call_event!(engine, problem, dispatch, id, t, value)

        elseif event == :station2call
            #once the ambulance arrives at the station change event to call2hospital
            println("Ambulance has arrived on the scene for event $id at time $t")
            station2call_event!(engine, problem, dispatch, id, t, value)
        elseif event == :call2hospital
            #once the ambulance arrives at the hospital it will hand over the patient and then
            #change the deploy back to its station
            println("Ambulance for event $id has arrived at the hospital at time $t")
            call2hosp_event!(engine, problem, dispatch, id, t, value)

        elseif event == :hospital2station
            #once the ambulance drops off the patient it will reroute to the station
            println("Ambulance for event $id is returning to the station at $t")
            return_event!(engine, problem, dispatch, id, t, value)

        elseif event == :done
            #increase the number of ambulances at that station now that it is done
            done_event!(engine, problem, dispatch, id, t, value)

        end
    end
    println("This is the Short Queue")
    while(~isempty(engine.shortqueue))
       (event, id, t, value) = dequeue!(engine.shortqueue)
        print("$event, $id, $t, $value")
        print("This is not empty")
    end
    engine.eventlog , engine.guiArray
end

function call_event!(
        ems::EMSEngine,
        problem::DispatchProblem,
        dispatch::DispatchModel,
        id::Int, # the id of the emergency call
        t::Int, # the time of the emergency call
        nbhd::Int; # the neighborhood the call is from
    )

    #check if there is an ambulance within the coverage matrix
    if sum(problem.deployment[problem.coverage[nbhd,:]]) == 0
        @assert false "$id: no ambulance reachable for call at $nbhd"
    #check if one of the ambulances within the coverage is available
    elseif sum(problem.available[problem.coverage[nbhd,:]]) > 0
        i = closest_available(dispatch, id, problem)
        @assert i > 0 "$id: dispatch from $i to nbhd $nbhd" # assume valid i (enforced by <if> condition)
        update_ambulances!(dispatch, i, -1)
        ems.eventlog[id, :dispatch_from] = i
        @assert problem.available[i] >= 0

        travel_time = ceil(Int, 60*problem.emergency_calls[id, Symbol("stn$(i)_min")])
        @assert travel_time >= 0
        if(t != problem.emergency_calls[id, :arrival_seconds])
            print("This is a shortfall")
        end
        ems.eventlog[id, :responsetime] = (travel_time + t - problem.emergency_calls[id,:arrival_seconds]) / 60 # minutes

        amb = respond_to!(dispatch, i, t)
        ems.eventlog[id, :ambulance] = amb

        event = gui_event("call responded", problem.emergency_calls[id, :neighborhood], i, id, amb,-1, problem.available[i] ,Dates.now())
        num_ambulances_array = Array{Integer}(undef, size(problem.available, 1))
        for i in 1:size(problem.available, 1)
            num_ambulances_array[i] = problem.available[i]
        end
        push!(ems.guiArray,event)
        push!(ems.guiArray,num_ambulances_array)
        enqueue!(ems.eventqueue, (:station2call, id, t + travel_time, amb), t + travel_time)
    #else queue it
    else
        event = gui_event("shortfall", problem.emergency_calls[id, :neighborhood], -1, id, -1,-1,-1 ,Dates.now())
        push!(ems.guiArray, event)
        enqueue!(ems.shortqueue, (:call, id, t, problem.emergency_calls[id, :neighborhood]), t)
        println(id, ": call from ", nbhd, " is a shortfall")
    end
end

function station2call_event!(
        ems::EMSEngine,
        problem::DispatchProblem,
        dispatch::DispatchModel,
        id::Int, # the id of the emergency call
        t::Int, # the time of the emergency call
        amb::Int
    )
    arriveatscene!(dispatch, amb, t)
    # time the ambulance spends at the scene
    scene_time = ceil(Int,60*0.4*rand(problem.turnaround)) # 60sec*0.4*mean(40minutes) ~ 15minutes
    ems.eventlog[id, :scenetime] = scene_time / 60 # minutes
    @assert scene_time > 0

    event = gui_event("ambulance arrived", problem.emergency_calls[id, :neighborhood], -1, id, amb, ceil(ems.eventlog[id, :responsetime]*60) ,-1 ,Dates.now())
    push!(ems.guiArray, event)
    enqueue!(ems.eventqueue, (:call2hospital, id, t + scene_time, amb), t + scene_time)
end

function call2hosp_event!(
        ems::EMSEngine,
        problem::DispatchProblem,
        dispatch::DispatchModel,
        id::Int, # the id of the emergency call
        t::Int, # the time of the emergency call
        amb::Int
    )
    h = let mintime = Inf, minindex = 0
        for h in 1:nrow(problem.hospitals)
            traveltime = problem.emergency_calls[id, Symbol("hosp$(h)_min")]
            if !isna(traveltime) && !(traveltime=="NA") && traveltime < mintime
                @assert traveltime >= 0
                minindex = h; mintime = traveltime
            end
        end
        minindex
    end
    @assert h != 0
    dispatch.hospital[amb] = ems.eventlog[id, :hospital] = h
    conveytime = 60*15 + ceil(Int, 60*problem.emergency_calls[id, Symbol("hosp$(h)_min")]) # ~20minutes
    ems.eventlog[id, :conveytime] = conveytime / 60 # minutes
    @assert conveytime >= 0 conveytime
    arriveathospital!(dispatch, amb, h, t)
    event = gui_event("ambulance hospital", problem.emergency_calls[id, :neighborhood], -1, id, amb, t,-1 ,Dates.now())
    push!(ems.guiArray, event)

    enqueue!(ems.eventqueue, (:hospital2station, id, t+conveytime, amb), t+conveytime)
end

function return_event!(
        ems::EMSEngine,
        problem::DispatchProblem,
        dispatch::DispatchModel,
        id::Int,
        t::Int,
        amb::Int
    )
    stn = returning_to!(dispatch, amb, t)
    h = dispatch.hospital[amb]
    returntime = ceil(Int,60*2*problem.hospitals[h, Symbol("stn$(stn)_min")]) # ~ 10minutes
    dispatch.hospital[amb] = 0
    ems.eventlog[id, :returntime] = returntime / 60 # minutes
    @assert returntime >= 0 returntime
    t_end = t + returntime
    event = gui_event("ambulance returning", problem.emergency_calls[id, :neighborhood], -1, id, amb, t,-1 ,Dates.now())
    push!(ems.guiArray, event)

    enqueue!(ems.eventqueue, (:done, id, t_end, amb), t_end)
end

function done_event!(
        ems::EMSEngine,
        problem::DispatchProblem,
        dispatch::DispatchModel,
        id::Int,
        t::Int,
        amb::Int
    )
    i = dispatch.assignment[amb]
    println("Ambulance from event $id has returned to station $i")
    @assert i > 0 "$id: dispatch from $i to nbhd $nbhd" # assume valid i (enforced by <if> condition)
    update_ambulances!(dispatch, i, 1)
    push!(dispatch.ambulances[i], amb)
    dispatch.status[amb] = :available
    ems.eventlog[id, :return_to] = id
    @assert problem.available[i] > 0
    event = gui_event("ambulance returned", problem.emergency_calls[id, :neighborhood], -1, id, amb, t,-1 ,Dates.now())
    push!(ems.guiArray, event)

    while(~isempty(ems.shortqueue))
        # record time difference
        (event, short_id, short_t, value) = dequeue!(ems.shortqueue)
        enqueue!(ems.eventqueue, (event, short_id, t , value), short_t)
    end
end

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
import Base.length
import JSON
export GRN, reset!, set_input!, get_output, step!, distance

struct GRN
    nin::Int64
    nout::Int64
    nreg::Int64
    ids::Array{Float64}
    enh::Array{Float64}
    inh::Array{Float64}
    beta::Float64
    delta::Float64
    weights::Array{Float64}
    cons::Array{Float64}
end

function get_weights(ids::Array{Float64}, enh::Array{Float64}, inh::Array{Float64},
                     beta::Float64)
    glen = length(ids)
    weights = zeros(length(ids), length(ids))
    for i in eachindex(ids)
        for j in eachindex(ids)
            # influence of i on j (based on excitatory and inhibitory proteins of i)
            weights[i, j] = (exp(-beta * abs(enh[i] - ids[j])) -
                             exp(-beta * abs(inh[i] - ids[j])))
        end
    end
    weights
end

function GRN(nin::Int64, nout::Int64, nreg::Int64, config::Config)
    grn_size = nin + nout + nreg
    ids = rand(grn_size)
    enh = rand(grn_size)
    inh = rand(grn_size)
    beta = rand() * (config.beta_max - config.beta_min) + config.beta_min
    delta = rand() * (config.delta_max - config.delta_min) + config.delta_min
    weights = get_weights(ids, enh, inh, beta)
    cons = ones(grn_size) ./ grn_size
    GRN(nin, nout, nreg, ids, enh, inh, beta, delta, weights, cons)
end

function GRN(nin::Int64, nout::Int64, nreg::Int64, ids::Array{Float64},
             enh::Array{Float64}, inh::Array{Float64},
             beta::Float64, delta::Float64)
    grn_size = nin + nout + nreg
    weights = get_weights(ids, enh, inh, beta)
    cons = ones(grn_size) ./ grn_size
    GRN(nin, nout, nreg, ids, enh, inh, beta, delta, weights, cons)
end

function GRN(grn::GRN)
    grn_size = grn.nin + grn.nout + grn.nreg
    cons = ones(grn_size) ./ grn_size
    weights = get_weights(grn.ids, grn.enh, grn.inh, grn.beta)
    GRN(grn.nin, grn.nout, grn.nreg, copy(grn.ids), copy(grn.enh),
        copy(grn.inh), grn.beta, grn.delta, weights, cons)
end

function GRN(json_repr::String)
    d = JSON.parse(json_repr)
    ids = Array{Float64}(d["ids"])
    enh = Array{Float64}(d["enh"])
    inh = Array{Float64}(d["inh"])
    GRN(d["nin"], d["nout"], d["nreg"], ids, enh, inh, d["beta"], d["delta"])
end

function length(grn::GRN)
    length(grn.ids)
end

function reset!(grn::GRN)
    grn.cons .= 1./length(grn.ids)
end

function set_input!(grn::GRN, input::Array{Float64})
    @assert length(input) == grn.nin
    # @assert all(input .>= 0.0)
    # @assert all(input .<= 1.0)
    for i in eachindex(input)
        grn.cons[i] = input[i]
    end
end

function get_output(grn::GRN)
    grn.cons[grn.nin+(1:grn.nout)]
end

function step!(grn::GRN)
    reg = deepcopy(grn.cons)
    reg[grn.nin+(1:grn.nout)] .= 0.0
    cons = max.(0.0, grn.delta / length(grn.ids) * (reg' * grn.weights))
    sumcons = sum(cons[(grn.nin+1):end])
    if sumcons > 0
        cons = cons ./ sumcons
    end
    for i in (grn.nin+1):length(grn.ids)
        grn.cons[i] = cons[i]
    end
end

function step!(grn::GRN, nsteps::Int64)
    for i in 1:nsteps
        step!(grn)
    end
end

function step!(grn::GRN, inputs::Array{Float64})
    set_input!(grn, inputs)
    step!(grn)
    get_output(grn)
end

function distance(g1::Array{Float64}, g2::Array{Float64})
    # TODO: align regulatory
    dist = 1.0
    if (length(g1) > 0) || (length(g2) > 0)
        if length(g1) > length(g2)
            dist = sum(abs.(g1 - [g2; zeros(length(g1)-length(g2))]))
        elseif length(g1) < length(g2)
            dist = sum(abs.(g2 - [g1; zeros(length(g2)-length(g1))]))
        else
            dist = sum(abs.(g1 .- g2))
        end
    end
    dist
end

function distance(g1::GRN, g2::GRN, config::Config)
    (distance(g1.ids, g2.ids) * config.dist_coef_ids +
     distance(g1.enh, g2.enh) * config.dist_coef_enh +
     distance(g1.inh, g2.inh) * config.dist_coef_inh +
     (abs(g1.beta - g2.beta) / (config.beta_max - config.beta_min)) +
     (abs(g1.delta - g2.delta) / (config.delta_max - config.delta_min)))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Foo
include("SubFoo1.jl")
include(joinpath("subdir", "SubFoo2.jl"))
import Bar, Baz, Qux

# This tests pkgdir and pathof when executed in toplevel
# on both the package itself and on its dependencies
@assert isdir(pkgdir(Foo))
@assert isdir(pkgdir(Bar))
@assert isfile(pathof(Foo))
@assert isfile(pathof(Bar))

this = "Foo1"
which = "path"
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Test
using Unfold
using StatsModels
using DataFrames
using Statistics
using Random

include("test_utilities.jl")

data, evts = loadtestdata("test_case_3a") #

data_r = reshape(data, (1, :))
data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (0, 0.05), sfreq = 10) # cut the data into epochs

##
f = @formula 0 ~ 1 + conditionA + continuousA # 1
m_mul = fit(Unfold.UnfoldModel, Dict(Any=>(f,times)), evts, data_e)

# test simple case
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_mul)
@test size(eff,1) == 2 # we specified 2 levels @ 1 time point
@test eff.conditionA ≈ [0.,1.] # we want to different levels
@test eff.yhat ≈ [2.0,5.0] # these are the perfect predicted values

# combination 2 levels /  6 values
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-2,0,2]),m_mul)
@test size(eff,1) == 6 # we want 6 values
@test eff.conditionA ≈ [0.,0.,0.,1.,1.,1.]
@test eff.continuousA ≈ [-2,0,2,-2,0,2.]

# testing typical value
eff_man = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[mean(evts.continuousA)]),m_mul)
eff_typ = Unfold.effects(Dict(:conditionA => [0,1]),m_mul)
@test eff_man.yhat ≈ eff_typ.yhat

## Testing Splines
f_spl = @formula 0 ~ 1 + conditionA + spl(continuousA, 4) # 1
m_mul_spl = fit(UnfoldModel, f_spl, evts, data_e, times)

eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_mul_spl)
@test size(eff,1) == 2 # we specified 2 levels @ 1 time point
@test eff.conditionA ≈ [0.,1.] # we want to different levels
@test eff.yhat ≈ [2.0,5.0] # these are the perfect predicted values

# combination 2 levels /  6 values
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-0.5,0,0.5]),m_mul_spl)
@test size(eff,1) == 6 # we want 6 values
@test eff.conditionA ≈ [0.,0.,0.,1.,1.,1.]
@test eff.continuousA ≈ [-0.5,0,0.5,-0.5,0,0.5]
@test eff.yhat ≈ [0,2,4,3,5,7]

# testing for safe predictions
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[2]),m_mul_spl)
@test all(ismissing.(eff.yhat ))

## Timeexpansion
data, evts = loadtestdata("test_case_3a") #
f = @formula 0 ~ 1 + conditionA + continuousA # 1
uf = fit(Unfold.UnfoldModel, Dict(Any=>(f,firbasis([0,0.1],10))), evts, data)
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),uf)
@test nrow(eff) == 4
@test eff.yhat ≈ [2., 2., 5., 5.]
@test eff.conditionA ≈ [0.,0.,1.,1.]
@test eff.continuousA ≈ [0.,0.,0.,0.]
## Two event Timeexpansion
data, evts = loadtestdata("test_case_4a") #
b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA")
b2 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisB")
f = @formula 0 ~ 1 # 1
m_tul = fit(UnfoldModel, Dict("eventA"=>(f,b1),"eventB"=>(f,b2)), evts, data,eventcolumn="type")

#p = predict(m_tul,DataFrame(:Cond => [0,1]))

eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_tul)
@test unique(eff.basisname)==["basisA","basisB"]
@test unique(eff.yhat) ≈ [2,3]
@test size(eff,1) == 2*2*20 # 2 basisfunctions, 2x conditionA, 1s a 20Hz

eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-1,0,1]),m_tul)
@test size(eff,1) == 2*6*20

## Different sized events + different Formulas
data, evts = loadtestdata("test_case_4a") #
evts[!,:continuousA] = rand(MersenneTwister(42),nrow(evts))
b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA")
b2 = firbasis(τ = (0.0, 0.5), sfreq = 20, name = "basisB")
f1 = @formula 0 ~ 1 # 1
f2 = @formula 0 ~ 1+continuousA # 1
m_tul = fit(UnfoldModel, Dict("eventA"=>(f1,b1),"eventB"=>(f2,b2)), evts, data,eventcolumn="type")
eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-1,0,1]),m_tul)
@test nrow(eff) == (length(b1.times)+length(b2.times))*6
@test sum(eff.basisname .== "basisA") == 120
@test sum(eff.basisname .== "basisB") == 66

## Test two channels
data, evts = loadtestdata("test_case_3a") #

data_r = repeat(reshape(data, (1, :)),3,1)
data_r[2,:] = data_r[2,:] .* 2

data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (0, 0.05), sfreq = 10) # cut the data into epochs

#
f = @formula 0 ~ 1 + conditionA + continuousA # 1
m_mul = fit(Unfold.UnfoldModel, Dict(Any=>(f,times)), evts, data_e)
m_tul = fit(Unfold.UnfoldModel, Dict(Any=>(f,firbasis([0,.05],10))), evts, data_r)
# test simple case
eff_m = Unfold.effects(Dict(:conditionA => [0,1,0,1],:continuousA =>[0]),m_mul)
eff_t = Unfold.effects(Dict(:conditionA => [0,1,0,1],:continuousA =>[0]),m_tul)

@test eff_m.yhat ≈ eff_t.yhat
@test length(unique(eff_m.channel)) == 3
@test eff_m[eff_m.channel .==1,:yhat] ≈ eff_m[eff_m.channel .==2,:yhat]./2
@test eff_m[eff_m.channel .==1,:yhat] ≈ [2,5,2,5.] # these are the perfect predicted values - note that we requested them twice


text:
@doc raw"""
    gauss_quadrature(q)

Return the Gauss-Legendre quadrature rule with `q` points (integrates polynomials of degree
$2q-1$ exactly).

A tuple `(x,w)` is returned.
"""
function gauss_quadrature(q)
    x = zeros(CeedScalar, q)
    w = zeros(CeedScalar, q)
    C.CeedGaussQuadrature(q, x, w)
    x, w
end

struct QuadratureMode{T} end
const Abscissa = QuadratureMode{:Abscissa}()
const AbscissaAndWeights = QuadratureMode{:AbscissaAndWeights}()

@doc raw"""
    lobatto_quadrature(q, mode::Mode=Abscissa)

Return the Gauss-Lobatto quadrature rule with `q` points (integrates polynomials of degree
$2q-3$ exactly).

If `mode` is `AbscissaAndWeights`, then both the weights and abscissa are returned as a
tuple `(x,w)`.

Otherwise, (if `mode` is `Abscissa`), then only the abscissa `x` are returned.
"""
function lobatto_quadrature(q, mode::Mode=Abscissa) where {Mode}
    return_weights = (mode == AbscissaAndWeights)
    x = zeros(CeedScalar, q)
    w = (return_weights) ? zeros(CeedScalar, q) : C_NULL
    C.CeedLobattoQuadrature(q, x, w)
    return_weights ? (x, w) : x
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
using Test
using Plots

include("../../../lib/DifferencesMatrix.jl")
using .DifferencesMatrix: backwardDiff, forwardDiff, centeredDiff

include("../../../lib/SpecialVectors.jl")
using .SpecialVectors: squares

include("../../../lib/KTBC.jl")
using .KTBC: CreateKTBC

# 8
dim = 5
forwardDiff(dim) - backwardDiff(dim)

backwardDiff(dim) * forwardDiff(dim)

backwardDiff(dim) * centeredDiff(dim)

# 10
forwardDiff(dim) * backwardDiff(dim)

function getH(n::Int)
    return 1 / (n + 1)
end

function getAnalogRange()
    return 0:10^(-2):1
end

# 13
function discretSolution(n::Int)
    h = getH(n)
    T = backwardDiff(n) * (-forwardDiff(n))
    u = h^2 * (T \ ones(n))
    return u
end

n = 7
analogSolution(x) = (1 / 2) * (1 - x^2)
steps = 0:1 / (n + 1):1

uₜ = analogSolution.([steps...])[2:(n + 1)]
plot(uₜ)

u₇ = discretSolution(n)
plot!(u₇)

e = uₜ - u₇
plot(e)

# 14
function cosAgainstDiscretCos(n::Int)
    x = getAnalogRange()

    h = getH(n)
    xᵢ = 0:h:1
    ω = 4 * π

    u = cos.(ω * x)
    uᵢ = cos.(ω * h * (1:n + 2))

    lineWidth = 2
    plot(x, u, label="analog cos", lw=lineWidth, title="Cos/DiscretCos Oscilations with N=$n")
    plot!(xᵢ, uᵢ, label="discret cos", lw=lineWidth)
end

cosAgainstDiscretCos(3)
cosAgainstDiscretCos(20)
cosAgainstDiscretCos(150)
cosAgainstDiscretCos(350)

# 16
function analogAndDiscretSolutions(n::Int)
    ω = 4 * π
    h = getH(n)

    rangeₐ = getAnalogRange()

    analogSolution(x) = ((ω)^-2) * cos(ω * x) - ((ω)^-2)
    u = analogSolution.(rangeₐ)

    rangeᵢ = h * (1:n)
    Kₙ, Tₙ, Bₙ, Cₙ = CreateKTBC(n)
    uᵢ = h^2 * (Kₙ \ cos.(ω * rangeᵢ))

    plot(rangeₐ, u, label="analog", title="Analog+Discret solutions with N=$n")
    plot!(rangeᵢ, uᵢ, label="discret")
end

analogAndDiscretSolutions(4)
analogAndDiscretSolutions(8)
analogAndDiscretSolutions(16)
analogAndDiscretSolutions(32)

# 18
function analogAndDiscretSolutions2(n::Int)
    h = getH(n)

    rangeₐ = getAnalogRange()

    analogSolution(x) = (x / 6) * (x^2 - 1)
    u = analogSolution.(rangeₐ)

    rangeᵢ = h * (1:n)
    Kₙ = CreateKTBC(n)[1]
    uᵢ = h^2 * (-Kₙ \ rangeᵢ)

    plot(rangeₐ, u, label="analog", title="Analog+Discret solutions with N=$n")
    plot!(rangeᵢ, uᵢ, label="discret")
end

analogAndDiscretSolutions2(4)
analogAndDiscretSolutions2(8)
analogAndDiscretSolutions2(32)
analogAndDiscretSolutions2(64)

# 19
function analogAndDiscretSolutions3(n::Int)
    h = getH(n)

    rangeₐ = getAnalogRange()

    cₐ = 1 / (ℯ - 1)
    analogSolution(x) = cₐ - cₐ * ℯ^x + x
    u = analogSolution.(rangeₐ)

    f = ones(n)
    ◬₊ = forwardDiff(n)
    ◬₋ = backwardDiff(n)
    ◬₀ = (◬₊ + ◬₋) / (2 * h)
    Kₙ = CreateKTBC(n)[1]
    K = Kₙ / (h^2)

    rangeᵢ = h * (1:n)
    uᵢ = round.((K + ◬₀) \ f, digits=4)
    u₊ = round.((K + ◬₊ / h) \ f, digits=4)

    plot(rangeₐ, u, label="analog", title="Analog+Discret solutions with N=$n")
    plot!(rangeᵢ, uᵢ, label="discret centered")
    plot!(rangeᵢ, u₊, label="discret uncentered")
end

analogAndDiscretSolutions3(4)
analogAndDiscretSolutions3(8)
analogAndDiscretSolutions3(16)
analogAndDiscretSolutions3(32)
analogAndDiscretSolutions3(64)


text:
abstract type Primitive end
abstract type Terminal <: Primitive end

struct PrimitiveFunction{T} <: Primitive
    func::T
    arity::Int
end

struct Variable{T} <: Terminal
    name::T
end

struct Const{T <: Real} <: Terminal
    value::T
end

struct RandConst <: Terminal
end

protected_div(x, y) = ifelse(y == 0, 1, x / y)
sigmoid(x) = 1 / (1 + exp(-x))

mutable struct PrimitiveSet
    functions::Vector{PrimitiveFunction}
    terminals::Vector{Terminal}

    terminal_ratio::Float64
    function_ratio::Float64

    PrimitiveSet() = new([], [], 0, 0)
end

arity(x::PrimitiveFunction) = x.arity
arity(x::Terminal) = 0

function update_ratio!(pset)
    function_ratio = length(pset.functions) / (length(pset.functions) + length(pset.terminals))
    pset.function_ratio = function_ratio
    pset.terminal_ratio = 1 - function_ratio
end

function add_function!(pset::PrimitiveSet, func::Function, arity::Int)
    prim = PrimitiveFunction(func, arity)
    push!(pset.functions, prim)
    update_ratio!(pset)
end

function add_const!(pset::PrimitiveSet, val)
    term = Const(val)
    push!(pset.terminals, term)
    update_ratio!(pset)
end

function add_rand!(pset)
    push!(pset.terminals, RandConst())
    update_ratio!(pset)
end

function add_variable!(pset, varname)
    var = Variable(varname)
    push!(pset.terminals, var)
    update_ratio!(pset)
end

struct GPNode{T <: Primitive}
    primitive::T
    children::Vector{GPNode}
end

GPNode(prim::PrimitiveFunction) = GPNode(prim, Vector{GPNode}(undef, prim.arity))
GPNode(prim::Terminal) = GPNode(prim, GPNode[])

dummy() = error("dummy node")
dummy_node() = GPNode(PrimitiveFunction(dummy, 0))

function Base.rand(x::Vector{Terminal})
    r = rand(1:length(x))
    if x[r] isa RandConst
        return Const(rand())
    else
        return x[r]
    end
end

function generate(pset, min_depth, max_depth, terminal_condition)
    local node
    if terminal_condition(pset, 0, min_depth, max_depth)
        prim = rand(pset.terminals)
        node = GPNode(prim)
    else
        prim = rand(pset.functions)
        node = GPNode(prim)
        _generate!(node, pset, 1, min_depth, max_depth, terminal_condition)
    end
    return node
end

function _generate!(node, pset, depth, min_depth, max_depth, terminal_condition)
    for i in 1:node.primitive.arity
        if terminal_condition(pset, depth, min_depth, max_depth)
            prim = rand(pset.terminals)
            node.children[i] = GPNode(prim)
        else
            prim = rand(pset.functions)
            node.children[i] = GPNode(prim)
        end
        if node.children[i].primitive isa PrimitiveFunction
            _generate!(node.children[i], pset, depth+1, min_depth, max_depth, terminal_condition)
        end
    end
end

cond_full(pset, depth, min_depth, max_depth) = depth == max_depth
cond_grow(pset, depth, min_depth, max_depth) = depth == max_depth || (depth >= min_depth && rand() < pset.terminal_ratio)

generate_full(pset, min_depth, max_depth) = generate(pset, max_depth, max_depth, cond_full)
generate_grow(pset, min_depth, max_depth) = generate(pset, min_depth, max_depth, cond_grow)

function ramped_half_and_half(cfg, data_train, data_test)
    entries = Vector{Entry}(undef, cfg.population_size)

    index = 1
    for i in 1:cfg.population_size÷2
        tree = generate_full(cfg.pset, cfg.min_depth, cfg.max_depth)
        entries[i] = Entry(tree, data_train, data_test, cfg, index=index)
        index += 1
    end

    for i in cfg.population_size÷2+1:cfg.population_size
        tree = generate_grow(cfg.pset, cfg.min_depth, cfg.max_depth)
        entries[i] = Entry(tree, data_train, data_test, cfg, index=index)
        index += 1
    end

    return entries
end

function evaluate(node::GPNode{PrimitiveFunction{T}}, vars) where T
    if node.primitive.arity == 2
        arg1 = evaluate(node.children[1], vars)
        arg2 = evaluate(node.children[2], vars)
        return node.primitive.func(arg1, arg2)
    elseif node.primitive.arity == 1
        arg = evaluate(node.children[1], vars)
        return node.primitive.func(arg)
    else
        error("functions with arity > 2 are not supported arity")
    end
end
evaluate(node::GPNode{Const{T}}, vars) where T <: Real =  node.primitive.value
evaluate(node::GPNode{Variable{T}}, vars) where T = vars[node.primitive.name]

# Interface to AbstractTrees.jl
function AbstractTrees.children(node::GPNode)
    if node.primitive isa PrimitiveFunction
        return node.children
    end
    return ()
end

function AbstractTrees.printnode(io::IO, node::GPNode)
    prim = node.primitive
    if prim isa PrimitiveFunction
        print(io, prim.func)
    elseif prim isa Variable
        print(io,prim.name)
    else
        print(io, prim.value)
    end
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
using Dates

abstract type Object <: Component end

include("../data/entities.jl")
struct Entity <: Object
    name::SubString{String}
end
Entity(name::String) = Entity(SubString(name))

struct LaTeXFragment <: Object
    contents::SubString{String}
    delimiters::Union{Tuple{SubString{String}, SubString{String}}, Nothing}
end
LaTeXFragment(contents::String, delimitors::Tuple{String,String}=("\\(", "\\)")) =
    LaTeXFragment(SubString(contents), SubString.(delimitors))

struct ExportSnippet <: Object
    backend::SubString{String}
    snippet::SubString{String}
end
ExportSnippet(backend::String, snippet::String) =
    ExportSnippet(SubString(backend), SubString(snippet))

mutable struct FootnoteReference{L <: Union{SubString{String}, Nothing},
                                 D <: Union{Vector{Object}, Nothing}} <: Object
    label::L
    definition::D
end
FootnoteReference(label::String, definition::Union{Vector{Object}, Nothing}=nothing) =
    FootnoteReference(SubString(label), definition)
FootnoteReference(::Nothing, definition::Vector{Object}) =
    FootnoteReference(nothing, definition)
FootnoteReference(::Nothing, definition::String) =
    FootnoteReference(nothing, Object[TextPlain(definition)])

mutable struct CitationReference <: Object
    prefix::Vector{Object}
    key::SubString{String}
    suffix::Vector{Object}
end
CitationReference(prefix::Vector{Object}, key::String, suffix::Vector{Object}=Object[]) =
    CitationReference(prefix, SubString(key), suffix)
CitationReference(prefix::String, key::String, suffix::String) =
    CitationReference(Object[TextPlain(prefix)], key, Object[TextPlain(suffix)])
CitationReference(prefix::String, key::String, ::Nothing=nothing) =
    CitationReference(Object[TextPlain(prefix)], key)
CitationReference(key::String, suffix::Vector{Object}=Object[]) =
    CitationReference(Object[], key, suffix)
CitationReference(::Nothing, key::String, suffix::String) =
    CitationReference(key, Object[TextPlain(suffix)])

mutable struct Citation <: Object
    style::Tuple{Union{SubString{String}, Nothing},
                 Union{SubString{String}, Nothing}}
    globalprefix::Vector{Object}
    citerefs::Vector{CitationReference}
    globalsuffix::Vector{Object}
end
Citation(style::String, gpref::Vector{Object},
         citerefs::Vector{CitationReference}, gsuf::Vector{Object}) =
    Citation((SubString(style), nothing), gpref, citerefs, gsuf)
Citation(style::Union{Tuple{String, Nothing}, Tuple{String}}, gpref::Vector{Object},
         citerefs::Vector{CitationReference}, gsuf::Vector{Object}) =
    Citation((SubString(style[1]), nothing), gpref, citerefs, gsuf)
Citation(style::Tuple{String, String}, gpref::Vector{Object},
         citerefs::Vector{CitationReference}, gsuf::Vector{Object}) =
    Citation((SubString(style[1]), SubString(style[2])), gpref, citerefs, gsuf)
Citation(style::Union{String, Tuple}, gpref::Vector{Object}, citerefs::Vector{CitationReference}) =
    Citation(style, gpref, citerefs, Object[])
Citation(gpref::Vector{Object}, citerefs::Vector{CitationReference}) =
    Citation((nothing, nothing), gpref, citerefs)
Citation(style::Union{String, Tuple}, citerefs::Vector{CitationReference}, gsuf::Vector{Object}=Object[]) =
    Citation(style, Object[], citerefs, gsuf)
Citation(citerefs::Vector{CitationReference}, gsuf::Vector{Object}=Object[]) =
    Citation((nothing, nothing), citerefs, gsuf)
Citation(style::Union{String, Tuple}, key::String) =
    Citation(style, Object[], CitationReference(key), Object[])
Citation(key::String) =
    Citation((nothing, nothing), CitationReference(key))

struct InlineBabelCall <: Object
    name::SubString{String}
    header::Union{SubString{String}, Nothing}
    arguments::Union{SubString{String}, Nothing}
end
InlineBabelCall(name::String, arguments::Union{SubString{String}, Nothing}=nothing) =
    InlineBabelCall(SubString(name), nothing, arguments)

struct InlineSourceBlock <: Object
    lang::SubString{String}
    options::Union{SubString{String}, Nothing}
    body::SubString{String}
end
InlineSourceBlock(lang::String, body::String) =
    InlineSourceBlock(SubString(lang), nothing, SubString(body))

struct LineBreak <: Object end

abstract type Link <: Object end

struct LinkPath <: Object
    protocol::Union{Symbol, SubString{String}}
    path::SubString{String}
end
LinkPath(protocol::Symbol, path::String) =
    LinkPath(protocol, SubString(path))
LinkPath(protocol::String, path::String) =
    LinkPath(SubString(protocol), SubString(path))

struct PlainLink <: Link
    path::LinkPath
end
PlainLink(protocol::Union{Symbol, String}, path::String) =
    PlainLink(LinkPath(protocol, path))

struct AngleLink <: Link
    path::LinkPath
end
AngleLink(protocol::Union{Symbol, String}, path::String) =
    AngleLink(LinkPath(protocol, path))

struct RegularLink <: Link
    path::LinkPath
    description::Union{Vector{Object}, Nothing}
end
RegularLink(protocol::Union{Symbol, String}, path::String, description::Union{Vector{Object}, Nothing}=nothing) =
    RegularLink(LinkPath(protocol, path), description)
RegularLink(protocol::Union{Symbol, String}, path::String, description::String) =
    RegularLink(LinkPath(protocol, path), Object[TextPlain(description)])

struct Macro <: Object
    name::SubString{String}
    arguments::Vector{SubString{String}}
end
Macro(name::String, arguments::Vector{String}=String[]) =
    Macro(SubString(name), SubString.(arguments))

mutable struct RadioTarget <: Object
    contents::Vector{Object}
end

mutable struct RadioLink <: Link
    radio::RadioTarget
end

mutable struct Target <: Object
    target::SubString{String}
end
Target(target::String) = Target(SubString(target))

abstract type StatisticsCookie <: Object end
struct StatisticsCookiePercent <: StatisticsCookie
    percentage::SubString{String}
end
struct StatisticsCookieFraction <: StatisticsCookie
    complete::Union{Integer, Nothing}
    total::Union{Integer, Nothing}
end
StatisticsCookie() = StatisticsCookieFraction(nothing, nothing)
StatisticsCookie(complete::Integer, total::Integer) =
    StatisticsCookieFraction(complete, total)
StatisticsCookie(percentage::String) =
    StatisticsCookiePercent(SubString(percentage))
StatisticsCookie(percentage::Float64) =
    StatisticsCookiePercent(SubString(string(round(Int, 100*percentage), '%')))

abstract type Script <: Object end
struct Subscript <: Script
    char::Char
    script::SubString{String}
end
Subscript(char::Char, script::String) = Subscript(char, SubString(script))
struct Superscript <: Script
    char::Char
    script::SubString{String}
end
Superscript(char::Char, script::String) = Superscript(char, SubString(script))

mutable struct TableCell <: Object
    contents::Vector{Object}
end
TableCell(text::String) = TableCell(Object[TextPlain(text)])

abstract type Timestamp <: Object end
mutable struct TimestampRepeaterOrDelay
    type::Symbol
    value::Real
    unit::Char
end
mutable struct TimestampDiary <: Timestamp
    sexp::SubString{String}
end
abstract type TimestampInstant <: Timestamp end
mutable struct TimestampActive <: TimestampInstant
    date::Date
    time::Union{Time, Nothing}
    repeater::Union{TimestampRepeaterOrDelay, Nothing}
    warning::Union{TimestampRepeaterOrDelay, Nothing}
    function TimestampActive(date::Date, time::Union{Time, Nothing}=nothing,
                             repeater::Union{TimestampRepeaterOrDelay, Nothing}=nothing,
                             warning::Union{TimestampRepeaterOrDelay, Nothing}=nothing)
        new(date, time, repeater, warning)
    end
end
TimestampActive(datetime::DateTime) =
    TimestampActive(Date(datetime), Time(datetime))
mutable struct TimestampInactive <: TimestampInstant
    date::Date
    time::Union{Time, Nothing}
    repeater::Union{TimestampRepeaterOrDelay, Nothing}
    warning::Union{TimestampRepeaterOrDelay, Nothing}
    function TimestampInactive(date::Date, time::Union{Time, Nothing}=nothing,
                               repeater::Union{TimestampRepeaterOrDelay, Nothing}=nothing,
                               warning::Union{TimestampRepeaterOrDelay, Nothing}=nothing)
        new(date, time, repeater, warning)
    end
end
TimestampInactive(datetime::DateTime) =
    TimestampInactive(Date(datetime), Time(datetime))
abstract type TimestampRange <: Timestamp end
mutable struct TimestampActiveRange <: TimestampRange
    start::TimestampActive
    stop::TimestampActive
end
mutable struct TimestampInactiveRange <: TimestampRange
    start::TimestampInactive
    stop::TimestampInactive
end

mutable struct TextMarkup{C <: Union{Vector{Object}, SubString{String}}} <: Object
    formatting::Symbol
    contents::C
end
function TextMarkup(formatting::Symbol, contents::String)
    if formatting in ('=', '~')
        TextMarkup(formatting, SubString(contents))
    else
        TextMarkup(formatting, TextPlain(contents))
    end
end

struct TextPlain{S <: Union{SubString{String}, String}} <: Object
    text::S
end


text:
"""
	TimeVaryingCBFQP <: Policy

Struct representing a Control Barrier Function based quadratic program (CBFQP).

# Fields
- `control`: function u = control(x) describing the control law.
"""
struct TimeVaryingCBFQP <: Policy
    control
end

"""
	(κ::TimeVaryingCBFQP)(x, t)

Solve CBF-QP at state x and time t.
"""
function (κ::TimeVaryingCBFQP)(x, t)
    return κ.control(x, t)
end

"""
    TimeVaryingCBFQP(Σ::ControlAffineSystem, CBF::ControlBarrierFunction, κ::TimeVaryingCLFQP)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, CBF::ControlBarrierFunction, κ::TimeVaryingCLFQP, A, b)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBF::HighOrderCBF, κ::TimeVaryingCLFQP)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBF::HighOrderCBF, κ::TimeVaryingCLFQP, A, b)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, CBFs::Vector{ControlBarrierFunction}, κ::TimeVaryingCLFQP)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, CBFs::Vector{ControlBarrierFunction}, κ::TimeVaryingCLFQP, A, b)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBFs::Vector{HighOrderCBF}, κ::TimeVaryingCLFQP)
    TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBFs::Vector{HighOrderCBF}, κ::TimeVaryingCLFQP, A, b)

Construct a TimeVaryingCBFQP.
"""
function TimeVaryingCBFQP(Σ::ControlAffineSystem, CBF::ControlBarrierFunction, κ::TimeVaryingCLFQP)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        @constraint(model, CBF.∇h(x) * (Σ.f(x) + Σ.g(x) * u) >= -CBF.α(CBF(x)))
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, CBF::ControlBarrierFunction, κ::TimeVaryingCLFQP, A, b)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        @constraint(model, CBF.∇h(x) * (Σ.f(x) + Σ.g(x) * u) >= -CBF.α(CBF(x)))
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBF::HighOrderCBF, κ::TimeVaryingCLFQP)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        @constraint(model, HOCBF.∇ψ(x) * (Σ.f(x) + Σ.g(x) * u) >= -HOCBF.α(HOCBF(x)))
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBF::HighOrderCBF, κ::TimeVaryingCLFQP, A, b)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        @constraint(model, HOCBF.∇ψ(x) * (Σ.f(x) + Σ.g(x) * u) >= -HOCBF.α(HOCBF(x)))
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, CBFs::Vector{ControlBarrierFunction}, κ::TimeVaryingCLFQP)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        for CBF in CBFs
            @constraint(model, CBF.∇h(x) * (Σ.f(x) + Σ.g(x) * u) >= -CBF.α(CBF(x)))
        end
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, CBFs::Vector{ControlBarrierFunction}, κ::TimeVaryingCLFQP, A, b)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        for CBF in CBFs
            @constraint(model, CBF.∇h(x) * (Σ.f(x) + Σ.g(x) * u) >= -CBF.α(CBF(x)))
        end
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBFs::Vector{HighOrderCBF}, κ::TimeVaryingCLFQP)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        for HOCBF in HOCBFs
            @constraint(model, HOCBF.∇ψ(x) * (Σ.f(x) + Σ.g(x) * u) >= -HOCBF.α(HOCBF(x)))
        end
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end

function TimeVaryingCBFQP(Σ::ControlAffineSystem, HOCBFs::Vector{HighOrderCBF}, κ::TimeVaryingCLFQP, A, b)
    function control(x, t)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u - κ(x,t)'u)
        for HOCBF in HOCBFs
            @constraint(model, HOCBF.∇ψ(x) * (Σ.f(x) + Σ.g(x) * u) >= -HOCBF.α(HOCBF(x)))
        end
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return TimeVaryingCBFQP(control)
end


text:
export Polygon, PolygonalPrism
export regpoly, isosceles

#= Polygon (for a base shape) =#

# Assume the followings for the polygon represented by Polygon:
# - The polygon is convex.
# - The vertices are listed in the counter-clockwise order around the origin.

mutable struct Polygon{K,K2,D,T} <: Shape{2,4,D,T}  # K2 = 2K
    v::SMatrix{K,2,T,K2}  # vertices
    n::SMatrix{K,2,T,K2}  # direction normals to edges
	l::SVector{2,T}		  # lower x,y bounds
	u::SVector{2,T}		  # upper x,y bounds
	sz::SVector{2,T}	  # "size" in x and y: sz = l-u
	rbnd::T				  # minimum distance before point is considered to be on boundary, rmin = √eps(T) * max(sz.data...)
    data::D  # auxiliary data
    # Polygon{K,K2,D,T}(v,n,data) where {K,K2,D,T<:Real} = new(v,n,data)  # suppress default outer constructor
end
# mutable struct Polygon{K,K2,D} <: Shape{2,4,D}  # K2 = 2K
#     v::SMatrix{K,2,Float64}  # vertices
#     n::SMatrix{K,2,Float64}  # direction normals to edges
#     data::D  # auxiliary data
#     Polygon{K,K2,D}(v,n,data) where {K,K2,D} = new(v,n,data)  # suppress default outer constructor
# end

function l_bnds_poly(v::SMatrix{K,2,T}) where {K,T<:Real}
	@tullio (min) res[b] :=  v[a,b]
	return SVector(res)
end
function u_bnds_poly(v::SMatrix{K,2,T}) where {K,T<:Real}
	@tullio (max) res[b] :=  v[a,b]
	return SVector(res)
end
# function Polygon(v::AbstractMatrix{T}, data::D=nothing) where {D,T<:Real}
#     K = size(v,1);

function ∆ϕ(∆v::SMatrix{K,2,<:Real}) where K
	# Calculate the increases in angle between neighboring edges.
	∆z = ∆v[:,1] + im * ∆v[:,2]  # SVector{K}: edge directions as complex numbers
	icurr = ntuple(identity, Val(K-1))
	inext = ntuple(x->x+1, Val(K-1))
	return angle.(∆z[SVector(inext)] ./ ∆z[SVector(icurr)])  # angle returns value between -π and π
end

function sort_v_if_needed(v::SMatrix{K,2,T}) where {K,T}
	w = v .- mean(v, dims=1)  # v in center-of-mass coordinates
	ϕ = mod.(atan.(w[:,2], w[:,1]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
	if !issorted(ϕ)	# TODO: make sort_verts shuffling fn with AD rules, currently unsorted verts would break differentiability
		# Do this only when ϕ is not sorted, because the following uses allocations.
		ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
		v = v[ind,:]  # SVector{K}: sorted v
	end
end

function n_norm(dv)
	local one_mone = [1, -1]
	@tullio n[i,j] := dv[i,3-j] * $one_mone[j] / sqrt( dv[i,1]^2 + dv[i,2]^2 ) nograd=one_mone
end

n_norm_fwd(dv) = Zygote.forwarddiff(n_norm,dv)

ChainRulesCore.@non_differentiable ∆ϕ(::Any)

function Polygon(v::SMatrix{K,2,T}, data::D=nothing) where {K,D,T<:Real}
    # Sort the vertices in the counter-clockwise direction

    # w = v .- mean(v, dims=1)  # v in center-of-mass coordinates
    # ϕ = mod.(atan.(w[:,2], w[:,1]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
    # if !issorted(ϕ)	# TODO: make sort_verts shuffling fn with AD rules, currently unsorted verts would break differentiability
    #     # Do this only when ϕ is not sorted, because the following uses allocations.
    #     ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
    #     v = v[ind,:]  # SVector{K}: sorted v
    # end
    ∆v = v - circshift(v,1)
    # Check all the angle increases are positive.  If they aren't, the polygon is not convex.
    # all(∆ϕ(∆v) .> 0) || throw("v = $v should represent vertices of convex polygon.")

	# outward normal directions to edges
	# one_mone = [	0. 1.
	# 				-1. 0. ]# SMatrix{2,2}(0., 1., -1., 0.)
    # n0 = ∆v * one_mone #[∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
	# @tullio nnorm[k] := ∆v[k,j]^2 |> sqrt
	# @tullio n[k,j] := n0[k,j] / nnorm[k] # n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize
	# n0 = [∆v[:,2] -∆v[:,1]]  # outward normal directions to edges
    # n = n0 ./ hypot.(n0[:,1],n0[:,2])  # normalize

	# n = n_norm(∆v)
	n = n_norm_fwd(∆v)
	l = l_bnds_poly(v)
	u = u_bnds_poly(v)
	sz = u-l
	rbnd = Zygote.@ignore(max(sz.data[1],sz.data[2])*Base.rtoldefault(T))
	# TODO: Ask why `bounds(::Polygon)` was prev. recomputed each time and why
	# the "size" (`sz`) in find_surfpt_nearby(Polygon was computed as abs.(l-u)?
	# Hopefully this doesn't break something or hurt performance
    return Polygon{K,2K,D,T}(v,SMatrix{K,2,T}(n),l,u,sz,rbnd,data) # Polygon{K,2K,D}(v,n,data)
end

function _∆xe_poly(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2})::SVector{K}  where {K} #,T<:Real}
	@tullio out[k] := n[k,a] * ( x[a] - v[k,a] ) # edge line eqns for a K-point Polygon{K} `s`
end
function ∆xe_poly(x::AbstractVector{<:Real},s::Polygon{K})::SVector{K} where K
	_∆xe_poly(x,s.v,s.n)
end
bounds(s::Polygon) = (s.l, s.u)

# function onbnd(Δxe::SVector{K},s::Polygon{K})::SVector{K,Bool} where K 		# fixed wrong Δ, keeping in case this breaks
function onbnd(∆xe::SVector{K},s::Polygon{K})::SVector{K,Bool} where K
	map(x->abs(x)≤s.rbnd,∆xe)
end
function onbnd(x::SVector{2},s::Polygon)::SVector{K,Bool} where K
	map(x->abs(x)≤s.rbnd,∆xe_poly(x,s))
end
function cout(x::SVector{2},s::Polygon)::Int
	mapreduce((a,b)->Int(a|b),+,onbnd(x,s),map(isposdef,∆xe_poly(x,s)))
end
function cout(∆xe::SVector{K,<:Real},obd::SVector{K,Bool})::Int where K
	mapreduce((a,b)->Int(a|b),+,obd,map(isposdef,∆xe))
end
function imin2(x::SVector{2},s::Polygon{K})::Tuple{Int,Int} where K
	∆xv = x' .- s.v
	l∆xv = hypot.(∆xv[:,1], ∆xv[:,2])
	imin = argmin(l∆xv)
	return  imin ,  mod1(imin-1,K) # imin,imin₋₁
end

function _∆x_poly(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2},kmax::Int)::SVector{2}  where {K} #,T<:Real}
	# @tullio Δx[i] := n[$kmax,a] * ( v[$kmax,a] - x[a] * n[$kmax,i]	# works but gradient doesn't vectorize
	( @tullio Δx_factor := n[$kmax,a] * ( v[$kmax,a] - x[a] ) )  * SVector(n[kmax,1],n[kmax,2])
end
∆x_poly(x::AbstractVector{<:Real},s::Polygon)::SVector{2} = _∆x_poly(x,s.v,s.n,argmax(_∆xe_poly(x,s.v,s.n)))

ChainRulesCore.@non_differentiable cout(::Any,::Any)
ChainRulesCore.@non_differentiable onbnd(::Any,::Any)
ChainRulesCore.@non_differentiable imin2(::Any,::Any)

# Polygon(v::AbstractMatrix{<:Real}, data=nothing) = (K = size(v,1); Polygon(SMatrix{K,2}(v), data))

Base.:(==)(s1::Polygon, s2::Polygon) = s1.v==s2.v && s1.n==s2.n && s1.data==s2.data  # assume sorted v
Base.isapprox(s1::Polygon, s2::Polygon) = s1.v≈s2.v && s1.n≈s2.n && s1.data==s2.data  # assume sorted v
Base.hash(s::Polygon, h::UInt) = hash(s.v, hash(s.n, hash(s.data, hash(:Polygon, h))))

Base.in(x::SVector{2,<:Real}, s::Polygon) = all(sum(s.n .* (x' .- s.v), dims=Val(2)) .≤ 0)
# Base.in(x::SVector{3,<:Real}, s::Polygon) = ( x2 = SVector(x[1], x[2]); all(sum(s.n .* (x2' .- s.v), dims=Val(2)) .≤ 0) )

function surfpt_nearby(x::SVector{2,T}, s::Polygon{K}) where {K,T<:Real}
    # Calculate the signed distances from x to edge lines.
	# ∆xe = ∆xe_poly(x,s)
	∆xe = _∆xe_poly(x,s.v,s.n)
	# Determine if x is outside of edges, inclusive.
	obd = Zygote.@ignore(onbnd(∆xe,s))
    # For x inside the polygon, it is easy to find the closest surface point: we can simply
    # pick the closest edge and find its point closest to x.
    # For x outside the polygon, there are many cases depending on how many edges x lies
    # outside of.  However, x that is sufficiently close to the polygon should be outside of
    # at most two edges.  (Outside two edges near the vertices, and one edge elsewhere.)
    # Therefore, we will consider only cout ≤ 2 below.
    co = Zygote.@ignore(cout(∆xe,obd))
    if co == 2  # x is outside two edges
        # We could choose to find ind corresponding to the two nonnegative ∆xe, but such an
        # operation leads to allocations because it does not create an SVector.  Instead,
        # find the closest vertex directly, which is the surface point we are looking for
        # for cout = 2.
        imin, imin₋₁ = imin2(x,s)
		surf = SVector{2}(s.v[imin,1],s.v[imin,2])
		@inbounds if obd[imin] && obd[imin₋₁]  # x is very close to vertex imin
            # nout = reinterpret(SVector{2,T}, normalize( [ s.n[imin,1]+s.n[imin₋₁,1],s.v[imin,2]+s.n[imin₋₁,2] ]))  #  s.n[imin,:] + s.n[imin₋₁,:]
			@inbounds nout = normalize( [ s.n[imin,1]+s.n[imin₋₁,1],s.n[imin,2]+s.n[imin₋₁,2] ])
		else
            # nout = reinterpret(SVector{2,T}, normalize( [ x[1] - s.v[imin,1] , x[2] - s.v[imin,2] ] ) )[1]
			@inbounds nout = normalize( [ x[1] - s.v[imin,1] , x[2] - s.v[imin,2] ] )
        end
    else  # cout ≤ 1 or cout ≥ 3
        # co ≤ 1 || @warn "x = $x is outside $cout edges: too far from polygon with vertices $(s.v); " *
        #                     "result could be inaccurate."
        # Choose the closest edge to x.
        # If cout = 0, all ∆xe are negative, so the largest ∆xe is the smallest in magnitude
        # and corresponds to the closest edge.
        # If cout = 1, all but one ∆xe are negative, so again the largest ∆xe is the only
        # nonnegative one and corresponds to the edge outside which x lies.
        # Even for cout = 3, this seems to be a reasonable choice from a simple geometric
        # argument.
        imax::Int = Zygote.@ignore(argmax(∆xe))
        surf = x + _∆x_poly(x,s.v,s.n,imax) # ∆x = (nmax⋅(vmax-x)) .* nmax
        nout = @inbounds SVector(s.n[imax,1],s.n[imax,2] )
    end

    return surf, nout
end

translate(s::Polygon{K,K2,D}, ∆::SVector{2,T}) where {K,K2,D,T<:Real} = Polygon{K,K2,D,T}(s.v .+ transpose(∆), s.n, s.data)

#= Factory methods =#
# Regular polygon
function regpoly(::Val{K},  # number of vertices
                 r::Real,  # distance between center and each vertex
                 θ::Real=π/2,  # angle from +x-direction towards first vertex; π/2 corresponds to +y-direction
                 c::SVector{2,T}=SVector(0.0,0.0),  # center location
                 data=nothing) where {K,T<:Real}
    ∆θ = 2π / K

    θs = θ .+ ∆θ .* SVector(ntuple(k->k-1, Val(K)))  # SVector{K}: angles of vertices
    v = c' .+ r .* [cos.(θs) sin.(θs)]  # SMatrix{K,2}: locations of vertices

    return Polygon(v, data)
end

regpoly(k::Integer,  # number of vertices
        r::Real,  # radius: distance from center to vertices
        θ::Real=π/2,  # angle of first vertex
        c::AbstractVector{<:Real}=[0.0,0.0],  # [x, y]: center of regular polygon
        data=nothing) =
   regpoly(Val(k), r, θ, SVector{2}(c), data)

# Isosceles triangle
function isosceles(base::NTuple{2,SVector{2,<:Real}},
                   h::Real,
                   data=nothing)
    m = (base[1] + base[2]) / 2  # midpoint of base
    bvec = normalize(base[2] - base[1])  # unit direction of base
    hvec = [-bvec[2], bvec[1]]  # unit direction of height
    p = m + h.*hvec  # apex

    v = [base[1] base[2] p]'  # vertices

    return Polygon(v, data)
end

isosceles(base::NTuple{2,AbstractVector{<:Real}},  # (end point 1, end point 2): two end points of base
          h::Real,  # height drawn normal to base; direction is such that base pt 1, base pt 2, apex are put in counter-clockwise order
          data=nothing) = isosceles(SVector{2}.(base), h, data)

# To-dos: parallegram, rhombus, ...

#= Polygonal prism =#
const PolygonalPrism{K,K2,T} = Prism{Polygon{K,K2,Nothing,T}}

# Below, if we called PolygonalPrism(c, ...) in the function body, it would call the inner
# constructor Prism{Polygon{K,K2,Nothing}}(c, ...) because PolygonalPrism = Prism{Polygon{K,K2,Nothing}},
# which is not what we want.
# To call the outer constructor of Prism, we should call Prism(c, ...) instead of PolygonalPrism(c, ...).
PolygonalPrism(c::SVector{3,<:Real},
               v::SMatrix{K,2,<:Real},  # 2D coordinates of base vertices in projected prism coordinates
               h::Real=Inf,
               a::SVector{3,<:Real}=SVector(0.0,0.0,1.0),
               data=nothing) where {K} =
    (â = normalize(a); Prism(c, Polygon(v), h, [orthoaxes(â)... â], data))

PolygonalPrism(c::AbstractVector{<:Real},  # center of prism
               v::AbstractMatrix{<:Real},  # vertices of base polygon
               h::Real=Inf,  # height of prism
               a::AbstractVector{<:Real}=[0.0,0.0,1.0],  # axis direction of prism
               data=nothing) =
    (K = size(v,1); PolygonalPrism(SVector{3}(c), SMatrix{K,2}(v), h, SVector{3}(a), data))

# Return the bounds of the center cut with respect to the prism center.
function bounds_ctrcut(s::PolygonalPrism{K}) where {K}
    ax = inv(s.p)  # prism axes: columns are not only unit vectors, but also orthogonal
    v = [s.b.v @SVector(zeros(K))]  # SMatrix{K,3}: 3D vectices in prism axis coordinates
    w = v * ax'  # SMatrix{K,3}: vertices in external coordinates (equivalent to (ax * v')')

    return minimum(w, dims=Val(1))[1,:], maximum(w, dims=Val(1))[1,:]
end


text:
using ApproxFun, Test
    import ApproxFun: Block

@testset "RaggedMatrix" begin
    cols=Int[rand(1:k+2) for k=1:5]
    B=ApproxFun.rrand(Float64,maximum(cols),cols)
    cols=Int[rand(1:k+2) for k=1:size(B,1)]
    A=ApproxFun.rrand(Float64,maximum(cols),cols)
    @test Matrix(A)*Matrix(B) ≈ Matrix(A*B)

    @test ApproxFun.RaggedMatrix(B) === B
    @test ApproxFun.RaggedMatrix{Float64}(B) === B
    @test Matrix(ApproxFun.RaggedMatrix{ComplexF64}(B)) == Matrix{ComplexF64}(Matrix(B))

    B = ApproxFun.brand(10,10,2,3)
    @test Matrix(B) == Matrix(ApproxFun.RaggedMatrix(B))
    @test ApproxFun.RaggedMatrix(B) == ApproxFun.RaggedMatrix{Float64}(B)
    @test ApproxFun.RaggedMatrix(ApproxFun.BandedMatrix{ComplexF64}(B)) == ApproxFun.RaggedMatrix{ComplexF64}(B)
end


text:
# Escrever sobre o arquivo
# io = open("crud.txt", "w")
# write(io, "\nEstou manipulando arquivo e aprendi o conceito de CRUD")

# Adicionar ao final do arquivo
io = open("crud.txt", "a")
write(io, "\nEstou manipulando arquivo e aprendi o conceito de CRUD")

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming