text:
###############################################################################
#
#  Isomorphisms
#
#
#  (C) 2022 Jeroen Hanselman
#
###############################################################################

export Isogeny

export isomorphism, is_isomorphic, isomorphism_data, isomorphism_to_isogeny, identity_map,
negation_map, automorphism_group_generators, rational_maps, transform_rstu

###############################################################################
#
#  Isomorphisms
#
###############################################################################

mutable struct EllCrvIso{T} <: Map{EllCrv, EllCrv, HeckeMap, EllCrvIso} where T<:RingElem
  header::MapHeader{EllCrv{T}, EllCrv{T}}
  domain::EllCrv{T}
  codomain::EllCrv{T}
  data::Tuple{T, T, T, T}

  function EllCrvIso(E::EllCrv{T}, data::Vector{T}) where T <:RingElem
    f = new{T}()
    f.domain = E

    if length(data)!= 4
      throw(DomainError(T, "Array needs to have length 4"))
    end
    r, s, t, u = data[1], data[2], data[3], data[4]
    f.data = (r, s, t, u)
    a1, a2, a3, a4, a6 = a_invars(E)

    a1new = divexact(a1 + 2*s, u)
    a2new = divexact(a2 - s*a1 + 3*r - s^2, u^2)
    a3new = divexact(a3 + r*a1 + 2*t, u^3)
    a4new = divexact(a4 - s*a3 + 2*r*a2 - (t + r*s)*a1 + 3*r^2 - 2*s*t, u^4)
    a6new = divexact(a6 + r*a4 + r^2*a2 + r^3 - t*a3 - t^2 - r*t*a1, u^6)

    F = EllipticCurve([a1new, a2new, a3new, a4new, a6new])
    f.codomain = F
    f.header = MapHeader(E, F)
    return f
  end

end

@doc Markdown.doc"""
    identity_map(E::EllCrv) -> EllCrvIso

Return the identity isomorphism on the elliptic curve $E$.
"""
function identity_map(E::EllCrv)
  return isomorphism(E, [0, 0, 0, 1])
end

@doc Markdown.doc"""
    negation_map(E::EllCrv) -> EllCrvIso

Return the negation isomorphism on the elliptic curve $E$.
"""
function negation_map(E::EllCrv)
  a1, a2, a3, a4, a6 = a_invars(E)
  return isomorphism(E, [0, -a1, -a3, -1])
end

@doc Markdown.doc"""
    isomorphism_data(f::EllCrvIso) -> (T, T, T, T)

Return the tuple $(r, s, t, u)$ that defines the isomorphism which is of the form
[x : y : 1] -> [(x - r)//u^2 : (y - s*(x-r) - t)//u^3 : 1]
"""
function isomorphism_data(f::EllCrvIso)
  return f.data
end

@doc Markdown.doc"""
    inv(f::EllCrvIso) -> EllCrvIso
Return the inverse of the isomorphism $f$.
"""
function inv(f::EllCrvIso{T}) where T<:RingElem
  r, s, t, u = isomorphism_data(f)
  newr = -r//u^2
  news = -s//u
  newt = (r*s-t)//u^3
  newu = 1//u

  g = EllCrvIso(codomain(f),[newr, news, newt, newu])
  return g
end

@doc Markdown.doc"""
    isomorphism_to_isogeny(f::EllCrvIso) -> Isogeny
Return the isomorphism $f$ as an isogeny.
"""
function isomorphism_to_isogeny(f::EllCrvIso)
  r, s, t, u = f.data
  E = f.domain
  phi = Isogeny(E)
  phi.codomain = f.codomain
  phi.degree = 1

  R = base_field(E)
  Rx, x = PolynomialRing(R, "x")
  phi.psi = one(Rx)

  phi.coordx = (x - r)//Rx(u^2)
  Rxy, y = PolynomialRing(Rx, "y")
  phi.coordy = (y - s*(x-r) - t)//Rxy(u^3)
  phi.header = f.header
  return phi
end

function compose(f::Isogeny, g::EllCrvIso)
  return compose(f, isomorphism_to_isogeny(g))
end

function compose(g::EllCrvIso, f::Isogeny)
  return compose(isomorphism_to_isogeny(g), f)
end

#Following Connell's Handbook for Elliptic Curves Chapter 4.4
@doc Markdown.doc"""
    is_isomorphic(E1::EllCrv, E2::EllCrv) -> Bool
Return true when $E1$ and $E2$ are isomorphic
"""
function is_isomorphic(E1::EllCrv{T}, E2::EllCrv{T}) where T
  K = base_field(E1)
  char = characteristic(K)
  j1 = j_invariant(E1)
  j2 = j_invariant(E2)

  #First check the j-invariants
  if j1 != j2
    return false
  end

  if char == 2
    E1, phi1 = simplified_model(E1)
    E2, phi2 = simplified_model(E2)
    a1, a2, a3, a4, a6 = a_invars(E1)
    _a1, _a2, _a3, _a4, _a6 = a_invars(E2)
    Rx, x = PolynomialRing(K, "x")
    # j-invariant non-zero
    if j1!=0
      f = x^2 + x + a2 + _a2
      return !isempty(roots(f))
    end
    # j-invariant is 0
    us = roots(x^3 - a3//_a3)
    for u in us
      g = x^4 + a3*x + a4 + u^4*_a4
      ss = roots(g)
      for s in ss
        h = x^2 + a3*x + s^6 + a4*s^2 + a6 + u^6*_a6
        ts = roots(h)
        if !isempty(ts)
          return true
        end
      end
    end
    return false
  end

  if char == 3 && j1 == 0
    E1, phi1 = simplified_model(E1)
    E2, phi2 = simplified_model(E2)
    a1, a2, a3, a4, a6 = a_invars(E1)
    _a1, _a2, _a3, _a4, _a6 = a_invars(E2)
    Rx, x = PolynomialRing(K, "x", cached = false)
    us = roots(x^4 - a4//_a4)
    for u in us
      g = x^3 + a4*x + a6 - u^6*_a6
      rs = roots(g)
      if !isempty(rs)
        return true
      end
    end
    return false
  end

  c4, c6 = c_invars(E1)
  _c4, _c6 = c_invars(E2)

  if j1!=0 && j1!=1728
    return issquare(c6//_c6)
  else
    Rx, x = PolynomialRing(K, "x")
    if j1 == 1728
      us = roots(x^4 - c4//_c4)
      return !isempty(us)
    end
    if j1 == 0
      us = roots(x^4 - c6//_c6)
      return !isempty(us)
    end
  end
  error("There is a bug in is_isomorphism")
end

#Following Connell's Handbook for Elliptic Curves Chapter 4.4
@doc Markdown.doc"""
    isomorphism(E1::EllCrv, E2::EllCrv) -> EllCrvIso
Return an isomorphism between $E1$ and $E2$ if they are isomorphic.
Otherwise return an error.
"""
function isomorphism(E1::EllCrv, E2::EllCrv)
  K = base_field(E1)
  char = characteristic(K)
  j1 = j_invariant(E1)
  j2 = j_invariant(E2)

  #First check the j-invariants
  j1 == j2 || error("Curves are not isomorphic")

  E1s, pre_iso = simplified_model(E1)
  E2s, _, post_iso = simplified_model(E2)
  a1, a2, a3, a4, a6 = a_invars(E1s)
  _a1, _a2, _a3, _a4, _a6 = a_invars(E2s)

  if char == 2
    Rx, x = PolynomialRing(K, "x")
    # j-invariant non-zero
    if j1!=0
      f = x^2 + x + a2 + _a2
      ss = roots(f)
      !isempty(ss) || error("Curves are not isomorphic")

      s = ss[1]
      phi = isomorphism(E1s, [0, s, 0, 1])
      F = codomain(phi)
      @assert F == E2s "There is a bug in isomorphism"
      return pre_iso*phi * post_iso
    end
    # j-invariant is 0
    us = roots(x^3 - a3//_a3)
    for u in us
      g = x^4 + a3*x + a4 + u^4*_a4
      ss = roots(g)
      for s in ss
        h = x^2 + a3*x + s^6 + a4*s^2 + a6 + u^6*_a6
        ts = roots(h)
        !isempty(ts) || error("Curves are not isomorphic")

        t = ts[1]
        phi = isomorphism(E1s, [s^2, s, t, u])
        F = codomain(phi)
        @assert F == E2s "There is a bug in isomorphism"
        return pre_iso * phi * post_iso
      end
    end
  end

  if char == 3 && j1 == 0
    Rx, x = PolynomialRing(K, "x")
    us = roots(x^4 - a4//_a4)
    for u in us
      g = x^3 + a4*x + a6 - u^6*_a6
      rs = roots(g)
      !isempty(rs) || error("Curves are not isomorphic")

      r = rs[1]
      phi = isomorphism(E1s, [r, 0, 0, u])
      F = codomain(phi)
      @assert F == E2s "There is a bug in isomorphism"
      return pre_iso * phi * post_iso
    end
  end

  #Characteristic != 2, 3 and characteristic 3 with j!= 0
  c4, c6 = c_invars(E1)
  _c4, _c6 = c_invars(E2)
  usq = (c6//_c6)//(c4//_c4)

  issquare(usq) || error("Curves are not isomorphic")
  u = sqrt(usq)
  phi = isomorphism(E1s, [0, 0, 0, u])
  F = codomain(phi)
  @assert F == E2s "There is a bug in isomorphism"
  return pre_iso * phi * post_iso
end
#=
function isomorphism(E::EllCrv, E2::EllCrv)
  char = characteristic(base_field(E))
  if char!= 2 && char!= 3
    if is_isomorphic(E, E2)
      a1, a2, a3 = a_invars(E)
      _a1, _a2, _a3 = a_invars(E2)

      c4, c6 = c_invars(E)
      _c4, _c6 = c_invars(E2)
      usq = (c6//_c6)//(c4//_c4)

      u = sqrt(usq)
      s = (_a1*u-a1)//2
      r = (_a2*u^2-a2+s*a1+s^2)//3
      t = (_a3*u^3-a3-r*a1)//2
      return isomorphism(E,[r,s,t,u])
    else
      throw(DomainError(E, "Curves are not isomorphic over the base field"))
    end
  else
    throw(DomainError(E, "Isomorphism test only implemented for characteristic not equal to 2 or 3"))
  end
end
=#

################################################################################
#
# T(r, s, t, u) transformation
#
################################################################################

# transformation T(r,s,t,u) as in Connell's handbook
@doc Markdown.doc"""
    transform_rstu(E::EllCrv, [r, s, t, u]::Vector{T})
    -> EllCrv, EllCrvIso, EllCrvIso
Return the transformation of E under the isomorphism given by
[(x - r)//u^2 : (y - s*(x-r) - t)//u^3 : 1], the isomorphism and the inverse of
the isomorphism
"""
function transform_rstu(E::EllCrv, T::Vector{S}) where S

  phi = isomorphism(E, T)

  return codomain(phi), phi, inv(phi)
end

@doc Markdown.doc"""
    isomorphism(E::EllCrv, [r, s, t, u]::Vector{T}) -> EllCrvIso
Return the isomorphism with domain E given by
[(x - r)//u^2 : (y - s*(x-r) - t)//u^3 : 1]. The codomain
is calculated automatically.
"""
function isomorphism(E::EllCrv{T}, data::Vector{T}) where T

  if length(data)!= 4
    throw(DomainError(data, "Array needs to have length 4"))
  end
  return EllCrvIso(E, data)
end

function isomorphism(E::EllCrv, data::Vector)

  if length(data)!= 4
    throw(DomainError(data, "Array needs to have length 4"))
   end

  K = base_field(E)
  data = map(K, data)
  return EllCrvIso(E, data)
end

function degree(f::EllCrvIso)
  return 1
end

@doc Markdown.doc"""
    image(f::EllCrvIso, P::EllCrvPt) -> EllCrvPt
Return the image of $P$ under the isomorphism $f$.
"""
function image(f::EllCrvIso, P::EllCrvPt)
  @assert domain(f) == parent(P)
  F = codomain(f)
  if !is_finite(P)
    return infinity(F)
  end
  r, s, t, u = isomorphism_data(f)
  xnew = divexact(1, u^2)*(P.coordx - r)
  ynew = divexact(1, u^3)*(P.coordy - s*u^2*xnew - t)
  Q = F([xnew, ynew])
  return Q
end

@doc Markdown.doc"""
    preiimage(f::EllCrvIso, P::EllCrvPt) -> EllCrvPt
Return the preimage of $P$ under the isomorphism $f$.
"""
function preimage(f::EllCrvIso, P::EllCrvPt)
  @assert codomain(f) == parent(P)
  E = domain(f)
  if !is_finite(P)
    return infinity(E)
  end
  r, s, t, u = isomorphism_data(f)
  xnew = u^2*P.coordx + r
  ynew = u^3*P.coordy + s*u^2*P.coordx + t
  Q = E([xnew, ynew])
  return Q
end

function compose(f::EllCrvIso, g::EllCrvIso)
  @req codomain(f) == domain(g) "Codomain and domain must be identitcal"
  r1, s1, t1, u1 = f.data
  r2, s2, t2, u2 = g.data

  r3 = r1 + u1^2*r2
  s3 = s1 + u1*s2
  t3 = t1 + u1^2*s1*r2 + u1^3*t2
  u3 = u1*u2

  return isomorphism(domain(f), [r3, s3, t3, u3])

end

@doc Markdown.doc"""
    automorphism_group_generators(E::EllCrv) -> Vector{EllCrvIso}
Return generators of the automorphism group of $E$.
"""
function automorphism_group_generators(E::EllCrv{T}) where {T}
  K = base_field(E)
  char = characteristic(K)
  j = j_invariant(E)

  #Group is Z/2Z
  if j!= 0 && j!= 1728
    aut = abelian_group(2)
    return [negation_map(E)]
  end

  Kx, x = PolynomialRing(K, cached = false)
  Es, pre_iso, post_iso = simplified_model(E)
  a1, a2, a3, a4, a6 = a_invars(Es)

  if char != 2 && char != 3
    if j == 1728
      f = x^2+1
      a = roots(f)
      if !isempty(roots(f))
        #Group is Z/4Z
        i = a[1]
        return [pre_iso * isomorphism(Es, [0, 0, 0, i]) * post_iso]
      end
    elseif j == 0
      f = x^2+x+1
      a = roots(f)
      if !isempty(roots(f))
        #Group is Z/6Z
        zeta3 = a[1]
        return [pre_iso * isomorphism(Es, [0, 0, 0, -zeta3]) * post_iso]
      end
    end
    #Group is Z/2Z
    return [negation_map(E)]
  end

  if char == 3 #And j-invariant is 0.
    us = roots(x^2 + 1)  #See if x^4 + 1 = (x^2 +1)*(x^2 -1) has a root that induces an element of order 4
    if !isempty(us)
      rs = roots(x^3 + a4*x + 2*a6)
      i = us[1]
      if !isempty(rs)
        # Group is dicyclic group of order 12. <a, b | a^6 = b^4 = id, a^3 = b^2, bab^(-1) = a^(-1)>
        r = rs[1]
        return [pre_iso * isomorphism(Es, [r, 0, 0, -1]) * post_iso, pre_iso * isomorphism(Es, [r, 0, 0, i]) * post_iso]
      else
        #Group is Z/4Z
        return [pre_iso * isomorphism(Es, [0, 0, 0, i]) * post_iso]
      end
    else
      rs = roots(x^2 + a4) #Now u^2 = 1, so x^3 + a4*x + a6 - u^6*a6 = x*(x^2+a4)
      if !isempty(rs)
        #Group is Z/6Z
        r = rs[1]
        return [pre_iso * isomorphism(Es, [0, 0, 0, -r]) * post_iso]
      end
    end
    #Group is Z/2Z
    return [negation_map(E)]
  end

  if char == 2 #And j-invariant is 0
    auts = EllCrvIso{T}[]
    us = roots(x^2 + x +1) #Check if there are us other than u = 1.
    if !isempty(us)
      u = us[1]
      g = x^4 + a3*x + (1+u)*a4
      ss = roots(g)
      for s in ss
        h = x^2 + a3*x + s^6 + a4*s^2
        ts = roots(h)
        for t in ts
          push!(auts, pre_iso * isomorphism(Es, [s^2, s, t, u]) * post_iso)
        end
      end
    end
    size = length(auts)
    if size == 8 #Group is SL(2,3)
      #Search for generators. One element of order 3 and one element of order 4 should suffice.
      #Element of order 3
      for phi in auts
        phi3 = phi * phi * phi
        if phi3 == identity_map(Es)
          g1 = phi
        end
      end
      #Element of order 4. Need to take u = 1.
      g = x^3 + a3*x
      s = roots(g)[1]
      h = x^2 + a3*x + s^6 + a4*s^2
      t = roots(h)[1]

      g2 = pre_iso * isomorphism(Es, [s^2, s, t, 1]) * post_iso

      return [g1, g2]
    elseif size == 2 #Group is Z/6Z
      g1 = auts[1]
      g2 = auts[2]
      if (g1 * g1 * g1) != identity_map(Es)
        return g1
      else
        return g2
      end
    end
  else
      #u = 1
    auts = EllCrvIso{T}[]
    g = x^3 + a3
    ss = roots(g)
    for s in ss
      h = x^2 + a3*x + s^6 + a4*s^2
      ts = roots(h)
      for t in ts
        push!(auts, pre_iso * isomorphism(Es, [s^2, s, t, 1]) * post_iso)
      end
    end
    size = length(auts)
    if size == 6 #Group isomorphic to Quaternions
      if auts[1] == inv(auts[2])
        return [auts[1], auts[3]]
      else
        return [auts[1], auts[2]]
      end
    elseif size == 2 # Group isomorphic to Z/4Z
      return [auts[1]]
    end
  #Group is Z/2Z
  return [negation_map(E)]
  end
end

@doc Markdown.doc"""
    rational_maps(f::EllCrvIso) -> Vector{Poly}
Return the rational maps defining the isomorphism.
"""
function rational_maps(f::EllCrvIso)
  K = base_field(domain(f))
  Kx, x = PolynomialRing(K, "x")
  Kxy, y = PolynomialRing(Kx, "y")
  r, s, t, u = f.data
  xnew = divexact(1, u^2)*(x - r)
  ynew = divexact(1, u^3)*(y - s*u^2*xnew - t)
  return [xnew, ynew]
end

@doc Markdown.doc"""
    ==(f::EllCrvIso, g::EllCrvIso) -> Bool

Return true if $f$ and $g$ define the same map over the same field.
"""
function ==(f::EllCrvIso, g::EllCrvIso) where T
  Ef = domain(f)
  Eg = domain(g)
  return f.data == g.data && Ef == Eg && base_field(Ef) == base_field(Eg)
end

################################################################################
#
#  Automorphism group
#
################################################################################

# Need a parent object for the morphisms
struct EllCrvIsoSet{T}
  crv::T
end

parent(f::EllCrvIso{T}) where {T} = EllCrvIsoSet{EllCrv{T}}(domain(f))

mutable struct EllCrvAutMap{S, T} <: Map{GrpGen, EllCrvIsoSet{T}, HeckeMap, EllCrvAutMap}
  G::GrpGen
  auts::Vector{T}
  header::MapHeader{GrpGen, EllCrvIsoSet{S}}

  function EllCrvAutMap(E::S, G::GrpGen, auts::Vector{T}) where {S, T}
    return new{S, T}(G, auts, MapHeader(G, EllCrvIsoSet{S}(E)))
  end
end

function image(f::EllCrvAutMap, g::GrpGenElem)
  return f.auts[g[]]
end

function preimage(f::EllCrvAutMap, m::EllCrvIso)
  @req parent(m) == codomain(f) "Isomorphism has wrong elliptic curve"
  for g in domain(f)
    if f(g) == m
      return g
    end
  end
  error("This should not happen")
end

function automorphism_group(E::EllCrv)
  gens = automorphism_group_generators(E)
  cl = closure(gens, *)
  G, Gtocl, cltoG = generic_group(cl, *)
  m = EllCrvAutMap(E, G, cl)
  return domain(m), m
end


text:
# calculate required particle num for kld sampling
# by wilson-hilferty transformation

module ParticleNumWh
  using Distributions

  function num(epsilon, delta, bin_num)
    return ceil(quantile(Chisq(bin_num-1), 1-delta)/(2*epsilon))
  end

  function num_wh(epsilon, delta, bin_num)
    dof = bin_num - 1
    z = quantile(Normal(), 1 - delta)
    return ceil(dof/(2*epsilon)*(1-2/(9*dof)+sqrt(2/(9*dof))*z)^3)
  end

  function main()
    bins = [2, 4, 8, 1000, 10000, 100000]
    for bin in bins
      println("bin:$(bin), ϵ=0.1, δ=0.01, chisq=$(num(0.1, 0.01, bin)), norm=$(num_wh(0.1, 0.01, bin))")
      println("bin:$(bin), ϵ=0.5, δ=0.01, chisq=$(num(0.5, 0.01, bin)), norm=$(num_wh(0.5, 0.01, bin))")
      println("bin:$(bin), ϵ=0.5, δ=0.05, chisq=$(num(0.5, 0.05, bin)), norm=$(num_wh(0.5, 0.05, bin))")
    end
  end
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using Base: @nloops, @ntuple
using Einsum
using LinearAlgebra
using TensorOperations

function computemarginal(d::DiscreteSpinGlass)
    # compute the marginal probiliby for each variable to take +1
    n = length(d)
    exp = quote
        Z = 0.0
        mP = zeros($n) # marginal probability
        @nloops $n i (d)->[-1.0,1.0] begin
            x = @ntuple $n i
            x = collect(x)
            local p = exp(Distributions.logpdf($d,x))
            term = p.*x
            term[term.<0] .= 0.0
            global mP .+= term
            global Z += p
        end
        return mP./Z, Z
    end
    eval(exp)
end

function get_neighbours(d::Union{DiscreteSpinGlass, ThirdOrderExpFamily})
    """
    Get list of neighbor ids for both factor nodes and variable nodes
    """
    fidx = d.edge_index
    vidx = Array{Array{Int64,1},1}([[] for _ in 1:d.nnode])
    for i in 1:d.nfactor
        for j in fidx[i]
            push!(vidx[j], i)
        end
    end
    return fidx, vidx
end

function get_factor_tensor(d::DiscreteSpinGlass)
    """
    get factor tensor on support [-1,1]^r where r is the rank of the tensor
    """
    res = []
    for i in 1:d.nfactor
        idx = d.edge_index[i]
        supp = [1.; -1.]
        alpha = d.edge_attr[i][1] # assume scalar
        if length(idx) == 1
            @tensor begin
                T[i] := supp[i]
            end
        elseif length(idx) == 2
            @tensor begin
                T[i,j] := supp[i] * supp[j]
            end
        elseif length(idx) == 3
            @tensor begin
                T[i,j,k] := supp[i] * supp[j] * supp[k]
            end
        end
        T = @. exp(-alpha * T)
        T ./= maximum(T)
        push!(res,T)
    end
    return res
end

function bp_3spin(fidx::AbstractVector, vidx::AbstractVector, ftensor::AbstractVector; niter::Integer=1, atol=1e-5)
    nf = length(fidx)
    nv = length(vidx)
    ## Initialzie messages
    mv2f = Dict{Tuple,Array}()
    mf2v = Dict{Tuple,Array}()
    for i in 1:nf
        for j in fidx[i]
            mf2v[(i,j)] = ones(2) / 2.  # 2 state spin
            mv2f[(j,i)] = ones(2) / 2.
        end
    end

    ## main loop
    isconverged=false
    for iiter in 1:niter
        mv2f_new = deepcopy(mv2f)
        mf2v_new = deepcopy(mf2v)
        # update variable-to-factor messages
        for j in 1:nv
            source_idx = combinations(1:length(vidx[j]),length(vidx[j])-1) |> collect |> reverse
            for (neighbor_id, i) in enumerate(vidx[j])
                if length(source_idx) > 1
                    mv2f_new[(j,i)] = foldl((x,y)->x.*y,[mf2v[(vidx[j][k],j)] for k in source_idx[neighbor_id]])
                end
            end
        end
        # update factor-to-variable messages
        for i in 1:nf
            factor_size = length(fidx[i])
            source_idx = combinations(1:factor_size,factor_size-1) |> collect |> reverse
            for (neighbor_id, j) in enumerate(fidx[i])
                if factor_size == 1
                    mf2v_new[(i,j)] = ftensor[i]
                else
                    perm = collect(1:factor_size )
                    perm[neighbor_id] = 1
                    perm[1] = neighbor_id
                    ft_perm = permutedims(ftensor[i],perm)
                    if factor_size == 2
                        T = ft_perm
                        k1 = source_idx[neighbor_id][1]
                        m1 = mv2f[(fidx[i][k1],i)]
                        @einsum tmp[a] := T[a,b] * m1[b]
                        mf2v_new[(i,j)] .= tmp
                    elseif factor_size == 3
                        T = ft_perm
                        k1 = source_idx[neighbor_id][1]
                        k2 = source_idx[neighbor_id][2]
                        m1 = mv2f[(fidx[i][k1],i)]
                        m2 = mv2f[(fidx[i][k2],i)]
                        @einsum tmp[a] := T[a,b,c] * m1[b] * m2[c]
                        mf2v_new[(i,j)] .= tmp
                    end
                    @assert all(tmp .>= 0) "$tmp\n$(T) \n$(m1)\n$(m2)"
                end
            end
        end
        # normalize message
        for k in keys(mv2f)
            mv2f_new[k] = mv2f_new[k] ./ sum(mv2f_new[k])
        end
        for k in keys(mf2v)
            mf2v_new[k] = mf2v_new[k] ./ sum(mf2v_new[k])
        end
        # check convergence
        # assert convergent if all message change < atol
        isconverged = true
        for k in keys(mv2f)
            if abs(mv2f[k][1]-mv2f_new[k][1]) > atol
                isconverged = false
                break
            end
        end
        for k in keys(mf2v)
            if abs(mf2v[k][1]-mf2v_new[k][1]) > atol
                isconverged = false
                break
            end
        end
        mv2f = mv2f_new
        mf2v = mf2v_new
    end

    ## Calculate marginal distributions
    # singleton marginals
    marginals = []
    for j in 1:nv
        incoming_messages = [mf2v[(i,j)] for i in vidx[j]]
        if length(incoming_messages)>1
            mP = foldl((x,y)->x.*y,incoming_messages)
        else
            mP = ones(2) / .2
        end
        mP ./= sum(mP)
        push!(marginals,mP[1])
    end

    return marginals, isconverged, mf2v, mv2f
end

function bp(d::DiscreteSpinGlass; niter::Integer=10, atol::Real=1e-5)
    if !has_singleton_factor(d)
        d = add_singleton_factor(d)
    end
    fidx, vidx = get_neighbours(d)
    ftensor = get_factor_tensor(d)
    mP_bp, isconverged, mf2v, mv2f = bp_3spin(fidx,vidx,ftensor, niter=niter, atol=atol)
    return mP_bp, isconverged
end


text:
# Tests for Factor Analysis
function t_fa()

N = 2	# variables
M = 250 # samples

matrixdata = randn(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(matrixdata)
A02 = j4pr.datacell(matrixdata, labels)
A03 = j4pr.datacell(matrixdata, multilabels)

for tr in [A01, A02, A03, matrixdata]
	for ts in [A01, A02, A03, matrixdata]
		Test.@test try
			Wfa = tr |> j4pr.fa()
			Wfar = j4pr.far(Wfa)
			ts |> Wfa
			(ts |> Wfa) |> Wfar
			true
		catch
			false
		end
	end
end

end


text:
using Test

@testset "tests" begin
  include(joinpath("advection", "runtests.jl"))
  include(joinpath("euler", "runtests.jl"))
end


text:
export @alchemy

macro alchemy(block)
    mol = repr(block)
    mol = replace(mol, ";"=>"\n")
    mol = strip(filter(c->!occursin(c,"{}():"), mol))
    mol = String(mol)
    mol = Molecules.parse_string(mol)
    quote
        plot($mol)
    end
end

macro alchemy(mol::Symbol)
    return :(Alchemy.plot($mol)) |> esc
end

macro alchemy(fpath::String)
    mol = Molecules.parse_file(fpath)
    return :(Alchemy.plot($mol)) |> esc
end

function plot(atoms::Vector{T}) where T <: Molecules.Atom
    S = Scene()
    plot!(S, atoms)
    return S
end

function plot!(S, atoms::Vector{T}) where T <: Molecules.Atom
    for i = eachindex(atoms)
        center = Point3f0(atoms[i].xyz)
        symbol = Molecules.symbol(atoms[i])
        r = get_atom_radius(symbol)
        c = get_atom_color(symbol)

        sphere = Sphere(center, r)

        mesh!(S, sphere, show_axis=false, color=c)
    end

    # Add chemical bonds for close atoms
    for i in eachindex(atoms)
        A1 = atoms[i]
        r1 = A1.xyz
        for j in (i+1):length(atoms)
            A2 = atoms[j]
            r2 = A2.xyz
            d = √sum((r1 .- r2).^2)
            if d < 1.55
                xyz = [Point3f0(A1.xyz)]
                uvw = [Point3f0(A2.xyz - A1.xyz)]
                arrows!(S, xyz, uvw, linewidth=0.1, arrowsize=0.0, linecolor=stegeman)
            end
        end
    end

    return S
end

max_stars_repo_name: gustavojra/Alchemy

max_issues_repo_name: gustavojra/Alchemy

max_forks_repo_name: gustavojra/Alchemy


text:
# runtest.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2009 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
using Ollam
using MNIST, Stage, LIBSVM, GZip
import Base: length, start, done, next
using Base.Test

# ----------------------------------------------------------------------------------------------------------------
# Helpers
# ----------------------------------------------------------------------------------------------------------------
immutable EachCol{T}
  matrix :: Array{T, 2}
end
length(e :: EachCol)      = size(e.matrix, 2)
start(e :: EachCol)       = 1
next(e :: EachCol, state) = (e.matrix[:, state], state + 1)
done(e :: EachCol, state) = state > length(e) ? true : false

type T
  name     :: AbstractString
  trainer  :: Function
  expected :: Float64
  testset  :: Matrix
end

# ----------------------------------------------------------------------------------------------------------------
# test hildreth
# ----------------------------------------------------------------------------------------------------------------
h = setup_hildreth(C = Inf, k = 2)
@expect abs(hildreth((SparseMatrixCSC)[ spzeros(1235, 1),
                sparsevec(Dict(470=>-1.0, 1070=>-1.0, 1231=>-1.0, 1232=>-1.0, 1233=>-1.0, 1234=>-1.0, 1235=>-1.0, 496=>1.0,
                           1058=>1.0, 1226=>1.0, 1227=>1.0, 1228=>1.0, 1229=>1.0, 1230=>1.0)) ], [ 0.0, 0.9733606705258293 ], h)[2] - 0.069526) < 0.00001

# ----------------------------------------------------------------------------------------------------------------
# test lazy maps
# ----------------------------------------------------------------------------------------------------------------
xxx = [1, 2, 3]
for i in lazy_map(f -> f+1, xxx)
  @expect i == xxx[i-1] + 1
end

# ----------------------------------------------------------------------------------------------------------------
# Setup buzz task
# ----------------------------------------------------------------------------------------------------------------
@info "reading Auto-MPG train/test data"

f = gzopen("test/auto-mpg.data.gz")
const mpg_data = readdlm(f)'
close(f)

srand(0)
const rand         = shuffle(collect(1:size(mpg_data, 2)))
const rall         = map(Float64,mpg_data[2:end-1, :])[:, rand]
const rall_truth   = vec(map(Float64,mpg_data[1, :]))[rand]
const norm_truth   = (rall_truth - minimum(rall_truth)) / (maximum(rall_truth) - minimum(rall_truth))
const m            = zeros(size(rall))
const s            = ones(size(rall))
const norm_all     = [ (rall[i, j] - m[i]) / s[i] for i = 1:size(rall, 1), j = 1:size(rall, 2) ]
const rtrain       = norm_all[:, 1:round(Int,size(rall, 2)*0.8)]
const rtrain_truth = norm_truth[1:round(Int,size(rall, 2)*0.8)]
const rtest        = norm_all[:, round(Int,size(rall, 2)*0.8)+1:end]
const rtest_truth  = norm_truth[round(Int,size(rall, 2)*0.8)+1:end]
@debug "size of train: $(size(rtrain)), test: $(size(rtest))"

# ----------------------------------------------------------------------------------------------------------------
# Regression tests
# ----------------------------------------------------------------------------------------------------------------
regression_tests = [
  T("perceptron",          (init) -> regress_perceptron(EachCol(rtrain), rtrain_truth, init; iterations = 400, average = false, learn_rate = 1e-7, C = 0.1), 0.106, rtest),
  T("averaged perceptron", (init) -> regress_perceptron(EachCol(rtrain), rtrain_truth, init; iterations = 400, average = true, learn_rate = 1e-7, C = 0.1),  0.114, rtest),
  T("MIRA",                (init) -> regress_mira(EachCol(rtrain), rtrain_truth, init; iterations = 400, average = false),                                   0.099, rtest),
  T("averaged MIRA",       (init) -> regress_mira(EachCol(rtrain), rtrain_truth, init; iterations = 400, average = true),                                    0.104, rtest),
]

# Ollam tests
for t in regression_tests
  init  = RegressionModel(length(rtrain[:, 1]))
  @timer "training $(t.name) model" model = t.trainer(init)
  res = sqrt(test_regression(model, EachCol(t.testset), rtest_truth))
  @info @sprintf("%s test set error rate: %7.3f", t.name, res)
  @expect abs(res - t.expected) < 0.001
end

# ----------------------------------------------------------------------------------------------------------------
# setup MNIST task
# ----------------------------------------------------------------------------------------------------------------
@info "reading MNIST train/test data"
const train_raw, train_truth    = traindata()
const test_raw, test_truth      = testdata()
const classes                   = Dict{Float64, Int32}()
const train                     = train_raw / 255.0
const test                      = test_raw / 255.0
# const train_deskew, test_deskew = begin
#   f   = gzopen("mnist-deslant.jd.gz")
#   trd = deserialize(f)
#   trt = deserialize(f)
#   tsd = deserialize(f)
#   tst = deserialize(f)
#   trd / 255.0, tsd / 255.0
# end

# ----------------------------------------------------------------------------------------------------------------
# Testing support
# ----------------------------------------------------------------------------------------------------------------
i = 1
for t in test_truth
  if !(t in keys(classes))
    classes[t] = i
    i += 1
  end
end
@info "truth set: $classes"

# ----------------------------------------------------------------------------------------------------------------
# Classification tests
# ----------------------------------------------------------------------------------------------------------------
classification_tests = [
  # deskew tests
  # T("5-best MIRA [deskew]",          (init) -> train_mira(EachCol(train_deskew), train_truth, init; k = 5, iterations = 30, average = false), 6.73, test_deskew),
  # T("averaged 5-best MIRA [deskew]", (init) -> train_mira(EachCol(train_deskew), train_truth, init; k = 5, iterations = 30, average = true),  4.54, test_deskew),
  # T("linear SVM [deskew]",           (init) -> train_svm(EachCol(train_deskew), train_truth; C = 0.001, iterations = 100),                    5.40, test_deskew),
  # T("linear libSVM",  (init) -> train_libsvm(EachCol(train_deskew), train_truth, C = 1.0, cache_size = 250.0, eps = 0.001, shrinking = true), 7.83, test_deskew),

  # baselines
  #T("perceptron",           (init) -> train_perceptron(EachCol(train), train_truth, init; iterations = 50, average = false),  12.74, test),
  T("averaged perceptron",  (init) -> train_perceptron(EachCol(train), train_truth, init; iterations = 50, average = true),   12.69, test),
  #T("1-best MIRA",          (init) -> train_mira(EachCol(train), train_truth, init; k = 1, iterations = 30, average = false), 14.26, test),
  T("averaged 1-best MIRA", (init) -> train_mira(EachCol(train), train_truth, init; k = 1, iterations = 30, average = true),   8.23, test),
  #T("5-best MIRA",          (init) -> train_mira(EachCol(train), train_truth, init; k = 5, iterations = 30, average = false), 11.46, test),
  T("averaged 5-best MIRA", (init) -> train_mira(EachCol(train), train_truth, init; k = 5, iterations = 30, average = true),   7.72, test),
  T("linear SVM",           (init) -> train_svm(EachCol(train), train_truth; C = 0.001, iterations = 100),                     8.73, test),
  # T("linear libSVM",        (init) -> train_libsvm(EachCol(train[:, 1:20000]), train_truth[1:20000],
  #                                                     C = 1.0, cache_size = 250.0, eps = 0.001, shrinking = true),             8.85, test),
]

# Ollam tests
for t in classification_tests
  init  = LinearModel(classes, length(train[:, 1]))
  @timer "training $(t.name) model" model = t.trainer(init)
  res = test_classification(model, EachCol(t.testset), test_truth) * 100.0
  @info @sprintf("%s test set error rate: %7.3f%%", t.name, res)
  @expect abs(res - t.expected) < 0.001
end

# baseline test with degree-3 RBF kernel
# @timer "libsvm 2nd-order rbf direct" model = svmtrain(train_truth, train, C = 10.0, verbose = true, shrinking = true)
# (predicted_labels, decision_values) = svmpredict(model, test)
# @info "test libsvm 2nd-order rbf direct: $((1.0 - mean(predicted_labels .== test_truth))*100.0)"
# @expect abs(((1.0 - mean(predicted_labels .== test_truth))*100.0) - 3.87) < 0.001


text:
##############################################################
#oxdlcccloloollllcccccclollcclodooodxkkxdddxxddollc:::::;;;;;#
#dkxdoooodddoolllcllloollllllodddddddxkOOkxxxxxxxxoccc:::c:::#
#xxdddxkkxdloollloddddoollol:;,;;lkOOkxxdddddxdxkkdollcccc:::#
#dxkxxxxkxxdooodxxxool,;,:;;'.',;;,cl;;coxxxxxkkkkkkxxdolcccc#
#kOkxkxxkkOkxxxkOOdl;.,,',',:;,;,,;;,::;;::oOkOOOOOkxdooollcc#
#OkkkkkO0KOOO0KKK0c:..'.',',,..,;,,;;;;,;;:l0KKXKKKK0OOkxdool#
#KKK00OkOkkO0KKKKdc;',',;::,',.;:,'..',;,'',;x0KXKKKK00OO0Okx#
#OOO000OOO0OO0kkd::;,,;,'';,.''.....',,,,..',:OOxk00K0K00K00O#
#K0000KKKKK00dc:;;,:.....';;;cc.....,;'.,,'',:;cc,dXXXXXKKK0O#
#XXXXXXXKKKKKd;',;;:..'',...',. .......,,'.',,;;,,lNNXXXXXKK0#
#XXNNNNNNXXkl..,,'',....,,'.',,.'.......;::,:;';,':x0XXKXXXXK#
#00000KKXXXd;'';,'.,;..''.;..';'',,''..,'',';;;'',,;l00XXXKK0#
#KKKKXXXKKKKo':od,.ll .,','. .''..':o:........',,,,cckKNNXXXX#
#KKKKKKKKKko::lc;,'.'.....',... .,;.'''.';'..''..,;coKXNNNNXX#
#XXXXXXXXkc,;,,cc:;,,:''......;..,;:,.,;l:c,;;.;c;lxXXNNNNNNN#
#XXXXXXXXX0cdcl:co,;;':ll,''..:.,'';:,.cxOxl,:c:0NWWWNXXXXXXN#
#NNNNNNNNNNNNXxkNNKOOXNNNX0K0kk.,o:clcxdONNXKoKNWNNNNNNNXXNNN#
#WWWWWWWWWWWWWWWWWWWWWWWWWWWWWK..NNNNNNWWWWNNNNNNNWNNNNNNNNNN#
#NNNNXXXXXXXKKKKKKK0kdxkxkxxdxl::dxdxxxxkOkkO0KXXXXXXKKKKKKKK#
#dddxxxxxxddddddxddddddddddddddddddooodoooolllllollllllllllll#
##############################################################

# BKTrees.jl - Burkhard-Keller trees implementation written
#              at 0x0α Research by Corneliu Cofaru, 2018

module BKTrees

using DataStructures: Deque
import Base: show

if VERSION < v"1.0"
    import Base: find
end

export BKTree, Node, add!, find

include("tree.jl")

end # module


text:
module CategoricalArrays
    export CategoricalPool, CategoricalValue
    export AbstractCategoricalArray, CategoricalArray,
           CategoricalVector, CategoricalMatrix
    export AbstractNullableCategoricalArray, NullableCategoricalArray,
           NullableCategoricalVector, NullableCategoricalMatrix

    export categorical, compact, droplevels!, levels, levels!, ordered, ordered!

    using Compat

    include("typedefs.jl")

    include("buildfields.jl")

    include("pool.jl")
    include("value.jl")

    include("array.jl")
    include("nullablearray.jl")

    if VERSION < v"0.5.0-dev"
        Base.convert{T,n,S}(::Type{Array{T}}, x::AbstractArray{S, n}) = convert(Array{T, n}, x)
        Base.convert{T,n,S}(::Type{Array{T,n}}, x::AbstractArray{S,n}) = copy!(Array{T}(size(x)), x)

        Base.convert{T}(::Type{Nullable   }, x::T) = Nullable{T}(x)

        Base.promote_op{S<:CategoricalValue, T<:CategoricalValue}(::typeof(@functorize(==)),
                                                                  ::Type{S}, ::Type{T}) = Bool
        Base.promote_op{S<:CategoricalValue, T<:CategoricalValue}(::typeof(@functorize(>)),
                                                                  ::Type{S}, ::Type{T}) = Bool
        Base.promote_op{S<:CategoricalValue, T<:CategoricalValue}(::typeof(@functorize(<)),
                                                                  ::Type{S}, ::Type{T}) = Bool
        Base.promote_op{S<:CategoricalValue, T<:CategoricalValue}(::typeof(@functorize(>=)),
                                                                  ::Type{S}, ::Type{T}) = Bool
        Base.promote_op{S<:CategoricalValue, T<:CategoricalValue}(::typeof(@functorize(<=)),
                                                                  ::Type{S}, ::Type{T}) = Bool
    end
end


text:
module SpaceJam

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Space Jam" => Ahorn.EntityPlacement(
        Maple.DreamBlock,
        "rectangle"
    )
)

Ahorn.nodeLimits(entity::Maple.DreamBlock) = 0, 1

Ahorn.minimumSize(entity::Maple.DreamBlock) = 8, 8
Ahorn.resizable(entity::Maple.DreamBlock) = true, true

function Ahorn.selection(entity::Maple.DreamBlock)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    nodes = get(entity.data, "nodes", ())
    if isempty(nodes)
        return Ahorn.Rectangle(x, y, width, height)

    else
        nx, ny = Int.(nodes[1])
        return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
    end
end

function renderSpaceJam(ctx::Ahorn.Cairo.CairoContext, x::Int, y::Int, width::Int, height::Int)
    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    Ahorn.drawRectangle(ctx, x, y, width, height, (0.0, 0.0, 0.0, 0.4), (1.0, 1.0, 1.0, 1.0))

    Ahorn.restore(ctx)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DreamBlock)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        renderSpaceJam(ctx, nx, ny, width, height)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DreamBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderSpaceJam(ctx, 0, 0, width, height)
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# auto-generated tests from julia-repl docstrings
using Test, LaurentPolynomials
function mytest(f::String,a::String,b::String)
  println(f," ",a)
  omit=a[end]==';'
  a=replace(a,"\\\\"=>"\\")
  a=repr(MIME("text/plain"),eval(Meta.parse(a)),context=:limit=>true)
  if omit a="nothing" end
  a=replace(a,r" *(\n|$)"s=>s"\1")
  a=replace(a,r"\n$"s=>"")
  b=replace(b,r" *(\n|$)"s=>s"\1")
  b=replace(b,r"\n$"s=>"")
  i=1
  while i<=lastindex(a) && i<=lastindex(b) && a[i]==b[i]
    i=nextind(a,i)
  end
  if a!=b print("exec=$(repr(a[i:end]))\nmanl=$(repr(b[i:end]))\n") end
  a==b
end
@testset verbose = true "Gapjm" begin
@testset "README.md" begin
@test mytest("README.md","Pol(:q)","Pol{Int64}: q")
@test mytest("README.md","@Pol q","Pol{Int64}: q")
@test mytest("README.md","Pol([1,2])","Pol{Int64}: 2q+1")
@test mytest("README.md","2q+1","Pol{Int64}: 2q+1")
@test mytest("README.md","Pol()","Pol{Int64}: q")
@test mytest("README.md","p=Pol([1,2,1],-1)","Pol{Int64}: q+2+q⁻¹")
@test mytest("README.md","q+2+q^-1","Pol{Int64}: q+2+q⁻¹")
@test mytest("README.md","valuation(p),degree(p)","(-1, 1)")
@test mytest("README.md","p[0], p[1], p[-1], p[10]","(2, 1, 1, 0)")
@test mytest("README.md","p[valuation(p):degree(p)]","3-element Vector{Int64}:\n 1\n 2\n 1")
@test mytest("README.md","p[begin:end]","3-element Vector{Int64}:\n 1\n 2\n 1")
@test mytest("README.md","coefficients(p)","3-element Vector{Int64}:\n 1\n 2\n 1")
@test mytest("README.md","Pol(1)","Pol{Int64}: 1")
@test mytest("README.md","convert(Pol{Int},1)","Pol{Int64}: 1")
@test mytest("README.md","scalar(Pol(1))","1")
@test mytest("README.md","convert(Int,Pol(1))","1")
@test mytest("README.md","Int(Pol(1))","1")
@test mytest("README.md","scalar(q+1)","nothing")
@test mytest("README.md","derivative(p)","Pol{Int64}: 1-q⁻²")
@test mytest("README.md","p=(q+1)^2","Pol{Int64}: q²+2q+1")
@test mytest("README.md","p=(q+1)^2","Pol{Int64}: q²+2q+1")
@test mytest("README.md","p/2","Pol{Float64}: 0.5q²+1.0q+0.5")
@test mytest("README.md","p//2","Pol{Rational{Int64}}: (1//2)q²+(1//1)q+1//2")
@test mytest("README.md","p(1//2)","9//4")
@test mytest("README.md","p(0.5)","2.25")
@test mytest("README.md","Pol([1,2,3],[2.0,1.0,3.0])","Pol{Float64}: 1.5q²-5.5q+6.0")
@test mytest("README.md","divrem(q^3+1,2q+1)","(0.5q²-0.25q+0.125, 0.875)")
@test mytest("README.md","divrem(q^3+1,2q+1//1)","((1//2)q²+(-1//4)q+1//8, 7//8)")
@test mytest("README.md","pseudodiv(q^3+1,2q+1)","(4q²-2q+1, 7)")
@test mytest("README.md","(4q^2-2q+1)*(2q+1)+7","Pol{Int64}: 8q³+8")
@test mytest("README.md","exactdiv(q+1,2.0)","Pol{Float64}: 0.5q+0.5")
@test mytest("README.md","a=1/(q+1)","Frac{Pol{Int64}}: 1/(q+1)")
@test mytest("README.md","Pol(2/a)","Pol{Int64}: 2q+2")
@test mytest("README.md","numerator(a)","Pol{Int64}: 1")
@test mytest("README.md","denominator(a)","Pol{Int64}: q+1")
@test mytest("README.md","m=[q+1 q+2;q-2 q-3]","2×2 Matrix{Pol{Int64}}:\n q+1  q+2\n q-2  q-3")
@test mytest("README.md","n=inv(Frac.(m))","2×2 Matrix{Frac{Pol{Int64}}}:\n (-q+3)/(2q-1)  (-q-2)/(-2q+1)\n (q-2)/(2q-1)   (q+1)/(-2q+1)")
@test mytest("README.md","map(x->x(1),n)","2×2 Matrix{Float64}:\n  2.0   3.0\n -1.0  -2.0")
@test mytest("README.md","map(x->x(1;Rational=true),n)","2×2 Matrix{Rational{Int64}}:\n  2//1   3//1\n -1//1  -2//1")
end
@testset "LaurentPolynomials.jl" begin
@test mytest("LaurentPolynomials.jl","gcd(2q+2,q^2-1)","Pol{Int64}: q+1")
@test mytest("LaurentPolynomials.jl","gcd(q+1//1,q^2-1//1)","Pol{Rational{Int64}}: (1//1)q+1//1")
@test mytest("LaurentPolynomials.jl","gcdx(q^3-1//1,q^2-1//1)","((1//1)q-1//1, 1//1, (-1//1)q)")
@test mytest("LaurentPolynomials.jl","powermod(q-1,3,q^2+q+1)","Pol{Int64}: 6q+3")
@test mytest("LaurentPolynomials.jl","p=Pol([1,1,1])","Pol{Int64}: q²+q+1")
@test mytest("LaurentPolynomials.jl","vals=p.(1:5)","5-element Vector{Int64}:\n  3\n  7\n 13\n 21\n 31")
@test mytest("LaurentPolynomials.jl","Pol(1:5,vals*1//1)","Pol{Rational{Int64}}: (1//1)q²+(1//1)q+1//1")
@test mytest("LaurentPolynomials.jl","Pol(1:5,vals*1.0)","Pol{Float64}: 1.0q²+1.0q+1.0")
end
end


text:
function resolve_import(x, state::State)
    if typof(x) === Using || typof(x) === Import
        u = typof(x) === Using
        i = 2
        n = length(x)

        root = par = getsymbolserver(state.server)
        while i <= n
            arg = x[i]
            if is_id_or_macroname(arg)
                if refof(arg) !== nothing
                    par = refof(arg)
                else
                    par = _get_field(par, arg, state)
                end
            elseif ispunctuation(arg) && kindof(arg) == CSTParser.Tokens.COMMA
                # end of chain, make available
                if i > 2
                    _mark_import_arg(x[i - 1], par, state, u)
                end
                par = root
            elseif isoperator(arg) && kindof(arg) == CSTParser.Tokens.COLON
                root = par
                if par !== nothing && i > 2 && isidentifier(x[i - 1]) && refof(x[i - 1]) === nothing
                    setref!(x[i - 1], par)
                end
            elseif ispunctuation(arg) && kindof(arg) == CSTParser.Tokens.DOT
                # dot between identifiers
                if par !== nothing && i > 2 && isidentifier(x[i - 1]) && refof(x[i - 1]) === nothing
                    setref!(x[i - 1], par)
                end
            elseif isoperator(arg) && kindof(arg) == CSTParser.Tokens.DOT
                # dot prexceding identifiser
                if par == root == getsymbolserver(state.server)
                    par = state.scope
                elseif par isa Scope && parentof(par) !== nothing
                    par = parentof(par)
                else
                    return
                end
            elseif !isoperator(arg)
                return
            end
            if i == n
                _mark_import_arg(x[i], par, state, u)
            end
            i += 1
        end
    end
end

function _mark_import_arg(arg, par, state, usinged)
    if par !== nothing && is_id_or_macroname(arg)
        if par isa Binding # mark reference to binding
            push!(par.refs, arg)
        end
        if par isa SymbolServer.VarRef
            par = SymbolServer._lookup(par, getsymbolserver(state.server), true)
            !(par isa SymbolServer.SymStore) && return
        end
        if bindingof(arg) === nothing
            if !hasmeta(arg)
                arg.meta = Meta()
            end
            arg.meta.binding = Binding(arg, par, _typeof(par, state), [], nothing, nothing)
        end

        if usinged
            if par isa SymbolServer.ModuleStore
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), par)
            elseif par isa Binding && par.val isa SymbolServer.ModuleStore
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), par.val)
            elseif par isa Binding && par.val isa EXPR && CSTParser.defines_module(par.val)
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), scopeof(par.val))
            elseif par isa Binding && par.val isa Binding && par.val.val isa EXPR && CSTParser.defines_module(par.val.val)
                add_to_imported_modules(state.scope, Symbol(valofid(arg)), scopeof(par.val.val))
            end
        end
    end
end

function add_to_imported_modules(scope::Scope, name::Symbol, val)
    if scope.modules isa Dict
        scope.modules[name] = val
    else
        modules = Dict(name => val)
    end
end
no_modules_above(s::Scope) = !CSTParser.defines_module(s.expr) || s.parent === nothing || no_modules_above(s.parent)
function get_named_toplevel_module(s::Scope, name::String)
    if CSTParser.defines_module(s.expr) && valofid(CSTParser.get_name(s.expr)) == name && no_modules_above(s)
        return s.expr
    elseif s.parent isa Scope
        return get_named_toplevel_module(s.parent, name)
    end
    return nothing
end
function _get_field(par, arg, state)
    arg_str_rep = CSTParser.str_value(arg)
    if par isa SymbolServer.EnvStore
        if (arg_scope = retrieve_scope(arg)) !== nothing && (tlm = get_named_toplevel_module(arg_scope, arg_str_rep)) !== nothing && hasbinding(tlm)
            return bindingof(tlm)
        elseif haskey(par, Symbol(arg_str_rep))
            return par[Symbol(arg_str_rep)]
        end
    elseif par isa SymbolServer.ModuleStore # imported module
        if Symbol(arg_str_rep) === par.name.name
            return par
        elseif haskey(par, Symbol(arg_str_rep))
            par = par[Symbol(arg_str_rep)]
            if par isa SymbolServer.VarRef # reference to dependency
                return SymbolServer._lookup(par, getsymbolserver(state.server), true)
            end
            return par
        end
        for used_module_name in par.used_modules
            used_module = maybe_lookup(par[used_module_name], state.server)
            if used_module !== nothing && isexportedby(Symbol(arg_str_rep), used_module)
                return used_module[Symbol(arg_str_rep)]
            end
        end
    elseif par isa Scope
        if scopehasbinding(par, arg_str_rep)
            return par.names[arg_str_rep]
        elseif par.modules !== nothing
            for used_module in values(par.modules)
                if used_module isa SymbolServer.ModuleStore && isexportedby(Symbol(arg_str_rep), used_module)
                    return maybe_lookup(used_module[Symbol(arg_str_rep)], state.server)
                elseif used_module isa Scope && scope_exports(used_module, arg_str_rep, state)
                    return used_module.names[arg_str_rep]
                end
            end
        end
    elseif par isa Binding
        if par.val isa Binding
            return _get_field(par.val, arg, state)
        elseif par.val isa EXPR && CSTParser.defines_module(par.val) && scopeof(par.val) isa Scope
            return _get_field(scopeof(par.val), arg, state)
        elseif par.val isa SymbolServer.ModuleStore
            return _get_field(par.val, arg, state)
        end
    end
    return
end


text:
module SigmaPoints
export MerweScaled, calculate_σ_points, num_σ_points

using LinearAlgebra: cholesky, det, Symmetric

"""
```julia
MerweScaled(
    ;n::Int64, α::Float64, β::Float64, κ::Float64, residual_x::Function = -
)
```

# Parameters for the Van der Merwe's σ-points
- `n::Int64`: Dimensionality of the state. `2n+1` will be generated
- `α::Float64`: Spread of the σ-points around the mean.
- `β::Float64`: Prior knowledge of the distribution of the mean.
- `κ::Float64`: Second scaling parameter.
- `Σ_w::Array{Float64, 1}`: Weights for each σ-point for the covariance.
- `m_w::Array{Float64, 1}`: Weights for each σ-point for the mean.
- `residual_x::Function = -`: Function to compute residual between states.
"""
struct MerweScaled
    n::Int64
    α::Float64
    β::Float64
    κ::Float64

    Σ_w::Array{Float64, 1}
    m_w::Array{Float64, 1}

    residual_x::Function
end

function MerweScaled(
    ;n::Int64, α::Float64, β::Float64, κ::Float64, residual_x::Function = -
)
    Σ_w, m_w = σ_weights(n=n, α=α, β=β, κ=κ)
    MerweScaled(n, α, β, κ, Σ_w, m_w, residual_x)
end

function Base.show(io::IO, merwe::MerweScaled)
    print(
        io, "Merwe Scaled:\n",
        "n: ", repr(MIME("text/plain"), merwe.n, context=io), "\n",
        "α: ", repr(MIME("text/plain"), merwe.α, context=io), "\n",
        "β: ", repr(MIME("text/plain"), merwe.β, context=io), "\n",
        "κ: ", repr(MIME("text/plain"), merwe.κ, context=io), "\n",
        "Σ_w: ", repr(MIME("text/plain"), merwe.Σ_w, context=io), "\n",
        "m_w: ", repr(MIME("text/plain"), merwe.m_w, context=io), "\n",
        "residual_x: ", repr(MIME("text/plain"), merwe.residual_x, context=io), "\n",
    )
end

"""
```julia
num_σ_points(σ_parameters::MerweScaled)::Int64
```

Return number of σ-points for given σ-parameters.
"""
num_σ_points(σ_parameters::MerweScaled)::Int64 = 2σ_parameters.n + 1

"""
```julia
σ_weights(
    ;n::Int64, α::Float64, β::Float64, κ::Float64
)::Tuple{Array{Float64, 1}, Array{Float64, 1}}
```
Compute weights for each σ-point for the covariance `P` and mean `x`.

# Arguments
- `n::Int64`: Dimensionality of the state. `2n+1` σ-points will be generated
- `α::Float64`: Spread of the σ-points around the mean.
- `β::Float64`: Prior knowledge of the distribution of the mean.
- `κ::Float64`: Second scaling parameter.
"""
function σ_weights(
    ;n::Int64, α::Float64, β::Float64, κ::Float64
)::Tuple{Array{Float64, 1}, Array{Float64, 1}}
    λ = α ^ 2 * (n + κ) - n
    c = 0.5 / (n + λ)

    Σ_w = [c for i = 1:2n + 1]
    m_w = copy(Σ_w)

    w = λ / (n + λ)
    Σ_w[1] = w + (1 - α ^ 2 + β)
    m_w[1] = w
    Σ_w, m_w
end

"""
```julia
calculate_σ_points(
    ;σ_parameters::MerweScaled, x::Array{Float64, 1}, P::Array{Float64, 2}
)::Array{Float64, 2}
```

Compute `2n+1 x n` Van der Merwe's σ-points
for an Unscented Kalman Filter (UKF) given the mean `x` and covariance `P`.

# Arguments
- `σ_parameters::MerweScaled`: Parameters of the Van der Merwe's σ-points.
- `x::Array{Float64, 1}`: Mean of the UKF.
- `P::Array{Float64, 2}`: Covariance of the `x` of the UKF.
"""
function calculate_σ_points(
    ;σ_parameters::MerweScaled, x::Array{Float64, 1}, P::Array{Float64, 2}
)::Array{Float64, 2}
    λ = σ_parameters.α ^ 2 * (σ_parameters.n + σ_parameters.κ) - σ_parameters.n
    U = cholesky(Symmetric(P * (λ + σ_parameters.n))).U

    Σ = zeros(Float64, 2σ_parameters.n + 1, σ_parameters.n)
    @inbounds begin
    Σ[1, :] = x
    for i = 1:σ_parameters.n
        Σ[i + 1, :] = σ_parameters.residual_x(x, -@view(U[i, :]))
        Σ[σ_parameters.n + i + 1, :] = σ_parameters.residual_x(x, @view(U[i, :]))
    end
    end
    Σ
end

end

max_stars_repo_name: pxl-th/KalmanFilters

max_issues_repo_name: pxl-th/KalmanFilters

max_forks_repo_name: pxl-th/KalmanFilters


text:
using BlockDiagonalMatrices
using BlockDiagonalMatrices: svd_blockwise
using LinearAlgebra
using Random
using Test

@testset "linalg.jl" begin
    rng = MersenneTwister(123456)
    N1, N2, N3 = 3, 4, 5
    N = N1 + N2 + N3
    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
    A = rand(rng, N, N + N1)
    B = rand(rng, N + N1, N + N2)
    A′, B′ = A', B'
    a = rand(rng, N)
    b = rand(rng, N + N1)

    @testset "mul!" begin
        c = similar(b1)
        d = similar(Matrix(b1))
        mul!(c, b1, b1)
        mul!(d, Matrix(b1), Matrix(b1))
        @test c ≈ d
        if VERSION ≥ v"1.3"
            mul!(c, b1, b1, 2.0, 3.0)
            mul!(d, Matrix(b1), Matrix(b1), 2.0, 3.0)
            @test c ≈ d
        end
    end

    @testset "Unary Linear Algebra" begin
        @testset "$f" for f in (adjoint, det, diag, eigvals, inv, pinv, svdvals, transpose, tr)
            @test f(b1) ≈ f(Matrix(b1))
        end

        @testset "permute=$p, scale=$s" for p in (true, false), s in (true, false)
            @testset "$g" for g in (eigmin, eigmax)
                # `b2` has real eigenvals, required for `eigmin`, `b1` has Complex eigenvals
                @test g(b2; permute=p, scale=s) ≈ g(Matrix(b2); permute=p, scale=s)
            end
            @testset "eigvals" begin
                result = eigvals(b2; permute=p, scale=s)
                expected = eigvals(Matrix(b2), permute=p, scale=s)
                @static if VERSION < v"1.2"
                    # the `eigvals` method we hit above did not sort real eigenvalues pre-v1.2
                    # but some `eigvals` methods did, so we always sort real eigenvalues.
                    @test sort!(result) ≈ sort!(expected)
                else
                    @test result ≈ expected
                end
            end
        end

        @testset "eigvals on LinearAlgebra types" begin
            # `eigvals` has different methods for different types, e.g. Hermitian
            b_herm = BlockDiagonal([Hermitian(rand(rng, 3, 3) + I) for _ in 1:3])
            @test eigvals(b_herm) ≈ eigvals(Matrix(b_herm))
            @test eigvals(b_herm, 1.0, 2.0) ≈ eigvals(Hermitian(Matrix(b_herm)), 1.0, 2.0)
        end

        # Requires a postive semi-definite matrix
        @testset "logdet" begin
            b̂1 = b1 * b1'
            @test logdet(b̂1) ≈ logdet(Matrix(b̂1))
        end
    end  # Unary

    @testset "Cholesky decomposition" begin
        X = [  4  12 -16
              12  37 -43
             -16 -43  98]
        U = [ 2.0 6.0 -8.0
              0.0 1.0  5.0
              0.0 0.0  3.0]

        B = BlockDiagonal([X, X])
        C = cholesky(B)
        @test C isa Cholesky{Float64, <:BlockDiagonal{Float64}}
        @test C.U ≈ cholesky(Matrix(B)).U
        @test C.U ≈ BlockDiagonal([U, U])
        @test C.L ≈ BlockDiagonal([U', U'])
        @test C.UL ≈ C.U
        @test C.uplo === 'U'
        @test C.info == 0

        M = BlockDiagonal(map(Matrix, blocks(C.L)))
        C = Cholesky(M, 'L', 0)
        @test C.U ≈ cholesky(Matrix(B)).U
        @test C.U ≈ BlockDiagonal([U, U])
        @test C.L ≈ BlockDiagonal([U', U'])
        @test C.UL ≈ C.L
        @test C.uplo === 'L'
        @test C.info == 0
    end  # Cholesky

    @testset "Singular Value Decomposition" begin
        X = [  4  12 -16
              12  37 -43
             -16 -43  98]
        B = BlockDiagonal([X, X])

        @testset "full=$full" for full in (true, false)

            @testset "svd_blockwise" begin
                U, S, Vt = svd_blockwise(B; full=full)
                F = SVD(U, S, Vt)
                @test B ≈ F.U * Diagonal(F.S) * F.Vt

                # Matrices should be BlockDiagonal
                @test F isa SVD{Float64, Float64, <:BlockDiagonal{Float64}}
                @test F.U isa BlockDiagonal
                @test F.V isa BlockDiagonal
                @test F.Vt isa BlockDiagonal

                # Should have same values, but not sorted so as to keep BlockDiagonal structure
                F_ = svd(Matrix(B), full=full)
                for fname in fieldnames(SVD)
                    @test sort(vec(getfield(F, fname))) ≈ sort(vec(getfield(F_, fname)))
                end
                # Singular values should be block-wise
                s = svd(X).S
                @test F.S == vcat(s, s)
            end

            @testset "svd" begin
                F = svd(B; full=full)
                F_ = svd(Matrix(B), full=full)

                @test F isa SVD
                @test B ≈ F.U * Diagonal(F.S) * F.Vt

                @test F == F_
                for fname in fieldnames(SVD)
                    @test getfield(F, fname) ≈ getfield(F_, fname)
                end

                # Singular values should be sorted in descending order
                @test F.S == sort(F.S, rev=true)
            end
        end
    end  # SVD
end


text:
localhost = filter(x -> x[1] == "localhost", worker_list)
remote_servers = filter(x -> x[1] != "localhost", worker_list)

length(localhost) > 0 && addprocs(localhost[1][2], exeflags=exeflags)
length(remote_servers) > 0 && addprocs(remote_servers, sshflags=sshflags, exeflags=exeflags)

# validate package versions
@everywhere function get_git_hash(path)
    cmd = `git -C $path rev-parse HEAD`
    (gethostname(), strip(read(cmd, String)))
end

function setup_julia_environment()
    local_githash = get_git_hash(JULIA_ENV)[2]
    for id in workers()
        remote_name, remote_githash = remotecall_fetch(get_git_hash, id, JULIA_ENV)
        @assert remote_githash == local_githash "Host: $remote_name has version mismatch." *
                                                    "Hash is '$remote_githash' instead of '$local_githash'."
    end
end
@everywhere using Pkg
setup_julia_environment()

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
using Coord
using Base.Test

include("includes.jl")
include("test_dynamics.jl")
include("test_encounter.jl")
include("test_environment.jl")
include("test_state.jl")
include("test_strategies.jl")
include("test_tables.jl")
include("test_uav.jl")

println("All tests pass!")


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

module ComputeManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_AvailabilitySetsApi.jl")
include("api_ImagesApi.jl")
include("api_SkusApi.jl")
include("api_UsageApi.jl")
include("api_VirtualMachineExtensionImagesApi.jl")
include("api_VirtualMachineExtensionsApi.jl")
include("api_VirtualMachineImagesApi.jl")
include("api_VirtualMachineScaleSetExtensionsApi.jl")
include("api_VirtualMachineScaleSetRollingUpgradesApi.jl")
include("api_VirtualMachineScaleSetVMsApi.jl")
include("api_VirtualMachineScaleSetsApi.jl")
include("api_VirtualMachineSizesApi.jl")
include("api_VirtualMachinesApi.jl")

# export models
export convert, AdditionalUnattendContent
export convert, ApiEntityReference
export convert, ApiError
export convert, ApiErrorBase
export convert, AvailabilitySetListResult
export convert, AvailabilitySetProperties
export convert, BootDiagnostics
export convert, BootDiagnosticsInstanceView
export convert, Caching
export convert, ComputeLongRunningOperationProperties
export convert, CreateOption
export convert, DataDisk
export convert, DataDiskImage
export convert, DiagnosticsProfile
export convert, DiskEncryptionSettings
export convert, DiskInstanceView
export convert, HardwareProfile
export convert, ImageDataDisk
export convert, ImageListResult
export convert, ImageOSDisk
export convert, ImageProperties
export convert, ImageStorageProfile
export convert, InnerError
export convert, InstanceViewStatus
export convert, KeyVaultKeyReference
export convert, KeyVaultSecretReference
export convert, LinuxConfiguration
export convert, ListUsagesResult
export convert, MaintenanceRedeployStatus
export convert, NetworkInterfaceReferenceProperties
export convert, NetworkProfile
export convert, OSDisk
export convert, OSDiskImage
export convert, OSProfile
export convert, OperationStatusResponse
export convert, Plan
export convert, PurchasePlan
export convert, Resource
export convert, ResourceSku
export convert, ResourceSkuCapabilities
export convert, ResourceSkuCapacity
export convert, ResourceSkuCosts
export convert, ResourceSkuRestrictions
export convert, ResourceSkusResult
export convert, RollingUpgradePolicy
export convert, RollingUpgradeProgressInfo
export convert, RollingUpgradeRunningStatus
export convert, RollingUpgradeStatusInfoProperties
export convert, Sku
export convert, SshConfiguration
export convert, SshPublicKey
export convert, StorageAccountType
export convert, StorageProfile
export convert, SubResource
export convert, SubResourceReadOnly
export convert, UpdateResource
export convert, UpgradePolicy
export convert, Usage
export convert, UsageName
export convert, VaultCertificate
export convert, VaultSecretGroup
export convert, VirtualHardDisk
export convert, VirtualMachineAgentInstanceView
export convert, VirtualMachineCaptureParameters
export convert, VirtualMachineCaptureResultProperties
export convert, VirtualMachineExtensionHandlerInstanceView
export convert, VirtualMachineExtensionImageProperties
export convert, VirtualMachineExtensionInstanceView
export convert, VirtualMachineExtensionProperties
export convert, VirtualMachineHealthStatus
export convert, VirtualMachineIdentity
export convert, VirtualMachineImageProperties
export convert, VirtualMachineInstanceView
export convert, VirtualMachineListResult
export convert, VirtualMachineProperties
export convert, VirtualMachineScaleSetDataDisk
export convert, VirtualMachineScaleSetExtensionListResult
export convert, VirtualMachineScaleSetExtensionProfile
export convert, VirtualMachineScaleSetExtensionProperties
export convert, VirtualMachineScaleSetIPConfigurationProperties
export convert, VirtualMachineScaleSetIdentity
export convert, VirtualMachineScaleSetInstanceView
export convert, VirtualMachineScaleSetInstanceViewStatusesSummary
export convert, VirtualMachineScaleSetListResult
export convert, VirtualMachineScaleSetListSkusResult
export convert, VirtualMachineScaleSetListWithLinkResult
export convert, VirtualMachineScaleSetManagedDiskParameters
export convert, VirtualMachineScaleSetNetworkConfigurationDnsSettings
export convert, VirtualMachineScaleSetNetworkConfigurationProperties
export convert, VirtualMachineScaleSetNetworkProfile
export convert, VirtualMachineScaleSetOSDisk
export convert, VirtualMachineScaleSetOSProfile
export convert, VirtualMachineScaleSetProperties
export convert, VirtualMachineScaleSetPublicIPAddressConfiguration
export convert, VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings
export convert, VirtualMachineScaleSetPublicIPAddressConfigurationProperties
export convert, VirtualMachineScaleSetSku
export convert, VirtualMachineScaleSetSkuCapacity
export convert, VirtualMachineScaleSetStorageProfile
export convert, VirtualMachineScaleSetUpdateIPConfigurationProperties
export convert, VirtualMachineScaleSetUpdateNetworkConfigurationProperties
export convert, VirtualMachineScaleSetUpdateNetworkProfile
export convert, VirtualMachineScaleSetUpdateOSDisk
export convert, VirtualMachineScaleSetUpdateOSProfile
export convert, VirtualMachineScaleSetUpdateProperties
export convert, VirtualMachineScaleSetUpdatePublicIPAddressConfiguration
export convert, VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties
export convert, VirtualMachineScaleSetUpdateStorageProfile
export convert, VirtualMachineScaleSetUpdateVMProfile
export convert, VirtualMachineScaleSetVMExtensionsSummary
export convert, VirtualMachineScaleSetVMInstanceIDs
export convert, VirtualMachineScaleSetVMInstanceRequiredIDs
export convert, VirtualMachineScaleSetVMInstanceView
export convert, VirtualMachineScaleSetVMListResult
export convert, VirtualMachineScaleSetVMProfile
export convert, VirtualMachineScaleSetVMProperties
export convert, VirtualMachineSize
export convert, VirtualMachineSizeListResult
export convert, VirtualMachineStatusCodeCount
export convert, WinRMConfiguration
export convert, WinRMListener
export convert, WindowsConfiguration
export convert, AvailabilitySet
export convert, Image
export convert, ImageReference
export convert, ManagedDiskParameters
export convert, NetworkInterfaceReference
export convert, RollingUpgradeStatusInfo
export convert, VirtualMachine
export convert, VirtualMachineCaptureResult
export convert, VirtualMachineExtension
export convert, VirtualMachineExtensionImage
export convert, VirtualMachineImageResource
export convert, VirtualMachineScaleSet
export convert, VirtualMachineScaleSetExtension
export convert, VirtualMachineScaleSetIPConfiguration
export convert, VirtualMachineScaleSetNetworkConfiguration
export convert, VirtualMachineScaleSetUpdate
export convert, VirtualMachineScaleSetUpdateIPConfiguration
export convert, VirtualMachineScaleSetUpdateNetworkConfiguration
export convert, VirtualMachineScaleSetVM
export convert, VirtualMachineImage

# export operations
export convert, AvailabilitySetsApi, ImagesApi, SkusApi, UsageApi, VirtualMachineExtensionImagesApi, VirtualMachineExtensionsApi, VirtualMachineImagesApi, VirtualMachineScaleSetExtensionsApi, VirtualMachineScaleSetRollingUpgradesApi, VirtualMachineScaleSetVMsApi, VirtualMachineScaleSetsApi, VirtualMachineSizesApi, VirtualMachinesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using DimensionalData, Test, Unitful
using Unitful: s

dim = Coord([(1.0,1.0,1.0), (1.0,2.0,2.0), (3.0,4.0,4.0), (1.0,3.0,4.0)], (X(), Y(), Z()))
da = DimArray(0.1:0.1:0.4, dim)
da2 = DimArray((0.1:0.1:0.4) * (1:1:3)', (dim, Ti(1s:1s:3s)))

@testset "regular indexing" begin
    @test da[Coord()] === da[Coord(:)] === da
    @test da[Coord([1, 2])] == [0.1, 0.2]
    @test da[Coord(4)] == 0.4
    @test da2[Coord(4), Ti(3)] ≈ 1.2
end

@testset "coord selector indexing" begin
    @test da[Coord(:, :, :)] == [0.1, 0.2, 0.3, 0.4]
    @test da[Coord(Between(1, 5), :, At(4.0))] == [0.3, 0.4]
    @test da[Coord(:, Between(1, 3), :)] == [0.1, 0.2]
    @test da2[Ti(1s), Coord(:, Between(1, 3), :)] ≈ [0.1, 0.2]
end

@testset "coord dimension indexing" begin
    @test da[Coord(Z(At(1.0)), Y(Between(1, 3)))] == [0.1]
end

@test index(da[Coord(:, Between(1, 3), :)], Coord) == [(1.0,1.0,1.0), (1.0,2.0,2.0)]

@test bounds(da) == (((1.0, 3.0), (1.0, 4.0), (1.0, 4.0)),)


text:
# error bars

#makie

using CairoMakie

f = Figure()
Axis(f[1, 1])

xs = 0:0.5:10
ys = 0.5 .* sin.(xs)

lowerrors = fill(0.1, length(xs))
higherrors = LinRange(0.1, 0.4, length(xs))

errorbars!(xs, ys, lowerrors, higherrors,
    color = range(0, 1, length = length(xs)),
    whiskerwidth = 10, linewidth=2, colormap = :tableau_hue_circle)

# plot position scatters so low and high errors can be discriminated
CairoMakie.scatter!(xs, ys, markersize = 4, color = :black)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using LinearAlgebra
using NSDEBase
using RungeKutta
using Test

u0 = [1.0, 1.0, 1.0]
tspan = (0.0, 1e-2)
function f!(du, u, t)
    du .= Matrix(1.0I, 3, 3) * u
end
problem = IVP(f!, u0, tspan)

@testset "Explicit" begin
    solver = ExplicitEuler(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = ExplicitMidpoint(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Heun2(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Ralston2(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Heun3(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Kutta3(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Ralston3(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = SSPRK3(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = RK4(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = Rule38(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = HeunEuler(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = F45(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = DP54(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = V65(h = 1e-2)
    @test solver isa ExplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
end

@testset "Implicit" begin
    solver = ImplicitEuler(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = ImplicitMidpoint(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = CrankNicolson(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = SDIRK3(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = GL4(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = GL6(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = LobattoIIIA4(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = LobattoIIIB2(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = LobattoIIIB4(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = LobattoIIIC2(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = LobattoIIIC4(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = RadauIA3(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = RadauIA5(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = RadauIIA3(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
    solver = RadauIIA5(h = 1e-2)
    @test solver isa ImplicitRungeKuttaSolver
    solution = solve(problem, solver)
    @test solution isa RungeKuttaSolution
end

@testset "Stability" begin
    @test ℛ(0.0, Euler().tableau) == 1.0
end


text:
using Distributions
using Base.Test

v = 7.0
S = eye(2)
S[1, 2] = S[2, 1] = 0.5

W = Wishart(v,S)
IW = InverseWishart(v,S)

for d in [W,IW]
    @test size(d) == size(rand(d))
    @test length(d) == length(rand(d))
    @test typeof(d)(params(d)...) == d
end

@test_approx_eq_eps mean(rand(W,100000)) mean(W) 0.1
@test_approx_eq_eps mean(rand(IW,100000)) mean(IW) 0.1

v = 3.0

@test_approx_eq_eps pdf(Wishart(v,S), S) 0.04507168 1e-8
@test_approx_eq_eps pdf(Wishart(v,S), inv(S)) 0.01327698 1e-8
@test_approx_eq_eps pdf(Wishart(v,inv(S)),S) 0.0148086 1e-8
@test_approx_eq_eps pdf(Wishart(v,inv(S)),inv(S)) 0.01901462 1e-8

@test_approx_eq_eps pdf(InverseWishart(v,S), S) 0.04507168 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,S), inv(S)) 0.006247377 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,inv(S)),S)  0.03147137 1e-8
@test_approx_eq_eps pdf(InverseWishart(v,inv(S)),inv(S)) 0.01901462 1e-8


text:
const CMat = Ptr{Cvoid}

"""
    A Petsc matrix, actually just a pointer to the actual C matrix
"""
struct PetscMat
    ptr::Ref{CMat}

    PetscMat() = new(Ref{CMat}())
end

# allows us to pass PetscMat objects directly into CMat ccall signatures
Base.cconvert(::Type{CMat}, mat::PetscMat) = mat.ptr[]

"""
MatSetValue(mat::PetscMat, i::PetscInt, j::PetscInt, v::PetscScalar, mode::InsertMode)

Wrapper to `MatSetValue`. Indexing starts at 0 (as in PETSc).

# Implementation
For an unknow reason, calling PETSc.MatSetValue leads to an "undefined symbol: MatSetValue" error.
So this wrapper directly call MatSetValues (anyway, this is what is done in PETSc...)
"""
function MatSetValue(mat::PetscMat, i::PetscInt, j::PetscInt, v::PetscScalar, mode::InsertMode)
    MatSetValues(mat, PetscIntOne, [i], PetscIntOne, [j], [v], mode)
end

function MatSetValue(mat::PetscMat, i::Integer, j::Integer, v::Number, mode::InsertMode)
    MatSetValue(mat, PetscInt(i), PetscInt(j), PetscScalar(v), mode)
end

"""
MatSetValues(mat::PetscMat, nI::PetscInt, I::Vector{PetscInt}, nJ::PetscInt, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode)

Wrapper to `MatSetValues`. Indexing starts at 0 (as in PETSc)
"""
function MatSetValues(mat::PetscMat, nI::PetscInt, I::Vector{PetscInt}, nJ::PetscInt, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode)
    error = ccall((:MatSetValues, libpetsc), PetscErrorCode,
        (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode),
        mat, nI, I, nJ, J, V, mode)
    @assert iszero(error)
end

function MatSetValues(mat::PetscMat, I::Vector{PetscInt}, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode)
    MatSetValues(mat, PetscInt(length(I)), I, PetscInt(length(J)), J, V, mode)
end

function MatSetValues(mat::PetscMat, I, J, V, mode::InsertMode)
    MatSetValues(mat, PetscInt.(I), PetscInt.(collect(J)), PetscScalar.(V), mode)
end

"""
    MatCreate(comm::MPI.Comm, mat::PetscMat)

Wrapper to MatCreate
"""
function MatCreate(comm::MPI.Comm, mat::PetscMat)
    error = ccall((:MatCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CMat}), comm, mat.ptr)
    @assert iszero(error)
end

function MatCreate(comm::MPI.Comm = MPI.COMM_WORLD)
    mat = PetscMat()
    MatCreate(comm, mat)
    return mat
end

"""
    MatSetSizes(mat::PetscMat, nrows_loc, ncols_loc, nrows_glo, ncols_glo)

Wrapper to MatSetSizes
"""
function MatSetSizes(mat::PetscMat, nrows_loc, ncols_loc, nrows_glo, ncols_glo)
    nr_loc = PetscInt(nrows_loc)
    nc_loc = PetscInt(ncols_loc)
    nr_glo = PetscInt(nrows_glo)
    nc_glo = PetscInt(ncols_glo)
    error = ccall((:MatSetSizes, libpetsc),
                PetscErrorCode,
                (CMat, PetscInt, PetscInt, PetscInt, PetscInt),
                mat, nr_loc, nc_loc, nr_glo, nc_glo
            )
    @assert iszero(error)
end

"""
    MatSetUp(mat::PetscMat)

Wrapper to MatSetUp
"""
function MatSetUp(mat::PetscMat)
    error = ccall((:MatSetUp, libpetsc), PetscErrorCode, (CMat,), mat)
    @assert iszero(error)
end

"""
    MatSetFromOptions(mat::PetscMat)

Wrapper to MatSetFromOptions
"""
function MatSetFromOptions(mat::PetscMat)
    error = ccall((:MatSetFromOptions, libpetsc), PetscErrorCode, (CMat,), mat)
    @assert iszero(error)
end

"""
    MatGetOwnershipRange(mat::PetscMat)

Wrapper to `MatGetOwnershipRange`

The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor.

# Warning
`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row
handled by the local processor.
"""
function MatGetOwnershipRange(mat::PetscMat)
    rstart = Ref{PetscInt}(0)
    rend = Ref{PetscInt}(0)

    error = ccall((:MatGetOwnershipRange, libpetsc), PetscErrorCode, (CMat, Ref{PetscInt}, Ref{PetscInt}), mat, rstart, rend)
    @assert iszero(error)

    return rstart[], rend[]
end

"""
    Wrapper to `MatAssemblyBegin`
"""
function MatAssemblyBegin(mat::PetscMat, type::MatAssemblyType)
    error = ccall((:MatAssemblyBegin, libpetsc), PetscErrorCode, (CMat, MatAssemblyType), mat, type)
    @assert iszero(error)
end

"""
    Wrapper to `MatAssemblyEnd`
"""
function MatAssemblyEnd(mat::PetscMat, type::MatAssemblyType)
    error = ccall((:MatAssemblyEnd, libpetsc), PetscErrorCode, (CMat, MatAssemblyType), mat, type)
    @assert iszero(error)
end

"""
    Wrapper to MatCreateVecs
"""
function MatCreateVecs(mat::PetscMat, vecr::PetscVec, veci::PetscVec)
    error = ccall((:MatCreateVecs, libpetsc), PetscErrorCode, (CMat, Ptr{CVec}, Ptr{CVec}), mat, vecr.ptr, veci.ptr)
    @assert iszero(error)
end

function MatCreateVecs(mat::PetscMat)
    vecr = PetscVec(); veci = PetscVec()
    MatCreateVecs(mat, vecr, veci)
    return vecr, veci
end

"""
    MatView(mat::PetscMat, viewer::PetscViewer = PetscViewerStdWorld())

Wrapper to `MatView`
"""
function MatView(mat::PetscMat, viewer::PetscViewer = PetscViewerStdWorld())
    error = ccall((:MatView, libpetsc), PetscErrorCode, (CMat, CViewer), mat, viewer)
    @assert iszero(error)
end

"""
    MatDestroy(mat::PetscMat)

Wrapper to MatDestroy
"""
function MatDestroy(mat::PetscMat)
    error = ccall((:MatDestroy, libpetsc), PetscErrorCode, (Ptr{CMat},), mat.ptr)
    @assert iszero(error)
end

"""
    MatCreateComposite(mat::PetscMat, smat::PetscMat)

Wrapper to `MatCompositeAddMat`
"""
# comm::MPI.Comm = MPI.COMM_WORLD
function MatCreateComposite(comm::MPI.Comm, nmat::PetscInt, mats::Vector{PetscMat}, mat::PetscMat)
    error = ccall((:MatCompositeAddMat, libpetsc), PetscErrorCode, (MPI.MPI_Comm, PetscInt, Ptr{CMat}, Ptr{CMat}),
                  comm, nmat, [x.ptr for x in mats], mat.ptr)
    @assert iszero(error)
end

function MatCreateComposite(comm::MPI.Comm, nmat::Integer, mats::Vector{PetscMat}, mat::PetscMat)
    MatCreateComposite(comm, PetscInt(nmat), mats, mat)
end

function MatCreateComposite(nmat::PetscInt, mats::Vector{PetscMat}, mat::PetscMat)
    MatCreateComposite(MPI.COMM_WORLD, nmat, mats, mat)
end

function MatCreateComposite(nmat::Integer, mats::Vector{PetscMat}, mat::PetscMat)
    MatCreateComposite(MPI.COMM_WORLD, PetscInt(nmat), mats, mat)
end


text:
module MonomialBasisTests

using Test
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials

xi = Point(2,3)
np = 5
x = fill(xi,np)

# order 0 degenerated case

order = 0
V = Float64
G = gradient_type(V,xi)
H = gradient_type(G,xi)
b = MonomialBasis{2}(V,order)
@test get_order(b) == 0
@test get_orders(b) == (0,0)

v = V[1.0,]
g = G[(0.0, 0.0),]
h = H[(0.0, 0.0, 0.0, 0.0),]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
Hbx = repeat(permutedims(h),np)
test_field(b,x,bx,grad=∇bx,hessian=Hbx)

# Real-valued Q space with isotropic order

order = 1
V = Float64
G = gradient_type(V,xi)
H = gradient_type(G,xi)
b = MonomialBasis{2}(V,order)

v = V[1.0, 2.0, 3.0, 6.0]
g = G[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (3.0, 2.0)]
h = H[(0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 1.0, 1.0, 0.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
Hbx = repeat(permutedims(h),np)
test_field(b,x,bx,grad=∇bx,hessian=Hbx)

# Real-valued Q space with an isotropic order

orders = (1,2)
V = Float64
G = gradient_type(V,xi)
b = MonomialBasis{2}(V,orders)

v = V[1.0, 2.0, 3.0, 6.0, 9.0, 18.0]
g = G[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (3.0, 2.0), (0.0, 6.0), (9.0, 12.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

# Vector-valued Q space with isotropic order

order = 1
V = VectorValue{3,Float64}
G = gradient_type(V,xi)
H = gradient_type(G,xi)
b = MonomialBasis{2}(V,order)

v = V[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],
      [2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0],
      [3.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 3.0],
      [6.0, 0.0, 0.0], [0.0, 6.0, 0.0], [0.0, 0.0, 6.0]]

g = G[[0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0],
      [0.0 0.0 0.0; 0.0 0.0 0.0], [1.0 0.0 0.0; 0.0 0.0 0.0],
      [0.0 1.0 0.0; 0.0 0.0 0.0], [0.0 0.0 1.0; 0.0 0.0 0.0],
      [0.0 0.0 0.0; 1.0 0.0 0.0], [0.0 0.0 0.0; 0.0 1.0 0.0],
      [0.0 0.0 0.0; 0.0 0.0 1.0], [3.0 0.0 0.0; 2.0 0.0 0.0],
      [0.0 3.0 0.0; 0.0 2.0 0.0], [0.0 0.0 3.0; 0.0 0.0 2.0]]

h = H[
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0),
      (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
Hbx = repeat(permutedims(h),np)
test_field(b,x,bx,grad=∇bx,hessian=Hbx)

# Vector-valued Q space with an-isotropic order

orders = (1,2)
V = VectorValue{2,Float64}
G = gradient_type(V,xi)
b = MonomialBasis{2}(V,orders)

v = V[
 (1.0, 0.0), (0.0, 1.0), (2.0, 0.0), (0.0, 2.0),
 (3.0, 0.0), (0.0, 3.0), (6.0, 0.0), (0.0, 6.0),
 (9.0, 0.0), (0.0, 9.0), (18.0, 0.0), (0.0, 18.0)]

g = G[
  (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0),
  (1.0, 0.0, 0.0, 0.0), (0.0, 0.0, 1.0, 0.0),
  (0.0, 1.0, 0.0, 0.0), (0.0, 0.0, 0.0, 1.0),
  (3.0, 2.0, 0.0, 0.0), (0.0, 0.0, 3.0, 2.0),
  (0.0, 6.0, 0.0, 0.0), (0.0, 0.0, 0.0, 6.0),
  (9.0, 12.0, 0.0, 0.0), (0.0, 0.0, 9.0, 12.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

# Real-valued P space

order = 1
V = Float64
G = gradient_type(V,xi)
filter = (e,o) -> sum(e) <= o
b = MonomialBasis{2}(V,order,filter)

v = V[1.0, 2.0, 3.0]
g = G[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

# Vector-valued P space

order = 1
V = VectorValue{3,Float64}
G = gradient_type(V,xi)
filter = (e,o) -> sum(e) <= o
b = MonomialBasis{2}(V,order,filter)

v = V[[1.0; 0.0; 0.0], [0.0; 1.0; 0.0], [0.0; 0.0; 1.0],
      [2.0; 0.0; 0.0], [0.0; 2.0; 0.0], [0.0; 0.0; 2.0],
      [3.0; 0.0; 0.0], [0.0; 3.0; 0.0], [0.0; 0.0; 3.0]]

g = G[[0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0],
      [0.0 0.0 0.0; 0.0 0.0 0.0], [1.0 0.0 0.0; 0.0 0.0 0.0],
      [0.0 1.0 0.0; 0.0 0.0 0.0], [0.0 0.0 1.0; 0.0 0.0 0.0],
      [0.0 0.0 0.0; 1.0 0.0 0.0], [0.0 0.0 0.0; 0.0 1.0 0.0],
      [0.0 0.0 0.0; 0.0 0.0 1.0]]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

order = 1
b = MonomialBasis{1}(Float64,order)
@test evaluate(b,Point{1,Float64}[(0,),(1,)]) == [1.0 0.0; 1.0 1.0]

b = MonomialBasis{0}(VectorValue{2,Float64},order)
@test evaluate(b,Point{0,Float64}[(),()]) == VectorValue{2,Float64}[(1.0, 0.0) (0.0, 1.0); (1.0, 0.0) (0.0, 1.0)]

end # module


text:
# Contains types for time axis manipulations
@def_structequal struct TcspcTimeAxis
    # Stores time bin information for the collected axes
    nbins::Unsigned
    tcal::AbstractFloat
    ψ::AbstractFloat
    # Define internal constructor
    function TcspcTimeAxis(tcal::AbstractFloat,nbins::Integer,
        ψ::AbstractFloat)
        # Round nbins to nearest power of two
        if nbins <= 0
            error("nbins ($nbins) must be greater than 0")
        end
        nbins = nextpow(2,nbins)
        if tcal < 0
            error("tcal must be larger than 0.0")
        end
        new(nbins,tcal,ψ)
    end
end

# Multiple dispatch with inputs of start, stop, and nbins
function TcspcTimeAxis(start::AbstractFloat,
    stop::AbstractFloat,nbins::Integer,
    ψ::AbstractFloat)
    if start >= stop
        error("start ($start) must be less than stop ($stop)")
    end
    # Determine tcal
    tcal = (stop - start) / nbins
    TcspcTimeAxis(tcal,nbins,ψ)
end

# Helper methods with time axis
function chan2time(tax::TcspcTimeAxis,chan::Integer)
    # Converts channel number to time point along axis
    if chan <= 0
        error("chan must be greater than 0")
    end
    # Return to middle of bin
    chan * (tax.tcal / 2)
end

function time2chan(tax::TcspcTimeAxis,t::AbstractFloat)
    if t < 0
        error("t must be greater than or equal to 0")
    end
    # Generate histogram of data
    h = fit(Histogram,t,time_bins(tax))
    # Find non-zero elements
    findfirst(h.weights)
end

time_bins(t::TcspcTimeAxis) = range(0,step=t.tcal,length=t.nbins+1)

t_vector(t::TcspcTimeAxis) = range(t.tcal,step=t.tcal,length=t.nbins)

function t_vector(tcal::AbstractFloat,nbins::Integer)
    tax = TcspcTimeAxis(tcal,nbins,44.7) # placeholder ψ
    range(tax.tcal,step=tax.tcal,length=nbins)
end

η(t::TcspcTimeAxis) = 3 * (cosd(t.ψ) ^ 2) - 1


text:
using SuiteSparseGraphBLAS.SparseArrayCompat
@testset "Matching SparseMatrixCSC" begin
    A = sprand(5, 5, 0.35)
    B = SparseMatrixGB(A)
    C = sprand(5, 5, 0.35)
    D = SparseMatrixGB(C)
    for i ∈ eachindex(A)
        @test A[i] == B[i]
    end
    # broadcasting tests
    @test A .+ C ≈ B .+ D
    @test A .- C ≈ B .- D
    @test A .* C ≈ B .* D
    @test min.(A, C) ≈ min.(B, D)
    @test max.(A, C) ≈ max.(B, D)

    # map tests
    @test map(one, A) == map(one, B)
    @test map(sin, A) ≈ map(sin, B)
    @test map(identity, A) == map(identity, B)
    @test map(cos, A) ≈ map(cos, B)
    @test map((x)-> 1.5*x^2 + 3.0, A) ≈ map((x)-> 1.5*x^2 + 3.0, B)
end


text:
qrf!(A,τ) = _qrf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
_qrf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
    LAPACK.geqrf!(A,τ)
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
    LAPACK.ormqr!('L','T',A,τ,B)
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
    LAPACK.ormqr!('L','C',A,τ,B)
apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)

qlf!(A,τ) = _qlf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
_qlf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
    LAPACK.geqlf!(A,τ)
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
    LAPACK.ormql!('L','T',A,τ,B)
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
    LAPACK.ormql!('L','C',A,τ,B)
apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)

function qr!(A::BlockBandedMatrix{T}) where T
    l,u = blockbandwidths(A)
    M,N = nblocks(A)
    bs = M < N ? BlockSizes((cumulsizes(blocksizes(A),1),)) : BlockSizes((cumulsizes(blocksizes(A),2),))
    τ = PseudoBlockVector{T}(undef, bs)
    for K = 1:min(N,M)
        KR = Block.(K:min(K+l,M))
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        qrf!(V,t)
        for J = K+1:min(K+u,N)
            apply_qr!(V, t, view(A,KR,Block(J)))
        end
    end
    QR(A,τ.blocks)
end

function ql!(A::BlockBandedMatrix{T}) where T
    l,u = blockbandwidths(A)
    M,N = nblocks(A)

    bs = if M < N
        throw(ArgumentError("Wide block-QL not implented"))
    else
        BlockSizes((cumulsizes(blocksizes(A),2),))
    end
    τ = PseudoBlockVector{T}(undef, bs)

    for K = N:-1:max(N - M + 1,1)
        μ = M+K-N
        KR = Block.(max(K-u,1):μ)
        V = view(A,KR,Block(K))
        t = view(τ,Block(K-N+min(M,N)))
        qlf!(V,t)
        for J = K-1:-1:max(K-l,1)
            apply_ql!(V, t, view(A,KR,Block(J)))
        end
    end
    QL(A,τ.blocks)
end

qr(A::BlockBandedMatrix) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
ql(A::BlockBandedMatrix) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))

qr(A::BandedBlockBandedMatrix) = qr(BlockBandedMatrix(A))
ql(A::BandedBlockBandedMatrix) = ql(BlockBandedMatrix(A))

function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
    Q = parent(adjQ)
    A = Q.factors
    l,u = blockbandwidths(A)
    N,M = nblocks(A)
    # impose block structure
    bs = BlockSizes((cumulsizes(blocksizes(A),1),))
    τ  = PseudoBlockArray(Q.τ, bs)
    B = PseudoBlockArray(Bin, bs)
    for K = 1:min(N,M)
        KR = Block.(K:min(K+l,N))
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        apply_qr!(V, t, view(B,KR))
    end
    Bin
end

function lmul!(adjQ::Adjoint{<:Any,<:QLPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
    Q = parent(adjQ)
    A = Q.factors
    l,u = blockbandwidths(A)
    N,M = nblocks(A)
    # impose block structure
    bs = BlockSizes((cumulsizes(blocksizes(A),1),))
    τ  = PseudoBlockArray(Q.τ, bs)
    B = PseudoBlockArray(Bin, bs)
    for K = N:-1:1
        KR = Block.(max(1,K-u):K)
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        apply_ql!(V, t, view(B,KR))
    end
    Bin
end

# avoid LinearALgebra Strided obsession

for Typ in (:StridedVector, :StridedMatrix, :AbstractVector, :AbstractMatrix)
    @eval function ldiv!(A::QR{<:Any,<:BlockSkylineMatrix}, B::$Typ)
        lmul!(adjoint(A.Q), B)
        M,N = nblocks(A.factors)
        MN = min(M,N)
        V = view(A.factors,Block.(1:MN), Block.(1:MN))
        materialize!(Ldiv(UpperTriangular(V), view(B,1:size(V,1),:)))
        B
    end
end

function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractVector)
    lmul!(adjoint(A.Q), B)
    materialize!(Ldiv(LowerTriangular(A.factors), B))
end

function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractMatrix)
    lmul!(adjoint(A.Q), B)
    materialize!(Ldiv(LowerTriangular(A.factors), B))
end

\(A::AbstractBlockBandedMatrix, b::AbstractVecOrMat) = qr(A)\b  # use QR because LU would be a _mess_ to implement


text:
@testset "dictionary difference" begin
    modvals = Dict(:x1=>1, :x2=>[0.0, 1.0])
    addvals = Dict()
    remvals = Dict(:z=>3)
    a = DictDifference(modvals, addvals, remvals)
    b = DictDifference(modvals, addvals, remvals)
    c = DictDifference(modvals, addvals, remvals)

    @testset "constructors" begin
        @test isa(a, DictDifference)
        @test a.modvals == modvals && a.addvals == addvals && a.remvals == remvals
    end

    @testset "equality operator" begin
        @test a == a
        @test a == b && b == a
        @test a == b && b == c && a == c
    end

    @testset "hash" begin
        @test hash(a) == hash(a)
        @test a == b && hash(a) == hash(b)
    end

    @testset "accessors" begin
        @test modified(a) == modvals
        @test added(a) == addvals
        @test removed(a) == remvals
    end

    @testset "difference" begin
        a = Dict(:x=>1)
        b = Dict(:x=>Dict(:s=>1, :t=>2))
        c = Dict(:x=>[1, 2])
        d = Dict(:x=>2, :y=>1)
        ev = []
        evi = Int[]
        ed = Dict()
        @test diff(a, a) == DictDifference(Dict(:x=>0), ed, ed)
        @test diff(b, b) == DictDifference(
            Dict(:x => DictDifference(Dict(:s=>0, :t=>0), ed, ed)), ed, ed
        )
        @test diff(c, c) == DictDifference(
            Dict(:x => VectorDifference([1, 2], evi, evi, [0, 0], ev, ev)), ed, ed
        )
        @test diff(a, d) == DictDifference(Dict(:x=>-1), Dict(:y=>1), ed)
        @test diff(d, a) == DictDifference(Dict(:x=>1), ed, Dict(:y=>1))
    end

    @testset "math" begin
        a = Dict(:x=>1, :y=>2)
        b = Dict(:y=>4, :z=>8)
        d = DictDifference(Dict(:y=>2), Dict(:z=>8), Dict(:x=>1))
        @test a + d == b
        @test d + a == b
    end
end


text:
# our iterator type
type OrderedSubsets{T}
    members::Vector{T}
    valueFunction
    pq::Collections.PriorityQueue
    started::Bool
end

"Enumerate all the subsets ordered by the provided value function."
function orderedsubsets{T}(members::Vector{T}, f, order=:descending)
    orderby = x->-f([x])
    pqOrder = Base.Order.Reverse
    if order == :ascending
        orderby = x->f([x])
        pqOrder = Base.Order.Forward
    end
    OrderedSubsets(
        sort(members, by=orderby),
        f,
        Collections.PriorityQueue(Vector{Int64}, Float64, pqOrder),
        false
    )
end

Base.eltype(it::OrderedSubsets) = Vector{eltype(it.members)}
function Base.length(it::OrderedSubsets)
    length(it.members) < 40 ? 2^length(it.members) : Inf
end

Base.start{T}(it::OrderedSubsets{T}) = nothing

function Base.next{T}(it::OrderedSubsets{T}, state)
    if !it.started
        it.started = true
        Collections.enqueue!(it.pq, [1], it.valueFunction([it.members[1]]))
        return T[],nothing
    end

    # get the next best subset
    nextSubset = Collections.dequeue!(it.pq)

    # move this subset to look for the next "last member"
    if nextSubset[end] < length(it.members)
        nextSubsetInc = copy(nextSubset)
        nextSubsetInc[end] += 1
        Collections.enqueue!(it.pq, nextSubsetInc, it.valueFunction(T[it.members[i] for i in nextSubsetInc]))
    end

    # add the grown version of this subset
    if nextSubset[end] < length(it.members)
        nextSubsetGrown = [nextSubset; nextSubset[end]+1]
        Collections.enqueue!(it.pq, nextSubsetGrown, it.valueFunction(T[it.members[i] for i in nextSubsetGrown]))
    end

    T[it.members[i] for i in nextSubset],nothing
end

Base.done(it::OrderedSubsets, state) = isempty(it.pq) && it.started

# our iterator type
type ESKernelSubsets
    itr
    complement::Bool
    variances::Vector{Float64}
    lastSubset
    M
    numDone
    numTotal
end

function eskernelSubsetWeight(x, M)
    s = length(x)
    if s == M || s == 0
        return 1e12
    end
    (M-1)/(binomial(M,s)*s*(M-s))
end

"Enumerate all the subsets in descending order by the provided value function."
function eskernelsubsets(members, variances)
    M = length(members)

    # ensure all the variances are unique
    variances = variances .+ collect(M:-1:1)/1e5

    function subsetValue(x)

        s = length(x)
        if s == M || s == 0
            return 1e12
        elseif s > M/2
            return -1
        elseif s == 1
            return 1e10 + minimum(variances[x])
        end

        w = (M-1)/(binomial(M,s)*s*(M-s))
        minimum(variances[x])*w
    end

    ESKernelSubsets(
        orderedsubsets(members, subsetValue),
        false,
        variances,
        Int64[],
        M,
        0,
        M < 40 ? 2^M : Inf
    )
end

Base.eltype(it::ESKernelSubsets) = Vector{Float32}
function Base.length(it::ESKernelSubsets)
    it.numTotal
end

function Base.start(it::ESKernelSubsets)
    start(it.itr)
end

function Base.next(it::ESKernelSubsets, state)
    if !it.complement
        subsetElements,tmp = next(it.itr, nothing)

        it.lastSubset = subsetElements
        it.complement = true
        out = zeros(Float32, it.M)
        out[subsetElements] = 1
    else
        out = ones(Float32, it.M)
        out[it.lastSubset] = 0
        it.complement = false
    end
    it.numDone += 1

    return (out,eskernelSubsetWeight(find(out), it.M)),nothing
end

Base.done(it::ESKernelSubsets, state) = it.numDone >= it.numTotal