text:
###############################################################################
#
#   fq_abs_series.jl : Power series over flint fmpz integers
#
###############################################################################

export fq_abs_series, FqAbsSeriesRing, PowerSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fq_abs_series)
   if iszero(a)
      return deepcopy(a)    # 0 + O(x^n)
   end
   prec = length(a) - 1
   prec < 0 && throw(DomainError())
   z = fq_abs_series(base_ring(a), Array{fq}(0), 0, prec)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FqAbsSeriesRing}) = fq_abs_series

parent_type(::Type{fq_abs_series}) = FqAbsSeriesRing

base_ring(R::FqAbsSeriesRing) = R.base_ring

var(a::FqAbsSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FqAbsSeriesRing) = R.prec_max

function normalise(a::fq_abs_series, len::Int)
   ctx = base_ring(a)
   if len > 0
      c = base_ring(a)()
      ccall((:fq_poly_get_coeff, :libflint), Void,
         (Ref{fq}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
          c, a, len - 1, ctx)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fq_poly_get_coeff, :libflint), Void,
            (Ref{fq}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
             c, a, len - 1, ctx)
      end
   end

   return len
end

function length(x::fq_abs_series)
   return ccall((:fq_poly_length, :libflint), Int,
                (Ref{fq_abs_series}, Ref{FqFiniteField}), x, base_ring(x))
end

precision(x::fq_abs_series) = x.prec

function coeff(x::fq_abs_series, n::Int)
   if n < 0
      return base_ring(x)()
   end
   z = base_ring(x)()
   ccall((:fq_poly_get_coeff, :libflint), Void,
         (Ref{fq}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
          z, x, n, base_ring(x))
   return z
end

zero(R::FqAbsSeriesRing) = R(0)

one(R::FqAbsSeriesRing) = R(1)

function gen(R::FqAbsSeriesRing)
   S = base_ring(R)
   z = fq_abs_series(S, [S(0), S(1)], 2, max_precision(R))
   z.parent = R
   return z
end

function deepcopy_internal(a::fq_abs_series, dict::ObjectIdDict)
   z = fq_abs_series(base_ring(a), a)
   z.prec = a.prec
   z.parent = parent(a)
   return z
end

function isgen(a::fq_abs_series)
   return precision(a) == 0 || ccall((:fq_poly_is_gen, :libflint), Bool,
                   (Ref{fq_abs_series}, Ref{FqFiniteField}), a, base_ring(a))
end

iszero(a::fq_abs_series) = length(a) == 0

isunit(a::fq_abs_series) = valuation(a) == 0 && isunit(coeff(a, 0))

function isone(a::fq_abs_series)
   return precision(a) == 0 || ccall((:fq_poly_is_one, :libflint), Bool,
                   (Ref{fq_abs_series}, Ref{FqFiniteField}), a, base_ring(a))
end

# todo: write an fq_poly_valuation
function valuation(a::fq_abs_series)
   for i = 1:length(a)
      if !iszero(coeff(a, i - 1))
         return i - 1
      end
   end
   return precision(a)
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FqAbsSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fq_abs_series}) = show_minus_one(fq)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fq_abs_series)
   z = parent(x)()
   ccall((:fq_poly_neg, :libflint), Void,
                (Ref{fq_abs_series}, Ref{fq_abs_series}, Ref{FqFiniteField}),
               z, x, base_ring(x))
   z.prec = x.prec
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fq_abs_series, b::fq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   prec = min(a.prec, b.prec)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fq_poly_add_series, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(a))
   return z
end

function -(a::fq_abs_series, b::fq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   prec = min(a.prec, b.prec)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fq_poly_sub_series, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(a))
   return z
end

function *(a::fq_abs_series, b::fq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(a)))
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.prec = prec
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:fq_poly_mullow, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(a))
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::fq, y::fq_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fq_poly_scalar_mul_fq, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series}, Ref{fq}, Ref{FqFiniteField}),
               z, y, x, base_ring(y))
   return z
end

*(x::fq_abs_series, y::fq) = y * x

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fq_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   z.prec = x.prec + len
   ccall((:fq_poly_shift_left, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, x, len, base_ring(x))
   return z
end

function shift_right(x::fq_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   if len >= xlen
      z.prec = max(0, x.prec - len)
   else
      z.prec = x.prec - len
      ccall((:fq_poly_shift_right, :libflint), Void,
            (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, x, len, base_ring(x))
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fq_abs_series, prec::Int)
   prec < 0 && throw(DomainError())
   if x.prec <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   ccall((:fq_poly_set_trunc, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, x, prec, base_ring(x))
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fq_abs_series, b::Int)
   b < 0 && throw(DomainError())
   if precision(a) > 0 && isgen(a) && b > 0
      return shift_left(a, b - 1)
   elseif length(a) == 1
      return parent(a)([coeff(a, 0)^b], 1, a.prec)
   elseif b == 0
      z = one(parent(a))
      set_prec!(z, precision(a))
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit !=0
         z = z*z
         if (UInt(bit) & b) != 0
            z *= a
         end
         bit >>= 1
      end
   end
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fq_abs_series, y::fq_abs_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   n = max(length(x), length(y))
   n = min(n, prec)
   return Bool(ccall((:fq_poly_equal_trunc, :libflint), Cint,
             (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               x, y, n, base_ring(x)))
end

function isequal(x::fq_abs_series, y::fq_abs_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || length(x) != length(y)
      return false
   end
   return Bool(ccall((:fq_poly_equal, :libflint), Cint,
             (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               x, y, length(x), base_ring(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fq_abs_series, y::fq)
   if length(x) > 1
      return false
   elseif length(x) == 1
      z = base_ring(x)()
      ccall((:fq_poly_get_coeff, :libflint), Void,
            (Ref{fq}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
             z, x, 0, base_ring(x))
      return z == y
   else
      return precision(x) == 0 || iszero(y)
   end
end

==(x::fq, y::fq_abs_series) = y == x

function ==(x::fq_abs_series, y::fmpz)
   if length(x) > 1
      return false
   elseif length(x) == 1
      z = base_ring(x)()
      ccall((:fq_poly_get_coeff, :libflint), Void,
            (Ref{fq}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
             z, x, 0, base_ring(x))
      return z == y
   else
      return precision(x) == 0 || iszero(y)
   end
end

==(x::fmpz, y::fq_abs_series) = y == x

==(x::fq_abs_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::fq_abs_series) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fq_abs_series, y::fq_abs_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   v2 = valuation(y)
   v1 = valuation(x)
   if v2 != 0
      if v1 >= v2
         x = shift_right(x, v2)
         y = shift_right(y, v2)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec, y.prec - v2 + v1)
   z = parent(x)()
   z.prec = prec
   ccall((:fq_poly_div_series, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, x, y, prec, base_ring(x))
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fq_abs_series, y::fq)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fq_poly_scalar_div_fq, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series}, Ref{fq}, Ref{FqFiniteField}),
               z, x, y, base_ring(x))
   return z
end

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fq_abs_series)
   a == 0 && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ccall((:fq_poly_inv_series, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               ainv, a, a.prec, base_ring(a))
   return ainv
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function fit!(z::fq_abs_series, n::Int)
   ccall((:fq_poly_fit_length, :libflint), Void,
         (Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
         z, n, base_ring(z))
   return nothing
end

function setcoeff!(z::fq_abs_series, n::Int, x::fq)
   ccall((:fq_poly_set_coeff, :libflint), Void,
                (Ref{fq_abs_series}, Int, Ref{fq}, Ref{FqFiniteField}),
               z, n, x, base_ring(z))
   return z
end

function mul!(z::fq_abs_series, a::fq_abs_series, b::fq_abs_series)
   lena = length(a)
   lenb = length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(z)))
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = min(lena + lenb - 1, prec)
   if lenz < 0
      lenz = 0
   end
   z.prec = prec
   ccall((:fq_poly_mullow, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               z, a, b, lenz, base_ring(z))
   return z
end

function addeq!(a::fq_abs_series, b::fq_abs_series)
   lena = length(a)
   lenb = length(b)
   prec = min(a.prec, b.prec)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   lenz = max(lena, lenb)
   a.prec = prec
   ccall((:fq_poly_add_series, :libflint), Void,
         (Ref{fq_abs_series}, Ref{fq_abs_series},
          Ref{fq_abs_series}, Int, Ref{FqFiniteField}),
               a, a, b, lenz, base_ring(a))
   return a
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fq_abs_series}, ::Type{T}) where {T <: Integer} = fq_abs_series

promote_rule(::Type{fq_abs_series}, ::Type{fq}) = fq_abs_series

promote_rule(::Type{fq_abs_series}, ::Type{fmpz}) = fq_abs_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FqAbsSeriesRing)()
   ctx = base_ring(a)
   z = fq_abs_series(ctx)
   z.prec = a.prec_max
   z.parent = a
   return z
end

function (a::FqAbsSeriesRing)(b::Integer)
   ctx = base_ring(a)
   if b == 0
      z = fq_abs_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_abs_series(ctx, [base_ring(a)(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FqAbsSeriesRing)(b::fmpz)
   ctx = base_ring(a)
   if b == 0
      z = fq_abs_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_abs_series(ctx, [base_ring(a)(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FqAbsSeriesRing)(b::fq)
   ctx = base_ring(a)
   if b == 0
      z = fq_abs_series(ctx)
      z.prec = a.prec_max
   else
      z = fq_abs_series(ctx, [b], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FqAbsSeriesRing)(b::fq_abs_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FqAbsSeriesRing)(b::Array{fq, 1}, len::Int, prec::Int)
   ctx = base_ring(a)
   z = fq_abs_series(ctx, b, len, prec)
   z.parent = a
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

function PowerSeriesRing(R::FqFiniteField, prec::Int, s::AbstractString; model=:capped_relative, cached = true)
   S = Symbol(s)

   if model == :capped_relative
      parent_obj = FqRelSeriesRing(R, prec, S, cached)
   elseif model == :capped_absolute
      parent_obj = FqAbsSeriesRing(R, prec, S, cached)
   end

   return parent_obj, gen(parent_obj)
end


text:
# Plot coordinates of the C alpha with best occupancy
@recipe function plot(residues::AbstractVector{PDBResidue})
    ca = CAmatrix(residues)
    ca[:,1], ca[:,2], ca[:,3]
end

ext: jl


text:
# Inverted Pendulum task for continous reinforcement learning as describe in XXX

@with_kw mutable struct InvertedPendulum <: MDP{Tuple{Float64,Float64},Float64}
    g::Float64 = 9.81
    m::Float64 = 2.
    l::Float64 = 8.
    M::Float64 = 0.5
    alpha::Float64 = 1/(m + M)
    dt::Float64 = 0.1
    discount::Float64 = 0.99
    cost::Float64 = -1.
end

actions(ip::InvertedPendulum) = [-50., 0., 50.]

initialstate(ip::InvertedPendulum) = ImplicitDistribution(rng -> ((rand(rng)-0.5)*0.1, (rand(rng)-0.5)*0.1, ))

function reward(ip::InvertedPendulum,
              s::Tuple{Float64,Float64},
              a::Float64,
              sp::Tuple{Float64,Float64})
    return isterminal(ip, sp) ? ip.cost : 0.0
end

discount(ip::InvertedPendulum) = ip.discount
isterminal(::InvertedPendulum, s::Tuple{Float64,Float64}) = abs(s[1]) > pi/2.

function dwdt(m::InvertedPendulum,th::Float64,w::Float64,u::Float64)
    num = m.g*sin(th)-m.alpha*m.m*m.l*(w^2)*sin(2*th)*0.5 - m.alpha*cos(th)*u
    den = (4/3)*m.l - m.alpha*m.l*(cos(th)^2)
    return num/den
end

# TODO
# function rk45(m::InvertedPendulum,s::Tuple{Float64,Float64},a::Float64)
#     k1 = dwdt(m,s[1],s[2],a)
#     #something...
# end

function euler(m::InvertedPendulum,s::Tuple{Float64,Float64},a::Float64)
    alph = dwdt(m,s[1],s[2],a)
    w_ = s[2] + alph*m.dt
    th_ = s[1] + s[2]*m.dt + 0.5*alph*m.dt^2
    if th_ > pi
        th_ -= 2*pi
    elseif th_ < -pi
        th_ += 2*pi
    end
    return (th_,w_)
end

function transition(ip::InvertedPendulum, s::Tuple{Float64,Float64}, a::Float64)
    ImplicitDistribution(s, a) do s, a, rng
        a_offset = 20*(rand(rng)-0.5)
        a_ = a + a_offset
        return euler(ip, s, a_)
    end
end

function convert_s(::Type{A}, s::Tuple{Float64,Float64}, ip::InvertedPendulum) where A<:AbstractArray
    v = copyto!(Array{Float64}(undef, 2), s)
    return v
end

function convert_s(::Type{Tuple{Float64,Float64}}, s::A, ip::InvertedPendulum) where A<:AbstractArray
    return (s[1], s[2])
end

# deprecated in POMDPs v0.9
POMDPs.initialstate_distribution(ip::InvertedPendulum) = ImplicitDistribution(rng -> ((rand(rng)-0.5)*0.1, (rand(rng)-0.5)*0.1, ))


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdstock

"""
    ''_antdstock(;kwargs...)

An AntdStock component.

Keyword arguments:
- `id` (String; optional)
- `annotations` (optional)
- `appendPadding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `data` (Array of Dict with Strings as keys and values of type String | Reals; required)
- `fallingFill` (String; optional)
- `height` (Real; optional)
- `label` (optional)
- `legend` (optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-CN', 'en-US'; optional)
- `meta` (optional)
- `padding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `renderer` (a value equal to: 'canvas', 'svg'; optional)
- `risingFill` (String; optional)
- `slider` (optional)
- `stockStyle` (optional): . stockStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `style` (Dict; optional)
- `tooltip` (optional)
- `width` (Real; optional)
- `xAxis` (optional)
- `xField` (String; required)
- `yAxis` (optional)
- `yField` (Array of Strings; required)
"""
function ''_antdstock(; kwargs...)
        available_props = Symbol[:id, :annotations, :appendPadding, :autoFit, :className, :data, :fallingFill, :height, :label, :legend, :loading_state, :locale, :meta, :padding, :renderer, :risingFill, :slider, :stockStyle, :style, :tooltip, :width, :xAxis, :xField, :yAxis, :yField]
        wild_props = Symbol[]
        return Component("''_antdstock", "AntdStock", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Entwine_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Entwine")
JLLWrappers.@generate_main_file("Entwine", UUID("9290252f-3c82-5b4b-93a0-1879861c8c5d"))
end  # module Entwine_jll


text:
using Cbc
using PowerSimulations
using PowerSystems
using DataStructures
using InfrastructureSystems
const IS = InfrastructureSystems
const PSI = PowerSimulations
const PSY = PowerSystems

include("../../../test/test_utils/get_test_data.jl")

abstract type TestOpProblem <: PSI.DecisionProblem end

system = build_c_sys5_re(; add_reserves=true)
solver = optimizer_with_attributes(Cbc.Optimizer)

devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalBasicDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticBranch),
    :TT => DeviceModel(TapTransformer, StaticBranch),
);
services = Dict{Symbol, ServiceModel}();

template = PSI.ProblemTemplate(CopperPlatePowerModel, devices, branches, services);

operation_problem = PSI.DecisionModel(TestOpProblem, template, system; optimizer=solver);

set_services_template!(
    operation_problem,
    Dict(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :Down_Reserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
    ),
)

op_results = solve!(operation_problem)
re_results = PSI.run_economic_dispatch(system; optimizer=solver, use_parameters=true)


text:
# developed with Julia 1.1.1
#
# functions to calibrate EMS control models

using EMSx

using JLD2
using ProgressMeter

function calibrate_sites(controller::EMSx.AbstractController,
    path_to_save_folder::String,
    path_to_price_csv_file::String,
    path_to_metadata_csv_file::String,
    path_to_train_data_folder::String)

    EMSx.make_directory(joinpath(path_to_save_folder, "value_functions"))
    prices = EMSx.load_prices(path_to_price_csv_file)
    sites = EMSx.load_sites(path_to_metadata_csv_file, nothing,
        path_to_train_data_folder, path_to_save_folder)

    elapsed = 0.0

    @showprogress for site in sites

        elapsed += @elapsed site_value_functions = calibrate_site(controller,
            site, prices)

    end

    println("Terminating model calibration in $(elapsed) seconds")

    return nothing

end

function calibrate_site(controller::EMSx.AbstractController, site::EMSx.Site,
    prices::EMSx.Prices)

    controller = EMSx.initialize_site_controller(controller, site, prices)

    timer = @elapsed value_function = compute_value_functions(controller)

    JLD2.save(joinpath(site.path_to_save_folder, "value_functions", site.id*".jld2"),
        Dict("value_function"=>value_function, "time"=>timer))

    return nothing

end


text:
using Documenter, AlphaVantage

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs"),
    repo   = "github.com/ellisvalentiner/AlphaVantage.jl.git",
    julia  = "0.6"
)


text:
## Initialize
import Pkg
Pkg.activate("/Users/arreyes/.julia/environments/v1.5/")
using Plots, BenchmarkTools
gr()
Pkg.activate("/Users/arreyes/.julia/dev/RamirezReyes_Yang_SpontaneousCyclogenesis/")
using AvailablePotentialEnergyFramework
using Statistics:mean
using NCDatasets
import ImageSegmentation: seeded_region_growing, labels_map, region_adjacency_graph
## Read test file
test_file_with_cyclones = "/Users/arreyes/Documents/Research/Developement/testNetCDF/f5e-4_2km_1000km_homoRad_2d_timesteps.nc"
test_file_with_cyclones_3d = "/Users/arreyes/Documents/Research/Developement/testNetCDF/f5e-4_2km_1000km_homoRad_3d_timesteps.nc"
ds = Dataset(test_file_with_cyclones);
ds_3d = Dataset(test_file_with_cyclones_3d);
pres = ds_3d["PP"][:,:,1,end] :: Array{Float32,2};
U = ds_3d["U"][:,:,:,end] :: Array{Float32,3};
V = ds_3d["V"][:,:,:,end] :: Array{Float32,3};
close(ds)
close(ds_3d)
##Algorithm has to be: Detect cyclone, recenter, detect again and recenter. It has to be done twice to avoid centering some areas that are on the edge of a cyclone and are still detected as a local maximum

windspeed = sqrt.(U.*U .+ V.*V);
pres_anomaly = pres .- mean(pres,dims=(1,2));

heatmap(permutedims(windspeed[:,:,1]),title="Wind speed field (m/s)")
p1 = heatmap(permutedims(pres_anomaly))
centers = AvailablePotentialEnergyFramework.findcyclonecenters_aspressureminima(pres_anomaly,-5,2000);
centers_x= [centers[i][1] for i in 1:length(centers)];
centers_y = [centers[i][2] for i in 1:length(centers)];
scatter!(p1,centers_x,centers_y,markersize=10,markershape=:cross,c=:red)

centers_and_labels,cyclones = AvailablePotentialEnergyFramework.detect_cyclones(pres,-5,2000)
labels_map2 = copy(labels_map(cyclones));
labels_map2[labels_map2 .== 1000] .= 0;
p2 = heatmap(permutedims(labels_map2),reuse=false,title="Identified and segmented cyclones")

G, vert_map = region_adjacency_graph(cyclones, (i,j)->true);

radiusbins = 1000:2000:512000;

presaverages = []
for i in 1:(length(radiusbins) - 1 )
push!(presaverages,
AvailablePotentialEnergyFramework.averageallindistance(radiusbins[i:i+1],pres,labels_map2.==1,centers_and_labels[1][1],2000))
end

plot(radiusbins[1:end-1]./1e3,presaverages./100,title="Azimuthally averaged surface pressure anomaly (hPa)")

windspeedaverages = AvailablePotentialEnergyFramework.averageallindistance(radiusbins[1:2],windspeed,labels_map2.==2,centers_and_labels[2][1],2000)
for i in 2:(length(radiusbins) - 1 )
    global windspeedaverages
    windspeedaverages = hcat(windspeedaverages,
        AvailablePotentialEnergyFramework.averageallindistance(radiusbins[i:i+1],windspeed,labels_map2.==2,centers_and_labels[2][1],2000))
end
heatmap(windspeedaverages,title = "Azimuthally averaged wind speed for a random cyclone (m/s)")
size(windspeedaverages)

presaverage = AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,pres,cyclones,centers_and_labels,2000)
windspeedaverage = AvailablePotentialEnergyFramework.azimuthalaverage_allcyclones(radiusbins,windspeed,cyclones,centers_and_labels,2000)
plot(radiusbins[1:end-1]./1e3,1:80,presaverage./100,title="Azimuthally averaged surface pressure anomaly (hPa) 3 cyclones 1 timestep")
heatmap(radiusbins[1:end-1]./1e3,1:80,windspeedaverage,title="Azimuthally averaged surface windspeed (m/s) 3 cyclones 1 timestep")


text:
export DummyRobot, DummyRobotParams

Base.@kwdef struct DummyRobotParams <: DeviceParams
  test::String = "" # WARNING: this is needed since otherwise the constructor cannot be called
  namedPositions::Dict{String, Vector{typeof(1.0u"mm")}} = Dict("origin" => [0,0,0]u"mm")
  coordinateSystem::ScannerCoordinateSystem = ScannerCoordinateSystem(3)
end

DummyRobotParams(dict::Dict) = params_from_dict(DummyRobotParams, prepareRobotDict(dict))

Base.@kwdef mutable struct DummyRobot <: Robot
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::DummyRobotParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}
  state::RobotState = INIT
  referenced::Bool = false
end

_isReferenced(rob::DummyRobot) = rob.referenced

_getPosition(rob::DummyRobot) = [1.0,0.0,0.0]u"mm"
dof(rob::DummyRobot) = 3
axisRange(rob::DummyRobot) = [[-Inf, Inf], [-Inf, Inf], [-Inf, Inf]]u"m"

function _moveAbs(rob::DummyRobot, pos::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  @debug "DummyRobot: movAbs pos=$pos"
end

function _moveRel(rob::DummyRobot, dist::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  @debug "DummyRobot: movRel dist=$dist"
end

function _enable(rob::DummyRobot)
  @debug "DummyRobot: enable"
end

function _disable(rob::DummyRobot)
  @debug "DummyRobot: disable"
end

function _reset(rob::DummyRobot)
  @debug "DummyRobot: reset"
end

function _setup(rob::DummyRobot)
  @debug "DummyRobot: setup"
end

function _doReferenceDrive(rob::DummyRobot)
  @debug "DummyRobot: Doing reference drive"
  rob.referenced = true
end


text:
include("../Utils.jl")

@doc """
An object of this class represents a point in parameter space.
There are functions defined to evaluate the log likelihood and
move around.
""" ->
type Particle
	params::Vector{Float64}
end

@doc """
A constructor. Makes params have length 20
""" ->
function Particle()
	return Particle(Array(Float64, (20, )))
end

@doc """
Generate params from the prior
""" ->
function from_prior!(particle::Particle)
	particle.params = -0.5 + rand(length(particle.params))
	return nothing
end

@doc """
Do a metropolis proposal. Return log(hastings factor for prior sampling)
""" ->
function perturb!(particle::Particle)
	i = rand(1:length(particle.params))
	particle.params[i] += randh()
	particle.params[i] = mod(particle.params[i] + 0.5, 1.0) - 0.5
	return 0.0
end

@doc """
Evaluate the log likelihood
""" ->
function log_likelihood(particle::Particle)
	logL1 = -length(particle.params)*0.5*log(2*pi*0.1^2)
	logL2 = -length(particle.params)*0.5*log(2*pi*0.01^2)
	for i in 1:length(particle.params)
		logL1 += -0.5*(particle.params[i]/0.1)^2
		logL2 += -0.5*((particle.params[i] - 0.031)/0.01)^2
	end
	return logsumexp([logL1, logL2 + log(100.0)])
end

@doc """
Convert to string, for output to sample.txt
"""
import Base.string
function string(particle::Particle)
	return join([string(signif(x, 6), " ") for x in particle.params])
end


text:
using EchelleInstruments
using Documenter

makedocs(;
    modules=[EchelleInstruments],
    authors="Eric Ford",
    repo="https://github.com/RvSpectML/EchelleInstruments.jl/blob/{commit}{path}#L{line}",
    sitename="EchelleInstruments.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RvSpectML.github.io/EchelleInstruments.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Instrument agnostic" => "agnostic.md",
        "Instruments" => [
            "NEID" => "neid.md",
            "EXPRES" => "expres.md"],
        "Index" => "longlist.md",
    ],
)

deploydocs(;
    repo="github.com/RvSpectML/EchelleInstruments.jl",
)

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
"""-----------------------------------------------------------------------------
    contracts a tensor represented with indices + values in a symmetric COO
  format. Rows of the indices correspond to sorted indices of a symmetric
  hyper-edge. Values correspond to the edge weights, and v is vector being
  contracted with the
-----------------------------------------------------------------------------"""
function tensor_vector_contraction(A::ThirdOrderSymTensor, x::Array{Float64,1})
    @assert length(x) == A.n

    y = zeros(Float64,A.n)
    edge_count = size(A.indices,1)
	j =1
    @inbounds for (i_1,i_2,i_3) in eachrow(A.indices)

        y[i_1] += 2 * x[i_2]*x[i_3]*A.values[j]
        y[i_2] += 2 * x[i_1]*x[i_3]*A.values[j]
        y[i_3] += 2 * x[i_1]*x[i_2]*A.values[j]

		j +=1
    end

    return y
end

function tensor_vector_contraction(A::UnweightedThirdOrderSymTensor, x::Array{Float64,1})
    @assert length(x) == A.n

    y = zeros(Float64,A.n)

    @inbounds for (i_1,i_2,i_3) in eachrow(A.indices)

        y[i_1] += 2 * x[i_2]*x[i_3]
        y[i_2] += 2 * x[i_1]*x[i_3]
        y[i_3] += 2 * x[i_1]*x[i_2]

    end

    return y
end

"""-----------------------------------------------------------------------------
    contracts a tensor represented with indices + values in a symmetric COO
  format with a vector to return a sparse matrix.
-----------------------------------------------------------------------------"""
function tri_sub_tensor(A::ThirdOrderSymTensor,x::Array{Float64,1})

	index = 0
	nnzs = size(A.indices,1)
	I = zeros(Int64,6*nnzs)
	J = zeros(Int64,6*nnzs)
	V = zeros(Float64,6*nnzs)

    @inbounds for nnz_index = 1:nnzs
		index += 6
		i,j,k = A.indices[nnz_index,:]
		val = A.values[nnz_index]

		I[index-5] = i
		J[index-5] = j
		I[index-4] = j
		J[index-4] = i
		V[index-5] = val*x[k]
		V[index-4] = val*x[k]

		I[index-3] = i
		J[index-3] = k
		I[index-2] = k
		J[index-2] = i
		V[index-3] = val*x[j]
		V[index-2] = val*x[j]

		I[index-1] = j
		J[index-1] = k
		I[index] = k
		J[index] = j
		V[index-1] = val*x[i]
		V[index] = val*x[i]

	end
	return sparse(I,J,V,A.n,A.n)
end

"""-----------------------------------------------------------------------------
    contracts a tensor represented with indices + values in a symmetric COO
  format with a vector to return a sparse matrix.
-----------------------------------------------------------------------------"""
function tri_sub_tensor(A::ThirdOrderSymTensor,x::Array{Float64,1},
    					I::Array{Int64,1},J::Array{Int64,1},V::Array{Float64,1})

	index = 0
	nnzs = size(A.indices,1)
    for nnz_index = 1:nnzs
		index += 6
		i,j,k = A.indices[nnz_index,:]
		val = A.values[nnz_index]

		I[index-5] = i
		J[index-5] = j
		I[index-4] = j
		J[index-4] = i
		V[index-5] = val*x[k]
		V[index-4] = val*x[k]

		I[index-3] = i
		J[index-3] = k
		I[index-2] = k
		J[index-2] = i
		V[index-3] = val*x[j]
		V[index-2] = val*x[j]

		I[index-1] = j
		J[index-1] = k
		I[index] = k
		J[index] = j
		V[index-1] = val*x[i]
		V[index] = val*x[i]

	end
	return sparse(I,J,V,A.n,A.n)
end

"""-----------------------------------------------------------------------------
    contracts a tensor represented with indices + values in a symmetric COO
  format with a vector to return a sparse matrix.
-----------------------------------------------------------------------------"""
function tri_sub_tensor(A::UnweightedThirdOrderSymTensor,x::Array{Float64,1})

	index = 0
	nnzs = 2*sum([length(x) for x in A_U_ten.indices])
	I = zeros(Int64,2*nnzs)
	J = zeros(Int64,2*nnzs)
	V = zeros(Float64,2*nnzs)

	for i in 1:A.n
		for (j,k) in A.indices[i]
			index += 2

			I[index-1] = j
			J[index-1] = k
			I[index] = k
			J[index] = j
			V[index-1] = x[i]
			V[index] = x[i]
		end
	end
	return sparse(I[1:index],J[1:index],V[1:index],A.n,A.n)
end

function kron_contract(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,
                       X::SparseMatrixCSC{Float64,Int64},rank::Int=minimum(size(X)))

	(U,S,VT),_ = svds(X,nsv = rank)
	println("input rank: ",rank)
	println("output rank:",length(S))

	#singular_indexes = [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*minimum(size(X))]

	#println("rank in kron contract is",length(singular_indexes))

    kron_contract(A,B,U,VT*diagm(S))
end

function kron_contract(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,
                       U::Array{Float64,2},V::Array{Float64,2})
    n,d1 = size(U)
    m,d2 = size(V)
    @assert d1 == d2
    @assert A.n == n
    @assert B.n == m

	#preallocate memory used for building sparse matrices
	A_nnzs = size(A.indices,1)
	A_I = zeros(Int64,6*A_nnzs)
	A_J = zeros(Int64,6*A_nnzs)
	A_V = zeros(Float64,6*A_nnzs)

	B_nnzs = size(B.indices,1)
	B_I = zeros(Int64,6*B_nnzs)
	B_J = zeros(Int64,6*B_nnzs)
	B_V = zeros(Float64,6*B_nnzs)

    result = zeros(A.n * B.n)

    @inbounds for i in 1:d1

	#	sub_A_i = tri_sub_tensor(A,U[:,i])
	#	sub_B_i = tri_sub_tensor(B,V[:,i])
        sub_A_i = tri_sub_tensor(A,U[:,i],A_I,A_J,A_V)
        sub_B_i = tri_sub_tensor(B,V[:,i],B_I,B_J,B_V)
        for j in 1:i

            A_update = (sub_A_i*U[:,j])
            B_update = (sub_B_i*V[:,j])

            if i == j
               result += kron(B_update,A_update)
            else
                result += kron(2*B_update,A_update)
            end
        end
    end

    return reshape(result, A.n,B.n)
end

function get_kron_contract_comps(A::Union{ThirdOrderSymTensor,UnweightedThirdOrderSymTensor},
    							 B::Union{ThirdOrderSymTensor,UnweightedThirdOrderSymTensor},
              			         U::Array{Float64,2},V::Array{Float64,2})

    n,d1 = size(U)
    m,d2 = size(V)
    @assert d1 == d2
    @assert A.n == n
    @assert B.n == m

#    result = zeros(A.n * B.n)

    ileave = (i,j) -> i + A.n*(j-1)

    max_rank = Int((d1+1)*d1/2)

	A_comps = zeros(n,max_rank)
	B_comps = zeros(m,max_rank)
	index = 1

	A_nnzs = size(A.indices,1)
	A_I = zeros(Int64,6*A_nnzs)
	A_J = zeros(Int64,6*A_nnzs)
	A_V = zeros(Float64,6*A_nnzs)

	B_nnzs = size(B.indices,1)
	B_I = zeros(Int64,6*B_nnzs)
	B_J = zeros(Int64,6*B_nnzs)
	B_V = zeros(Float64,6*B_nnzs)

    @inbounds for i in 1:d1

        sub_A_i = tri_sub_tensor(A,U[:,i],A_I,A_J,A_V)
        sub_B_i = tri_sub_tensor(B,V[:,i],B_I,B_J,B_V)

        for j in 1:i

			if i == j
				A_comps[:,index] = (sub_A_i*U[:,j])
				B_comps[:,index] = (sub_B_i*V[:,j])
			else

				A_comps[:,index] = sqrt(2)*(sub_A_i*U[:,j])
				B_comps[:,index] = sqrt(2)*(sub_B_i*V[:,j])
			end
			index += 1
        end
    end

    return A_comps, B_comps
end

function implicit_contraction(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,x::Array{Float64,1})

	error("Code is Broken, must be fixed")
    @assert length(x) == A.n*B.n
    m = A.n
    n = B.n
	X = reshape(x,A.n,B.n)
    Y = similar(X)

    #ileave = (i,j) -> i + m*(j-1)
	i = 1
	j = 1
    @inbounds for (i_1,i_2,i_3) in eachrow(A.indices)

		A_val = A.values[i]
		i += 1
		j = 1
#        for (i_1,i_2,i_3) in permutations(A.indices[i,:])
		for (j_1,j_2,j_3) in eachrow(B.indices)

			#i_1,i_2,i_3
			@inbounds Y[i_1,j_1] += 2*A_val*B.values[j]*X[i_2,j_2]*X[i_3,j_3]
			@inbounds Y[i_2,j_2] += 2*A_val*B.values[j]*X[i_1,j_1]*X[i_3,j_3]
			@inbounds Y[i_3,j_3] += 2*A_val*B.values[j]*X[i_1,j_1]*X[i_2,j_2]

			#i_1,i_3,i_2
			@inbounds Y[i_1,j_1] += 2*A_val*B.values[j]*X[i_3,j_2]*X[i_2,j_3]
			@inbounds Y[i_3,j_2] += 2*A_val*B.values[j]*X[i_1,j_1]*X[i_2,j_3]
			@inbounds Y[i_2,j_3] += 2*A_val*B.values[j]*X[i_1,j_1]*X[i_3,j_2]

			#i_2,i_1,i_3
			@inbounds Y[i_2,j_1] += 2*A_val*B.values[j]*X[i_1,j_2]*X[i_3,j_3]
			@inbounds Y[i_1,j_2] += 2*A_val*B.values[j]*X[i_2,j_1]*X[i_3,j_3]
			@inbounds Y[i_3,j_3] += 2*A_val*B.values[j]*X[i_2,j_1]*X[i_1,j_2]

			#i_2,i_3,i_1
			@inbounds Y[i_2,j_1] += 2*A_val*B.values[j]*X[i_3,j_2]*X[i_1,j_3]
			@inbounds Y[i_3,j_2] += 2*A_val*B.values[j]*X[i_2,j_1]*X[i_1,j_3]
			@inbounds Y[i_1,j_3] += 2*A_val*B.values[j]*X[i_2,j_1]*X[i_3,j_2]

			#i_3,i_2,i_1
			@inbounds Y[i_3,j_1] += 2*A_val*B.values[j]*X[i_2,j_2]*X[i_1,j_3]
			@inbounds Y[i_2,j_2] += 2*A_val*B.values[j]*X[i_3,j_1]*X[i_1,j_3]
			@inbounds Y[i_1,j_3] += 2*A_val*B.values[j]*X[i_3,j_1]*X[i_2,j_2]

			#i_3,i_1,i_2
			@inbounds Y[i_3,j_1] += 2*A_val*B.values[j]*X[i_1,j_2]*X[i_2,j_3]
			@inbounds Y[i_1,j_2] += 2*A_val*B.values[j]*X[i_3,j_1]*X[i_2,j_3]
			@inbounds Y[i_2,j_3] += 2*A_val*B.values[j]*X[i_3,j_1]*X[i_1,j_2]

		end
		j += 1
  #      end
    end

    return Y[:]
end

function impTTVnodesym(A::UnweightedThirdOrderSymTensor,
                       B::UnweightedThirdOrderSymTensor,x::Vector{Float64})

	X = reshape(x,A.n,B.n)
	Y = similar(X)
	Y .= 0
	@inbounds for a = 1:A.n
		for b = 1:B.n
			newval = 0.0
			for (jp,kp) in A.indices[a]
				@simd for pair in B.indices[b]
					j,k = pair
					@inbounds newval += X[jp,j]*X[kp,k]+X[kp,j]*X[jp,k]
				end
			end
			Y[a,b] = 2*newval
		end
	end
	y = Y[:]
	return y

end

function impTTVnodesym(nG::Int,nH::Int,x::Vector{Float64},Gti, Hti)

	X = reshape(x,nG,nH)
	#Xt = copy(X')
	Y = similar(X)
	Y .= 0
	@inbounds for g = 1:nG
		for h = 1:nH
			newval = 0.0
			for (jp,kp) in Gti[g]
				@simd for pair in Hti[h]
					j,k = pair
					@inbounds newval += X[jp,j]*X[kp,k]+X[kp,j]*X[jp,k]
				end
			end
			Y[g,h] = 2*newval
		end
	end
	y = Y[:]
	return y
end

#TODO: move to appropriate file
function SSHOPM(A::ThirdOrderSymTensor,β::Float64, max_iter::Int, tol::Float64,
	            x_0::Array{Float64,1}=ones(A.n);update_user::Int=-1)

	#normalize init vector
	x_k = copy(x_0)
	x_k ./=norm(x_0)

    lambda_k = Inf
    i = 1

    while true

		x_k_1 = tensor_vector_contraction(A,x_k)
		lambda_k_1 = x_k_1'*x_k

		if β != 0.0
			x_k_1 .+= β*x_k
        end

		x_k_1 ./= norm(x_k_1)

		res = abs(lambda_k_1 - lambda_k)
		if update_user != -1 && i % update_user == 0
			println("iteration $(i)    λ: $(lambda_k_1)   |λ_k_1 - λ_k| = $res")
		end

        if res < tol || i >= max_iter
            return x_k_1, lambda_k_1
        else
			lambda_k = copy(lambda_k_1)
			x_k = copy(x_k_1)
            i += 1
        end

    end

end


text:
# -*- coding: utf-8 -*-
include("gerda_data_processing/init.jl")
include("gerda_data_processing/fct.jl")
include("fct-sim.jl")

sim_file        = "/remote/ceph2/group/gerda/data/simulation/gerda-mage-sim/gedet/intrinsic_bege/0nbb/edep/raw-gedet-intrinsic_bege-0nbb-edep-000.root"
config_dir      = "/res/Impurity_Scan/config-dep/"
plots_base_path = "pulses/sim/raw_0vbb/plots/"
base_path_raw   = "pulses/sim/raw_0vbb/"
base_path       = "pulses/sim/0vbb/"
n_sim_events    = 2500
sampling_time   = 1u"ns"
hits_threshold  = 0.005; # MeV
generate_wf    = false;
generate_cl_wf = true;

function simulate_waveforms_safe(events, output_filename, sim)
    wf_p = []
    wf_n = []
    skip = []
    no_skip = []
    @showprogress for evt in events
        edep =  evt.edep
        pos = []
        for sp in evt.pos
            push!(pos, T[sp[1].val, sp[2].val, sp[3].val] ./ 1000)
        end
        pos = CartesianPoint.(pos)
        e = Event(pos, edep)
        try simulate!(e, sim)
            catch
            push!(skip, findfirst(x->x == evt, events))
            continue
        end
        push!(no_skip, findfirst(x->x == evt, events))
        push!(wf_p, e.waveforms[1])
        push!(wf_n, e.waveforms[2])
    end
#     return events
    tt = Table(
        evtno = Int32.(vcat(events.evtno[no_skip], events.evtno[no_skip])),
        event_mult = Int32.(vcat(events.event_mult[no_skip], events.event_mult[no_skip])),
        detno = VectorOfArrays(vcat(events.detno[no_skip], events.detno[no_skip])),
        hits_totnum = Int32.(vcat(events.hits_totnum[no_skip], events.hits_totnum[no_skip])),
        edep = VectorOfArrays(vcat(events.edep[no_skip], events.edep[no_skip])),
        pos = VectorOfArrays(vcat(events.pos[no_skip], events.pos[no_skip])),
        waveform = ArrayOfRDWaveforms(Array{typeof(wf_p[1]), 1}(vcat(wf_p, wf_n)))
    );
    if length(skip) > 0
        output_filename = split(output_filename, ".h5")[1] * "_skipped-" * string(length(skip)) * ".h5"
    end
    HDF5.h5open(output_filename, "w") do h5f
        LegendHDF5IO.writedata(h5f, "generated_waveforms", tt)
    end
end

# Search for available channels
channels = []
for c in 6:1:12
    files = glob("pulses/sim/raw_0vbb/raw*/"*string(c)*"-"*channel_to_name[c]*"/");
    if length(files) > 0 && channel_to_bege[c]
        push!(channels, c)
    end
end

function check_charge_deposits(positions, sim)
    good_event = true
    for pos in positions
        point = CartesianPoint{Float32}(pos[1].val, pos[2].val, pos[3].val) ./1000
        if !in(point, sim.detector.semiconductors[1].geometry)
            good_event = false
        end
    end
    return good_event
end

for ch in channels
    Base.run(`clear`)
    @info("Start ch"*string(ch)*" | "*channel_to_name[ch])
    @info(">---------------------------------------------------------------<")
    files = glob(base_path_raw*"*/*-"*channel_to_name[ch]*"/*.h5")
    p = Progress(length(files), dt=0.5,
                    barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                    barlen=10)
    for file in files
        # Check whether the file has already been simulated
        events = h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end

        last_file = glob(base_path_raw*"*"*channel_to_name[ch]*"/"*channel_to_name[ch]*"-cl-wf-*"*split(basename(file), ".")[1]*"_*")
        if length(last_file) == 0
            existance_check = false
        else
            last_file       = last_file[end]
            last_event      = parse(Int64,split(split(last_file, "_evts_")[2],"-")[end])
            existance_check = isfile(last_file) && last_event == size(events,1) && stat(last_file).size > 0
        end
        if existance_check == false
            # Read the SSD simulation file of depleted detector
            config_file = glob("res/Impurity_Scan/config-dep/"*channel_to_name[ch]*"*.config")[1]
            sim_output  = glob("res/Impurity_Scan/output/"*channel_to_name[ch]*"*/" )[1]
            # Open *.h5 file containing detector specific events
            sim = @suppress readcreateh5(config_file, sim_output);

            #
            ## Check for charge deposit location and filter if necessary
            t = events |> @filter(check_charge_deposits(_.pos, sim)) |> Table
            if size(events,1) != size(t,1)
                HDF5.h5open(file, "w") do h5f
                    LegendHDF5IO.writedata( h5f, "data", Table(evtno = t.evtno, event_mult = t.event_mult, detno = VectorOfArrays(t.detno), hits_totnum = t.hits_totnum,
                                                                            edep = VectorOfArrays(t.edep), pos = VectorOfArrays(t.pos)))
                end
                events = h5open(file, "r") do h5f
                    LegendHDF5IO.readdata(h5f, "data")
                end
            end

            isdir(base_path_raw*string(ch)*"-"*channel_to_name[ch]) ? "Directory exists" : mkpath(base_path_raw*string(ch)*"-"*channel_to_name[ch])
            file_id = findall(x->x == file, files)[1]
            output_dir       = base_path_raw*string(ch)*"-"*channel_to_name[ch]
            if generate_wf
                output_basename  = channel_to_name[ch]*"-wf-"*split(basename(file), ".")[1]
                i = 1
                j = n_sim_events
                while i <= size(events,1)
                    if j > size(events,1)
                        j = size(events,1)
                    end
                    step_rng = i:1:j
                    str_rng = lpad(i, 6, "0") * "-" * lpad(j, 6, "0")
                    output_filename = joinpath(output_dir, output_basename) * "_evts_" * str_rng * ".h5"
                    @info(str_rng)
                    try @suppress SSD.simulate_waveforms(events[step_rng], sim, output_dir, output_basename * "_evts_" * str_rng, chunk_n_physics_events = 10000, Δt = sampling_time);
                    catch
                        @suppress_err simulate_waveforms_safe(events[step_rng], output_filename)
                    end

                    i += n_sim_events
                    j += n_sim_events
                end
#                 t_with_waveforms = @suppress SSD.simulate_waveforms(events, sim, output_dir, output_basename, chunk_n_physics_events = n_sim_events, Δt = sampling_time);
            end

            if generate_cl_wf
                events_clustered = SSD.cluster_detector_hits(events, 0.2u"mm")
                output_basename  = channel_to_name[ch]*"-cl-wf-"*split(basename(file), ".")[1]
                i = 1
                j = n_sim_events
                while i <= size(events,1)
                    if j > size(events,1)
                        j = size(events,1)
                    end
                    step_rng = i:1:j
                    str_rng = lpad(i, 6, "0") * "-" * lpad(j, 6, "0")
                    output_filename = joinpath(output_dir, output_basename) * "_evts_" * str_rng * ".h5"

                    @info(str_rng)
                    try @suppress SSD.simulate_waveforms(events_clustered[step_rng], sim, output_dir, output_basename * "_evts_" * str_rng, chunk_n_physics_events = 10000, Δt = sampling_time);
                    catch
                        @info("Safe simulation")
                        @suppress_err simulate_waveforms_safe(events_clustered[step_rng], output_filename, sim)
                    end

                    i += n_sim_events
                    j += n_sim_events
                end

#                 t_clustered_with_waveforms = try @suppress SSD.simulate_waveforms(events_clustered, sim, output_dir, output_basename, chunk_n_physics_events = n_sim_events, Δt = sampling_time);
#                     catch; end
            end
        end
        next!(p)
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Dynare_preprocessor_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Dynare_preprocessor")
JLLWrappers.@generate_main_file("Dynare_preprocessor", UUID("8bf4d177-eddb-5a20-998e-55ecc2283fe8"))
end  # module Dynare_preprocessor_jll


text:
using BasicBlockRewriter
using Documenter

DocMeta.setdocmeta!(BasicBlockRewriter, :DocTestSetup, :(using BasicBlockRewriter); recursive=true)

makedocs(;
    modules=[BasicBlockRewriter],
    authors="Tim Holy <tim.holy@gmail.com> and contributors",
    repo="https://github.com/timholy/BasicBlockRewriter.jl/blob/{commit}{path}#{line}",
    sitename="BasicBlockRewriter.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://timholy.github.io/BasicBlockRewriter.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/timholy/BasicBlockRewriter.jl",
    devbranch="main",
)


text:
struct OrnsteinUhlenbeck{TB, Tβ, Ta} <: Evolution
    B::TB
    β::Tβ
    a::Ta
end
b(t, x, OU::OrnsteinUhlenbeck) = OU.B*x + OU.β
riccatti(t, P, OU::OrnsteinUhlenbeck) = OU.B*P + P*OU.B + OU.a

struct Euler{T} <: Evolution
    M::T
end

function evolve(M::Euler{<:OrnsteinUhlenbeck}, (s, u)::Pair, t)
    x, P = u
    t => Gaussian(x + b(s, x, M.M)*dt, P + riccatti(t, P, M.M)*dt)
end

function dyniterate(M::Euler{<:OrnsteinUhlenbeck}, (s, u)::Pair, t)
    x, P = u
    dub(t => Gaussian(x + b(s, x, M.M)*dt, P + riccatti(t, P, M.M)*dt))
end

function evolve(M::OrnsteinUhlenbeck, (s, u)::Pair, t)
    x, P = u
    t => Gaussian(x + b(s, x, M.M)*dt, P + riccatti(t, P, M.M)*dt)
end


text:
module CellIntegration

using Gridap
using Gridap.Kernels: IntegrateKernel

export integrate

function integrate(
  cellfun::CellMap{Point{D},1,T,N},
  phi::CellGeomap{D,Z},
  quad::CellQuadrature{D}) where {D,Z,T,N}
  z = coordinates(quad)
  w = weights(quad)
  f = evaluate(cellfun,z)
  j = evaluate(gradient(phi),z)
  m = meas(j)
  k = _prepare_kernel(f,m,w)
  apply(k,f,m,w)
end

function _prepare_kernel(f::CellArray{T,N},m::CellArray{J},w::CellArray{W}) where {T,N,J,W}
  S = Base._return_type(*,(T,J,W))
  IntegrateKernel(Val(N),S)
end

function integrate(
  cellfun::CellMap{Point{D},1},trian::Triangulation{D,Z},quad::CellQuadrature{D}) where {D,Z}
  phi = CellGeomap(trian)
  integrate(cellfun,phi,quad)
end

function integrate(
  fun::Function,trian::Triangulation{D,Z},quad::CellQuadrature{D}) where {D,Z}
  phi = CellGeomap(trian)
  cellfun = compose(fun,phi)
  integrate(cellfun,phi,quad)
end

end # module CellIntegration


text:
function encode(t::Tokenizer, doc::Vector)
    unk, space, lf = t.dict["UNKNOWN"], t.dict[" "], t.dict["\n"]
    chars = Int[]
    ranges = UnitRange{Int}[]
    pos = 1
    for sent in doc
        for (word,tag) in sent
            for c in tag
                c == '_' && continue
                if c == 'S' # space
                    push!(chars, space)
                elseif c == 'N' # newline
                    push!(chars, lf)
                end
                pos += 1
            end
            for c in word
                push!(chars, push!(t.dict,string(c)))
            end
            push!(ranges, pos:pos+length(word)-1)
            pos += length(word)
        end
    end
    chars, ranges
end

function flatten(data::Vector)
    res = Int[]
    for x in data
        append!(res, x)
    end
    res
end

function train(t::Tokenizer, nepochs::Int, traindata::Vector, testdata::Vector)
    function conv(data)
        chars, ranges = encode(t, data)
        tags = encode(t.tagset, ranges)
        # mini-batch
        chars_batch = [chars[i:min(i+1000-1,length(chars))] for i=1:1000:length(chars)]
        tags_batch = [tags[i:min(i+1000-1,length(tags))] for i=1:1000:length(tags)]
        chars_batch, tags_batch
    end
    train_x, train_y = conv(traindata)
    test_x, test_y = conv(testdata)

    opt = SGD(0.001)
    for epoch = 1:nepochs
        println("epoch: $(epoch)")
        loss = fit(train_x, train_y, t.model, crossentropy, opt)
        println("loss: $(loss)")

        test_z = map(test_x) do x
            argmax(t.model(x).data, 1)
        end
        y, z = flatten(test_y), flatten(test_z)
        ranges_y, ranges_z = decode(t.tagset, y), decode(t.tagset, z)
        fscore(ranges_y, ranges_z)

        acc = accuracy(flatten(test_y), flatten(test_z))

        println("test acc.: $(acc)")
        println("")
    end
end

function fscore(golds::Vector{UnitRange{Int}}, preds::Vector{UnitRange{Int}})
    g = Set(golds)
    p = Set(preds)
    inter = intersect(g, p)
    correct = length(inter)

    prec = correct / length(p)
    recall = correct / length(g)
    fvalue = 2*recall*prec/(recall+prec)
    println("precision: $(prec)")
    println("recall: $(recall)")
    println("f-value: $(fvalue)")
end

function accuracy(golds::Vector{Int}, preds::Vector{Int})
    @assert length(golds) == length(preds)
    correct = 0
    total = 0
    for i = 1:length(golds)
        golds[i] == preds[i] && (correct += 1)
        total += 1
    end
    correct / total
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
import LessOLS: normal_sampler, ols, show, equation

# create sample with known properties
ns = normal_sampler(a=0, b=20, β_0=5, β=[1, 2.5], sd_e=0.5)
sample1 = ns(100000)

# attempt a linear model
lm1 = ols(sample1, intercept=false)
show(lm1)

# check the coefficients are about right (for a model with zero intercept)
function almost_equals(a, b, threshold)
    diff = map(abs, (a-b))
    return all(diff .< threshold)
end
@assert almost_equals(lm1.beta, [1.2143, 2.713], 0.01)


text:
using AdvancedHMC,
    ForwardDiff,
    LinearAlgebra,
    MCMCDiagnosticTools,
    Optim,
    Pathfinder,
    Random,
    Statistics,
    StatsFuns,
    Test

Random.seed!(0)

struct RegressionModel{X,Y}
    x::X
    y::Y
end

function (m::RegressionModel)(θ)
    logσ = θ[1]
    σ = exp(logσ)
    α = θ[2]
    β = θ[3:end]
    x = m.x
    y = m.y
    lp = normlogpdf(σ) + logtwo
    lp += normlogpdf(α)
    lp += sum(normlogpdf, β)
    y_hat = muladd(x, β, α)
    lp += sum(eachindex(y_hat, y)) do i
        return normlogpdf(y_hat[i], σ, y[i])
    end
    return lp
end

function mean_and_mcse(f, θs)
    zs = map(f, θs)
    ms = mean(zs)
    ses = map(mcse, eachrow(reduce(hcat, zs)))
    return ms, ses
end

function compare_estimates(f, xs1, xs2, α=0.05)
    nparams = length(first(xs1))
    α /= nparams  # bonferroni correction
    p = α / 2
    m1, s1 = mean_and_mcse(f, xs1)
    m2, s2 = mean_and_mcse(f, xs2)
    zs = @. (m1 - m2) / sqrt(s1^2 + s2^2)
    @test all(norminvcdf(p) .< zs .< norminvccdf(p))
end

@testset "AdvancedHMC integration" begin
    A = Diagonal(rand(5))
    B = randn(5, 2)
    D = exp(Symmetric(randn(2, 2)))
    M⁻¹ = Pathfinder.WoodburyPDMat(A, B, D)
    A2 = Diagonal(rand(5))
    B2 = randn(5, 2)
    D2 = exp(Symmetric(randn(2, 2)))
    M⁻¹2 = Pathfinder.WoodburyPDMat(A2, B2, D2)
    ℓπ(x) = -sum(abs2, x) / 2
    ∂ℓπ∂θ(x) = -x

    @testset "RankUpdateEuclideanMetric" begin
        metric = Pathfinder.RankUpdateEuclideanMetric(M⁻¹)
        @test metric.M⁻¹ === M⁻¹
        metric_dense = AdvancedHMC.DenseEuclideanMetric(Symmetric(Matrix(M⁻¹)))
        h = AdvancedHMC.Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
        h_dense = AdvancedHMC.Hamiltonian(metric_dense, ℓπ, ∂ℓπ∂θ)
        r = randn(5)
        θ = randn(5)

        metric2 = Pathfinder.RankUpdateEuclideanMetric(3)
        @test metric2.M⁻¹ ≈ I
        @test size(metric2) == (3, 3)
        metric2 = Pathfinder.RankUpdateEuclideanMetric((4,))
        @test metric2.M⁻¹ ≈ I
        @test size(metric2) == (4, 4)
        metric2 = Pathfinder.RankUpdateEuclideanMetric(Float32, (4,))
        @test metric2.M⁻¹ ≈ I
        @test size(metric2) == (4, 4)
        @test eltype(metric2.M⁻¹) === Float32

        @test size(metric) == (5, 5)
        @test AdvancedHMC.renew(metric, M⁻¹2).M⁻¹ === M⁻¹2
        @test sprint(show, metric) == "RankUpdateEuclideanMetric(diag=$(diag(metric.M⁻¹)))"
        @test AdvancedHMC.neg_energy(h, r, θ) ≈ AdvancedHMC.neg_energy(h_dense, r, θ)
        @test AdvancedHMC.∂H∂r(h, r) ≈ AdvancedHMC.∂H∂r(h_dense, r)
        compare_estimates(
            identity,
            [rand(metric) for _ in 1:10_000],
            [rand(metric_dense) for _ in 1:10_000],
        )
    end

    @testset "sample" begin
        ndraws = 1_000
        nadapts = 500
        nparams = 5
        x = 0:0.01:1
        y = sin.(x) .+ randn.() .* 0.2 .+ x
        X = [x x .^ 2 x .^ 3]
        θ₀ = randn(nparams)
        ℓπ = RegressionModel(X, y)

        metric = DiagEuclideanMetric(nparams)
        hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
        ϵ = find_good_stepsize(hamiltonian, θ₀)
        integrator = Leapfrog(ϵ)
        proposal = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
        adaptor = StanHMCAdaptor(
            MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator)
        )
        samples1, stats1 = sample(
            hamiltonian,
            proposal,
            θ₀,
            ndraws,
            adaptor,
            nadapts;
            drop_warmup=true,
            progress=false,
        )

        result_pf = pathfinder(ℓπ; dim=5, optimizer=Optim.LBFGS(; m=6))

        @testset "Initial point" begin
            metric = DiagEuclideanMetric(nparams)
            hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
            ϵ = find_good_stepsize(hamiltonian, θ₀)
            integrator = Leapfrog(ϵ)
            proposal = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
            adaptor = StepSizeAdaptor(0.8, integrator)
            samples2, stats2 = sample(
                hamiltonian,
                proposal,
                result_pf.draws[:, 1],
                ndraws,
                adaptor,
                nadapts;
                drop_warmup=true,
                progress=false,
            )
            compare_estimates(identity, samples2, samples1)
        end

        @testset "Initial point and metric" begin
            metric = DiagEuclideanMetric(diag(result_pf.fit_distribution.Σ))
            hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
            ϵ = find_good_stepsize(hamiltonian, θ₀)
            integrator = Leapfrog(ϵ)
            proposal = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
            adaptor = StepSizeAdaptor(0.8, integrator)
            samples3, stats3 = sample(
                hamiltonian,
                proposal,
                result_pf.draws[:, 1],
                ndraws,
                adaptor,
                nadapts;
                drop_warmup=true,
                progress=false,
            )
            compare_estimates(identity, samples3, samples1)
        end

        @testset "Initial point and final metric" begin
            metric = Pathfinder.RankUpdateEuclideanMetric(result_pf.fit_distribution.Σ)
            hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)
            ϵ = find_good_stepsize(hamiltonian, θ₀)
            integrator = Leapfrog(ϵ)
            proposal = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
            adaptor = StepSizeAdaptor(0.8, integrator)
            samples4, stats4 = sample(
                hamiltonian,
                proposal,
                result_pf.draws[:, 1],
                ndraws,
                adaptor,
                nadapts;
                drop_warmup=true,
                progress=false,
            )
            compare_estimates(identity, samples4, samples1)
        end
    end
end


text:
using MakeCshared
using Test

@testset "MakeCshared" begin

println("# printout function check")
@test printout() == nothing

println("# calc function check")
@test calc(4) == 16
@test calc(8) == 64
@test_throws MethodError calc("c")

println("# hello function check")
@test hello("Test") == "Hello Test World !!"
@test_throws MethodError hello(3)

end


text:
for N in [Float64, Rational{Int}, Float32]
    B = BallInf(ones(N, 2), N(3))
    H = Hyperrectangle(ones(N, 2), ones(N, 2))
    E = EmptySet{N}(2)

    # intersection of two sets
    I = Intersection(B, H)

    # swap
    I2 = swap(I)
    @test I.X == I2.Y && I.Y == I2.X

    # dim
    @test dim(I) == 2

    # support vector (currently throws an error)
    @test_throws ErrorException σ(ones(N, 2), I)

    # membership
    @test ones(N, 2) ∈ I && N[5, 5] ∉ I

    # boundedness (more tests in Float64-only section)
    @test isbounded(I)
    @test isbounded(Singleton(N[1]) ∩ HalfSpace(N[1], N(1)))

    # emptiness of intersection
    @test !isempty_known(I)
    @test !isempty(I)
    @test isempty_known(I)
    @test !isempty(I)

    # =================
    # IntersectionArray
    # =================

    # relation to base type (internal helper functions)
    @test LazySets.array_constructor(Intersection) == IntersectionArray
    @test LazySets.is_array_constructor(IntersectionArray)

    # intersection of an array of sets
    IA = IntersectionArray([B, H])

    # dim
    @test dim(IA) == 2

    # support vector (currently throws an error)
    @test_throws ErrorException σ(ones(N, 2), IA)

    # boundedness
    @test isbounded(IA)
    @test isbounded(IntersectionArray([Singleton(N[1]), HalfSpace(N[1], N(1))]))
    # the following tests crash because ρ(::IntersectionArray) is not implemented yet
    @test_throws ErrorException isbounded(IntersectionArray([HalfSpace(N[1], N(1)), HalfSpace(N[1], N(-1))]))
    @test_throws ErrorException !isbounded(IntersectionArray([HalfSpace(ones(N, 2), N(1)), HalfSpace(ones(N, 2), N(-1))]))

    # isempty
    @test_throws MethodError isempty(IA)

    # membership
    @test ones(N, 2) ∈ IA && N[5, 5] ∉ IA

    # array getter
    v = Vector{LazySet{N}}()
    @test array(IntersectionArray(v)) ≡ v

    # constructor with size hint and type
    IntersectionArray(10, N)

    # constraints_list for polytopic intersection

    # ================
    # common functions
    # ================

    # absorbing element
    @test absorbing(Intersection) == absorbing(IntersectionArray) == EmptySet
    @test I ∩ E == E ∩ I == IA ∩ E == E ∩ IA == E ∩ E == E
end

# ======================
# Tests for Float64 only
# ======================
for N in [Float64]
    # constraints_list for polytopic intersection
    B = BallInf(ones(N, 2), N(3))
    H = Hyperrectangle(ones(N, 2), ones(N, 2))
    I = B ∩ H
    IA = IntersectionArray([B, H])
    clist1 = constraints_list(I)
    clist2 = constraints_list(IA)
    @test ispermutation(clist1, clist2) &&
          ispermutation(clist1, [HalfSpace(N[1, 0], N(2)),
                                 HalfSpace(N[0, 1], N(2)),
                                 HalfSpace(N[-1, 0], N(0)),
                                 HalfSpace(N[0, -1], N(0))])

    # HalfSpace vs. Ball1 intersection
    X = Ball1(zeros(2), N(1));
    d = normalize(N[1, 0])

    # flat intersection at x = 1
    H = HalfSpace(N[-1, 0], N(-1)); # x >= 1

    # default algorithm
    @test ρ(d, X ∩ H) == ρ(d, H ∩ X) == N(1)

    # intersection at x = 0
    H = HalfSpace(N[1, 0], N(0)); # x <= 0

    # default algorithm
    @test ρ(d, X ∩ H) < 1e-6 && ρ(d, H ∩ X) < 1e-6

    # specify  line search algorithm
    @test ρ(d, X ∩ H, algorithm="line_search") < 1e-6 &&
        ρ(d, H ∩ X, algorithm="line_search") < 1e-6

    # boundedness
    @test isbounded(HalfSpace(N[1], N(1)) ∩ HalfSpace(N[-1], N(1)))
    @test !isbounded(HalfSpace(ones(N, 2), N(1)) ∩ HalfSpace(-ones(N, 2), N(1)))

    # HalfSpace vs. Ball2 intersection
    B2 = Ball2(zeros(2), N(1));
    @test ρ(d, B2 ∩ H) < 1e-6 && ρ(d, H ∩ B2) < 1e-6

    # Ball1 vs. Hyperplane intersection
    H = Hyperplane(N[1, 0], N(0.5)); # x = 0.5
    @test isapprox(ρ(d, X ∩ H, algorithm="line_search"), N(0.5), atol=1e-6)
    # For the projection algorithm, if the linear map is taken lazily we can use Ball1
    @test isapprox(ρ(d, X ∩ H, algorithm="projection", lazy_linear_map=true), N(0.5), atol=1e-6)
    # But the default is to take the linear map concretely; in this case, we *may*
    # need Polyhedra (in the general case), for the concrete linear map. As a valid workaround
    # if we don't want to load Polyhedra here is to convert the given set to a polygon in V-representation
    @test isapprox(ρ(d, convert(VPolygon, X) ∩ H, algorithm="projection", lazy_linear_map=false), N(0.5), atol=1e-6)

    # =====================
    # concrete operations
    # =====================
    cap =  HPolytope([HalfSpace(N[1], N(1))]) ∩ HPolytope([HalfSpace(N[-1], N(1))])  # x <= 1 && x >= -1
    p = linear_map(reshape([N(1/2)], 1, 1), cap)
    @test N[-0.5] ∈ p && N[0.5] ∈ p && (N[1.0] ∉ p || N[1.0] ∉ p)
end


text:
export debounce_time

import Base: close
import Base: show

# TODO: Work in progress
# TODO: Untested and undocumented

debounce_time(due_time::Int) = DebounceTimeOperator(due_time)

struct DebounceTimeOperator <: InferableOperator
    due_time :: Int
end

function on_call!(::Type{L}, ::Type{L}, operator::DebounceTimeOperator, source) where L
    return proxy(L, source, DebounceTimeProxy(operator.due_time))
end

operator_right(operator::DebounceTimeOperator, ::Type{L}) where L = L

struct DebounceTimeProxy <: ActorSourceProxy
    due_time :: Int
end

actor_proxy!(::Type{L}, proxy::DebounceTimeProxy, actor::A)   where { L, A } = DebounceTimeActor{L, A}(proxy.due_time, actor)
source_proxy!(::Type{L}, proxy::DebounceTimeProxy, source::S) where { L, S } = DebounceTimeObservable{L, S}(source)

struct DebounceTimeCompletionException <: Exception end
struct DebounceTimeCancellationException <: Exception end

mutable struct DebounceTimeActor{L, A} <: Actor{L}
    is_cancelled  :: Bool
    is_completed  :: Bool
    due_time      :: Int
    actor         :: A
    last_received :: Union{Nothing, L}
    condition     :: Condition

    DebounceTimeActor{L, A}(due_time::Int, actor::A) where { L, A } = begin
        self = new(false, false, due_time, actor, nothing, Condition())

        @async begin
            try
                if self.is_completed && !self.is_cancelled
                    complete!(self.actor)
                else
                    while !self.is_cancelled && !self.is_completed
                        wait(self.condition)
                        sleep(self.due_time / MILLISECONDS_IN_SECOND)
                        if !self.is_cancelled
                            next!(self.actor, self.last_received)
                            if self.is_completed
                                complete!(self.actor)
                            end
                        end
                    end
                end
            catch err
                if err isa DebounceTimeCompletionException
                    complete!(self.actor)
                elseif !(err isa DebounceTimeCancellationException)
                    error!(self.actor, err)
                end
            end
        end

        return self
    end
end

function on_next!(actor::DebounceTimeActor{L}, data::L) where L
    actor.last_received = data
    notify(actor.condition)
end

function on_error!(actor::DebounceTimeActor, err)
    actor.is_cancelled = true
    actor.is_completed = true
    close(actor, DebounceTimeCancellationException())
    error!(actor.actor, err)
end

function on_complete!(actor::DebounceTimeActor)
    actor.is_completed = true
    close(actor)
end

Base.close(actor::DebounceTimeActor, excp = DebounceTimeCompletionException()) = notify(actor.condition, excp; error = true)

@subscribable struct DebounceTimeObservable{L, S} <: Subscribable{L}
    source :: S
end

function on_subscribe!(observable::DebounceTimeObservable, actor::DebounceTimeActor)
    return DebounceTimeSubscription(actor, subscribe!(observable.source, actor))
end

struct DebounceTimeSubscription <: Teardown
    actor
    subscription
end

as_teardown(::Type{<:DebounceTimeSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::DebounceTimeSubscription)
    if !subscription.actor.is_cancelled
        close(subscription.actor, DebounceTimeCancellationException())
    end
    subscription.actor.is_cancelled = true
    unsubscribe!(subscription.subscription)
    return nothing
end

Base.show(io::IO, ::DebounceTimeOperator)              = print(io, "DebounceTimeOperator()")
Base.show(io::IO, ::DebounceTimeProxy)                 = print(io, "DebounceTimeProxy()")
Base.show(io::IO, ::DebounceTimeActor{L})      where L = print(io, "DebounceTimeActor($L)")
Base.show(io::IO, ::DebounceTimeObservable{L}) where L = print(io, "DebounceTimeObservable($L)")
Base.show(io::IO, ::DebounceTimeSubscription)          = print(io, "DebounceTimeSubscription()")


text:
module ImplicitEquations

import IntervalArithmetic: Interval, diam, isempty
using RecipesBase
using CommonEq

include("predicates.jl")
include("intervals.jl")
include("tupper.jl")
include("asciigraphs.jl")
include("plot_recipe.jl")

export screen, I_
export asciigraph
export Lt, ≪, Le, ≦, Eq, ⩵, Ne, ≶, ≷, Ge, ≧, Gt, ≫
export Neq

#export GRAPH, OInterval
#export Region, compute, negate_op
#export TRUE, FALSE, MAYBE

end # module


text:
type SubtreeMutation <: Mutation
  stage::AbstractString
  representation::KozaTreeRepresentation
  rate::Float
  prob_terminal::Float
  min_depth::Int
  max_depth::Int
  max_tree_depth::Int

  SubtreeMutation(stage::AbstractString, rep::KozaTreeRepresentation, rate::Float, pt::Float, mind::Int, maxd::Int) =
    new(stage, rep, rate, pt, mind - 1, maxd - 1, rep.max_depth)
end

"""
Provides a definition for a subtree mutation operator.
"""
type SubtreeMutationDefinition <: MutationDefinition
  """
  The name of the stage that this operator works on.
  """
  stage::AbstractString

  """
  The mutation rate for this operator determines the probability that a given
  individual will be subject to subtree mutation.
  """
  rate::Float

  """
  The probability of selecting a termainl when generating within the subtree.
  """
  p_terminal::Float

  """
  The minimum depth of the generated sub-tree.
  """
  min_depth::Int

  """
  The maximum depth of the generated sub-tree.
  """
  max_depth::Int

  SubtreeMutationDefinition() = new("", 0.01, 0.5, 1, 5)
end

"""
TODO: Description of sub-tree mutation operator.
"""
function subtree_mutation(def::Function)
  d = SubtreeMutationDefinition()
  def(d)
  d
end

"""
Composes a subtree mutation operator from a provided definition.
"""
function compose!(def::SubtreeMutationDefinition, sp::Species)
  def.stage = def.stage == "" ? genotype(sp).label : def.stage
  r = sp.stages[def.stage].representation
  SubtreeMutation(def.stage, r, def.rate, def.p_terminal, def.min_depth, def.max_depth)
end

num_inputs(o::SubtreeMutation) = 1
num_outputs(o::SubtreeMutation) = 1

"""
Performs subtree mutation on a provided Koza tree.
"""
function operate!{T <: KozaTree}(
  o::SubtreeMutation,
  stage::IndividualStage{T}
)
  # Enforce the mutation rate.
  rand() >= o.rate && return

  # Select a random node in the tree.
  tree = get(stage)
  n1 = sample(tree)
  d = depth(n1)

  # Choose the depth to grow the tree till.
  td = min(o.max_tree_depth, d + rand(o.min_depth:o.max_depth))

  # Grow the sub-tree to a depth between its current and maximum depth
  # and insert it into the tree at the position of the selected node.
  n2 = build(KozaGrowBuilder,
    o.representation,
    n1.parent,
    d,
    td,
    o.prob_terminal)
  swap!(n1.parent, n1, n2)

  return
end


text:
using Test
using SolidStateDetectors

import SolidStateDetectors.ConstructiveSolidGeometry as CSG
import SolidStateDetectors.ConstructiveSolidGeometry: Geometry

T = Float64

example_primitive_dir = joinpath(@__DIR__, "../../examples/example_primitive_files")
@testset "Test primitive read-in" begin
    @testset "Box" begin
        box_widths = Geometry(T, joinpath(example_primitive_dir, "Box.yaml"))
        box_halfwidths = Geometry(T, joinpath(example_primitive_dir, "Box_halfwidths.yaml"))
        box_hXhYhZ = Geometry(T, joinpath(example_primitive_dir, "Box_hXhYhZ.yaml"))
        box_old = Geometry(T, joinpath(example_primitive_dir, "Box_old.yaml"))
        @test box_widths isa CSG.Box
        @test box_widths == box_halfwidths == box_hXhYhZ == box_old
    end

    @testset "Cone" begin
        cone = Geometry(T, joinpath(example_primitive_dir, "Cone.yaml"))
        cone_old = Geometry(T, joinpath(example_primitive_dir, "Cone_old.yaml"))
        @test cone isa CSG.Cone{T, <:Any, <:Tuple}
        @test cone == cone_old

        cone_tube = Geometry(T, joinpath(example_primitive_dir, "Cone_tube.yaml"))
        @test cone_tube isa CSG.VaryingTube{T}
    end

    @testset "Ellipsoid" begin
        ellipsoid_full_sphere = Geometry(T, joinpath(example_primitive_dir, "Ellipsoid_full_sphere.yaml"))
        @test ellipsoid_full_sphere isa CSG.FullSphere{T}
    end

    @testset "RegularPrism" begin
        hexagon = Geometry(T, joinpath(example_primitive_dir, "RegularPrism_hexagon.yaml"))
        @test hexagon isa CSG.HexagonalPrism{T}
    end

    @testset "Torus" begin
        torus = Geometry(T, joinpath(example_primitive_dir, "Torus.yaml"))
        @test torus isa CSG.Torus{T}
    end
end


text:
# Pure tests of functions defined in linalg.jl
using JLD, BenchmarkTools

"""
Build companion matrix from a vector.
"""
function compan(c::AbstractVector{T}) where T
    N = length(c)
    A = diagm(-1 => ones(T, N-1))
    @inbounds for i in 1:N
        A[1,i] = - c[i]/c[1]
    end
    A
end

# generate a nasty ill-conditioned matrix
N = 27
(@isdefined IC) || (const IC = compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) )))

# B-chain in the region where stabilization is necessary
B_SVD_U = load("data/linalg.jld", "B_SVD_U")
B_SVD_D = load("data/linalg.jld", "B_SVD_D")
B_SVD_Vt = load("data/linalg.jld", "B_SVD_Vt")

B_QR_U = load("data/linalg.jld", "B_QR_U")
B_QR_D = load("data/linalg.jld", "B_QR_D")
B_QR_T = load("data/linalg.jld", "B_QR_T")

@testset "Linear Algebra" begin

    @testset "Decompositions" begin

        @testset "QR aka UDT" begin
            Random.seed!(1234)
            x = rand(10,10)

            u,d,t = decompose_udt(x)

            @test isapprox(u * u', I) # unitarity of u
            @test sort!(count.(iszero,eachcol(t))) == 0:9 # triangularity of t (istriu(t) == false due to pivoting)
            @test isapprox(d, [2.0718140718296794, 1.3713719555859616, 1.1527980473310226, 1.1359672199189663,
             0.9230156804883671, 0.7796143878032683, 0.5299931859601968, 0.43829746451880564,
              0.36921251585078685, 0.13351187833543765])
            @test isapprox(u * Diagonal(d) * t, x)

            D = zero(d);
            u, t = decompose_udt!(copy(x), D)

            @test isapprox(u * u', I) # unitarity of u
            @test sort!(count.(iszero,eachcol(t))) == 0:9 # triangularity of t (istriu(t) == false due to pivoting)
            @test isapprox(D, [2.0718140718296794, 1.3713719555859616, 1.1527980473310226, 1.1359672199189663,
             0.9230156804883671, 0.7796143878032683, 0.5299931859601968, 0.43829746451880564,
              0.36921251585078685, 0.13351187833543765])
            @test isapprox(u * Diagonal(D) * t, x)
        end

        @testset "SVD aka UDV" begin
            Random.seed!(1234)
            x = rand(10,10)

            u,d,vt = decompose_udv(x)

            @test isapprox(u * u', I) # unitarity of u
            @test isapprox(vt * vt', I) # unitarity of v
            @test isapprox(u * Diagonal(d) * vt, x)
            @test isapprox(d, [4.69085124499646, 1.3700826637207126, 1.2218064335214107, 1.069747771166964,
             0.8846933503887872, 0.7163644709197314, 0.5338397904461203, 0.4033960707621825,
              0.2653967719523671, 0.10076327879773647])

            u2, d2, vt2 = decompose_udv!(x)

            @test isapprox(u, u2)
            @test isapprox(d, d2)
            @test isapprox(vt, vt2)
        end
    end

    @testset "Other" begin
        # Diagonal(I)
        @test Diagonal(I) == true
        @test Diagonal(3*I) == 3
        @test Diagonal(4.2*I) == 4.2

        # eye!
        x = rand(4,4)
        eye!(x)
        @test x == I

        # Kronecker Delta
        @test kd(1,3) == false
        @test kd(1,1) == true
        @test kd(1.2,1.3) == false
        @test kd(1.2,1.2) == true
        @test kd(1.,1) == true
        @test kd(3+im*4, 3+im*4) == true
        @test kd(3+im*4, 3+im*2) == false
    end

    @testset "Safe multiply UDT and UDV" begin
        # UDT - calculating square of ill-conditioned slice matrix
        u = B_QR_U; d = B_QR_D; t = B_QR_T;
        # naive = (u * Diagonal(d) * t)^2
        # u_naive, d_naive, t_naive = decompose_udt(naive)
        U,D,T = multiply_safely(u,d,t,u,d,t)
        @test isapprox(maximum(D), 5.8846316709257896e16)
        @test isapprox(minimum(D), 1.9123571535539083e-24)

        # TODO: UDV
    end

    @testset "UDT and UDV inversion business (well-conditioned)" begin
        let mc = mc_from_inxml("parameters/O3_generic_small_system.in.xml")
            Random.seed!(1234)
            # X = rand(10,10);
            # Y = rand(10,10);
            X = similar(mc.s.greens)
            Y = similar(mc.s.greens)
            rand!(X)
            rand!(Y)
            inv_I_plus_X = inv(I + X)
            inv_X_plus_Y = inv(X+Y)
            inv_X = inv(X)
            res = similar(X)

            # ----------- UDT / QR ---------------------
            @testset "UDT" begin
                u,d,t = decompose_udt(X);
                uy,dy,ty = decompose_udt(Y);

                # inv_udt
                @test isapprox(inv_udt(u,d,t), inv_X)
                @test isapprox(inv_udt(u,d,t), inv_X)
                @test inv_udt!(res,u,d,t) == nothing
                @test isapprox(res, inv_X)

                # inv_one_plus_udt
                @test isapprox(inv_one_plus_udt(u,d,t), inv_I_plus_X)
                @test inv_one_plus_udt!(mc, res, u,d,t) == nothing
                @test isapprox(res, inv_I_plus_X)

                @test inv_one_plus_udt_scalettar!(mc, res, u,d,t) == nothing
                @test isapprox(res, inv_I_plus_X)

                # UDT_to_mat!
                @test isapprox(UDT_to_mat(u,d,t), X)
                @test UDT_to_mat!(res, u,d,t) == nothing
                @test isapprox(res, X)

                @test isapprox(UDT_to_mat(u,d,t, invert=true), inv_X)
                @test UDT_to_mat!(res, u,d,t, invert=true) == nothing
                @test isapprox(res, inv_X)

                # inv_sum_udts
                @test UDT_to_mat!(res, inv_sum_udts(u,d,t,uy,dy,ty)...) == nothing
                @test isapprox(res, inv_X_plus_Y)
                @test inv_sum_udts!(mc, res, u,d,t,uy,dy,ty) == nothing
                @test isapprox(res, inv_X_plus_Y)

                @test UDT_to_mat!(res, inv_sum_udts_scalettar(u,d,t,uy,dy,ty)...) == nothing
                @test isapprox(res, inv_X_plus_Y)

                @test inv_sum_udts_scalettar!(mc, res, u,d,t,uy,dy,ty) == nothing
                @test isapprox(res, inv_X_plus_Y)
            end

            # ----------- UDV / SVD ---------------------
            @testset "UDV" begin
                u,d,v = decompose_udv(X);
                uy,dy,vy = decompose_udv(Y);

                # inv_udv
                @test isapprox(inv_udv(u,d,v), inv_X)
                @test isapprox(inv_udv(u,d,v), inv_X)
                @test inv_udv!(res,u,d,v) == nothing
                @test isapprox(res, inv_X)

                # inv_one_plus_udv
                @test isapprox(inv_one_plus_udv(u,d,v), inv_I_plus_X)
                @test isapprox(inv_one_plus_udv_alt(u,d,v), inv_I_plus_X)
                @test isapprox(inv_one_plus_udv_scalettar(u,d,v), inv_I_plus_X)

                @test inv_one_plus_udv_scalettar!(mc, res, u,d,v) == nothing
                @test isapprox(res, inv_I_plus_X)

                # UDV_to_mat!
                @test isapprox(UDV_to_mat(u,d,v), X)
                @test UDV_to_mat!(res, u,d,v) == nothing
                @test isapprox(res, X)

                @test isapprox(UDV_to_mat(u,d,v, invert=true), inv_X)
                @test UDV_to_mat!(res, u,d,v, invert=true) == nothing
                @test isapprox(res, inv_X)

                # inv_sum_udvs
                @test UDV_to_mat!(res, inv_sum_udvs(u,d,v,uy,dy,vy)...) == nothing
                @test isapprox(res, inv_X_plus_Y)
                @test inv_sum_udvs!(mc, res, u,d,v,uy,dy,vy) == nothing
                @test isapprox(res, inv_X_plus_Y)
            end
        end
    end
end

nothing

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc