text:
module MultiplesOfPi

using FastRationals: FastRational

export PiExpTimes
export Pi

"""
	PiExpTimes(x, n)

Construct a number that behaves as `x*π^n` for a real `x` and an integer `n`

# Examples
```jldoctest
julia> PiExpTimes(2, 3)
2Pi^3
```

See also: [`Pi`](@ref)
"""
struct PiExpTimes{T<:Real,N<:Real} <: AbstractIrrational
	x :: T
	n :: N
end

PiExpTimes{T,N}(x::Real) where {T,N} = PiExpTimes{T,N}(T(x), zero(N))
PiExpTimes{T,N}(p::PiExpTimes{S}) where {T,S,N} = PiExpTimes{T,N}(T(p.x), N(p.n))

PiExpTimes{PiExpTimes{T}}(p::PiExpTimes{T}, n::Real) where {T<:Real} = PiExpTimes{PiExpTimes{T}, typeof(n)}(p, n)
PiExpTimes{PiExpTimes}(p::PiExpTimes, n::Real) = PiExpTimes{PiExpTimes, typeof(n)}(p, n)
PiExpTimes{Irrational{:π}}(p::Irrational{:π}, n::Real) = PiExpTimes{Irrational{:π}, typeof(n)}(p, n)

PiExpTimes{T}(x::Real, n::Real = 0) where {T<:Real} = PiExpTimes{T,typeof(n)}(T(x), n)
PiExpTimes{T}(::Irrational{:π}, n::Real) where {T<:Real} = PiExpTimes{T,typeof(n)}(T(1), n + 1)
PiExpTimes{T}(p::PiExpTimes, n::Real) where {T<:Real} = PiExpTimes{T,typeof(n)}(p.x, p.n + n)

PiExpTimes(::Irrational{:π}, n::Real = 0) = PiExpTimes{Int}(1, n + one(n))
# avoid nesting
PiExpTimes(p::PiExpTimes, n::Real = 0) = PiExpTimes(p.x, p.n + n)
PiExpTimes(x::Real) = PiExpTimes(x, 0)

"""
	Pi

The number `PiExpTimes(1, 1)`, numerically equivalent to `pi`.
Using `Pi` instead of `pi` often produces results that avoid floating-point inaccuracies.

# Examples
```jldoctest
julia> sin(Pi)
0.0

julia> sin(π)
1.2246467991473532e-16

julia> exp(im*Pi) + 1 == 0
true

julia> exp(im*π) + 1 == 0
false

julia> (2//3)Pi + Pi == (5//3)Pi
true

julia> (2//3)π + π == (5//3)π
false
```
"""
const Pi = PiExpTimes(1, 1)

function Base.:(<)(p::PiExpTimes, q::PiExpTimes)
	p.x == q.x ? p.n < q.n :
	p.n == q.n ? p.x < q.x :
	float(p) < float(q)
end
function Base.:(==)(p::T, q::T) where {T<:PiExpTimes}
	iszero(p.x) && iszero(q.x) ? true :
	(p.n == q.n && p.x == q.x) ? true :
	false
end

# Compare nested types by simplifying and promoting to a common type
function Base.:(==)(p::PiExpTimes, q::PiExpTimes)
	p2, q2 = promote(simplify(p), simplify(q))
	p2 == q2
end

Base.:(==)(p::PiExpTimes, x::Real) = p == PiExpTimes(x, 0)
Base.:(==)(x::Real, p::PiExpTimes) = p == x

simplify(p) = p
simplify(p::PiExpTimes{Irrational{:π}}) = PiExpTimes(1, p.n + one(p.n))
simplify(p::PiExpTimes{<:PiExpTimes}) = simplify(p.x) * Pi^(p.n)

for f in [:(<), :(==)]
	@eval Base.$f(p::PiExpTimes, ::Irrational{:π}) = $f(p, Pi)
	@eval Base.$f(::Irrational{:π}, p::PiExpTimes) = $f(Pi, p)
end

Base.promote_rule(::Type{Irrational{:π}}, ::Type{PiExpTimes{T,N}}) where {T,N} = PiExpTimes{T,N}
Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{Irrational{:π}}) where {T,N} = PiExpTimes{T,N}

Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{Bool}) where {T,N} = PiExpTimes{T,N}
Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{S}) where {T,S<:AbstractIrrational,N} = PiExpTimes{promote_type(S,T),N}
Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{S}) where {T,S<:Real,N} = PiExpTimes{promote_type(S,T),N}
Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{Float16}) where {T,N} = PiExpTimes{promote_type(Float16,T),N}
Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{Float32}) where {T,N} = PiExpTimes{promote_type(Float32,T),N}

Base.promote_rule(::Type{BigFloat}, ::Type{PiExpTimes{T,N}}) where {T,N} = PiExpTimes{promote_type(BigFloat,T),N}

Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{PiExpTimes{S,M}}) where {T,S,N,M} = PiExpTimes{promote_type(T, S),promote_type(M, N)}

Base.promote_rule(::Type{PiExpTimes{T,N}}, ::Type{Complex{S}}) where {T<:Real,S<:Real,N} = Complex{PiExpTimes{promote_type(T, S),N}}

Base.AbstractFloat(p::PiExpTimes{T,N}) where {T,N} = promote_type(T,N,Float64)(p)

for T in (:Float64,:Float32,:Float16)
	@eval begin
		function Base.$T(p::PiExpTimes)
			n = p.n
			if iszero(n)
				return $T(p.x)
			else
				return $T(p.x*float(pi)^n)
			end
		end
	end
end

# BigFloats need to be handled differently to gain accuracy
# Convert π to BigFloat before multiplying
function Base.BigFloat(p::PiExpTimes)
	n = p.n
	if iszero(n)
		return BigFloat(p.x)
	else
		return BigFloat(p.x*BigFloat(π)^n)
	end
end

Base.iszero(p::PiExpTimes{<:Real,Int}) = iszero(p.x)
Base.iszero(p::PiExpTimes) = iszero(float(p))
Base.isfinite(p::PiExpTimes) = isfinite(float(p))
Base.isnan(p::PiExpTimes) = isnan(float(p))

# Unfortunately Irrational numbers do not have a multiplicative identity of the same type,
# so we make do with something that works
Base.one(::Type{<:PiExpTimes}) = true
Base.isone(p::PiExpTimes) = (x = simplify(p); isone(x.x) && iszero(x.n))

Base.zero(p::PiExpTimes{T,N}) where {T,N} = PiExpTimes{T,N}(zero(T), p.n)
Base.zero(::Type{PiExpTimes{T,N}}) where {T,N} = PiExpTimes{T,N}(zero(T), zero(N))

Base.sign(p::PiExpTimes) = sign(p.x)
Base.signbit(p::PiExpTimes) = signbit(p.x)

# Define trigonometric functions

# this is type-unstable, but Union splitting does its magic
@inline dividebypiexp(p::PiExpTimes, n) = p.n == n ? p.x : p.x * π^float(p.n-n)
@inline dividebypi(p) = dividebypiexp(p, 1)

@inline _check_exp_apply(f, fzeroexp, p) = iszero(p.n) ? fzeroexp(p.x) : f(dividebypi(p))
Base.sin(p::PiExpTimes) = _check_exp_apply(sinpi, sin, p)
# Base.sin(p::PiExpTimes) = iszero(p.n) ? sin(p.x) : sinpi(dividebypi(p))
Base.cos(p::PiExpTimes) = _check_exp_apply(cospi, cos, p)
if VERSION > v"1.6.0"
	Base.sincos(p::PiExpTimes) = _check_exp_apply(sincospi, sincos, p)
else
	Base.sincos(p::PiExpTimes) = (sin(p),cos(p))
end

function Base.cis(p::PiExpTimes)
	s, c = sincos(p);
	Complex(c, s)
end

@inline function Base.tan(p::PiExpTimes{T}) where {T}
	Tf = float(T)
	if isone(simplify(p).n)
		r = abs(rem(p.x, one(p.x)))
		iszero(r) && return copysign(zero(Tf), p.x)
		2r == one(r) && return copysign(Tf(Inf), p.x)
		4r == one(r) && return copysign(one(Tf), p.x)
	end
	tan(float(p))
end

function Base.tanh(z::Complex{<:PiExpTimes})
	iszero(real(z)) && return im*tan(imag(z))
	tanh(float(z))
end

Base.sinc(p::PiExpTimes) = sinc(float(p))

# Arithmetic operators

function match_exponent(p1, p2)
	minexp = min(p1.n, p2.n)
	x = dividebypiexp(p1, minexp)
	y = dividebypiexp(p2, minexp)
	x, y, minexp
end

function Base.:(+)(p1::PiExpTimes, p2::PiExpTimes)
	x, y, exponent = match_exponent(p1, p2)
	PiExpTimes(x + y, exponent)
end

Base.:(-)(p::PiExpTimes) = PiExpTimes(-p.x, p.n)

function Base.:(-)(p1::PiExpTimes, p2::PiExpTimes)
	x, y, exponent = match_exponent(p1, p2)
	PiExpTimes(x - y, exponent)
end

Base.inv(p::PiExpTimes) = PiExpTimes(inv(p.x), -p.n)

function Base.:(*)(p1::PiExpTimes, p2::PiExpTimes)
	PiExpTimes(p1.x*p2.x, p1.n+p2.n)
end
function Base.:(/)(p1::PiExpTimes, p2::PiExpTimes)
	PiExpTimes(p1.x/p2.x, p1.n-p2.n)
end

for f in [:(+), :(-), :(*), :(/)]
	@eval Base.$f(p::PiExpTimes, ::Irrational{:π}) = $f(p, Pi)
	@eval Base.$f(::Irrational{:π}, p::PiExpTimes) = $f(Pi, p)
end

Base.:(*)(b::Bool, p::PiExpTimes) = ifelse(b, p, flipsign(zero(p), p.x))
Base.:(/)(b::Bool, p::PiExpTimes) = PiExpTimes(b / p.x, -p.n)

Base.:(^)(p::PiExpTimes, q::PiExpTimes) = PiExpTimes(float(p.x)^oftype(float(p.x), q), q*p.n)
Base.:(^)(::Irrational{:π}, p::PiExpTimes) = Pi^p
Base.:(^)(p::PiExpTimes, ::Irrational{:π}) = PiExpTimes(p.x^π, Pi*p.n)
_default_exp(p, n) = PiExpTimes(p.x^n, n*p.n)
# These methods need to be defined to handle ambiguities
Base.:(^)(p::PiExpTimes, n::AbstractIrrational) = _default_exp(p, n)
Base.:(^)(p::PiExpTimes, n::Rational) = _default_exp(p, n)
Base.:(^)(p::PiExpTimes, n::Integer) = _default_exp(p, n)
Base.:(^)(p::PiExpTimes, n::Real) = _default_exp(p, n)

# Rational divide
Base.:(//)(p::PiExpTimes, q::PiExpTimes) = PiExpTimes(p.x//q.x, p.n - q.n)
Base.:(//)(p::PiExpTimes, y::Real) = PiExpTimes(p.x//y, p.n)
Base.:(//)(p::PiExpTimes, ::Irrational{:π}) = PiExpTimes(p.x//1, p.n-1)
Base.:(//)(y::Real, p::PiExpTimes) = PiExpTimes(y//p.x, -p.n)
Base.:(//)(::Irrational{:π}, p::PiExpTimes) = PiExpTimes(1//p.x, -p.n+1)

_half(n::Integer) = FastRational(n, 2)
_half(n::Rational) = n//2
_half(n) = n/2
_third(n::Integer) = FastRational(n, 3)
_third(n::Rational) = n//3
_third(n) = n/3

Base.sqrt(p::PiExpTimes) = PiExpTimes(sqrt(p.x), _half(p.n))
Base.cbrt(p::PiExpTimes) = PiExpTimes(cbrt(p.x), _third(p.n))

# Pretty-printing

function Base.show(io::IO, p::PiExpTimes)
	if !isfinite(p)
		print(io, float(p))
		return
	end
	x = p.x

	expstr(p) = isone(p.n) ? "Pi" : "Pi^"*string(p.n)

	tostr(x,p) = isone(x) ? expstr(p) : string(x)*"*"*expstr(p)
	tostr(x::Integer,p) = isone(x) ? expstr(p) : string(x)*expstr(p)
	tostr(x::AbstractFloat,p) = isone(x) ? expstr(p) :
								isinf(x) || isnan(x) ? string(x) :
								string(x)*"*"*expstr(p)
	tostr(x::Rational,p) = "("*string(x)*")"*expstr(p)

	str = iszero(x) ? string(x) : tostr(x,p)
	print(io, str)
end

end # module


text:
######################################################################
# DeterminantalPointProcesses.jl
# Determinantal Point Processes in Julia
# http://github.com/theogf/DeterminantalPointProcesses.jl
# MIT Licensed
######################################################################

__precompile__(true)

module DeterminantalPointProcesses

using Distributed
using Distributions: Distributions, pdf, logpdf
using LinearAlgebra
using Random: Random, rand, bitrand, AbstractRNG, MersenneTwister, GLOBAL_RNG
using Requires
using SharedArrays
import Base: rand

export
    # point process types and aliases
    DeterminantalPointProcess,
    DPP,
    kDeterminantalPointProcess,
    kDPP,

    # mehtods
    logpdf,             # log probability mass function
    pdf,                # probability mass function
    rand,               # generate samples
    randmcmc            # generate samples using MCMC

function __init__()
    @require KernelFunctions="ec8451be-7e33-11e9-00cf-bbf324bd1392" include("kernelcompat.jl")
end
### source files

# Types
include("types.jl")

# pdf and rand methods
include("prob.jl")
include("rand.jl")

# utilities
include("utils.jl")

end # module

ext: jl


text:
# # [Getting Started: Optimize!](@id Optimize)
#
# This example illustrates how to set up and solve optimization problems and how
# to further get data from the algorithm using [`DebugOptions`](@ref) and
# [`RecordOptions`](@ref)
#
# To start from the quite general case: A __Solver__ is an algorithm that aims
# to solve
#
# $\operatorname*{argmin}_{x∈\mathcal M} f(x)$
#
# where $\mathcal M$ is a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) and
# $f\colon\mathcal M \to \mathbb R$ is the cost function.
#
# In `Manopt.jl` a __Solver__ is an algorithm that requires a [`Problem`](@ref)
# `p` and [`Options`](@ref) `o`. While former contains __static__ data,
# most prominently the manifold $\mathcal M$ (usually as `p.M`) and the cost
# function $f$ (usually as `p.cost`), the latter contains __dynamic__
# data, i.e. things that usually change during the algorithm, are allowed to
# change, or specify the details of the algorithm to use. Together they form a
# `plan`. A `plan` uniquely determines the algorithm to use and provide all
# necessary information to run the algorithm.
#
# ## Example
#
# A gradient plan consists of a [`GradientProblem`](@ref) with the fields `M`,
# cost function $f$ as well as `gradient` storing the gradient function
# corresponding to $f$. Accessing both functions can be done directly but should
# be encapsulated using [`get_cost`](@ref)`(p,x)` and [`get_gradient`](@ref)`(p,x)`,
# where in both cases `x` is a point on the [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`.
# Second, the [`GradientDescentOptions`](@ref) specify that the algorithm to solve
# the [`GradientProblem`](@ref) will be the [gradient
# descent](https://en.wikipedia.org/wiki/Gradient_descent) algorithm. It requires
# an initial value `o.x0`, a [`StoppingCriterion`](@ref) `o.stop`, a
# [`Stepsize`](@ref) `o.stepsize` and a retraction `o.retraction` and it
# internally stores the last evaluation of the gradient at `o.∇` for convenience.
# The only mandatory parameter is the initial value `x0`, though the defaults for
# both the stopping criterion ([`StopAfterIteration`](@ref)`(100)`) as well as the
# stepsize ([`ConstantStepsize`](@ref)`(1.)` are quite conservative, but are
# chosen to be as simple as possible.
#
# With these two at hand, running the algorithm just requires to call `xOpt = solve(p,o)`.
#
# In the following two examples we will see, how to use a higher level interface
# that allows to more easily activate for example a debug output or record values during the iterations
#
# ## The given Dataset
#
exportFolder = joinpath(@__DIR__,"..","..","docs","src","assets","images","tutorials") #src
using Manopt, Manifolds
using Random, Colors
# For a persistent random set we use
n = 100
σ = π/8
M = Sphere(2)
x = 1/sqrt(2)*[1., 0., 1.]
Random.seed!(42)
data = [exp(M,x,random_tangent(M,x,Val(:Gaussian),σ)) for i ∈ 1:n ]
nothing #hide
# and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)
black = RGBA{Float64}(colorant"#000000")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
nothing #hide
#
# Then our data looks like
#
asymptote_export_S2_signals(exportFolder*"/startDataAndCenter.asy"; #src
    points = [ [x], data], #src
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]), #src
    dotSize = 3.5, cameraPosition = (1.,.5,.5) #src
) #src
render_asymptote(exportFolder*"/startDataAndCenter.asy"; render = 2) #src
#md # ```julia
#md # asymptote_export_S2_signals("startDataAndCenter.asy";
#md #     points = [ [x], data],
#md #     colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
#md #     dotSize = 3.5, cameraPosition = (1.,.5,.5)
#md # )
#md # render_asymptote("startDataAndCenter.asy"; render = 2)
#md # ```
#md #
#md # ![The data of noisy versions of $x$](../assets/images/tutorials/startDataAndCenter.png)
#
# ## Computing the Mean
#
# To compute the mean on the manifold we use the characterization, that the
# Euclidean mean minimizes the sum of squared distances, and end up with the
# following cost function. Its minimizer is called
# [Riemannian Center of Mass](https://arxiv.org/abs/1407.2087).
#
F = y -> sum(1/(2*n) * distance.(Ref(M),Ref(y),data).^2)
∇F = y -> sum(1/n*∇distance.(Ref(M),data,Ref(y)))
nothing #hide
#
# note that the [`∇distance`](@ref) defaults to the case `p=2`, i.e. the
# gradient of the squared distance. For details on convergence of the gradient
# descent for this problem, see [[Afsari, Tron, Vidal, 2013](#AfsariTronVidal2013)]
#
# The easiest way to call the gradient descent is now to call
# [`steepest_descent`](@ref)
xMean = steepest_descent(M,F,∇F,data[1])
nothing; #hide
# but in order to get more details, we further add the `debug=` options, which
# act as a [decorator pattern](https://en.wikipedia.org/wiki/Decorator_pattern)
# using the [`DebugOptions`](@ref) and [`DebugAction`](@ref)s. The latter store
# values if that's necessary, for example for the [`DebugChange`](@ref) that prints
# the change during the last iteration. The following debug prints
#
# `# i | x: | Last Change: | F(x): ``
#
# as well as the reason why the algorithm stopped at the end.
# Here, the format shorthand and the [`DebugFactory`] are used, which returns a
# [`DebugGroup`](@ref) of [`DebugAction`](@ref) performed each iteration and the stop,
# respectively.
xMean = steepest_descent(M,F,∇F,data[1];
   debug = [:Iteration," | ", :x, " | ", :Change, " | ", :Cost, "\n", :Stop]
)
nothing #hide
#
asymptote_export_S2_signals(exportFolder*"/startDataCenterMean.asy"; #src
    points = [ [x], data, [xMean] ], #src
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange]), #src
    dotSize = 3.5, cameraPosition = (1.,.5,.5) #src
) #src
render_asymptote(exportFolder*"/startDataCenterMean.asy"; render = 2) #src
#md # ```julia
#md # asymptote_export_S2_signals("startDataCenterMean.asy";
#md #     points = [ [x], data, [xMean] ],
#md #     colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange]),
#md #     dotSize = 3.5, cameraPosition = (1.,.5,.5)
#md # )
#md # render_asymptote("startDataCenterMean.asy"; render = 2)
#md # ```
#md #
#md # ![The resulting mean (orange)](../assets/images/tutorials/startDataCenterMean.png)
#
# ## Computing the Median
#
# Similar to the mean you can also define the median as the minimizer of the
# distances, see for example [[Bačák, 2014](#Bačák2014)], but since
# this problem is not differentiable, we employ the Cyclic Proximal Point (CPP)
# algorithm, described in the same reference. We define
#
F2 = y -> sum( 1/(2*n) * distance.(Ref(M),Ref(y),data))
proxes = Function[ (λ,y) -> prox_distance(M,λ/n,di,y,1) for di in data ]
nothing #hide
# where the `Function` is a helper for global scope to infer the correct type.
#
# We then call the [`cyclic_proximal_point`](@ref) as
o = cyclic_proximal_point(M,F2,proxes,data[1];
    debug = [:Iteration," | ", :x, " | ", :Change, " | ", :Cost, "\n", 50, :Stop],
    record = [:Iteration, :Change, :Cost],
    return_options = true
)
xMedian = get_solver_result(o)
values = get_record(o)
nothing # hide
# where the differences to [`steepest_descent`](@ref) are as follows
#
# * the third parameter is now an Array of proximal maps
# * debug is reduces to only every 50th iteration
# * we further activated a [`RecordAction`](@ref) using the `record=` optional
#   parameter. These work very similar to those in debug, but they
#   collect their data in an array. The high level interface then returns two
#   variables; the `values` do contain an array of recorded
#   datum per iteration. Here a Tuple containing the iteration, last change and
#   cost respectively; see [`RecordGroup`](@ref), [`RecordIteration`](@ref),
#   [`RecordChange`](@ref), [`RecordCost`](@ref) as well as the [`RecordFactory`](@ref)
#   for details. The `values` contains hence a tuple per iteration,
#   that itself consists of (by order of specification) the iteration number,
#   the last change and the cost function value.
#
# These recorded entries read
values
# The resulting median and mean for the data hence are
#
asymptote_export_S2_signals(exportFolder*"/startDataCenterMedianAndMean.asy"; #src
    points = [ [x], data, [xMean], [xMedian] ], #src
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange, TolVibrantMagenta]), #src
    dotSize = 3.5, cameraPosition = (1.,.5,.5) #src
) #src
render_asymptote(exportFolder*"/startDataCenterMedianAndMean.asy"; render = 2) #src
#md # ```julia
#md # asymptote_export_S2_signals("startDataCenterMean.asy";
#md #     points = [ [x], data, [xMean], [xMedian] ],
#md #     colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange, TolVibrantMagenta]),
#md #     dotSize = 3.5, cameraPosition = (1.,.5,.5)
#md # )
#md # render_asymptote("startDataCenterMedianAndMean.asy"; render = 2)
#md # ```
#md #
#md # ![The resulting mean (orange) and median (magenta)](../assets/images/tutorials/startDataCenterMedianAndMean.png)
#
# ## Literature
#
# ```@raw html
# <ul>
# <li id="Bačák2014">[<a>Bačák, 2014</a>]
#   Bačák, M: <emph>Computing Medians and Means in Hadamard Spaces.</emph>,
#   SIAM Journal on Optimization, Volume 24, Number 3, pp. 1542–1566,
#   doi: <a href="https://doi.org/10.1137/140953393">10.1137/140953393</a>,
#   arxiv: <a href="https://arxiv.org/abs/1210.2145">1210.2145</a>.</li>
#   <li id="AfsariTronVidal2013">[<a>Afsari, Tron, Vidal, 2013</a>]
#    Afsari, B; Tron, R.; Vidal, R.: <emph>On the Convergence of Gradient
#    Descent for Finding the Riemannian Center of Mass</emph>,
#    SIAM Journal on Control and Optimization, Volume 51, Issue 3,
#    pp. 2230–2260.
#    doi: <a href="https://doi.org/10.1137/12086282X">10.1137/12086282X</a>,
#    arxiv: <a href="https://arxiv.org/abs/1201.0925">1201.0925</a></li>
# </ul>
# ```


text:
module IsPower

using Hecke, InteractiveUtils
import Nemo

function Hecke.roots_upper_bound(f::fmpz_poly)
  a = coeff(f, degree(f))
  return max(fmpz(1), maximum([ceil(fmpz, abs(coeff(f, i)//a)) for i=0:degree(f)]))
end

function Hecke.roots_upper_bound(f::fmpq_poly)
  a = coeff(f, degree(f))
  return max(fmpq(1), maximum([abs(coeff(f, i)//a) for i=0:degree(f)]))
end

function ispower_mod_p(a::nf_elem, i::Int)

  p = 2^10
  K = parent(a)
  f = K.pol
  first = true
  cnt = 0
  local opt
  while true
    p = next_prime(p)
    if gcd(p-1, i) > 1
      continue
    end
    lp = factor(f, GF(p))
    if any(x->x>1, values(lp.fac))
      continue
    end
    if any(x->degree(x) > 100, keys(lp.fac))
      @show p, "deg", maximum(degree, keys(lp.fac))
      continue
    end
    if length(keys(lp.fac)) > 30
      @show p, "len", length(keys(lp.fac))
      continue
    end
    if first
      opt = (p, length(lp))
      first = false
    else
      if length(lp) < opt[2]
        opt = (p, length(lp))
        cnt = 0
      else
        cnt += 1
      end
    end
    if cnt > 2*i
      break
    end
  end
  println("using ", opt[1], " with ", opt[2], " local factors")

  C = Hecke.qAdicConj(K, opt[1])
  p = opt[1]

  local power_sum
  if false
    r_pr = 20
    local con_r
    while true
      con_r = map(abs, conjugates(a, r_pr))
      map(contains_zero, con_r)
      if any(x->contains_zero(x), con_r)
        r_pr += 20
      else
        con_r = [root(x, i) for x = con_r]
        break
      end
    end
    power_sum = l-> Hecke.upper_bound(sum(x^l for x = con_r), fmpz)
  else
    B = ceil(fmpz, roots_upper_bound(parent(a).pol))
    c = fmpz(0)
    for j=a.elem_length-1:-1:0
      c = B*c+ceil(fmpz, abs(coeff(a, j)))
    end
    c = root(c, i)
    power_sum = l-> degree(parent(a))*c^l
  end

  @show bd = map(power_sum, 1:15)
  pr = clog(bd[end], p) + 145
  println("using a precision of ", pr)
  con_p = conjugates(a, C, pr,all = false, flat = false)
  @assert all(x->degree(minpoly(ResidueField(parent(x))[2](x))) == degree(parent(x)), con_p)

  con_pr = [roots(x, i) for x = con_p] #select primes to minimize this

  @show "#roots per local factor"
  @show map(length, con_pr)
  #TODO: if there is only one root/local => lll not neccessary,
  #      go directly to CRT
  con_pr_j = [[one(parent(x)) for x = y] for y = con_pr]
  no_fac = sum(map(length, con_pr))
  j = 0
  trafo = identity_matrix(FlintZZ, no_fac)
  no_rt = no_fac
  while true
    j += 1
    if j >= length(bd)
      bd = vcat(bd, map(power_sum, length(bd)+1:length(bd)+5))
    end

    con_pr_j = [ con_pr[i] .* con_pr_j[i] for i=1:length(con_pr)]
    @assert con_pr_j[1][1] == con_pr[1][1]^j

    data = matrix([reduce(vcat, [map(x -> lift(trace(x)), y) for y = con_pr_j])])
    data = sub(trafo, 1:no_rt, 1:no_fac)*data
    k = clog(bd[j], p)
    @assert k < pr
    pk = fmpz(p)^(pr-k)
    map_entries!(x->rem(div(x, fmpz(p)^(k)), pk), data, data)
    if iszero(data)
      println("nothing new...")
      continue
    end
    trafo = hcat(trafo, data)
    data = zero_matrix(FlintZZ, 1, ncols(trafo))
    data[1, end] = pk
    trafo = vcat(trafo, data)
    #= roots: products of the local roots that ar small
       corresponds to sums of local power-sums that are small
       the "correct" one will have no_fac ones in the front
         rear  (a+b) < bd => (a+b)/c < bd/c
         round((a+b)/c) - round(a/c) - round(b/c)
    =#
    no_rt, trafo = lll_with_removal(trafo, fmpz(p)^2*fmpz(2*no_fac)) #THINK
    trafo = sub(trafo, 1:no_rt, 1:ncols(trafo))
    d = Dict{fmpz_mat, Array{Int, 1}}()
    for l=1:no_fac
      k = trafo[:, l]
      if haskey(d, k)
        push!(d[k], l)
      else
        d[k] = [l]
      end
    end
    println("partition: ", values(d))
    if all(x->length(x) >= length(con_pr), values(d))
      # we have a candidate!
      @show "success"
      m = Base.minimum(map(length, values(d)))
      pa = d[findfirst(x->length(x) == m, d)]
      co = reduce(vcat, con_pr)
      mp = Dict{Int, Int}()
      i = 1
      s = 1
      for c = con_pr
        for x = c
          mp[s] = i
          s += 1
        end
        i += 1
      end
      @show mp, pa
      mC = completions(C)
      @show [mp[x] for x = pa]
      mo = factors(C)
      va = zeros(K, length(mo))
      for x = pa
        y = Hecke.mod_sym(preimage(mC[mp[x]], co[x]), fmpz(p)^pr)
        va[mp[x]] = y
      end
#      va = [Hecke.mod_sym(preimage(mC[mp[x]], co[x]), fmpz(p)^pr) for x = pa]
      @assert length(mo) == length(va)
      pk = fmpz(C.C.H.p)^Int(C.C.H.prev)
      res = crt(map(Hecke.Globals.Zx, va), mo, C.C.H, pk)
      return res(gen(K)), pk
    end
    if j > 40 error("") end
  end
end

function Hecke.crt(A::Array{<:NfAbsOrdElem, 1}, I::Array{<:NfAbsOrdIdl, 1})
  while length(I) > 1
    II = typeof(I[1])[]
    AA = typeof(A[1])[]
    i = 1
    while i< length(I)
      push!(II, I[i]*I[i+1])
      push!(AA, crt(A[i], I[i], A[i+1], I[i+1]))
      i += 2
    end
    if i == length(I)
      push!(II, I[end])
      push!(AA, A[end])
    end
    I = II
    A = AA
  end
end

function Base.lastindex(M::MatElem, i::Int)
  i == 1 && return nrows(M)
  i == 2 && return ncols(M)
  error("illegal dimensino")
end

function Hecke.ispower(a::qadic, i::Int)
  @assert i>0
  if i==1
    return true, a
  end

  k, mk = ResidueField(parent(a))
  v = valuation(a)
  if v % i != 0
    return false, a
  end
  rv = divexact(v, i)
  a = divexact(a, uniformizer(parent(a))^v)

  fl, bp = ispower(mk(a), i)
  if !fl
    return false, a
  end

  b = preimage(mk, bp)

  p = precision(a)
  while p > 0
    b = b - (b^i-a)//(i*b^(i-1))
    p = div(p, 2)
  end
  return true, b*uniformizer(parent(a))^rv
end

function Hecke.roots(a::qadic, i::Int)
  @assert i>0
  if i==1
    return a
  end

  k, mk = ResidueField(parent(a))
  v = valuation(a)
  if v % i != 0
    error("elem not an $i-th power")
  end
  rv = divexact(v, i)
  a = divexact(a, uniformizer(parent(a))^v)

  fl, bp = ispower(mk(a), i)
  if !fl
    error("elem not an $i-th power")
  end

  local zeta
  j = gcd(size(k)-1, i)
  if j != 1
    while true
      zeta = rand(k)^(divexact(size(k)-1, j))
      if all(l-> !isone(zeta^l), 1:j-1)
        break
      end
    end
  else
    zeta = k(1)
  end

  b = preimage(mk, bp)
  z = preimage(mk, zeta)

  p = precision(a)
  while p > 0
    b = b - (b^i-a)//(i*b^(i-1))
    if !isone(zeta)
      z = z - (z^j-one(parent(z)))//(j*z^(j-1))
    end
    p = div(p, 2)
  end
  b *= uniformizer(parent(a))^rv
  if !isone(zeta)
    return [b*z^l for l=0:Int(j)-1]
  else
    return [b]
  end
end

function completions(C::qAdicConj)
  mC = []
  K = C.K
  a = gen(K)
  z = conjugates(a, C, 5, all = false)
  for t = z
    _, x = completion(K, t)
    push!(mC, x)
  end
  return mC
end

function factors(C::qAdicConj)
  return [x[1] for x = Hecke.factor_mod_pk(Array, C.C.H, Int(C.C.H.N))]
end

function Base.getindex(H::Hecke.HenselCtx, s::Symbol, i::Int)
  r = H.r
  @assert 1<= i<=2*r-2
  if s == :l
    return unsafe_load(H.link, i)
  elseif s == :v
    f = Hecke.Globals.Zx()
    @GC.preserve f H ccall((:fmpz_poly_set, Nemo.libflint), Cvoid, (Ref{fmpz_poly}, Ptr{Hecke.fmpz_poly_raw}), f, H.v+(i-1)*sizeof(Hecke.fmpz_poly_raw))
    return f
  elseif s == :w
    f = Hecke.Globals.Zx()
    @GC.preserve f H ccall((:fmpz_poly_set, Nemo.libflint), Cvoid, (Ref{fmpz_poly}, Ptr{Hecke.fmpz_poly_raw}), f, H.w+(i-1)*sizeof(Hecke.fmpz_poly_raw))
    return f
  end
end

function Hecke.mod_sym!(f::fmpz_poly, n::fmpz)
  for i=degree(f):-1:0
    setcoeff!(f, i, Hecke.mod_sym(coeff(f, i), n))
  end
  return f
end

function Hecke.crt(v::Array{fmpz_poly, 1}, m::Array{fmpz_poly, 1}, H::Hecke.HenselCtx, pk::fmpz = fmpz(H.p)^Int(H.prev))
  if length(v) == 1
    return v[1]
  end
  p = fmpz(H.p)
  r = H.r
  res = fmpz_poly[]
  for i = 1:2*r-2
    @show mu = H[:l, i]
    if mu < 0
      push!(res, divrem(v[-mu], H[:v, -mu])[2])
      mod_sym!(res[end], pk)
      continue
    end
    f = H[:v, mu+1]
    g = H[:v, mu+2]
    t = f*H[:w, mu+1] + g*H[:w, mu+2] - 1
#    @show map(x->valuation(coeff(t, x), p), 0:degree(t))
#    @assert all(x->valuation(coeff(t, x), p) > 1, 0:degree(t))
    push!(res, res[mu+2]*f*H[:w, mu+1] + res[mu+1]*g*H[:w, mu+2])
    mod_sym!(res[end], pk)
    res[end] = divrem(res[end], H[:v, length(res)])[2]
    mod_sym!(res[end], pk)
#    t = divrem(res[end], f)[2] - res[mu+1];
#    @show map(x->valuation(coeff(t, x), p), 0:degree(t))
#
#    t = divrem(res[end], g)[2] - res[mu+2];
#    @show map(x->valuation(coeff(t, x), p), 0:degree(t))
  end
  mu = 2*r-4
  f = H[:v, mu+1]
  g = H[:v, mu+2]
#    t = f*H[:w, mu+1] + g*H[:w, mu+2] - 1
#    @show map(x->valuation(coeff(t, x), p), 0:degree(t))
  push!(res, res[mu+2]*f*H[:w, mu+1] + res[mu+1]*g*H[:w, mu+2])
  mod_sym!(res[end], pk)

  f = deepcopy(H.f)
  f.parent = res[end].parent
  res[end] = divrem(res[end], f)[2]
  mod_sym!(res[end], pk)

  return res[end]
end

function Hecke.degree(R::qadic) #TODO XXXX
  return R.length
end

function Hecke.ideal(R::NfAbsOrd, f::MapFromFunc{AnticNumberField,FlintQadicField})
  K = nf(R)
  Cq = codomain(f)
  if degree(Cq) == degree(K)
    return prime(Cq)^precision(Cq)*R
  end
  p = prime(Cq)
  C = qAdicConj(K, Int(p))
  r = roots(C.C, 5)
  fa = f(gen(K))
  k = precision(fa)
  lf = Hecke.factor_mod_pk(C.C.H, k)
  for fp = keys(lf)
    kfa = fp(setprecision(fa, 2))
    if iszero(kfa)
      s = R(fp(gen(K)))
      if valuation(norm(s), p) > k
        s += p^k
      end
      i = ideal(R, p^k, s)
      i.gens_normal = p
      i.norm = (p^k)^degree(Cq)
      i.minimum = p^k
      return i
    end
  end
end

end


text:
Base.@kwdef struct LucasKanade
    iterations::Int64 = 30
    window_size::Int64 = 9
    pyramid_levels::Int64 = 3
    eigenvalue_threshold::Float64 = 1e-4
    ϵ::Float64 = 1e-2
end

function optflow!(
    displacement, first_pyramid, second_pyramid, points, algorithm::LucasKanade,
)
    has_enough_layers =
        length(first_pyramid.layers) > algorithm.pyramid_levels &&
        length(second_pyramid.layers) > algorithm.pyramid_levels
    has_enough_layers || throw("Not enough layers in pyramids.")

    n_points = length(points)
    status = trues(n_points)
    n_good = n_points

    window = algorithm.window_size
    mode = BSpline(Linear())

    for level in (algorithm.pyramid_levels + 1):-1:1
        level_resolution = axes(first_pyramid.layers[level])
        # Interpolate layer to get sub-pixel precision.
        # We never go out-of-bound, so there is no need to extrapolate.
        # We can safely use `interpolate!` function (thus saving some memory), because
        # when doing linear interpolation, there is no modification to the input array.
        # See: https://github.com/JuliaMath/Interpolations.jl/blob/master/src/b-splines/prefiltering.jl#L33
        interploated_layer = interpolate!(@inbounds(second_pyramid.layers[level]), mode)

        Threads.@threads for n in 1:n_points
            @inbounds status[n] || continue

            point = get_pyramid_coordinate(@inbounds(points[n]), level)
            offsets = get_offsets(point, point, window, level_resolution)
            grid = get_grid(point, offsets)

            G_inv, min_eigenvalue = compute_spatial_gradient(first_pyramid, grid, level)
            if min_eigenvalue < algorithm.eigenvalue_threshold
                @inbounds status[n] = false
                n_good -= 1
                continue
            end

            pyramid_contribution = SVector{2, Float64}(0.0, 0.0)
            for _ in 1:algorithm.iterations
                putative_flow = @inbounds displacement[n] + pyramid_contribution
                putative_correspondence = point + putative_flow
                if !lies_in(level_resolution, putative_correspondence)
                    @inbounds status[n] = false
                    n_good -= 1
                    break
                end

                new_offsets = get_offsets(
                    point, putative_correspondence, window, level_resolution)
                # Recalculate gradient only if the offset changes.
                if new_offsets != offsets
                    offsets = new_offsets
                    grid = get_grid(point, offsets)
                    G_inv, min_eigenvalue = compute_spatial_gradient(
                        first_pyramid, grid, level)
                    if min_eigenvalue < algorithm.eigenvalue_threshold
                        @inbounds status[n] = false
                        n_good -= 1
                        break
                    end
                else
                    offsets = new_offsets
                    grid = get_grid(point, offsets)
                end

                estimated_flow = compute_flow_vector(
                    putative_correspondence,
                    first_pyramid, interploated_layer, level,
                    grid, offsets, G_inv)

                # Epsilon termination criteria.
                abs(estimated_flow[1]) < algorithm.ϵ &&
                    abs(estimated_flow[2]) < algorithm.ϵ && break

                pyramid_contribution += estimated_flow
                in_bounds = lies_in(
                    level_resolution, putative_correspondence + estimated_flow)
                if !in_bounds
                    @inbounds status[n] = false
                    n_good -= 1
                    break
                end
            end
            @inbounds status[n] || continue
            @inbounds displacement[n] += pyramid_contribution
            level > 1 && (@inbounds displacement[n] *= 2.0;)
        end
    end

    displacement, status, n_good
end

function compute_partial_derivatives(Iy, Ix; kwargs...)
    Iyy = typeof(Iy)(undef, size(Iy))
    Ixx = typeof(Iy)(undef, size(Iy))
    Iyx = typeof(Iy)(undef, size(Iy))
    compute_partial_derivatives!(Iyy, Ixx, Iyx, Iy, Ix; kwargs...)
end

function compute_partial_derivatives!(
    Iyy, Ixx, Iyx, Iy, Ix;
    squared = typeof(Iy)(undef, size(Iy)),
    filtered = typeof(Iy)(undef, size(Iy)), σ = 4.0,
)
    kernel_factors = get_kernel(σ, 2)

    squared .= Iy .* Iy
    imfilter!(filtered, squared, kernel_factors)
    integral_image!(Iyy, filtered)

    squared .= Ix .* Ix
    imfilter!(filtered, squared, kernel_factors)
    integral_image!(Ixx, filtered)

    squared .= Iy .* Ix
    imfilter!(filtered, squared, kernel_factors)
    integral_image!(Iyx, filtered)

    Iyy, Ixx, Iyx
end

function integral_image!(integral, img)
    sd = coords_spatial(img)
    cumsum!(integral, img; dims=sd[1])
    for i ∈ 2:length(sd)
        cumsum!(integral, integral; dims=sd[i])
    end
    integral
end

function _compute_spatial_gradient(
    grid, Iyy_integral, Iyx_integral, Ixx_integral,
)
    sum_Iyy = boxdiff(Iyy_integral, grid[1], grid[2])
    sum_Ixx = boxdiff(Ixx_integral, grid[1], grid[2])
    sum_Iyx = boxdiff(Iyx_integral, grid[1], grid[2])
    SMatrix{2, 2, Float64}(sum_Iyy, sum_Iyx, sum_Iyx, sum_Ixx)
end

function compute_spatial_gradient(pyramid::LKPyramid, grid, level)
    G = _compute_spatial_gradient(
        grid, pyramid.Iyy[level], pyramid.Iyx[level], pyramid.Ixx[level])
    U, S, V = svd2x2(G)
    G_inv = pinv2x2(U, S, V)
    min_eigenvalue = min(S[1, 1], S[2, 2]) / prod(length.(grid))

    G_inv, min_eigenvalue
end

function prepare_linear_system(corresponding_point, A, Iy, Ix, offsets, B)
    P, Q = size(A)
    by, bx = 0.0, 0.0

    @inbounds for q in 1:Q, p in 1:P
        r = corresponding_point[1] + offsets[1][p]
        c = corresponding_point[2] + offsets[2][q]

        δI = A[p, q] - B(r, c)
        by += δI * Iy[p, q]
        bx += δI * Ix[p, q]
    end

    SVector{2, Float64}(by, bx)
end

function compute_flow_vector(
    corresponding_point,
    first_pyramid::LKPyramid, etp, level,
    grid, offsets, G_inv,
)
    b = prepare_linear_system(
        corresponding_point,
        view(first_pyramid.layers[level], grid[1], grid[2]),
        view(first_pyramid.Iy[level], grid[1], grid[2]),
        view(first_pyramid.Ix[level], grid[1], grid[2]),
        offsets, etp)
    G_inv * b
end

@inline lies_in(area, point) = @inbounds(
    first(area[1]) ≤ point[1] ≤ last(area[1]) &&
    first(area[2]) ≤ point[2] ≤ last(area[2]))

"""
# Arguments
- `level`: Level of the pyramid in `[1, levels]` range.
"""
@inline get_pyramid_coordinate(point, level) = floor.(Int, point ./ 2 ^ (level - 1))

function get_offsets(point, new_point, window, image_axes)
    rows, cols = image_axes

    up = floor(Int, min(window, min(point[1], new_point[1]) - first(rows)))
    down = floor(Int, min(window, last(rows) - max(point[1], new_point[1])))
    left = floor(Int, min(window, min(point[2], new_point[2]) - first(cols)))
    right = floor(Int, min(window, last(cols) - max(point[2], new_point[2])))

    (-up:down, -left:right)
end

@inline get_grid(point, offsets) = @inbounds (
    (point[1] + offsets[1][begin]):(point[1] + offsets[1][end]),
    (point[2] + offsets[2][begin]):(point[2] + offsets[2][end]))


text:
module Fit

# Simulate a dataset and test functions: loglikelihood!() and fit!()
using MixedModelsBLB, LinearAlgebra, Random, Distributions, Test

Random.seed!(123)
N = 100; p = 4; q = 2; reps = 5
βtrue  = [0.1; 6; -3; 1]
σ²true = 1.5
σtrue  = sqrt(σ²true)
Σtrue  = [1.0 0; 0 3.0] # Matrix(Diagonal([2.0; 1.2; 1.0]))
Ltrue  = Matrix(cholesky(Symmetric(Σtrue)).L)
obsvec = Vector{blblmmObs{Float64}}(undef, N)
X = Matrix{Float64}(undef, reps, p)
X[:, 1] = ones(reps)
Z = Matrix{Float64}(undef, reps, q)
Z[:, 1] = ones(reps)
storage_q = Vector{Float64}(undef, q)
re_storage = Vector{Float64}(undef, q)
y = Vector{Float64}(undef, reps)
fenames = vcat("Intercept", "x" .* string.([1:1:(p-1);]))
renames = ["Intercept", "z1"]

# Generate data
for i in 1:N
    randn!(y) # y = standard normal error
    # first column intercept, remaining entries iid std normal
    @views randn!(X[:, 2:p]) #Distributions.rand!(Normal(), X[:, 2:p])
    BLAS.gemv!('N', 1., X, βtrue, σtrue, y) # y = Xβ + σtrue * standard normal error
    randn!(storage_q)
    BLAS.gemv!('N', 1., Ltrue, storage_q, 0., re_storage)
    # first column intercept, remaining entries iid std normal
    @views randn!(Z[:, 2:q]) #Distributions.rand!(Normal(), Z[:, 2:q])
    BLAS.gemv!('N', 1., Z, re_storage, 1., y) # y = Xβ + Zα + error
    # y = X * βtrue .+ Z * (Ltrue * randn(q)) .+ σtrue * randn(ns[i])
    # form a blblmmObs instance
    obsvec[i] = blblmmObs(copy(y), copy(X), copy(Z))
end

# form a LmmModel instance
lmm = blblmmModel(obsvec, fenames, renames, N)
obsN = lmm.data[N]

# Data was simulated correctly
@test obsN.y[1] ≈ -13.248956453781018 atol = 1e-5

@testset "loglikelihood!()" begin
    @test loglikelihood!(obsN, βtrue, [σ²true], Ltrue, true, true) ≈ -10.753503805549846 atol = 1e-5

    init_ls!(lmm)
    @test loglikelihood!(lmm, true, true) ≈ -1090.2939385079767 atol = 1e-5
end

@testset "fit!()" begin
    fit!(lmm)
    @test loglikelihood!(lmm, true, false) ≈ -994.0524142638102 atol = 1e-5

    @test lmm.β[1] ≈ 0.20516618221633012 atol = 1e-5

    @test norm(lmm.∇β) < 1e-5
    @test norm(lmm.∇σ²) < 1e-5
    @test norm(lmm.∇L) < 1e-5
end

end


text:
using Test

import MathOptInterface
const MOI = MathOptInterface

@testset "mixed_complementarity" begin
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    )
    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock) -> MOI.Utilities.mock_optimize!(
            mock, config.optimal_status, [2.8, 0.0, 0.8, 1.2]
        )
    )
    MOI.Test.mixed_complementaritytest(mock, config)
end

@testset "math_program_complementarity_constraints" begin
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    )
    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock) -> MOI.Utilities.mock_optimize!(
            mock, config.optimal_status, [1.0, 0.0, 3.5, 0.0, 0.0, 0.0, 3.0, 6.0]
        )
    )
    MOI.Test.math_program_complementarity_constraintstest(mock, config)
end


text:
@testset "Inverter Components" begin
    converter = AvgCnvFixedDC(690.0, 2750000.0) #S_rated goes in Watts
    @test converter isa PowerSystems.DynamicComponent
    dc_source = FixedDCSource(600.0) #Not in the original data, guessed.
    @test dc_source isa PowerSystems.DynamicComponent
    filter = LCLFilter(0.08, 0.003, 0.074, 0.2, 0.01)
    @test filter isa PowerSystems.DynamicComponent
    pll = PLL(500.0, 0.084, 4.69)
    @test pll isa PowerSystems.DynamicComponent
    virtual_H = VirtualInertia(2.0, 400.0, 20.0, 2 * pi * 50.0)
    @test virtual_H isa PowerSystems.DeviceParameter
    Q_control = ReactivePowerDroop(0.2, 1000.0)
    @test Q_control isa PowerSystems.DeviceParameter
    outer_control = VirtualInertiaQdroop(virtual_H, Q_control)
    @test outer_control isa PowerSystems.DynamicComponent
    vsc = CombinedVIwithVZ(0.59, 736.0, 0.0, 0.0, 0.2, 1.27, 14.3, 0.0, 50.0, 0.2)
    @test vsc isa PowerSystems.DynamicComponent
end

@testset "Dynamic Inverter" begin
    nodes_OMIB = [
        Bus(
            1, #number
            "Bus 1", #Name
            "REF", #BusType (REF, PV, PQ)
            0, #Angle in radians
            1.06, #Voltage in pu
            (min = 0.94, max = 1.06), #Voltage limits in pu
            69,
            nothing,
            nothing,
        ), #Base voltage in kV
        Bus(2, "Bus 2", "PV", 0, 1.045, (min = 0.94, max = 1.06), 69, nothing, nothing),
    ]

    converter = AvgCnvFixedDC(
        138.0, #Rated Voltage
        100.0,
    ) #Rated MVA

    branch_OMIB = [PSY.Line(
        "Line1", #name
        true, #available
        0.0, #active power flow initial condition (from-to)
        0.0, #reactive power flow initial condition (from-to)
        Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
        0.01, #resistance in pu
        0.05, #reactance in pu
        (from = 0.0, to = 0.0), #susceptance in pu
        18.046, #rate in MW
        1.04,
    )]  #angle limits (-min and max)

    dc_source = FixedDCSource(1500.0) #Not in the original data, guessed.

    filt = LCLFilter(
        0.08, #Series inductance lf in pu
        0.003, #Series resitance rf in pu
        0.074, #Shunt capacitance cf in pu
        0.2, #Series reactance rg to grid connection (#Step up transformer or similar)
        0.01,
    ) #Series resistance lg to grid connection (#Step up transformer or similar)

    pll = PLL(
        500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
        0.084, #k_p: PLL proportional gain
        4.69,
    ) #k_i: PLL integral gain

    virtual_H = VirtualInertia(
        2.0, #Ta:: VSM inertia constant
        400.0, #kd:: VSM damping coefficient
        20.0, #kω:: Frequency droop gain in pu
        2 * pi * 50.0,
    ) #ωb:: Rated angular frequency

    Q_control = ReactivePowerDroop(
        0.2, #kq:: Reactive power droop gain in pu
        1000.0,
    ) #ωf:: Reactive power cut-off low pass filter frequency

    outer_control = VirtualInertiaQdroop(virtual_H, Q_control)

    vsc = CombinedVIwithVZ(
        0.59, #kpv:: Voltage controller proportional gain
        736.0, #kiv:: Voltage controller integral gain
        0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
        0.0, #rv:: Virtual resistance in pu
        0.2, #lv: Virtual inductance in pu
        1.27, #kpc:: Current controller proportional gain
        14.3, #kiv:: Current controller integral gain
        0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
        50.0, #ωad:: Active damping low pass filter cut-off frequency
        0.2,
    ) #kad:: Active damping gain

    test_inverter = DynamicInverter(
        2, #number
        "DARCO", #name
        nodes_OMIB[1], #bus location
        1.0, #ω_ref
        1.02, #V_ref
        0.5, #P_ref
        0.0, #Q_ref
        100.0, #MVABase
        converter, #Converter
        outer_control, #OuterControl
        vsc, #Voltage Source Controller
        dc_source, #DC Source
        pll, #Frequency Estimator
        filt,
    ) #Output Filter

    @test test_inverter isa PowerSystems.Component

    sys = PSY.System(100)
    for bus in nodes_OMIB
        PSY.add_component!(sys, bus)
    end
    for lines in branch_OMIB
        PSY.add_component!(sys, lines)
    end
    PSY.add_component!(sys, test_inverter)

    @test collect(PSY.get_components(DynamicInverter, sys))[1] == test_inverter

end


text:
using COSMO
using Convex, SparseArrays, LinearAlgebra

export filterConvexPrices, isArbitrageFree

function isArbitrageFree(strikes::Vector{T}, callPrices::Vector{T}, forward::T)::Tuple{Bool,Int} where {T}
    for (xi, yi) in zip(strikes, callPrices)
        if yi < max(forward - xi, 0)
            return (false, i)
        end
    end

    for i = 2:length(callPrices)-1
        s0 = (callPrices[i] - callPrices[i-1]) / (strikes[i] - strikes[i-1])
        s1 = (callPrices[i] - callPrices[i+1]) / (strikes[i] - strikes[i+1])
        if s0 <= -1
            return (false, i)
        end
        if s0 >= s1
            return (false, i)
        end
        if s1 >= 0
            return (false, i + 1)
        end
    end

    return (true, -1)
end

#create closest set of arbitrage free prices
function filterConvexPrices(
    strikes::Vector{T},
    callPrices::Vector{T}, #undiscounted!
    weights::Vector{T},
    forward::T;
    tol = 1e-8,
)::Tuple{Vector{T},Vector{T},Vector{T}} where {T}
    if isArbitrageFree(strikes, callPrices, forward)[1]
        return strikes, callPrices, weights
    end
    n = length(callPrices)
    z = Variable(n)
    G = spzeros(T, 2 * n, n)
    h = zeros(T, 2 * n)
    for i = 2:n-1
        dym = (strikes[i] - strikes[i-1])
        dy = (strikes[i+1] - strikes[i])
        G[i, i-1] = -1 / dym
        G[i, i] = 1 / dym + 1 / dy
        G[i, i+1] = -1 / dy
    end
    G[1, 1] = 1 / (strikes[2] - strikes[1])
    G[1, 2] = -G[1, 1]
    G[n, n] = 1 / (strikes[n] - strikes[n-1])
    G[n, n-1] = -G[n, n]
    for i = 1:n
        h[i] = -tol
        G[n+i, i] = -1
        h[n+i] = -max(forward - strikes[i], 0) - tol
    end
    h[1] = 1 - tol
    W = spdiagm(weights)
    problem = minimize(square(norm(W * (z - callPrices))), G * z <= h)
    #solve!(problem, () -> SCS.Optimizer(verbose = 0))
    Convex.solve!(problem, () -> COSMO.Optimizer(verbose = false, eps_rel = 1e-8, eps_abs = 1e-8))
    #println("problem status is ", problem.status, " optimal value is ", problem.optval)
    strikesf = strikes
    pricesf = Convex.evaluate(z)
    return strikesf, pricesf, weights
end

#transform to X coordinate for collocation
function makeXFromUndiscountedPrices(strikesf::Vector{T}, pricesf::Vector{T}; slopeTolerance = 1e-8) where {T}
    return makeXFromUndiscountedPrices(strikesf, pricesf, s -> -norminv(-s), slopeTolerance)
end

function makeXFromUndiscountedPrices(strikesf::Vector{T}, pricesf::Vector{T}, invphi, slopeTolerance = 1e-8) where {T}
    n = length(strikesf)

    pif = ones(T, n)
    for i = 2:n-1
        dxi = strikesf[i+1] - strikesf[i]
        dxim = strikesf[i] - strikesf[i-1]
        dzi = (pricesf[i+1] - pricesf[i]) / dxi
        dzim = (pricesf[i] - pricesf[i-1]) / dxim
        s = (dxim * dzi + dxi * dzim) / (dxim + dxi)
        previous = pif[i-1] < slopeTolerance ? pif[i-1] : pif[i-1] - slopeTolerance
        pif[i] = min(-s , previous)
    end
    dzi = (pricesf[n] - pricesf[n-1]) / (strikesf[n] - strikesf[n-1])
    pif[n] = -dzi
    dzim = -pif[n-2]
    if dzi * dzim < 0 || dzi < dzim || abs(dzi) < slopeTolerance
        pif = pif[1:n-1]
        strikesf = strikesf[1:n-1]
    end
    dzim = (pricesf[2] - pricesf[1]) / (strikesf[2] - strikesf[1])
    pif[1] = -dzim
    if dzim <= -1 + slopeTolerance
        pif = pif[2:end]
        strikesf = strikesf[2:e, d]
    end
    xif = invphi.(-pif)
    return strikesf, pif, xif
end


text:
using MacroTools: flatten, isexpr, prewalk, postwalk

"""
    @staticexpand @a_macro begin
        # some code...
        @static if condition
            # conditional code
        end
        # more code...
    end

Pre-processes a macro expression, evaluating any `@static if condition` block and retaining
or dropping the inner expression based on the result. The inner expression must be valid
code given the surrounding context, as if the `@static` conditional statement did not
exist.
"""
macro staticexpand(expr::Expr)
    sentinel = gensym("sentinel")
    isexpr(expr, :macrocall) || error("Expected a macro to expand")
    # Evaluate the static statements
    expr′ = prewalk(expr) do ex
        # Search for @static macro
        isexpr(ex, :macrocall) || return ex
        ex = ex::Expr
        ex.args[1] === Symbol("@static") || return ex
        # Found it --- just let @static run and do its thing
        ret = Base.macroexpand(__module__, ex)
        # Replace a `nothing` result with a private sentinel so that we can accurately
        # filter it back out.
        return ret === nothing ? sentinel : ret
    end
    # Remove any instances of our sentinel result
    expr′ = postwalk(expr′) do ex
        isexpr(ex, :block) || return ex
        ex = ex::Expr
        filter!(x -> x !== sentinel, ex.args)
        return ex
    end
    # Then merge the true branches' blocks with the surrounding conditionals
    expr′ = postwalk(flatten, expr′)
    return esc(expr′)
end


text:
module Pages

using HTTP, JSON3, Documenter

export Endpoint, HTTP, JSON3, Documenter
export GET, HEAD, POST, PUT, DELETE, CONNECT, OPTIONS, TRACE, PATCH

include("Endpoints.jl"); using .Endpoints
include("Callbacks.jl"); using .Callbacks
include("Server.jl"); using .Server
include("Messaging.jl"); using .Messaging
const broadcast = Messaging.broadcast
include("DOM.jl"); using .DOM
include("Docs.jl"); using .Docs
include("Examples/Examples.jl"); using .Examples

end


text:
struct DFNGammaStrategy <: SVDD.InitializationStrategyGamma
    min_γ
    max_γ
end

DFNGammaStrategy() = DFNGammaStrategy(0.0, 100.0)

function dfn(nn, gamma)
    near = mean(exp(-gamma*nn_dist[2]) for nn_dist in nn[2])
    far = mean(exp(-gamma*nn_dist[end]) for nn_dist in nn[2])
    return 2*near - 2*far
end

function _dfn(data)
    kdtree = KDTree(data)
    nn = knn(kdtree, data, size(data, 2), true)
    return ((gamma) -> dfn(nn, gamma))
end

function SVDD.calculate_gamma(model, strategy::DFNGammaStrategy)
    f(gamma) = - _dfn(model.data)(gamma)
    res = Optim.optimize(f, strategy.min_γ, strategy.max_γ, abs_tol=0.1)
    return res.minimizer
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
"""
    smooth_vars_and_write_to_netcdf!(output_file,input_file,vars_to_smooth,window_h,window_t)
Take a netcdf file, a list of variables and two integers and performs a moving mean smoothing of these variables. Write the smoothed fields into a netcdf file called output_file
"""
function smooth_vars_and_write_to_netcdf!(output_file,input_file,vars_to_smooth,window_h,window_t)

    ds_orig = Dataset(input_file)
    ds_new = Dataset(output_file,"c", attrib = OrderedDict(
        "history"                   => "Filtered data from file: $input_file, using $window_h points in the space dimensions and $window_t points in the time dimension"
    ))
    # Dimensions
    for (dim,dim_size) in ds_orig.dim
        @info "Writing $dim coordinate in netcdf file"
        flush(stdout)
        ds_new.dim[dim] = dim_size # unlimited dimension
        v = ds_orig[dim]
        create_var = defVar(ds_new,dim, eltype(v.var), dimnames(v),attrib = v.var.attrib)
        create_var[:] = v[:]
    end
    # Declare variables
    for current_var in vars_to_smooth
        @info "Processing $current_var"
        flush(stdout)
        v = ds_orig[current_var]
        create_var = defVar(ds_new,current_var, eltype(v.var), dimnames(v),attrib = v.var.attrib)
        create_var[:] = filter_array(v[:],window_h,window_t)
        @info "Finished processing $current_var"
        flush(stdout)
    end
    close(ds_new)
    close(ds_orig)
end

const dimensions_of_APE_variables = Dict(
    "rho0"    => ("z","t"),
    "N2"      => ("z","t"),
    "B"       => ("x","y","z","t"),
    "QCONVEC" => ("x","y","z","t"),
    "QRAD"    => ("x","y","z","t"),
    "APE"     => ("t",),
    "KE"      => ("t",),
    "APERAD"  => ("t",),
    "APEDIA"  => ("t",),
    "APEWN2"  => ("t",),
    "APEUb2"  => ("t",),
    "APEVb2"  => ("t",),
    "APErate" => ("t",),
    "APESURF" => ("t",)
)

"""
Create NetCDF for storing ape budget variables
"""
function create_APE_netcdf(filename,var_size)

    Dataset(filename,"c") do ds

        # Dimensions

        ds.dim["x"] = var_size[1]
        ds.dim["y"] = var_size[2]
        ds.dim["z"] = var_size[3]
        ds.dim["time"] = Inf # unlimited dimension

        # Declare variables

        ncx = defVar(ds,"x", Float32, ("x",), deflatelevel=1,  attrib = OrderedDict(
            "units"                     => "m",
        ))

        ncy = defVar(ds,"y", Float32, ("y",), deflatelevel=1, attrib = OrderedDict(
            "units"                     => "m",
        ))

        ncz = defVar(ds,"z", Float32, ("z",), deflatelevel=1,  attrib = OrderedDict(
            "units"                     => "m",
            "long_name"                 => "height",
        ))

        nctime = defVar(ds,"time", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "units"                     => "d",
            "long_name"                 => "time",
        ))

        ncAPE = defVar(ds,"APE", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy",
            "units"                     => "J/m^2     ",
        ))

        ncKE = defVar(ds,"KE", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Kinetic energy",
            "units"                     => "J/m^2     ",
        ))
        ncAPERAD = defVar(ds,"APERAD", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy production by radiative processes",
            "units"                     => "W/m^2     ",
        ))
        ncAPEDIA = defVar(ds,"APEDIA", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy production by radiative processes other than radiative and surface fluxes",
            "units"                     => "W/m^2     ",
        ))
        ncAPEWN2 = defVar(ds,"APEWN2", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy conversion into kinetic energy",
            "units"                     => "W/m^2     ",
        ))
        ncAPEUB2 = defVar(ds,"APEUb2", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy advection by the U component of the wind",
            "units"                     => "W/m^2     ",
        ))
        ncAPEVB2 = defVar(ds,"APEVb2", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy advection by the V component of the wind",
            "units"                     => "W/m^2     ",
        ))
        ncAPESURF = defVar(ds,"APESURF", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy production by surface fluxes",
            "units"                     => "W/m^2     ",
        ))
        ncAPErate = defVar(ds,"APErate", Float32, ("time",), deflatelevel=1,  attrib = OrderedDict(
            "long_name"                 => "Available Potential Energy rate of change",
            "units"                     => "W/m^2     ",
        ))

        ncrho0 = defVar(ds,"rho0", Float32, ("z","time"), deflatelevel=1, shuffle=true, attrib = OrderedDict(
            "units"                     => "kg/m^3",
            "long_name"                 => "Horizontally averaged density",
        ))

        ncN2 = defVar(ds,"N2", Float32, ("z","time"), deflatelevel=1, shuffle=true, attrib = OrderedDict(
            "units"                     => "1/s^2",
            "long_name"                 => "Brunt-Vaisala frequency",
        ))

        ncQCONVEC = defVar(ds,"QCONVEC", Float32, ("x", "y", "z", "time"), deflatelevel=1, shuffle=true, attrib = OrderedDict(
            "long_name"                 => "Convective heating departure from horizontal mean",
            "units"                     => "m/s^3       ",
        ))

        ncQRAD = defVar(ds,"QRAD", Float32, ("x", "y", "z", "time"), deflatelevel=1, shuffle=true, attrib = OrderedDict(
            "long_name"                 => "Radiative heating rate departure from horizontal mean",
            "units"                     => "m/s^3     ",
        ))

        ncB = defVar(ds,"B", Float32, ("x", "y", "z", "time"), deflatelevel=1, shuffle=true, attrib = OrderedDict(
            "long_name"                 => "Buoyancy rate departure from horizontal mean",
            "units"                     => "m/s^2     ",
        ))

        for itime = 1:var_size[end]
            nctime[itime] = itime * 3600.0
        end

    end

end

"""
Create NetCDF for storing ape budget variables
"""
function set_netcdf_var!(filename,var,data)
    @info "Started writing: $var into file: $filename"
    n_dimensions = ndims(data)
    last_dimension_size = last(size(data))
    Dataset(filename,"a") do ds
        #for ind in 1:last_dimension_size
        #    #@show ndims(data),last_dimension_size, ind,  selectdim(data,n_dimensions,ind)
        #    @info "Writing step: " ind
        #    selectdim(ds[var],n_dimensions,ind)[:] =  selectdim(data,n_dimensions,ind)
        #end
        if n_dimensions == 1
            ds[var][:] = data[:]
        elseif n_dimensions == 2
            ds[var][:,:] = data[:,:]
        elseif n_dimensions == 3
            ds[var][:,:,:] = data[:,:,:]
        elseif n_dimensions == 4
            ds[var][:,:,:,:] = data[:,:,:,:]
        end
        @info "Finished writing: $var into file: $filename"
    end
end


text:
"""
Overload of LightGraphs.bfs_parents that provides an alternative outneighbors function.

See [`LightGraphs.bfs_parents`](@ref)
"""
LightGraphs.bfs_parents(g::AbstractGraph, s::Integer; dir = :out) =
    (dir == :out) ?
    LightGraphs._bfs_parents(g, s, outneighbors_ranked) :
    LightGraphs._bfs_parents(g, s, inneighbors)

"""
For the given graph and vertex, return the list of the vertex's out neighbors
ranked by edge weight between vertex and the neighbor (ascending order by default).
This function is invoked by [`bfs_parents`](@ref) and is not intended to be called directly.
    @ TODO #17 Benchmark this to see if using views would be faster.
"""
function outneighbors_ranked(g, v; order=:asc)
    alln = collect(outneighbors(g,v))
    length(alln) == 1 && return alln
    W = LightGraphs.weights(g)
    T = eltype(W)
    rw = T[]
    for n in alln
        w = W[v,n]
        push!(rw, w)
    end
    idx = Array{Int}(undef, length(rw))
    sortperm!(idx, rw; rev = (order == :desc))
    return alln[idx]
end


text:
###################################
## Bayesian Poisson Mixture Model

import StatsFuns.logsumexp
include("sampler.jl")

####################
## Types
type Gam
    # Parameters of Gamma distribution
    a::Vector{Float64}
    b::Float64
end

type BPMM
    # Parameters of Bayesian Poisson Mixture Model
    D::Int
    K::Int
    alpha::Vector{Float64}
    cmp::Vector{Gam}
end

type Poi
    # Parameters of Poisson Distribution
    lambda::Vector{Float64}
end

type PMM
    # Parameters of Poisson Mixture Model
    D::Int
    K::Int
    phi::Vector{Float64}
    cmp::Vector{Poi}
end

####################
## Common functions
function sample_PMM(bpmm::BPMM)
    cmp = Vector{Poi}()
    for c in bpmm.cmp
        lambda = Vector{Float64}()
        for d in 1 : bpmm.D
            push!(lambda, gamma_sample(c.a[d], c.b))
        end
        push!(cmp, Poi(lambda))
    end
    phi = dirichlet_sample(bpmm.alpha)
    return PMM(bpmm.D, bpmm.K, phi, cmp)
end

function sample_data(pmm::PMM, N::Int)
    X = zeros(pmm.D, N)
    Z = categorical_sample(pmm.phi, N)
    for n in 1 : N
        k = indmax(Z[:, n])
        for d in 1 : pmm.D
            X[d,n] = poisson_sample(pmm.cmp[k].lambda[d])
        end
    end
    return X, Z
end

function init_Z(X::Matrix{Float64}, bpmm::BPMM)
    N = size(X, 2)
    K = bpmm.K
    Z = categorical_sample(ones(K)/K, N)
    return Z
end

function calc_bound(X::Matrix{Float64}, pri::BPMM, pos::BPMM)
    ln_expt_Z = update_Z(pos, X)
    expt_Z = exp(ln_expt_Z)
    K, N = size(expt_Z)
    D = size(X, 1)

    expt_ln_lambda = zeros(D, K)
    expt_lambda = zeros(D, K)
    expt_lik = 0
    for k in 1 : K
        expt_ln_lambda[:,k] = digamma(pos.cmp[k].a) - log(pos.cmp[k].b)
        expt_lambda[:,k] = pos.cmp[k].a / pos.cmp[k].b
        for n in 1 : N
            expt_lik += expt_Z[k,n] * (X[:, n]' * expt_ln_lambda[:,k]
                                       - sum(expt_lambda[:,k]) - sum(lgamma(X[:,n]+1)))[1]
        end
    end

    expt_ln_pZ = sum(expt_Z' * (digamma(pos.alpha) - digamma(sum(pos.alpha))))
    expt_ln_qZ = sum(expt_Z .* ln_expt_Z)

    KL_lambda = 0
    for k in 1 : K
        KL_lambda += (sum(pos.cmp[k].a)*log(pos.cmp[k].b) - sum(pri.cmp[k].a)*log(pri.cmp[k].b)
                      - sum(lgamma(pos.cmp[k].a)) + sum(lgamma(pri.cmp[k].a))
                      + (pos.cmp[k].a - pri.cmp[k].a)' * expt_ln_lambda[:,k]
                      + (pri.cmp[k].b - pos.cmp[k].b) * sum(expt_lambda[:,k])
                      )[1]
    end
    KL_pi = (lgamma(sum(pos.alpha)) - lgamma(sum(pri.alpha))
             - sum(lgamma(pos.alpha)) + sum(lgamma(pri.alpha))
             + (pos.alpha - pri.alpha)' * (digamma(pos.alpha) - digamma(sum(pos.alpha)))
             )[1]

    VB = expt_lik + expt_ln_pZ - expt_ln_qZ - (KL_lambda + KL_pi)
    return VB
end

function add_stats(bpmm::BPMM, X::Matrix{Float64}, Z::Matrix{Float64})
    D = bpmm.D
    K = bpmm.K
    sum_Z = sum(Z, 2)
    alpha = [bpmm.alpha[k] + sum_Z[k] for k in 1 : K]
    cmp = Vector{Gam}()

    XZ = X*Z';
    for k in 1 : K
        a = [(bpmm.cmp[k].a[d] + XZ[d,k])::Float64 for d in 1 : D]
        b = bpmm.cmp[k].b + sum_Z[k]
        push!(cmp, Gam(a, b))
    end
    return BPMM(D, K, alpha, cmp)
end

remove_stats(bpmm::BPMM, X::Matrix{Float64}, Z::Matrix{Float64}) = add_stats(bpmm, X, -Z)

####################
## used for Variational Inference
function update_Z(bpmm::BPMM, X::Matrix{Float64})
    D, N = size(X)
    K = bpmm.K
    ln_expt_Z = zeros(K, N)
    tmp = zeros(K)

    sum_digamma_tmp = digamma(sum(bpmm.alpha))
    for k in 1 : K
        tmp[k] = - sum(bpmm.cmp[k].a) / bpmm.cmp[k].b
        tmp[k] += digamma(bpmm.alpha[k]) - sum_digamma_tmp
    end
    ln_lambda_X = [X'*(digamma(bpmm.cmp[k].a) - log(bpmm.cmp[k].b)) for k in 1 : K]
    for n in 1 : N
        tmp_ln_pi =  [tmp[k] + ln_lambda_X[k][n] for k in 1 : K]
        ln_expt_Z[:,n] = tmp_ln_pi - logsumexp(tmp_ln_pi)
    end
    return ln_expt_Z
end

function winner_takes_all(Z::Matrix{Float64})
    Z_ret = zeros(size(Z))
    for n in 1 : size(Z_ret, 2)
        idx = indmax(Z[:,n])
        Z_ret[idx,n] = 1
    end
    return Z_ret
end

####################
## used for Gibbs Sampling
function sample_Z_GS(pmm::PMM, X::Matrix{Float64})
    D, N = size(X)
    K = pmm.K
    Z = zeros(K, N)

    tmp = [-sum(pmm.cmp[k].lambda) + log(pmm.phi[k]) for k in 1 : K]
    ln_lambda_X = [X'*log(pmm.cmp[k].lambda) for k in 1 : K]
    for n in 1 : N
        tmp_ln_phi = [(tmp[k] + ln_lambda_X[k][n])::Float64 for k in 1 : K]
        tmp_ln_phi = tmp_ln_phi - logsumexp(tmp_ln_phi)
        Z[:,n] = categorical_sample(exp(tmp_ln_phi))
    end
    return Z
end

####################
## used for Collapsed Gibbs Sampling
function calc_ln_NB(Xn::Vector{Float64}, gam::Gam)
    ln_lkh = [(gam.a[d]*log(gam.b)
               - lgamma(gam.a[d])
               + lgamma(Xn[d] + gam.a[d])
               - (Xn[d] + gam.a[d])*log(gam.b + 1)
               )::Float64 for d in 1 : size(Xn, 1)]
    return sum(ln_lkh)
end

function sample_Zn(Xn::Vector{Float64}, bpmm::BPMM)
    ln_tmp = [(calc_ln_NB(Xn, bpmm.cmp[k]) + log(bpmm.alpha[k])) for k in 1 : bpmm.K]
    ln_tmp = ln_tmp -  logsumexp(ln_tmp)
    Zn = categorical_sample(exp(ln_tmp))
    return Zn
end

function sample_Z_CGS(Z::Matrix{Float64}, X::Matrix{Float64}, bpmm::BPMM)
    D, N = size(X)
    K = size(Z, 1)
    for n in randperm(N)
        # remove
        bpmm = remove_stats(bpmm, X[:,[n]], Z[:,[n]])
        # sample
        Z[:,n] = sample_Zn(X[:,n], bpmm)
        # insert
        bpmm = add_stats(bpmm, X[:,[n]], Z[:,[n]])
    end
    return Z, bpmm
end

####################
## Algorithm main
function learn_BPMM_VI(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    ## Variational Inference
    ## for Bayesian Poisson Mixture Model

    # initialisation
    expt_Z = init_Z(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, expt_Z)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # E-step
        expt_Z = exp(update_Z(bpmm, X))
        # M-step
        bpmm = add_stats(prior_bpmm, X, expt_Z)
        # calc VB
        VB[i] = calc_bound(X, prior_bpmm, bpmm)
    end

    # assign binary values
    #Z = winner_takes_all(expt_Z)
    Z = expt_Z
    return Z, bpmm, VB
end

function learn_BPMM_GS(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    ## Gibbs Sampling
    ## for Bayesian Poisson Mixture Model

    # initialisation
    Z = init_Z(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, Z)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # sample parameters
        pmm = sample_PMM(bpmm)
        # sample latent variables
        Z = sample_Z_GS(pmm, X)
        # update current model
        bpmm = add_stats(prior_bpmm, X, Z)
        # calc VB
        VB[i] = calc_bound(X, prior_bpmm, bpmm)
    end

    return Z, bpmm, VB
end

function learn_BPMM_CGS(X::Matrix{Float64}, prior_bpmm::BPMM, max_iter::Int)
    ## Collapsed Gibbs Sampling
    ## for Bayesian Poisson Mixture Model

    # initialisation
    Z = init_Z(X, prior_bpmm)
    bpmm = add_stats(prior_bpmm, X, Z)
    VB = NaN * zeros(max_iter)

    # inference
    for i in 1 : max_iter
        # directly sample Z
        Z, bpmm = sample_Z_CGS(Z, X, bpmm)
        # calc VB
        VB[i] = calc_bound(X, prior_bpmm, bpmm)
    end

    return Z, bpmm, VB
end

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
module Job

using ..Db: ENV_CONN, mangrove_reviews, DbInput, insert
using ..Model: get_chains, subs, qualities, kids, neutralities
using Turing: mean, group

function store_inferred(conn_string::String = ENV_CONN)
    chains = conn_string |> mangrove_reviews |> get_chains
    sub_qualities = mean(group(chains, :qualities))
    kid_neutralities = mean(group(chains, :biases))
    input = DbInput(
        subs(sub_qualities),
        qualities(sub_qualities),
        kids(kid_neutralities),
        neutralities(kid_neutralities)
    )
    insert(conn_string, input)
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
# ------------------- CLASSES ------------------------
#include("defArc.jl")

type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    flow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64, freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0.)

#include("fitTraffic.jl")

## Solve an inverse tarffic problem over polynomials
## of degree at most d
## optionally use a regularizer from the poly kernel

using JuMP
using Gurobi
#using LinAlg

polyEval(coeffs, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

polyEval(coeffs::Array{Float64, 1}, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

bpacost(flow::Float64, capacity::Float64, freeflowtime::Float64) = freeflowtime*(1 + .15 * (flow/capacity)^4)
bpacost(flow::Float64, arc) = bpacost(flow, arc.capacity, arc.freeflowtime)
bpacost(arc::Arc) = bpacost(arc.flow, arc)

# --------------------  FUNCTIONS  -------------------

function setUpFitting(deg::Int, c::Float64)
    m = 0
	m = Model(solver=GurobiSolver(OutputFlag=false))

	@variable(m, coeffs[1:deg+1])
	@variable(m, Calphas[1:deg+1])

	#build the graham matrix; cf. Ref. [21] (Regularization Networks and Support Vector Machines), page 47
	samples = linspace(0, 1, deg + 1)
	k(x,y) = (c + x*y)^deg
	K = [ k(x,y) for x = samples, y=samples]
	K = convert(Array{Float64, 2}, K)
	#assert(rank(K) == deg+1)

	C = chol(K + 1e-6* eye(deg+1))
	for i=1:deg + 1
		@constraint(m, polyEval(coeffs, samples[i]) == sum{C[j, i] * Calphas[j], j=1:deg+1})
	end

	@variable(m, reg_term >= 0)
	reg_term_ = QuadExpr(Calphas[:], Calphas[:], ones(deg+1), AffExpr())

	@constraint(m, reg_term >= reg_term_)

	return m, coeffs, reg_term
end

function fixCoeffs(m, fcoeffs, coeffs)
	for (fc, c) in zip(fcoeffs, coeffs[:])
		@constraint(m, fc == c)
	end
end

function addResid(m, coeffs, ys, demands, arcs, scaling)
	@variable(m, resid)
	@variable(m, dual_cost)
	@variable(m, primal_cost)

	@constraint(m, dual_cost == sum{demands[(s,t)] * (ys[(s,t), t] - ys[(s,t), s]), (s,t)=keys(demands)})
	@constraint(m, primal_cost == sum{a.flow * a.freeflowtime * polyEval(coeffs, a.flow/a.capacity), a=values(arcs)})

	@constraint(m, resid >= (dual_cost - primal_cost) / scaling )
	@constraint(m, resid >= (primal_cost - dual_cost) / scaling )
	return resid
end

function addIncreasingCnsts(m, coeffs, arcs; TOL=0.)
	sorted_flows = sort([a.flow / a.capacity for a in values(arcs)])
	@constraint(m, polyEval(coeffs, 0) <= polyEval(coeffs, sorted_flows[1]))
	for i = 2:length(sorted_flows)
		@constraint(m, polyEval(coeffs, sorted_flows[i-1]) <= polyEval(coeffs, sorted_flows[i]) + TOL)
	end
    @constraint(m, coeffs[1] == 1)  # enforce g(0) = 1
end

#equates the total cost of the network to the true total cost
function normalize(m, coeffs, tot_true_cost::Float64, arcs)
	@constraint(m,
		sum{a.freeflowtime * a.flow * polyEval(coeffs, a.flow / a.capacity), a=values(arcs)} == tot_true_cost)
end

function normalize(m, coeffs, scaled_flow::Float64, cost::Float64)
	@constraint(m, polyEval(coeffs, scaled_flow) == cost)
end

function normalize(m, coeffs, scaled_flows::Array{Float64, 1}, avgCost::Float64)
    @constraint(m, sum{polyEval(coeffs, f), f=scaled_flows} == avgCost * length(scaled_flows))
end

function addNetworkCnsts(m, coeffs, demands, arcs, numNodes)
	@variable(m, ys[keys(demands), 1:numNodes])
	for k = keys(arcs)
		a = arcs[k]
		rhs = a.freeflowtime * polyEval(coeffs, a.flow/a.capacity)
		for od in keys(demands)
			@constraint(m, ys[od, k[2]] - ys[od, k[1]] <= rhs)
		end
	end
	return ys
end

function read_demand_file(file_path)
    file = open(file_path)
    demands = Dict()
    n = n_nodes  # number of nodes
    for i = 1:n
        demands[(i,i)] = 0.0
    end
    for line in eachline(file)
        OD_demand = split(line, ",")
        key, value = (parse(Int, OD_demand[1]), parse(Int, OD_demand[2])), parse(Float64, split(OD_demand[3], "\n")[1])
        demands[key] = value
    end
    close(file)
    return demands
end

using JSON
function read_flow_after_cons(file_path)

    flow_after_conservation = readstring(file_path);
    flow_after_conservation = replace(flow_after_conservation, "NaN", 0);
    flow_after_conservation = JSON.parse(flow_after_conservation);
    return flow_after_conservation
end

function read_link_day_min_dict(file_path)
    link_day_minute_Apr_dict = readstring(file_path);
    link_day_minute_Apr_dict = replace(link_day_minute_Apr_dict, "NaN", 0);
    link_day_minute_Apr_dict = JSON.parse(link_day_minute_Apr_dict);
    return link_day_minute_Apr_dict
end

##########
#Fitting Funcs
##########

function train(indices, lam::Float64, deg::Int, c::Float64, demand_data, flow_data, arcs, vArcs, demands; fcoeffs=nothing)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))
    m, coeffs, reg_term = setUpFitting(deg, c)

    addIncreasingCnsts(m, coeffs, arcs, TOL=1e-8)  #uses the original obs flows

    avgCost = mean( [bpacost(a.flow, a.capacity, 1.0) for a in values(arcs)] )
    normalize(m, coeffs, [a.flow / a.capacity for a in values(arcs)], avgCost)

    resids = Variable[]

    for i = indices
        #copy the flow data over to the arcs, demand data to demands (slow)
        for (ix, a) in enumerate(vArcs)
            a.flow = flow_data[ix, i]
        end
        for odpair in keys(demands)
            demands[odpair] = demand_data[odpair][i]
        end

        #Dual Feasibility
        ys = addNetworkCnsts(m, coeffs, demands, arcs, numNodes)

        #add the residual for this data point
        push!(resids, addResid(m, coeffs, ys, demands, arcs, 1e6))
    end

    if fcoeffs != nothing
        fixCoeffs(m, fcoeffs, coeffs)
    end
    @objective(m, Min, sum{resids[i], i = 1:length(resids)} + lam*reg_term)
    solve(m)
    #println(getObjectiveValue(m) - lam * getValue(reg_term) )
    return [getvalue(coeffs[i]) for i =1:length(coeffs)]
end

# --------------------  PARAMETERS -------------------

#Importing parameters
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport parameters_julia

out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
year = parameters_julia.year
insts = parameters_julia.instances_ID
deg_grid = parameters_julia.deg_grid
c_grid = parameters_julia.c_grid
lamb_grid = parameters_julia.lamb_grid
#println(deg_grid)
week_day_Apr_list_1 = parameters_julia.week_day_list_1
week_day_Apr_list_2 = parameters_julia.week_day_list_2
week_day_Apr_list_3 = parameters_julia.week_day_list_3

n_nodes = parameters_julia.n_zones
month = parameters_julia.month
month_w = parameters_julia.month_w
year = parameters_julia.year

using Graphs
#for inst in insts
    inst = "NT"
	# Assing file names
	link_day_min_path = out_dir * "link_min_dict" * files_ID * ".json"
    #flow_path =  out_dir * "flows_after_QP" * files_ID * "_" * inst * ".json"
    flow_path =  out_dir * "density_links" * files_ID * "_" * inst * ".json"
    save_file_path = out_dir * "coeffs_dict_" * month_w * "_" * inst * ".json"
    OD_demand_path = out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * "full" * "_weekday_" * inst * files_ID * ".txt"

    demands = read_demand_file(OD_demand_path)
    flow_after_conservation = read_flow_after_cons(flow_path)
    link_day_minute_Apr_dict = read_link_day_min_dict(link_day_min_path)
    key1 = [uppercase(key1) for key1 in keys(flow_after_conservation)]
    key2 = key1[1]
    arc_ids = [uppercase(ids) for ids in keys(flow_after_conservation[key2])]
    arcs_1 = Dict[]
    for j in week_day_Apr_list_1
        arcs_1_ = Dict()
        for i in arc_ids
            key = "link_$(i)_$(year)_$(month)_$(j)"
            initNode = link_day_minute_Apr_dict[key]["init_node"]
            termNode = link_day_minute_Apr_dict[key]["term_node"]
            capacity = link_day_minute_Apr_dict[key]["capac_$inst"]
            freeflowtime = link_day_minute_Apr_dict[key]["free_flow_time"]
            flow = link_day_minute_Apr_dict[key]["avg_flow_$inst"]
            arcs_1_[(initNode, termNode)] = Arc(initNode, termNode, capacity, freeflowtime, flow)
        end
        push!(arcs_1, arcs_1_)
    end

    arcs_2 = Dict[]
    for j in week_day_Apr_list_2
        arcs_2_ = Dict()
        for i in arc_ids
            key = "link_$(i)_$(year)_$(month)_$(j)"
            initNode = link_day_minute_Apr_dict[key]["init_node"]
            termNode = link_day_minute_Apr_dict[key]["term_node"]
            capacity = link_day_minute_Apr_dict[key]["capac_$inst"]
            freeflowtime = link_day_minute_Apr_dict[key]["free_flow_time"]
            flow = link_day_minute_Apr_dict[key]["avg_flow_$inst"]
            arcs_2_[(initNode, termNode)] = Arc(initNode, termNode, capacity, freeflowtime, flow)
        end
        push!(arcs_2, arcs_2_)
    end
    arcs_3 = Dict[]
    for j in week_day_Apr_list_3
        arcs_3_ = Dict()
        for i in arc_ids
            key = "link_$(i)_$(year)_$(month)_$(j)"
            initNode = link_day_minute_Apr_dict[key]["init_node"]
            termNode = link_day_minute_Apr_dict[key]["term_node"]
            capacity = link_day_minute_Apr_dict[key]["capac_$inst"]
            freeflowtime = link_day_minute_Apr_dict[key]["free_flow_time"]
            flow = link_day_minute_Apr_dict[key]["avg_flow_$inst"]
            arcs_3_[(initNode, termNode)] = Arc(initNode, termNode, capacity, freeflowtime, flow)
        end
        push!(arcs_3, arcs_3_)
    end
    ##########
    # Set up demand data and flow data
    ##########
    numData = length(arcs_1)
    sigma = .0
    flow_data_1 = Array(Float64, length(arcs_1[1]), numData)
    flow_data_2 = Array(Float64, length(arcs_2[1]), numData)
    flow_data_3 = Array(Float64, length(arcs_3[1]), numData)

    demand_data = Dict()

    numNodes = maximum(map(pair->pair[1], keys(demands)))
    g = simple_inclist(numNodes, is_directed=true)
    vArcs = Arc[]
    for arc in values(arcs_1[1])
        add_edge!(g, arc.initNode, arc.termNode)
        push!(vArcs, arc)
    end

    for iRun = 1:numData
        for odpair in keys(demands)
            if ! haskey(demand_data, odpair)

                demand_data[odpair] = [demands[odpair], ]
            else
                push!(demand_data[odpair], demands[odpair])
            end
        end

        #flow_data[:, iRun] = [a.flow::Float64 for a in vArcs]
        flow_data_1[:, iRun] = [a.flow::Float64 for a in values(arcs_1[iRun])]
        flow_data_2[:, iRun] = [a.flow::Float64 for a in values(arcs_2[iRun])]
        flow_data_3[:, iRun] = [a.flow::Float64 for a in values(arcs_3[iRun])]
    end

    # Train

    coeffs_dict_Apr_PM = Dict()
   # println(deg_grid)
    coeffs_keys = Vector()
   #coeffs_keys_2 = []
   #coeffs_keys_3 = []
    for deg in deg_grid
        for c in c_grid
            for lam in lamb_grid
                #println(numData)

                coeffs_dict_Apr_1 = train(1:numData,lam,deg,c,demand_data,flow_data_1,arcs_1[1], vArcs, demands)
                #println(coeffs_dict_Apr_1[1])
                #println(isa(coeffs_dict_Apr_1[1], Number))

                if isnan(coeffs_dict_Apr_1[1]) == false
                    coeffs_dict_Apr_PM[(deg,c,lam,1)] = coeffs_dict_Apr_1
                    append!(coeffs_keys, ["(" * string(deg) * ", " *  string(c) * ", " * string(lam) * ", 1)"])
                end

                coeffs_dict_Apr_2 = train(1:numData,lam,deg,c,demand_data,flow_data_1,arcs_2[1], vArcs, demands)
                if isnan(coeffs_dict_Apr_2[1]) == false
                    coeffs_dict_Apr_PM[(deg,c,lam,2)] = coeffs_dict_Apr_2
                    append!(coeffs_keys, ["(" * string(deg) * ", " *  string(c) * ", " * string(lam) * ", 2)"])
                end

                coeffs_dict_Apr_3 = train(1:numData,lam,deg,c,demand_data,flow_data_1,arcs_3[1], vArcs, demands)
                if isnan(coeffs_dict_Apr_3[1]) == false
                    coeffs_dict_Apr_PM[(deg,c,lam,3)] = coeffs_dict_Apr_3
                    append!(coeffs_keys, ["(" * string(deg) * ", " *  string(c) * ", " * string(lam) * ", 3)"])
                end

                #coeffs_dict_Apr_2 = train(1:numData,lam,deg,c,demand_data,flow_data_2,arcs_1[1], vArcs, demands)
                #coeffs_dict_Apr_3 = train(1:numData,lam,deg,c,demand_data,flow_data_3,arcs_1[1], vArcs, demands)

                #coeffs_dict_Apr_PM[(deg,c,lam,2)] = train(1:numData,lam,deg,c,demand_data,flow_data_2,arcs_2[1], vArcs, demands)
                #coeffs_dict_Apr_PM[(deg,c,lam,3)] = train(1:numData,lam,deg,c,demand_data,flow_data_3,arcs_3[1], vArcs, demands)

            end
        end
    end

    outfile = open(save_file_path, "w")
    JSON.print(outfile, coeffs_dict_Apr_PM)
    close(outfile)

    outfile = open(out_dir * "coeffs_keys_" * month_w * "_" * inst * ".json", "w")
    JSON.print(outfile, coeffs_keys)
    close(outfile)
#end

max_stars_repo_path: 01. InvOpt_salo_f.jl

max_issues_repo_path: 01. InvOpt_salo_f.jl

max_forks_repo_path: 01. InvOpt_salo_f.jl


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

mutable struct IoK8sApiCoreV1PodDNSConfig <: SwaggerModel
    nameservers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nameservers
    options::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodDNSConfigOption} } # spec name: options
    searches::Any # spec type: Union{ Nothing, Vector{String} } # spec name: searches

    function IoK8sApiCoreV1PodDNSConfig(;nameservers=nothing, options=nothing, searches=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("nameservers"), nameservers)
        setfield!(o, Symbol("nameservers"), nameservers)
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        validate_property(IoK8sApiCoreV1PodDNSConfig, Symbol("searches"), searches)
        setfield!(o, Symbol("searches"), searches)
        o
    end
end # type IoK8sApiCoreV1PodDNSConfig

const _property_map_IoK8sApiCoreV1PodDNSConfig = Dict{Symbol,Symbol}(Symbol("nameservers")=>Symbol("nameservers"), Symbol("options")=>Symbol("options"), Symbol("searches")=>Symbol("searches"))
const _property_types_IoK8sApiCoreV1PodDNSConfig = Dict{Symbol,String}(Symbol("nameservers")=>"Vector{String}", Symbol("options")=>"Vector{IoK8sApiCoreV1PodDNSConfigOption}", Symbol("searches")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodDNSConfig }) = collect(keys(_property_map_IoK8sApiCoreV1PodDNSConfig))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodDNSConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodDNSConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodDNSConfig }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodDNSConfig[property_name]

function check_required(o::IoK8sApiCoreV1PodDNSConfig)
    true
end

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


text:
## %PREDICTMD_GENERATED_BY%

import PredictMDExtra
PredictMDExtra.import_all()

import PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
import CSVFiles
import DataFrames
import FileIO
import RDatasets
import Random
import Test
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "boston_housing", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin data preprocessing code

Random.seed!(999)

df = RDatasets.dataset("MASS", "Boston")

## PredictMD requires that you provide your data in a DataFrame.

## If your data are in a CSV file (e.g. "data.csv"), load them into
## a DataFrame named `df` with:
## df = DataFrames.DataFrame(CSVFiles.load("data.csv"; type_detect_rows = 10_000))

## If your data are in a gzipped CSV file (e.g. "data.csv.gz"), load them into
## a DataFrame named `df` with:
## df = DataFrames.DataFrame(CSVFiles.load(CSVFiles.File(CSVFiles.format"CSV", "data.csv.gz"); type_detect_rows = 10_000))

## If your data are in some other format, use the appropriate Julia package to
## load your data into a DataFrame named `df`.

# PREDICTMD IF INCLUDE TEST STATEMENTS
df1 = DataFrames.DataFrame()
df1[:x] = randn(5)
df1_filename = joinpath(PredictMD.maketempdir(), "df1.csv")
CSVFiles.save(df1_filename, df1)
df2 = DataFrames.DataFrame(CSVFiles.load(df1_filename; type_detect_rows = 10_000))
Test.@test( all(df1[:x] .== df2[:x]) )
df3 = DataFrames.DataFrame()
df3[:y] = randn(5)
df3_filename = joinpath(PredictMD.maketempdir(), "df3.csv.gz")
CSVFiles.save(CSVFiles.File(CSVFiles.format"CSV", df3_filename), df3)
df4 = DataFrames.DataFrame(CSVFiles.load(CSVFiles.File(CSVFiles.format"CSV", df3_filename); type_detect_rows = 10_000))
Test.@test( all(df3[:y] .== df4[:y]) )
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

categorical_feature_names = Symbol[]
continuous_feature_names = Symbol[
    :Crim,
    :Zn,
    :Indus,
    :Chas,
    :NOx,
    :Rm,
    :Age,
    :Dis,
    :Rad,
    :Tax,
    :PTRatio,
    :Black,
    :LStat,
    ]
categorical_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "categorical_feature_names.jld2",
    )
continuous_feature_names_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "continuous_feature_names.jld2",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("", categorical_feature_names_filename)
@debug("", isfile(categorical_feature_names_filename))
@debug("", dirname(categorical_feature_names_filename))
@debug("", isdir(dirname(categorical_feature_names_filename)))
@debug("", joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
@debug("", isdir(joinpath(PROJECT_OUTPUT_DIRECTORY, "data")))
@debug("", PROJECT_OUTPUT_DIRECTORY)
@debug("", isdir(PROJECT_OUTPUT_DIRECTORY))
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

FileIO.save(
    categorical_feature_names_filename,
    "categorical_feature_names",
    categorical_feature_names,
    )
FileIO.save(
    continuous_feature_names_filename,
    "continuous_feature_names",
    continuous_feature_names,
    )
feature_names = vcat(categorical_feature_names, continuous_feature_names)

single_label_name = :MedV

continuous_label_names = Symbol[single_label_name]
categorical_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

df = df[:, vcat(feature_names, label_names)]
DataFrames.dropmissing!(df; disallowmissing=true,)
PredictMD.shuffle_rows!(df)

PredictMD.fix_column_types!(
    df;
    categorical_feature_names = categorical_feature_names,
    continuous_feature_names = continuous_feature_names,
    categorical_label_names = categorical_label_names,
    continuous_label_names = continuous_label_names,
    )
PredictMD.check_column_types(
    df;
    categorical_feature_names = categorical_feature_names,
    continuous_feature_names = continuous_feature_names,
    categorical_label_names = categorical_label_names,
    continuous_label_names = continuous_label_names,
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
Test.@test PredictMD.check_no_constant_columns(df)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

features_df = df[feature_names]
labels_df = df[label_names]

DataFrames.describe(logger_stream, labels_df[single_label_name])

(trainingandtuning_features_df,
    trainingandtuning_labels_df,
    testing_features_df,
    testing_labels_df,) = PredictMD.split_data(
        features_df,
        labels_df,
        0.75,
        )
(training_features_df,
    training_labels_df,
    tuning_features_df,
    tuning_labels_df,) = PredictMD.split_data(
        trainingandtuning_features_df,
        trainingandtuning_labels_df,
        2/3,
        )

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
FileIO.save(trainingandtuning_features_df_filename, trainingandtuning_features_df)
FileIO.save(trainingandtuning_labels_df_filename, trainingandtuning_labels_df)
FileIO.save(testing_features_df_filename, testing_features_df)
FileIO.save(testing_labels_df_filename, testing_labels_df)
FileIO.save(training_features_df_filename, training_features_df)
FileIO.save(training_labels_df_filename, training_labels_df)
FileIO.save(tuning_features_df_filename, tuning_features_df)
FileIO.save(tuning_labels_df_filename, tuning_labels_df)

# PREDICTMD IF INCLUDE TEST STATEMENTS
temp_dir = mktempdir()
atexit(() -> rm(temp_dir; force = true, recursive = true))
temp_gz_filename_testing_features_df = joinpath(temp_dir, "temp_filename_testing_features_df.csv.gz")
temp_gz_filename_testing_labels_df = joinpath(temp_dir, "temp_filename_testing_labels_df.csv.gz")

original_testing_features_df = deepcopy(testing_features_df)
original_testing_labels_df = deepcopy(testing_labels_df)

CSVFiles.save(CSVFiles.File(CSVFiles.format"CSV", temp_gz_filename_testing_features_df), original_testing_features_df)
CSVFiles.save(CSVFiles.File(CSVFiles.format"CSV", temp_gz_filename_testing_labels_df), original_testing_labels_df)

roundtrip_testing_features_df = DataFrames.DataFrame(CSVFiles.load(CSVFiles.File(CSVFiles.format"CSV", temp_gz_filename_testing_features_df); type_detect_rows = 10_000))
roundtrip_testing_labels_df = DataFrames.DataFrame(CSVFiles.load(CSVFiles.File(CSVFiles.format"CSV", temp_gz_filename_testing_labels_df); type_detect_rows = 10_000))

for column in names(roundtrip_testing_features_df)
    for i = 1:size(roundtrip_testing_features_df, 1)
        if ismissing(roundtrip_testing_features_df[i, column])
            Test.@test ismissing(original_testing_features_df[i, column])
        else
            Test.@test roundtrip_testing_features_df[i, column] == original_testing_features_df[i, column]
        end
    end
end

for column in names(roundtrip_testing_labels_df)
    for i = 1:size(roundtrip_testing_labels_df, 1)
        if ismissing(roundtrip_testing_labels_df[i, column])
            Test.@test ismissing(original_testing_labels_df[i, column])
        else
            Test.@test roundtrip_testing_labels_df[i, column] == original_testing_labels_df[i, column]
        end
    end
end

rm(temp_dir; force = true, recursive = true)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End data preprocessing code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "boston_housing", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 32
modulus = 4800
base_frequency = 125000000
samples_per_period = div(modulus, dec)
periods_per_frame = 10

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)

modeDAC(rp, "STANDARD")

frequencyDAC(rp, 1, 1, base_frequency / modulus)
println(" frequency = $(frequencyDAC(rp,1,1))")

signalTypeDAC(rp, 1 , "SINE")
amplitudeDAC(rp, 1, 1, 0.5)
offsetDAC(master(rp), 1, 0)
phaseDAC(rp, 1, 1, 0.0 )

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "INTERNAL")
masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)
uFirstPeriod = readFrames(rp, 0, 1)
sleep(0.1)
fr = currentFrame(rp)
@show fr

uCurrentPeriod = readFrames(rp, fr, 1)

figure(1)
clf()
plot(vec(uFirstPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,1,:,:]))
legend(("first period", "current period"))
savefig("images/simple.png")

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBlocks.ConstGate

@testset "construct" begin
    @test_throws AssertionError put(2, 1 => swap(2, 1, 2))
end

@testset "apply!" begin
    n = 6
    Reg = rand_state(n)

    pb = PutBlock{n}(X, (3,))
    rb = repeat(n, X, (3,))
    @test apply!(copy(Reg), pb) ≈ apply!(copy(Reg), rb)
    @test pb |> applymatrix ≈ mat(pb)
    pb = PutBlock{n}(rot(X, 0.3), (3,))
    @test pb |> applymatrix ≈ mat(pb)

    pb = PutBlock{n}(rot(CNOT, 0.3), (6, 3))
    @test pb |> applymatrix ≈ mat(pb)
    pb = PutBlock{n}(matblock(mat(rot(CNOT, 0.3)) |> Matrix), (6, 3))
    @test pb |> applymatrix ≈ mat(pb)

    pb = PutBlock{n}(rot(X, 0.3), (3,))
    @test pb |> applymatrix ≈ mat(pb)
    pb = PutBlock{n}(matblock(mat(rot(X, 0.3)) |> Matrix), (3,))
    @test pb |> applymatrix ≈ mat(pb)

    Cb = control(n, (3,), 5 => X)
    pb = PutBlock{n}(CNOT, (3, 5))
    @test apply!(copy(Reg), Cb) ≈ apply!(copy(Reg), pb)

    blks = [control(2, 1, 2 => Z)]
    @test (chsubblocks(pb, blks) |> subblocks .== blks) |> all

    pb = PutBlock{1000}(X, (3,))
    @test pb |> ishermitian
    @test pb |> isunitary
    @test pb |> isreflexive
end

@testset "test swap gate" begin
    include("swap_gate.jl")
end

@testset "rotation gate" begin
    reg = rand_state(5)
    @test apply!(copy(reg), put(5, 2 => Rx(0.3))) |> state ≈ mat(put(5, 2 => Rx(0.3))) * reg.state
    @test apply!(copy(reg), put(5, 2 => Ry(0.3))) |> state ≈ mat(put(5, 2 => Ry(0.3))) * reg.state
    @test apply!(copy(reg), put(5, 2 => Rz(0.3))) |> state ≈ mat(put(5, 2 => Rz(0.3))) * reg.state
end


text:
module Graphs

# modules
export
    PathDirections

# types
export
    DirectedGraph,
    SpanningTree,
    TreePath,
    Edge,
    Vertex

# functions
export
    vertextype,
    edgetype,
    vertices,
    edges,
    source,
    target,
    out_edges,
    in_edges,
    out_neighbors,
    in_neighbors,
    num_vertices,
    num_edges,
    add_vertex!,
    add_edge!,
    remove_vertex!,
    remove_edge!,
    rewire!,
    replace_edge!,
    reindex!,
    root,
    edge_to_parent,
    edges_to_children,
    tree_index,
    ancestors,
    lowest_common_ancestor,
    subtree_vertices,
    direction,
    directions

using Base.Iterators: flatten
import SparseArrays: sparsevec

include("abstract.jl")
include("directed_graph.jl")
include("spanning_tree.jl")
include("tree_path.jl")
include("edge_vertex_wrappers.jl")

end # module