text:
export istotally_real, istotally_complex, conjugates, conjugates_real,
       conjugates_complex, conjugates_log, complex_conjugation

################################################################################
#
#  Totally real/complex predicates
#
################################################################################

@doc Markdown.doc"""
    istotally_real(K::NumberField) -> Bool

Returns true if and only if $K$ is totally real, that is, if all roots of the
defining polynomial are real.
"""
function istotally_real(K::NumField)
  return iszero(signature(K)[2])
end

@doc Markdown.doc"""
    istotally_complex(K::AnticNumberField) -> Bool

Returns true if and only if $K$ is totally complex, that is, if all roots of the
defining polynomial are not real.
"""
function istotally_complex(K::NumField)
  return iszero(signature(K)[1])
end

################################################################################
#
#  Conjugates and real embeddings
#
################################################################################

@doc Markdown.doc"""
    conjugates(x::nf_elem, abs_tol::Int) -> Vector{acb}

Compute the conjugates of `x` as elements of type `acb`.
Recall that we order the complex conjugates
$\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
$\sigma_{i}(x) = \overline{sigma_{i + s}(x)}$ for $r + 1 \leq i \leq r + s$.

Every entry `y` of the vector returned satisfies
`radius(real(y)) < 2^-abs_tol` and `radius(imag(y)) < 2^-abs_tol` respectively.
"""
function conjugates(x::NumFieldElem, abs_tol::Int = 32, T = arb)
  if T === arb
    return conjugates_arb(x, abs_tol)
  else
    error("Cannot return conjugates as type $T")
  end
end

# This is for quick and dirty computations
function __conjugates_arb(x::nf_elem, prec::Int = 32)
  K = parent(x)
  d = degree(K)
  r1, r2 = signature(K)
  conjugates = Array{acb}(undef, r1 + 2*r2)

  c = conjugate_data_arb_roots(K, -1)

  CC = AcbField(prec, cached = false)
  RR = ArbField(prec, cached = false)

  xpoly = arb_poly(parent(K.pol)(x), prec)

  for i in 1:r1
    o = RR()
    ccall((:arb_poly_evaluate, libarb), Nothing,
          (Ref{arb}, Ref{arb_poly}, Ref{arb}, Int),
           o, xpoly, c.real_roots[i], prec)

    if !isfinite(o)
      error("oops")
    end
    conjugates[i] = CC(o)
  end

  for i in 1:r2
    tacb = CC()
    ccall((:arb_poly_evaluate_acb, libarb), Nothing,
          (Ref{acb}, Ref{arb_poly}, Ref{acb}, Int),
           tacb, xpoly, c.complex_roots[i], prec)

    if !isfinite(tacb)
      error("oops")
    end

    conjugates[r1 + i] = tacb
    conjugates[r1 + i + r2] = conj(conjugates[r1 + i])
  end

  return conjugates
end

function conjugates_arb(x::nf_elem, abs_tol::Int = 32)
  K = parent(x)
  d = degree(K)
  r1, r2 = signature(K)
  conjugates = Array{acb}(undef, r1 + 2*r2)
  target_tol = abs_tol
  abs_tol = Int(floor(abs_tol * 1.1))

  while true
    prec_too_low = false
    c = conjugate_data_arb_roots(K, abs_tol)

    if abs_tol < 0
      error("Precision overflow in conjugates_arb_log")
    end

    CC = AcbField(abs_tol, cached = false)
    RR = ArbField(abs_tol, cached = false)

    xpoly = arb_poly(parent(K.pol)(x), abs_tol)

    for i in 1:r1
      o = RR()
      ccall((:arb_poly_evaluate, libarb), Nothing,
            (Ref{arb}, Ref{arb_poly}, Ref{arb}, Int),
             o, xpoly, c.real_roots[i], abs_tol)

      if !isfinite(o) || !radiuslttwopower(o, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
        break
      end
      conjugates[i] = CC(o)
    end

    if prec_too_low
      continue
    end

    for i in 1:r2
      tacb = CC()
      ccall((:arb_poly_evaluate_acb, libarb), Nothing,
            (Ref{acb}, Ref{arb_poly}, Ref{acb}, Int),
             tacb, xpoly, c.complex_roots[i], abs_tol)

      if !isfinite(tacb) || !radiuslttwopower(tacb, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
        break
      end

      conjugates[r1 + i] = tacb
      conjugates[r1 + i + r2] = conj(conjugates[r1 + i])
    end

    if prec_too_low
      continue
    end

    for i in 1:d
      expand!(conjugates[i], -target_tol)
    end
    return conjugates
  end
end

@doc Markdown.doc"""
    conjugates_arb_real(x::nf_elem, abs_tol::Int) -> Vector{arb}

Compute the real conjugates of `x` as elements of type `arb`.

Every entry `y` of the array returned satisfies
`radius(y) < 2^-abs_tol`.
"""
function conjugates_real(x::nf_elem, abs_tol::Int = 32, T = arb)
  if T === arb
    return conjugates_arb_real(x, abs_tol)
  else
    error("Cannot return real conjugates as type $T")
  end
end

function conjugates_arb_real(x::nf_elem, abs_tol::Int = 32)
  r1, r2 = signature(parent(x))
  c = conjugates_arb(x, abs_tol)
  z = Array{arb}(undef, r1)

  for i in 1:r1
    z[i] = real(c[i])
  end

  return z
end

@doc Markdown.doc"""
    conjugates_complex(x::nf_elem, abs_tol::Int) -> Vector{acb}

Compute the complex conjugates of `x` as elements of type `acb`.
Recall that we order the complex conjugates
$\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
$\sigma_{i}(x) = \overline{sigma_{i + s}(x)}$ for $r + 1 \leq i \leq r + s$.

Every entry `y` of the array returned satisfies
`radius(real(y)) < 2^-abs_tol` and `radius(imag(y)) < 2^-abs_tol`.
"""
function conjugates_complex(x::nf_elem, abs_tol::Int = 32, T = arb)
  if T === arb
    return conjugates_arb_complex(x, abs_tol)
  else
    error("Cannot return real conjugates as type $T")
  end
end

function conjugates_arb_complex(x::nf_elem, abs_tol::Int)
  r1, r2 = signature(parent(x))
  c = conjugates_arb(x, abs_tol)
  z = Vector{acb}(undef, r2)

  for i in (r1 + 1):(r1 + r2)
    z[i - r1] = c[i]
  end

  return z
end

################################################################################
#
#  Logarithms of conjugates
#
################################################################################

@doc Markdown.doc"""
    conjugates_arb_log(x::nf_elem, abs_tol::Int) -> Array{arb, 1}

Returns the elements
$(\log(\lvert \sigma_1(x) \rvert),\dotsc,\log(\lvert\sigma_r(x) \rvert),
\dotsc,2\log(\lvert \sigma_{r+1}(x) \rvert),\dotsc,
2\log(\lvert \sigma_{r+s}(x)\rvert))$ as elements of type `arb` with radius
less then `2^-abs_tol`.
"""
function conjugates_log(x::nf_elem, abs_tol::Int = 32, T = arb)
  if T === arb
    return conjugates_arb_log(x, abs_tol)
  else
    throw(error("Cannot return real conjugates as type ", T))
  end
end

function conjugates_arb_log(x::nf_elem, abs_tol::Int)
  K = parent(x)
  c = conjugate_data_arb_roots(K, abs_tol)
  r1 = length(c.real_roots)
  r2 = length(c.complex_roots)
  d = degree(K)
  target_tol = abs_tol

  # TODO: Replace this using multipoint evaluation of libarb
  z = Array{arb}(undef, r1 + r2)
  while true
    prec_too_low = false
    c = conjugate_data_arb_roots(K, abs_tol)
    if abs_tol > 2^20
      error("Something wrong in conjugates_arb_log")
    end
    xpoly = arb_poly(parent(K.pol)(x), abs_tol)
    RR = ArbField(abs_tol, cached = false)
    for i in 1:r1
      o = RR()
      ccall((:arb_poly_evaluate, libarb), Nothing,
            (Ref{arb}, Ref{arb_poly}, Ref{arb}, Int),
            o, xpoly, c.real_roots[i], abs_tol)
      abs!(o, o)
      log!(o, o)
      z[i] = o

      if !isfinite(z[i]) || !radiuslttwopower(z[i], -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
        break
      end
    end

    if prec_too_low
      continue
    end

    CC = AcbField(abs_tol, cached = false)

    tacb = CC()
    for i in 1:r2
      oo = RR()
      ccall((:arb_poly_evaluate_acb, libarb), Nothing,
            (Ref{acb}, Ref{arb_poly}, Ref{acb}, Int),
            tacb, xpoly, c.complex_roots[i], abs_tol)
      abs!(oo, tacb)
      log!(oo, oo)
      mul2exp!(oo, oo, 1)
      z[r1 + i] = oo

      if !isfinite(z[r1 + i]) || !radiuslttwopower(z[r1 + i], -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
        break
      end
    end

    if prec_too_low
      continue
    end

    for i in 1:r1 + r2
      zz = deepcopy(z[i])
      expand!(z[i], -target_tol)
    end
    return z
  end
end

function conjugates_arb_log(x::nf_elem, R::ArbField)
  z = conjugates_arb_log(x, R.prec)
  return map(R, z)
end

################################################################################
#
#  Minkowski map
#
################################################################################

@doc Markdown.doc"""
    minkowski_map(a::nf_elem, abs_tol::Int) -> Array{arb, 1}

Returns the image of $a$ under the Minkowski embedding.
Every entry of the array returned is of type `arb` with radius less then
`2^(-abs_tol)`.
"""
function minkowski_map(a::T, abs_tol::Int = 32) where T <: NumFieldElem
  z = _minkowski_map_and_apply(a, abs_tol, identity)
  return z
end

# The following function computes the minkowski_map, applies G to the output.
# G mus be a function (::Vector{arb}, abs_tol::Int) -> Bool, *
# where the first return value indicates if the result is good enough
function _minkowski_map_and_apply(a, abs_tol, G, work_tol::Int = abs_tol)
  K = parent(a)
  A = Array{arb}(undef, absolute_degree(K))
  c = conjugates_arb(a, work_tol)::Vector{acb}
  r, s = signature(K)

  for i = 1:r
    @assert isreal(c[i])
    A[i] = real(c[i])
    if !radiuslttwopower(A[i], -abs_tol)
      error("this should not happen")
    end
  end

  if work_tol > 2^18 || abs_tol > 2^18
    throw(error("asdsd"))
  end

  #R = ArbField(precision(parent(c[1])), false)
  R = ArbField(2 * work_tol, cached = false)
  sqrt2 = sqrt(R(2))

  for i in 1:s
    t = c[r + i]
    A[r + 2*i - 1] = sqrt2 * real(t)
    A[r + 2*i] = sqrt2 * imag(t)
    if !radiuslttwopower(A[r + 2*i], -abs_tol)
      return _minkowski_map_and_apply(a, abs_tol, G, Int(floor(work_tol * 2)))
    end
  end

  if typeof(G) === typeof(identity)
    return A
  else
    b, B = G(A, abs_tol)
    if b
      return B
    else
      return _minkowski_map_and_apply(a, abs_tol, G, 2*work_tol)
    end
  end
end

################################################################################
#
#  T_2
#
################################################################################

function t2(x::S, abs_tol::Int = 32, T = arb) where S <: NumFieldElem
  if T === arb
    g = function(w, abs_tol)
      z = mapreduce(y -> y^2, +, w)
      return radiuslttwopower(z, -abs_tol), z
    end
    c = _minkowski_map_and_apply(x, abs_tol, g)
    return c
  else
    throw(NotImplemented())
  end
end

############################################################################
#
#  Signs of real embeddings
#
############################################################################

#@doc Markdown.doc"""
##    _signs(a::nf_elem) -> Array{Int, 1}
#> For a non-zero element $a$ return the signs of all real embeddings.
#"""
function _signs(a::nf_elem)
  if iszero(a)
    error("element must not be zero")
  end

  p = 16
  r1, r2 = signature(parent(a))

  if r1 == 0
    return Int[]
  end

  s = Array{Int}(undef, r1)
  while true
    c = conjugates(a, p)
    done = true
    for i=1:r1
      if contains(reim(c[i])[1], 0)
        p *= 2
        done = false
        break
      end
      s[i] = ispositive(real(c[i])) ? 1 : -1
    end
    if done
      return s
    end
  end
end

#@doc Markdown.doc"""
##    signs(a::FacElem{nf_elem, AnticNumberField}) -> Array{Int, 1}
#> For a non-zero element $a$ in factored form,
#> return the signs of all real embeddings.
#"""
function _signs(a::FacElem{nf_elem, AnticNumberField})
  r1, r2 = signature(base_ring(a))
  if r1 == 0
    return Int[]
  end

  s = ones(Int, r1)

  for (k, e) = a.fac
    if iseven(e)
      continue
    end
    s .*= _signs(k)
  end
  return s
end

@doc Markdown.doc"""
    complex_conjugation(K::AnticNumberField)

Given a totally complex normal number field, this function returns an
automorphism which is the restriction of complex conjugation at one embedding.
"""
function complex_conjugation(K::AnticNumberField; auts::Vector{NfToNfMor} = NfToNfMor[])
  if !isempty(auts)
    A = auts
  else
    A = automorphisms(K)
    if length(A) < degree(K)
      error("Number field must be normal")
    end
  end
  a = gen(K)
  d = length(A)
  !istotally_complex(K) && error("Number field must be totally complex")
  #First, quick and dirty. If only one automorphism works, then we return it
  p = 32
  while true
    c = conjugates(a, p)
    for i = 1:d
      if !isinvolution(A[i])
        continue
      end
      cc = conj.(conjugates(image_primitive_element(A[i]), p))
      for k = 1:d
        if overlaps(c[k], cc[k])
          found = true
          for j = 1:d
            if j == k
              continue
            end
            if overlaps(c[j], cc[k])
              found = false
              break
            end
          end
          if found
            return A[i]
          end
        end
      end
    end
    p = 2 * p
    if p > 2^18
      error("Precision too high in complex_conjugation")
    end
  end
  error("something went wrong!")
end

function _find_complex_conjugation(K::AnticNumberField, A::Vector{NfToNfMor})
  a = gen(K)
  #First, quick and dirty. If only one automorphism works, then we return it
  p = 32
  while true
    c = conjugates(a, p)
    overlap = false
    for i = 1:length(A)
      if !isinvolution(A[i])
        continue
      end
      cc = conj.(conjugates(image_primitive_element(A[i]), p))
      for k = 1:degree(K)
        if overlaps(c[k], cc[k])
          overlap = true
          found = true
          for j = 1:degree(K)
            if j == k
              continue
            end
            if overlaps(c[j], cc[k])
              found = false
              break
            end
          end
          if found
            return true, A[i]
          end
        end
      end
    end
    if !overlap
      break
    end
    if p > 2^18
      error("Precision too high in complex_conjugation")
    end
  end
  return false, A[1]
end

function iscomplex_conjugation(f::NfToNfMor)
  K = domain(f)
  @assert K == codomain(f)
  !istotally_complex(K) && error("Number field must be totally complex")
  p = 32
  d = degree(K)
  a = gen(K)
  img_a = image_primitive_element(f)
  while true
    c = conjugates(a, p)
    cc = conj.(conjugates(img_a, p))
    for i = 1:degree(K)
      if !overlaps(c[i], cc[i])
        return false
      end
    end
    #Now I need to assure that the precision is enough.
    found = true
    for j in 1:d-1
      if overlaps(c[d], cc[j])
        found = false
        break
      end
    end
    if found
      return true
    end
    p = 2 * p
    if p > 2^18
      error("Precision too high in complex_conjugation")
    end
  end

end


text:
using StanModels, StanSample, MCMCChains, CSV

df = CSV.read(stanmodels_path("..", "data", "WaffleDivorce.csv"), DataFrame)
StanModels.scale!(df, [:MedianAgeMarriage])

# Define the Stan language model

m5_1s = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] divorce; // Predictor
 vector[N] median_age; // Outcome
}

parameters {
 real a; // Intercept
 real bA; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
  # priors
  a ~ normal(10, 10);
  bA ~ normal(0, 1);
  sigma ~ uniform(0, 10);

  # model
  divorce ~ normal(a + bA*median_age , sigma);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m_5_1s = SampleModel("m5.1s", m5_1s);

# Input data for cmdstan

m5_1_data = Dict("N" => length(df[!, :Divorce]), "divorce" => df[!, :Divorce],
    "median_age" => df[!, :MedianAgeMarriage_s]);

# Sample using cmdstan

rc = stan_sample(m_5_1s, data=m5_1_data);

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.69 0.22  9.34 10.03  2023    1
bA    -1.04 0.21 -1.37 -0.71  1882    1
sigma  1.51 0.16  1.29  1.79  1695    1
"

# Describe the draws
if success(rc)
  chn = read_samples(m_5_1s; output_format=:mcmcchains)
  show(chn)
end

ext: jl


text:
using Documenter
using BioCCP

makedocs(sitename="BioCCP",
    format = Documenter.HTML(),
    modules=[BioCCP],
    pages=Any["BioCCP"=> "man/BioCCP.md"])

#=
deploydocs(
            repo = "github.com/kirstvh/BioCCP.jl.git",
        )

=#

max_forks_repo_name: kirstvh/BioCCP


text:
suite["core"] = BenchmarkGroup()

s = BenchmarkGroup()
suite["core"]["edges"] = s

function b1(g)
    i = 0
    for e in edges(g)
        i += 1
    end
    return i
end

for G in GLIST
    for n in [10,100]
        g = CompleteGraph(n)
        s["fullyconn-$G-$n"] = @benchmarkable b1($g)
    end
end

for G in DGLIST
    for n in [10,100]
        g = CompleteDiGraph(n)
        s["fullyconn-$G-$n"] = @benchmarkable b1($g)
    end
end


text:
using
    MathOptInterface,
    BenchmarkTools,
    Random,
    IterTools,
    GaussianProcesses,
    Plots
    Random;

func(x) = (x[1]-2)^2*x[1]+x[2];

numSamps = 5
numPlot = 20
sampGrid = Vector{Float64}[]
plotGrid = Vector{Float64}[]

for p in product(LinRange(0,3,numSamps), LinRange(0,3,numSamps))
    push!(sampGrid, collect(p))
end

for p in product(LinRange(0,3,numPlot), LinRange(0,3,numPlot))
    push!(plotGrid, collect(p))
end

trueFunc = func.(plotGrid)

mz = MeanZero()
mc = MeanConst(1.)
#ml = MeanLin()
mp = MeanPoly([2. 3.; 1. 2.])
kern = SE(0.0, 0.0)
sampGridT = reduce(hcat, sampGrid)
sampVals = func.(sampGrid)+0.1*randn(numSamps^2)
gp = GP(sampGridT, sampVals, mz, kern)
GaussianProcesses.optimize!(gp)

attemptGrid = zeros(numPlot,numPlot)
for i in 1:length(plotGrid)
    attemptGrid[Int(trunc((i-1)/numPlot))+1, i-Int(trunc((i-1)/numPlot))*numPlot] = predict_f(gp, hcat(plotGrid[i]))[1][1]
end
trueGrid = zeros(numPlot,numPlot)
for i in 1:length(plotGrid)
    trueGrid[Int(trunc((i-1)/numPlot))+1, i-Int(trunc((i-1)/numPlot))*numPlot] = func(plotGrid[i])
end

plot(trueGrid, st=:surface, label = "True Function")
plot!(attemptGrid, st=:surface, label = "GP Attempt")

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using Base.Test
using StrPack

@struct type A
    a::UInt8
end

@struct type B
    a::Array{UInt8,2}(2,2)
end

@struct type C
    a::String(5)
end

@struct type D
    a::C
end

@struct type E
    a::Array{C,1}(2)
end

@struct immutable F
    a::Array{Float64,2}(3,2)
end

abstract type abstractG; end

@struct type G1 <: abstractG
  a::UInt8
end

@struct immutable G2 <: abstractG
  a::UInt8
end

@struct type Hvl_t
    len::Csize_t
    p::Ptr{Void}
end

function roundtrip(a)
    ios = IOBuffer()
    pack(ios, a)
    seek(ios, 0)
    isequal(a, unpack(ios, typeof(a)))
end

@test roundtrip(A(0xbc))

@test !roundtrip(B(zeros(UInt8, 0, 0)))
@test !roundtrip(B([0xab 0xcd]))
@test roundtrip(B([0x01 0x02; 0x03 0x04]))
@test !roundtrip(B([0x01 0x02 0xab; 0x03 0x04 0xcd; 0xac 0xbd 0xef]))

@test roundtrip(C(""))
@test roundtrip(C("a"))
@test roundtrip(C("asdfg"))
@test !roundtrip(C("asdfgh"))

@test roundtrip(D(C("")))

@test !roundtrip(E([C("a")]))
@test roundtrip(E([C("a"); C("b")]))

@test roundtrip(F([1. 2.; 3. 4.; 5. 6.]))

@test roundtrip(G1(UInt8(0)))
@test roundtrip(G2(UInt8(0)))

@test roundtrip(Hvl_t(UInt64(0), ccall(:jl_environ, Ptr{Void}, (Int,), 0)))


text:
function inser-sort!(A)
    if length(A) < 1 # check if A is empty
        return "A < 1"
    end
    for i in 2:length(A) # Julia uses 1 based indexing, no point at starting at 1 wchich has 1 element
        key = A[i] # Set key equal to the element we're currently checking
        j = i-1 # j will be the index checking the elements before i
        while j > 0 && key < A[j] # loop until value is less than the current element
            A[j+1] = A[j]
            j = j - 1
        end
        A[j+1] = key # set equal to the last index that was bigger than current
    end
end

A = [0,1,2,4,2]
#test push

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
using Luxor, Colors

# n-gons
@png begin
    for i in 23:-1:3
        sethue(i/20, 0.7, .99)
        ngonside(O, 75, i, 0, :fill)
        sethue("navy")
        ngonside(O, 75, i, 0, :stroke)
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# mdp step simulator and stepthrough
let
    mdp = GridWorld()
    solver = RandomSolver(MersenneTwister(2))
    policy = solve(solver, mdp)
    sim = StepSimulator("s,sp,r,a,ai", rng=MersenneTwister(3), max_steps=100)
    n_steps = 0
    for (s, sp, r, a, ai) in simulate(sim, mdp, policy)
        @test isa(s, state_type(mdp))
        @test isa(sp, state_type(mdp))
        @test isa(r, Float64)
        @test isa(a, action_type(mdp))
        @test isa(ai, Void)
        n_steps += 1
    end
    @test n_steps <= 100

    n_steps = 0
    for s in stepthrough(mdp, policy, "s", rng=MersenneTwister(4), max_steps=100)
        @test isa(s, state_type(mdp))
        n_steps += 1
    end
    @test n_steps <= 100
end

# pomdp step simulator and stepthrough
let
    mdp = BabyPOMDP()
    policy = FeedWhenCrying()
    up = PrimedPreviousObservationUpdater(true)
    sim = StepSimulator("s,sp,r,a,b,ui,i,ai", rng=MersenneTwister(3), max_steps=100)
    n_steps = 0
    for (s, sp, r, a, b, ui, i, ai) in simulate(sim, mdp, policy, up)
        @test isa(s, state_type(mdp))
        @test isa(sp, state_type(mdp))
        @test isa(r, Float64)
        @test isa(a, action_type(mdp))
        @test isa(b, Bool)
        @test ui == nothing
        @test ai == nothing
        @test ui == nothing
        n_steps += 1
    end
    @test n_steps == 100

    n_steps = 0
    for r in stepthrough(mdp, policy, "r", rng=MersenneTwister(4), max_steps=100)
        @test isa(r, Float64)
        @test r <= 0
        n_steps += 1
    end
    @test n_steps == 100
end

# example from stepthrough documentation
let
    pomdp = BabyPOMDP()
    policy = RandomPolicy(pomdp)

    for (s, a, o, r) in stepthrough(pomdp, policy, "s,a,o,r", max_steps=10)
        println("in state $s")
        println("took action $o")
        println("received observation $o and reward $r")
    end
end


text:
module InformationTools
using Distributions
export omega2, omega2!, fval, fval!, qrfactchk, calculate_f, calculate_omega2, computessr,
       computessr!, poisson_deviance_grid, poisson_deviance_grid_pval, poisson_deviance_pval,
       makeperms

function check_ss_effect(x)
    if x < -3e-9
        error("ss_effect = $x < -3e-9")
    end
end

function omega2!(out, fullrss::AbstractArray, fulldfe, reducedrss::AbstractArray, reduceddfe, n)
    broadcast!(omega2, out, fullrss, fulldfe, reducedrss, reduceddfe, n)
end

function omega2(fullrss, fulldfe, reducedrss, reduceddfe, n)
    ss_effect = reducedrss - fullrss
    check_ss_effect(ss_effect)
    df_effect = reduceddfe - fulldfe
    ms_effect = ss_effect/df_effect
    ms_error = fullrss/fulldfe
    (df_effect*(ms_effect - ms_error))/(df_effect*ms_effect + (n - df_effect)*ms_error)
end

function fval!(out, fullrss::AbstractArray, fulldfe, reducedrss::AbstractArray, reduceddfe)
    broadcast!(fval, out, fullrss, fulldfe, reducedrss, reduceddfe)
end

function fval(fullrss, fulldfe, reducedrss, reduceddfe)
    ss_effect = reducedrss - fullrss
    check_ss_effect(ss_effect)
    (ss_effect/(reduceddfe - fulldfe))./(fullrss/fulldfe)
end

qrfactchk(X) = ((r = rank(X)) == size(X, 2) || error("rank is $(r) but size is $(size(X, 2))"); qrfact(X))

function calculate_f(rast, dmat::LinAlg.QRCompactWY, reduced_dmat::LinAlg.QRCompactWY; stat::Symbol=:omega2)
    reducedrss = computessr(reduced_dmat, rast)
    fullrss = computessr(dmat, rast)
    fulldfe = size(rast, 1) - size(dmat, 2)
    reduceddfe = size(rast, 1) - size(reduced_dmat, 2)
    return fval.(fullrss, fulldfe, reducedrss, reduceddfe)
end

function calculate_omega2(rast, dmat::LinAlg.QRCompactWY, reduced_dmat::LinAlg.QRCompactWY; stat::Symbol=:omega2)
    reducedrss = computessr(reduced_dmat, rast)
    fullrss = computessr(dmat, rast)
    fulldfe = size(rast, 1) - size(dmat, 2)
    reduceddfe = size(rast, 1) - size(reduced_dmat, 2)
    return omega2.(fullrss, fulldfe, reducedrss, reduceddfe, size(rast, 1))
end

computessr(qrf::LinAlg.QRCompactWY, y::StridedMatrix) = vec(sumabs2(view(qrf[:Q]'y, size(qrf, 2)+1:size(qrf, 1), :), 1))
computessr!(qrf::LinAlg.QRCompactWY, y::StridedMatrix) =
    vec(sumabs2(view(LAPACK.gemqrt!('L', 'T', qrf.factors, qrf.T, y), size(qrf, 2)+1:size(qrf, 1), :), 1))
function computessr!(out::StridedVector, qrf::LinAlg.QRCompactWY, y::StridedMatrix)
    sumabs2!(reshape(out, 1, length(out)), view(LAPACK.gemqrt!('L', 'T', qrf.factors, qrf.T, y), size(qrf, 2)+1:size(qrf, 1), :))
    out
end

function cumsum2!(X)
    for j = 2:size(X, 2), i = 1:size(X, 1)
        @inbounds X[i, j] = X[i, j-1] + X[i, j]
    end
    X
end

@inline xlogy{T<:Real}(x::T, y::T) = x > zero(T) ? x * log(y) : zero(x)
function poisson_deviance_grid(Y, ningroup, minlength=1)
    @assert sum(ningroup) == size(Y, 1)
    cs = sum(Y, 1)
    cumsum2!(cs)
    groupcs = fill(-1, length(ningroup), size(Y, 2))
    lls = fill(NaN, size(Y, 2), size(Y, 2))
    for itp = 1:size(Y, 2)
        iidx = 0
        for igrp = 1:length(ningroup)
            grpend = iidx+ningroup[igrp]
            s = 0
            @simd for i = iidx+1:grpend
                @inbounds s += Y[i, itp]
            end
            groupcs[igrp, itp] = s
            iidx = grpend
        end
    end
    cumsum2!(groupcs)

    for istart = 1:size(Y, 2), iend = istart+minlength:size(Y, 2)
        groupll = 0.0
        iidx = 0
        for igrp = 1:length(ningroup)
            k = groupcs[igrp, iend] - (istart == 1 ? 0 : groupcs[igrp, istart-1])
            groupll += -k + xlogy(float(k), k/(ningroup[igrp]*(iend-istart+1)))
        end
        k = cs[iend] - (istart == 1 ? 0 : cs[istart-1])
        ll = -k + xlogy(float(k), k/(size(Y, 1)*(iend-istart+1)))
        lls[iend, istart] = groupll - ll
    end

    return lls
end

function poisson_deviance_grid_pval(Y, ningroup, perm, minlength=1)
    @assert sum(ningroup) == size(Y, 1)
    cs = sum(Y, 1)
    cumsum2!(cs)
    groupcs = fill(-1, length(ningroup), size(Y, 2))

    truevals = poisson_deviance_grid(Y, ningroup, minlength)
    trueval = maximum(truevals)
    isnan(trueval) && return NaN

    ngt = 0
    for ibs = 1:size(perm, 2)
        for itp = 1:size(Y, 2)
            iidx = 0
            for igrp = 1:length(ningroup)
                grpend = iidx+ningroup[igrp]
                s = 0
                @simd for i = iidx+1:grpend
                    @inbounds s += Y[perm[i, ibs], itp]
                end
                groupcs[igrp, itp] = s
                iidx = grpend
            end
        end
        cumsum2!(groupcs)

        bestll = -Inf
        for istart = 1:size(Y, 2), iend = istart+minlength:size(Y, 2)
            groupll = 0.0
            for igrp = 1:length(ningroup)
                k = groupcs[igrp, iend] - (istart == 1 ? 0 : groupcs[igrp, istart-1])
                groupll += -k + xlogy(float(k), k/(ningroup[igrp]*(iend-istart+1)))
            end
            k = cs[iend] - (istart == 1 ? 0 : cs[istart-1])
            ll = -k + xlogy(float(k), k/(size(Y, 1)*(iend-istart+1)))
            newll = groupll - ll
            bestll = max(newll, bestll)
        end

        ngt += bestll >= trueval
    end

    return ngt/size(perm, 2)
end

function poissonll(y, ningroup, perm, ibs)
    iidx = 0
    tk = 0
    groupll = 0.0
    for igrp = 1:length(ningroup)
        grpend = iidx+ningroup[igrp]
        k = 0
        @simd for i = iidx+1:grpend
            @inbounds k += y[perm[i, ibs]]
        end
        groupll += -k + xlogy(float(k), k/ningroup[igrp])
        tk += k
        iidx = grpend
    end
    groupll - (-tk + xlogy(float(tk), tk/length(y)))
end

function poisson_deviance_pval(y, ningroup, perm)
    @assert sum(ningroup) == size(y, 1)
    trueval = poissonll(y, ningroup, 1:length(y), 1)
    isnan(trueval) && return NaN

    ngt = 0
    for ibs = 1:size(perm, 2)
        ll = poissonll(y, ningroup, perm, ibs)
        ngt += ll >= trueval
    end

    return ngt/size(perm, 2)
end

function makeperms(n, nperm)
    perms = zeros(Int, n, nperm)
    inds = collect(1:n)
    for iperm = 1:nperm
        perms[:, iperm] = shuffle!(inds)
    end
    perms
end
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# Convenience functions for initializing AtomsBase systems.

# See NBodySimulator.generate_bodies_in_cell_nodes
function generate_atoms_in_cubic_cell(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, L::Unitful.Length, reference_temp::Unitful.Temperature)
    average_velocity = √(u"k" * reference_temp / elements[symbol].atomic_mass)
    velocities = average_velocity * randn(Float64, (3, n))
    particles = Vector{AtomsBase.Atom}(undef, n)

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x = dL/2:dL:L, y = dL/2:dL:L, z = dL/2:dL:L
        if count > n
            break
        end
        particles[count] = AtomsBase.Atom(symbol, SVector{3}(x, y, z), SVector{3}(velocities[:, count]))
        count += 1
    end

    bounding_box = SVector{3}(SVector{3}.(eachrow(L * I(3))))
    periodic_system(particles, bounding_box)
end


text:
module LinQuadInterpMod

import FEM: AbstractInterpolator, Point2, FENode2, Vertex4, inv2x2t!, det2x2,
       dNmatrix, Jmatrix, dNdxmatrix, Nvec, get_area

export LinQuadInterp

immutable LinQuadInterp <: AbstractInterpolator
    N::Vector{Float64}
    dN::Matrix{Float64}
    dNdx::Matrix{Float64}
    J::Matrix{Float64}
end

function LinQuadInterp()
    N = zeros(4)
    dN = zeros(4,2)
    dNdx = Array(Float64, 4, 2)
    J = Array(Float64, 2, 2)
    LinQuadInterp(N, dN, dNdx, J)
end
get_allocated_N(i::LinQuadInterp) = i.N
get_allocated_dN(i::LinQuadInterp) = i.dN
get_allocated_dNdx(i::LinQuadInterp) = i.dNdx
get_allocated_J(i::LinQuadInterp) = i.J

# Shape functions in local coords
function Nvec(interp::LinQuadInterp, loc_coords::Point2)

    ξ = loc_coords[1]
    η = loc_coords[2]

    N = get_allocated_N(interp)

    N[1] = (1 + ξ) * (1 + η) * 0.25
    N[2] = (1 - ξ) * (1 + η) * 0.25
    N[3] = (1 - ξ) * (1 - η) * 0.25
    N[4] = (1 + ξ) * (1 - η) * 0.25

    return N
end

function dNmatrix(interp::LinQuadInterp, loc_coords::Point2)
    ξ = loc_coords[1]
    η = loc_coords[2]

    dN = get_allocated_dN(interp)

    dN[1,1] =  (1 + η) * 0.25
    dN[1,2] =  (1 + ξ) * 0.25

    dN[2,1] = -(1 + η) * 0.25
    dN[2,2] =  (1 - ξ) * 0.25

    dN[3,1] = -(1 - η) * 0.25
    dN[3,2] = -(1 - ξ) * 0.25

    dN[4,1] =  (1 - η) * 0.25
    dN[4,2] = -(1 + ξ) * 0.25

    return dN
end

function Jmatrix(interp::LinQuadInterp,
                 vertices::Vertex4, nodes::Vector{FENode2},
                 dN::Matrix{Float64})

    J = get_allocated_J(interp)
    fill!(J, 0.0)
    for row in 1:size(dN, 1)
        x = nodes[vertices[row]].coords.x
        y = nodes[vertices[row]].coords.y

        J[1, 1] += dN[row, 1] * x
        J[1, 2] += dN[row, 1] * y
        J[2, 1] += dN[row, 2] * x
        J[2, 2] += dN[row, 2] * y
    end
    return J
end

function dNdxmatrix(interp::LinQuadInterp, local_coords::Point2,
                    vertices::Vertex4, nodes::Vector{FENode2})
        dN = dNmatrix(interp, local_coords)
        J = Jmatrix(interp, vertices, nodes, dN)
        A_mul_B!(interp.dNdx, dN, inv2x2t!(J))
        return interp.dNdx
end

end # module


text:
# run this script as `julia evaluate_performance_single.jl res.bson` or `julia evaluate_performance_single.jl dir`
# in the second case, it will recursively search for all compatible files in subdirectories
target = ARGS[1]

using DrWatson
@quickactivate
using EvalMetrics
using FileIO
using BSON
using DataFrames
using ValueHistories
using LinearAlgebra

function compute_stats(f::String)
	data = load(f)
	println(abspath(f))
	#map(p->println("$(p[1]) = $(p[2])"), collect(pairs(data[:parameters])))
	scores_labels = [(data[:val_scores], data[:val_labels]), (data[:tst_scores], data[:tst_labels])]
	setnames = ["validation", "test"]

	results = []
	for (scores, labels) in scores_labels
		scores = vec(scores)
		roc = EvalMetrics.roccurve(labels, scores)
		auc = EvalMetrics.auc_trapezoidal(roc...)
		prc = EvalMetrics.prcurve(labels, scores)
		auprc = EvalMetrics.auc_trapezoidal(prc...)

		t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
		cm5 = ConfusionMatrix(labels, scores, t5)
		tpr5 = EvalMetrics.true_positive_rate(cm5)
		f5 = EvalMetrics.f1_score(cm5)

		push!(results, [auc, auprc, tpr5, f5])
	end

	DataFrame(measure = ["AUC", "AUPRC", "TPR@5", "F1@5"], validation = results[1], test = results[2])
end

function query_stats(target::String)
	if isfile(target)
		try
			println(compute_stats(target))
		catch e
			@info "$target not compatible"
		end
	else
		query_stats.(joinpath.(target, filter(x->!(occursin("model_", string(x))), readdir(target))))
	end
end

query_stats(target)


text:
using Oceananigans.Architectures
using Oceananigans.Architectures: device
import Oceananigans.Architectures: architecture
using CUDA, CUDA.CUSPARSE
using KernelAbstractions: @kernel, @index

using LinearAlgebra, SparseArrays, IncompleteLU
using SparseArrays: nnz

import LinearAlgebra.ldiv!

"""
`JacobiPreconditioner`
    stores only the diagonal `d` of D⁻¹ where `D = diag(A)`
    is applied to `r` with a vector multiplication `d .* r`

`ILUFactorization`
    stores two sparse lower and upper trianguilar matrices `L` and `U` such that `LU ≈ A`
    is applied to `r` with `forward_substitution!(L, r)` followed by `backward_substitution(U, r)`
    constructed with `ilu(A, τ = drop_tolerance)`

`SparseInversePreconditioner`
    stores a sparse matrix `M` such that `M ≈ A⁻¹`
    is applied to `r` with a matrix multiplication `M * r`
    constructed with
    `simplified_inverse_preconditioner(A)`
        -> assumes that the sparsity of `M` is the same as the sparsity of `A`
    `sparse_approximate_preconditioner(A, ε = tolerance, nzrel = relative_maximum_number_of_elements)`
        -> starts constructing the sparse inverse of A from identity matrix until, either a tolerance (ε) is met or nnz(M) = nzrel * nnz(A)

The suggested preconditioners are

on the `CPU`
`ilu()` (superior to everything always and in every situation!)

on the `GPU`
`sparse_inverse_preconditioner()` (if `Δt` is constant)
`simplified_inverse_preconditioner()` or `JacobiPreconditioner` (if Δt is variable)

`ilu()` cannot be used on the GPU because preconditioning the solver with a direct LU (or Choleski) type
of preconditioner would require too much computation for the `ldiv!(P, r)` step completely hindering the performances

"""
validate_settings(T, arch, settings)                                 = settings
validate_settings(::Val{:Default}, arch, settings)                   = arch isa CPU ? (τ = 0.001,) : (ε = 0.1, nzrel = 1.0)
validate_settings(::Val{:SparseInverse}, arch, settings::Nothing)    = (ε = 0.1, nzrel = 1.0)
validate_settings(::Val{:ILUFactorization}, arch, settings::Nothing) = (τ = 0.001, )

validate_settings(::Val{:ILUFactorization}, arch, settings) = haskey(settings, :τ) ?
                                                                     settings :
                                                                     throw(ArgumentError("τ has to be specified for ILUFactorization"))
validate_settings(::Val{:SparseInverse}, arch, settings)    = haskey(settings, :ε) && haskey(settings, :nzrel) ?
                                                                     settings :
                                                                     throw(ArgumentError("both ε and nzrel have to be specified for SparseInverse"))

function build_preconditioner(::Val{:Default}, matrix, settings)
    default_method = architecture(matrix) isa CPU ? :ILUFactorization : :SparseInverse
    return build_preconditioner(Val(default_method), matrix, settings)
end

build_preconditioner(::Val{:None},              A, settings) = Identity()
build_preconditioner(::Val{:Jacobi},            A, settings) = JacobiPreconditioner(A)
build_preconditioner(::Val{:SparseInverse},     A, settings) = sparse_inverse_preconditioner(A, ε = settings.ε, nzrel = settings.nzrel)
build_preconditioner(::Val{:SimplifiedInverse}, A, settings) = simplified_inverse_preconditioner(A)
build_preconditioner(::Val{:ILUFactorization},  A, settings) = ilu(A, τ = settings.τ)

@inline architecture(::CuSparseMatrixCSC) = GPU()
@inline architecture(::SparseMatrixCSC)   = CPU()

struct JacobiPreconditioner{V}
    invdiag::V
end

function JacobiPreconditioner(A::AbstractMatrix)
    M    = size(A, 1)
    arch = architecture(A)
    dev  = device(arch)

    invdiag = arch_array(arch, zeros(eltype(A), M))
    col, row, val = unpack_constructors(arch, constructors(arch, A))

    loop! = _get_inv_diag!(dev, 256, M)
    event = loop!(invdiag, col, row, val; dependencies=Event(dev))
    wait(dev, event)

    return JacobiPreconditioner(invdiag)
end

@kernel function _multiply_in_place!(u, invdiag, v)
    i = @index(Global, Linear)
    @inbounds u[i] = invdiag[i] * v[i]
end

function  LinearAlgebra.ldiv!(u, precon::JacobiPreconditioner, v)
    invdiag = precon.invdiag
    arch = architecture(invdiag)
    dev  = device(arch)

    M = length(invdiag)
    loop! = _multiply_in_place!(dev, 256, M)
    event = loop!(u, invdiag, v; dependencies=Event(dev))
    wait(dev, event)
end

function  LinearAlgebra.ldiv!(precon::JacobiPreconditioner, v)
    ldiv!(v,precon,v)
end

abstract type AbstractInversePreconditioner{M} end

function  LinearAlgebra.ldiv!(u, precon::AbstractInversePreconditioner, v)
    mul!(u, matrix(precon), v)
end

function  LinearAlgebra.ldiv!(precon::AbstractInversePreconditioner, v)
    mul!(v, matrix(precon), v)
end

struct SparseInversePreconditioner{M} <: AbstractInversePreconditioner{M}
    Minv :: M
end

@inline matrix(p::SparseInversePreconditioner)  = p.Minv

function simplified_inverse_preconditioner(A::AbstractMatrix)

    arch                  = architecture(A)
    constr                = deepcopy(constructors(arch, A))
    colptr, rowval, nzval = copy_unpack_constructors(arch, constr)
    dev                   = device(arch)

    M = size(A, 1)

    diag = arch_array(arch, zeros(eltype(nzval), M))

    loop! = _get_diag!(dev, 256, M)
    event = loop!(diag, colptr, rowval, nzval; dependencies=Event(dev))
    wait(dev, event)

    loop! = _initialize_simplified_inverse_preconditioner!(dev, 256, M)
    event = loop!(nzval, colptr, rowval, diag; dependencies=Event(dev))
    wait(dev, event)

    constr_new = (colptr, rowval, nzval)

    Minv = arch_sparse_matrix(arch, constructors(arch, M, constr_new))

    return SparseInversePreconditioner(Minv)
end

@kernel function _initialize_simplified_inverse_preconditioner!(nzval, colptr, rowval, diag)
    col = @index(Global, Linear)

    for idx = colptr[col] : colptr[col+1] - 1
        if rowval[idx] == col
            nzval[idx] = diag[col]
        else
            nzval[idx] = - nzval[idx] * diag[rowval[idx]]
        end
    end
end

function sparse_inverse_preconditioner(A::AbstractMatrix; ε, nzrel)

   # let's choose an initial sparsity => diagonal
   A_cpu    = arch_sparse_matrix(CPU(), A)
   Minv_cpu = spai_preconditioner(A_cpu, ε = ε, nzrel = nzrel)

   Minv = arch_sparse_matrix(architecture(A), Minv_cpu)
   return SparseInversePreconditioner(Minv)
end


text:
# least squares penalty

"""
  LeastSquares(A, b, λ=1.0)

Returns the function `f(x) = (λ/2)⋅||Ax-b||^2`.
"""

type LeastSquares{R <: RealOrComplex} <: ProximableFunction
  A::AbstractArray{R,2}
  b::AbstractArray{R,1}
  Atb::AbstractArray{R,1}
  lambda::Real
  gamma::Real
  S::AbstractArray{R,2}
  U::Union{LinAlg.Cholesky, SparseArrays.CHOLMOD.Factor}
  function LeastSquares(A::AbstractArray{R,2}, b::AbstractArray{R,1}, lambda::Real)
    if size(A, 1) != length(b)
      error("A and b have incompatible dimensions")
    end
    if lambda <= 0
      error("lambda must be positive")
    end
    if size(A,1) >= size(A,2)
      new(A, b, A'*b, lambda, -1, A'*A)
    else
      new(A, b, A'*b, lambda, -1, A*A')
    end
  end
end

LeastSquares{R <: RealOrComplex}(A::AbstractArray{R,2}, b::AbstractArray{R,1}, lambda::Real=1.0) =
  LeastSquares{R}(A, b, lambda)

function (f::LeastSquares{R}){R <: RealOrComplex}(x::AbstractArray{R,1})
  return (f.lambda/2)*vecnorm(f.A*x - f.b, 2)^2
end

function factor_step!{R <: RealOrComplex}(f::LeastSquares{R}, gamma::Real)
  # factor step, two cases: (1) tall A, (2) fat A
  lamgam = f.lambda*gamma
  if issparse(f.A)
    f.U = cholfact(f.S; shift=1.0/lamgam)
  else
    f.U = cholfact(f.S + I/lamgam)
  end
  f.gamma = gamma
end

function prox!{R <: RealOrComplex}(f::LeastSquares{R}, x::AbstractArray{R,1}, y::AbstractArray{R,1}, gamma::Real=1.0)
  # if gamma different from f.gamma then call factor_step!
  if gamma != f.gamma
    factor_step!(f, gamma)
  end
  lamgam = f.lambda*gamma
  # solve step, two cases: (1) tall A, (2) fat A
  q = f.Atb + x/lamgam
  if size(f.A,1) >= size(f.A,2)
    y[:] = f.U\q
  else
    y[:] = lamgam*(q - (f.A'*(f.U\(f.A*q))))
  end
  return (f.lambda/2)*norm(f.A*y-f.b, 2)^2
end

fun_name(f::LeastSquares) = "least-squares penalty"
fun_dom{R <: Real}(f::LeastSquares{R}) = "AbstractArray{Real,1}"
fun_dom{C <: Complex}(f::LeastSquares{C}) = "AbstractArray{Complex,1}"
fun_expr(f::LeastSquares) = "x ↦ (λ/2)||A*x - b||^2"
fun_params(f::LeastSquares) = string("λ = $(f.lambda), A = ", typeof(f.A), " of size ", size(f.A), ", b = ", typeof(f.b), " of size ", size(f.b))

function prox_naive{R <: RealOrComplex}(f::LeastSquares, x::AbstractArray{R,1}, gamma::Real=1.0)
  lamgam = f.lambda*gamma
  y = (f.A'*f.A + I/lamgam)\(f.Atb + x/lamgam)
  fy = (f.lambda/2)*norm(f.A*y-f.b)^2
  return y, fy
end


text:
# x25519.jl gives some modular arithmetics modulo 2^255-19 and the point
# addition and scalar multiplication on curve25519. The computation is
# implemented with the intention of implementing it on top of an FPGA.

# Carry saved adder.
function csa(a,b,c,carry_in=0)
    s = (a ⊻ b) ⊻ c
    c = (((a & b) | (a & c) | (b & c)) << 1) | carry_in
    return (s, c)
end

# Modular addition of two numbers in 0:2^255-20 modulo p=2^255-19.
function addmod_25519(a, b)
#    return (a + b) % (BigInt(2)^255-19)
    zs, zc = csa(a, b, ~(BigInt(2)^255-19) & (BigInt(2)^255-1), 1)
    sel = ((zs+zc) >> 255) != 0 ? 1 : 0
    return (sel == 1) ? (zs+zc) & (BigInt(2)^255-1) : a+b
end

# Modular subtraction of two numbers in 0:2^255-20 modulo p=2^255-19.
function submod_25519(a, b)
#    t = (a - b) % (BigInt(2)^255 - 19)
#    return (t < 0) ? t + (BigInt(2)^255 - 19) : t
     zs, zc = csa(a, ~b & (BigInt(2)^255-1), 0)
     ws, wc = csa(a, ~b & (BigInt(2)^255-1), (BigInt(2)^255-19))
     sel = ((zs+zc+1) >> 255) != 0 ? 1 : 0
     return (sel == 1) ? (zs+zc+1) & (BigInt(2)^255-1) : (ws+wc+1) & (BigInt(2)^255-1)
end

# Modular multiplication of two numbers in 0:2^255-20 modulo p=2^255-19.
# with radix-8 interleave.
function multmod_25519(x, y)
#    return (x*y) % (BigInt(2)^255 - 19)
    p = BigInt(2)^255 - 19
    lut1 = [(0,0), (y,0), (y<<1, 0), (y<<1, y),
            (y<<2, 0), (y<<2, y), (y<<2, y<<1),
            csa(y<<2,y<<1,y)]
    lut2 = [i*8*19 for i=0:23]
    s = c = n = BigInt(0)
    for i in 85:-1:1
        ms, mc = lut1[(x >> (3*i-3)) + 1]
        x = x % BigInt(2)^(3*i-3)
        s1 = 8*(s % BigInt(2)^255)
        c1 = 8*(c % BigInt(2)^255)
        s2 = (s1 ⊻ ms) ⊻ c1
        c2 = ((s1 & ms) | (s1 & c1) | (ms & c1)) << 1
        s3 = (s2 ⊻ mc) ⊻ c2
        c3 = ((s2 & mc) | (s2 & c2) | (mc & c2)) << 1
        s = (s3 ⊻ n) ⊻ c3
        c = ((s3 & n) | (s3 & c3) | (n & c3)) << 1
        n = lut2[(s >> 255) + (c >> 255) + 1]
    end
    # TODO better reduction
    # s c = (s[254:0]+19*(s>>255), c[254:0]+19*(c>>255))
    s, c = csa((s % BigInt(2)^255), (c % BigInt(2)^255), n>>3)
    n = lut2[(s >> 255) + (c >> 255) + 1]
    #println("1st reduce $(s>>253) $(c>>253)")
    s, c = csa((s % BigInt(2)^255), (c % BigInt(2)^255), n>>3)
    z = s + c
    #println("2nd reduce $(s>>253) $(c>>253) $(z>>253)")
    if (z > BigInt(2)^255 - 19)
        z -= BigInt(2)^255 - 19
    end
    return z
end

# Point addition of two points on twisted Edwards curve 25519 where
# points are represented with (x, y, t, z) cordinates. If affine variable
# is true, then the result is affinized i.e. z=1.
function point_add(x1,y1,t1,z1,x2,y2,t2,z2; affine=false)
    r1 = submod_25519(y1, x1)
    r2 = submod_25519(y2, x2)

    r3 = addmod_25519(y1, x1)
    r4 = addmod_25519(y2, x2)

    r5 = multmod_25519(r1, r2)
    r6 = multmod_25519(r3, r4)

    r7 = multmod_25519(t1, t2)
    r8 = multmod_25519(z1, z2)

    k = 16295367250680780974490674513165176452449235426866156013048779062215315747161

    r7 = multmod_25519(k, r7)
    r8 = addmod_25519(r8, r8)

    r1 = submod_25519(r6, r5)
    r2 = submod_25519(r8, r7)

    r3 = addmod_25519(r8, r7)
    r4 = addmod_25519(r6, r5)

    x3 = multmod_25519(r1, r2)
    y3 = multmod_25519(r3, r4)

    t3 = multmod_25519(r1, r4)
    z3 = multmod_25519(r2, r3)

    if (affine == true)
        zinv = invmod_25519_M(z3)
        x3 = multmod_25519(x3, zinv)
        y3 = multmod_25519(y3, zinv)
        t3 = multmod_25519(x3, y3)
        z3 = BigInt(1)
    end

    return (x3, y3, t3, z3)
end

# Load precomputed radix-16 lookup table of the scalar multiplications of
# the base point on curve25519.
# The nth line of each files is the x,y,t-coordinates of (i*16^j)B expressed in
# 64 hexadecimal digits where i = n %16, j = n ÷ 16 and B is the base point.
function load_precomp_points()
    fx=open("../x25519/x_precomp_32k.dat", "r")
    fy=open("../x25519/y_precomp_32k.dat", "r")
    ft=open("../x25519/t_precomp_32k.dat", "r")
    xlines=readlines(fx)
    ylines=readlines(fy)
    tlines=readlines(ft)
    global x_precomp_data=map(x -> parse(BigInt,x,base=16), xlines)
    global y_precomp_data=map(x -> parse(BigInt,x,base=16), ylines)
    global t_precomp_data=map(x -> parse(BigInt,x,base=16), tlines)
    close(fx)
    close(fy)
    close(ft)
end

# Scalar multiplication of the base point on curve25519 with the 16-radix
# precomputed lookup table of scalar multiplications.
# Call load_precomp_points() before calling this function.
function scalarmultB(k)
    kl = [(k >> (4*i)) % 16 for i in 0:63]

    # k0 B from lut
    px = x_precomp_data[1+kl[1]]
    py = y_precomp_data[1+kl[1]]
    pt = t_precomp_data[1+kl[1]]
    pz = BigInt(1)
    for i in 2:64
        # k_n 16^n B from lut
        qx = x_precomp_data[1+(i-1)*16+kl[i]]
        qy = y_precomp_data[1+(i-1)*16+kl[i]]
        qt = t_precomp_data[1+(i-1)*16+kl[i]]
	qz = BigInt(1)
        # p = p + q
        #println("$(i): px $(px) py $(py) qx $(qx) qy $(qy)")
	px, py, pt, pz = point_add(px, py, pt, pz, qx, qy, qt, qz,
                                   affine=(i == 64))
    end
    return px, py, pz
end

function modpow2_n(x, n)
    for _ in 1:n
      x = multmod_25519(x,x)
    end
    return x
end

# Modular inverse module 2^255-19 by the Fermat's little theorem
function invmod_25519(a)
      a2   = multmod_25519(a,a)
      a4   = multmod_25519(a2,a2)
      a8   = multmod_25519(a4,a4)
      a9   = multmod_25519(a8,a)
      a11  = multmod_25519(a9,a2)
      a22  = multmod_25519(a11,a11)
      a_2_5  = multmod_25519(a22,a9)
      b_2_10 = modpow2_n(a_2_5,5)
      a_2_10 = multmod_25519(b_2_10,a_2_5)
      b_2_20 = modpow2_n(a_2_10,10)
      a_2_20 = multmod_25519(b_2_20,a_2_10)
      b_2_40 = modpow2_n(a_2_20,20)
      a_2_40 = multmod_25519(b_2_40,a_2_20)
      b_2_50 = modpow2_n(a_2_40,10)
      a_2_50 = multmod_25519(b_2_50,a_2_10)
      b_2_100 = modpow2_n(a_2_50,50)
      a_2_100 = multmod_25519(b_2_100,a_2_50)
      b_2_200 = modpow2_n(a_2_100,100)
      a_2_200 = multmod_25519(b_2_200,a_2_100)
      b_2_250 = modpow2_n(a_2_200,50)
      a_2_250 = multmod_25519(b_2_250,a_2_50)
      b_2_255 = modpow2_n(a_2_250,5)
      inv = multmod_25519(b_2_255,a11)
      return inv
end

function even(x)
    return ((x & 1) == 0) ? true : false
end

function signbit(x)
    return (x < 0) ? 1 : 0
end

function bitsize(M)
    if (M < 0)
        M = -M
    end
    return length(string(M,base=2))
end

# Montgomery modular inverse.
# Input X in 1:M-1 and M
# Output Lrs in 1:M-1 where Lrs = Xinv 2^n mod M
#  if real_inverse is false, then n = bitsize(M) else n = 0
#  i.e. when real_inverse is true, return the real inverse instead of
#  the almost inverse.
function inv_montgomery(X, M; real_inverse=false)
    # Phase1
    k = -bitsize(M)
    Luv = BigInt(0)
    Ruv = BigInt(X) << 1
    Lrs = BigInt(0)
    Rrs = BigInt(1)
    Luv = (Luv>>1)+Ruv
    Ruv = M
    Lrs = Lrs + Rrs
    Rrs = 0
    while true
        # println("Luv $(Luv) Ruv $(Ruv) Lrs $(Lrs) Rrs $(Rrs) k $(k)")
        SLuv, SRuv = signbit(Luv), signbit(Ruv)
        if (even(Luv>>1))
            if Luv == 0 # SLuv == signbit(-Luv)
                break
            end
            Luv = Luv>>1
            Rrs = Rrs<<1
            k = k+1
        else
            tmpuv = Luv>>1
            tmprs = Lrs
            Lrs = Lrs + Rrs
            if (SLuv ⊻ SRuv) == 1
                Luv = (Luv>>1)+Ruv
            else
                Luv = (Luv>>1)-Ruv
            end
            k = k+1
            ctrl = ((~SLuv & ~SRuv) | (~SLuv & SRuv)) & 1
            if signbit(Luv) == ctrl
                Ruv = tmpuv
                Rrs = tmprs<<1
            else
                Rrs = Rrs<<1
            end
        end
    end
    # Here we have Lrs == M
    Lrs = Lrs-Rrs
    Rrs = M
    if signbit(Lrs) == 1
        Lrs = Lrs + Rrs
    end

    # Phase 2
    if (real_inverse)
        k = k + bitsize(M)
    end

    while (k != 0)
        # println("Phase2 Lrs $(Lrs) Rrs $(Rrs) k $(k)")
        k = k - 1
        if even(Lrs)
            Lrs = Lrs>>1
        else
            Lrs = (Lrs+Rrs)>>1
        end
    end
    return Lrs
end

# Modular inverse modulo 2^255-19 with Montgomery modular inverse.
function invmod_25519_M(x)
    P = BigInt(2)^255-19
    iM = inv_montgomery(x, P, real_inverse=true)
    # iM = inv_montgomery(x, P)
    # i2n = inverse of 2^255 mod P = inverse of 19 mod P
    # i2n = 21330121701610878104342023554231983025602365596302209165163239159352418617876
    # return multmod_25519(iM, i2n)
    return iM
end


text:
using LinearAlgebra
using Random

using Jutils.Transforms

@testset "Shift" begin
    Random.seed!(201808071455)

    pts = rand(Float64, 4)
    shift = rand(Float64, 4)
    ptworkspace = rand(Float64, 4)
    ptoutput = rand(Float64, 4)
    mxworkspace = rand(Float64, 4, 4)
    mxoutput = rand(Float64, 4, 4)

    @test applytrans((Shift(shift),), pts, ptworkspace, ptoutput) == pts + shift
    @test applytrans_grad((Shift(shift),), pts, ptworkspace, ptoutput, mxworkspace, mxoutput) == Matrix(1.0I, 4, 4)
end

@testset "Updim" begin
    Random.seed!(201808151553)

    pts = rand(Float64, 4)
    new = rand()
    ptworkspace = rand(Float64, 5)
    ptoutput = rand(Float64, 5)

    @test applytrans((Updim{4,1}(new),), pts, ptworkspace, ptoutput) == [new, pts...]
    @test applytrans((Updim{4,2}(new),), pts, ptworkspace, ptoutput) == [pts[1], new, pts[2:4]...]
    @test applytrans((Updim{4,3}(new),), pts, ptworkspace, ptoutput) == [pts[1:2]..., new, pts[3:4]...]
    @test applytrans((Updim{4,4}(new),), pts, ptworkspace, ptoutput) == [pts[1:3]..., new, pts[4]]
    @test applytrans((Updim{4,5}(new),), pts, ptworkspace, ptoutput) == [pts..., new]

    ptworkspace = rand(Float64, 2)
    ptoutput = rand(Float64, 2)
    mxworkspace = rand(Float64, 2, 2)
    mxoutput = rand(Float64, 2, 2)

    pts = rand(Float64, 1)
    @test applytrans_grad((Updim{1,1}(0.0),), pts, ptworkspace, ptoutput, mxworkspace, mxoutput) == [0.0 1.0; 1.0 0.0]
    @test applytrans_grad((Updim{1,2}(0.0),), pts, ptworkspace, ptoutput, mxworkspace, mxoutput) == [1.0 0.0; 0.0 -1.0]
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
using Plots, Printf, LinearAlgebra

@views av(A) = 0.5 .* (A[1:end-1] .+ A[2:end])

@views function diffusion_1D(; implicit=false, do_viz=true)
    # Physics
    lx     = 10.0       # domain size
    ttot   = 1.0        # total simulation time
    dt     = 0.2        # physical time step
    # Numerics
    nx     = 512        # numerical grid resolution
    tol    = 1e-8       # tolerance
    itMax  = 1e4        # max number of iterations
    nout   = 10         # tol check
    CFL    = 0.8        # CFL number: 1.0 implicit, 0.8 explicit
    ε      = 1e-2       # small number to avoid division by 0
    # Derived numerics
    dx     = lx / nx    # grid size
    Vpdτ   = CFL * dx
    xc     = LinRange(dx/2, lx - dx/2, nx)
    # Array allocation
    qHx    = zeros(nx-1)
    qHx2   = zeros(nx-1)
    ResH   = zeros(nx-2)
    Re     = zeros(nx  )
    θr_dτ  = zeros(nx-1)
    dτ_ρ   = zeros(nx-2)
    # Initial condition
    H0     = exp.(-(xc .- lx/2).^2)
    Hold   = ones(nx) .* H0
    H      = ones(nx) .* H0
    D      = H.^3
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot / dt))
    # Physical time loop
    while it < nt
        iter = 0; err = 2 * tol
        # Pseudo-transient iteration
        while err > tol && iter < itMax
            # Diffusion coefficient
            D          .= H.^3
            # Assign iter params
            Re         .= π .+ sqrt.(π^2 .+ lx^2 ./ max.(D,ε) ./ dt)
            θr_dτ      .= lx ./ Vpdτ ./ av(Re)
            dτ_ρ       .= Vpdτ .* lx ./ max.(D[2:end-1],ε) ./ Re[2:end-1]
            # PT updates
            if implicit
                qHx        .= (qHx .* θr_dτ .- av(D) .* diff(H) ./ dx) ./ (1.0 .+ θr_dτ)
                H[2:end-1] .= (H[2:end-1] .+  dτ_ρ .* (Hold[2:end-1] ./ dt .- diff(qHx) ./ dx)) ./ (1.0 .+ dτ_ρ ./ dt)
            else
                qHx        .= qHx        .+  1.0 ./ θr_dτ .* ( .-qHx .-av(D) .* diff(H) ./ dx )
                H[2:end-1] .= H[2:end-1] .+  dτ_ρ  .* ( .-(H[2:end-1] .- Hold[2:end-1]) ./ dt .- diff(qHx) ./ dx)
            end
            iter += 1
            # Check error (explicit residual)
            if iter % nout == 0
                qHx2 .= .-av(D) .* diff(H) ./ dx
                ResH .= .-(H[2:end-1] .- Hold[2:end-1]) ./ dt .- diff(qHx2) ./ dx
                err = norm(ResH) / sqrt(length(ResH))
            end
        end
        # Update H
        Hold .= H
        ittot += iter; it += 1; t += dt
    end
    @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, nx, ittot)
    # Visualise
    if do_viz plot(xc, Array(H0), linewidth=3); display(plot!(xc, Array(H), legend=false, framestyle=:box, linewidth=3, xlabel="lx", ylabel="H", title="nonlinear diffusion (nt=$it, iters=$ittot)")) end
    return
end

diffusion_1D()


text:
# Adapted from:
# Example 3.99 of
# Blekherman, G.; Parrilo, P. & Thomas, R.
# Semidefinite Optimization and Convex Algebraic Geometry
# Society for Industrial and Applied Mathematics, 2012

@testset "[BPT12] Example 3.99 with $(factory.constructor)" for factory in sdp_factories
    @polyvar x y

    m = SOSModel(factory)

    @variable m α

    @constraint(m, x^2 + α*y <= 10, domain = @set x^2 + y^2 == 1)

    @objective m Max α

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

    @test JuMP.value(α) ≈ 6 atol=1e-6

    # FIXME JuMP does not support modification once loaded on jump/moi yet
#    @objective m Min α
#
#    JuMP.optimize!(m)
#
#    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
#
#    @test JuMP.value(α) ≈ -6 atol=1e-6
end


text:
using Kinetic, ProgressMeter, Plots, LinearAlgebra, JLD2, Flux
using Flux: onecold, @epochs
using ProgressMeter: @showprogress

cd(@__DIR__)
@load "data.jld2" X1 Y1 X2 Y2

device = cpu

X1 = Float32.(X1) |> device
Y1 = Float32.(Y1) |> device
X2 = Float32.(X2) |> device
Y2 = Float32.(Y2) |> device

nn = Chain(
    RNN(7, 28, relu),
    #Dense(7, 28, relu),
    RNN(28, 56, relu),
    RNN(56, 28, relu),
    Dense(28, 2),
)

data = Flux.Data.DataLoader((X1, Y1), shuffle = true) |> device
ps = params(nn)
sqnorm(x) = sum(abs2, x)
loss(x, y) = sum(abs2, nn(x) - y) / size(x, 2) #+ 1e-6 * sum(sqnorm, ps)
cb = () -> println("loss: $(loss(X1, Y1))")
opt = ADAM()

for epoch in 1:10
    println("epoch: $epoch, loss: $(loss(X1, Y1))")
    gs = Flux.gradient(ps) do
        loss(X1, Y1)
    end
    Flux.Optimise.update!(opt, ps, gs)
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
@testset "709.to-lower-case.jl" begin
    @test to_lower_case("QWaqwfkjAvaOdEakgbFbk") == "qwaqwfkjavaodeakgbfbk"
end


text:
using ViewerGL
using Distributed
if nprocs() == 1
    addprocs(3)
end
using LarSurf

block_size = [20, 20, 20]
LarSurf.lsp_setup(block_size)

segmentation = LarSurf.generate_truncated_sphere(10)

V, FV = LarSurf.lsp_get_surface(segmentation)
FVtri = LarSurf.triangulate_quads(FV)
objlines = LarSurf.Lar.lar2obj(V, FVtri, "tetris_tri.obj")

ViewerGL.VIEW([
    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

Vs = LarSurf.Smoothing.smoothing_FV_taubin(V, FV, 0.4, -0.3, 50)

ViewerGL.VIEW([
    ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

objlines = LarSurf.Lar.lar2obj(Vs, FVtri, "tetris_tri_taubin.obj")


text:
# Lucas Ondel, 2021

# Set of standards labels used across corpora
const LABELS = Dict(
    :speechunit => "speech-unit",
    :nonspeechunit => "non-speech-unit",
    :sil => "<sil>",
    :unk => "<unk>",
)

# Standard file names
const LOCATIONS = Dict(
    :wavs => "wav.scp",
    :local => ".local",
    :lang => "lang",
    :speakers => "speakers",
    :utt2spk => "uttids_speakers",
    :uttids => "uttids",
    :ali => "ali",
    :trans => "trans",
    :words => "words",
    :units => "units",
    :lexicon => "lexicon",
)

# Non-speech units (silence, breath, ...)
const NON_SPEECH_UNITS = Set([
    LABELS[:sil],
    LABELS[:unk]
])

max_issues_repo_name: lucasondel/SpeechDatasets

max_forks_repo_name: lucasondel/SpeechDatasets


text:
using RPOMDPToolbox
using RPOMDPModels
using RobustValueIteration
using Base.Test
const RPBVI = RobustValueIteration

unc = 0.1
prob = SimpleBaby2RPOMDP(-5.0, -10.0, 0.9, unc)
b0 = DiscreteBelief(prob, [0.5, 0.5])
@test pdf(b0,:hungry) == 0.5
@test pdf(b0,:full) == 0.5

# testing updater initialization
bs = [[b, 1-b] for b in 0.0:0.05:1.0]
solver = RPBVISolver(beliefpoints = bs, max_iterations = 100)
policy = RPBVI.solve(solver, prob)
up = RobustUpdater(prob, policy.alphas)
isd = initial_state_distribution(prob)
b4 = initialize_belief(up, isd)
@test pdf(b4, :hungry) == pdf(isd, :hungry)
@test pdf(b4, :full) == pdf(isd, :full)

# testing update function; if we feed baby, it won't be hungry
a = :feed
o = :crying
b4p = update(up, b4, a, o)
@show b4p.b
@test pdf(b4p,:hungry) ≈ 0.0 atol = 1e-3
@test pdf(b4p,:full) ≈ 1.0 atol = 1e-3

a = :nothing
o = :crying
b4p = update(up, b4, a, o)
@show b4p.b

@test isapprox(pdf(b4p,:hungry), 0.470588, atol=1e-4)
@test isapprox(pdf(b4p,:full), 0.52941, atol=1e-4)


text:
using Statistics
using LinearAlgebra
using PyPlot
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../test/useful_routines.jl")
include("transformation-elasticity-solver.jl")

function lame_lambda(k, m)
    return k - 2m / 3
end

function bulk_modulus(l, m)
    return l + 2m / 3
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

polyorder = 3
nelmts = 33
penaltyfactor = 1e2

folderpath = "examples/potential/cube/"
filename =
    "polyorder-" * string(polyorder) * "-nelmts-" * string(nelmts) * ".png"

width = 1.0              # mm
K1, K2 = 247.0, 192.0    # GPa
mu1, mu2 = 126.0, 87.0   # GPa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 3.93e3           # Kg/m^3
rho2 = 3.68e3           # Kg/m^3
V01 = 1.0 / rho1
V02 = 1.0 / rho2

ΔG0 = -14351.0/0.147

theta0 = -0.067
transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

meshwidth = [width, width]
numqp = required_quadrature_order(polyorder) + 2

dx = width / nelmts
penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

basis = TensorProductBasis(2, polyorder)
corner = [0.8,0.8]

mesh, cellquads, facequads, interfacequads = construct_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    corner,
    numqp,
)

nodaldisplacement = free_slip_bc_nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    penalty,
)

refseedpoints, spatialseedpoints, seedcellids =
    CutCellDG.seed_zero_levelset_with_interfacequads(interfacequads, mesh)
normals = CutCellDG.collect_interface_normals(interfacequads, mesh)

spatialpoints = spatialseedpoints[1, :, :]
productreferencepoints = refseedpoints[1, :, :]
parentreferencepoints = refseedpoints[2, :, :]
productreferencecellids = seedcellids[1, :]
parentreferencecellids = seedcellids[2, :]

angularposition = angular_position(spatialpoints)
sortidx = sortperm(angularposition)

parentreferencepoints = parentreferencepoints[:, sortidx]
productreferencepoints = productreferencepoints[:, sortidx]
parentreferencecellids = parentreferencecellids[sortidx]
productreferencecellids = productreferencecellids[sortidx]
spatialpoints = spatialpoints[:, sortidx]
normals = normals[:, sortidx]
angularposition = angularposition[sortidx]

parentstrain = parent_strain_at_reference_points(
    nodaldisplacement,
    basis,
    parentreferencepoints,
    parentreferencecellids,
    mesh,
)
parentstress = parent_stress(parentstrain, stiffness)

productstrain = product_strain_at_reference_points(
    nodaldisplacement,
    basis,
    theta0,
    productreferencepoints,
    productreferencecellids,
    mesh,
)
productstress = product_stress(productstrain, stiffness, theta0)

parentstrainenergy = V02 * strain_energy(parentstress, parentstrain)
productstrainenergy = V01 * strain_energy(productstress, productstrain)

parentradialtraction = CutCellDG.traction_force_at_points(parentstress, normals)
parentsrr = CutCellDG.traction_component(parentradialtraction, normals)

productradialtraction =
    CutCellDG.traction_force_at_points(productstress, normals)
productsrr = CutCellDG.traction_component(productradialtraction, normals)

parentdilatation = CutCellDG.dilatation(parentstrain)
productdilatation = CutCellDG.dilatation(productstrain)

parentcompwork = V02 * (1 .+ parentdilatation) .* parentsrr
productcompwork = V01 * (1 .+ productdilatation) .* productsrr

parentpotential = parentstrainenergy - parentcompwork
productpotential = productstrainenergy - productcompwork

potentialdifference = (productpotential - parentpotential) * 1e9 / abs(ΔG0)

fig, ax = PyPlot.subplots()
ax.plot(angularposition, potentialdifference, color = "black")
ax.grid()
ax.set_ylim(0.0,0.45)
ax.set_xlabel("Angular Position (deg)")
ax.set_ylabel(L"([\Phi] - [G0])/[G0]")
fig.tight_layout()
fig
fig.savefig(folderpath * filename)


text:
type InventoryModel <: RustModel
    dims::NTuple{5,Int}             #dd,dc,ds,dy,dr=dims
    pricelevel::Float64
    spotbuy::Float64
    wholesaleprices::Vector{Float64}
    uaa::Array{Float64,2}
    ubb::Array{Float64,1}

    ### RustModel interface
    beta::Float64
    u::Array{Float64,2}                 # u[k,a]: flow utilities
    v::Array{Float64,2}                 # v[k,a]: choice-specific discounted value
    V::Array{Float64,1}                 # V[k]: expected value
    pi::Array{Float64,3}                # pi[k,k',a]: condition state transitions
    p::Array{Float64,2}                 # p[k,a]: conditional choice probabilities (CCP)
end

# dcs2k(dims::Tuple,id,ic,is)=sub2ind(dims,id,ic,is)
# k2dcs(dims::Tuple,ik)=ind2sub(dims,ik)

#id = 1+ demand
#ic = wholesaleprice category
#is = 1 + stock at the beginning of the period
#iy = 1 + sales
#ir = 1 + stock next perido (decision variable)
#iq (implicit) = orders = sales + stock next period - stock today = (iy-1) + (ir - 1) - (is - 1)
#I can "sell" above demand for 0 profit
#I can "sell" above stock by paying a high spotprice (should only happen with extremely high random utility shocks)
#I can destroy goods for free with iq=iy+ir-is <0
uconstant(model::InventoryModel,id,ic,is,iy,ir)=
    model.pricelevel*(min(id-1,iy-1)-model.wholesaleprices[ic]*max(iy-1+ir-is,0)-model.spotbuy*max(iy-is,0))
ualpha(model::InventoryModel,id,ic,is,iy,ir)=-model.pricelevel*(is-1)
ubeta(model::InventoryModel,id,ic,is,iy,ir)=-model.pricelevel*1(iy-1+ir-is>0)

function ua(model::InventoryModel,id,ic,is,iy,ir,itheta)
    itheta==1?ualpha(model::InventoryModel,id,ic,is,iy,ir):ubeta(model::InventoryModel,id,ic,is,iy,ir)
end

function inventorymodel(dims,pricelevel,spotbuy,wholesaleprices)
    (dd,dc,ds)=dims
    dtheta=2
    dy,dr=ds,ds
    dk=dd*dc*ds
    da=dy*dr
    model=InventoryModel(
        (dd,dc,ds,dy,dr),
        pricelevel,
        spotbuy,
        wholesaleprices,
        zeros(dk*da,dtheta),  #uaa
        zeros(dk*da),         #ubb
        0.95,
        zeros(dk,da),       #u
        zeros(dk,da),       #v
        zeros(dk),          #V
        zeros(dk,dk,da),    #pi
        zeros(dk,da)        #p
        )
    uaa=[ua(model,id,ic,is,iy,ir,itheta) for id=1:dd,ic=1:dc,is=1:ds,iy=1:dy,ir=1:dr,itheta=1:dtheta]
    model.uaa=reshape(uaa,dk*da,dtheta)
    ubb=[uconstant(model,id,ic,is,iy,ir) for id=1:dd,ic=1:dc,is=1:ds,iy=1:dy,ir=1:dr]
    model.ubb=vec(ubb)
    model
end

function u!(model::InventoryModel,theta)
    model.u[:]=model.uaa*theta+model.ubb
end

function pi!(model::InventoryModel,pidc::AbstractMatrix)
    dd,dc,ds,dy,dr=model.dims
    dk=dd*dc*ds
    da=dy*dr
    for ik=1:dk
        idc=ind2sub((dd*dc,ds),ik)[1]
        for jk=1:dk
            jdc,js=ind2sub((dd*dc,ds),jk)
            for ia=1:da
                iy,ir=ind2sub((dy,dr),ia)
                model.pi[ik,jk,ia]=pidc[idc,jdc]*1(ir==js)
            end
        end
    end
end

# inventorymodel((dd,dc,ds),pricelevel,spotbuy,wholesaleprices)
# im1(pricelevel)=inventorymodel((5,3,10),pricelevel,10.0,[0.9,1.0,1.1])

function pidc1(model::InventoryModel)
    dd,dc,ds,dy,dr=model.dims
    pid=Tridiagonal(fill(0.1,dd-1),fill(0.8,dd),fill(0.1,dd-1))
    @assert dc==3
    pid[1,2]=0.2
    pid[dd,dd-1]=0.2
    pic=[0.1 0.9 0; 0.05 0.9 0.05; 0 0.9 0.1]
    kron(pic,pid)    #reverse order is on purpose the get the right layout
end

function pidc2(model::InventoryModel)
    dd,dc,ds,dy,dr=model.dims
    @assert dd==5
    @assert dc==3
    pid=[
    0.95 0.0 0.0 0.0 0.05;
    0.0 0.0 0.0 0.0 1.0;
    0.0 0.0 0.0 0.0 1.0;
    0.0 0.0 0.0 0.0 1.0;
    0.05 0.0 0.0 0.0 0.95;
    ]
    pic=[0.1 0.9 0; 0.05 0.9 0.05; 0 0.9 0.1]
    kron(pic,pid)    #reverse order is on purpose the get the right layout
end

function detensor(model::InventoryModel)
    dd,dc,ds,dy,dr=model.dims
    pidc=reshape(model.pi[:,:,1],(dd,dc,ds,dd,dc,ds))
    pic=sum(pidc[1,:,1,:,:,:],[2,4])[:,1,:,1]
    pid=sum(pidc[:,1,1,:,:,:],[3,4])[:,:,1,1]
    pid,pic
end

function twostep(model::InventoryModel)
    "todo"
    # V=(I-model.beta*model.pi[:,:,1])\(model.u[:,1]-log(model.p[:,1]))
    # thetau2=(I-model.beta*model.pi[:,:,2])*V+log(model.p[:,2])
    # model.u2\thetau2
end

#################### Identification

#a=[ A_1  I-betaPi_1; ... ; A_da  I-betaPi_da ]

# function idrank(model::InventoryModel,eta)
#     dk,da=size(model.u)
#     dtheta=size(model.uaa,2)
#     b=[1(ind2sub((dk,da),ika)[1]==jk)-model.beta*model.pi[ind2sub((dk,da),ika)[1],jk,ind2sub((dk,da),ika)[2]] for ika=1:dk*da,jk=1:dk]
#     a=hcat(sparse(model.uaa),sparse(b))
#     aa=full(a'*a)
#     ai=inv(aa)[1:dtheta,:]
#     function ff(eta1)
#         pi!(model,eta1)
#         piv=vcat([model.pi[:,:,ia]*model.V for ia=1:da]...)
#         model.beta*ai*(a'*piv)
#     end
#     jpi=Calculus.finite_difference_jacobian(ff,eta)
#     jpi
# end


text:
# C:\Users\wilhe\AppData\Local\Programs\Julia-1.6.0\bin\julia.exe --track-allocation=user
# include("C:\\Users\\wilhe\\Desktop\\Package Development\\EAGO.jl\\benchmark\\memory_allocation_analysis.jl"); exit()
#using Coverage; m = analyze_malloc("C:\\Users\\wilhe\\Desktop\\Package Development\\EAGO.jl\\src")
#filter!(x -> x.bytes > 0, m)

using JuMP, EAGO

m = Model(EAGO.Optimizer)
set_optimizer_attribute(m, "output_iterations", 1)
set_optimizer_attribute(m, "iteration_limit", 5)
set_optimizer_attribute(m, "verbosity", 0)

# Define bounded variables
xL = [10.0; 0.0; 0.0; 0.0; 0.0; 85.0; 90.0; 3.0; 1.2; 145.0]
xU = [2000.0; 16000.0; 120.0; 5000.0; 2000.0; 93.0; 95.0; 12.0; 4.0; 162.0]
@variable(m, xL[1] <= x[1] <= xU[1])
@NLobjective(m, Max, -2.0 * x[1])

JuMP.optimize!(m)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# ---
# title: 1338. Reduce Array Size to The Half
# id: problem1338
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/reduce-array-size-to-the-half/description/>
# hidden: true
# ---
#
# Given an array `arr`.  You can choose a set of integers and remove all the
# occurrences of these integers in the array.
#
# Return _the minimum size of the set_ so that **at least** half of the integers
# of the array are removed.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,3,3,3,5,5,5,2,2,7]
#     Output: 2
#     Explanation: Choosing {3,7} will make the new array [5,5,5,2,2] which has size 5 (i.e equal to half of the size of the old array).
#     Possible sets of size 2 are {3,5},{3,2},{5,2}.
#     Choosing set {2,7} is not possible as it will make the new array [3,3,3,3,5,5,5] which has size greater than half of the size of the old array.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [7,7,7,7,7,7]
#     Output: 1
#     Explanation: The only possible set you can choose is {7}. This will make the new array empty.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1,9]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: arr = [1000,1000,3,7]
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: arr = [1,2,3,4,5,6,7,8,9,10]
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `arr.length` is even.
#   * `1 <= arr[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
using CGAL
using Test

using CGAL.Cxx

cxx"""
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3          Point_3;
typedef K::Vector_3         Vector_3;
typedef K::Triangle_3       Triangle_3;
typedef K::Tetrahedron_3    Tetrahedron_3;
"""

indexable_to_array(v::Cxx.CxxCore.CppValue, n) = [icxx"auto r =$(v)[$i]; r;" for i = 0:(n-1)]
indexable_to_array(v::Cxx.CxxCore.CppPtr, n) = [icxx"auto r =(*($(v)))[$i]; r;" for i = 0:(n-1)]

@testset "Vector_3 - construct/add" begin
    pv1 = @cxxnew Vector_3(10, 20, 30)
    pv2 = @cxxnew Vector_3(1, 2, 3)

    v3 = icxx"*$pv1 + *$pv2;"

    @test 11 == icxx"auto r = $(v3).x(); r;"
    @test 22 == icxx"auto r = $(v3).y(); r;"
    @test 33 == icxx"auto r = $(v3).z(); r;"

    @test indexable_to_array(pv1, 3) == [10, 20, 30]
    @test indexable_to_array(pv2, 3) == [1, 2, 3]
    @test indexable_to_array(v3, 3) == [11, 22, 33]
end

@testset "Triangle_3" begin
    tri1 = icxx"Triangle_3(Point_3(1,-1,0), Point_3(-1,-1,0), Point_3(0,1,0));"
    tri2 = icxx"Triangle_3(Point_3(1,-1,1), Point_3(-1,-1,1), Point_3(0,1,1));"
    tri3 = icxx"Triangle_3(Point_3(1,0,-1), Point_3(-1,0,-1), Point_3(0,0,1));"
    @test icxx"do_intersect($tri1, $tri1);" == true
    @test icxx"do_intersect($tri1, $tri2);" == false
    @test icxx"do_intersect($tri1, $tri3);" == true

    @test (@cxx CGAL::do_intersect(tri1, tri3)) == true

    @test icxx"$tri1.squared_area();" == 4
end

@testset "Tetrahedron_3" begin
    tet1 = icxx"Tetrahedron_3(Point_3(1,-1,0), Point_3(-1,-1,0), Point_3(0,1,0), Point_3(0,0,1));"
    tri4 = icxx"Triangle_3(Point_3(1,-1,0.5), Point_3(-1,-1,0.5), Point_3(0,1,0.5));"
    tri5 = icxx"Triangle_3(Point_3(1,-1,10.0), Point_3(-1,-1,10.0), Point_3(0,1,10.0));"
    @test (@cxx CGAL::do_intersect(tet1, tri4)) == true
    @test (@cxx CGAL::do_intersect(tet1, tri5)) == false
end


text:
module MosekExtraMathProgTests

using MathProgBase
using MathProgBase.SolverInterface
using Mosek
using FactCheck

solver = MosekSolver(QUIET=true)

facts("[mathprogextra]") do

    # context("addelmtest") do
    #     objtol = 1e-7
    #     primaltol = 1e-6
    #     m = MathProgBase.model(Mosek.MosekSolver())

    #     loadproblem!(m,
    #                  sparse([1,1,1],[1,2,3],[1.0,1.0,1.0]), # A
    #                  Float64[-1.0,-1.0,-1.0],Float64[1.0,1.0,1.0], # lbx,ubx
    #                  zeros(Float64,3), # c
    #                  Float64[1.0], Float64[1.0], # lbc, ubc
    #                  :Max)

    #     j0  = MathProgBase.addvar!(m, 0.0, Inf, 0.0)
    #     i0  = MathProgBase.addconstr!(m, [4], [1.0], 1.0, 1.0)
    #     k0 = MathProgBase.addquadconstr!(m, [], [], [1,2,3], [1,2,3], [-1.0,1.0,1.0], '<', 0.0)

    #     j1  = MathProgBase.addvar!(m, 0.0, Inf, 0.0)
    #     i1  = MathProgBase.addconstr!(m, [5], [1.0], 1.0, 1.0)
    #     k1 = MathProgBase.addquadconstr!(m, [], [], [1,3], [2,3], [-1.0,1.0], '<', 0.0)

    #     println("j0 = ",j0,", i1 = ",i0,", k0 = ",k0)
    #     println("j1 = ",j1,", i1 = ",i1,", k1 = ",k1)

    #     Mosek.writedata(m.task,"mathprogtestextra.opf")
    #     println("End")

    #     #MathProgBase.optimize!(m)
    #     #stat=MathProgBase.status(m) # = :Unbounded
    # end

    context("dualsigntest") do
        duals=true
        # Problem 4 - lo1 from MOSEK docs
        # Property: All duals are non-zero
        # min      [ 3 1 5 1 ] * [ x y z w ]
        #  st [35]-[ 3 1 2   ]               ZERO
        #     [15]-[ 2 1 3 1 ] * [ x y z w ] NONPOS
        #     [25]-[   2   3 ]               NONNEG
        #     [10]-[   1     ]               NONPOS
        #     x,y,z,w > 0
        # Test: Signs of duals for ZERO, NONNEG, NONPOS
        println("Problem 4")
        duals = true
        if duals
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
                                      [ -3.0, -1.0, -5.0, -1.0 ],
                                      [  3.0   1.0   2.0   0.0 ;
                                       2.0   1.0   3.0   1.0 ;
                                       0.0   2.0   0.0   3.0 ;
                                       0.0   1.0   0.0   0.0 ],
                                      [ 30.0, 15.0, 25.0, 10.0 ],
            Any[(:Zero,1),(:NonPos,2),(:NonNeg,3),(:NonNeg,4)],
            [(:NonNeg,1:4)])
            MathProgBase.optimize!(m)
            @fact MathProgBase.status(m) --> :Optimal
            d = MathProgBase.getdual(m)
            # d[1] is free

            @fact d[2] --> less_than(1e-6)
            @fact d[3] --> greater_than(-1e-6)
            @fact d[4] --> greater_than(-1e-6)
        end
    end

    # context("cqo1test1") do
    #     s=MosekSolver()
    #     duals=true
    #     # Problem cqo1 from MOSEK docs
    #     # min  x4 + 0.5 * x5 + x6
    #     # s.t. x1 + x2 + 2 x3= 1
    #     #      x4^2  > x1^2 + x2^2
    #     #      x5*x6 > x3^2
    #     #      x1,x2,x3,x4,x5,x6 > 0

    #     # Input as conic, add conic constraints as quadratic constraints
    #     let m = MathProgBase.Conic(s)
    #         MathProgBase.loadproblem!(m,
    #                                   #  x1  x2  x3  x4  x5  x6
    #                                   [ 0.0 0.0 0.0 1.0 0.5 1.0 ], # c
    #                                   [ 1.0 1.0 2.0 0.0 0.0 0.0 ], # A
    #                                   [ 1.0 ], # b
    #                                   [ (:Zero,1) ],   # constr domain
    #                                   [ (:NonNeg,1:6) ]) # var    domain
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 4,1,2 ], Int32[4,1,2], Float64[-1.0, 1.0, 1.0], '<', 0.0 )
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 5,3   ], Int32[6,3  ], Float64[-1.0, 1.0     ], '<', 0.0 )

    #         MathProgBase.optimize!(m)

    #         @fact MathProgBase.status(m) --> :Optimal

    #         pobj = MathProgBase.getobjval(m)

    #         @fact pobj --> roughly(7.07106782e-01, 1e-6)

    #         xx = MathProgBase.getsolution(m)
    #         xc = MathProgBase.getconstrsolution(m)

    #         @fact (pobj-dot(xx,[ 0.0,0.0,0.0,1.0,0.5,1.0 ])) --> roughly(0.0,1e-6)
    #     end
    # end
    # context("cqo1test2") do
    #     s=MosekSolver()
    #     duals=true

    #     # Input as linear, add conic constraints as quadratic constraints
    #     let m = MathProgBase.LinearQuadraticModel(s)
    #         MathProgBase.loadproblem!(m,
    #                                   [ 1.0 1.0 2.0 0.0 0.0 0.0 ], # A
    #                                   [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ], # blx
    #                                   [ Inf, Inf, Inf, Inf, Inf, Inf ], # bux
    #                                   [ 0.0, 0.0, 0.0, 1.0, 0.5, 1.0 ], # c
    #                                   [ 1.0 ], # blc
    #                                   [ 1.0 ], # buc
    #                                   :Min )
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 4,1,2 ], Int32[4,1,2], Float64[-1.0, 1.0, 1.0], '<', 0.0 )
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 5,3   ], Int32[6,3  ], Float64[-1.0, 1.0     ], '<', 0.0 )

    #         MathProgBase.optimize!(m)

    #         @fact MathProgBase.status(m) --> :Optimal

    #         pobj = MathProgBase.getobjval(m)

    #         @fact pobj --> roughly(7.07106782e-01, 1e-6)

    #         xx = MathProgBase.getsolution(m)
    #         xc = MathProgBase.getconstrsolution(m)
    #         @fact (pobj-dot(xx,[ 0.0,0.0,0.0,1.0,0.5,1.0 ])) --> roughly(0.0,1e-6)

    #     end
    # end

    # context("cqo1test4") do
    #     s=MosekSolver()
    #     duals=true
    #     # Input as linear, add conic constraints as quadratic constraints
    #     let m = MathProgBase.model(s)
    #         MathProgBase.loadproblem!(m,
    #                                   [ 1.0 1.0 2.0 0.0 0.0 0.0 ], # A
    #                                   [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ], # blx
    #                                   [ Inf, Inf, Inf, Inf, Inf, Inf ], # bux
    #                                   [ 0.0, 0.0, 0.0, 1.0, 0.5, 1.0 ], # c
    #                                   [ 1.0 ], # blc
    #                                   [ 1.0 ], # buc
    #         :Min )
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 2,1,4 ], Int32[2,1,4], Float64[ 1.0, 1.0,-1.0], "<", 0.0 )
    #         MathProgBase.addquadconstr!(m, [],[], Int32[ 3,5   ], Int32[3,6  ], Float64[ 1.0,-1.0     ], "<", 0.0 )

    #         MathProgBase.optimize!(m)

    #         @fact MathProgBase.status(m) --> :Optimal

    #         pobj = MathProgBase.getobjval(m)

    #         @fact pobj --> roughly(7.07106782e-01, 1e-6)

    #         xx = MathProgBase.getsolution(m)
    #         xc = MathProgBase.getconstrsolution(m)

    #         @fact (pobj-dot(xx,[ 0.0,0.0,0.0,1.0,0.5,1.0 ])) --> roughly(0.0,1e-6)
    #     end
    # end

    context("qcqo1test") do
        duals=true
        # Problem qcqo1 from MOSEK docs
        # min     - x2      + x1^2 -     x1*x3+ 0.1 x2^2 +     x3^2

        # s.t. x1 + x2 + x3 - x1^2 + 0.2 x1*x3 -    x2^2 - 0.1 x3^2 >= 1.0
        #      x1,x2,x3 >= 0

        m = MathProgBase.LinearQuadraticModel(solver)
        MathProgBase.loadproblem!(m,
                                  zeros(Float64,(0,3)), # A
                                  [ 0.0, 0.0, 0.0 ], # blx
                                  [ Inf, Inf, Inf ], # bux
                                  [ 0.0, -1.0, 0.0 ], # c
                                  Float64[], # blc
                                  Float64[], # buc
                                  :Min)
        MathProgBase.setquadobj!(m,[1,1,2,3],[1,3,2,3], [2.0,-1.0,0.2,2.0]) # specify lower triangular
        MathProgBase.addquadconstr!(m,  # specify quadratid _terms_ (not triangular)
                                    [1,2,3], [1.0,1.0,1.0], # linear part
                                    [1,1,2,3], [1,3,2,3], [-1.0,0.2,-1.0,-0.1],
                                    '>',
                                    1.0)
        MathProgBase.optimize!(m)

        @fact MathProgBase.status(m) --> :Optimal

        pobj = MathProgBase.getobjval(m)

        @fact pobj --> roughly(-4.91766743e-01,1e-6)

        xx = MathProgBase.getsolution(m)

        @fact pobj - (-xx[2]+xx[1]^2-xx[1]*xx[3]+0.1*xx[2]^2+xx[3]^2) --> roughly(0.0,1e-6)

        #xc = MathProgBase.getquadconstrsolution(m)
        #@test_approx_eq_eps xc[1] (xx[1]+xx[2]+xx[3] - xx[1]^2 + 0.2*xx[1]*xx[3] - xx[2]^2 - 0.1 * xx[3]^2) 1e-6
    end

    context("modifymodeltest") do
        mmin = MathProgBase.LinearQuadraticModel(solver)
        mmax = MathProgBase.LinearQuadraticModel(solver)

        MathProgBase.loadproblem!(mmin,
                     [  1.0   0.0   0.0   0.0 ;
                        0.0   1.0   0.0   0.0 ;
                        0.0   0.0   1.0   0.0 ;
                        0.0   0.0   0.0   1.0 ],
                     [ -Inf, -Inf, -Inf, -Inf ], # blx
                     [  Inf,  Inf,  Inf,  Inf ], # bux
                     [  1.0,  1.0,  1.0,  1.0 ], # c
                     [  0.0,  0.0,  0.0,  0.0 ], # blc
                     [  0.0,  0.0,  0.0,  0.0 ], # buc
                     :Min)

        MathProgBase.loadproblem!(mmax,
                     [  1.0   0.0   0.0   0.0 ;
                        0.0   1.0   0.0   0.0 ;
                        0.0   0.0   1.0   0.0 ;
                        0.0   0.0   0.0   1.0 ],
                     [ -Inf, -Inf, -Inf, -Inf ], # blx
                     [  Inf,  Inf,  Inf,  Inf ], # bux
                     [  1.0,  1.0,  1.0,  1.0 ], # c
                     [  0.0,  0.0,  0.0,  0.0 ], # blc
                     [  0.0,  0.0,  0.0,  0.0 ], # buc
                     :Max)

        writeproblem(mmax,"testmax_orig.opf")
        # test: modify buc
        setconstrUB!(mmax,[ 1.0, 2.0, 3.0, 4.0 ]);
        MathProgBase.optimize!(mmax)
        @fact MathProgBase.status(mmax) --> :Optimal
        xx = MathProgBase.getsolution(mmax)

        @fact xx[1] --> roughly(1.0, 1e-8)
        @fact xx[2] --> roughly(2.0, 1e-8)
        @fact xx[3] --> roughly(3.0, 1e-8)
        @fact xx[4] --> roughly(4.0, 1e-8)

        writeproblem(mmax,"testmax.opf")
        #writeproblem(mmax,"testmin.task")

        # test: modify blc
        setconstrLB!(mmin,[ -1.0, -2.0, -3.0, -4.0 ]);
        MathProgBase.optimize!(mmin)
        @fact MathProgBase.status(mmin) --> :Optimal
        xx = MathProgBase.getsolution(mmin)

        @fact xx[1] --> roughly(-1.0, 1e-8)
        @fact xx[2] --> roughly(-2.0, 1e-8)
        @fact xx[3] --> roughly(-3.0, 1e-8)
        @fact xx[4] --> roughly(-4.0, 1e-8)

        # test: modify bux
        setvarUB!(mmax,[ 0.5, 1.5, 2.5, 3.5 ]);
        MathProgBase.optimize!(mmax)
        @fact MathProgBase.status(mmax) --> :Optimal
        xx = MathProgBase.getsolution(mmax)

        @fact xx[1] --> roughly(0.5, 1e-8)
        @fact xx[2] --> roughly(1.5, 1e-8)
        @fact xx[3] --> roughly(2.5, 1e-8)
        @fact xx[4] --> roughly(3.5, 1e-8)

        # test: modify blx
        setvarLB!(mmin,[ -0.5, -1.5, -2.5, -3.5 ]);
        MathProgBase.optimize!(mmin)
        @fact MathProgBase.status(mmin) --> :Optimal
        xx = MathProgBase.getsolution(mmin)

        @fact xx[1] --> roughly(-0.5, 1e-8)
        @fact xx[2] --> roughly(-1.5, 1e-8)
        @fact xx[3] --> roughly(-2.5, 1e-8)
        @fact xx[4] --> roughly(-3.5, 1e-8)
    end
end

end # module


text:
struct PixPos
     pix::UInt8
     pos::NamedTuple{(:x, :y, :z),Tuple{Float64, Float64, Float64}}
end
PixPos(pix::UInt8, pos::Vector{Float64}) = PixPos(pix, (x=pos[1], y=pos[2], z=pos[3]))
PixPos(pix::UInt8, pos::Array{Float64,1}) = PixPos(pix, (x=pos[1], y=pos[2], z=pos[3]))

@testset "Homography to Pose" begin

    fx = 1000.
    fy = 1000.
    cx = 320.
    cy = 240.

    K = [fx 0  cx;
          0 fy cy]

    C = [K; 0 0 1]

    detector = AprilTagDetector()
    detector.quad_decimate = 1.0 #NOTE see line 84 in runtest.jl

    #_____________________________________________________________________

    #test offset only
    test_cTw = [[1.0  0.0  0.0  -100.0;
                 0.0  1.0  0.0  -100.0;
                 0.0  0.0  1.0  1500.0;
                 0.0  0.0  0.0     1.0]]

    #test offset only
    push!(test_cTw, [1.0  0.0  0.0   320.0;
                     0.0  1.0  0.0  -200.0;
                     0.0  0.0  1.0  1500.0;
                     0.0  0.0  0.0     1.0])

    #test horizontal angle
    push!(test_cTw, [ 0.707107  0.0  0.707107  0.0;
                      0.0       1.0  0.0       0.0;
                     -0.707107  0.0  0.707107  1500.0;
                      0.0       0.0  0.0       1.0])

    #test offset and horizontal angle
    push!(test_cTw, [ 0.707107  0.0  0.707107  100.0;
                      0.0       1.0  0.0       100.0;
                     -0.707107  0.0  0.707107  1000.0;
                      0.0       0.0  0.0       1.0])

    #test vertical angle
    push!(test_cTw, [ 1.0  0.0       0.0       0.0;
                      0.0  0.707107  0.707107  0.0;
                      0.0 -0.707107  0.707107  1500.0;
                      0.0  0.0       0.0       1.0])

    #test 2 angles
    push!(test_cTw, [0.707107   0.0       -0.707107     0.0;
                     0.5        0.707107   0.5        -60.6602;
                     0.5       -0.707107   0.5       2060.66;
                     0.0        0.0        0.0          1.0])

    for idx = 1:length(test_cTw)
        cTw = test_cTw[idx]

        back = ones(UInt8, 301,301)*0x40
        back[11:290,11:290] =  ones(UInt8, 280,280)*0x80

        #tag 0 draw in camera frame
        #_____________________________________________________________________
        tag0 = kron(reinterpret(UInt8,getAprilTagImage(0)), ones(UInt8, 20,20))
        back[51:250,51:250] = tag0'

        x = collect(-150.:150)
        y = collect(-150.:150)
        z = 0.

        tag0cloud = Array{PixPos,2}(undef, 301,301)
        for i = 1:301, j = 1:301
           tag0cloud[i,j] = PixPos(back[i,j],[x[i]; y[j]; z])
        end

        #Projection
        #_____________________________________________________________________

        projImg = zeros(Float32, 480,640);
        cnt = zeros(Int, 480,640);

        for i = 1:301, j = 1:301
            tc = tag0cloud[i,j]
            pos = C*cTw[1:3,1:4]*[tc.pos.x; tc.pos.y; tc.pos.z; 1]
            pos /= pos[3]
            u = round(Int,pos[2])
            v = round(Int,pos[1])
            if 1 < u < 480 && 1 < v < 640
                cnt[u,v] += 1f0
                n = cnt[u,v]
                projImg[u,v] =  projImg[u,v]* (n-1)/n + tc.pix/255f0/n
            end
        end

        tags = detector(Gray{N0f8}.(projImg))

        pose = homographytopose(tags[1].H, fx, fy, cx, cy, taglength = 160.)
        # display(pose)

        if idx != 2 #FIXME? skip this test on 2 (bigger offsets like this one fails test), homographytopose is just not accurate enough
            @test all(isapprox.(pose[1:3,1:3], cTw[1:3,1:3], atol = 0.05))
            @test all(isapprox.(pose[1:3,4], cTw[1:3,4], atol = 10.))
        end
        pose,err1 = tagOrthogonalIteration(tags[1], fx, fy, cx, cy, taglength = 160.)
        @test all(isapprox.(pose[1:3,1:3], cTw[1:3,1:3], atol = 0.05))
        @test all(isapprox.(pose[1:3,4], cTw[1:3,4], atol = 10.))

    end

end


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using Plots
using Printf
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity, SurfaceTKEFlux, MixingLength
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize

Nz = 64
Ex, Ey = (1, 3)
sz = ColumnEnsembleSize(Nz=Nz, ensemble=(Ex, Ey))

ensemble_grid = RectilinearGrid(size = sz,
                                       halo = ColumnEnsembleSize(Nz=1),
                                       z = (-128, 0),
                                       topology = (Flat, Flat, Bounded))

default_closure = CATKEVerticalDiffusivity()
closure_ensemble = [default_closure for i = 1:Ex, j = 1:Ey]

Qᵇ = 1e-7
Qᵇ_ensemble = [Qᵇ for i = 1:Ex, j = 1:Ey]

closure_ensemble[1, 1] = CATKEVerticalDiffusivity(mixing_length = MixingLength(Cᴷc⁻=1.0, Cᴸᵇ=0.0))
closure_ensemble[1, 2] = CATKEVerticalDiffusivity(mixing_length = MixingLength(Cᴷc⁻=1.0, Cᴸᵇ=2.0))
closure_ensemble[1, 3] = CATKEVerticalDiffusivity(mixing_length = MixingLength(Cᴷc⁻=1.0, Cᴸᵇ=4.0))

b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ_ensemble))

model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
                                    closure = closure_ensemble,
                                    boundary_conditions = (; b=b_bcs),
                                    tracers = (:b, :e),
                                    buoyancy = BuoyancyTracer())

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

z = znodes(model.tracers.b)

simulation = Simulation(model, Δt = 1minute/2, stop_time = 0.0)

b = model.tracers.b
bz = Field(∂z(b))

function column_bz(j)
    compute!(bz)
    return view(interior(bz), 1, j, :)
end

Nt = 10
hs = zeros(Ey, Nt)
times = zeros(Nt)
for n = 1:Nt
    simulation.stop_time += 4hour
    run!(simulation)

    h = [-z[argmax(column_bz(j))] for j = 1:3]

    hs[:, n] .= h
    times[n] = model.clock.time
end

b1 = view(interior(b), 1, 1, :)
b2 = view(interior(b), 1, 2, :)
b3 = view(interior(b), 1, 3, :)

compute!(bz)
zbz = znodes(bz)

bz1 = view(interior(bz), 1, 1, :)
bz2 = view(interior(bz), 1, 2, :)
bz3 = view(interior(bz), 1, 3, :)

t = model.clock.time
b_plot = plot(b1, z, linewidth = 2, label = @sprintf("b1 t = %s", prettytime(t)), xlabel = "Buoyancy", ylabel = "z", legend=:bottomright)
plot!(b_plot, b2, z, linewidth = 2, linestyle=:dash, label = @sprintf("b2 t = %s", prettytime(t)))
plot!(b_plot, b3, z, linewidth = 2, linestyle=:dot, label = @sprintf("b3 t = %s", prettytime(t)))

bz_plot = plot(bz1, zbz, linewidth = 2, label = @sprintf("bz1 t = %s", prettytime(t)), xlabel = "Buoyancy", ylabel = "z", legend=:bottomright)
plot!(bz_plot, bz2, zbz, linewidth = 2, linestyle=:dash, label = @sprintf("bz2 t = %s", prettytime(t)))
plot!(bz_plot, bz3, zbz, linewidth = 2, linestyle=:dot, label = @sprintf("bz3 t = %s", prettytime(t)))

bbz_plot = plot(b_plot, bz_plot, layout=(1, 2), size=(1200, 600))

display(bbz_plot)


text:
using LinearAlgebra
using OpenCL

const sum_kernel = "
   __kernel void sum(__global const float *a,
                     __global const float *b,
                     __global float *c)
    {
      int gid = get_global_id(0);
      c[gid] = a[gid] + b[gid];
    }
"

const len = 50_000
#const len = 1

a = rand(Float32, len)
b = rand(Float32, len)

device, ctx, queue = cl.create_compute_context()
println(device)

a_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=a)
b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b)
c_buff = cl.Buffer(Float32, ctx, :w, length(a))

p = cl.Program(ctx, source=sum_kernel) |> cl.build!
k = cl.Kernel(p, "sum")

queue(k, size(a), nothing, a_buff, b_buff, c_buff)

r = cl.read(queue, c_buff)

if isapprox(norm(r - (a+b)), zero(Float32))
    println("Success!")
else
    println("Norm should be 0.0f")
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Plane(p, u, v)

A plane embedded in R³ passing through point `p`,
defined by non-parallel vectors `u` and `v`.

    Plane(p, n)

Alternatively specify point `p` and a given normal
vector `n` to the plane.
"""
struct Plane{T} <: Primitive{3,T}
  p::Point{3,T}
  u::Vec{3,T}
  v::Vec{3,T}
end

function Plane{T}(p::Point{3,T}, n::Vec{3,T}) where {T}
  # origin of coordinate system
  o = Vec{3,T}(zero(T), zero(T), zero(T))

  uv = Vec{3,T}[]
  for i in 1:3
    # subtract projection of Euclidean basis onto normal
    e = Vec{3,T}(ntuple(j -> j == i ? one(T) : zero(T), 3))
    v = e - (e⋅n)/(n⋅n) * n

    # check if it is a valid vector
    isapprox(v, o, atol=atol(T)) || push!(uv, v)

    # we need two vectors
    length(uv) == 2 && break
  end

  Plane{T}(p, uv[1], uv[2])
end

Plane(p::Point{3,T}, n::Vec{3,T}) where {T} = Plane{T}(p, n)

Plane(p::Tuple, u::Tuple, v::Tuple) = Plane(Point(p), Vec(u), Vec(v))

Plane(p::Tuple, n::Tuple) = Plane(Point(p), Vec(n))

paramdim(::Type{<:Plane}) = 2

isconvex(::Type{<:Plane}) = true

boundary(::Plane) = nothing

(p::Plane)(s, t) = p.p + s*p.u + t*p.v

"""
    origin(plane)

Return the origin of the `plane`.
"""
origin(p::Plane) = p.p

"""
    normal(plane)

Normal vector to the `plane`.
"""
normal(p::Plane) = p.u × p.v


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

struct AdminTeamsAdminsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminTeamsAdminsList(_api::AdminTeamsAdminsApi, team_id::String; limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.teams.admins.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List all of the admins on a given workspace.
Param: team_id::String (required)
Param: limit::Int32
Param: cursor::String
Return: Dict{String, Any}
"""
function adminTeamsAdminsList(_api::AdminTeamsAdminsApi, team_id::String; limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsAdminsList(_api, team_id; limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsAdminsList(_api::AdminTeamsAdminsApi, response_stream::Channel, team_id::String; limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsAdminsList(_api, team_id; limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminTeamsAdminsList


text:
# Code taken from FluxML/Optimisers.jl

const Param{T<:Number} = Union{AbstractArray{T},T}

_apply(opt, x, x̄, state) = apply(opt, x, x̄, state)
_apply(opt, x, x̄, ::Nothing) = apply(opt, x, x̄)

# Immutable updates
update!(opt, x, xb, state=state_init(opt, x)) = _update!(opt, x, xb, state)
update(opt, x, xb, state=state_init(opt, x))  = _update(opt, x, xb, state)

function _update(opt, x::Param, x̄::Param, state = nothing)
  Δ, state = _apply(opt, x, x̄, state)
  return x .- Δ, state
end

# Mutable updates

# Figure out if we can do in-place
inplace(x, y) = false
inplace(x, y::Nothing) = true
inplace(x::AbstractArray, x̄::AbstractArray) = true
inplace(x, x̄::NamedTuple) = all(inplace(getfield(x, f), getfield(x̄, f)) for f in fieldnames(typeof(x̄)))

function _update!(opt, x::AbstractArray{<:Number}, x̄::AbstractArray, state = nothing)
  Δ, state = _apply(opt, x, x̄, state)
  #println("fr = ", x)
  #println("to = ", Δ)
  x .-= Δ
  #println("after = ", x)
  return state
end

#function _update!(opt, x::Tuple, x̄::Tuple, state=nothing)
function _update!(opt, x, x̄::Tuple, state=nothing)
  for f in propertynames(x̄)
    f̄ = getindex(x̄, f)

    f̄ === nothing || _update!(opt, getindex(x, f), f̄, state)
  end
  return state
end

function _update!(opt, x, x̄::NamedTuple, state=nothing)
  for f in propertynames(x̄)
    f == :tuple_all_weights && continue
    f̄ = getproperty(x̄, f)

    f̄ === nothing || _update!(opt, getfield(x, f), f̄, state)
  end
  return state
end

function _update!(opt, x, x̄, state=nothing)
  for f in propertynames(x̄)
    if f isa Int
      f̄ = getindex(x̄, f)
      f̄ === nothing || _update!(opt, getindex(x, f), f̄, state)
    else
      f̄ = getproperty(x̄, f)
      f̄ === nothing || _update!(opt, getfield(x, f), f̄, state)
    end
  end
  return state
end


text:
orig_stdout = STDOUT
sr, sw = redirect_stdout()

println("hello")
println("world")

redirect_stdout(orig_stdout)

println(sr)