text:
"""
This  package,  which  depends  on  no other package, implements univariate
Laurent  polynomials  (type  `Pol{T}`),  and  univariate rational fractions
(type `Frac{Pol{T}}`).

The  initial motivation was to have a simple way to port GAP polynomials to
Julia. The reasons for still having my own package are multiple:

  - I need  to have  a simple  and flexible  interface, which  I hope this
    provides.
  - I need my polynomials to behave  well when coefficients are in a ring,
    in which case I use pseudo-division and subresultant gcd.
  - For my uses my polynomials are several times faster than those in the
    package `Polynomials`.
  - I need my polynomials to  work as well as possible with coefficients of
    types  `T` where elements have a `zero`  method but `T` itself does not
    have  one (because `T` does not  contain the necessary information). An
    example is modular arithmetic with a `BigInt` modulus which thus cannot
    be part of the type.
  - a final justification is that `LaurentPolynomials` is used by `PuiseuxPolynomials`.

Laurent polynomials have the parametric type `Pol{T}`, where `T`is the type
of   the  coefficients.  They  are  constructed   by  giving  a  vector  of
coefficients  of  type  `T`,  and  a  valuation  (an  `Int`).  We call true
polynomials those whose valuation is `≥0`.

There  is a  current variable  name (a  `Symbol`) used to print polynomials
nicely  at  the  repl  or  in  IJulia  or  Pluto.  This name can be changed
globally,  or just changed for printing a given polynomial. But polynomials
do not record individually with which symbol they should be printed.

# Examples
```julia-repl
julia> Pol(:q) # define symbol used for printing and return Pol([1],1)
Pol{Int64}: q

julia> @Pol q  # same as q=Pol(:q)  useful to start session with polynomials
Pol{Int64}: q

julia> Pol([1,2]) # valuation is taken to be 0 if omitted
Pol{Int64}: 2q+1

julia> 2q+1       # same polynomial
Pol{Int64}: 2q+1

julia> Pol()   # omitting all arguments gives Pol([1],1)
Pol{Int64}: q

julia> p=Pol([1,2,1],-1) # here the valuation is specified to be -1
Pol{Int64}: q+2+q⁻¹

julia> q+2+q^-1 # same polynomial
Pol{Int64}: q+2+q⁻¹
```

```julia-rep1
julia> print(p) # if not nice printing give an output which can be read back
Pol([1, 2, 1],-1)

# change the variable for printing just this time
julia> print(IOContext(stdout,:limit=>true,:varname=>"x"),p)
x+2+x⁻¹

julia> print(IOContext(stdout,:TeX=>true),p) # TeXable output (used in Pluto, IJulia)
q+2+q^{-1}
```

A  polynomial can be  taken apart with  the functions `valuation`, `degree`
and `getindex`. An index `p[i]` gives the coefficient of degree `i` of `p`.

```julia-repl
julia> valuation(p),degree(p)
(-1, 1)

julia> p[0], p[1], p[-1], p[10]
(2, 1, 1, 0)

julia> p[valuation(p):degree(p)]
3-element Vector{Int64}:
 1
 2
 1

julia> p[begin:end]  # the same as the above line
3-element Vector{Int64}:
 1
 2
 1

julia> coefficients(p)  # the same again
3-element Vector{Int64}:
 1
 2
 1
```

A  polynomial  is  a  *scalar*  if  the  valuation  and degree are `0`. The
function  `scalar` returns the constant coefficient  if the polynomial is a
scalar, and `nothing` otherwise.

```julia-repl
julia> Pol(1)
Pol{Int64}: 1

julia> convert(Pol{Int},1) # the same thing
Pol{Int64}: 1

julia> scalar(Pol(1))
1

julia> convert(Int,Pol(1)) # the same thing
1

julia> Int(Pol(1))         # the same thing
1

julia> scalar(q+1) # nothing; convert would give an error
```

In arrays `Pol{T}` of different types `T` are promoted to the same type `T`
(when  the `T`  involved have  a promotion)  and a  number is promoted to a
polynomial.

Usual  arithmetic (`+`, `-`,  `*`, `^`, `/`,  `//`, `one`, `isone`, `zero`,
`iszero`,  `==`) works. Elements  of type `<:Number`  or of type  `T` for a
`Pol{T}`   are  considered  as   scalars  for  scalar   operations  on  the
coefficients.

```julia-repl
julia> derivative(p)
Pol{Int64}: 1-q⁻²

julia> p=(q+1)^2
Pol{Int64}: q²+2q+1

julia> p/2
Pol{Float64}: 0.5q²+1.0q+0.5

julia> p//2
Pol{Rational{Int64}}: (1//2)q²+(1//1)q+1//2

julia> p(1//2) # value of p at 1//2
9//4

julia> p(0.5)
2.25

# interpolation: find p taking values [2.0,1.0,3.0] at [1,2,3]
julia> Pol([1,2,3],[2.0,1.0,3.0])
Pol{Float64}: 1.5q²-5.5q+6.0
```

Polynomials  are scalars  for broadcasting.  They can  be sorted (they have
`cmp`   and  `isless`  functions  which   compare  the  valuation  and  the
coefficients), they can be keys in a `Dict` (they have a `hash` function).

The  functions  `divrem`,  `div`,  `%`,  `gcd`,  `gcdx`,  `lcm`, `powermod`
operate  between  true  polynomials  over  a  field,  using  the polynomial
division.  Over a ring it is better  to use `pseudodiv` and `srgcd` instead
of  `divrem`  and  `gcd`  (by  default  `gcd`  between  integer polynomials
delegates to `srgcd`).

`exactdiv`  does  division  (over  a  field  or  a  ring) when it is exact,
otherwise gives an error.

```julia-repl
julia> divrem(q^3+1,2q+1) # changes coefficients to field elements
(0.5q²-0.25q+0.125, 0.875)

julia> divrem(q^3+1,2q+1//1) # case of coefficients already field elements
((1//2)q²+(-1//4)q+1//8, 7//8)

julia> pseudodiv(q^3+1,2q+1) # pseudo-division keeps the ring
(4q²-2q+1, 7)

julia> (4q^2-2q+1)*(2q+1)+7 # but multiplying back gives a multiple of the polynomial
Pol{Int64}: 8q³+8

julia> exactdiv(q+1,2.0) # exactdiv(q+1,2) would give an error
Pol{Float64}: 0.5q+0.5
```

Finally,   `Pol`s  have   methods  `conj`,   `adjoint`  which   operate  on
coefficients,  methods `positive_part`,  `negative_part` and  `bar` (useful
for  Kazhdan-Lusztig  theory)  and  a  method  `randpol`  to produce random
polynomials.

Inverting polynomials is a way to get a rational fraction  `Frac{Pol{T}}`.
Rational fractions are normalized so that the numerator and denominator are
true  polynomials prime to each other.  They have the arithmetic operations
`+`,  `-` ,  `*`, `/`,  `//`, `^`,  `inv`, `one`, `isone`, `zero`, `iszero`
(which can operate between a `Pol` or a `Number` and a `Frac{Pol{T}}`).

```julia-repl
julia> a=1/(q+1)
Frac{Pol{Int64}}: 1/(q+1)

julia> Pol(2/a) # convert back to `Pol`
Pol{Int64}: 2q+2

julia> numerator(a)
Pol{Int64}: 1

julia> denominator(a)
Pol{Int64}: q+1

julia> m=[q+1 q+2;q-2 q-3]
2×2 Matrix{Pol{Int64}}:
 q+1  q+2
 q-2  q-3

julia> n=inv(Frac.(m)) # convert to rational fractions to invert the matrix
2×2 Matrix{Frac{Pol{Int64}}}:
 (-q+3)/(2q-1)  (-q-2)/(-2q+1)
 (q-2)/(2q-1)   (q+1)/(-2q+1)

julia> map(x->x(1),n) # evaluate at 1 the inverse matrix
2×2 Matrix{Float64}:
  2.0   3.0
 -1.0  -2.0

julia> map(x->x(1;Rational=true),n) # evaluate at 1 using //
2×2 Matrix{Rational{Int64}}:
  2//1   3//1
 -1//1  -2//1
```

Rational fractions are also scalars for broadcasting and can be sorted
(have `cmp` and `isless` methods).
"""
module LaurentPolynomials
export degree, valuation, Pol, derivative, shift, positive_part, negative_part,
       bar, derivative, srgcd, Frac, @Pol, scalar, coefficients, exactdiv, root,
       randpol, pseudodiv

const varname=Ref(:x)

struct Pol{T}
  c::Vector{T}
  v::Int
  # Unexported inner constructor that bypasses all checks
  global Pol_(c::AbstractVector{T},v::Integer) where T=new{T}(c,v)
end

"""
  `Pol(c::AbstractVector,v::Integer=0;check=true,copy=true)`

  Make a polynomial of valuation `v` with coefficients `c`.

  Unless  `check` is `false`  normalize the result  by making sure that `c`
  has  no leading or  trailing zeroes (do  not set `check=false` unless you
  are sure this is already the case).

  Unless  `copy=false` the  contents of  `c` are  copied (you  can gain one
  allocation by setting `copy=false` if you know the contents can be shared)
"""
function Pol(c::AbstractVector{T},v::Integer=0;check=true,copy=true)where T
  if check # normalize c so there are no leading or trailing zeroes
    b=findfirst(!iszero,c)
    b===nothing && return Pol_(T[],0)
    e=findlast(!iszero,c)
    if b!=1 || e!=length(c) || copy return Pol_(view(c,b:e),v+b-1) end
  end
  Pol_(copy ? Base.copy(c) : c,v)
end

Pol()=Pol_([1],1)
Pol(a::T) where T<:Number=iszero(a) ? Pol_(T[],0) : Pol_([a],0)

"""
 `Pol(t::Symbol)`

 Sets the name of the variable for printing `Pol`s to `t`, and returns
 the polynomial of degree 1 equal to that variable.
"""
function Pol(t::Symbol)
  varname[]=t
  Pol()
end

"""
 `@Pol q`

 is equivalent to `q=Pol(:q)` excepted it creates `q` in the global scope of
 the current module, since it uses `eval`.
"""
macro Pol(t)
  if !(t isa Symbol) error("usage: @Pol <variable name>") end
  Base.eval(Main,:($t=Pol($(Core.QuoteNode(t)))))
end

Base.broadcastable(p::Pol)=Ref(p)
Base.copy(p::Pol)=Pol_(copy(p.c),p.v)

degree(a::Number)=0 # convenient
degree(p::Pol)=length(p.c)-1+p.v
valuation(p::Pol)=p.v
coefficients(p::Pol)=p.c
Base.lastindex(p::Pol)=degree(p)
Base.firstindex(p::Pol)=valuation(p)
@inbounds Base.getindex(p::Pol{T},i::Integer) where T=
i in firstindex(p):lastindex(p) ? p.c[i-p.v+1] : iszero(p) ? zero(T) : zero(p.c[1])

Base.getindex(p::Pol,i::AbstractVector{<:Integer})=getindex.(Ref(p),i)

Base.convert(::Type{Pol{T}},a::Number) where T=iszero(a) ? Pol_(T[],0) :
                                                           Pol_([T(a)],0)

(::Type{Pol{T}})(a) where T=convert(Pol{T},a)

# like convert for vectors does not always make a copy
Base.convert(::Type{Pol{T}},p::Pol{T1}) where {T,T1}=Pol_(Vector{T}(p.c),p.v)

function Base.promote_rule(a::Type{Pol{T1}},b::Type{Pol{T2}})where {T1,T2}
  Pol{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{Pol{T1}},b::Type{T2})where {T1,T2<:Number}
  Pol{promote_type(T1,T2)}
end

Base.isinteger(p::Pol)=isinteger(scalar(p))

function scalar(p::Pol{T})where T
  if iszero(p) T(0)
  elseif iszero(valuation(p)) && iszero(degree(p)) p[0]
  end
end

function Base.convert(::Type{T},p::Pol) where T<:Number
  res=scalar(p)
  if res===nothing throw(InexactError(:convert,T,p)) end
  convert(T,res)
end

(::Type{T})(p::Pol) where T<:Number=convert(T,p)

Base.cmp(a::Pol,b::Pol)=cmp((a.c,a.v),(b.c,b.v))
Base.isless(a::Pol,b::Pol)=cmp(a,b)==-1
Base.hash(a::Pol, h::UInt)=hash(a.v,hash(a.c,h))

(p::Pol{T})(x) where T=iszero(p) ? zero(T) : evalpoly(x,p.c)*x^p.v

# efficient p↦ qˢ p
shift(p::Pol{T},s) where T=Pol_(p.c,p.v+s)

positive_part(p::Pol)=p.v>=0 ? copy(p) : Pol(view(p.c,1-p.v:length(p.c)),0)

negative_part(p::Pol)=degree(p)<=0 ? copy(p) : Pol(view(p.c,1:1-p.v),p.v)

# q↦ q⁻¹ on p
bar(p::Pol)=Pol_(reverse(p.c),-degree(p))

Base.:(==)(a::Pol, b::Pol)= a.c==b.c && a.v==b.v
Base.:(==)(a::Pol,b)= (iszero(a) && iszero(b))||(b!==nothing && scalar(a)==b)
Base.:(==)(b,a::Pol)= a==b

Base.one(a::Pol{T}) where T=Pol_([iszero(a) ? one(T) : one(a.c[1])],0)
Base.one(::Type{Pol{T}}) where T=Pol_([one(T)],0)
Base.one(::Type{Pol})=one(Pol{Int})
Base.isone(a::Pol)=!iszero(a) && scalar(a)==1
Base.zero(::Type{Pol{T}}) where T=Pol_(T[],0)
Base.zero(::Type{Pol})=zero(Pol{Int})
Base.zero(a::Pol{T}) where T=Pol_(T[],0)
Base.iszero(a::Pol)=isempty(a.c)
# next 4 stuff to make transpose and inv using LU work (abs is stupid)
Base.conj(p::Pol{T}) where T=Pol_(conj.(p.c),p.v)
Base.abs(p::Pol)=p
Base.adjoint(a::Pol)=conj(a)
Base.transpose(a::Pol)=a

ismonomial(p::Pol)=length(p.c)==1

function Base.show(io::IO, ::MIME"text/html", a::Pol)
  print(io, "\$")
  show(IOContext(io,:TeX=>true),a)
  print(io, "\$")
end

function Base.show(io::IO, ::MIME"text/plain", a::Pol)
  if !haskey(io,:typeinfo)
    print(io,typeof(a),": ")
    io=IOContext(io,:typeinfo=>typeof(a))
  end
  show(io,a)
end

# 3 next methods copied from Util.jl in order to have a self-contained file.

# determines which coefficients should be bracketed for unambiguous display
function bracket_if_needed(c::String)
  if match(r"^[-+]?([^-+*/]|√-|{-)*(\(.*\))?$",c)!==nothing c
  else "("*c*")"
  end
end

function format_coefficient(c::String)
  if c=="1" ""
  elseif c=="-1" "-"
  else bracket_if_needed(c)
  end
end

const supvec=['⁰','¹','²','³','⁴','⁵','⁶','⁷','⁸','⁹']

function stringexp(io::IO,n::Integer)
  if isone(n) ""
  elseif get(io,:TeX,false)
    n in 0:9 ? "^"*string(n) : "^{"*string(n)*"}"
  elseif get(io,:limit,false)
    res=Char[]
    if n<0 push!(res,'⁻'); n=-n end
    for i in reverse(digits(n)) push!(res,supvec[i+1]) end
    String(res)
  else "^"*string(n)
  end
end

function Base.show(io::IO,p::Pol)
  if !get(io,:limit,false) && !get(io,:TeX,false)
    if ismonomial(p) && isone(p.c[1]) && p.v==1 print(io,"Pol()")
    else print(io,"Pol(",p.c)
      if !iszero(p.v) print(io,",",p.v) end
      print(io,")")
    end
  elseif iszero(p) print(io,"0")
  else
    var=string(get(io,:varname,varname[]))
    for deg in degree(p):-1:valuation(p)
      c=p[deg]
      if iszero(c) continue end
      c=repr(c; context=IOContext(io,:typeinfo=>typeof(c)))
      if !iszero(deg)
        c=format_coefficient(c)*var*stringexp(io,deg)
      end
      if c[1]!='-' && deg!=degree(p) c="+"*c end
      print(io,c)
    end
  end
end

function Base.:*(a::Pol{T1},b::Pol{T2})where {T1,T2}
  T=promote_type(T1,T2)
  if iszero(a) || iszero(b) return zero(Pol{T}) end
  # below not zero(T) for types T like Mod{T1} which have no zero method
  res=fill(zero(a.c[1]*b.c[1]),length(a.c)+length(b.c)-1)
  for i in eachindex(a.c), j in eachindex(b.c)
@inbounds res[i+j-1]+=a.c[i]*b.c[j]
  end
  Pol_(res,a.v+b.v)
end

Base.:*(a::Pol, b::Number)=Pol(a.c.*b,a.v)
Base.:*(a::Pol{T}, b::T) where T=Pol(a.c.*b,a.v;copy=false)
Base.:*(b::Number, a::Pol)=a*b
Base.:*(b::T, a::Pol{T}) where T=a*b

Base.:^(a::Pol, n::Real)=isinteger(n) ? a^Int(n) : root(a,1//n)

Base.:^(a::Pol, n::Integer)=ismonomial(a) ? Pol_([a.c[1]^n],n*a.v) :
         n>=0 ? Base.power_by_squaring(a,n) :
                Base.power_by_squaring(inv(a),-n)

function Base.:+(a::Pol,b::Pol)
  if iszero(a) return b elseif iszero(b) return a end
  z=zero(a.c[1]+b.c[1])
  d=b.v-a.v
  if d>=0
    res=fill(z,max(length(a.c),d+length(b.c)))
  @inbounds view(res,eachindex(a.c)).=a.c
  @inbounds view(res,d.+eachindex(b.c)).+=b.c
    Pol(res,a.v;copy=false)
  else
    res=fill(z,max(length(a.c)-d,length(b.c)))
  @inbounds view(res,eachindex(a.c).-d).=a.c
  @inbounds view(res,eachindex(b.c)).+=b.c
    Pol(res,b.v;copy=false)
  end
end

function Base.:-(a::Pol,b::Pol)
  if iszero(a) return -b elseif iszero(b) return a end
  z=zero(a.c[1]+b.c[1])
  d=b.v-a.v
  if d>=0
    res=fill(z,max(length(a.c),d+length(b.c)))
  @inbounds view(res,eachindex(a.c)).=a.c
  @inbounds view(res,d.+eachindex(b.c)).-=b.c
    Pol(res,a.v;copy=false)
  else
    res=fill(z,max(length(a.c)-d,length(b.c)))
  @inbounds view(res,eachindex(a.c).-d).=a.c
  @inbounds view(res,eachindex(b.c)).-=b.c
    Pol(res,b.v;copy=false)
  end
end

Base.:+(a::Pol, b::Number)=a+Pol(b)
Base.:+(b::Number, a::Pol)=Pol(b)+a
Base.:-(a::Pol)=Pol_(-a.c,a.v)
Base.:-(b::Number, a::Pol)=Pol(b)-a
Base.:-(a::Pol,b::Number)=a-Pol(b)
Base.div(a::Pol,b::Number)=Pol(div.(a.c,b),a.v;copy=false)

exactdiv(a,b)=a/b  # generic version for fields
function exactdiv(a::Integer,b::Integer) # define for integral domains
  (d,r)=divrem(a,b)
  if !iszero(r) error(b," does not exactly divide ",a) end
  d
end

function coeffexactdiv(a::Pol,b)
  if isone(b) return a end
  c=exactdiv.(a.c,b)
  Pol_(c,a.v)
end
exactdiv(a::Pol,b::Number)=coeffexactdiv(a,b)

Base.:/(p::Pol,q::Number)=Pol_(p.c./q,p.v)
Base.://(p::Pol,q::Number)=Pol_(p.c.//q,p.v)

derivative(a::Pol)=Pol([(i+a.v-1)*v for (i,v) in enumerate(a.c)],a.v-1,copy=false)

"""
`divrem(a::Pol, b::Pol)`

`a` and `b` should be true polynomials (nonnegative valuation).
Computes  `(q,r)` such  that `a=q*b+r`  and `degree(r)<degree(b)`.
Type stable if the coefficients of `b` are in a field.
"""
function Base.divrem(a::Pol, b::Pol)
  if iszero(b) throw(DivideError) end
  if degree(b)>degree(a) return (zero(a),a) end
  if a.v<0 || b.v<0 error("arguments should be true polynomials") end
  d=inv(b.c[end])
  z=zero(a.c[1]+b.c[1]+d)
  r=fill(z,1+degree(a))
  view(r,a.v+1:length(r)).=a.c
  q=fill(z,length(r)-degree(b))
  for i in length(r):-1:degree(b)+1
    if iszero(r[i]) c=zero(d)
    else c=r[i]*d
         view(r,i-length(b.c)+1:i) .-= c .* b.c
    end
    q[i-degree(b)]=c
  end
  Pol(q),Pol(r)
end

function exactdiv(a::Pol,b::Pol)
  if isone(b) || iszero(a) return a end
  if iszero(b) throw(DivideError) end
  d=a.v-b.v
  if !iszero(a.v) a=shift(a,-a.v) end
  if !iszero(b.v) b=shift(b,-b.v) end
  if degree(b)>degree(a) error(b," does not exactly divide ",a) end
  z=zero(a.c[1]+b.c[1])
  r=fill(z,1+degree(a))
  view(r,a.v+1:length(r)).=a.c
  q=fill(z,length(r)-degree(b))
  for i in length(r):-1:degree(b)+1
    c=exactdiv(r[i],b.c[end])
    view(r,i-length(b.c)+1:i) .-= c .* b.c
    q[i-length(b.c)+1]=c
  end
  if !iszero(r) error(b," does not exactly divide ",a) end
  res=Pol(q,d)
end

"""
`pseudodiv(a::Pol, b::Pol)`

pseudo-division  of `a` by `b`.  If `d` is the  leading coefficient of `b`,
computes   `(q,r)`   such   that   `d^(degree(a)+1-degree(b))a=q*b+r`   and
`degree(r)<degree(b)`. Does not do division so works over any ring.
For true polynomials (errors if the valuation of `a` or of `b` is negative).

See Knuth AOCP2 4.6.1 Algorithm R
"""
function pseudodiv(a::Pol, b::Pol)
 if isone(b) || iszero(a) return a,zero(a) end
  if iszero(b) throw(DivideError) end
  d=b.c[end]
  if degree(a)<degree(b) return (zero(a),d^(degree(a)+1-degree(b))*a) end
  if a.v<0 || b.v<0 error("arguments should be true polynomials") end
  z=zero(a.c[1]+b.c[1])
  r=fill(z,1+degree(a))
  view(r,a.v+1:length(r)).=a.c
  q=fill(z,length(r)-degree(b))
  for i in length(r):-1:degree(b)+1
    c=r[i]
    r.*=d
    q.*=d
    if !iszero(c)
      for j in eachindex(b.c) r[j+i-length(b.c)]-=c*b.c[j] end
    end
    q[i-degree(b)]=c
  end
  Pol(q),Pol(r)
end

"""
`srgcd(a::Pol,b::Pol)`

sub-resultant gcd: gcd of polynomials over a unique factorization domain

See Knuth AOCP2 4.6.1 Algorithm C
"""
function srgcd(a::Pol,b::Pol)
  if degree(b)>degree(a) a,b=b,a end
  if iszero(b) return a end
  ca=gcd(a.c);a=coeffexactdiv(a,ca)
  cb=gcd(b.c);b=coeffexactdiv(b,cb)
  d=gcd(ca,cb)
  g=1
  h=1
  while true
    δ=degree(a)-degree(b)
    q,r=pseudodiv(a,b)
    if iszero(r)
      cb=gcd(b.c);b=coeffexactdiv(b,cb)
      return isone(d) ? b : Pol_(b.c .*d,b.v)
    elseif degree(r)==0
      return Pol_([d],0)
    end
    a=b
    gh=g*h^δ
    b=coeffexactdiv(r,gh)
    g=a[end]
    if δ>0 h=exactdiv(g^δ,h^(δ-1)) end
  end
end

Base.gcd(p::Pol{<:Integer},q::Pol{<:Integer})=srgcd(p,q)
Base.gcd(v::AbstractArray{<:Pol})=reduce(gcd,v)
Base.lcm(p::Pol,q::Pol)=exactdiv(p*q,gcd(p,q))
Base.lcm(m::AbstractArray{<:Pol})=reduce(lcm,m)

Base.div(a::Pol, b::Pol)=divrem(a,b)[1]
Base.:%(a::Pol, b::Pol)=divrem(a,b)[2]
Base.:%(a::Pol{<:Integer}, b::Pol{<:Integer})=divrem(a,b)[2]

"""
`gcd(p::Pol,  q::Pol)` computes the  `gcd` of the  polynomials. It uses the
subresultant algorithms for the `gcd` of integer polynomials.

```julia-repl
julia> gcd(2q+2,q^2-1)
Pol{Int64}: q+1

julia> gcd(q+1//1,q^2-1//1)
Pol{Rational{Int64}}: (1//1)q+1//1
```
"""
function Base.gcd(p::Pol,q::Pol)
  if degree(q)>degree(p)
    p,q=q,p
  end
  p,q=promote(p,q)
  while !iszero(q)
    q=q/q.c[end]
    (q,p)=(divrem(p,q)[2],q)
  end
  p*inv(p.c[end])
end

"""
  `gcdx(a::Pol,b::Pol)`

for  polynomials  over  a  field  returns `d,u,v`  such  that `d=ua+vb` and
`d=gcd(a,b)`.

```julia-repl
julia> gcdx(q^3-1//1,q^2-1//1)
((1//1)q-1//1, 1//1, (-1//1)q)
```
"""
function Base.gcdx(a::Pol, b::Pol)
  a,b=promote(a, b)
  # a0, b0=a, b
  s0, s1=one(a), zero(a)
  t0, t1=s1, s0
  # The loop invariant is: s0*a0 + t0*b0 == a
  x,y=a,b
  while y != 0
    q,q1=divrem(x, y)
    x, y=y, q1
    s0, s1=s1, s0 - q*s1
    t0, t1=t1, t0 - q*t1
  end
  (x, s0, t0)./x[end]
end

"""
`powermod(p::Pol, x::Integer, q::Pol)` computes ``p^x \\pmod m``.
```julia-repl
julia> powermod(q-1,3,q^2+q+1)
Pol{Int64}: 6q+3
```
"""
function Base.powermod(p::Pol, x::Integer, q::Pol)
  x==0 && return one(q)
  b=p%q
  t=prevpow(2, x)
  r=one(q)
  while true
    if x>=t
     r=(r*b)%q
      x-=t
    end
    t >>>= 1
    t<=0 && break
    r=(r*r)%q
  end
  r
end

"""
`randpol(T,d)`

random polynomial of degree `d` with coefficients from `T`
"""
randpol(T,d::Integer)=Pol(rand(T,d+1))

"""
`Pol(x::AbstractVector,y::AbstractVector)`

Interpolation:  find a `Pol` (of  nonnegative valuation) of smallest degree
taking  values `y` at points  `x`. The values `y`  should be in a field for
the function to be type stable.

```julia-repl
julia> p=Pol([1,1,1])
Pol{Int64}: q²+q+1

julia> vals=p.(1:5)
5-element Vector{Int64}:
  3
  7
 13
 21
 31

julia> Pol(1:5,vals*1//1)
Pol{Rational{Int64}}: (1//1)q²+(1//1)q+1//1

julia> Pol(1:5,vals*1.0)
Pol{Float64}: 1.0q²+1.0q+1.0
```
"""
function Pol(pts::AbstractVector,vals::AbstractVector)
  vals=copy(vals)
  a=map(eachindex(pts))do i
    for k in i-1:-1:1
      if pts[i]==pts[k] error("interpolating points must be distinct") end
      vals[k]=(vals[k+1]-vals[k])/(pts[i]-pts[k])
    end
    vals[1]
  end
  p=Pol_([a[end]],0)
  for i in length(pts)-1:-1:1
    p=p*(Pol()-pts[i])+a[i]
  end
  p
end

Base.denominator(p::Pol)=iszero(p) ? 1 : lcm(denominator.(p.c))
Base.numerator(p::Pol{<:Rational{T}}) where T=convert(Pol{T},p*denominator(p))

function root(x::Pol,n::Union{Integer,Rational{<:Integer}}=2)
  n=Int(n)
  if iszero(x) return x end
  if !ismonomial(x) || !iszero(x.v%n) error("root($x,$n) not implemented") end
  c=x.c[1]
  Pol_([isone(c) ? c : root(c,n)],div(x.v,n))
end

root(x::AbstractFloat,n=2)=x^(1/n)
#---------------------- Frac-------------------------------------
## cannot use Rational{T} since it forces T<:Integer
struct Frac{T}
  num::T
  den::T
  # Unexported inner constructor that bypasses all checks
  global Frac_(num::T,den::T) where T=new{T}(num,den)
end

Base.numerator(a::Frac)=a.num
Base.denominator(a::Frac)=a.den

function Base.convert(::Type{Frac{T}},p::Frac{T1}) where {T,T1}
  Frac_(convert(T,p.num),convert(T,p.den))
end

(::Type{Frac{T}})(a::Frac{T}) where T=a
(::Type{Frac{T}})(a::Number) where T=convert(Frac{T},a)

Base.broadcastable(p::Frac)=Ref(p)

Base.copy(a::Frac)=Frac_(a.num,a.den)
Base.one(a::Frac)=Frac_(one(a.num),one(a.den))
Base.one(::Type{Frac{T}}) where T =Frac_(one(T),one(T))
Base.isone(a::Frac)=a.num==a.den
Base.zero(::Type{Frac{T}}) where T =Frac_(zero(T),one(T))
Base.zero(a::Frac)=Frac_(zero(a.num),one(a.num))
Base.iszero(a::Frac)=iszero(a.num)
# next 3 methods are to make inv using LU work (abs is stupid)
Base.abs(p::Frac)=p
Base.conj(p::Frac)=Frac_(conj(p.num),conj(p.den))
Base.adjoint(a::Frac)=conj(a)
Base.:(==)(a::Frac,b::Frac)=a.num==b.num && a.den==b.den
Base.cmp(a::Frac,b::Frac)=cmp((a.num,a.den),(b.num,b.den))
Base.isless(a::Frac,b::Frac)=cmp(a,b)==-1

function Base.show(io::IO, ::MIME"text/plain", a::Frac)
  if !haskey(io,:typeinfo)
    print(io,typeof(a),": ")
    io=IOContext(io,:typeinfo=>typeof(a))
  end
  show(io,a)
end

function Base.show(io::IO,a::Frac)
  if !get(io, :limit, false) && !get(io, :TeX, false)
    print(io,"Frac(",a.num,",",a.den,")")
    return
  end
  if haskey(io,:typeinfo) && isone(a.den)
    print(io,a.num)
    return
  end
  n=sprint(show,a.num; context=io)
  print(io,bracket_if_needed(n))
  n=sprint(show,a.den; context=io)
  print(io,"/",bracket_if_needed(n))
end

Base.inv(a::Frac)=Frac_(a.den,a.num)

Base.:*(a::Frac,b::Frac)=Frac(a.num*b.num,a.den*b.den)
Base.://(a::Frac,b::Frac)=Frac(a.num*b.den,a.den*b.num)
Base.:/(a::Frac,b::Frac)=a//b

Base.:^(a::Frac, n::Integer)= n>=0 ? Base.power_by_squaring(a,n) :
                              Base.power_by_squaring(inv(a),-n)
Base.:+(a::Frac,b::Frac)=Frac(a.num*b.den+a.den*b.num,a.den*b.den)
Base.:-(a::Frac)=Frac_(-a.num,a.den)
Base.:-(a::Frac,b::Frac)=Frac(a.num*b.den-a.den*b.num,a.den*b.den)
Base.:+(a::Frac{<:T},b::Union{Number,T}) where T=+(promote(a,b)...)
Base.:+(b::Union{Number,T},a::Frac{<:T}) where T=+(promote(a,b)...)
Base.:-(a::Frac{<:T},b::Union{Number,T}) where T=-(promote(a,b)...)
Base.:-(b::Union{Number,T},a::Frac{<:T}) where T=-(promote(b,a)...)
Base.:*(a::Frac{<:T},b::T) where T=Frac(a.num*b,a.den)
Base.:*(b::T,a::Frac{<:T}) where T=Frac(a.num*b,a.den)

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

# make both pols, one of valuation 0
function make_positive(a::Pol,b::Pol)
  v=a.v-b.v
  shift(a,max(v,0)-a.v),shift(b,-min(v,0)-b.v)
end

Frac(a::Pol,b::Pol;prime=false)=Frac(promote(a,b)...;prime)

"""
`Frac(a::Pol,b::Pol;prime=false)

Polynomials  `a` and `b` are promoted to same coefficient type, and checked
for  being true polynomials (otherwise they are both multiplied by the same
power  of  the  variable  so  they  become  true  polynomials),  and unless
`prime=true` they are checked for having a non-trivial `gcd`.
"""
function Frac(a::T,b::T;prime=false)::Frac{T} where T<:Pol
  if iszero(a) return Frac_(a,one(a))
  elseif iszero(b) error("division by 0")
  end
  a,b=make_positive(a,b)
  if !(prime || ismonomial(b) || ismonomial(a))
    d=gcd(a,b)
    a,b=exactdiv(a,d),exactdiv(b,d)
  end
  if scalar(b)==-1 a,b=(-a,-b) end
  Frac_(a,b)
end

function Frac(a::Pol{<:Rational},b::Pol{<:Rational};k...)
  Frac(numerator(a)*denominator(b),numerator(b)*denominator(a);k...)
end

function Pol(p::Frac{<:Pol})
  if ismonomial(p.den)==1
    if isone(p.den.c[1]^2) return Pol(p.num.c .*p.den.c[1],p.num.v-p.den.v)
    else return Pol(p.num.c .//p.den.c[1],p.num.v-p.den.v)
    end
  end
  error("cannot convert ",p," to Pol")
end

Base.convert(::Type{Pol{T}},p::Frac) where {T}=convert(Pol{T},Pol(p))

function Base.convert(::Type{Frac{T}},p::Pol) where T
  Frac(convert(T,p),one(T);prime=true)
end

function Base.convert(::Type{Frac{T}},p::Number) where {T<:Pol}
  Frac_(convert(T,p),convert(T,1))
end

function Base.promote_rule(a::Type{T1},b::Type{Frac{T2}})where {T1<:Pol,T2<:Pol}
  Frac{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{T1},b::Type{Frac{T2}})where {T1<:Number,T2<:Pol}
  Frac{promote_type(T1,T2)}
end

function Frac(a::Pol)
  if a.v>0 return Frac_(a,one(a)) end
  Frac(a,one(a);prime=true)
end

function Base.://(a::Pol,b::Pol)
  if ismonomial(b) return Pol(a.c.//b.c[1],a.v-b.v) end
  Frac(a,b)
end

bestinv(x)=isone(x) ? x : isone(-x) ? x : inv(x)

function Base.inv(p::Pol)
  if ismonomial(p) return Pol_([bestinv(p.c[1])],-p.v) end
  Frac_(make_positive(one(p),p)...)
end

Base.:/(p::Pol,q::Pol)=p*inv(q)

Base.://(a::Frac{<:Pol},b::Pol)=Frac(a.num,a.den*b)
Base.://(a::Frac{<:Pol},b::Number)=Frac(a.num,a.den*b;prime=true)
Base.:/(a::Frac{<:Pol},b::Union{Number,Pol})=a//b
Base.://(a::Union{Number,Pol},b::Frac{<:Pol})=a*inv(b)
Base.:/(a::Union{Number,Pol},b::Frac{<:Pol})=a//b
Base.://(p::Number,q::Pol)=Frac(Pol(p),q;prime=true)
Base.:/(p::Number,q::Pol)=p*inv(q)

Base.:*(a::Frac{<:Pol},b::Number)=Frac(a.num*b,a.den;prime=true)
Base.:*(b::Number,a::Frac{<:Pol})=a*b

(p::Frac{<:Pol})(x;Rational=false)=Rational ? p.num(x)//p.den(x) : p.num(x)/p.den(x)
# @btime inv(Frac.([q+1 q+2;q-2 q-3])) setup=(q=Pol())
# 32.556 μs (938 allocations: 57.22 KiB)
end


text:
function bubbleSort(lst)
    len = length(lst)
    for i in 1:len-1
        j = len
        while j > i
            if lst[j] < lst[j-1]
                lst[j], lst[j-1] = lst[j-1], lst[j]
            end
            j -= 1
        end
    end
    lst
end

lst = [8, 11, 2, 1, 4, 19, 9, 14, 13, 5, 6, 22]
lst = bubbleSort(lst)
println(lst)

ext: jl

max_stars_repo_name: Bingwen-Hu/hackaway

max_issues_repo_name: Bingwen-Hu/hackaway

max_forks_repo_name: Bingwen-Hu/hackaway


text:
module GHWT_2d

include("utils.jl")

using ..GraphSignal, ..GraphPartition, ..BasisSpecification, ..GHWT, LinearAlgebra, SparseArrays

include("common.jl")

export ghwt_analysis_2d, ghwt_bestbasis_2d, ghwt_synthesis_2d

"""
    dmatrix = ghwt_analysis_2d(matrix::Array{Float64,2}, GProws::GraphPart, GPcols::GraphPart)
For a matrix, equipped with row and column weight recursive partitionings
and weight matrices, generate the redundant matrix of GHWT
expansion coefficients.

### Input Arguments
   `matrix`              the matrix to be analyzed
   `GProws`              the recursive partitioning on the rows
   `GPcols`              the recursive partitioning on the columns

### Output Arguments
   `dmatrix`                the GHWT dictionary expansion coefficients

Copyright 2019 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_analysis_2d(matrix::Array{Float64,2}, GProws::GraphPart, GPcols::GraphPart)
    if isempty(GProws.tag)
        ghwt_core!(GProws)
    end

    if isempty(GPcols.tag)
        ghwt_core!(GPcols)
    end

    matrix_re = matrix[GProws.ind, GPcols.ind]

    (N, jmax_row) = size(GProws.rs)
    N = N-1

    # expand on the row direction
    (frows, fcols) = size(matrix_re)
    dmatrix = zeros((N,jmax_row,fcols))
    dmatrix[:,jmax_row,:] = matrix_re

    ghwt_core!(GProws, dmatrix)

    dmatrix = reshape(dmatrix, (frows*jmax_row, fcols))'

    # expand on the column direction
    (N, jmax_col) = size(GPcols.rs)
    N = N - 1
    dmatrix2 = zeros(N,jmax_col,size(dmatrix,2))
    dmatrix2[:, jmax_col, :] = dmatrix
    ghwt_core!(GPcols, dmatrix2)

    dmatrix = reshape(dmatrix2, (fcols*jmax_col, frows*jmax_row))'

    return Array{Float64,2}(dmatrix)
end

"""
    dvec, BSrows, BScols = ghwt_2d_bestbasis(matrix::Array{Float64,2},GProws::GraphPart,GPcols::GraphPart,
    costfun_rows::Any = 1.0, costfun_cols::Any = 1.0, flatten_rows::Any = 1.0, flatten_cols::Any = 1.0)

For a matrix, equipped with row and column weight recursive partitionings
and weight matrices, generate the (non-redundant) matrix of GHWT
expansion coefficients, using the best basis algorithm to select best bases for the rows and columns

### Input Argument
*   `matrix`              the matrix to be analyzed
*   `GProws`              the recursive partitioning on the rows
*   `GPcols`              the recursive partitioning on the columns
*   `costfun_rows`        the cost functional to be used for the rows
*   `costfun_cols`        the cost functional to be used for the columns
*   `flatten_rows`        the method for flattening vector-valued data to scalar-valued data for the rows
*   `flatten_cols`        the method for flattening vector-valued data to scalar-valued data for the columns

### Output Argument
*   `dvec`                the GHWT expansion coefficients (not redundant)
*   `BSrows`              the best basis on the rows
*   `BScols`              the best basis on the columns

Copyright 2019 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_bestbasis_2d(matrix::Array{Float64,2},GProws::GraphPart,GPcols::GraphPart;
    costfun_rows::Any = 1.0, costfun_cols::Any = 1.0, flatten_rows::Any = 1.0, flatten_cols::Any = 1.0)

    # generate GraphSig objects for the matrix
    rows,cols = size(matrix)
    Grows = GraphSig(spzeros(rows,rows),f = matrix)
    Gcols = GraphSig(spzeros(cols,cols),f = Matrix{Float64}(matrix'))

    # analyze the data matrix using the rows
    dmatrix = ghwt_analysis!(Grows, GP = GProws)

    # analyze the data matrix using the columns
    dmatrix2 = ghwt_analysis!(Gcols, GP = GPcols)

    # find the row and column best bases
    dvec,BSrows = ghwt_bestbasis(dmatrix,GProws,cfspec = costfun_rows,flatten = flatten_rows)
    ~,BScols = ghwt_bestbasis(dmatrix2,GPcols,cfspec = costfun_cols,flatten = flatten_cols)

    # analyze the row best-basis coefficient matrix using the columns
    Gcols.f = Matrix{Float64}(dvec')
    dmatrix = ghwt_analysis!(Gcols, GP = GPcols)

    # extract the col best-basis coefficients
    dvec = Matrix{Float64}(dmatrix2dvec(dmatrix, GPcols, BScols)')

    return dvec, BSrows, BScols
end

"""
    dvec = ghwt_2d_dmatrix2dvec(matrix::Array{Float64,2},BSrows::BasisSpec,BScols::BasisSpec)
For a matrix, equipped with row and column weight recursive partitionings
and weight matrices, generate the (non-redundant) matrix of GHWT
expansion coefficients, using the best basis algorithm to select best bases for the rows and columns

### Input Argument
*   `matrix`              the matrix to be analyzed
*   `GProws`              the recursive partitioning on the rows
*   `GPcols`              the recursive partitioning on the columns

### Output Argument
*   `dvec`                the GHWT expansion coefficients (not redundant)

Copyright 2019 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_2d_haar(matrix::Array{Float64,2},GProws::GraphPart,GPcols::GraphPart)

    # generate GraphSig objects for the matrix
    rows,cols = size(matrix)
    Grows = GraphSig(spzeros(rows,rows),f = matrix)
    Gcols = GraphSig(spzeros(cols,cols),f = Matrix{Float64}(matrix'))

    # analyze the data matrix using the rows
    dmatrix = ghwt_analysis!(Grows, GP = GProws)

    # analyze the data matrix using the columns
    #dmatrix2 = ghwt_analysis!(Gcols, GP = GPcols)

    # find the row and column best bases
    #BH=bs_haar(GP)
    #dvec=dmatrix2dvec(dc2f, GP, BH)
    BSrows = bs_haar(GProws)
    dvec = dmatrix2dvec(dmatrix, GProws, BSrows)
    BScols = bs_haar(GPcols)
    #dvec,BSrows = ghwt_bestbasis(dmatrix,GProws,cfspec = costfun_rows,flatten = flatten_rows)
    #~,BScols = ghwt_bestbasis(dmatrix2,GPcols,cfspec = costfun_cols,flatten = flatten_cols)

    # analyze the row best-basis coefficient matrix using the columns
    Gcols.f = Matrix{Float64}(dvec')
    dmatrix = ghwt_analysis!(Gcols, GP = GPcols)

    # extract the col best-basis coefficients
    dvec = Matrix{Float64}(dmatrix2dvec(dmatrix, GPcols, BScols)')

    return dvec, BSrows, BScols
end

"""
matrix = ghwt_synthesis_2d(dvec::Matrix{Float64},GProws::GraphPart,GPcols::GraphPart,BSrows::BasisSpec,BScols::BasisSpec)

Given a vector of GHWT expansion coefficients and info about the graph
partitioning and the choice of basis, reconstruct the signal

### Input Arguments
*   `dvec`        the expansion coefficients corresponding to the chosen basis
*   `GProws`      the recursive partitioning on the rows
*   `GPcols`      the recursive partitioning on the columns
*   `BSrows`      the best basis on the rows
*   `BScols`      the best basis on the columns

### Output Arguments
*   matrix      the reconstructed matrix

Copyright 2019 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""

function ghwt_synthesis_2d(dvec::Matrix{Float64},GProws::GraphPart,GPcols::GraphPart,BSrows::BasisSpec,BScols::BasisSpec)
    # reconstruct
    dvec = ghwt_synthesis(Matrix{Float64}(dvec'), GPcols, BScols)
    matrix = ghwt_synthesis(Matrix{Float64}(dvec'), GProws, BSrows)
    return matrix
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CustomResourceSubresourceScale defines how to serve the scale subresource for CustomResources.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale(;
        labelSelectorPath=nothing,
        specReplicasPath=nothing,
        statusReplicasPath=nothing,
    )

    - labelSelectorPath::String : labelSelectorPath defines the JSON path inside of a custom resource that corresponds to Scale `status.selector`. Only JSON paths without the array notation are allowed. Must be a JSON Path under `.status` or `.spec`. Must be set to work with HorizontalPodAutoscaler. The field pointed by this JSON path must be a string field (not a complex selector struct) which contains a serialized label selector in string form. More info: https://kubernetes.io/docs/tasks/access-kubernetes-api/custom-resources/custom-resource-definitions#scale-subresource If there is no value under the given path in the custom resource, the `status.selector` value in the `/scale` subresource will default to the empty string.
    - specReplicasPath::String : specReplicasPath defines the JSON path inside of a custom resource that corresponds to Scale `spec.replicas`. Only JSON paths without the array notation are allowed. Must be a JSON Path under `.spec`. If there is no value under the given path in the custom resource, the `/scale` subresource will return an error on GET.
    - statusReplicasPath::String : statusReplicasPath defines the JSON path inside of a custom resource that corresponds to Scale `status.replicas`. Only JSON paths without the array notation are allowed. Must be a JSON Path under `.status`. If there is no value under the given path in the custom resource, the `status.replicas` value in the `/scale` subresource will default to 0.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale <: SwaggerModel
    labelSelectorPath::Any # spec type: Union{ Nothing, String } # spec name: labelSelectorPath
    specReplicasPath::Any # spec type: Union{ Nothing, String } # spec name: specReplicasPath
    statusReplicasPath::Any # spec type: Union{ Nothing, String } # spec name: statusReplicasPath

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale(;labelSelectorPath=nothing, specReplicasPath=nothing, statusReplicasPath=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale, Symbol("labelSelectorPath"), labelSelectorPath)
        setfield!(o, Symbol("labelSelectorPath"), labelSelectorPath)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale, Symbol("specReplicasPath"), specReplicasPath)
        setfield!(o, Symbol("specReplicasPath"), specReplicasPath)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale, Symbol("statusReplicasPath"), statusReplicasPath)
        setfield!(o, Symbol("statusReplicasPath"), statusReplicasPath)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale = Dict{Symbol,Symbol}(Symbol("labelSelectorPath")=>Symbol("labelSelectorPath"), Symbol("specReplicasPath")=>Symbol("specReplicasPath"), Symbol("statusReplicasPath")=>Symbol("statusReplicasPath"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale = Dict{Symbol,String}(Symbol("labelSelectorPath")=>"String", Symbol("specReplicasPath")=>"String", Symbol("statusReplicasPath")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale)
    (getproperty(o, Symbol("specReplicasPath")) === nothing) && (return false)
    (getproperty(o, Symbol("statusReplicasPath")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceScale }, name::Symbol, val)
end


text:
import ONNXNaiveNASflux: fluxlayers, sources, actfuns, invariantops, pseudotransparentops, optype, nodes
using ONNXNaiveNASflux.NaiveNASflux

# Logging to avoid CI timeouts
@info "  Test padding and sources"

@testset "Read padding" begin
    import ONNXNaiveNASflux: prev

    @test prev(2) == 2
    @test prev([1,2]) == [1,2]
    @test prev([1,2,3,4]) == [2,4,1,3]
    @test prev([1,2,3,4,5,6]) == [3,6,2,5,1,4]
end

@testset "Sources" for tc in
    (
    (name="test_constant", ninputs=0, noutputs=1),
    )

    model, gb, inputs, outputs = prepare_node_test(tc.name, tc.ninputs, tc.noutputs)

    @testset "$(tc.name) op $(optype(node))" for node in nodes(gb)
        @test haskey(sources, optype(node))
        res = sources[optype(node)](node.attribute, params(node)...)

        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end

    @testset "$(tc.name) graph" begin
        cg = load(model)
        res = cg()
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]

        # Also test that it we get the same thing by serializing and then deserializing
        io = PipeBuffer()
        save(io, cg)
        cg = load(io)
        res = cg()
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end
end

# For testing since ONNX states that recurrent layers take 3D input while flux uses
# an Array of 2D Arrays
function (l::Flux.Recur)(x::AbstractArray{T, 3}) where T
    # ONNX shape for RNNs inputs is [seq_length, batch_size, input_size]
    # ONNX.jl reverses this to [input_size, batch_size, seq_length]
    # Unstacking it to a sequence of [input_size, batch_size]
    inseq =Flux.unstack(x, 3)
    out = nothing
    for inpt in inseq
         out = l(inpt)
     end
    # Just to turn it back to ONNX shape.
    # In the testdata only the last output in the sequence is present in the reference
    return reshape(out, size(out)..., 1)
end

@info "  Test Flux layers"

@testset "Fluxlayer $(tc.name)" for tc in
    (
    (name="test_averagepool_1d_default", ninputs=1, noutputs=1),
    # (name="test_averagepool_2d_ceil", ninputs=1, noutputs=1), Not supported!
    (name="test_averagepool_2d_default", ninputs=1, noutputs=1),
    #(name="test_averagepool_2d_pads", ninputs=1, noutputs=1), Not supported!
    (name="test_averagepool_2d_strides", ninputs=1, noutputs=1),
    (name="test_averagepool_3d_default", ninputs=1, noutputs=1),
    (name="test_basic_conv_with_padding", ninputs=2, noutputs=1),
    (name="test_basic_conv_without_padding", ninputs=2, noutputs=1),
    (name="test_batchnorm_epsilon", ninputs=5, noutputs=1),
    (name="test_batchnorm_example", ninputs=5, noutputs=1),
    (name="test_conv_with_strides_and_asymmetric_padding", ninputs=2, noutputs=1),
    (name="test_conv_with_strides_no_padding", ninputs=2, noutputs=1),
    (name="test_conv_with_strides_padding", ninputs=2, noutputs=1),
    (name="test_dropout_default", ninputs=1, noutputs=1),
    (name="test_dropout_random", ninputs=1, noutputs=1),
    #(name="test_gemm_all_attributes", ninputs=3, noutputs=1), Not supported!
    (name="test_gemm_alpha", ninputs=3, noutputs=1),
    (name="test_gemm_beta", ninputs=3, noutputs=1),
    #(name="test_gemm_default_matrix_bias", ninputs=3, noutputs=1), Not supported!
    (name="test_gemm_default_no_bias", ninputs=2, noutputs=1),
    (name="test_gemm_default_scalar_bias", ninputs=3, noutputs=1),
    (name="test_gemm_default_single_elem_vector_bias", ninputs=3, noutputs=1),
    (name="test_gemm_default_vector_bias", ninputs=3, noutputs=1),
    (name="test_gemm_default_zero_bias", ninputs=3, noutputs=1),
    #(name="test_gemm_transposeA", ninputs=3, noutputs=1), Not supported!
    (name="test_gemm_transposeB", ninputs=3, noutputs=1),
    (name="test_lstm_defaults", ninputs=3, noutputs=1),
    (name="test_lstm_with_initial_bias", ninputs=4, noutputs=1),
    # (name="test_lstm_with_peepholes", ninputs=8, noutputs=1), Not supported!
    (name="test_maxpool_1d_default", ninputs=1, noutputs=1),
    #(name="test_maxpool_2d_ceil", ninputs=1, noutputs=1), Not supported!
    (name="test_maxpool_2d_default", ninputs=1, noutputs=1),
    #(name="test_maxpool_2d_dilations", ninputs=1, noutputs=1), Not supported!
    #(name="test_maxpool_2d_pads", ninputs=1, noutputs=1), Not supported!
    (name="test_maxpool_2d_strides", ninputs=1, noutputs=1),
    (name="test_maxpool_3d_default", ninputs=1, noutputs=1),
    (name="test_maxpool_3d_default", ninputs=1, noutputs=1),
    (name="test_rnn_seq_length", ninputs=4, noutputs=1),
    )

    model, gb, inputs, outputs = prepare_node_test(tc.name, tc.ninputs, tc.noutputs)

    @testset "$(tc.name) op $(optype(node))" for node in nodes(gb)
        @test haskey(fluxlayers, optype(node))
        op = fluxlayers[optype(node)](node.attribute, params(node)...)

        res = op(Float32.(inputs[1]))
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end

    @testset "$(tc.name) graph" begin
        cg = load(model)
        res = cg(Float32.(inputs[1]))
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]

        # Also test that it we get the same thing by serializing and then deserializing
        io = PipeBuffer()
        save(io, cg)
        cg = load(io)
        res = cg(Float32.(inputs[1]))
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end
end

@info "  Test Flux activation functions"

@testset "Activation functions $(tc.name)" for tc in
    (
    (name="test_elu", ninputs=1, noutputs=1),
    (name="test_elu_default", ninputs=1, noutputs=1),
    (name="test_elu_example", ninputs=1, noutputs=1),
    (name="test_relu", ninputs=1, noutputs=1),
    (name="test_selu", ninputs=1, noutputs=1),
    (name="test_selu_default", ninputs=1, noutputs=1),
    (name="test_selu_example", ninputs=1, noutputs=1),
    )

    model, gb, inputs, outputs = prepare_node_test(tc.name, tc.ninputs, tc.noutputs)

    @testset "$(tc.name) op $(optype(node))" for node in nodes(gb)
        @test haskey(actfuns, optype(node))
        op = actfuns[optype(node)](node.attribute, params(node)...)
        @test op.(inputs[1]) ≈ outputs[1]

        @test haskey(invariantops, optype(node))
        bcop = invariantops[optype(node)](node.attribute, params(node)...)
        @test bcop(inputs[1]) ≈ outputs[1]

    end
end

@info "  Test stateless ops"

@testset "Invariant op $(tc.name)" for tc in
    (
    (name="test_flatten_axis0", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_axis1", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_axis2", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_axis3", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_default_axis", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_negative_axis1", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_negative_axis2", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_flatten_negative_axis3", ninputs=1, noutputs=1, fd=pseudotransparentops),
    (name="test_globalaveragepool", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_globalaveragepool_precomputed", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_globalmaxpool", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_globalmaxpool_precomputed", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_default_axes_keepdims_example", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_default_axes_keepdims_random", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_do_not_keepdims_example", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_do_not_keepdims_random", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_keepdims_example", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_keepdims_random", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_negative_axes_keepdims_example", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reduce_mean_negative_axes_keepdims_random", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_reshape_extended_dims", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_negative_dim", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_negative_extended_dims", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_one_dim", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_reduced_dims", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_reordered_all_dims", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_reordered_last_dims", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_zero_and_negative_dim", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_reshape_zero_dim", ninputs=2, noutputs=1, fd=pseudotransparentops),
    (name="test_softmax_axis_0", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_axis_1", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_axis_2", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_default_axis", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_example", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_large_number", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_softmax_negative_axis", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_squeeze", ninputs=1, noutputs=1, fd=invariantops),
    (name="test_squeeze_negative_axes", ninputs=1, noutputs=1, fd=invariantops)
    )

    model, gb, inputs, outputs = prepare_node_test(tc.name, tc.ninputs, tc.noutputs)

    @testset "$(tc.name) op $(optype(node))" for node in nodes(gb)
        @test haskey(tc.fd, optype(node))
        op = tc.fd[optype(node)](node.attribute, params(node)...)
        res = op(inputs[1])
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end

    @testset "$(tc.name) graph" begin
        cg = load(model, size(inputs[1]))
        res = cg(inputs[1])
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]

        # Also test that it we get the same thing by serializing and then deserializing
        io = PipeBuffer()
        save(io, cg)
        cg = load(io, size(inputs[1]))
        res = cg(inputs[1])
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end
end

@testset "Vertex $(tc.name)" for tc in
    (
    (name="test_add", ninputs=2, noutputs=1),
    #(name="test_add_bcast", ninputs=2, noutputs=1), # Op is supported, but we get the wrong idea about what type of inputvertex to create from 3D input
    (name="test_concat_1d_axis_0", ninputs=2, noutputs=1),
    (name="test_concat_1d_axis_negative_1", ninputs=2, noutputs=1),
    (name="test_concat_2d_axis_0", ninputs=2, noutputs=1),
    (name="test_concat_2d_axis_1", ninputs=2, noutputs=1),
    (name="test_concat_2d_axis_negative_1", ninputs=2, noutputs=1),
    (name="test_concat_2d_axis_negative_2", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_0", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_1", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_2", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_negative_1", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_negative_2", ninputs=2, noutputs=1),
    (name="test_concat_3d_axis_negative_3", ninputs=2, noutputs=1),
    (name="test_mul", ninputs=2, noutputs=1),
    )

    model, gb, inputs, outputs = prepare_node_test(tc.name, tc.ninputs, tc.noutputs)

    @testset "$(tc.name) graph" begin
        cg = load(model)
        res = cg(inputs[1:length(cg.inputs)]...)
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]

        # Also test that it we get the same thing by serializing and then deserializing
        io = PipeBuffer()
        save(io, cg)
        cg = load(io)
        res = cg(inputs[1:length(cg.inputs)]...)
        @test size(res) == size(outputs[1])
        @test res ≈ outputs[1]
    end
end

@testset "Deserialize with inputs" begin
    using NaiveNASflux: GenericFlux2D, layertype

    function sumgraph()
        ivs = denseinputvertex.(["in1", "in2"], 4)
        g_org = CompGraph(ivs, "out" >> ivs[1] + ivs[2])
        pb = PipeBuffer()
        save(pb, g_org, "in1" => missing, "in2" => missing)
        return pb
    end

    insize(t::Tuple) = ONNXNaiveNASflux.int_size(t[NaiveNASflux.actdim(length(t))])
    insize(p::Pair) = p |> last |> insize
    @testset "Input format $inshapes" for inshapes in (
        ((4,1), (4,1)),
        ("in1" => (4,1), "in2" => (4,1)),
        ((4,missing), (4, :B)),
        ((:I, 3), (:I, 4))
    )
        expsizes =  insize.(inshapes) |> collect
        g_new = if any(==(0), expsizes)
            @test_logs (:warn, r"No valid input sizes") load(sumgraph(), inshapes...)
        else
            load(sumgraph(), inshapes...)
        end
        @test nout.(g_new.inputs) == expsizes
        @test layertype.(g_new.inputs) == [GenericFlux2D(), GenericFlux2D()]
    end

    inshape(t::Tuple) = t |> length |> ONNXNaiveNASflux.guess_layertype
    @testset "Mixshape format $inshapes" for inshapes in (
        ((1,1,5,1), (5,1)),
        ((5,1), (1,1,5,1)),
    )
        g_new = load(sumgraph(), inshapes...)
        @test nout.(g_new.inputs) == [5, 5]
        @test layertype.(g_new.inputs) == inshape.(inshapes |> collect)
    end

    @testset "Malformed input $inshapes" for inshapes in (
        ((4,1), (4,1), (4,1)),
        ("in1" => (4,1), "in2" => (4,1), "in2" => (4,1)),
        ("in1" => (4,1), "notin2" => (4,1))
    )
        @test_throws AssertionError load(sumgraph(), inshapes...)
    end
end

@testset "Deserialize with merging" begin
    function remodel(f, args...)
        pb = PipeBuffer()
        save(pb, f, args...)
        return load(pb, args...)
    end

    @testset "Merge activation function" begin
        m = remodel(Dense(3,4, relu), (3, missing))
        @test nvertices(m) == 2
        @test layer(m.outputs[1]).σ == relu
    end

    @testset "Merge Reshape and $gp" for gp in (
        ONNXNaiveNASflux.globalmeanpool,
        ONNXNaiveNASflux.globalmaxpool
    )
        m = remodel(Chain(
            Conv((3,3), 3 => 3),
            x -> gp(x, xf -> reshape(xf, 3, :)),
        ), (4, 4, 3, missing))

        @test nvertices(m) == 3
    end

    @testset "Merge constant" begin
        m = remodel(x -> relu.(x) .+ (3))
        @test nvertices(m) == 3
        @test m([-1,1]) == [3, 4]
    end
end

@testset "Infer CompGraph shapes" begin
    function remodel(f, args...;kwargs...)
        pb = PipeBuffer()
        save(pb, f, args...)
        return load(pb,; kwargs...)
    end

    @testset "Simple Dense" begin
        v1 = denseinputvertex("v1", 3)
        v2 = fluxvertex("v2", Dense(nout(v1), 4), v1)
        g = remodel(CompGraph(v1, v2), missing)

        @test nout(inputs(g)[]) == nout(v1)
    end

    @testset "After invariant" begin
        v1 = denseinputvertex("v1", 3)
        v2 = invariantvertex("v2", identity, v1)
        v3 = fluxvertex("v3", Dense(nout(v2), 4), v2)
        g = remodel(CompGraph(v1, v3), missing)

        @test nout(inputs(g)[]) == nout(v1)
    end

    @testset "Can't infer after concat" begin
        # I suppose in this case we could infer it as we know nin(v3) = 2 * nout(v1)
        # Seems like too much of an edge case to be worth considering though
        v1 = denseinputvertex("v1", 3)
        v2 = concat("v2", v1, v1)
        v3 = fluxvertex("v3", Dense(nout(v2), 4), v2)
        g = @test_logs (:warn, r"No valid input sizes") remodel(CompGraph(v1, v3), (missing, :B))

        @test nout(inputs(g)[]) == 0
    end

    @testset "Flatten$label" for (label, layerfun, exputilsize) in
        (
        ("", identity, 1),
        ("with ActivationContribution", ActivationContribution, 60),
        )
        using ONNXNaiveNASflux: Flatten, create_vertex_default, defaultutility

        v1 = conv2dinputvertex("v1", 3)
        v2 = fluxvertex("v2", Conv((2,2), 3=>5), v1)
        v3 = absorbvertex("v3", Flatten(-1), v2)

        g = remodel(CompGraph(v1, v3), (5,4,3,:B); vfun=(args...) -> create_vertex_default(args...; layerfun))

        @test nout(inputs(g)[]) == nout(v1)
        @test nout(g[end]) == 60
        @test length(defaultutility(g[end])) == exputilsize
    end
end


text:
# This file is part of Kpax3. License is MIT.

#=
R = [1; 1; 2; 2; 3; 1]
k = length(unique(R))

priorR = Kpax3.EwensPitman(settings.α, settings.θ)
priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

ss = UInt8[1; 2; 3]
obj = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)

M = -Inf
for s1 in ss, s2 in ss, s3 in ss, s4 in ss
  S = [s1; s2; s3; s4]
  lp = Kpax3.logcondpostS(S, obj.cl, obj.k, obj.v, obj.n1s, priorC)

  if lp > M
    M = lp
  end
end

trueSp = zeros(Float64, 3, 4)
for s1 in ss, s2 in ss, s3 in ss, s4 in ss
  S = [s1; s2; s3; s4]
  lp = Kpax3.logcondpostS(S, obj.cl, obj.k, obj.v, obj.n1s, priorC)

  trueSp[s1, 1] += exp(lp - M)
  trueSp[s2, 2] += exp(lp - M)
  trueSp[s3, 3] += exp(lp - M)
  trueSp[s4, 4] += exp(lp - M)
end

trueSp = exp(M + log(trueSp))

for i in 1:12
  @printf("%.100f\n", trueSp[i])
end

Kpax3.rpostpartitioncols!(obj.C, obj.cl, obj.k, obj.v, obj.n1s, priorC);
obj.C[obj.cl[1:obj.k], :]

obj.C[obj.cl[1:obj.k], :] = Ctest1

p = exp(Kpax3.logcondpostC(obj.C, obj.cl, obj.k, obj.v, obj.n1s, priorC))
@printf("%.100f\n", p)
=#

function test_partition_cols_constructor()
  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  n1s = Float64[1; 5; 3; 3; 2; 2; 3; 3; 1; 3; 2; 1; 3; 2; 4; 2; 3; 3]

  r1 = 2.0
  r2 = 100.0

  A1 = zeros(Float64, 4, m)
  A1[1, :] .= (r1 + 1.0) .* (n1s .+ 0.5) ./ (n + 1)
  A1[2, :] .= 1.0
  A1[3, :] .= 1.0
  A1[4, :] .= r1

  B1 = zeros(Float64, 4, m)
  B1[1, :] .= (r1 + 1.0) .- A1[1, :]
  B1[2, :] .= 1.0
  B1[3, :] .= r1
  B1[4, :] .= 1.0

  A2 = zeros(Float64, 4, m)
  A2[1, :] .= n1s .+ 0.5
  A2[2, :] .= 1.0
  A2[3, :] .= 1.0
  A2[4, :] .= r2

  B2 = zeros(Float64, 4, m)
  B2[1, :] .= n .- n1s .+ 0.5
  B2[2, :] .= 1.0
  B2[3, :] .= r2
  B2[4, :] .= 1.0

  for k in 1:n
    for γ in ([1.0; 0.0; 0.0], [0.0; 1.0; 0.0], [0.0; 0.0; 1.0], [0.4; 0.3; 0.3], [0.5; 0.3; 0.2], [0.7; 0.2; 0.1], [0.1; 0.1; 0.1], [0.3; 0.1; 0.1], [0.0; 0.2; 0.1])
      x1 = Kpax3.AminoAcidPriorCol(data, γ, r1)
      x2 = Kpax3.AminoAcidPriorCol(data, γ, r2)

      @test isapprox(x1.logγ[1], log(γ[1] / sum(γ)), atol=ε)
      @test isapprox(x2.logγ[2], log(γ[2] / sum(γ)), atol=ε)
      @test isapprox(x2.logγ[3], log(γ[3] / sum(γ)), atol=ε)

      @test isapprox(x1.logω[k][1], log(1.0 - 1.0 / k), atol=ε)
      @test isapprox(x2.logω[k][2], log(1.0 / k), atol=ε)

      for b in 1:m, s in 1:4
        @test isapprox(x1.A[s, b], A1[s, b], atol=ε)
        @test isapprox(x1.B[s, b], B1[s, b], atol=ε)

        @test isapprox(x2.A[s, b], A2[s, b], atol=ε)
        @test isapprox(x2.B[s, b], B2[s, b], atol=ε)
      end
    end
  end

  nothing
end

test_partition_cols_constructor()

function test_partition_cols_functions()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  data = UInt8[1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 0 0 1;
               1 0 0 0 1 1]

  (m, n) = size(data)

  settings = Kpax3.KSettings(ifile, ofile, maxclust=6, maxunit=1)

  R = [1; 1; 2; 2; 3; 1]
  k = length(unique(R))

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)

  C = zeros(UInt8, settings.maxclust, m)

  cs = (UInt8[1; 1; 1], UInt8[2; 2; 2], UInt8[3; 3; 3], UInt8[3; 3; 4], UInt8[3; 4; 3], UInt8[4; 3; 3], UInt8[3; 4; 4], UInt8[4; 3; 4], UInt8[4; 4; 3], UInt8[4; 4; 4])

  logp1 = zeros(Float64, (2 + 2^k)^m)
  logp2 = zeros(Float64, (2 + 2^k)^m)
  logp3 = zeros(Float64, (2 + 2^k)^m)

  l = 0

  for c1 in cs, c2 in cs, c3 in cs, c4 in cs
    l += 1

    tmp = hcat(c1, c2, c3, c4)

    state.C[state.cl[1], :] = C[1, :] = tmp[1, :]
    state.C[state.cl[2], :] = C[2, :] = tmp[2, :]
    state.C[state.cl[3], :] = C[3, :] = tmp[3, :]

    logp1[l] = Kpax3.logpriorC(state.C, state.cl, k, priorC)
    logp2[l] = Kpax3.logpriorC(C, k, priorC)
    logp3[l] = Kpax3.logcondpostC(state.C, state.cl, k, state.v, state.n1s, priorC)
  end

  M = maximum(logp1)
  p1 = exp(M + log(sum(exp.(logp1 .- M))))

  M = maximum(logp2)
  p2 = exp(M + log(sum(exp.(logp2 .- M))))

  M = maximum(logp3)
  p3 = exp(M + log(sum(exp.(logp3 .- M))))

  @test isapprox(p1, 1.0, atol=ε)
  @test isapprox(p2, 1.0, atol=ε)
  @test isapprox(p3, 1.0, atol=ε)

  ss = UInt8[1; 2; 3]
  logp4 = zeros(Float64, 3^m)
  l = 0

  for s1 in ss, s2 in ss, s3 in ss, s4 in ss
    l += 1
    S = [s1; s2; s3; s4]
    logp4[l] = Kpax3.logcondpostS(S, state.cl, k, state.v, state.n1s, priorC)
  end

  M = maximum(logp4)
  p4 = exp(M + log(sum(exp.(logp4 .- M))))

  @test isapprox(p4, 1.0, atol=ε)

  nothing
end

test_partition_cols_functions()

function test_partition_cols_simulations()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  data = UInt8[1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 0 0 1;
               1 0 0 0 1 1]

  (m, n) = size(data)

  settings = Kpax3.KSettings(ifile, ofile, maxclust=1, maxunit=1)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  state = Kpax3.AminoAcidState(data, [1; 1; 2; 2; 3; 1], priorR, priorC, settings)

  N = 1000000

  Sp = zeros(Float64, 3, m)
  Cp = zero(Float64)
  Ctest = UInt8[2 1 4 1; 2 1 3 1; 2 1 3 1]

  trueSp = hcat([0.514127006003589515081841909704962745308876037597656250000,
                 0.389747596850238686716494385109399445354938507080078125000,
                 0.096125397146171798201663705185637809336185455322265625000],
                [0.469929717191958307154209251166321337223052978515625000000,
                 0.356242671218090556362056986472452990710735321044921875000,
                 0.173827611589951080972582531103398650884628295898437500000],
                [0.400605930740164040138040491001447662711143493652343750000,
                 0.402853999936236717438475807284703478217124938964843750000,
                 0.196540069323599214667908086084935348480939865112304687500],
                [0.653732577496086042501133306359406560659408569335937500000,
                 0.344352962631683023886353112175129354000091552734375000000,
                 0.001914459872230926673619677558235707692801952362060546875])

  trueCp = 0.023444249569984754177909280770109035074710845947265625

  for t in 1:N
    Kpax3.rpostpartitioncols!(state.C, state.cl, state.k, state.v, state.n1s, priorC)

    for b in 1:m
      if state.C[state.cl[1], b] == 0x01
        @test state.C[state.cl[2], b] == 0x01
        @test state.C[state.cl[3], b] == 0x01

        Sp[1, b] += 1.0
      elseif state.C[state.cl[1], b] == 0x02
        @test state.C[state.cl[2], b] == 0x02
        @test state.C[state.cl[3], b] == 0x02

        Sp[2, b] += 1.0
      else
        Sp[3, b] += 1.0
      end
    end

    if all(state.C[state.cl[1:state.k], :] .== Ctest)
      Cp += 1.0
    end
  end

  Sp /= N

  @test isapprox(Sp, trueSp, rtol=0.01)

  Cp /= N

  @test isapprox(Cp, trueCp, atol=0.005)

  nothing
end

test_partition_cols_simulations()


text:
"""
    TraceLogger{M,V}
Logs a move, movedelta, acceptance and timestamp every time Δ is called. Must be used in a Cassette-context MetaheuristicLogging.
"""
struct TraceLogger{M,V}
    moves::Vector{M}
    deltas::Vector{V}
    accepted::Vector{Bool}
    times::Vector{UInt64}
end

TraceLogger(::Type{M}, ::Type{V}) where {M,V} = TraceLogger(M[], V[], Bool[], UInt64[])

function log!(l::TraceLogger, move, Δv,a, t)
    push!(l.times, t)
    push!(l.moves, move)
    push!(l.accepted, a)
    push!(l.deltas, Δv)
    nothing
end

function Base.empty!(t::TraceLogger)
    empty!(t.moves)
    empty!(t.accepted)
    empty!(t.deltas)
    empty!(t.times)
    return t
end

Tables.istable(::TraceLogger) = true
Tables.rowaccess(::Type{TraceLogger{M,V}}) where {M,V} = true
Tables.columnaccess(::Type{TraceLogger{M,V}}) where {M,V} = false

function Tables.rows(t::TraceLogger{M,V}) where{M,V}
    map( x-> (iteration = x[1], move = x[2], Δv = x[3], accepted = x[4], Δt = x[5]), Iterators.zip(eachindex(t.moves), t.moves, t.deltas, t.accepted, Float64.(t.times)))
end


text:
module TCX
using EzXML, Dates, DataFrames, Geodesy, Mocking
import Base.show

export parse_tcx_dir, parse_tcx_file, getActivityType, getDataFrame, getDistance, getDistance2, getDuration, getAverageSpeed, getAveragePace

const OK = 200
const CLIENT_ERROR = 400
const CLIENT_TCX_ERROR = 401
const NOT_FOUND = 404
const SERVER_ERROR = 500

struct TrackPoint
    Time::DateTime
    Latitude::Float64
    Longtitude::Float64
    HeartRateBpm::Int32
    AltitudeMeter::Float64
    DistanceMeter::Float64
end

struct TCXRecord
    Id::DateTime
    Name::String
    ActivityType::String
    DistanceStatic::Float64
    DurationStatic::Float64
    HeartRate::Int32
    TrackPoints::Array{TrackPoint}
end

function parse_tcx(tcxdoc::EzXML.Document)
    root_element = root(tcxdoc)
    # Check if TCX
    if nodename(root_element) != "TrainingCenterDatabase"
        return CLIENT_TCX_ERROR, nothing
    end

    activities = findnode(root_element, "Activities")
    activity = findnode(activities, "Activity")
    aType = activity["Sport"]
    id = findnode(activity, "Id")
    xid = nodecontent(id)
    aId = convertToDateTime(xid)
    lap = findnode(activity, "Lap")
    aName = nodecontent(lap)
    xtime = findnode(lap, "TotalTimeSeconds")
    aTime = parse(Float64, nodecontent(xtime))
    xDistance= findnode(lap, "DistanceMeters")
    aDistance = nodecontent0(Float64, xDistance)
    xbpm = findnode(lap, "AverageHeartRateBpm", false)
    aHeartRateBpm = (xbpm === nothing) ? 0 : nodecontent0(Int32, firstelement(xbpm))
    tp_Points = findnode(lap, "Track")
    aTrackPoints = Array{TrackPoint, countelements(tp_Points)}[]
    for tp in elements(tp_Points)
        xtime = nodecontent(findnode(tp, "Time"))
        tp_time = convertToDateTime(xtime)
        position = findnode(tp, "Position", false)
        if position == nothing
            tp_lat = tp_lon = 0.0
        else
            tp_lat = nodecontent0(Float64, findnode(position, "LatitudeDegrees"))
            tp_lon = nodecontent0(Float64, findnode(position, "LatitudeDegrees"))
        end
        xbpm = findnode(tp, "HeartRateBpm", false)
        tp_bpm = (xbpm === nothing) ? 0 : nodecontent0(Int32, firstelement(xbpm))
        tp_dist = nodecontent0(Float32, findnode(tp, "DistanceMeters", false))
        tp_alt = nodecontent0(Float64, findnode(tp, "AltitudeMeters", false))

        aTrackPoints = vcat(aTrackPoints, TrackPoint(tp_time, tp_lat, tp_lon, tp_bpm, tp_alt, tp_dist))
    end

    return OK, TCXRecord(aId, aName, aType, aDistance, aTime, aHeartRateBpm, aTrackPoints)
end

function parse_tcx_str(str::String)
    try
        status, parsed_tcx = parse_tcx(EzXML.parsexml(str))
        warn_on_tcx_error(status, str, false)
        return status, parsed_tcx
    catch e
        if isa(e, EzXML.XMLError)
            @error "Invalid XML string: $str"
            return CLIENT_ERROR, nothing
        end
    end
end

function parse_tcx_file(file::String)
    file_path = abspath(file)
    if isfile(file_path) == false
        return NOT_FOUND, nothing
    end
    xmldoc = try @mock EzXML.readxml(file_path)
    catch e
       if isa(e, EzXML.XMLError)
           # Not a valid XML document
           @warn "Invalid XML document: $file_path"
           return CLIENT_ERROR, nothing
       else
           throw(e)
       end
    end

    status, parsed_tcx = parse_tcx(xmldoc)
    warn_on_tcx_error(status, file_path, true)

    return status, parsed_tcx
end

function warn_on_tcx_error(status::Int, thing::String, isFile::Bool)
    if status == CLIENT_TCX_ERROR
        @warn "Invalid TCX $(isFile ? document : string): $(thing)"
    end
end

# Returns nodecontent(node), or zero::dType if node==nothing.
function nodecontent0(dType, node)
    if node !== nothing
        return parse(dType, nodecontent(node))
    else
        return dType(0)
    end
end

# Returns the first occurance of a child node of node
# named name. Returns nothing if there
# is no such node. Issue a warning if warn is true.
function findnode(node::EzXML.Node, name::String, warn::Bool=true)
    n = firstnode(node)
    while(nodename(n) !== name && hasnextnode(n))
        n = nextnode(n)
    end
    if nodename(n) == name
        return n
    else
        if warn
            @warn "Can't find $name"
        end
        return nothing
    end
end

function parse_tcx_dir(path::String)
    if ispath(path) == false
        @warn "Invalid path: $path"
        return SERVER_ERROR, nothing
    end

    tcxArray = Array{TCXRecord}[]
    searchdir(path, key) = filter(x->occursin(key, x), readdir(path))

    for f in searchdir(path, ".tcx")
        err, tcx = parse_tcx_file(joinpath(path, f))
        if err == OK
            tcxArray = vcat(tcxArray, tcx)
        end
    end

    if length(tcxArray) > 0
        return OK, tcxArray
    else
        return NOT_FOUND, nothing
    end
end

function getActivityType(record::TCXRecord)
    return record.ActivityType
end

function getDataFrame(record::TCXRecord)
    return DataFrame(record.TrackPoints)
end

function getDataFrame(tcxArray::Array{Any, 1})
    aTP = Array{TrackPoint}[]
    for t in tcxArray
        aTP = vcat(aTP, t.TrackPoints)
    end
    return DataFrame(aTP)
end

function getDistance(record::TCXRecord)
    return record.DistanceStatic
end

function getDistance2(record::TCXRecord)
    total_distance = 0
    df = getDataFrame(record)
    num_of_rows = size(df, 1)
    for i in 1:num_of_rows
        if i < num_of_rows
            total_distance += distance(
                                       LLA(df[i, :Latitude], df[i, :Longtitude], df[i, :AltitudeMeter]),
                                       LLA(df[i+1, :Latitude], df[i+1, :Longtitude], df[i+1, :AltitudeMeter])
               )
        end
    end
    return total_distance
end

function getAverageSpeed(record::TCXRecord)
    return (record.DistanceStatic /1000) / (record.DurationStatic / 3600)  # km/h
end

function getAveragePace(record::TCXRecord)
    return (record.DurationStatic / 60) / (record.DistanceStatic / 1000) # min/km
end

function getDuration(record::TCXRecord)
    return record.DurationStatic
end

#=
= Converts a datetime string into the proper datetime based on string length.
=
= Will assume that an ArgumentError is due to
= https://github.com/JuliaLang/julia/issues/23049 and will attempt to work
= around this.
=#
function convertToDateTime(datestr::String)::DateTime
    m = match(r"\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(Z?|\.\d{1,3}Z?)", datestr)
    format_prefix = "yyyy-mm-ddTHH:MM:SS"
    if m === nothing
        msg = "'$(datestr)' is improperly formatted. Must be in the form "
        msg = msg * "'$(format_prefix)Z' or '$(format_prefix).sssZ'"
        throw(ArgumentError(msg))
    else
        suffix = replace(m.captures[1], r"\d" => "s")
        try
            return @mock DateTime(m.match, format_prefix * suffix)
        catch e
            if isa(e, ArgumentError)
                # OK! FINE! NO Z THEN!
                return DateTime(m.match[1:end-1], format_prefix * (suffix[1:end-1]))
            else
                throw(e)
            end
        end
    end
end

Base.show(io::IO, tcx::TCXRecord) = print(io, "$(tcx.ActivityType) $(tcx.DistanceStatic/1000) km at $(tcx.Id) for $(tcx.DurationStatic) seconds.")
end #module_end


text:
# ---
# title: 1277. Count Square Submatrices with All Ones
# id: problem1277
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/count-square-submatrices-with-all-ones/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix of ones and zeros, return how many **square**
# submatrices have all ones.
#
#
#
# **Example 1:**
#
#
#
#     Input: matrix =
#     [
#       [0,1,1,1],
#       [1,1,1,1],
#       [0,1,1,1]
#     ]
#     Output: 15
#     Explanation:
#     There are **10** squares of side 1.
#     There are **4** squares of side 2.
#     There is  **1** square of side 3.
#     Total number of squares = 10 + 4 + 1 = **15**.
#
#
# **Example 2:**
#
#
#
#     Input: matrix =
#     [
#       [1,0,1],
#       [1,1,0],
#       [1,1,0]
#     ]
#     Output: 7
#     Explanation:
#     There are **6** squares of side 1.
#     There is **1** square of side 2.
#     Total number of squares = 6 + 1 = **7**.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 300`
#   * `1 <= arr[0].length <= 300`
#   * `0 <= arr[i][j] <= 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Base.Iterators: drop, take
using IterTools: iterated, imap

"""
    fixedpoint(f, guess, stopfun)

return the result of applying `guess = f(guess)`
until convergence. Convergence is decided by applying
`stopfun(guess)` which returns a Boolean.
"""
function fixedpoint(f, guess, stopfun)
    for state in iterated(f, guess)
        stopfun(state) && return state
    end
end

mutable struct StopFunction{T,S}
    olderror::T
    counter::Int
    tol::S
    maxit::Int
    minit::Int
end

"""
    (st::StopFunction)(state)

stopfunction for vumps, returning true if error is smaller than expectation or the maximum
number of iterations is reached. Implemented as a closure since it needs to remember
the last error it saw for comparison.
"""
function (st::StopFunction)(state)
    st.counter += 1
    st.counter > st.maxit - 1 && return true

    error = state[2]
    error <= st.tol && st.counter > st.minit - 1 && return true
    st.olderror = error

    return false
end


text: my_f(x,y) = 2x+3y


text:
for (T1, T2) in ((:T, :(float(T))), (:(Union{Missing,T}), :(Union{Missing,float(T)})))
  @eval begin

    # unweighted windowed function application

    function rolling(fun2::Function, data1::V, data2::V, windowspan::Int) where {T, V<:AbstractVector{$T1}}
        nvals  = nrolled(min(length(data1),length(data2)), windowspan)
        offset = windowspan - 1
        result = zeros($T2, nvals)

        @inbounds for idx in eachindex(result)
            result[idx] = fun2( view(data1, idx:idx+offset), view(data2, idx:idx+offset) )
        end

        return result
    end

    # weighted windowed function application
    #=
    function rolling(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::F) where
                     {T, N<:Number, V<:AbstractVector{T}, F<:Vector{N}}

        length(weighting) != windowspan &&
            throw(WeightsError(length(weighting), windowspan))

        nvals  = nrolled(min(length(data1), length(data2)), windowspan)
        offset = windowspan - 1
        result = zeros(float(T), nvals)

        @inbounds for idx in eachindex(result)
            result[idx] = fun2( view(data1, idx:idx+offset) .* weighting,  view(data2, idx:idx+offset) .* weighting )
        end

        return result
    end
    =#

    # weighted windowed function application
    function rolling(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::F) where
                        {T, V<:AbstractVector{$T1}, N<:Number, F<:Vector{N}}

        length(weighting) != windowspan &&
           throw(WeightsError(length(weighting), windowspan))

        nvals  = nrolled(min(length(data1), length(data2)), windowspan)
        offset = windowspan - 1
        result = zeros($T2, nvals)
        curwin1 = zeros($T2, windowspan)
        curwin2 = zeros($T2, windowspan)
        v = view(weighting,1:length(weighting))

        @inbounds for idx in eachindex(result)
           for i=1:windowspan
               j = i + idx - 1
               curwin1[i] = data1[j] * v[i]
               curwin2[i] = data2[j] * v[i]
           end
           result[idx] = fun2( curwin1, curwin2 )
        end

        return result
    end

    # weighted windowed function application
    function rolling(fun2::Function, data1::V, data2::V, windowspan::Int, weighting1::F, weighting2::F) where
                        {T, V<:AbstractVector{$T1}, N<:Number, F<:Vector{N}}

        length(weighting1) != windowspan &&
           throw(WeightsError(length(weighting1), windowspan))
        length(weighting2) != windowspan &&
           throw(WeightsError(length(weighting2), windowspan))

        nvals  = nrolled(min(length(data1), length(data2)), windowspan)
        offset = windowspan - 1
        result = zeros($T2, nvals)
        curwin1 = zeros($T2, windowspan)
        curwin2 = zeros($T2, windowspan)
        v1 = view(weighting1,1:windowspan)
        v2 = view(weighting2,1:windowspan)

        @inbounds for idx in eachindex(result)
           for i=1:offset
               j = i + idx - 1
               curwin1[i] = data1[j] * v1[i]
               curwin2[i] = data2[j] * v2[i]
           end
           result[idx] = fun2( curwin1, curwin2 )
        end

        return result
    end

    rolling(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::W) where
                    {T, V<:AbstractVector{$T1}, W<:AbstractWeights} =
        rolling(fun2, data1, data2, windowspan, weighting.values)

    # unweighted windowed function tapering

    function tapers2(fun2::Function, data1::V, data2::V) where {T, V<:AbstractVector{$T1}}
        nvals  = min(length(data1), length(data2))
        result = zeros($T2, nvals)

        @inbounds for idx in nvals:-1:1
            result[idx] = fun2( view(data1, 1:idx), view(data2, 1:idx) )
        end

        return result
    end

    function tapers2(fun2::Function, data1::V, data2::V, ntocopy::Int) where {T, V<:AbstractVector{$T1}}
        nvals  = min(length(data1), length(data2))
        ntocopy = min(nvals, max(0, ntocopy))

        result = zeros($T2, nvals)
        if ntocopy > 0
           result[1:ntocopy] = data[1:ntocopy]
        end
        ntocopy += 1

        @inbounds for idx in nvals:-1:ntocopy
            result[idx] = fun2( view(data1, 1:idx), view(data2, 1:idx) )
        end

        return result
    end

    function tapers2(fun2::Function,  data1::V, data2::V, trailing_data::V) where
            {T, V<:AbstractVector{$T1}}

        ntrailing = axes(trailing_data)[1].stop
        nvals  = min(length(data1), length(data2)) + ntrailing
        result = zeros($T2, nvals)

        result[1:ntrailing] = trailing_data[1:ntrailing]
        ntrailing += 1

        @inbounds for idx in nvals:-1:ntrailing
            result[idx] = fun2( view(data1, 1:idx),  view(data2, 1:idx) )
        end

        return result
    end

    # weighted windowed function tapering

    function tapers2(fun2::Function, data1::V, data2::V, weighting::F) where
                     {T, N<:Number, V<:AbstractVector{$T1}, F<:Vector{N}}

        nvals  = min(length(data1), length(data2))
        nweighting = length(weighting)

        nweighting == nvals || throw(WeightsError(length(weighting), nvals))

        result = zeros($T2, nvals)

        @inbounds for idx in nvals:-1:1
            wts = normalize(view(weighting, (nweighting-idx+1):nweighting))
            result[idx] = fun2( view(data1, 1:idx) .* wts, view(data2, 1:idx) .* wts   )
        end

        return result
    end

    tapers2(fun2::Function, data1::V, data2::V, weighting::W) where
                    {T, V<:AbstractVector{$T1}, W<:AbstractWeights} =
        tapers2(fun2, data1, data2, weighting.values)

  end
end


text:
# DAY 2 PROBLEM 2 ADVENT OF CODE
# Get the final position of the sub starting from (0,0)

commands = """forward 2
down 4
down 1
down 4
forward 3
down 6
down 5
forward 3
forward 8
down 2
down 3
up 8
down 5
up 7
down 7
forward 5
up 2
down 6
forward 7
forward 1
forward 2
forward 7
up 7
forward 6
down 3
down 1
up 9
down 2
up 1
down 1
up 6
forward 6
down 7
forward 6
up 1
down 6
forward 2
up 7
forward 4
forward 8
forward 7
down 7
forward 8
down 1
down 6
down 7
forward 4
down 3
up 7
down 5
down 9
up 8
up 4
down 2
down 3
up 7
forward 6
forward 6
forward 8
forward 2
up 5
down 8
down 3
down 3
down 4
down 9
down 6
up 6
forward 4
down 6
forward 3
forward 3
down 4
down 8
down 2
up 5
up 5
forward 3
forward 5
down 7
forward 6
forward 9
forward 8
forward 2
down 3
down 3
down 7
down 1
down 1
down 1
down 2
down 8
down 6
forward 6
up 1
up 6
down 7
down 1
forward 1
up 2
up 8
up 8
forward 2
down 1
down 8
down 7
down 1
forward 1
down 9
up 3
down 3
forward 2
down 3
up 6
down 2
forward 7
down 9
down 6
down 1
forward 6
down 4
down 1
down 3
forward 3
down 5
forward 9
down 5
down 7
up 8
forward 8
forward 8
down 6
down 1
forward 8
down 4
up 4
up 4
up 2
forward 2
forward 2
down 1
up 8
down 1
down 7
forward 5
down 9
down 2
up 3
down 1
down 5
forward 6
down 7
up 3
forward 7
down 4
down 3
forward 4
up 8
down 4
forward 4
forward 2
forward 5
down 5
up 2
forward 4
down 4
forward 6
down 4
forward 1
down 5
forward 2
forward 2
down 8
forward 4
forward 7
down 3
up 3
forward 2
forward 6
forward 8
down 2
forward 4
down 2
up 9
down 9
down 2
forward 5
up 4
forward 2
down 2
down 3
forward 1
down 2
forward 8
forward 8
down 4
forward 6
down 3
down 3
down 5
forward 8
forward 4
forward 1
up 4
up 2
forward 8
down 8
forward 2
forward 6
up 1
up 5
forward 2
forward 4
forward 7
forward 8
forward 2
down 3
down 1
down 9
down 6
up 5
up 6
forward 6
down 3
down 2
down 1
forward 5
forward 2
forward 7
down 8
down 7
forward 7
up 8
forward 7
down 1
up 4
forward 9
forward 4
forward 1
down 3
down 9
down 7
forward 1
down 3
forward 3
down 4
down 7
forward 4
up 6
down 8
up 1
forward 6
forward 1
down 7
down 8
up 9
up 4
down 3
down 7
forward 8
up 2
up 6
forward 8
down 1
up 4
up 4
forward 8
down 2
down 4
down 3
forward 5
down 8
forward 1
down 2
forward 9
forward 3
up 6
down 6
forward 6
forward 4
down 6
down 3
down 3
forward 6
down 5
up 4
down 9
down 3
down 6
up 9
forward 6
down 2
forward 7
up 8
down 3
down 7
down 9
forward 6
down 1
forward 2
down 1
down 3
down 3
forward 5
forward 2
up 5
forward 4
up 7
down 9
forward 7
forward 3
down 6
forward 1
down 1
up 8
down 9
up 3
down 7
up 9
forward 7
down 7
down 9
forward 9
forward 7
up 9
down 7
down 2
down 7
up 2
down 3
down 9
down 6
forward 7
forward 8
forward 8
forward 6
forward 9
forward 4
down 4
down 5
down 7
forward 6
forward 2
forward 4
forward 9
down 4
forward 6
down 7
up 1
down 7
forward 9
forward 7
down 4
down 3
up 6
forward 8
forward 7
down 8
forward 4
up 6
up 4
forward 9
forward 4
forward 4
forward 7
down 1
up 6
forward 8
forward 3
up 6
forward 4
down 1
up 2
forward 1
down 5
forward 5
up 4
down 6
down 3
up 8
forward 9
down 2
forward 4
forward 8
down 9
forward 5
forward 2
down 9
down 8
forward 8
down 7
up 6
forward 1
up 9
up 3
forward 9
down 6
forward 9
down 3
down 3
forward 7
forward 5
down 8
down 9
down 3
down 6
up 8
down 9
forward 8
down 7
down 5
down 1
up 4
down 9
forward 6
forward 9
up 6
up 4
forward 3
forward 2
forward 1
down 1
down 2
forward 8
up 6
forward 5
up 4
down 1
forward 5
down 3
down 6
up 7
forward 2
forward 6
forward 7
forward 4
down 5
down 4
forward 4
down 6
up 2
up 2
forward 7
forward 3
down 8
down 1
down 8
forward 7
forward 7
up 5
forward 4
up 8
down 9
down 4
down 4
forward 5
down 1
forward 2
down 6
up 4
down 8
down 1
down 9
down 5
up 5
forward 4
down 2
down 8
down 4
forward 4
forward 5
down 8
up 9
forward 7
forward 6
down 8
down 3
up 7
down 7
forward 2
forward 5
forward 7
down 9
up 1
down 6
down 2
forward 6
forward 3
forward 3
up 9
forward 4
down 5
down 7
forward 8
forward 6
forward 5
down 9
down 5
down 1
down 7
forward 9
forward 8
down 2
down 4
down 1
up 5
up 5
forward 5
down 3
down 1
forward 8
up 9
up 3
down 3
up 3
up 5
forward 8
down 3
up 3
down 9
up 6
up 8
forward 5
up 2
down 6
forward 3
down 2
down 4
forward 9
forward 6
forward 3
up 5
down 9
down 7
forward 9
forward 7
forward 5
up 5
up 1
down 6
forward 4
forward 4
down 7
down 1
up 3
forward 6
forward 4
down 1
forward 5
forward 3
forward 1
forward 3
up 3
up 9
down 7
down 4
forward 8
down 8
down 3
up 2
down 8
forward 5
down 7
forward 6
down 9
up 5
forward 4
down 2
forward 6
down 8
down 7
forward 8
forward 5
down 2
forward 7
forward 5
forward 7
down 8
forward 5
down 8
down 6
down 7
down 9
forward 9
down 6
forward 8
up 6
up 1
down 5
forward 1
forward 7
up 2
up 5
up 6
down 5
down 5
forward 7
down 9
down 2
forward 9
forward 3
down 5
up 2
up 8
forward 5
forward 8
up 1
forward 3
forward 1
up 4
forward 1
down 9
down 6
forward 1
down 4
down 4
forward 9
down 3
up 6
down 3
forward 6
forward 6
down 3
forward 6
down 3
down 1
forward 3
down 7
up 9
forward 1
down 7
down 2
up 8
down 1
down 9
down 1
down 4
down 6
down 3
down 7
down 2
down 9
down 2
forward 4
up 3
down 4
up 4
down 1
forward 5
forward 7
down 7
forward 9
forward 6
down 8
forward 6
forward 7
up 3
down 3
up 6
forward 7
up 4
forward 4
down 1
up 8
forward 7
down 2
up 6
forward 1
forward 3
up 9
up 8
up 5
forward 7
up 5
down 6
forward 7
forward 7
down 4
down 3
forward 2
down 8
up 9
up 6
forward 7
forward 5
down 9
down 2
up 5
down 3
down 3
up 5
down 8
forward 7
down 4
down 2
up 9
down 5
down 8
down 5
down 6
forward 9
down 3
down 8
forward 3
down 1
down 9
forward 1
down 3
up 9
up 3
forward 8
up 2
down 4
up 5
up 4
down 9
down 5
up 3
forward 2
down 8
forward 8
forward 7
up 4
down 9
down 6
up 1
forward 9
up 8
forward 4
up 3
down 4
up 2
up 7
down 2
forward 3
down 8
down 9
up 7
up 8
forward 3
forward 1
forward 7
forward 5
forward 9
forward 2
up 1
down 1
up 4
forward 1
up 9
forward 7
forward 2
down 6
down 5
forward 9
forward 4
down 6
down 6
up 8
down 3
up 8
down 3
forward 2
down 1
down 1
forward 5
down 1
forward 9
up 8
forward 2
down 5
up 8
up 8
forward 8
forward 8
forward 3
forward 2
forward 8
forward 9
forward 8
forward 6
forward 4
up 7
forward 9
forward 8
up 7
forward 6
forward 9
forward 8
down 7
forward 9
down 4
down 1
up 1
up 9
forward 2
down 6
down 2
down 8
down 6
up 8
forward 7
up 9
forward 5
forward 4
forward 8
up 4
forward 4
up 6
forward 7
forward 1
up 8
down 6
forward 7
forward 3
forward 2
down 4
forward 4
down 7
down 6
down 2
up 3
up 5
down 7
down 9
up 8
down 1
up 1
down 8
up 8
forward 8
down 6
down 1
down 6
forward 3
down 9
down 5
up 3
down 1
down 1
forward 4
down 4
up 3
forward 8
up 4
down 3
down 5
down 3
forward 6
forward 3
down 2
forward 9
forward 3
forward 2
down 2
forward 6
down 1
down 1
forward 5
forward 4
forward 6
down 7
forward 7
forward 3
forward 1
up 3
down 6
forward 1
up 9
forward 9
forward 5
forward 3
forward 3
down 3
up 8
forward 5
up 6
forward 2
down 7
forward 2
forward 8
forward 8
forward 3
up 9
down 5
down 3
forward 7
up 9
forward 4
down 1
down 3
down 5
down 2
forward 9
up 6
down 3
down 7
down 3
up 6
forward 3
down 4
forward 2
down 8
down 2
forward 7
down 2
down 9
forward 1
down 1
down 9
down 6
forward 5
down 1
up 1
forward 5
forward 4
forward 9
down 3
forward 3
forward 5
down 9
forward 9
down 8
down 2
forward 1
up 1
down 5
forward 2
up 9
forward 9
forward 7
forward 9
forward 3
down 7
forward 2
down 4
up 3
down 7
down 6
forward 2
down 2
forward 8
up 9
down 1
forward 7
down 8
forward 3
down 2
down 5
down 5
down 3
forward 1
up 9
up 9
down 8
down 6
up 7
forward 7
down 4
forward 6
down 9
up 5
up 6
forward 4
forward 1
forward 1
down 7
down 8
down 2
down 4
down 3
up 8
down 3
forward 3
forward 8
up 3
down 2
forward 4
down 3
forward 5
up 1
down 9
down 1
down 4
forward 3
forward 6
forward 7
forward 2
forward 9
forward 1
forward 6
forward 7
forward 2
up 1
down 6
down 1
forward 6
down 6
down 5
forward 1"""

"""
    parse_commands(commands<:AbstractString) ::Vector{Vector{<:AbstractString}}

Get the list of commands and their parameter.
"""
function parse_commands(comms::AbstractString) ::Vector{Vector{AbstractString}}
    return [split(x, " ") for x in split(comms, "\n")]
end

"""
    execute_commands(comms::Vector{Vector{AbstractString}}) ::Vector{Integer}

Get the final position of the sub.
"""
function execute_commands(comms::Vector{Vector{AbstractString}}) ::Vector{Integer}
    # Add third item for aim.
    position::Vector{Integer} = [0, 0, 0]
    for (command, param) in comms
        param = parse(Int32, param)
        if (command == "forward")
            position[1] += param
            position[2] += (position[3] * param)
        # up and down now change aim.
        elseif (command == "up")
            position[3] -= param
        elseif (command == "down")
            position[3] += param
        end
    end
    return position[1:2]
end

commands |> parse_commands |> execute_commands |> prod


text:
@doc raw"""
    GeneralLinear{n,𝔽} <:
        AbstractGroupManifold{𝔽,MultiplicationOperation,DefaultEmbeddingType}

The general linear group, that is, the group of all invertible matrices in ``𝔽^{n×n}``.

The default metric is the left-``\mathrm{GL}(n)``-right-``\mathrm{O}(n)``-invariant metric
whose inner product is
```math
⟨X_p,Y_p⟩_p = ⟨p^{-1}X_p,p^{-1}Y_p⟩_\mathrm{F} = ⟨X_e, Y_e⟩_\mathrm{F},
```
where ``X_p, Y_p ∈ T_p \mathrm{GL}(n, 𝔽)``,
``X_e = p^{-1}X_p ∈ 𝔤𝔩(n) = T_e \mathrm{GL}(n, 𝔽) = 𝔽^{n×n}`` is the corresponding
vector in the Lie algebra, and ``⟨⋅,⋅⟩_\mathrm{F}`` denotes the Frobenius inner product.

By default, tangent vectors ``X_p`` are represented with their corresponding Lie algebra
vectors ``X_e = p^{-1}X_p``.
"""
struct GeneralLinear{n,𝔽} <:
       AbstractGroupManifold{𝔽,MultiplicationOperation,DefaultGroupDecoratorType} end

GeneralLinear(n, 𝔽::AbstractNumbers=ℝ) = GeneralLinear{n,𝔽}()

function allocation_promotion_function(::GeneralLinear{n,ℂ}, f, ::Tuple) where {n}
    return complex
end

function check_point(G::GeneralLinear, p; kwargs...)
    mpv = check_point(decorated_manifold(G), p; kwargs...)
    mpv === nothing || return mpv
    detp = det(p)
    if iszero(detp)
        return DomainError(
            detp,
            "The matrix $(p) does not lie on $(G), since it is not invertible.",
        )
    end
    return nothing
end
check_point(::GT, ::Identity{GT}; kwargs...) where {GT<:GeneralLinear} = nothing
function check_point(G::GeneralLinear, e::Identity; kwargs...)
    return DomainError(e, "The identity element $(e) does not belong to $(G).")
end

function check_vector(G::GeneralLinear, p, X; kwargs...)
    mpv = check_vector(decorated_manifold(G), p, X; kwargs...)
    mpv === nothing || return mpv
    return nothing
end

decorated_manifold(::GeneralLinear{n,𝔽}) where {n,𝔽} = Euclidean(n, n; field=𝔽)

default_metric_dispatch(::GeneralLinear, ::EuclideanMetric) = Val(true)
default_metric_dispatch(::GeneralLinear, ::LeftInvariantMetric{EuclideanMetric}) = Val(true)

distance(G::GeneralLinear, p, q) = norm(G, p, log(G, p, q))

@doc raw"""
    exp(G::GeneralLinear, p, X)

Compute the exponential map on the [`GeneralLinear`](@ref) group.

The exponential map is
````math
\exp_p \colon X ↦ p \operatorname{Exp}(X^\mathrm{H}) \operatorname{Exp}(X - X^\mathrm{H}),
````

where ``\operatorname{Exp}(⋅)`` denotes the matrix exponential, and ``⋅^\mathrm{H}`` is
the conjugate transpose. [^AndruchowLarotondaRechtVarela2014][^MartinNeff2016]

[^AndruchowLarotondaRechtVarela2014]:
    > Andruchow E., Larotonda G., Recht L., and Varela A.:
    > “The left invariant metric in the general linear group”,
    > Journal of Geometry and Physics 86, pp. 241-257, 2014.
    > doi: [10.1016/j.geomphys.2014.08.009](https://doi.org/10.1016/j.geomphys.2014.08.009),
    > arXiv: [1109.0520v1](https://arxiv.org/abs/1109.0520v1).
[^MartinNeff2016]:
    > Martin, R. J. and Neff, P.:
    > “Minimal geodesics on GL(n) for left-invariant, right-O(n)-invariant Riemannian metrics”,
    > Journal of Geometric Mechanics 8(3), pp. 323-357, 2016.
    > doi: [10.3934/jgm.2016010](https://doi.org/10.3934/jgm.2016010),
    > arXiv: [1409.7849v2](https://arxiv.org/abs/1409.7849v2).
"""
exp(::GeneralLinear, p, X)

function exp!(G::GeneralLinear, q, p, X)
    expX = exp(X)
    if isnormal(X; atol=sqrt(eps(real(eltype(X)))))
        return compose!(G, q, p, expX)
    end
    compose!(G, q, expX', exp(X - X'))
    compose!(G, q, p, q)
    return q
end
function exp!(::GeneralLinear{1}, q, p, X)
    p1 = p isa Identity ? p : p[1]
    q[1] = p1 * exp(X[1])
    return q
end
function exp!(G::GeneralLinear{2}, q, p, X)
    if isnormal(X; atol=sqrt(eps(real(eltype(X)))))
        return compose!(G, q, p, exp(SizedMatrix{2,2}(X)))
    end
    A = SizedMatrix{2,2}(X')
    B = SizedMatrix{2,2}(X) - A
    compose!(G, q, exp(A), exp(B))
    compose!(G, q, p, q)
    return q
end

function get_coordinates(
    ::GeneralLinear{n,ℝ},
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {n}
    return vec(X)
end

function get_coordinates!(
    ::GeneralLinear{n,ℝ},
    Xⁱ,
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {n}
    return copyto!(Xⁱ, X)
end

function get_vector(
    ::GeneralLinear{n,ℝ},
    p,
    Xⁱ,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {n}
    return reshape(Xⁱ, n, n)
end

function get_vector!(
    ::GeneralLinear{n,ℝ},
    X,
    p,
    Xⁱ,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {n}
    return copyto!(X, Xⁱ)
end

function group_exp!(::GeneralLinear{1}, q, X)
    q[1] = exp(X[1])
    return q
end
group_exp!(::GeneralLinear{2}, q, X) = copyto!(q, exp(SizedMatrix{2,2}(X)))

function group_log!(::GeneralLinear{1}, X::AbstractMatrix, p::AbstractMatrix)
    X[1] = log(p[1])
    return X
end

inner(::GeneralLinear, p, X, Y) = dot(X, Y)

invariant_metric_dispatch(::GeneralLinear, ::LeftAction) = Val(true)

inverse_translate_diff(::GeneralLinear, p, q, X, ::LeftAction) = X
inverse_translate_diff(::GeneralLinear, p, q, X, ::RightAction) = p * X / p

function inverse_translate_diff!(G::GeneralLinear, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, inverse_translate_diff(G, p, q, X, conv))
end

# find sU for s ∈ S⁺ and U ∈ U(n, 𝔽) that minimizes ‖sU - p‖²
function _project_Un_S⁺(p)
    n = LinearAlgebra.checksquare(p)
    F = svd(p)
    s = mean(F.S)
    U = F.U * F.Vt
    return rmul!(U, s)
end

@doc raw"""
    log(G::GeneralLinear, p, q)

Compute the logarithmic map on the [`GeneralLinear(n)`](@ref) group.

The algorithm proceeds in two stages. First, the point ``r = p^{-1} q`` is projected to the
nearest element (under the Frobenius norm) of the direct product subgroup
``\mathrm{O}(n) × S^+``, whose logarithmic map is exactly computed using the matrix
logarithm. This initial tangent vector is then refined using the
[`NLsolveInverseRetraction`](@ref).

For `GeneralLinear(n, ℂ)`, the logarithmic map is instead computed on the realified
supergroup `GeneralLinear(2n)` and the resulting tangent vector is then complexified.

Note that this implementation is experimental.
"""
log(::GeneralLinear, p, q)

function log!(G::GeneralLinear{n,𝔽}, X, p, q) where {n,𝔽}
    pinvq = inverse_translate(G, p, q, LeftAction())
    𝔽 === ℝ && det(pinvq) ≤ 0 && throw(OutOfInjectivityRadiusError())
    e = Identity(G, pinvq)
    if isnormal(pinvq; atol=sqrt(eps(real(eltype(pinvq)))))
        log_safe!(X, pinvq)
    else
        # compute the equivalent logarithm on GL(dim(𝔽) * n, ℝ)
        # this is significantly more stable than computing the complex algorithm
        Gᵣ = GeneralLinear(real_dimension(𝔽) * n, ℝ)
        pinvqᵣ = realify(pinvq, 𝔽)
        Xᵣ = realify(X, 𝔽)
        eᵣ = Identity(Gᵣ, pinvqᵣ)
        log_safe!(Xᵣ, _project_Un_S⁺(pinvqᵣ))
        inverse_retraction = NLsolveInverseRetraction(ExponentialRetraction(), Xᵣ)
        inverse_retract!(Gᵣ, Xᵣ, eᵣ, pinvqᵣ, inverse_retraction)
        unrealify!(X, Xᵣ, 𝔽, n)
    end
    translate_diff!(G, X, p, e, X, LeftAction())
    return X
end
function log!(::GeneralLinear{1}, X, p, q)
    p1 = p isa Identity ? p : p[1]
    X[1] = log(p1 \ q[1])
    return X
end

manifold_dimension(G::GeneralLinear) = manifold_dimension(decorated_manifold(G))

LinearAlgebra.norm(::GeneralLinear, p, X) = norm(X)

project(::GeneralLinear, p) = p
project(::GeneralLinear, p, X) = X

project!(::GeneralLinear, q, p) = copyto!(q, p)
project!(::GeneralLinear, Y, p, X) = copyto!(Y, X)

Base.show(io::IO, ::GeneralLinear{n,𝔽}) where {n,𝔽} = print(io, "GeneralLinear($n, $𝔽)")

translate_diff(::GeneralLinear, p, q, X, ::LeftAction) = X
translate_diff(::GeneralLinear, p, q, X, ::RightAction) = p \ X * p

function translate_diff!(G::GeneralLinear, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, translate_diff(G, p, q, X, conv))
end

vector_transport_to(::GeneralLinear, p, X, q, ::ParallelTransport) = X

vector_transport_to!(::GeneralLinear, Y, p, X, q, ::ParallelTransport) = copyto!(Y, X)


text:
module SinKing

include("./types.jl")
include("./network.jl")
include("./signals.jl")

module AgentParts
include("./agent_parts/iz_neuron.jl")
# include("./agent_parts/conductance_injection.jl")
# include("./agent_parts/markram_transmitter.jl")
include("./agent_parts/lif_neuron.jl")
include("./agent_parts/dc.jl")
end

module Agents
#include("./agents/iz_cond_mrkrm.jl")
include("./agents/lif_simple.jl")
include("./agents/dc_source.jl")
include("./agents/delayer.jl")
end

end # module


text:
list_test_str = split(
"",
'\n'
)

list_str = let input = @__DIR__() * "/day_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list)
end

testit("Test Part 1", part1(list_test_str), ==(nothing))
testit("Challenge Part 1", part1(list_str), ==(nothing))

# Part 2

function part2(list)
end

testit("Test Part 2", part2(list_test_str), ==(nothing))
testit("Challenge Part 2", part2(list_str), ==(nothing))


text:
using BinaryBuilder

name = "LibOSXUnwind"
version = v"0.0.6"

# Collection of sources required to build libosxunwind
sources = [
    ArchiveSource("https://github.com/JuliaLang/libosxunwind/archive/v$(version).tar.gz",
                  "c3943fdd063c2c8c249778326cbfecdc62c804a0509c4ae4799604e1817f9058"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libosxunwind*/

EXTRA_CFLAGS="-ggdb3 -O2"

FLAGS=(
    CC="$CC"
    CXX="$CXX"
    CFLAGS="${CFLAGS} -ggdb3 -O2"
    CXXFLAGS="${CXXFLAGS} -ggdb3 -O2"
    PREFIX="${prefix}"
)

# When all you have is a hammer...
make -j${nproc} "${FLAGS[@]}"

# Manual installation as the osxunwind `Makefile` doesnt' even know how to do this
mkdir -p ${libdir}
cp libosxunwind.dylib ${libdir}/
cp libosxunwind.a ${libdir}/
cp -aR include ${prefix}/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(Sys.isdarwin, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libosxunwind", :libosxunwind),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil