text:
module MPolyFact

using Hecke
import Hecke: Nemo, @vprint, @hassert, @vtime

add_verbose_scope(:AbsFact)
add_assert_scope(:AbsFact)

Hecke.example("mfactor.jl")

function Hecke.norm(f::MPolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(Kx)
  n = nvars(Kx)
  Qx, x = PolynomialRing(QQ, [String(x) for x= symbols(Kx)])
  Qxy, y = PolynomialRing(Qx, "y")
  gg = [MPolyBuildCtx(Qx) for i=1:degree(K)]
  for (c, e) = zip(coeffs(f), exponent_vectors(f))
    for i=0:degree(K)-1
      d = coeff(c, i)
      if !iszero(d)
        push_term!(gg[i+1], d, e)
      end
    end
  end
  g = Qxy(map(finish, gg))
  return resultant(g, K.pol(y))
end

function Hecke.lead(f::fmpq_mpoly)
  return first(coeffs(f))
end

function Hecke.ismonic(f::fmpq_mpoly)
  return isone(lead(f))
end

#dodgy
function (k::Nemo.GaloisField)(a::fmpq)
  return k(numerator(a))//k(denominator(a))
end
function (k::Nemo.FqNmodFiniteField)(a::fmpq)
  return k(numerator(a))//k(denominator(a))
end
function (R::FmpzMPolyRing)(f::fmpq_mpoly)
  return map_coeffs(ZZ, f, parent = R)
end

function Hecke.representation_matrix(a::ResElem{<:PolyElem})
  R = parent(a)
  S = base_ring(base_ring(R))
  n = degree(modulus(R))
  m = zero_matrix(S, n, n)
  for i=0:n-1
    m[1, i+1] = coeff(a.data, i)
  end
  for j=2:n
    a *= gen(R)
    for i=0:n-1
      m[j, i+1] = coeff(a.data, i)
    end
  end
  return m
end

function Hecke.preimage(phi::Nemo.FinFieldMorphism, x::FinFieldElem)
  return preimage_map(phi)(x)
end

#should be in AA/ Nemo
function Nemo.canonical_unit(a::SeriesElem)
  iszero(a) && return one(parent(a))
  v = valuation(a)
  v == 0 && return a
  v > 0 && return shift_right(a, v)
  return shift_left(a, -v)
end

#TODO: this is for rings, not for fields, maybe different types?
function Base.gcd(a::T, b::T) where {T <: SeriesElem}
  iszero(a) && iszero(b) && return a
  iszero(a) && return gen(parent(a))^valuation(b)
  iszero(b) && return gen(parent(a))^valuation(a)
  return gen(parent(a))^min(valuation(a), valuation(b))
end

function Base.lcm(a::T, b::T) where {T <: SeriesElem}
  iszero(a) && iszero(b) && return a
  iszero(a) && return a
  iszero(b) && return b
  return gen(parent(a))^max(valuation(a), valuation(b))
end

Hecke.inv(phi :: Nemo.FinFieldMorphism) = preimage_map(phi)

#= AbsFact: strategy
 f in Q[x,y] = Q[x][y], not necc. monic, irre. over Q

 find evaluation x <- x+ly s.th. g = f(0, y) is square-free
 find p s.th g mod p is square-free and such that the lcm of the factors
   is small
 compute the roots up to a bound B over the F_q splitting field
   as power series
 compute the coeffs. of f (as elements in Q[x]) as power series
=#
mutable struct RootCtx
  f::fmpq_mpoly
  R::Array{<:SeriesElem, 1}
  o::Array{<:SeriesElem, 1} #1/f'(r)
  RP::Array{Array{<:SeriesElem , 1}, 1}  #root powers
  t::Int

  all_R::Array{<:SeriesElem, 1} #all roots - if different from R

  function RootCtx(f::fmpq_mpoly, r::Array{<:RingElem, 1}, t::Int = 0)
    @assert nvars(parent(f)) == 2

    s = parent(r[1])

    S = PowerSeriesRing(s, 2, "s")[1]
    l = new()
    l.f = f
    l.R = [S(x) for x = r]
    for i=1:length(r)
      set_precision!(l.R[i], 1)
    end
    g = map_coeffs(parent(r[1]), f)
    tt = gen(S) - t
    set_precision!(tt, 1)
    l.o = [inv(evaluate(derivative(g, 1), [x, tt])) for x = l.R]
    l.t = t

    l.RP = [[one(S) for x = r], l.R]
    return l
  end
end

function root(R::RootCtx, i::Int, j::Int)
  if precision(R.R[1]) != precision(R.RP[1][1])
    o = one(parent(R.R[1]))
    R.RP = [[set_precision(o, precision(R.R[1])) for x = R.R], copy(R.R)]
  end
  if length(R.RP) > j
#    @assert R.RP[j+1][i] == R.R[i]^j
    return (R.RP[j+1][i])
  end
  while length(R.RP) <= j+1
    push!(R.RP, R.RP[2] .* R.RP[end])
  end

  s = R.RP[j+1][i]
#  @assert s == R.R[i]^j
  return (s)
end

#TODO: in Nemo, rename to setprecision
#      fix/report series add for different length
function set_precision(a::SeriesElem, i::Int)
  b = deepcopy(a)
  set_precision!(b, i)
  return b
end

function newton_lift!(R::RootCtx)

  #TODO: given that f might be sparse, do NOT compute all powers
  #      of the roots, only those needed - and this in an "optimized"
  #      way
  S = parent(R.R[1])
  T = base_ring(S)
#  set_precision!(S, 2*precision(S))
  S.prec_max = 2*precision(R.R[1])+1

  t = gen(S) - R.t

  for i = 1:length(R.R)
    a = R.R[i]
    o = R.o[i]
    set_precision!(a, 2*precision(a))
    set_precision!(o, precision(a))
  end
  set_precision!(t, precision(R.R[1])+1)

  for i=1:length(R.R)
    a = R.R[i]
    o = R.o[i]
    ev_f = zero(S)
    ev_fs = zero(S)

    if precision(R.R[1]) > 2
      for j=1:length(R.f)
        e = exponent_vector(R.f, j)
        c = coeff(R.f, j)
        if e[1] > 0
          ev_fs += S(T(c*e[1]))*root(R, i, e[1] - 1) * t^e[2]
        end
      end

      o = R.o[i] = o*(2-o*ev_fs)
    end

    for j=1:length(R.f)
      e = exponent_vector(R.f, j)
      c = coeff(R.f, j)
      ev_f += S(T(c))*t^e[2] * root(R, i, e[1])
    end
    R.R[i] = a - ev_f*o
#    @assert evaluate(R.f, [a, t]) == ev_f
  end
end

function roots(f::fmpq_mpoly; p::Int=2^25, pr::Int = 2)
  @assert nvars(parent(f)) == 2
  #requires f to be irred. over Q - which is not tested
  #requires f(x, 0) to be same degree and irred. - should be arranged by the absolute_bivariate_factorisation

  #f in Qxy
  Zx = Hecke.Globals.Zx
  ff = map_coeffs(ZZ, f)
  #TODO: 0 might not be a good evaluation point...
  #f needs to be irreducible over Q and g square-free
  g = evaluate(ff, [gen(Zx), Zx(0)])
  @assert degree(g) == degree(f, 1)

  d = degree(g)
  best_p = p
  pc = 0
  while true
    p = next_prime(p)
    gp = factor(g, GF(p))
    if any(x->x>1, values(gp.fac))
      @vprint :AbsFact 1 "not squarefree mod $p\n"
      continue
    end
    e = lcm([degree(x) for x = keys(gp.fac)])
    if e < d
      d = e
      best_p = p
      pc = 0
    else
      pc += 1
    end

    if e == 1 || pc > div(degree(g), 2)
      @vprint :AbsFact 1 "using $best_p of degree $d\n"
      p = best_p
      break
    end
  end
  F = FiniteField(p, d)[1]
  @hassert :AbsFact 1 !iszero(discriminant(map_coeffs(F, g)))
  @vtime :AbsFact 2 r = Set(Hecke.roots(g, F))
  @assert length(r) == degree(g)
  #use action of Frobenius to lift less roots!!!
  rr = Tuple{eltype(r), Int}[]
  while length(r) > 0
    s = pop!(r)
    d = degree(minpoly(s))
    push!(rr, (s, d))
    for i=1:d-1
      s = frobenius(s)
      pop!(r, s)
    end
  end
  @vprint :AbsFact 1 "need to seriously lift $(length(rr)) elements\n"

  R = RootCtx(f, [x[1] for x = rr])
  for i=1:pr
    newton_lift!(R)
  end

  @vprint :AbsFact 1 "now Frobenius action...\n"

  S = typeof(R.R[1])[]
  all_o = []
  for i=1:length(rr)
    push!(S, R.R[i])
    o = [length(S)]
    T = S[end]
    for i=1:rr[i][2]-1
      T = deepcopy(T)
      for j=0:T.length-1
        setcoeff!(T, j, frobenius(coeff(T, j)))
      end
      push!(S, T)
      push!(o, length(S))
    end
    push!(all_o, o)
  end
  @vprint :AbsFact 2 "orbits: $all_o\n"
  R.all_R = S
  return R
end

function more_precision(R::RootCtx)
  @vtime :AbsFact 2 newton_lift!(R)
  S = typeof(R.R[1])[]
  for r = R.R
    push!(S, r)
    s = coeff(r, 0)
    @assert !iszero(s)
    @assert valuation(r) == 0
    d = degree(minpoly(s))
    T = r
    for i=1:d-1
      T = deepcopy(T)
      for j=0:T.length-1
        setcoeff!(T, j, frobenius(coeff(T, j)))
      end
      push!(S, T)
    end
  end
  R.all_R = S
end

#check with Nemo/ Dan if there are better solutions

function Hecke.leading_coefficient(f::MPolyElem, i::Int)
  g = MPolyBuildCtx(parent(f))
  d = degree(f, i)
  for (c, e) = zip(coeffs(f), exponent_vectors(f))
    if e[i] == d
      e[i] = 0
      push_term!(g, c, e)
    end
  end
  return finish(g)
end

function make_monic(f::MPolyElem, i::Int)
  d = degree(f, i)
  cf = [MPolyBuildCtx(parent(f)) for j=0:d]
  for (c, e) = zip(coeffs(f), exponent_vectors(f))
    a = e[i]
    e[i] = 0
    push_term!(cf[a+1], c, e)
  end
  df = map(finish, cf)
  for j=0:d-1
    df[j+1] *= df[d+1]^(d-1-j)
  end
  df[d+1] = one(parent(f))
  return sum(gen(parent(f), i)^j*df[j+1] for j=0:d)
end

function Hecke.content(f::MPolyElem, i::Int)
  d = degree(f, i)
  cf = [MPolyBuildCtx(parent(f)) for j=0:d]
  for (c, e) = zip(coeffs(f), exponent_vectors(f))
    a = e[i]
    e[i] = 0
    push_term!(cf[a+1], c, e)
  end
  df = map(finish, cf)
  return reduce(gcd, df)
end

function Hecke.coefficients(f::MPolyElem, i::Int)
  d = degree(f, i)
  cf = [MPolyBuildCtx(parent(f)) for j=0:d]
  for (c, e) = zip(coeffs(f), exponent_vectors(f))
    a = e[i]
    e[i] = 0
    push_term!(cf[a+1], c, e)
  end
  return map(finish, cf)
end

function combination(RC::RootCtx)
  #R is a list of roots, ie. power series over F_q (finite field)
  f = RC.f
  R = RC.all_R
  Ft = parent(R[1])
  t = gen(Ft)
  n = precision(R[1])
  @assert all(x->precision(x) == n, R)

  #ps = [[div(x^i % tn, td^i) for i = 1:n] for x = R]

  F = base_ring(Ft)
  k = degree(F)

  p = characteristic(F)
  Fp = GF(p)

  #the paper
  # https://www.math.univ-toulouse.fr/~cheze/facto_abs.m
  #seems to say that a combination works iff the sum of the
  # deg 2 terms vanish (for almost all choices of specialisations)
  #this would be even easier
  #also: it would make the lifting above trivial (precision 2)
  #deal: if monic in x of degree n and deg_y(coeff of x^(n-1)) = r,
  # then for the 1st power sum we get r as a degree bound. In the paper
  # r == 1...
  # reasoning: sum of second highest terms of factor == second highest
  # of poly => bound
  # similar would be for other power sums - if I'd compute them.
  #all vanishing results in the paper do not depend on the assumption
  #D. Rupprecht / Journal of Symbolic Computation 37 (2004) 557–574
  # doi:10.1016/S0747-7171(02)00011-1
  # https://core.ac.uk/download/pdf/82425943.pdf
  #
  # This allows (should allow) to only use the 1st power sum (trace)
  # which makes denominators easier
  #
  #for non-monic: think about interaction with any_order
  #
  # instead of LLL use direct rref over GF(q)
  # for non-monics, the generalised equation order might be better than
  # trying to scale
  #
  j = 0
  nn = matrix(Fp, 0, length(R), [])
  d = degree(f, 2)+degree(f, 1)
  lc = leading_coefficient(f, 1)
  d += degree(lc, 2)

  ld = evaluate(map_coeffs(x->F(ZZ(x)), lc), [set_precision(Ft(0), n), set_precision(gen(Ft), n)])
  @assert precision(ld) >= n
  R = R .* ld

  pow = 1
  bad = 0
  last_rank = length(R)
  while true
    j += 1
    while pow*d+j >= n
      @vprint :AbsFact 1 "need more precicsion: $n ($d, $pow, $j)\n"
      more_precision(RC)
      R = RC.all_R
      n = precision(R[1])
      if false && n > 570 #too small - but a safety valve
        error("too much n")
      end
      set_precision!(ld, n)
      R = R .* ld
      @assert precision(R[1]) >= n
    end

    mn = matrix([[Fp(coeff(coeff(x^pow, pow*d+j), lk)) for lk = 0:k-1] for x = R])

    if false && iszero(mn)
      @vprint :AbsFact 2 "found zero column, disgarding\n"
      bad += 1
      if bad > max(2, div(length(R), 2))
        pow += 1
        @vprint :AbsFact 1 "increasing power to $pow\n"
        j = 0
        bad = 0
      end
      continue
    else
      @vprint :AbsFact 2 "found non zero column\n"
    end

    nn = vcat(nn, mn)

    ke = kernel(nn)
    @vprint :AbsFact 1 "current kernel dimension: $(ke[1])\n"
    if last_rank == ke[1]
      bad += 1
      if bad > max(2, div(length(R), 2))
        pow += 1
        @vprint :AbsFact 1 "increasing power to $pow\n"
        j = 0
        bad = 0
        continue
      end
    else
      bad = 0
      last_rank = ke[1]
    end
    if ke[1] == 0 || mod(length(R), ke[1]) != 0
      continue
    end
    m = ke[2]
    z = m'*m
    if z != div(length(R), ke[1])
      @vprint :AbsFact 1 "not a equal size partition\n"
      continue
    end
    return m'
  end
end

function field(RC::RootCtx, m::MatElem)
  R = RC.all_R
  P = RC.f

  #we have roots, we need to combine roots for each row in m where the entry is pm 1
  #the coeffs then live is a number field, meaning that the elem sym functions or
  #the power sums will be needed
  #the field degree apears to be nrows(m)...

  #need primitive element, can use power sums up to #factors

  #we will ONLY find one factor, the others are galois conjugate
  #the field here is not necc. normal

  d_f = div(ncols(m), nrows(m))
  @assert ncols(m) == length(R)

  Ft = parent(R[1])
  F = base_ring(Ft)
  t = gen(Ft)
  d = precision(R[1])

  tf = divexact(total_degree(P), nrows(m))

  #TODO given that all powers are used, build them up properly
  #TODO use Frobenius (again) to save on multiplications
  #TODO invest work in one factor only - need only powers of the roots involved there
  #     the other factor is then just a division away
  #     if complete orbits are combined, use the trace (pointwise) rather than powers
  @vprint :AbsFact 2 "combining: $([findall(x->!iszero(x), collect(m[i, :])) for i=1:nrows(m)])\n"
  RP = [[set_precision(x, tf+2) for x = R]]
  @vtime :AbsFact 2 for j=2:d_f
    push!(RP, RP[1] .* RP[end])
  end
  @vtime :AbsFact 2 el = [[sum(RP[j][i] for i=1:ncols(m) if m[lj, i] != 0) for j=1:d_f] for lj=1:nrows(m)]

  #now find the degree where the coeffs actually live:
  k = 1
  for x = el
    for y = x
      d = degree(minpoly(coeff(y, valuation(y))))
      k = max(d, k)
    end
  end

  @vprint :AbsFact 1 "target field has (local) degree $k\n"

  Qq = QadicField(characteristic(F), k, 10)
  k, mk = ResidueField(Qq)
  #TODO: debug and fix in Nemo
  F.overfields = Dict{Int64,Array{Nemo.FinFieldMorphism,1}}()
  F.subfields = Dict{Int64,Array{Nemo.FinFieldMorphism,1}}()
  k.overfields = Dict{Int64,Array{Nemo.FinFieldMorphism,1}}()
  k.subfields = Dict{Int64,Array{Nemo.FinFieldMorphism,1}}()

  phi = embed(k, F)

  kt, t = PolynomialRing(k)
  kXY, (X, Y) = PolynomialRing(k, ["X", "Y"])

  el = [Hecke.power_sums_to_polynomial([kt([preimage(phi, coeff(x, i)) for i=0:precision(x)])(Y) for x = y])(X) for y = el]
  #assuming one coeff is primtive...
  #the first coeff has to be primitive as the other could be re-computed by lifting, hence cannot
  #enlarge the field
  j = 1
  local s
  while true
    s = Set([(coeff(x, j)) for x = el])
    if length(s) == length(el)
      break
    end
    j += 1
  end

  @vprint :AbsFact 2 "$j-th coeff is primitive\n"
  @vprint :AbsFact 1 "hopefully $(length(el)) degree field\n"

  QqXY, _ = PolynomialRing(Qq, 2)

  el = [map_coeffs(x->preimage(mk, x), y, parent = QqXY) for y = el]

  #at this point, hopefully, prod(el) = P/lead

  pr = 10
  while true
    pr *= 2
    @vprint :AbsFact 1  "using p-adic precision of $pr\n"
    setprecision!(Qq, pr+1)
    el = [map_coeffs(x->setprecision(x, pr), y, parent = QqXY) for y = el]

    #TODO: allow continue to lift
    @vtime :AbsFact 1 el = Main.MFactor.lift_prime_power(P, el, [0], [degree(P, 2)], pr)

    pk = prime(Qq)^pr
    p = [coeff(sum(coeff(x, j)^l for x = el), 0) for l=1:length(el)]
    p = map(lift, p)
    p = map(x->QQ(Hecke.mod_sym(x, pk)), p)
#    @show p = map(x->rational_reconstruction(x, pk), p) #all arranged to be integral, so farey not used here
#    if !all(x->x[1], p)
#      @show "reco failed"
#      continue
#    end
#    p = [x[2]//x[3] for x = p]

    pp = Hecke.power_sums_to_polynomial(p)
    if any(x->!isone(denominator(x)), coefficients(pp))
      @vprint :AbsFact 2 "poly not integral, increasing p-adic precision\n"
      continue
    end
    k, a = number_field(Hecke.power_sums_to_polynomial(p))

    @vprint :AbsFact 1  "using as number field: $k\n"

    m = matrix([[(coeff(x, j)^l) for x = el] for l=0:degree(k)-1])
    kx, x = k["x"]
    kX, (X, Y) = PolynomialRing(k, ["X", "Y"])
    B = MPolyBuildCtx(kX)
    for j=1:length(el[1])
      n = matrix([[coeff(x, j)] for x = el])
      s = solve(m, n')
      @assert all(x->iszero(coeff(s[x, 1], 1)), 1:degree(k))
      s = [rational_reconstruction(lift(coeff(s[i, 1], 0))% pk, pk) for i=1:degree(k)]
      if !all(x->x[1], s)
        break
      end
      push_term!(B, k([x[2]//x[3] for x = s]), exponent_vector(el[1], j))
    end
    q = finish(B)
    if length(q) < length(el[1])
      continue
    end
    return q
  end
end

function absolute_bivariate_factorisation(f::fmpq_mpoly)
  d = degree(f, 1)
  R = parent(f)
  x, y = gens(R)

  if degree(f, 2) < d
    @vprint :AbsFact 1 "swapping variables to have smaller degree\n"
    f = evaluate(f, [y, x])
    a, ca = absolute_bivariate_factorisation(f)
    S = parent(a)
    X, Y = gens(S)
    return evaluate(a, [Y, X]), evaluate(ca, [Y, X])
  end

  if degree(f, 2) == d && !isone(leading_coefficient(f, 1)) && isone(leading_coefficient(f, 2))
    @vprint :AbsFact 1 "swapping variables to be monic\n"
    f = evaluate(f, [y, x])
    a, ca = absolute_bivariate_factorisation(f)
    S = parent(a)
    X, Y = gens(S)
    return evaluate(a, [Y, X]), evaluate(ca, [Y, X])
  end

  Qt, t = PolynomialRing(QQ, cached = false)
  s =-1
  while true
    s += 1
    @vprint :AbsFact 1 "substitution to $s\n"
    z = evaluate(f, [t, Qt(s)])
    if degree(z) == d && issquarefree(z)
      break
    end
  end
  ff = evaluate(f, [x, y+s])
  gg = make_monic(ff, 1)
  r = roots(gg)
  z = combination(r)
  if nrows(z) == 1
    return f, one(parent(f))
  end
  a = field(r, z)
  S = parent(a)
  X, Y = gens(S)
  l = leading_coefficient(ff, 1)
  if !isone(l)
    ll = evaluate(l, [S(0), Y])
    a = evaluate(a, [X*ll, Y])
    a = divexact(a, ll^(degree(a, 1)-1))
    a = divexact(a, content(a, 1))
  end
  a = evaluate(a, [X, Y-s])
  return a, divexact(map_coeffs(base_ring(a), f, parent =parent(a)), a)
end

end

module MPolyLift
using Hecke

mutable struct LiftPairCtx
  f::MPolyElem{qadic}
  g::MPolyElem{qadic}
  h::MPolyElem{qadic}
  a::MPolyElem{qadic}
  b::MPolyElem{qadic}
  pr_p ::Int
  pr_y ::Int

  function LiftPairCtx(g::MPolyElem{qadic}, h::MPolyElem{qadic})
    return LiftPairCtx(g*h, g, h)
  end
  function LiftPairCtx(f::MPolyElem{qadic}, g::MPolyElem{qadic}, h::MPolyElem{qadic})
    r = new()
    r.f = f
    R = base_ring(g)
    Rx = parent(g)
    k, mk = ResidueField(R)

    gp = map_coeffs(mk, g)
    hp = map_coeffs(mk, h, parent = parent(gp))

    kx, x = PolynomialRing(k, cached = false)
    _, ap, bp = gcdx(evaluate(gp, [x, 0*x]), evaluate(hp, [x, 0*x]))
    r.a = setprecision(map_coeffs(x->preimage(mk, x), ap)(gen(Rx, 1)), 1)
    r.b = setprecision(map_coeffs(x->preimage(mk, x), bp)(gen(Rx, 1)), 1)

    r.g = divrem(setprecision(g, 1), [gen(Rx, 2)])[2]
    r.h = divrem(setprecision(h, 1), [gen(Rx, 2)])[2]

    r.pr_p = r.pr_y = 1

    return r
  end
end

function lift_y(r::LiftPairCtx)
  y = gen(parent(r.f), 2)
  yk = y^r.pr_y
  fgh = divexact(r.f - r.g*r.h, yk)
  G = divrem(fgh*r.b, [r.g])[2]*yk+r.g
  H = divrem(fgh*r.a, [r.h])[2]*yk+r.h
  ogh = divexact(1-r.a*G - r.b*H, yk)
  B = divrem(ogh*r.b, [r.g])[2]*yk+r.b
  A = divrem(ogh*r.a, [r.h])[2]*yk+r.a
  yk *= yk
  r.g = divrem(G, yk)[2]
  r.h = divrem(H, yk)[2]
  r.a = divrem(A, yk)[2]
  r.b = divrem(B, yk)[2]
  r.pr_y *= 2
end

function lift_p(r::LiftPairCtx)
  p = Hecke.uniformizer(base_ring(r.f))
  pk = p^r.pr_p
  r.pr_p *= 2
  y = gen(parent(r.f), 2)
  yk = y^r.pr_y
  r.g = setprecision(r.g, r.pr_p)
  r.h = setprecision(r.h, r.pr_p)
  r.a = setprecision(r.a, r.pr_p)
  r.b = setprecision(r.b, r.pr_p)

  fgh = divexact(r.f - r.g*r.h, pk)
  G = divrem(divrem(divrem(fgh*r.b, [yk])[2], r.g)[2]*pk+r.g, [yk])[2]
  H = divrem(divrem(divrem(fgh*r.a, [yk])[2], r.h)[2]*pk+r.h, [yk])[2]
  ogh = divexact(divrem(1-r.a*G - r.b*H, [yk])[2], pk)
  B = divrem(divrem(ogh*r.b, [yk])[2], r.g)[2]*pk+r.b
  A = divrem(divrem(ogh*r.a, [yk])[2], r.h)[2]*pk+r.a
  r.g = G
  r.h = H
  r.a = A
  r.b = B
end

mutable struct LiftCtx
  f::MPolyElem{qadic}
  C::Array{LiftPairCtx, 1}
  pr_p::Int
  pr_y::Int
  function LiftCtx(f::MPolyElem{qadic}, g::Array{<:MPolyElem{qadic}, 1})
    C = Array{LiftPairCtx, 1}()
    for i=1:div(length(g), 2)
      push!(C, LiftPairCtx(g[2*i-1], g[2*i]))
    end
    i = 1
    if isodd(length(g))
      push!(C, LiftPairCtx(g[end], C[1].f))
      i += 1
    end
    while 2*i <= length(C)
      push!(C, LiftPairCtx(C[2*i-1].f, C[2*i].f))
      i += 1
    end

    r = new()
    r.f = f
    @show C[end].f
    @show f
    @assert degree(C[end].f, 1) == degree(f, 1)
    C[end].f = f
    r.C = C
    r.pr_p = r.pr_y = 1
    return r
  end
end

function lift_y(R::LiftCtx)
  C = R.C
  i = length(C)
  j = i-1

  while i >= 1
    lift_y(C[i])
    if j>= 2
      C[j].f = C[i].h
      C[j-1].f = C[i].g
    elseif j >= 1
      C[j].f = C[i].h
    end
    j -= 2
    i -= 1
  end
  R.pr_y *= 2
end

function lift_p(R::LiftCtx)
  C = R.C
  i = length(C)
  j = i-1

  while i >= 1
    lift_p(C[i])
    if j>= 2
      C[j].f = C[i].h
      C[j-1].f = C[i].g
    elseif j >= 1
      C[j].f = C[i].h
    end
    j -= 2
    i -= 1
  end
  R.pr_p *= 2
end

end

#= revised strategy until I actually understand s.th. better
 - assumming f is monic in y, irreducible over Q and f(0) is squarefree
 - find a prime s.th. f(0, y) is square-free with a small degree splitting field
 - compute roots in F_q[[x]] (finite field!)
 - use combine to find the combinations giving the proper factor
 - find the (small) field the factor lives over
 - lift the factor in the q-adic field
 - find the number field

 TODO:
  bounds on the precisions (poly prec)
  non-monic: use "any_order", the generizlised equation order to produce
    integral power sums -> make_monic
  find good evaluation points
  more variables
  more rings

example:

Qxy, (y, x) = PolynomialRing(QQ, ["y", "x"])
include("/home/fieker/Downloads/n60s3.m");
  #from  https://www.math.univ-toulouse.fr/~cheze/n60s3.m

  r = MPolyFact.roots(P)
  c = MPolyFact.combination(r)
  q = MPolyFact.field(r, c)

  from the Rupprecht paper, but the 3rd is boring (probably wrong)

  y^9+3*x^5*y^6+5*x^4*y^5+3*x^10*y^3-3*x^6*y^3+5*x^9*y^2+x^15

  y^5*x^5+9*x^8*y^4-6*x^14*y^2-18*x^10*y^6-18*x^6*y^10+x^20+5*x^16*y^4+10*x^12*y^8+10*x^8*y^12+5*y^16*x^4+9*y^8*x^4-6*y^14*x^2+y^20

  -125685*x+151959*x^8+917230*x^6+8717398*y^5*x^5+5108544*x^8*y^4-1564434*x^5+7744756*x^5*y^6+306683*x^3*y^6+413268*x^4*y^6+9081976*x^6*y^6+1317780*x^6*y^5+76745*x^4*y^5-15797040*x^7*y^5+99348*x^3*y^5+4106178*x^6*y^4+2010995*x^4*y^4-11264228*x^7*y^4-12465712*x^5*y^4+40908*x^2*y^4+404227*x^3*y^4-9204694*x^7*y^3-49266*x^2*y^3-3500343*x^4*y^3+1512264*x^3*y^3+6405504*x^8*y^3+9879662*x^6*y^3-3821606*x^5*y^3-592704*x^9*y^3-8503779*x^5*y^2-783216*x^9*y^2+10608275*x^6*y^2+574917*x^2*y^2-10143*x*y^2+5943180*x^4*y^2-3295022*x^3*y^2+3452692*x^8*y^2-6432756*x^7*y^2-344988*x^9*y+67473*x*y+2548458*x^4*y-2646351*x^7*y+1059606*x^8*y-3698541*x^5*y-491400*x^2*y+430155*x^3*y+4011984*x^6*y+1530912*x^4+617526
=#


text:
using Test, Random
using SignalBroadcastingSim

@testset "start" begin

Random.seed!(0)

@test rand(Int) == -4635026124992869592
@test 1+2 == 3

end

ext: jl


text:
include("gemm.jl")

using LinearAlgebra

let
    nsamples = 10
    ns = 48:48:1500
    for iter = 1:length(ns)
        n = m = k = ns[iter]
        C = rand(n,m)
        A = rand(n,k)
        B = rand(k,m)

        # julia mul!
        t_mul! = floatmax(Float64)
        for s = 1:nsamples
            t = @elapsed mul!(C,A,B)
            t_mul! = min(t,t_mul!)
        end
        GFLOPS_mul![iter] = 2*n*m*k / t_mul! / 1e9

        t_ijp! = floatmax(Float64)
        for s = 1:nsamples
            t = @elapsed mygemm_ijp!(C,A,B)
            t_ijp! = min(t,t_ijp!)
        end

    end
    df = DataFrame(m = ns,

    )
end


text:
fmatime = matime = 1.0;
for i in 1:5 matime=min(matime, (@timed muladd(3.0,2.0,1.0))[2]) end;
for i in 1:5 fmatime=min(fmatime, (@timed fma(3.0,2.0,1.0))[2]) end;
const HasFMAonChip = fmatime < 2*matime;


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--optimize", "-O"
        arg_type = Int
        metavar = "{0,1,2,3}"
        range_tester = (x->x ∈ (0, 1, 2, 3))
        default = 2
        help = "set the optimization level"
    end
end

function main()
    parsed_args = parse_args(s)
    for (k, v) ∈ parsed_args
        show("$k => $v")
    end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
Round to the nearest integer.
If the number is halfway between two integers,
then round down. (See definition 2.8 (pg 27) of Bremner.)

Motivation: The built-in round function sends -0.5 to 0 instead of -1.
This function correctly handles negative numbers per the examples in Bremner.
"""
function round_ties_down{T<:AbstractFloat}(x::T)
    if abs(floor(x) - x) == 0.5
        round(Int, x, RoundDown)
    else
        round(Int, x)
    end
end

round_ties_down{T<:Integer}(x::T) = x


text:
using .CLocales
using .LC

set_locale!(LocaleId("C"), LC.ALL)
const cl = locale()
const ALLCAT = collect(LC.ALL)

@test cl != Ptr{Nothing}(0)

@testset "clocale identifiers C for $sym"  for sym in ALLCAT
    @test clocale_id(sym) == "C"
end

for cat in ALLCAT
    ENV[string("LC_", cat)] = ""
end

clocid = "fr_FR.utf8"
ENV["LANG"] = clocid
ENV["LC_ALL"] = ""

set_locale!(LocaleId(""), LC.ALL)

@testset "clocale identifiers from envionment for $sym" for sym in ALLCAT
    @test clocale_id(sym) == clocid
end

@testset "clocale codesets from environment for $sym" for sym in ALLCAT
    try
        nli = eval(CLocales, Meta.parse(string(sym, "_CODESET")))
        @test nl_langinfo(nli) == "UTF-8"
    end
end

@testset "setting and querying clocalesi prog" begin
    @test clocale_id(LC.NUMERIC) == "fr_FR.utf8"
    set_locale!(LocaleId("C"))
    set_locale!(LocaleId("tr_TR"), LC.MESSAGES)
    set_locale!(LocaleId("en-us"), LC.MONETARY)
    @test clocale_id(LC.CTYPE) == "C"
    @test clocale_id(LC.NUMERIC) == "C"
    @test clocale_id(LC.MESSAGES) == "tr_TR.utf8"
    @test clocale_id(LC.MONETARY) == "en_US.utf8"
end

ENV["LC_ALL"] = "de_DE.utf8"
set_locale!(LocaleId(""), LC.MESSAGES | LC.NAME)

@testset "setting and querying clocales ENV" begin
    @test clocale_id(LC.NUMERIC) == "C"
    @test clocale_id(LC.MESSAGES) == "de_DE.utf8"
    @test clocale_id(LC.NAME) == "de_DE.utf8"
    @test clocale_id(LC.MONETARY) == "en_US.utf8"
end

import .CLocales:

CTYPE_CODESET,

RADIXCHAR, THOUSEP, THOUSANDS_SEP, GROUPING, NUMERIC_CODESET,

ABDAY, DAY, ABMON, MON, AM_STR, PM_STR, D_T_FMT, D_FMT, T_FMT, T_FMT_AMPM, ERA,
ERA_YEAR, ERA_D_FMT, ALT_DIGITS, ERA_D_T_FMT, ERA_T_FMT, TIME_CODESET,

COLLATE_CODESET,

YESEXPR, NOEXPR, YESSTR, NOSTR, MESSAGES_CODESET,

INT_CURR_SYMBOL, CURRENCY_SYMBOL, MON_DECIMAL_POINT, MON_THOUSANDS_SEP, MON_GROUPING,
POSITIVE_SIGN, NEGATIVE_SIGN, INT_FRAC_DIGITS, FRAC_DIGITS, P_CS_PRECEDES, P_SEP_BY_SPACE,
N_CS_PRECEDES, N_SEP_BY_SPACE, P_SIGN_POSN, N_SIGN_POSN, CRNCYSTR, INT_P_CS_PRECEDES,
INT_P_SEP_BY_SPACE, INT_N_CS_PRECEDES, INT_N_SEP_BY_SPACE, INT_P_SIGN_POSN, INT_N_SIGN_POSN,
MONETARY_CODESET,

PAPER_HEIGHT, PAPER_WIDTH, PAPER_CODESET,

NAME_FMT, NAME_GEN, NAME_MR, NAME_MRS, NAME_MISS, NAME_MS, NAME_CODESET,

ADDRESS_POSTAL_FMT, ADDRESS_COUNTRY_NAME, ADDRESS_COUNTRY_POST, ADDRESS_COUNTRY_AB2,
ADDRESS_COUNTRY_AB3, ADDRESS_COUNTRY_CAR, ADDRESS_COUNTRY_NUM, ADDRESS_COUNTRY_ISBN,
ADDRESS_LANG_NAME, ADDRESS_LANG_AB, ADDRESS_LANG_LIB, ADDRESS_LANG_TERM, ADDRESS_CODESET,

TELEPHONE_TEL_INT_FMT, TELEPHONE_TEL_DOM_FMT, TELEPHONE_INT_SELECT, TELEPHONE_INT_PREFIX,
TELEPHONE_CODESET,

MEASUREMENT, MEASUREMENT_CODESET,

IDENTIFICATION_TITLE, IDENTIFICATION_SOURCE, IDENTIFICATION_ADDRESS, IDENTIFICATION_CONTACT,
IDENTIFICATION_EMAIL, IDENTIFICATION_TEL, IDENTIFICATION_FAX, IDENTIFICATION_LANGUAGE,
IDENTIFICATION_TERRITORY, IDENTIFICATION_AUDIENCE, IDENTIFICATION_APPLICATION, IDENTIFICATION_ABBREVIATION,
IDENTIFICATION_REVISION, IDENTIFICATION_DATE, IDENTIFICATION_CATEGORY, IDENTIFICATION_CODESET

const DE_ITEMS = [
    RADIXCHAR => ",",
    THOUSEP => ".",
    THOUSANDS_SEP => ".",
    GROUPING => 3,

    ABDAY(2) => "Mo",
    DAY(1) => "Sonntag",
    DAY(7) => "Samstag",
    ABMON(6) => "Jun",
    MON(6) => "Juni",
    AM_STR => "",
    PM_STR => "",
    D_T_FMT => "%a %d %b %Y %T %Z",
    D_FMT => "%d.%m.%Y",
    T_FMT => "%T",
    T_FMT_AMPM => "",
    ERA => "",
    ERA_YEAR => "",
    ERA_D_FMT => "",
    ALT_DIGITS => "",
    ERA_D_T_FMT => "",
    ERA_T_FMT => "",

    YESEXPR => "^[jJyY].*",
    NOEXPR => "^[nN].*",

    INT_CURR_SYMBOL => "EUR ",
    CURRENCY_SYMBOL => "€" ,
    MON_DECIMAL_POINT => ",",
    MON_THOUSANDS_SEP => ".",
    MON_GROUPING => 3,
    POSITIVE_SIGN => "",
    NEGATIVE_SIGN => "-",
    INT_FRAC_DIGITS => 2,
    FRAC_DIGITS => 2,
    P_CS_PRECEDES => 0,
    P_SEP_BY_SPACE => 1,
    N_CS_PRECEDES => 0,
    N_SEP_BY_SPACE => 1,
    P_SIGN_POSN => 1 ,
    N_SIGN_POSN => 1,
    CRNCYSTR => "+€",
    INT_P_CS_PRECEDES => 0,
    INT_P_SEP_BY_SPACE => 1,
    INT_N_CS_PRECEDES => 0,
    INT_N_SEP_BY_SPACE => 1,
    INT_P_SIGN_POSN => 1,
    INT_N_SIGN_POSN => 1,

    PAPER_HEIGHT => 297,
    PAPER_WIDTH => 210,

    NAME_FMT => "%d%t%g%t%m%t%f",
    NAME_GEN => "",
    NAME_MR => "Herr",
    NAME_MRS => "Frau",
    NAME_MISS => "Fräulein",
    NAME_MS  => "Frau",

    ADDRESS_POSTAL_FMT => "%f%N%a%N%d%N%b%N%s %h %e %r%N%z %T%N%c%N",
    ADDRESS_COUNTRY_NAME => "Deutschland",
    ADDRESS_COUNTRY_POST => "D",
    ADDRESS_COUNTRY_AB2 => "DE",
    ADDRESS_COUNTRY_AB3 => "DEU",
    ADDRESS_COUNTRY_CAR => "D",
    ADDRESS_COUNTRY_NUM => 276,
    ADDRESS_COUNTRY_ISBN => "3",
    ADDRESS_LANG_NAME => "Deutsch",
    ADDRESS_LANG_AB => "de",
    ADDRESS_LANG_LIB => "deu",
    ADDRESS_LANG_TERM => "ger",

    TELEPHONE_TEL_INT_FMT => "+%c %a %l",
    TELEPHONE_TEL_DOM_FMT => "%A %l",
    TELEPHONE_INT_SELECT => "00",
    TELEPHONE_INT_PREFIX => "49",

    MEASUREMENT => 1,

    IDENTIFICATION_TITLE => "German locale for Germany",
    IDENTIFICATION_SOURCE => "Free Software Foundation, Inc.",
    IDENTIFICATION_ADDRESS => "http://www.gnu.org/software/libc/",
    IDENTIFICATION_CONTACT => "",
    IDENTIFICATION_EMAIL => "bug-glibc-locales@gnu.org",
    IDENTIFICATION_TEL => "",
    IDENTIFICATION_FAX => "",
    IDENTIFICATION_LANGUAGE => "German",
    IDENTIFICATION_TERRITORY => "Germany",
    IDENTIFICATION_AUDIENCE => "",
    IDENTIFICATION_APPLICATION => "",
    IDENTIFICATION_ABBREVIATION => "",
    IDENTIFICATION_REVISION => "1.0",
    IDENTIFICATION_DATE => "2000-06-24",
    IDENTIFICATION_CATEGORY => "de_DE:2000",

]

set_locale!(LocaleId("de-de")) # set all categories
@testset "nl_langinfo queries $(category(p.first))-$(offset(p.first)) => '$(p.second)'" for p in DE_ITEMS
    item, res = p
    @test nl_langinfo(item) == res
end


text:
using Registrator.CommentBot: make_trigger, parse_comment

using Test

@testset "Server" begin

@testset "Trigger comment" begin
    trigger = make_trigger(Dict("trigger" => "@JuliaRegistrator"))
    @test match(trigger, "@JuliaRegistrator hi") !== nothing
    @test match(trigger, "@juliaregistrator hi") !== nothing
    @test match(trigger, "@Zuliaregistrator hi") === nothing
    @test match(trigger, "@zuliaregistrator hi") === nothing
end

@testset "parse_comment" begin
    @test parse_comment("register()") == ("register", Dict{Symbol,String}())
    @test parse_comment("approved()") == ("approved", Dict{Symbol,String}())

    @test parse_comment("register(target=qux)") == ("register", Dict(:target=>"qux"))
    @test parse_comment("register(target=qux, branch=baz)") == ("register", Dict(:branch=>"baz",:target=>"qux"))
    @test parse_comment("register(target=foo, branch=bar)") == ("register", Dict(:branch=>"bar",:target=>"foo"))

    @test parse_comment("register(branch=\"release-1.0\")") == ("register", Dict(:branch=>"release-1.0"))

    @test parse_comment("register") == ("register", Dict{Symbol, String}())
    @test parse_comment("register branch=foo target=bar") == ("register", Dict(:branch => "foo", :target => "bar"))
    @test parse_comment("register branch = foo target = bar") == ("register", Dict(:branch => "foo", :target => "bar"))
    @test parse_comment("register branch=foo, target=bar") == ("register", Dict(:branch => "foo", :target => "bar"))

    @test parse_comment("register(branch=foo)\nfoobar\"") == ("register", Dict(:branch => "foo"))

    @test parse_comment("register branch=foo branch=bar") == (nothing, nothing)
end

end


text:
const v_e = 2550  # rocket exhaust velocity [m/s]
const M = 250000  # rocket dry mass [kg]

function rocket_fuel(v)
    m_fuel = 0

    # Your code goes here!

    return m_fuel
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
"""
A contiguous memory.

This type works like a `Vector` method.
"""
struct Memory
    ptr::Ptr{UInt8}
    size::UInt
end

Memory(data::ByteData) = Memory(pointer(data), sizeof(data))
Base.length(mem::Memory) = mem.size
Base.lastindex(mem::Memory) = Int(mem.size)
Base.checkbounds(mem::Memory, i::Integer) = (1 ≤ i ≤ lastindex(mem)) || throw(BoundsError(mem, i))
Base.getindex(mem::Memory, i::Integer) = (@boundscheck checkbounds(mem, i); unsafe_load(mem.ptr, i))
Base.setindex!(mem::Memory, val::UInt8, i::Integer) = (@boundscheck checkbounds(mem, i); unsafe_store!(mem.ptr, val, i))


text:
"""
    hwe(b::Bgen, v::Variant; T=Float32, decompressed=nothing)
    hwe(p::Preamble, d::Vector{UInt8}, idx::Vector{<:Integer}, layout::UInt8,
        rmask::Union{Nothing, Vector{UInt16}})
Hardy-Weinberg equilibrium test for diploid biallelic case
"""
function hwe(p::Preamble, d::Vector{UInt8}, startidx::Integer, layout::UInt8,
    rmask::Union{Nothing, Vector{UInt16}})
    @assert layout == 2 "hwe only supported for layout 2"
    @assert p.bit_depth == 8 && p.max_probs == 3 && p.max_ploidy == p.min_ploidy
    idx1 = startidx

    # "counts" times 255.
    n00 = 0
    n01 = 0
    n11 = 0
    if p.n_samples >= 16
        @inbounds for n in 1:16:(p.n_samples - p.n_samples % 16)
            idx_base = idx1 + ((n-1) >> 1) << 2
            if rmask !== nothing
                rs = vload(Vec{16,UInt16}, rmask, n)
                if sum(rs) == 0
                    continue
                end
            end
            r = reinterpret(Vec{16, UInt16}, vload(Vec{32, UInt8}, d, idx_base))
            first  = (r & mask_odd)
            second = (r & mask_even) >> 8
            third  = 0x00ff - first - second
            if rmask !== nothing
                first = first * rs
                second = second * rs
                third = third * rs
            end
            n00 += sum(first)
            n01 += sum(second)
            n11 += sum(third)
        end
    end
    rem = p.n_samples % 16
    if rem != 0
        @inbounds for n in ((p.n_samples - rem) + 1) : p.n_samples
            rmask !== nothing && rmask[n] == 0 && continue
            idx_base = idx1 + ((n - 1) << 1)
            n00 += d[idx_base]
            n01 += d[idx_base + 1]
            n11 += 255 - d[idx_base] - d[idx_base + 1]
        end
    end
    @inbounds for n in p.missings
        rmask !== nothing && rmask[n] == 0 && continue
        idx_base = idx1 + ((n - 1) << 1)
        n00 -= d[idx_base]
        n01 -= d[idx_base + 1]
        n11 -= 255 - d[idx_base] - d[idx_base + 1]
    end
    n00 *= one_255th
    n01 *= one_255th
    n11 *= one_255th
    return hwe(n00, n01, n11)
end

@inline ccdf_chisq_1(x) = gamma_inc(convert(typeof(x), 1/2), x/2, 0)[2]
"""
    hwe(n00, n01, n11)
Hardy-Weinberg equilibrium test. `n00`, `n01`, `n11` are counts of homozygotes
and heterozygoes respectively. Output is the p-value of type Float64.
"""
function hwe(n00::Real, n01::Real, n11::Real)
    n = n00 + n01 + n11
    n == 0 && return 1.0
    p0 = (n01 + 2n00) / 2n
    (p0 ≤ 0.0 || p0 ≥ 1.0) && return 1.0
    p1 = 1 - p0
    # Pearson's Chi-squared test
    e00 = n * p0 * p0
    e01 = 2n * p0 * p1
    e11 = n * p1 * p1
    ts = (n00 - e00)^2 / e00 + (n01 - e01)^2 / e01 + (n11 - e11)^2 / e11
    #pval = ccdf(Chisq(1), ts)
    pval = ccdf_chisq_1(ts)

    # TODO Fisher exact test
    return pval
end

"""
    maf(b::Bgen, v::Variant; T=Float32, decompressed=nothing)
    maf(p::Preamble, d::Vector{UInt8}, idx::Vector{<:Integer}, layout::UInt8,
        rmask::Union{Nothing, Vector{UInt16}})
Minor-allele frequency for diploid biallelic case
"""
function maf(p::Preamble, d::Vector{UInt8}, startidx::Integer, layout::UInt8,
    rmask::Union{Nothing, Vector{UInt16}})
    @assert layout == 2 "maf only supported for layout 2"
    @assert p.bit_depth == 8 && p.max_probs == 3 && p.max_ploidy == p.min_ploidy
    idx1 = startidx
    # "counts" times 255.
    dosage_total = 0
    if p.n_samples >= 16
        @inbounds for n in 1:16:(p.n_samples - p.n_samples % 16)
            idx_base = idx1 + ((n-1) >> 1) << 2
            if rmask !== nothing
                rs = vload(Vec{16,UInt16}, rmask, n)
                if sum(rs) == 0
                    continue
                end
            end
            r = reinterpret(Vec{16, UInt16}, vload(Vec{32, UInt8}, d, idx_base))
            first  = (r & mask_odd)  << 1
            second = (r & mask_even) >> 8
            s = first + second
            if rmask !== nothing
                s = s * rs
            end

            dosage_total += sum(s)
        end
    end
    rem = p.n_samples % 16
    if rem != 0
        @inbounds for n in ((p.n_samples - rem) + 1) : p.n_samples
            rmask !== nothing && rmask[n] == 0 && continue
            idx_base = idx1 + ((n - 1) << 1)
            dosage_total += 2 * d[idx_base] + d[idx_base + 1]
        end
    end
    @inbounds for n in p.missings
        rmask !== nothing && rmask[n] == 0 && continue
        idx_base = idx1 + ((n - 1) << 1)
        dosage_total -= 2 * d[idx_base] + d[idx_base + 1]
    end
    dosage_total *= one_255th
    dosage_total /= (rmask !== nothing ? sum(rmask) : p.n_samples - length(p.missings))
    dosage_total < 1.0 ? dosage_total / 2 : 1 - dosage_total / 2
end

"""
    info_score(b::Bgen, v::Variant; T=Float32, decompressed=nothing)
    info_score(p::Preamble, d::Vector{UInt8}, idx::Vector{<:Integer}, layout::UInt8,
        rmask::Union{Nothing, Vector{UInt16}})
Information score of the variant.
"""
function info_score(p::Preamble, d::Vector{UInt8}, startidx::Integer, layout::UInt8,
    rmask::Union{Nothing, Vector{UInt16}})
    @assert layout == 2 "info_score only supported for layout 2"
    @assert p.bit_depth == 8 && p.max_probs == 3 && p.max_ploidy == p.min_ploidy
    @assert length(p.missings) == 0 "current implementation does not allow missingness"
    idx1 = startidx
    # "counts" times 255.
    samples_cum = 0
    mean_cum = 0.0
    sumsq_cum = 0.0
    dosage_sum = 0.0
    dosage_sumsq = 0.0
    rs = Vec{16,UInt16}((1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
    rs_float = Vec{16, Float32}((1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
    rs_float_ones = rs_float
    rs_floatarr = ones(Float32, 16)
    if p.n_samples >= 16
        @inbounds for n in 1:16:(p.n_samples - p.n_samples % 16)
            idx_base = idx1 + ((n-1) >> 1) << 2
            if rmask !== nothing
                rs = vload(Vec{16,UInt16}, rmask, n)
                if sum(rs) == 0
                    continue
                elseif sum(rs) == 16
                    rs_float = rs_float_ones
                else
                    rs_floatarr .= @view(rmask[n:n+15])
                    rs_float = vload(Vec{16, Float32}, rs_floatarr, 1)
                end
            end
            r = reinterpret(Vec{16, UInt16}, vload(Vec{32, UInt8}, d, idx_base))
            first  = (r & mask_odd)  << 1
            second = (r & mask_even) >> 8
            s = first + second
            if rmask !== nothing
                s = s * rs
            end
            dosage_float = one_255th * convert(
                Vec{16, Float32}, s)
            samples_new = sum(rs)
            sum_new = sum(dosage_float)
            samples_prev = samples_cum
            samples_cum = samples_cum + samples_new
            mean_new = sum_new / samples_new
            diff = dosage_float - mean_new
            if rmask !== nothing
                diff *= rs_float
            end
            sumsq_new = sum(diff ^ 2)
            delta = mean_new - mean_cum
            mean_cum += delta * samples_new / samples_cum
            sumsq_cum += sumsq_new + delta ^ 2 * samples_prev * samples_new / samples_cum
        end
    end
    rem = p.n_samples % 16
    if rem != 0
        @inbounds for n in ((p.n_samples - rem) + 1) : p.n_samples
            rmask !== nothing && rmask[n] == 0 && continue
            idx_base = idx1 + ((n - 1) << 1)
            dosage = one_255th * (2 * d[idx_base] + d[idx_base + 1])
            delta = dosage - mean_cum
            samples_cum += 1
            mean_cum += delta / samples_cum
            delta2 = dosage - mean_cum
            sumsq_cum += delta * delta2
        end
    end
    n_samples = rmask !== nothing ? sum(rmask) : p.n_samples
    p = mean_cum / 2
    v = sumsq_cum / (n_samples - 1)
    v / (2p * (1-p))
end

function counts!(p::Preamble, d::Vector{UInt8}, startidx::Integer, layout::UInt8,
    rmask::Union{Nothing, Vector{UInt16}}; r::Union{Nothing,Vector{<:Integer}}=nothing, dosage::Bool=true)
    if dosage
        @assert layout == 2 "hwe only supported for layout 2"
        @assert p.bit_depth == 8 && p.max_probs == 3 && p.max_ploidy == p.min_ploidy
        idx1 = startidx
        if r !== nothing
            @assert length(r) == 512
            fill!(r, 0)
        else
            r = zeros(UInt, 512)
        end
        rs = Vec{16,UInt16}((1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
        if p.n_samples >= 16
            @inbounds for n in 1:16:(p.n_samples - p.n_samples % 16)
                idx_base = idx1 + ((n-1) >> 1) << 2
                if rmask !== nothing
                    rs = vload(Vec{16,UInt16}, rmask, n)
                    if sum(rs) == 0
                        continue
                    end
                end
                q = reinterpret(Vec{16, UInt16}, vload(Vec{32, UInt8}, d, idx_base))
                first  = (q & mask_odd)  << 1
                second = (q & mask_even) >> 8
                s = first + second
                if rmask !== nothing
                    s = s * rs
                end
                @inbounds for i in 1:16
                    if rs[i] != 0
                        ss = s[i]
                        r[ss + 1] += 1
                    end
                end
            end
        end
        rem = p.n_samples % 16
        if rem != 0
            @inbounds for n in ((p.n_samples - rem) + 1) : p.n_samples
                rmask !== nothing && rmask[n] == 0 && continue
                idx_base = idx1 + ((n - 1) << 1)
                r[2 * d[idx_base] + d[idx_base + 1] + 1] += 1
            end
        end
        # subtract back counts for missings
        @inbounds for n in p.missings
            rmask !== nothing && rmask[n] == 0 && continue
            idx_base = idx1 + ((n - 1) << 1)
            r[2 * d[idx_base] + d[idx_base + 1] + 1] -= 1
        end
        r[end] = length(p.missings)
    else
        @error "counts for non-dosage not implemented yet"
    end
    r
end

for ftn in [:maf, :hwe, :info_score, :counts!]
    @eval begin
        function $(ftn)(b::Bgen, v::Variant; T=Float32, decompressed=nothing,
            is_decompressed=false, rmask=nothing, kwargs...)
            io, h = b.io, b.header
            if (decompressed !== nothing && !is_decompressed) ||
                (decompressed === nothing && (v.genotypes === nothing ||
                v.genotypes.decompressed === nothing))
                decompressed = decompress(io, v, h; decompressed=decompressed)
            else
                decompressed = v.genotypes.decompressed
            end
            startidx = 1
            if v.genotypes === nothing
                p = parse_preamble(decompressed, h, v)
                v.genotypes = Genotypes{T}(p, decompressed)
            else
                p = v.genotypes.preamble
            end
            if h.layout == 2
                startidx += 10 + h.n_samples
            end

            $(ftn)(p, decompressed, startidx, h.layout, rmask; kwargs...)
        end
    end
end


text:
@testset "SOCP Atoms: " for solver in solvers
    if can_solve_socp(solver)
        @testset "norm 2 atom" begin
            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            p = Problem{BigFloat}(:minimize, norm2(A * x + b))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.64888 atol=TOL
            @test evaluate(norm2(A * x + b)) ≈ 0.64888 atol=TOL

            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            lambda = 1
            p = Problem{BigFloat}(:minimize, norm2(A * x + b) + lambda * norm2(x), x >= 1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 14.9049 atol=TOL
            @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=TOL

            x = Variable(2)

            p = Problem{BigFloat}(:minimize, norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1)
            @test vexity(p) == ConvexVexity()

            solve!(p, solver)
            @test p.optval ≈ 14.9049 atol=TOL
            @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=TOL

            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            lambda = 1
            p = Problem{BigFloat}(:minimize, norm2(A * x + b) + lambda * norm_1(x), x >= 1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 15.4907 atol=TOL
            @test evaluate(norm2(A * x + b) + lambda * norm_1(x)) ≈ 15.4907 atol=TOL
        end

        @testset "frobenius norm atom" begin
            m = Variable(4, 5)
            c = [m[3, 3] == 4, m >= 1]
            p = Problem{BigFloat}(:minimize, norm(vec(m), 2), c)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ sqrt(35) atol=TOL
            @test evaluate(norm(vec(m), 2)) ≈ sqrt(35) atol=TOL
        end

        @testset "quad over lin atom" begin
            x = Variable(3, 1)
            A = [2 -3 5; -2 9 -3; 5 -8 3]
            b = [-3; 9; 5]
            c = [3 2 4]
            d = -3
            p = Problem{BigFloat}(:minimize, quadoverlin(A*x + b, c*x + d))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 17.7831 atol=TOL
            @test (evaluate(quadoverlin(A * x + b, c * x + d)))[1] ≈ 17.7831 atol=TOL
        end

        @testset "sum squares atom" begin
            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            p = Problem{BigFloat}(:minimize, sumsquares(A*x + b))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.42105 atol=TOL
            @test (evaluate(sumsquares(A * x + b)))[1] ≈ 0.42105 atol=TOL
        end

        @testset "square atom" begin
            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            p = Problem{BigFloat}(:minimize, sum(square(A*x + b)))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.42105 atol=TOL
            @test evaluate(sum(square(A * x + b))) ≈ 0.42105 atol=TOL

            x = Variable(2, 1)
            A = [1 2; 2 1; 3 4]
            b = [2; 3; 4]
            expr = A * x + b
            p = Problem{BigFloat}(:minimize, sum(dot(^)(expr,2))) # elementwise ^
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.42105 atol=TOL
            @test evaluate(sum(broadcast(^, expr, 2))) ≈ 0.42105 atol=TOL

            p = Problem{BigFloat}(:minimize, sum(dot(*)(expr, expr))) # elementwise *
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.42105 atol=TOL
            @test evaluate(sum((dot(*))(expr, expr))) ≈ 0.42105 atol=TOL
        end

        @testset "inv pos atom" begin
            x = Variable(4)
            p = Problem{BigFloat}(:minimize, sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 2 atol=TOL
            @test evaluate(sum(invpos(x))) ≈ 2 atol=TOL

            x = Variable(3)
            p = Problem{BigFloat}(:minimize, sum(dot(/)([3,6,9], x)), x<=3)
            solve!(p, solver)
            @test x.value ≈ fill(3.0, (3, 1)) atol=TOL
            @test p.optval ≈ 6 atol=TOL
            @test evaluate(sum((dot(/))([3, 6, 9], x))) ≈ 6 atol=TOL

            x = Variable()
            p = Problem{BigFloat}(:minimize, sum([3,6,9]/x), x<=3)
            solve!(p, solver)
            @test x.value ≈ 3 atol=TOL
            @test p.optval ≈ 6 atol=TOL
            @test evaluate(sum([3, 6, 9] / x)) ≈ 6 atol=TOL
        end

        @testset "geo mean atom" begin
            x = Variable(2)
            y = Variable(2)
            p = Problem{BigFloat}(:minimize, geomean(x, y), x >= 1, y >= 2)
            # not DCP compliant
            @test vexity(p) == ConcaveVexity()
            p = Problem{BigFloat}(:maximize, geomean(x, y), 1 < x, x < 2, y < 2)
            # Just gave it a vector as an objective, not okay
            @test_throws Exception solve!(p, solver)

            p = Problem{BigFloat}(:maximize, sum(geomean(x, y)), 1 < x, x < 2, y < 2)
            solve!(p, solver)
            @test p.optval ≈ 4 atol=TOL
            @test evaluate(sum(geomean(x, y))) ≈ 4 atol=TOL
        end

        @testset "sqrt atom" begin
            x = Variable()
            p = Problem{BigFloat}(:maximize, sqrt(x), 1 >= x)
        end

        @testset "quad form atom" begin
            x = Variable(3, 1)
            A = [0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
            p = Problem{BigFloat}(:minimize, quadform(x, A), [x >= 1])
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 6.1464 atol=TOL
            @test (evaluate(quadform(x, A)))[1] ≈ 6.1464 atol=TOL

            x = Variable(3, 1)
            A = -1.0*[0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
            c = [3 2 4]
            p = Problem{BigFloat}(:maximize, c*x , [quadform(x, A) >= -1])
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 3.7713 atol=TOL
            @test (evaluate(quadform(x, A)))[1] ≈ -1 atol=TOL
        end

        @testset "huber atom" begin
            x = Variable(3)
            p = Problem{BigFloat}(:minimize, sum(huber(x, 1)), x >= 2)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 9 atol=TOL
            @test evaluate(sum(huber(x, 1))) ≈ 9 atol=TOL
        end

        @testset "rational norm atom" begin
            A = [1 2 3; -1 2 3]
            b = A * ones(3)
            x = Variable(3)
            p = Problem{BigFloat}(:minimize, norm(x, 4.5), [A * x == b])
            @test vexity(p) == ConvexVexity()
            # Solution is approximately x = [1, .93138, 1.04575]
            solve!(p, solver)
            @test p.optval ≈ 1.2717 atol=TOL
            @test evaluate(norm(x, 4.5)) ≈ 1.2717 atol=TOL
        end

        @testset "rational norm dual norm" begin
            v = [0.463339, 0.0216084, -2.07914, 0.99581, 0.889391]
            x = Variable(5)
            q = 1.379;  # q norm constraint that generates many inequalities
            qs = q / (q - 1);  # Conjugate to q
            p = Problem{BigFloat}(:minimize, x' * v)
            p.constraints += (norm(x, q) <= 1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)  # Solution is -norm(v, q / (q - 1))
            @test p.optval ≈ -2.144087 atol=TOL
            @test sum(evaluate(x' * v)) ≈ -2.144087 atol=TOL
            @test evaluate(norm(x, q)) ≈ 1 atol=TOL
            @test sum(evaluate(x' * v)) ≈ -(sum(abs.(v) .^ qs) ^ (1 / qs)) atol=TOL
        end

        @testset "rational norm atom sum" begin
            A = [-0.719255  -0.229089
                -1.33632   -1.37121
                0.703447  -1.4482]
            b = [-1.82041, -1.67516, -0.866884]
            q = 1.5
            xvar = Variable(2)
            p = Problem{BigFloat}(:minimize, .5 * sumsquares(xvar) + norm(A * xvar - b, q))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            # Compute gradient, check it is zero(ish)
            x_opt = xvar.value
            margins = A * x_opt - b
            qs = q / (q - 1);  # Conjugate
            denom = sum(abs.(margins).^q)^(1/qs)
            g = x_opt + A' * (abs.(margins).^(q-1) .* sign.(margins)) / denom
            @test p.optval ≈ 1.7227 atol=TOL
            @test norm(g, 2) ^ 2 ≈ 0 atol=TOL
        end

        @testset "norm consistent with Base for matrix variables" begin
            A = randn(4, 4)
            x = Variable(4, 4)
            x.value = A
            # Matrix norm
            @test evaluate(opnorm(x)) ≈ opnorm(A) atol=TOL
            @test evaluate(opnorm(x, 1)) ≈ opnorm(A, 1) atol=TOL
            @test evaluate(opnorm(x, 2)) ≈ opnorm(A, 2) atol=TOL
            @test evaluate(opnorm(x, Inf)) ≈ opnorm(A, Inf) atol=TOL
            # Vector norm
            # TODO: Once the deprecation for norm on matrices is removed, remove the `vec` calls
            @test evaluate(norm(vec(x), 1)) ≈ norm(vec(A), 1) atol=TOL
            @test evaluate(norm(vec(x), 2)) ≈ norm(vec(A), 2) atol=TOL
            @test evaluate(norm(vec(x), 7)) ≈ norm(vec(A), 7) atol=TOL
            @test evaluate(norm(vec(x), Inf)) ≈ norm(vec(A), Inf) atol=TOL
        end

        @testset "Fixed and freed variables" begin
            @testset "fix and free addition" begin
                x = Variable()
                y = Variable()

                p = Problem{BigFloat}(:minimize, x+y, x>=0, y>=0)
                solve!(p, solver)
                @test p.optval ≈ 0 atol=TOL

                y.value = 4
                fix!(y)
                solve!(p, solver)
                @test p.optval ≈ 4 atol=TOL

                free!(y)
                solve!(p, solver)
                @test p.optval ≈ 0 atol=TOL
            end

            @testset "fix multiplication" begin
                a = [1,2,3,2,1]
                x = Variable(length(a))
                gamma = Variable(Positive())
                fix!(gamma, 0.7)

                p = Problem{BigFloat}(:minimize, norm(x-a) + gamma*norm(x[1:end-1] - x[2:end]))
                solve!(p, solver)
                o1 = p.optval
                # x should be very close to a
                @test o1 ≈ 0.7 * norm(a[1:end - 1] - a[2:end]) atol=TOL
                # increase regularization
                fix!(gamma, 1.0)
                solve!(p, solver)
                o2 = p.optval
                # x should be very close to mean(a)
                @test o2 ≈ norm(a .- mean(a)) atol=TOL

                @test o1 <= o2
            end
        end
    end
end


text:
# pyarg_tuples[i] is a pointer to a python tuple of length i-1
# const pyarg_tuples = Vector{PyPtr}(32)
const pyarg_tuples = PyPtr[]

"""
Low-level version of `pycall!(ret, o, ...; kwargs...)`
Sets `ret.o` to the result of the call, and returns `ret::PyObject`
"""
function _pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, args, kwargs)
    if isempty(kwargs)
        kw = C_NULL
    else
        kw = PyObject(Dict{AbstractString, Any}([Pair(string(k), v) for (k, v) in kwargs]))
    end
    _pycall!(ret, o, args, length(args), kw)
end

"""
Low-level version of `pycall!(ret, o, ...)` for when `kw` is already in python
friendly format but you don't have the python tuple to hold the arguments
(`pyargsptr`). Sets `ret.o` to the result of the call, and returns `ret::PyObject`.
"""
function _pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, args,
  nargs::Int=length(args), kw::Union{Ptr{Cvoid}, PyObject}=C_NULL)
    # pyarg_tuples[i] is a pointer to a python tuple of length i-1
    for n in length(pyarg_tuples):nargs
        push!(pyarg_tuples, @pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,), n))
    end
    check_pyargsptr(nargs)
    pyargsptr = pyarg_tuples[nargs+1]
    return _pycall!(ret, pyargsptr, o, args, nargs, kw) #::PyObject
end

"""
Handle the situation where a callee had previously called incref on the
arguments tuple that was passed to it. We need to hold the only reference to the
arguments tuple, since setting a tuple item is only allowed when there is just
one reference to the tuple (tuples are supposed to be immutable in Python in all
other cases). Note that this actually happens when creating new builtin
exceptions, ref: https://github.com/python/cpython/blob/480ab05d5fee2b8fa161f799af33086a4e68c7dd/Objects/exceptions.c#L48
OTOH this py"def foo(*args): global z; z=args" doesn't trigger this.
Fortunately, this check for ob_refcnt is fast - only a few cpu clock cycles.
"""
function check_pyargsptr(nargs::Int)
    if nargs > 0 && unsafe_load(pyarg_tuples[nargs+1]).ob_refcnt > 1
        pydecref_(pyarg_tuples[nargs+1])
        pyarg_tuples[nargs+1] =
            @pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,), nargs)
    end
end

"""
Low-level version of `pycall!(ret, o, ...)` for when `kw` is already in python
friendly format and you have the python tuple to hold the arguments (`pyargsptr`).
Sets the tuple's values to the python version of your arguments, and calls the
function. Sets `ret.o` to the result of the call, and returns `ret::PyObject`.
"""
function _pycall!(ret::PyObject, pyargsptr::PyPtr, o::Union{PyObject,PyPtr},
  args, nargs::Int=length(args), kw::Union{Ptr{Cvoid}, PyObject}=C_NULL)
    pysetargs!(pyargsptr, args, nargs)
    return __pycall!(ret, pyargsptr, o, kw) #::PyObject
end

"""
```
pysetargs!(pyargsptr::PyPtr, args...)
```
Convert `args` to `PyObject`s, and set them as the elements of the Python tuple
pointed to by `pyargsptr`
"""
function pysetargs!(pyargsptr::PyPtr, args, N::Int)
    for i = 1:N
        pysetarg!(pyargsptr, args[i], i)
    end
end

"""
```
pysetarg!(pyargsptr::PyPtr, arg, i::Integer=1)
```
Convert `arg` to a `PyObject`, and set it as the `i-1`th element of the Python
tuple pointed to by `pyargsptr`
"""
function pysetarg!(pyargsptr::PyPtr, arg, i::Integer=1)
    pyarg = PyObject(arg)
    pyincref(pyarg) # PyTuple_SetItem steals the reference
    @pycheckz ccall((@pysym :PyTuple_SetItem), Cint,
                     (PyPtr,Int,PyPtr), pyargsptr, i-1, pyarg)
end

"""
Lowest level version of  `pycall!(ret, o, ...)`, assumes `pyargsptr` and `kw`
have all their args set to Python values, so we can just call the function `o`.
Sets `ret.o` to the result of the call, and returns `ret::PyObject`.
"""
function __pycall!(ret::PyObject, pyargsptr::PyPtr, o::Union{PyObject,PyPtr},
  kw::Union{Ptr{Cvoid}, PyObject})
    sigatomic_begin()
    try
        retptr = @pycheckn ccall((@pysym :PyObject_Call), PyPtr, (PyPtr,PyPtr,PyPtr), o,
                        pyargsptr, kw)
        pyincref_(retptr)
        pydecref(ret)
        ret.o = retptr
    finally
        sigatomic_end()
    end
    return ret #::PyObject
end

"""
```
pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, returntype::Type, args...; kwargs...)
```
Set `ret` to the result of `pycall(o, returntype, args...; kwargs)` and return
`convert(returntype, ret)`.
Avoids allocating an extra PyObject for `ret`. See `pycall` for other details.
"""
pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...) =
    return convert(returntype, _pycall!(ret, o, args, kwargs))

function pycall!(ret::PyObject, o::T, returntype::Type{PyObject}, args...; kwargs...) where
  T<:Union{PyObject,PyPtr}
    return _pycall!(ret, o, args, kwargs)
end

pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, ::Type{PyAny}, args...; kwargs...) =
    return convert(PyAny, _pycall!(ret, o, args, kwargs))

"""
```
pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...)
```
Call the given Python function (typically looked up from a module) with the
given args... (of standard Julia types which are converted automatically to the
corresponding Python types if possible), converting the return value to
returntype (use a returntype of PyObject to return the unconverted Python object
reference, or of PyAny to request an automated conversion)
"""
pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...) =
    return convert(returntype, _pycall!(PyNULL(), o, args, kwargs))::returntype

pycall(o::Union{PyObject,PyPtr}, ::Type{PyAny}, args...; kwargs...) =
    return convert(PyAny, _pycall!(PyNULL(), o, args, kwargs))

(o::PyObject)(args...; kwargs...) =
    return convert(PyAny, _pycall!(PyNULL(), o, args, kwargs))

PyAny(o::PyObject) = convert(PyAny, o)

"""
    @pycall func(args...)::T

Convenience macro which turns `func(args...)::T` into pycall(func, T, args...)
"""
macro pycall(ex)
    if !(isexpr(ex,:(::)) && isexpr(ex.args[1],:call))
        throw(ArgumentError("Usage: @pycall func(args...)::T"))
    end
    func = ex.args[1].args[1]
    args, kwargs = ex.args[1].args[2:end], []
    if isexpr(args[1],:parameters)
        kwargs, args = args[1], args[2:end]
    end
    T = ex.args[2]
    :(pycall($(map(esc,[kwargs; func; T; args])...)))
end


text:
using Documenter
using LogNormals
using Distributions, StatsBase, StatsPlots, Missings, MissingStrategies

push!(LOAD_PATH,"../src/")
# need to add Statistics and Distributions to Project.toml in docs/
#DocMeta.setdocmeta!(LogNormals, :DocTestSetup, :(using Statistics,Distributions,LogNormals); recursive=true)
#DocMeta.setdocmeta!(Distributions, :DocTestSetup, :(using Statistics,Distributions,LogNormals); recursive=true)
makedocs(sitename="LogNormals.jl",
         doctest=false,
         pages = [
            "Home" => "index.md",
            "LogNormal properties" => "lognormalprops.md",
            "Fit to statistic" => "fitstats.md",
            "sem correlated" => "semcor.md",
         ],
         modules = [LogNormals],
         format = Documenter.HTML(prettyurls = false)
)
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/bgctw/LogNormals.jl.git",
    devbranch = "main"
)


text:
using Dates
import Hyperscript as HS

node = HS.m

# March 5, 2019
dfmt(d) = Dates.format(d, "U d, yyyy")

function hfun_dfmt(p::Vector{String})
    d = getlvar(Symbol(p[1]))
    return dfmt(d)
end

# ===============================================
# Logic to show the list of tags for a page
# ===============================================

function hfun_page_tags()
    tags = get_page_tags()
    base = globvar(:tags_prefix)
    return join(
        (
            node("span", class="tag",
                node("a", href="/$base/$id/", name)
            )
            for (id, name) in tags
        ),
        node("span", class="separator", "•")
    )
end

# ===============================================
# Logic to show the list of tags for a page
# ===============================================

function hfun_taglist()
    return hfun_list_posts(getlvar(:tag_name))
end

# ===============================================
# Logic to retrieve posts in posts/ and display
# them as a list sorted by anti-chronological
# order.
#
# Assumes that 'date' and 'title' are defined for
# all posts.
# ===============================================

function hfun_list_posts(t::String)
    return string(
        node("ul",
                (
                    node("li",
                        node("span", class="date", Dates.format(p.date, "U d, yyyy") * "  — "),
                        node("a", class="title", href=p.href, p.title)
                    )
                    for p in get_posts(t)
                )...
            )
        )
end
hfun_list_posts() = hfun_list_posts("")

function get_posts(t::String)
    # find all valid "posts/xxx.md" files, exclude the index which is where
    # the post-list gets placed
    paths = joinpath.(
        "posts",
        filter!(
            p -> endswith(p, ".md") && p != "index.md",
            readdir("posts")
        )
    )
    # for each of those posts, retrieve date and title, both are expected
    # to be there
    posts = [
        (
            date  = getvarfrom(:date, rp),
            title = getvarfrom(:title, rp),
            href  = "/$(splitext(rp)[1])",
            tags  = get_page_tags(rp)
        )
        for rp in paths
    ]
    sort!(posts, by = x -> x.date, rev=true)
    if !isempty(t)
        filter!(
            p -> t in values(p.tags),
            posts
        )
    end
    return posts
end


text:
using RobRecSolver
using Base.Test

tests = [
    "initial_scenario_tests",
    "minimum_knapsack_problem_tests",
    "minimum_assignment_problem_tests",
    "incremental_problem_tests",
    "recoverable_problem_tests",
    "evaluation_problem_tests",
    "adversarial_problem_tests",
    "selection_lower_bound_tests",
    "lagrangian_lower_bound_tests"
    ]

@time for test in tests
    println("Starting $(test)...")
    include("$(test).jl")
    println("$(test) finished.")
end


text:
################################################################################
## Kullback-Leibler - KL
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate(div::KL, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("First array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    for i = eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        if ai < 0 || bi < 0
            r = Inf
            break
        end
        ui = ai / bi
        r += ai * log(ui) - ai + bi
    end
    return r
end

function evaluate(div::KL, a::AbstractVector{T}) where T <: AbstractFloat
    r = zero(T)
    onet = one(T)
    for i = eachindex(a)
        @inbounds ai = a[i]
        r += xlogx(ai) - ai + onet
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient(div::KL, a::T, b::T) where T <: AbstractFloat
    if b <= 0
        u = convert(T, -Inf)
    elseif a > 0 && b > 0
        u = log(a/b)
    else
        u = convert(T, -Inf)
    end
    return u
end

function gradient(div::KL, a::T) where T <: AbstractFloat
    if a > 0
        u = log(a)
    else
        u = convert(T, -Inf)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian(div::KL, a::T, b::T) where T <: AbstractFloat
    if b==0
        u = convert(T, Inf)
    elseif a > 0 && b > 0
        u = 1.0/a
    else
        u = convert(T, Inf)
    end
    return u
end

function hessian(div::KL, a::T) where T <: AbstractFloat
    if a > 0
        u = 1.0/a
    else
        u = convert(T, Inf)
    end
    return u
end

################################################################################
## Modified Kullback-Leibler - MKL
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate(div::MKL, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    ϑ   = div.ϑ
    f0, f1, f2, uϑ = div.m
    @inbounds for i = eachindex(a,b)
        ai = a[i]
        bi = a[i]
        ui = ai/bi
        if ui >= uϑ
            r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui > 0 && ui < uϑ
            r += (ai*log(ui) - ai + bi)
        else
            r = oftype(a, Inf)
        end
    end
    return r
end

function evaluate(div::MKL, a::AbstractVector{T}) where T <: AbstractFloat
    r = zero(T)
    f0, f1, f2, uϑ = div.m
    @inbounds for i = eachindex(a)
        ui = a[i]
        if ui >= uϑ
            r += f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2
        elseif ui>0 && ui<uϑ
            r += ui*log(ui) - ui + 1.0
        else
            r = oftype(ui, Inf)
            break
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient(div::MKL, a::T) where T <: AbstractFloat
    f0, f1, f2, uϑ = div.m
    if a >= uϑ
        u =  f1 + f2*(a-uϑ)
    else
        u = gradient(div.d, a)
    end
    return u
end

function gradient(div::MKL, a::T, b::T) where T <: AbstractFloat
    f0, f1, f2, uϑ = div.m
    ui = a/b
    if ui > uϑ
        u = f1 + f2*(ui-uϑ)
    else
        u = gradient(div.d, a, b)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian(div::MKL, a::T) where T <: AbstractFloat
    f0, f1, f2, uϑ = div.m
    if a >= uϑ
        u  = f2
    else
        u = hessian(div.d, a)
    end
    return u
end

function hessian(div::MKL, a::T, b::T) where T <: AbstractFloat
    f0, f1, f2, uϑ = div.m
    ui = a/b
    if ui >= uϑ
        u = f2/b
    else
        u = hessian(div.d, a, b)
    end
    return u
end

################################################################################
## Fully Modified Kullback-Leibler - FMKL
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate(div::FMKL, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("First array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m

    @inbounds for i = eachindex(a,b)
        ai = a[i]
        bi = a[i]
        ui = ai/bi

        if ui >= uϑ
            r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui <= uφ
            r += (g0 + g1*(ui-uφ) + .5*g2*(ui-uφ)^2)*bi
        else
            r += (ai*log(ui) - ai + bi)
        end
    end
    return r
end

function evaluate(div::FMKL, a::AbstractVector{T}) where T <: AbstractFloat
    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    @inbounds for i = eachindex(a,b)
        ai = a[i]
        if ai >= uϑ
            r += (f0 + f1*(ai-uϑ) + .5*f2*(ai-uϑ)^2)*bi
        elseif ai <= uφ
            r += (g0 + g1*(ai-uφ) + .5*g2*(ai-uφ)^2)*bi
        else
            r += (ai*log(ai) - ai + 1.0)
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient(div::FMKL, a::T) where T <: AbstractFloat
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u =  f1 + f2*(a-uϑ)
    elseif a <= uφ
        u =  g1 + g2*(a-uφ)
    else
        u = gradient(div.d, a)
    end
    return u
end

function gradient(div::FMKL, a::T, b::T) where T <: AbstractFloat
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    ui = a/b
    if ui > uϑ
        u = f1 + f2*(ui-uϑ)
    elseif ui <= uφ
        u = g1 + g2*(ui-uφ)
    else
        u = gradient(div.d, a, b)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian(div::FMKL, a::T) where T <: AbstractFloat
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u  = f2
    elseif a <= uφ
        u = g2
    else
        u = hessian(div.d, a)
    end
    return u
end

function hessian(div::FMKL, a::T, b::T) where T <: AbstractFloat
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    ui = a/b
    if ui >= uϑ
        u = f2/b
    elseif ui <= uφ
        u = g2/b
    else
        u = hessian(div.d, a, b)
    end
    return u
end


text:
include("differencial.jl")
"""
`RK4(F,A,dt)`
Runge-Kutta法で刻み時間`dt`秒後のFの値を計算する関数
## Arguments
- `F::Any`: 計算したい時間`t`の関数の`t`秒時点の値 関数`dif`で微分値を求められるもの
- `A::Any`:微分計算に必要な定数の組
- `dt::Number`  :刻み時間
## Returns
- `dt`秒後のFの値
"""
function RK4(F,A,dt)
    f=dif(F,A)
    k1=dt.*f
    k2=dt.*dif((F.+k1),A)
    k3=dt.*dif((F.+0.5.*(k1.+k2)),A)
    k4=dt.*dif((F.+k1.+k3),A)
    return F.+(k1.+2.0.*k2.+2.0.*k3.+k4)./6
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The clnlbeam problem

# Based on an AMPL model by Hande Y. Benson
#
# Copyright (C) 2001 Princeton University
# All Rights Reserved

# Permission to use, copy, modify, and distribute this software and its
# documentation for any purpose and without fee is hereby granted, provided that
# the above copyright notice appear in all copies and that the copyright notice
# and this permission notice appear in all supporting documentation.

# Source: H. Maurer and H.D. Mittelman, "The non-linear beam via optimal control
# with bound state variables", Optimal Control Applications and Methods 12, pp.
# 19-31, 1991.

using JuMP
import Ipopt
import Test  #src

function example_clnlbeam()
    N = 1000
    h = 1/N
    alpha = 350
    model = Model(Ipopt.Optimizer)
    @variables(model, begin
           -1 <= t[1:(N + 1)] <= 1
        -0.05 <= x[1:(N + 1)] <= 0.05
                 u[1:(N + 1)]
    end)
    @NLobjective(
        model,
        Min,
        sum(
            0.5 * h * (u[i + 1]^2 + u[i]^2) +
            0.5 * alpha * h * (cos(t[i + 1]) + cos(t[i]))
            for i = 1:N
        ),
    )
    @NLconstraint(
        model,
        [i = 1:N],
        x[i + 1] - x[i] - 0.5 * h * (sin(t[i + 1]) + sin(t[i])) == 0,
    )
    @constraint(
        model,
        [i = 1:N],
        t[i + 1] - t[i] - 0.5 * h * u[i + 1] - 0.5 * h * u[i] == 0,
    )
    optimize!(model)
    println("""
    termination_status = $(termination_status(model))
    primal_status      = $(primal_status(model))
    objective_value    = $(objective_value(model))
    """)
    Test.@test termination_status(model) == MOI.LOCALLY_SOLVED  #src
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT  #src
    Test.@test objective_value(model) ≈ 350.0  #src
    return
end

example_clnlbeam()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distributions
using Test

# #1474
@testset "median and quantile" begin
    for (a, b) in ((-8, -6), (-8, -5), (-5, 1), (-5, 2), (-2, 3), (-2, 4), (2, 4), (2, 5))
        d = DiscreteUniform(a, b)

        for x in a:b
            @test @inferred(quantile(d, cdf(d, x))) === x
        end

        @test @inferred(median(d)) === quantile(d, 1//2)
    end
end

@testset "type stability regression tests" begin
    @test @inferred(mgf(DiscreteUniform(2, 5), 0//1)) === 1.0
    @test @inferred(cf(DiscreteUniform(1, 3), 0//1)) === 1.0 + 0.0im
end


text:
#=
We need to track the current OpenGL context.
Since we can't do this via pointer identity  (OpenGL may reuse the same pointers)
We go for this slightly ugly version.
=#
const ACTIVE_OPENGL_CONTEXT = Base.RefValue{Any}(nothing)

function current_context()
    isassigned(ACTIVE_OPENGL_CONTEXT) || error("No active context")
    ctx = ACTIVE_OPENGL_CONTEXT[]
    ctx === nothing && error("No active context")
    return ctx
end

function is_current_context(x)
    return x == ACTIVE_OPENGL_CONTEXT[]
end

function native_context_alive(x)
    error("Not implemented for $(typeof(x))")
end

"""
Is current context & is alive
"""
is_context_active(x) = is_current_context(x) && context_alive(x)

"""
Has context been destroyed or is it still living?
"""
context_alive(x) = native_context_alive(x)

function native_switch_context!(x)
    error("Not implemented for $(typeof(x))")
end

"""
Invalidates the current context
"""
function switch_context!()
    # for reverting to no context
    ACTIVE_OPENGL_CONTEXT[] = nothing
end

"""
Switches to a new context `x`. Is a noop if `x` is already current
"""
function switch_context!(x)
    if !is_current_context(x)
        ACTIVE_OPENGL_CONTEXT[] = x
        native_switch_context!(x)
    end
end


text:
using GaussQuadrature, GLMakie

function cubed_sphere_warp(
    a,
    b,
    c,
    R = max(abs(a), abs(b), abs(c)),
)

    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        ζ1 = sR / sqrt(X^2 + Y^2 + 1)
        ζ2, ζ3 = X * ζ1, Y * ζ1
        ζ1, ζ2, ζ3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : formulas for Face I from Ronchi, Iacono, Paolucci (1996)
        #              but for us face II of the developed net of the cube
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : formulas for Face II from Ronchi, Iacono, Paolucci (1996)
        #              but for us face III of the developed net of the cube
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : formulas for Face III from Ronchi, Iacono, Paolucci (1996)
        #              but for us face IV of the developed net of the cube
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : formulas for Face IV from Ronchi, Iacono, Paolucci (1996)
        #              but for us face I of the developed net of the cube
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : formulas for Face V from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : formulas for Face VI from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubed_sphere_warp(::EquiangularCubedSphere): $a, $b, $c")
    end

    return x1, x2, x3
end

N = 5
ξ¹, ω¹ = GaussQuadrature.legendre(N, GaussQuadrature.both)
ξ² = copy(ξ¹)
ξ³ = copy(ξ¹)

Ne = 12
vertices = collect(range(-1, 1, length = Ne + 1))
x¹ = zeros(N, Ne)
x² = copy(x¹)
x³ = copy(x¹)

for e in 1:Ne
    @. x¹[:, e] = (ξ¹ + 1) / 2 * (vertices[e+1] - vertices[e]) + vertices[e]
    @. x²[:, e] = (ξ¹ + 1) / 2 * (vertices[e+1] - vertices[e]) + vertices[e]
    @. x³[:, e] = (ξ¹ + 1) / 2 * (vertices[e+1] - vertices[e]) + vertices[e]
end

Nex = Ney = Ne
face1 = [(x¹[1], x²[j, ex], x³[k, ey]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface1 = [cubed_sphere_warp(face1[i]...) for i in eachindex(face1)]

face2 = [(x¹[end], x²[j, ex], x³[k, ey]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface2 = [cubed_sphere_warp(face2[i]...) for i in eachindex(face2)]

face3 = [(x¹[j, ey], x²[1], x³[k, ex]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface3 = [cubed_sphere_warp(face3[i]...) for i in eachindex(face1)]

face4 = [(x¹[j, ey], x²[end], x³[k, ex]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface4 = [cubed_sphere_warp(face4[i]...) for i in eachindex(face1)]

face5 = [(x¹[j, ex], x²[k, ey], x³[1]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface5 = [cubed_sphere_warp(face5[i]...) for i in eachindex(face1)]

face6 = [(x¹[j, ex], x²[k, ey], x³[end]) for j in 1:N, k in 1:N, ex in 1:Nex, ey in 1:Ney]
warpedface6 = [cubed_sphere_warp(face6[i]...) for i in eachindex(face1)]

# Now Plot
lw = 1  # linewidth
ms = 25 # marker size

fig = Figure(resolution = (1400, 1200))
ax = LScene(fig[1, 1])
# fig, ax, sc = scatter(warpedface1, color = :red, markersize = 0.0, show_axis = false)
# scatter(warpedface1, color = :red, markersize = 0.0, show_axis = false)
r = range(-1, 1, length = Ne + 1)

xa = zeros(Ne + 1, Ne + 1)
xb = zeros(Ne + 1, Ne + 1)
xc = zeros(Ne + 1, Ne + 1)

xa .= r
xb .= r'
xc .= 1

a = [cubed_sphere_warp(a, b, c)[1] for (a, b, c) in zip(xa, xb, xc)] .* 1.00
b = [cubed_sphere_warp(a, b, c)[2] for (a, b, c) in zip(xa, xb, xc)] .* 1.00
c = [cubed_sphere_warp(a, b, c)[3] for (a, b, c) in zip(xa, xb, xc)] .* 1.00

wireframe!(ax,
    a, b, c,
    show_axis = false,
    linewidth = lw)

wireframe!(ax,
    -c, a, b,
    show_axis = false,
    linewidth = lw * 5)

scatter!(ax, warpedface1, color = :red, markersize = ms * 0.75,)
scatter!(ax, warpedface3, color = :blue, markersize = ms,)

rotate_cam!(fig.scene.children[1], (π/12, π, 0))

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
@testset "partition" begin
    for n in 1:17
        for p in 1:n
            Is = [partition(n, p, i) for i in 1:p]
            @test vcat(Is...) == collect(1:n)
            @test maximum(length, Is) <= minimum(length, Is) + 1
        end
    end
end

@testset "translate_partition" begin
    n = 17
    for p in 1:n
        for i in 1:p
            for q in 1:n
                j = translate_partition(n, p, q, i)
                @test first(partition(n, p, i)) in partition(n, q, j)
            end
        end
    end
end

@testset "align_partitions" begin
    n = 10
    p, q, i = 4, 5, 1
    @test align_partitions(n, p, q, i) == 1
    p, q, i = 4, 5, 2
    @test align_partitions(n, p, q, i) == 2
    p, q, i = 4, 5, 3
    @test align_partitions(n, p, q, i) == 2
    p, q, i = 4, 5, 4
    @test align_partitions(n, p, q, i) == 2
    p, q, i = 4, 2, 2
    @test align_partitions(n, p, q, i) == 1
    p, q, i = 4, 2, 4
    @test align_partitions(n, p, q, i) == 2

    n = 17
    for p in 1:n
        for i in 1:p
            for q in 1:n
                j = align_partitions(n, p, q, i)
                ip = translate_partition(n, q, p, j)
                @test first(partition(n, p, ip)) == first(partition(n, q, j))
            end
        end
    end
end


text:
include("rotate.jl")

include("clockwise.jl")
include("counter.jl")


text:
module Example

using Twinkle, Twinkle.FlexUI
using Colors
using StaticArrays
using GeometryBasics
import Luxor

println("Start [$(splitext(basename(@__FILE__))[1])]")

app = App()
prop!(app, :title, "Twinkle Example - Luxor")

#---------------------------------------------------------------
# Define Variables
#---------------------------------------------------------------

angle = addVariable!(app, Variable(name="angle", type="flota64",value=8))
count = addVariable!(app, Variable(name="count", type="flota64",value=16))

image = addVariable!(app, Variable(
    name="image",
    type="image",
    value="",
))

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------

ui = VContainer(
    H1Label("Value Controls"),
    Slider(
        text="Angle (π over value)",
        trailing_text="[\$()]",
        min=2,
        max=16,
        value=8,
        variable="angle"
    ),
    Slider(
        text="Count",
        trailing_text="[\$()]",
        min=1,
        max=30,
        value=16,
        variable="count"
    ),
    H1Label("Static Image Example"),
    Image(
        source="\$(image)",
        # height="50%",
    ),
    H1Label("Image Viewer allowing Pan (left-drag) and Zoom (wheel)"),
    PanZoom(
        style="width: 100%; height=400px;",
        content = Container(
            direction = "row warp",
            children = [
                Image(
                    source="\$(image)",
                ),
            ]
        ),
    ),

)
FlexUI.controls!(app, ui)

#---------------------------------------------------------------
# the render function
#---------------------------------------------------------------
function render()
    Luxor.Drawing(600, 400, :svg)
    Luxor.origin()
    Luxor.background("white")

    for θ in range(0, step=π/angle[], length=count[])
        Luxor.gsave()
        Luxor.scale(0.2)
        Luxor.rotate(θ)
        Luxor.translate(350, 0)
        Luxor.julialogo(action=:fill, bodycolor=Luxor.randomhue())
        Luxor.grestore()
    end

    Luxor.gsave()
    Luxor.scale(0.3)
    Luxor.juliacircles()
    Luxor.grestore()

    Luxor.translate(150, -150)
    Luxor.scale(0.3)
    Luxor.julialogo()
    Luxor.finish()

    image["svg"] = Luxor.preview()
end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
using LogDensityProblems, DynamicHMC
using Distributions, Parameters, Random
using PositiveFactorizations, Calculus
using Flux

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/gali_bayesian.jl"))
include(string(pwd(),"/mcmc/priors_v2.jl"))
include(string(pwd(),"/mcmc/foos.jl"))

yy,shocks = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,500)

prob = DSGE_Model(yy[:,2],1/3)

prob((bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

t = problem_transform(prob)

P = TransformedLogDensity(t,prob)

aa = [0.99,6,2/3,1,1,1,1.5,0.5/4,0.5]

true_vals = TransformVariables.inverse(t,(bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

LogDensityProblems.logdensity(P,true_vals)

hes = -Calculus.hessian(x->LogDensityProblems.logdensity(P,x),true_vals)

hes_inv = inv(hes)

hes_inv = cholesky(Positive,hes_inv)

hes_inv = hes_inv.L*hes_inv.L'

isposdef(hes_inv)

dist = MvNormal(hes_inv)

N = 10000

dados = zeros(Float32,10,N)

i = 1

while i <= N
    x = rand(dist)
    y = LogDensityProblems.logdensity(P,x)
    if y < -90000
        continue
    else
        dados[1:9,i] = x
        dados[10,i] = y
        global i += 1
    end
    println(i)
end

dados_orig = dados
dados[10,:] = log.(-dados[10,:])
dados = Flux.normalise(dados)

model_flux = Chain(
    Dense(9,20,relu),
    Dense(20,50,relu),
    Dense(50,50,relu),
    Dense(50,1)
)

function make_minibatch(X, Y, siz)
    indexes = sample(1:size(Y,1),siz)
    X_batch = X[:,indexes]
    Y_batch = Y[indexes]
    return (X_batch, Y_batch), indexes
end

#mini_batch_test = [make_minibatch(dados_x,dados_y,100) for i in 1:100]

loss_flux(x,y) = Flux.mse(model_flux(x),y)

optim_flux = Flux.RMSProp(0.00001)

epochs_max = 100
batch_size = 625
patience = 5

ofs_prop = 0.2

ofs_size = Int(ceil(batch_size*ofs_prop))

acc(x,y) = mean([(model_flux(x)[1,i]-y[i])^2 for i in 1:length(y)])

final_countdown = 0
epoch = 2
last_update = 0

while epoch <= epochs_max&&final_countdown <= patience
    batch,id = make_minibatch(dados[1:9,:],dados[10,:],batch_size)
    batch = [batch for i in 1:batch_size]

    train_batch = batch[1:(ofs_size)]

    train_batch_ids = id[1:(ofs_size)]
    val_batch_ids = id[(ofs_size+1):batch_size]

    loss_in_bef = acc(dados[1:9,train_batch_ids],dados[10,train_batch_ids])*1000
    loss_val_bef = acc(dados[1:9,val_batch_ids],dados[10,val_batch_ids])*1000

    Flux.train!(loss_flux,Flux.params(model_flux),train_batch,optim_flux)

    loss_in_aft = acc(dados[1:9,train_batch_ids],dados[10,train_batch_ids])*1000
    loss_val_aft = acc(dados[1:9,val_batch_ids],dados[10,val_batch_ids])*1000

    delta_loss_val = loss_val_aft - loss_val_bef
    delta_loss_in = loss_in_aft - loss_in_bef

    println("Epoch ", epoch, " ΔLoss out-of-sample *1000 ", delta_loss_val," ΔLoss in-sample *1000 ", delta_loss_in)

    if delta_loss_val > 1e-4
        if last_update != epoch-1
            global final_countdown = 0
        end
        global final_countdown += 1
        global last_update = epoch
    end
        global epoch += 1
end

idxs = sample(1:size(dados,2),500)

scatter(reshape(model_flux(dados[1:9,idxs]),500), label = "Fit")
scatter!(dados[10,idxs], label = "Data")

scatter(reshape(model_flux(dados[1:9,idxs]),500),dados[10,idxs])

acc(dados[1:9,idxs],dados[10,idxs])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
"""

$(TYPEDEF)

Structure that contains some data required to compute the linked cells.

$(TYPEDFIELDS)

"""
struct Box{FloatVector,IntVector}
    sides::FloatVector
    nc::IntVector
    l::FloatVector
    lcell::Int
end

function Box(lcell::Int, sides::AbstractVector, cutoff::Float64)
    # Compute the number of cells in each dimension
    nc = SVector{3,Int}(max.(1, trunc.(Int, sides / (cutoff / lcell))))
    l = SVector{3,Float64}(sides ./ nc)
    Box(sides, nc, l, lcell)
end


text:
module MinComposite

using LinearAlgebra
using OptimPackNextGen
using Printf

import ..Utils: affine_model!

export ProblemParams, SubproblemOptParams, value_function!, min_composite

struct ProblemParams
    c::Function
    ∇c_adjoint::Function
    gauge_projection::Function
    dual_gauge::Union{Function, Nothing}
end

mutable struct IntermediateVars{T}
    r::Vector{T}
    z::Vector{T}
    tau::Real
end

struct SubproblemOptParams{T}
    mem::Integer
    lower::Union{T,Vector{T},Nothing}
    upper::Union{T,Vector{T},Nothing}
    ftol::Tuple{Real,Real}
    maxiter::Integer
    maxeval::Integer
    function SubproblemOptParams{T}(mem::Integer=10,
                                    lower::Union{T,Vector{T},Nothing}=nothing,
                                    upper::Union{T,Vector{T},Nothing}=nothing,
                                    ftol::Tuple{Real,Real}=(0.0,1e-3),
                                    maxiter::Integer=20,
                                    maxeval::Integer=100) where {T}
        mem > 0 || error("mem must be > 0, got $(mem)")
        ftol[1] ≥ 0.0 || error("ftol[1] must be >= 0, got $(ftol[1])")
        ftol[2] ≥ 0.0 || error("ftol[2] must be >= 0, got $(ftol[2])")
        maxiter > 0 || error("maxiter must be > 0, got $(maxiter)")
        maxeval > 0 || error("maxeval must be > 0, got $(maxeval)")
        new(mem,lower,upper,ftol,maxiter,maxeval)
    end
end

function value_function_objective!(x, df, params::ProblemParams, info::IntermediateVars{T})::Float64 where {T}
    params.c(info.r, x)
    params.gauge_projection(info.z, info.r, info.tau)
    info.r .= info.r .- info.z
    params.∇c_adjoint(df, x, info.r)
    return 0.5*norm(info.r)^2
end

function value_function!(x::Vector{T}, params::ProblemParams, info::IntermediateVars, opt_params::SubproblemOptParams{T}) where {T}
    obj! = (x,g)->value_function_objective!(x,g,params,info)
    x .= vmlmb(obj!, x, mem=opt_params.mem, ftol=opt_params.ftol, maxiter=opt_params.maxiter, maxeval=opt_params.maxeval)

    f = obj!(x,zeros(size(x)))
    if params.dual_gauge != nothing
        df = -params.dual_gauge(info.r)::T
    else
        df = convert(0.0, T)
    end
    return f, df
end

"""
min_composite(A, b, gauge_projection, x0, dims, opt_params, dual_gauge)

Value function approach to solve problems of the form
    min_{x} h(A*x - b)

for a non-smooth, often convex, function h(z) with easy to project on level sets.
"""
function min_composite(A::Matrix{T},
                       b::Vector{T},
                       gauge_projection::Function,
                       opt_params::SubproblemOptParams{T},
                       dual_gauge::Union{Function, Nothing}=nothing;
                       num_outer_iter::Integer = 20,
                       x0::Union{Vector{T},Nothing}=nothing,
                       initial_tau::Real=0.0,
                       use_newton::Bool=true,
                       verbose::Bool=false) where {T}
    c = (y,x)->affine_model!(y,A,x,b)
    ∇c_adjoint = (y,x,r)->mul!(y,A',r)
    return min_composite(c, ∇c_adjoint, gauge_projection, size(A), opt_params, dual_gauge,
                        num_outer_iter=num_outer_iter, x0=x0, initial_tau=initial_tau, use_newton=use_newton, verbose=verbose)
end

"""
min_composite(c, ∇c_adjoint, gauge_projection, x0, dims, opt_params, dual_gauge)

Value function approach to solve problems of the form
    min_{x} h(c(x))

"""
function min_composite(c::Function,
                       ∇c_adjoint::Function,
                       gauge_projection::Function,
                       dims::Tuple{Integer,Integer},
                       opt_params::SubproblemOptParams{T},
                       dual_gauge::Union{Function, Nothing}=nothing;
                       num_outer_iter::Integer=20,
                       x0::Union{Vector{T}, Nothing}=nothing,
                       initial_tau::Real=0.0,
                       use_newton::Bool=true,
                       verbose::Bool=false) where {T}
    if dual_gauge == nothing
        use_newton = false
    end
    problem_params = ProblemParams(c,∇c_adjoint,gauge_projection,dual_gauge)
    m,n = dims
    if x0 == nothing
        x = zeros(dims[2])
    else
        x = copy(x0)
    end
    τ1 = initial_tau
    info = IntermediateVars{T}(zeros(m),zeros(m),initial_tau)
    vf! = function(x,τ)
        info.tau = τ
        f,df = value_function!(x,problem_params, info, opt_params)
        return f, df
    end

    for i in 1:num_outer_iter
        if use_newton
            f1, df1 = vf!(x,τ1)
            verbose && @printf("itr %3d - τ: %3.3e f: %3.3e df: %3.3e\n",i,τ1,f1,df1)
            τ1 = τ1 - f1/df1
        else
            if i==1
                f1, df1 = vf!(x,τ1)
                τ2 = max((1.001)*τ1,τ1 + 0.001)
            end
            verbose && @printf("itr %3d - τ: %3.3e f: %3.3e \n",i,τ1,f1)
            f2, df2 = vf!(x,τ2)
            τ_new = τ2 - f2*(τ2 - τ1)/(f2 - f1)
            τ1, f1, df1 = τ2, f2, df2
            τ2 = τ_new
        end
    end
    return x
end

end


text:
function mre(predictions, actualdata)
	#Julia: simple version of mre( without checking the third parameter)
	#%if the third parameter is provided as true, then we will first convert
	#%everything into integers!
	#if nargin == 3 && varargin{1}
	#    predictions = ceil(predictions);
	#    actualdata = ceil(actualdata);
	#end;

	#% mean absolute error

	out = zeros(1,size(predictions,2));
	for i=1:size(predictions,2)
		if i > size(actualdata,2)
			continue
		end
		act = actualdata[:,i];
		pre = predictions[:,i];
		nonz = find(act);

		#%My own definition
		out[i] = mean(abs(pre[nonz]-act[nonz])./abs(act[nonz]));
		#%Weka!!
		#%avgAct = mean(act);
		#%out(i) = sum(abs(pre-act)) / sum(abs(avgAct-act));
	end
	out
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer