text:
export SupervisedModel
export iterations, cost, fitness, predict, predictProb, accuracy
export stop!, state, remember!, history, trainingCurve, name
export @remember!, @history

import StatsBase.predict

function defaultCallback()
end

# ==========================================================================
# Interface for all models

cost(model::Any) = @_not_implemented
cost(model::Any, data::DataSource) = @_not_implemented
fitness(model::Any) = -1 * cost(model)
fitness(model::Any, data::DataSource) = -1 * cost(model, data)
#predict(model::Any, data::DataSource) = @_not_implemented # Statsbase takes care of this
predictProb(model::Any, data::DataSource) = @_not_implemented
accuracy(model::Any, data::DataSource) = @_not_implemented
train!(callback::Function, model::Any, data::DataSource, solver::AbstractSolver, args...; nargs...) = @_not_implemented

# ==========================================================================
# Use default callback if no other is specified

train!(model::Any, data::DataSource, solver::AbstractSolver, args...; nargs...) =
  train!(defaultCallback, model, data, solver, args...; nargs...)

# ==========================================================================
# Embedding supervised learning model for convenience

type SupervisedModel{T<:Any}
  _state::State
  _iterations::Int
  _native::T
  _history::TrainingHistory{Int}
end

function SupervisedModel{T<:Any}(model::T)
  SupervisedModel{T}(:uninitialized, zero(Int), model, TrainingHistory(Int))
end

# ==========================================================================
# Macros for simpler function defintion

macro _static_always(f, command, args...)
  esc(:(($f{T<:Any})(model::SupervisedModel{T},$(args...)) = $command))
end

macro _static_training_only(f, command, args...)
  msg = string(f, " only defined during model training")
  esc(:(($f{T<:Any})(model::SupervisedModel{T},$(args...)) = :training <= model._state <= :trystop ? $command : throw(StateError($msg))))
end

macro _static_initialized_only(f, command, args...)
  msg = string(f, " not defined for uninitialized or erroneous models")
  esc(:(($f{T<:Any})(model::SupervisedModel{T},$(args...)) = model._state < :initialized ? throw(StateError($msg)) : $command))
end

macro _delegate_initialized_only(f, args...)
  vars = map(x->x.args[1], args)
  msg = string(f, " not defined for uninitialized or erroneous models")
  esc(:(($f{T<:Any})(model::SupervisedModel{T},$(args...)) = model._state < :initialized ? throw(StateError($msg)) : $f(model._native,$(vars...))))
end

# ==========================================================================
# Define the local and delegating functions

@_static_always(name, string(typeof(model._native)))
@_static_always(state, model._state) # :untrained, :training, :trained
@_static_always(iterations, model._iterations)
@_static_training_only(stop!, model._state = :trystop)
@_static_training_only(remember!, push!(model._history, model._iterations, key, value), key::Symbol, value)
@_static_initialized_only(history, get(model._history, key), key::Symbol)
@_static_initialized_only(trainingCurve, get(model._history, :trainingCurve))
@_delegate_initialized_only(cost)
@_delegate_initialized_only(cost, data::DataSource)
@_delegate_initialized_only(fitness)
@_delegate_initialized_only(fitness, data::DataSource)
@_delegate_initialized_only(predict, data::DataSource)
@_delegate_initialized_only(predictProb, data::DataSource)
@_delegate_initialized_only(accuracy, data::DataSource)

macro remember!(model, func)
  key = string(func)
  esc(:(@remember!($model, symbol($key), $func)))
end

macro history(model, func)
  key = string(func)
  esc(:(@get($model, symbol($key))))
end

# ==========================================================================
# Train method

function train!{T<:Any, D<:LabeledDataSource}(
    userCallback::Function,
    model::SupervisedModel{T},
    data::D,
    solver::AbstractSolver,
    args...;
    max_iter::Int = 1000,
    break_every::Int = -1, # -1 = automatic, 0 = off
    nargs...)
  @assert max_iter > 0
  break_every = break_every < 0 ? safeFloor(max_iter / 10) : break_every
  enable_callback = break_every > 0

  # Define callback for the native model
  function localCallback()
    model._iterations += 1
    if model._iterations % break_every == 0
      remember!(model, :trainingCurve, cost(model)) # This should always be trivial cheap
      userCallback()
      model._state._idx == @stateidx(:trystop) ? :stop : nothing
    end
  end

  # If callback is disabled only pass the empty function
  cb = enable_callback ? localCallback : defaultCallback

  # Train the native model
  model._state = :training
  try
    itersOld = model._iterations
    itersNew = train!(cb, model._native, data, solver, args...; max_iter=max_iter, nargs...)
    model._iterations = itersOld + itersNew
    model._state = :trained
  catch ex
    print("An exception was thrown while training: ", ex)
    model._state = :error
  end
  model
end

ext: jl


text:
export AbstractPolynomial

const SymbolLike = Union{AbstractString,Char,Symbol}

"""
    AbstractPolynomial{T}

An abstract container for various polynomials.

# Properties
- `coeffs` - The coefficients of the polynomial
- `var` - The indeterminate of the polynomial
"""
abstract type AbstractPolynomial{T} end

# We want  ⟒(P{α…,T}) = P{α…}; this default
# works for most cases
⟒(P::Type{<:AbstractPolynomial}) = constructorof(P)

# convert `as` into polynomial of type P based on instance, inheriting variable
# (and for LaurentPolynomial the offset)
_convert(p::P, as) where {P <: AbstractPolynomial} = ⟒(P)(as, p.var)

"""
    Polynomials.@register(name)

Given a polynomial with `name`, creates some common convenience constructors and conversions to minimize code required for implementation of a new polynomial type.

# Example
```julia
struct MyPolynomial{T} <: AbstractPolynomial{T} end

Polynomials.@register MyPolynomial
```

# Implementations
This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`).

For constructors, it implements the shortcut for `MyPoly(...) = MyPoly{T}(...)`, singleton constructor `MyPoly(x::Number, ...)`,  conversion constructor `MyPoly{T}(n::S, ...)`, and `variable` alternative  `MyPoly(var=:x)`.
"""
macro register(name)
    poly = esc(name)
    quote
        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
        Base.convert(P::Type{<:$poly}, p::$poly{T}) where {T} = P(coeffs(p), p.var)
        Base.promote(p::P, q::Q) where {T, P <:$poly{T}, Q <: $poly{T}} = p,q
        Base.promote_rule(::Type{<:$poly{T}}, ::Type{<:$poly{S}}) where {T,S} =
            $poly{promote_type(T, S)}
        Base.promote_rule(::Type{<:$poly{T}}, ::Type{S}) where {T,S<:Number} =
            $poly{promote_type(T, S)}
        $poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
            $poly{T}(coeffs, Symbol(var))
        $poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S<:Number} =
            $poly(T.(x), Symbol(var))
        function $poly(coeffs::G, var::SymbolLike=:x) where {G}
            !Base.isiterable(G) && throw(ArgumentError("coeffs is not iterable"))
            $poly(collect(coeffs), var)
        end
        $poly{T}(n::S, var::SymbolLike = :x) where {T, S<:Number} =
            n *  one($poly{T}, Symbol(var))
        $poly(n::S, var::SymbolLike = :x)  where {S  <: Number} = n * one($poly{S}, Symbol(var))
        $poly{T}(var::SymbolLike=:x) where {T} = variable($poly{T}, Symbol(var))
        $poly(var::SymbolLike=:x) = variable($poly, Symbol(var))
    end
end

macro registerN(name, params...)
    poly = esc(name)
    αs = tuple(esc.(params)...)
    quote
        Base.convert(::Type{P}, q::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = q
        Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q
        Base.promote(p::P, q::Q) where {$(αs...),T, P <:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = p,q
        Base.promote_rule(::Type{<:$poly{$(αs...),T}}, ::Type{<:$poly{$(αs...),S}}) where {$(αs...),T,S} =
            $poly{$(αs...),promote_type(T, S)}
        Base.promote_rule(::Type{<:$poly{$(αs...),T}}, ::Type{S}) where {$(αs...),T,S<:Number} =
            $poly{$(αs...),promote_type(T,S)}

        function $poly{$(αs...),T}(x::AbstractVector{S}, var::SymbolLike = :x) where {$(αs...),T,S}
            $poly{$(αs...),T}(T.(x), Symbol(var))
        end
        $poly{$(αs...)}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {$(αs...),T} =
            $poly{$(αs...),T}(coeffs, Symbol(var))
        $poly{$(αs...),T}(n::Number, var::SymbolLike = :x) where {$(αs...),T} = n*one($poly{$(αs...),T}, Symbol(var))
        $poly{$(αs...)}(n::Number, var::SymbolLike = :x) where {$(αs...)} = n*one($poly{$(αs...)}, Symbol(var))
        $poly{$(αs...),T}(var::SymbolLike=:x) where {$(αs...), T} = variable($poly{$(αs...),T}, Symbol(var))
        $poly{$(αs...)}(var::SymbolLike=:x) where {$(αs...)} = variable($poly{$(αs...)}, Symbol(var))
    end
end

# deprecated. If desired,  replace with  @registerN  type  parameters... macro
# Macros to register POLY{α, T} and POLY{α, β, T}
macro register1(name)
    @warn "@register1 is deprecated use @registerN"
    poly = esc(name)
    quote
        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
        Base.promote(p::P, q::Q) where {α,T, P <:$poly{α,T}, Q <: $poly{α,T}} = p,q
        Base.promote_rule(::Type{<:$poly{α,T}}, ::Type{<:$poly{α,S}}) where {α,T,S} =
            $poly{α,promote_type(T, S)}
        Base.promote_rule(::Type{<:$poly{α,T}}, ::Type{S}) where {α,T,S<:Number} =
            $poly{α,promote_type(T,S)}
        function $poly{α,T}(x::AbstractVector{S}, var::SymbolLike = :x) where {α,T,S}
            $poly{α,T}(T.(x), Symbol(var))
        end
        $poly{α}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {α,T} =
            $poly{α,T}(coeffs, Symbol(var))
        $poly{α,T}(n::Number, var::SymbolLike = :x) where {α,T} = n*one($poly{α,T}, Symbol(var))
        $poly{α}(n::Number, var::SymbolLike = :x) where {α} = n*one($poly{α}, Symbol(var))
        $poly{α,T}(var::SymbolLike=:x) where {α, T} = variable($poly{α,T}, Symbol(var))
        $poly{α}(var::SymbolLike=:x) where {α} = variable($poly{α}, Symbol(var))
    end
end

# Macro to register POLY{α, β, T}
macro register2(name)
    @warn "@register2  is deprecated use @registerN"
    poly = esc(name)
    quote
        Base.convert(::Type{P}, p::P) where {P<:$poly} = p
        Base.promote(p::P, q::Q) where {α,β,T, P <:$poly{α,β,T}, Q <: $poly{α,β,T}} = p,q
        Base.promote_rule(::Type{<:$poly{α,β,T}}, ::Type{<:$poly{α,β,S}}) where {α,β,T,S} =
            $poly{α,β,promote_type(T, S)}
        Base.promote_rule(::Type{<:$poly{α,β,T}}, ::Type{S}) where {α,β,T,S<:Number} =
            $poly{α,β,promote_type(T, S)}
        $poly{α,β}(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {α,β,T} =
            $poly{α,β,T}(coeffs, Symbol(var))
        $poly{α,β,T}(x::AbstractVector{S}, var::SymbolLike = :x) where {α,β,T,S<:Number} = $poly{α,β,T}(T.(x), var)
        $poly{α,β,T}(n::Number, var::SymbolLike = :x) where {α,β,T} = n*one($poly{α,β,T}, Symbol(var))
        $poly{α,β}(n::Number, va::SymbolLiker = :x) where {α,β} = n*one($poly{α,β}, Symbol(var))
        $poly{α,β,T}(var::SymbolLike=:x) where {α,β, T} = variable($poly{α,β,T}, Symbol(var))
        $poly{α,β}(var::SymbolLike=:x) where {α,β} = variable($poly{α,β}, Symbol(var))
    end
end


text:
using ParallelStencil
using LinearAlgebra

include("part2_utils.jl")

""" Computes matrix-vector product (∇^2 - c) * x in a matrix-free way, without using shared memory. """
@parallel_indices (ix, iy) function matrix_free_matvec_prod!(T, hx, hy, c, dT2)
    if (1 < ix < size(T, 1) && 1 < iy < size(T, 2))
        dT2[ix, iy] = (( T[ix+1, iy] - 2*T[ix, iy] + T[ix-1, iy] ) / hx^2 +
                       ( T[ix, iy+1] - 2*T[ix, iy] + T[ix, iy-1] ) / hy^2) - c * T[ix, iy]
    end
    return nothing
end

""" Computes matrix-vector product (∇^2 - c) * x in a matrix-free way, using shared memory. """
@parallel_indices (ix, iy) function matrix_free_matvec_prod_shmem!(T_, hx, hy, c, dT2)
    tx = @threadIdx().x + 1
    ty = @threadIdx().y + 1
    T = @sharedMem(eltype(T_), (@blockDim().x+2, @blockDim().y+2))
    if (1 <= ix <= size(T_, 1) && 1 <= iy <= size(T_, 2))
        T[tx, ty] = T_[ix, iy]
    end
    @sync_threads()
    if (1 < ix < size(T_, 1) && 1 < iy < size(T_, 2))
        if (@threadIdx().x == 1)             T[tx-1, ty] = T_[ix-1, iy] end
        if (@threadIdx().x == @blockDim().x) T[tx+1, ty] = T_[ix+1, iy] end
        if (@threadIdx().y == 1)             T[tx, ty-1] = T_[ix, iy-1] end
        if (@threadIdx().y == @blockDim().y) T[tx, ty+1] = T_[ix, iy+1] end
        @sync_threads()
        dT2[ix, iy] = (( T[tx+1, ty] - 2*T[tx, ty] + T[tx-1, ty] ) / hx^2 +
                       ( T[tx, ty+1] - 2*T[tx, ty] + T[tx, ty-1] ) / hy^2) - c * T[tx, ty]
    end
    return nothing
end

""" Computes matrix-free matrix-vector product (∇^2 - c) * x, based on an execution policy. """
function matrix_free_matvec_prod_wrapper!(p, hx, hy, c, p_hat; execution_policy=parallel_shmem)
    nx, ny = size(p)
    if execution_policy == parallel_shmem
        threads = (32, 8)
        blocks = ((nx, ny) .- 1) .÷ (threads) .+ 1
        shmem = prod(threads.+2) * sizeof(Float64)
        @parallel blocks threads shmem=shmem matrix_free_matvec_prod_shmem!(p, hx, hy, c, p_hat)
    elseif execution_policy == parallel
        @parallel matrix_free_matvec_prod!(p, hx, hy, c, p_hat)
    elseif execution_policy == serial
        error()
    else
        error()
    end
    @synchronize()
end

""" Unpreconditioned matrix-free Conjugate Gradient method for solving (∇^2 - c) * x = b. """
function cg!(x_in, b, hx, hy, c, tol, Nmax; execution_policy=parallel_shmem, verbose=false)
    nx, ny = size(b)
    normb  = norm(b)
    tolb   = tol * normb
    r      = copy(b)
    p      = copy(r)
    p_hat  = copy(r)
    x      = @zeros(nx, ny)
    normr  = Inf
    rho    = sum(r.*r)

    for i = 1:Nmax
        # p_hat .= A*p
        matrix_free_matvec_prod_wrapper!(p, hx, hy, c, p_hat, execution_policy=execution_policy)
        alpha = rho / sum(p.*p_hat)
        x .+= alpha * p
        r .-= alpha * p_hat
        normr = norm(r)
        if verbose
            println("$(i): $(normr/normb)")
        end
        if (normr < tolb)
            if verbose
                println("CG converged in $(i) iterations :)")
            end
            break
        end
        rho_old = rho
        rho = sum(r.*r)
        beta = rho / rho_old
        p .= r + beta * p
    end

    x_in .= x

    return sqrt(sum(r.^2) / (nx * ny))
end


text:
const twopi  = Double64(6.283185307179586, 2.4492935982947064e-16)
const onepi  = Double64(3.141592653589793, 1.2246467991473532e-16)
const halfpi = Double64(1.5707963267948966, 6.123233995736766e-17)
const qrtrpi = Double64(0.7853981633974483, 3.061616997868383e-17)
const sixteenthpi = Double64(0.19634954084936207, 7.654042494670958e-18)
const thirtysecondpi = Double64(0.09817477042468103, 3.827021247335479e-18)
const threesixteenthpi = Double64(0.5890486225480862, 2.296212748401287e-17)

atanxy(x::T, y::T) where {T<:Real} = atan(y, x)

#=
     sin(a) from the Taylor series.
     Assumes |a| <= pi/32.
=#
function sin_taylor(a::Double64)
    iszero(a) && return(a)

    x = -square(a)
    r = a
    for i = 3:2:nused_inv_fact
        r = r * x
        t = r * inv_fact[i]
        a = a + t
    end

    return a
end

#=
   1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! ...
=#
function cos_taylor(a::Double64)
    iszero(a) && return(one(Double64))

    x2 = square(a)
    r = one(a)
    a = one(a)
    for i = 2:4:(nused_inv_fact-2)
        r = r * x2
        t = r * inv_fact[i]
        a = a - t
        r = r * x2
        t = r * inv_fact[i+2]
        a = a + t
    end

    return a
end

function sincos_taylor(a::Double64)
    if iszero(a)
        return a, (one(Double64))
    end
    s = sin_taylor(a)
    c = cos_taylor(a)
    return s,c
end

function tan_taylor(a::Double64)
    iszero(a) && return(zero(Double64))
    s = sin_taylor(a)
    c = cos_taylor(a)
    return s/c
end

function csc_taylor(a::Double64)
    return inv(sin_taylor(a))
end

function sec_taylor(a::Double64)
    return inv(cos_taylor(a))
end

function cot_taylor(a::Double64)
    iszero(a) && return(zero(Double64))
    s = sin_taylor(a)
    c = cos_taylor(a)
    return c/s
 end

function index_npio32(x::DoubleFloat{T}) where {T<:IEEEFloat}
    x < npio32[1] && return 1
    x >= npio32[end] && return length(npio32)
    result = 1
    while x >= npio32[result]
        result += 1
    end
    return max(1,result-1)
end

#=
   sin(a+b) = sin(a)*cos(b) + cos(a)*sin(b)
   cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
=#

@inline function sin_circle(x::DoubleFloat{T}) where {T<:IEEEFloat}
    idx = index_npio32(x)
    pipart = npio32[idx]
    rest = x - pipart
    sin_part = sin_npio32[idx]
    cos_part = cos_npio32[idx]
    sin_rest, cos_rest = sincos_taylor(rest)
    result1 = sin_part * cos_rest
    result2 = cos_part * sin_rest
    result  = result1 + result2
    return result
end

@inline function cos_circle(x::DoubleFloat{T}) where {T<:IEEEFloat}
    idx = index_npio32(x)
    pipart = npio32[idx]
    rest = x - pipart
    sin_part = sin_npio32[idx]
    cos_part = cos_npio32[idx]
    sin_rest, cos_rest = sincos_taylor(rest)
    result1 = cos_part * cos_rest
    result2 = sin_part * sin_rest
    result  = result1 - result2
    return result
end

@inline function tan_circle(x::DoubleFloat{T}) where {T<:IEEEFloat}
    idx = index_npio32(x)
    pipart = npio32[idx]
    rest = x - pipart
    sin_part = sin_npio32[idx]
    cos_part = cos_npio32[idx]
    sin_rest, cos_rest = sincos_taylor(rest)
    sin_result1 = sin_part * cos_rest
    sin_result2 = cos_part * sin_rest
    sin_result  = sin_result1 + sin_result2
    cos_result1 = cos_part * cos_rest
    cos_result2 = sin_part * sin_rest
    cos_result  = cos_result1 - cos_result2
    return sin_result / cos_result
end

function sincos_circle(x::DoubleFloat{T}) where {T<:IEEEFloat}
    idx = index_npio32(x)
    pipart = npio32[idx]
    rest = x - pipart
    sin_part = sin_npio32[idx]
    cos_part = cos_npio32[idx]
    sin_rest, cos_rest = sincos_taylor(rest)
    s1 = sin_part * cos_rest
    s2 = cos_part * sin_rest
    s  = s1 + s2
    c1 = cos_part * cos_rest
    c2 = sin_part * sin_rest
    c  = c1 - c2
    return s, c
end

@inline function sin_kernel(x::DoubleFloat{T}) where {T<:IEEEFloat}
    signbit(x) && return -sin(abs(x))
    iszero(x) && return zero(typeof(x))
    !isfinite(x) && return nan(typeof(x))
    if x >= twopi
       x = mod2pi(x)
    end
    if x >= onepi
       z = -sin_circle(x_minus_onepi(x))
    elseif x >= halfpi
       z = cos_circle(x_minus_halfpi(x))
    elseif x <= thirtysecondpi
       z = sin_taylor(x)
    else
       z = sin_circle(x)
    end
    return z
end

@inline function cos_kernel(x::DoubleFloat{T}) where {T<:IEEEFloat}
    signbit(x) && return cos(abs(x))
    iszero(x) && return one(typeof(x))
    !isfinite(x) && return nan(typeof(x))
    if x >= twopi
       x = mod2pi(x)
    end
    if x >= onepi
       z = -cos_circle(x_minus_onepi(x))
    elseif x >= halfpi
       z = -sin_circle(x_minus_halfpi(x))
    elseif x <= thirtysecondpi
       z = cos_taylor(x)
    else
       z = cos_circle(x)
    end
    return z
end

function cos(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("cos(x) only defined for finite x"))
    return abs(x.hi) < 6.28125 ? cos_kernel(x) : DoubleFloat{T}(cos(Quadmath.Float128(x)))
end

function sin(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("sin(x) only defined for finite x"))
    return abs(x.hi) < 6.28125 ? sin_kernel(x) : DoubleFloat{T}(sin(Quadmath.Float128(x)))
end

Base.sincos(x::DoubleFloat) = (sin(x), cos(x))

#=
function tangent(x::T) where {T}
    signbit(x) && return -tangent(abs(x))
    x > halfpi(T) && return tangent(DoubleFloats.modhalfpi(x))
    x > qrtrpi(T) && return inv(tangent_0_qrtrpi(halfpi(T)-x))
    return sin(x)/cos(x)
end
function tangent_0_qrtrpi(x)
   c = cos(x)                #  c = cos(x); cc = c*c; return sqrt((1-cc)/cc)
   s = sqrt(1 - c*c)
   return s/c
end

=#
function tan(x::Double64)
    isnan(x) && return x
    isinf(x) && throw(DomainError("tan(x) only defined for finite x"))
    abs(HI(x)) >= 0.36815538909255385 && return Double64(tan(Float128(x)))  # (15/128 * pi)

    abs(mod1pi(x-Double64(pi)/2)) <= eps(one(DoubleFloat{Float64})) && return DoubleFloat{Float64}(Inf)
    iszero(x) && return zero(typeof(x))
    !isfinite(x) && return nan(typeof(x))
    signbit(x) && return -tan(-x)
    HI(x) <= 2.0e-12 && return x

    y = mod1pi(x)                        # 0 <= y < pi
    if y >= halfpi
        y = x_minus_onepi(y)                # -pi/2 < y < 0
        return tan(y)
    elseif y >= qrtrpi
        y = x_minus_qrtrpi(y)       # 0 < y < pi/4
        t = tan(y)
        return (1+t)/(1-t)
    elseif y >= threesixteenthpi
        y = -x_minus_qrtrpi(y)     #   0 < y < pi/16
        t = tan(y)
        return (1-t)/(1+t)
    end
    return tan_circle(y)               # 0 <= y < 3pi/16 [(3/4 * pi/4)< 0.5891]
end

const tan0qrtrpi_numercoeffs = [
 Double64(-4.589387262410812e-34, 3.615269061456329e-50),
 Double64(-1.1277602868617984, -3.7260835757356473e-17),
 Double64(0.023504022820282806, -8.687970367035411e-19),
 Double64(0.15800109650013386, 1.0039888332689995e-17),
 Double64(-0.0032234179678582767, 7.13543487768582e-20),
 Double64(-0.00485067399529607, 2.5007870884389995e-19),
 Double64(9.183636558700604e-5, 6.505326980348073e-21),
 Double64(4.3378563109937034e-5, -1.5862612770718012e-21),
 Double64(-6.679371552585612e-7, 2.526644583588966e-24),
 Double64(-8.978969477885054e-8, 4.807483176770941e-24),
 Double64(6.62853386447739e-10, 1.0925131977072127e-26)
];

const tan0qrtrpi_denomcoeffs = [
 Double64(-1.1277602868617984, -3.726083575735698e-17),
 Double64(0.023504022820282806, -8.68797036610413e-19),
 Double64(0.5339211921207333, 5.0215742527305713e-17),
 Double64(-0.011058092241285879, 2.163933343013038e-19),
 Double64(-0.03245636645396739, -2.0950660440965625e-18),
 Double64(0.0006439974033112582, -4.907277021564753e-20),
 Double64(0.0005359286750031091, 3.4429445937395886e-20),
 Double64(-9.392512261553479e-6, -6.384836458590209e-22),
 Double64(-2.4709845162823393e-6, -1.6760465043339782e-22),
 Double64(3.015269279339435e-8, -8.41236508210075e-25),
 Double64(1.5859442637424446e-9, 6.30689666197131e-26)
];

const tan0qrtrpi_numerpoly = Poly(tan0qrtrpi_numercoeffs);
const tan0qrtrpi_denompoly = Poly(tan0qrtrpi_denomcoeffs);

function tan0qrtrpi(x::Double64)
     numer = polyval(tan0qrtrpi_numerpoly, x)
     denom = polyval(tan0qrtrpi_denompoly, x)
     return numer/denom
end

function csc(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("csc(x) only defined for finite x"))
    return inv(sin(x))
end

function sec(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("sec(x) only defined for finite x"))
    return inv(cos(x))
end

function cot(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("cot(x) only defined for finite x"))
    abs(mod1pi(x)) <= eps(one(DoubleFloat{T})) && return DoubleFloat{T}(Inf)
    return inv(tan(x))
end

function sinpi(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("sinpi(x) only defined for finite x"))
    return DoubleFloat{T}(sinpi(Quadmath.Float128(x)))
    #=
    y = Double64(x)
    hi,lo = mul322(pi_1o1_t64, HILO(y))
    y = Double64(hi, lo)
    z = sin(y)
    return DoubleFloat{T}(z)
    =#
end

function cospi(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("cospi(x) only defined for finite x"))
    return DoubleFloat{T}(cospi(Quadmath.Float128(x)))
    #=
     y = Double64(x)
    hi,lo = mul322(pi_1o1_t64, HILO(y))
    y = Double64(hi, lo)
    z = cos(y)
    return DoubleFloat{T}(z)
    =#
end

function tanpi(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && throw(DomainError("tanpi(x) only defined for finite x"))
    return sinpi(x)/cospi(x)
    #=
    y = Double64(x)
    hi,lo = mul322(pi_1o1_t64, HILO(y))
    y = Double64(hi, lo)
    z = tan(y)
    return DoubleFloat{T}(z)
    =#
end


text:
#=Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние прямоугольные перегородки
в пятой задаче перегородки были изолированы-теперь перегородки могут быть любыми

РЕЗУЛЬТАТ: Робот - в исходном положении и в углах поля стоят маркеры=#
function mark_angles_but_harder(r::Robot)
    steps=0
    num_steps=[]
    while isborder(r,Sud)==false && isborder(r,West)==false
        push!(num_steps,get_num_movements!(r,West))
        push!(num_steps,get_num_movements!(r,Sud))
        steps=steps+2 #сумма отрезков(образуют изгибы)
    end
    for side in (Nord,Ost,Sud,West)
        get_num_movements!(r,side)
        putmarker!(r)
    end
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for _ in 1:num_steps[steps]
            move!(r,side)
        end
        steps=steps-1
    end
end

function get_num_movements!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
"constructor: univariate Particles"
function Particles{T <: AbstractFloat}(x::Vector{T}, w::Vector{Float64})
    p = 1
    n = length(x)
    @assert n == length(w)
    isMultivariate = false
    isNormalized = (sum(w) == 1.0)
    @assert all(w .>= 0)
    return Particles(n, p, isMultivariate, isNormalized, x, w)
end

"constructor: multivariate Particles"
function Particles{T <: AbstractFloat}(x::Matrix{T}, w::Vector{Float64})
    p, n = size(x)
    @assert n == length(w)
    isMultivariate = true
    isNormalized = (sum(w) == 1.0)
    @assert all(w .>= 0)
    return Particles(n, p, isMultivariate, isNormalized, x, w)
end

function normalize!(p::Particles)
    p.isNormalized = true
    p.w = p.w / sum(p.w)
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
function rand(::Type{T}, mod::S, n::Int64) where {T<:Union{DNASeq, RNASeq}, S <: NASM}
  return rand(T, Weights(_π(mod)), n)
end

function simulate!(root_seq::T,
                   tree::Tree,
                   mod::S,
                   site_rates::Vector{Float64}) where {T <: Union{DNASeq, RNASeq}, S <: NASM}
  # Simulation order
  visit_order = reverse(postorder(tree))
  # Sequence length
  len = length(site_rates)
  node_data = Dict{Int64, T}()
  node_data[visit_order[1]] = root_seq
  # Error checking
  if length(root_seq) != len
    throw(ErrorException("Dimension of root sequence must match length of site rates"))
  end
  # Iterate through nodes
  for i in visit_order[2:end]
    source = tree.branches[tree.nodes[i].in[1]].source
    source_seq = node_data[source]
    branch_length = tree.branches[tree.nodes[i].in[1]].length
    wv = [Weights(P(mod, branch_length * site_rates[j]) * source_seq.data[:, j]) for j = 1:len]
    node_data[i] = rand(T, wv, checkinput=false)
  end
  return node_data
end

function simulate(::Type{T},
                  tree::Tree,
                  mod::S,
                  site_rates::Vector{Float64}) where {T <: Union{DNASeq, RNASeq}, S <: NASM}
  # Simulation order
  visit_order = reverse(postorder(tree))
  # Sequence length
  len = length(site_rates)
  node_data = Dict{Int64, T}()
  # Generate root sequence
  node_data[visit_order[1]] = rand(T, mod, len)
  # Iterate through nodes
  for i in visit_order[2:end]
    source = tree.branches[tree.nodes[i].in[1]].source
    source_seq = node_data[source]
    branch_length = tree.branches[tree.nodes[i].in[1]].length
    wv = [Weights(P(mod, branch_length * site_rates[j]) * source_seq.data[:, j]) for j = 1:len]
    node_data[i] = rand(T, wv, checkinput=false)
  end
  return node_data
end

function simulate(::Type{T},
                   tree::Tree,
                   mod::S,
                   n::Int64) where {T <: Union{DNASeq, RNASeq}, S <: NASM}
  # Simulation order
  visit_order = reverse(postorder(tree))
  node_data = Dict{Int64, T}()
  # Generate root sequence
  node_data[visit_order[1]] = rand(T, mod, n)
  # Iterate through nodes
  for i in visit_order[2:end]
    source = tree.branches[tree.nodes[i].in[1]].source
    source_seq = node_data[source]
    branch_length = tree.branches[tree.nodes[i].in[1]].length
    pmat = P(mod, branch_length)
    wv = [Weights(pmat * source_seq.data[:, j]) for j = 1:n]
    node_data[i] = rand(T, wv, checkinput=false)
  end
  return node_data
end


text:
using Documenter, SphericalFunctions

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

makedocs(
    sitename="SphericalFunctions.jl",
    modules = [SphericalFunctions],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),  # Use clean URLs, unless built as a "local" build
        edit_link = "main",  # Link out to "main" branch on github
        canonical = "https://moble.github.io/SphericalFunctions.jl/stable/",
    ),
    pages = [
        "Introduction" => "index.md",
        "Basics" => "manual.md",
        "Utilities" => "utilities.md",
    ],
    # doctest = false
)

deploydocs(
    repo="github.com/moble/SphericalFunctions.jl.git",
    devbranch="main"
)


text:
function NetworkLoad(testsystem)
   PUcorrection = 1
   if testsystem == "AP15busDN"
      usingMatPower = 0
   elseif testsystem == "IEEE33busDN" || testsystem == "IEEE33busDN2" || testsystem == "ISONE8busTN"
      usingMatPower = 0
   elseif testsystem == "R2-25.00-1"
      usingMatPower = 0
      #PUcorrection = 0.01000
   else
      usingMatPower = 1
      PUcorrection = 100
      if testsystem =="case141_OPF"
         PUcorrection = 10
      end
   end

   if usingMatPower == 0
      filename_Node = pwd()"/data/$testsystem/Node.csv"
      filename_Generator = pwd()"/data/$testsystem/Generator.csv"
      filename_Line = pwd()"/data/$testsystem/Line.csv"

      busmat = readcsv(filename_Node, header=true)[1]
      buses = Bus[]
      for i in 1:size(busmat,1)
         bindex = busmat[i,1]
         btype = 0
         Pd = busmat[i,2]
         Qd = busmat[i,3]
         Gs = busmat[i,6]
         Bs = busmat[i,7]
         area = 0
         Vm = 0
         Va = 0
         baseKV = busmat[i,8]
         bzone = 0
         Vmax = busmat[i,4]
         Vmin = busmat[i,5]
         b = Bus(bindex, btype, Pd, Qd, Gs, Bs, area, Vm, Va, baseKV, bzone, Vmax, Vmin)
         push!(buses, b)
      end

      genmat = readcsv(filename_Generator, header=true)[1]
      generators = Generator[]
      for i in 1:size(genmat,1)
            gindex = genmat[i,1]
            gtype = genmat[i,2]
            location = genmat[i,3]
            Pg = genmat[i,4]
            Qg = genmat[i,5]
            Qmax = genmat[i,6]
            Qmin = genmat[i,7]
            Vg = genmat[i,8]
            mBase = genmat[i,9]
            status = genmat[i,10]
            Pmax = genmat[i,11]
            Pmin = genmat[i,12]
            cost = [genmat[i,13], genmat[i,14], genmat[i,15]]
            SUcost = genmat[i,16]
            SDcost = genmat[i,17]
            RU = genmat[i,18]
            RD = genmat[i,19]
            UPtime = genmat[i,20]
            DNtime = genmat[i,21]
            g = Generator(gindex, gtype, location, Pg, Qg, Qmax, Qmin, Vg, mBase, status, Pmax, Pmin, cost, SUcost, SDcost, RU, RD, UPtime, DNtime)
            push!(generators, g)
      end

      branchmat = readcsv(filename_Line, header=true)[1]
      lines = Line[]
      for i in 1:size(branchmat,1)
         lindex = i
         fbus = Int(branchmat[i,2])
         tbus = Int(branchmat[i,1])
         r = branchmat[i,3]
         x = branchmat[i,4]
         b = branchmat[i,5]
         u = branchmat[i,6]
         push!(buses[tbus].children, fbus)#children
         push!(buses[fbus].ancestor, tbus)#ancestor
         push!(buses[tbus].inline, lindex)
         push!(buses[fbus].outline, lindex)
         l = Line(lindex, fbus, tbus, r, x, b, u)
         push!(lines,l)
      end
      for l in 1:length(lines)
         for k in 1:length(lines)
            if lines[l].fbus == lines[k].tbus
               push!(lines[l].ancestor, k)
            elseif lines[l].tbus == lines[k].fbus
               push!(lines[l].children, k)
            end
         end
      end
      datamat = Dict()
      datamat["bus"] = busmat
      datamat["branch"] = branchmat
      datamat["gen"] = genmat
   elseif usingMatPower == 1
      ## MatPowerCases load ### ex) testsystem = "case6ww" /testsystem = "case118"
      mpc = loadcase(testsystem)
      datamat = mpc
      #Sbase = mpc["baseMVA"] * 1e6 # MVA to VA
      #Vbase = mpc["bus"][:,10] * 1e3 # kV to V
      #Zbase = Vbase^2 / Sbase

      buses = Bus[]
      for i in 1:size(mpc["bus"],1)
         bindex = mpc["bus"][i,1]
         btype = mpc["bus"][i,2]
         Pd = mpc["bus"][i,3]
         Qd = mpc["bus"][i,4]
         Gs = mpc["bus"][i,5]
         Bs = mpc["bus"][i,6]
         area = mpc["bus"][i,7]
         Vm = mpc["bus"][i,8]
         Va = mpc["bus"][i,9]
         baseKV = mpc["bus"][i,10]
         bzone = mpc["bus"][i,11]
         Vmax = mpc["bus"][i,12]
         Vmin = mpc["bus"][i,13]
         b = Bus(bindex, btype, Pd, Qd, Gs, Bs, area, Vm, Va, baseKV, bzone, Vmax, Vmin)
         push!(buses, b)
      end

      generators = Generator[]
      for i in 1:size(mpc["gen"],1)
         gindex = i
         gtype = "NotDefined"
         location = mpc["gen"][i,1]
         Pg = mpc["gen"][i,2]
         Qg = mpc["gen"][i,3]
         Qmax = mpc["gen"][i,4]
         Qmin = mpc["gen"][i,5]
         Vg = mpc["gen"][i,6]
         mBase = mpc["gen"][i,7]
         status = mpc["gen"][i,8]
         Pmax = mpc["gen"][i,9]
         Pmin = mpc["gen"][i,10]
         cost = mpc["gencost"][i,5:end]
         SUcost = mpc["gencost"][i,2]
         SDcost = mpc["gencost"][i,3]
         RU = 0
         RD = 0
         UPtime = 0
         DNtime = 0
         g = Generator(gindex, gtype, location, Pg, Qg, Qmax, Qmin, Vg, mBase, status, Pmax, Pmin, cost, SUcost, SDcost, RU, RD, UPtime, DNtime)
         push!(generators, g)
      end

      lines = Line[]
      for i in 1:size(mpc["branch"],1)
         lindex = i
         fbus = Int(mpc["branch"][i,2]) #definition of fbus is based on  {Low, 2013} which is opposite of Matpower
         tbus = Int(mpc["branch"][i,1])
         r = mpc["branch"][i,3]/PUcorrection
         x = mpc["branch"][i,4]/PUcorrection
         b = mpc["branch"][i,5]
         u = 10000
         push!(buses[tbus].children, fbus)#children
         push!(buses[fbus].ancestor, tbus)#ancestor
         push!(buses[tbus].inline, lindex)
         push!(buses[fbus].outline, lindex)
         l = Line(lindex, fbus, tbus, r, x, b, u)
         push!(lines,l)
      end
      #case1
      if testsystem == "case141_OPF"
         for l in 1:6
            lines[l].u = 50
         end
         for l in 7:length(lines)
            lines[l].u = 5.0
         end
      end
      #case2
      if testsystem == "case141_OPF"
         for l in 15:31
            #lines[l].u = 1.4
         end
         for l in 46:51
            #lines[l].u = 1.25
         end
      end

      for l in 1:length(lines)
         for k in 1:length(lines)
            if lines[l].fbus == lines[k].tbus
               push!(lines[l].ancestor, k)
            elseif lines[l].tbus == lines[k].fbus
               push!(lines[l].children, k)
            end
         end
      end
   end
   return buses, lines, generators, datamat
end


text:
@testset "mala" begin

    # smoke test a function without retval gradient
    @gen function foo()
        x = @trace(normal(0, 1), :x)
        return x
    end

    (trace, _) = generate(foo, ())
    (new_trace, accepted) = mala(trace, select(:x), 0.1)

    # smoke test a function with retval gradient
    @gen (grad) function foo()
        x = @trace(normal(0, 1), :x)
        return x
    end

    (trace, _) = generate(foo, ())
    (new_trace, accepted) = mala(trace, select(:x), 0.1)
end


text:
module InertiaWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractForceModel = @jimport org.orekit.forces.AbstractForceModel
const Class = @jimport java.lang.Class
const Field = @jimport org.hipparchus.Field
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeDerivativesEquations = @jimport org.orekit.propagation.numerical.FieldTimeDerivativesEquations
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const InertialForces = @jimport org.orekit.forces.inertia.InertialForces
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const TimeDerivativesEquations = @jimport org.orekit.propagation.numerical.TimeDerivativesEquations
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export InertialForces
export acceleration
export add_contribution
export depends_on_position_only
export equals
export get_class
export get_events_detectors
export get_field_events_detectors
export get_parameter_driver
export get_parameters
export get_parameters_drivers
export hash_code
export init
export is_supported
export notify
export notify_all
export to_string
export wait

include("inertial_forces.jl")

end


text:
import Pkg
Pkg.instantiate()

Pkg.add("Statistics")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("StatsBase")

using DataFrames, CSV, Statistics, StatsBase
data = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

# describe(data[:ATT])

for col in eachcol(data);
    if (eltype(col[2]) <: Real) println(describe(col[2]));
    end;
end;


text:
module Sass

using libsass_jll

include("version.jl")
include("context.jl")
include("options.jl")
include("compiler.jl")
include("julian_api.jl")

examplefolder = normpath(joinpath(@__DIR__, "..", "test", "css"))

end # module


text:
validate(square) = (square < 1 || square > 64) && throw(DomainError("Invalid square"))

"""Calculate the number of grains on square `square`."""
function on_square(square)
    validate(square)
    Int128(2)^(square-1)
end

"""Calculate the total number of grains after square `square`."""
function total_after(square)
    validate(square)
    (Int128(2)^square)-1
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Flux, DiffEqFlux, DiffEqSensitivity
using DifferentialEquations
using DiffEqSensitivity
using DiffEqNoiseProcess
using Plots
using Printf, BSON
using QuantumOptics
using Statistics
using Zygote
using StaticArrays
using LinearAlgebra
using JLD  #saving models

#################################################
#read parameters from external file
include("parameters_qb.jl")
using Main.parameters

struct Parameters
	dim::Int64
	w::Float32
	force_mag::Float32
	n_steps::Int64
	dt::Float32
	n_substeps::Int64
	gamma::Float32
end

MyParameters=Parameters(
	parameters.parameters_para["N"],
	parameters.parameters_para["w"],
	parameters.parameters_para["force_mag"],
	parameters.parameters_para["max_episode_steps"],
	parameters.parameters_para["dt"],
	parameters.parameters_para["n_substeps"],
	parameters.parameters_para["gamma"])

# nr of parallel realizations
n_par = 64

# loss function hyperparameters
C1 = 0.8f0 # evolution state fid
C2 = 0.001f0 # action amplitudes
C3 = 1.8f0*MyParameters.n_steps/50 #enhance last 50 steps

########################
using Random
Random.seed!(2)

################################
# Hamiltonian and evolution
basis = SpinBasis(1//2) #max n=dim-1
σ_z = sigmaz(basis)
σ_x = sigmax(basis)
σ_p = sigmap(basis)
σ_m = sigmam(basis)
σ_pm = σ_p*σ_m

H_0pre = MyParameters.w/2.0f0*σ_z #drift Hamiltonian
H_1pre = σ_x #control field

#to the dense regime, all real, otherwise Flux.train! complains
H_0 = real(Array{Float32}(H_0pre.data))
H_1 = real(Array{Float32}(H_1pre.data))

Arrayσ_p = real(Array{Float32}(σ_p.data))
Arrayσ_m = real(Array{Float32}(σ_m.data))
Arrayσ_pm = real(Array{Float32}(σ_pm.data))

###################
# DRIFT TERM FOR SDE
function qb_dynamics_dt!(du,u, α, t)  #du/u have a standard dimension
	ψRe = u[1:MyParameters.dim]   #(dim,)
	ψIm = u[MyParameters.dim+1:end]

	HRe = H_0.+α[1]*H_1 #size (dim,dim) ,alpha is just a single number

	Reρ = ψRe.*transpose(ψRe)+ψIm.*transpose(ψIm)
	ex_x = real(sum(diag((Arrayσ_p+Arrayσ_m)*Reρ)))

	HIm = -Arrayσ_pm/2 +ex_x*Arrayσ_m
	du[1:MyParameters.dim] = dψRe = HIm*ψRe+HRe*ψIm; #size dim)
	du[MyParameters.dim+1:end] = dψIm = HIm*ψIm-HRe*ψRe;
end

###################
# DIFFUSION TERM FOR SDE
function qb_dynamics_dW!(du,u, α, t)  #last action IS NOT STORED IN u
	ψRe = u[1:MyParameters.dim]   #(dim)
	ψIm = u[MyParameters.dim+1:end]

	HRe = Arrayσ_m

	du[1:MyParameters.dim,:] = dψRe = HRe*ψRe; #size (dim)
	du[MyParameters.dim+1:end,:] = dψIm = HRe*ψIm;
end

################################################
#model
@info("Constructing the model...") #input: current state
state_1 = Dense(2*MyParameters.dim, 256,relu ,initb = Flux.glorot_uniform)
state_2 = Dense(256, 128,relu, initb = Flux.glorot_uniform)
state_3 = Dense(128, 64,relu, initb = Flux.glorot_uniform)
state_4 = Dense(64, 1,softsign,initb = Flux.glorot_uniform)

model = Chain(state_1, state_2, state_3,state_4)
p1, re = Flux.destructure(model)

###############################################
# initial state anywhere on the Bloch sphere
u0 = Array{Float32,2}(undef, 2*MyParameters.dim,n_par)
fill!(u0,0.0f0)
u0[2,:] .= 1.0f0 #down state

function prepare_initial!(u0) #random position on the Bloch sphere
	fill!(u0,0.0f0)
	theta = acos.(2*rand(n_par).-[1])  #uniform sampling for cos(theta) between -1 and 1
	phi = rand(n_par)*2*pi
	#real parts
	u0[1,:] += cos.(theta/2)
	u0[2,:] += sin.(theta/2).*cos.(phi)
	#imag parts
	#u0[3,:].+=0
	u0[4,:] += sin.(theta/2).*sin.(phi)
	#for sure normalize initial state
	norm_factor = sqrt.(sum(u0[1:2*MyParameters.dim,:].^2,dims=1))
	u0[1:2*MyParameters.dim,:] = u0[1:2*MyParameters.dim,:]./norm_factor
	return u0
end

##########################
# target state
# ψtar = |up>
ut_complex = Array{Float32,1}(undef, MyParameters.dim)
fill!(ut_complex,0.0f0)
ut_complex[1,:] .= 1.0f0
Re_ut = real(Array{Float32}(ut_complex))
Im_ut = imag(Array{Float32}(ut_complex))
#To visualize
Fock_target = SVector{MyParameters.dim}(Re_ut.^2+Im_ut.^2)

############################################
#Mutate vector outside gradients so Zygote will not complain!
# Necessary for collectiong results during the training
mut(A,index,x) = (A[index] = x)
mut_row(A,index,x) = (A[index,:]=x)
mut_vec(A,x) = (A.=x)
#tell Zygote not to take gradients of these
Zygote.@nograd mut
Zygote.@nograd mut_row
Zygote.@nograd mut_vec

##################################
# time range for the solver
t_interval = round(MyParameters.n_substeps*MyParameters.dt,digits=5)
tspan = (0.0f0,t_interval)

#########################################
#Static Arrays to collect results
mean_fid_store = zeros(MVector{MyParameters.n_steps+1,Float32})
std_fid_store = zeros(MVector{MyParameters.n_steps+1,Float32})
single_traj_fid_store = zeros(MMatrix{MyParameters.n_steps+1,n_par,Float32})
mean_action_store = zeros(MVector{MyParameters.n_steps,Float32})
std_action_store = zeros(MVector{MyParameters.n_steps,Float32})
single_traj_action_store = zeros(MMatrix{MyParameters.n_steps,n_par,Float32})
Fock_end_example = zeros(MVector{MyParameters.dim,Float32})

###########################
# compute loss
###########################

#normalize the state using callback function
condition(u,t,integrator) = true
function affect!(integrator)
	integrator.u = integrator.u/norm(integrator.u)
end
cb = DiscreteCallback(condition,affect!,save_positions = (false,false))

#creating the noise process
CreateGrid(W1) =  (NoiseGrid(Array{Float32}((0.0:MyParameters.dt:(t_interval+MyParameters.dt))),W1))
Zygote.@nograd CreateGrid #avoid taking grads of this function

# set scalar random process
W = sqrt(MyParameters.dt)*randn(Float32,MyParameters.n_substeps+1) #for 1 trajectory
W1 = cumsum([zero(MyParameters.dt); W[1:end-1]], dims=1)
NG = CreateGrid(W1)
# define SDE problem
prepare_initial!(u0)
α = zeros(Float32,1,n_par)
prob = SDEProblem{true}(qb_dynamics_dt!,qb_dynamics_dW!, u0[:,1], tspan,α[1,1], noise=NG )

function loss_along_trajectory(p1)
	#initial values
	loss = 0.0f0
	u = u0
	α = zeros(Float32,1,n_par) #start with zero actions

	function prob_func(prob, i, repeat)
		#prepare tge vector of Wiener Process
		W = sqrt(MyParameters.dt)*randn(Float32,MyParameters.n_substeps+1) #for 1 trajectory
		W1 = cumsum([zero(MyParameters.dt); W[1:end-1]], dims=1)
		NG = CreateGrid(W1) # EM and RKMil
		remake(prob,p = α[i:i],u0 = u[:,i],noise = NG)
	end

	#initial state
	Re_uj = @view u[1:MyParameters.dim,:]
	Im_uj = @view u[MyParameters.dim+1:end,:]
	fid = zeros(Float32,1,n_par)
	fid += abs2.(sum(Re_ut.*Re_uj,dims=1))
#	fid += abs2.(sum(Im_ut.*Im_uj,dims=1))  #commented because Im_ut=0, for general target uncomment
#	fid += abs2.(sum(Im_ut.*Re_uj,dims=1))
	fid += abs2.(sum(Re_ut.*Im_uj,dims=1))
#	fid += 2*(sum(Re_ut.*Re_uj,dims=1)).*sum(Im_ut.*Im_uj,dims=1)
#	fid -= 2*(sum(Re_ut.*Im_uj,dims=1)).*sum(Im_ut.*Re_uj,dims=1)
	mut_row(single_traj_fid_store,1,fid)
	mut(mean_fid_store,1,mean(fid))
	mut(std_fid_store,1,std(fid))

	for j in 1:MyParameters.n_steps
		#update action from NN
		α = re(p1)(u).*MyParameters.force_mag #dim (1,n_par)
		#define ensemble problem
		ensembleprob = EnsembleProblem(prob,prob_func = prob_func)

		u = Array(solve(ensembleprob, RKMil(),
		ensemblealg = EnsembleThreads(), # EnsembleCPUArray(), EnsembleDistributed()
		sensealg = ForwardDiffSensitivity(),
		dt = MyParameters.dt,# save_everystep=false,
		saveat = [t_interval],callback = cb, adaptive = false, trajectories = n_par, batch_size = n_par,
		save_noise = false,
		timeseries_errors = false,weak_timeseries_errors = false,weak_dense_errors = false
		))[:,1,:]

		Re_uj = @view u[1:MyParameters.dim,:]
		Im_uj = @view u[MyParameters.dim+1:end,:]
		fid = zeros(Float32,1,n_par)
		fid += abs2.(sum(Re_ut.*Re_uj,dims=1)) #average over n_par
	#	fid += abs2.(sum(Im_ut.*Im_uj,dims=1))
	#	fid +=a bs2.(sum(Im_ut.*Re_uj,dims=1))
		fid += abs2.(sum(Re_ut.*Im_uj,dims=1))
	#	fid += 2*(sum(Re_ut.*Re_uj,dims=1)).*sum(Im_ut.*Im_uj,dims=1)
	#	fid -= 2*(sum(Re_ut.*Im_uj,dims=1)).*sum(Im_ut.*Re_uj,dims=1)
		loss += C1*MyParameters.gamma^j*(1-mean(fid))

		#store data
		mut(mean_fid_store,j+1,mean(fid))
		mut(std_fid_store,j+1,std(fid))
		mut_row(single_traj_fid_store,j+1,fid)
		mut(mean_action_store,j,mean(α))
		mut(std_action_store,j,std(α))
		mut_row(single_traj_action_store,j,α)
		#punish large actions--note, that for j we are pointing to the j-1st action!
		loss += C2*MyParameters.gamma^(j)*(mean(abs2.(α))) #mimic...sum max valus
   		#emphasize the main interval
   		if j>(MyParameters.n_steps-50)
  			loss += C3*MyParameters.gamma^j*(1-mean(fid))
   			mut_vec(Fock_end_example,Re_uj[:,1].^2+Im_uj[:,1].^2);
   		end
	end
	return loss
end

###################################
# training loop parameters
epochs = 3000
println("total epochs: ",epochs)
training = zeros(epochs)
maxgrads = zeros(epochs)
some_nans = zeros(epochs)

data = Iterators.repeated((), epochs)
opt = ADAM(0.0001)

###
#grad clipping
function clip(x)
    min.(max.(x,-40), 40)
end

using DelimitedFiles  #save txt files
## Save model
using BSON: @save
function qb_train!(loss, p1, data, opt,u0)
	ps = Flux.params(p1)
	iter = 0
	for d in data
		iter += 1
		prepare_initial!(u0)  #different initial states!
		#ini_fid = abs2.(sum(Re_ut.*u0[1:2],dims=1))+abs2.(sum(Re_ut.*u0[3:4],dims=1)) #initial mean fidelity
		@show iter
		# @show mean(ini_fid)
		@time gs = gradient(ps) do
			training_loss = loss(Zygote.hook(clip,p1))
			mut(training,iter,training_loss)
			println("loss: ",training_loss)
			return training_loss
		end
		maxgrads[iter] = maximum(abs.(gs[p1]))
		some_nans[iter] = sum(isnan.(gs[p1]))
		println("is nan: ",sum(isnan.(gs[p1])))
		println("max grad: ",maximum(abs.(gs[p1])))
		if iter%1 == 0
			fig1 = plot( [1:MyParameters.n_steps+1,1:MyParameters.n_steps+1],
			  [mean_fid_store mean_fid_store],
			    fillrange=[mean_fid_store.-std_fid_store mean_fid_store.+std_fid_store], fillalpha=0.3, c=:blue,
				xlabel = "steps", ylim=(0,1),xlim=(0,MyParameters.n_steps), lw = 1.5, title="Fidelity" ,legend=false)
			fig1=plot!(1:MyParameters.n_steps+1,single_traj_fid_store[:,1], c=:black  )
			fig2 = plot( [1:MyParameters.n_steps,1:MyParameters.n_steps],
			  [mean_action_store mean_action_store],
			    fillrange = [mean_action_store.-std_action_store mean_action_store.+std_action_store], fillalpha=0.3, c=:orange,
				xlabel = "steps", ylim=(-MyParameters.force_mag,MyParameters.force_mag), xlim=(0,MyParameters.n_steps), lw = 1.5, title="Action" ,legend=false)
			fig2=plot!(1:MyParameters.n_steps,single_traj_action_store[:,1], c=:black  )
		# uncomment to display every epoch performance
		#	display(plot(fig1, fig2, layout = (1, 2), legend = false,size=(800,360)));
			if iter == 1 || iter%100 == 0
				png(string("Figures/Figure_",iter)) #for saving figs
				epoch_output = hcat(1:MyParameters.n_steps,mean_fid_store[1:end-1],single_traj_fid_store[1:end-1,1],mean_action_store,single_traj_action_store[:,1]) #save results for epoch
				writedlm(string("Data/Epoch",iter,".txt"), epoch_output)
			end
			println("iter: ", iter)
		end
		if iter%1000 == 0
			@save string("Data/model_state-in_parallel_",iter,".bson")  p1 re opt
		end
		println("++++++++++++++++++++++")
		Flux.Optimise.update!(opt, ps, gs)
  	end
	@save string("Data/model_state-in_parallel_END.bson")  p1 re opt
end

#training
Random.seed!(10)
qb_train!(loss_along_trajectory, p1, data, opt,u0)

#plotting of the loss function
plot(training, xlabel = "epochs", ylabel="loss", lw = 1.5, title = "Training")
savefig("Figures/Loss.pdf")
#save training a txt
writedlm("Data/Loss.txt", training)


text:
mutable struct CutCollection{T <: AbstractFloat}
    cuts::Vector{SparseOptimalityCut{T}}
    q::T
    considered::Int
    id::Int

    function CutCollection(::Type{T}, id::Integer = 1) where T <: AbstractFloat
        new{T}(Vector{SparseOptimalityCut{T}}(), zero(T), 0, id)
    end
end

collection_size(collection::CutCollection) = length(collection.cuts)
considered(collection::CutCollection) = collection.considered

function aggregate(collection::CutCollection)
    return aggregate(collection.cuts, collection.id)
end

function renew!(collection::CutCollection{T}, id::Integer) where T <: AbstractFloat
    empty!(collection.cuts)
    collection.q = zero(T)
    collection.considered = 0
    collection.id = id
    return nothing
end

function add_to_collection!(collection::CutCollection, cut::HyperPlane, x::AbstractVector)
    collection.considered += 1
    return nothing
end

function add_to_collection!(collection::CutCollection, cut::HyperPlane{OptimalityCut}, x::AbstractVector)
    push!(collection.cuts, cut)
    collection.q += cut(x)
    collection.considered += 1
    return nothing
end


text:
"""
    function dl_xys_mj()
---
This function is to download all 新语丝·月刊.
"""
function dl_xys_mj()
    url = "http://www.xys.org/xys/magazine/GB"
    for year in 1994:2021
        page = parsehtml(read(download("$url/$year/"), String))
        for k in eachmatch(sel"li", page.root)
            file = Gumbo.text(k)
            if file[1:3] == "xys"
                download("$url/$year/$file", "tmp.txt") # man change save path
                run(`iconv -c -f gbk -t utf8 tmp.txt -o $file`)
            end
        end
    end
    #run(`tar jcvf xys.tar.bz2 xys*.txt`)
end


text:
"""
The value target is the discounted root value of the search tree td_steps into the
future, plus the discounted sum of all rewards until then.
"""
function compute_target_value(history::GameHistory, index::Int)::Float32
    bootstrap_index = index + conf.td_steps
    if bootstrap_index < length(history.root_values)
        root_values = (isnothing(history.reanalysed_predicted_root_values) ? history.root_values : history.reanalysed_predicted_root_values )
        last_step_value = (history.to_play_history[bootstrap_index] == history.to_play_history[index] ? root_values[bootstrap_index] : -root_values[bootstrap_index] )
        value = last_step_value * conf.discount^conf.td_steps
		for (i, reward) in enumerate(history.reward_history[index:bootstrap_index])
			# The value is oriented from the perspective of the current player
			value += (history.to_play_history[index] == history.to_play_history[index + i] ? reward : -reward) * conf.discount^i
		end
    else
        value = 0.0f0
    end
	# @info "compute_target_value" index bootstrap_index conf.td_steps size(history.reward_history)
    return value
end

"""
Generate targets for every unroll steps.
"""
function make_target(history::GameHistory, state_index::Int)::Tuple{Vector{Float32},Vector{Float32},Array{Float32,2},Vector{Int}}
    target_values, target_rewards, target_policies, actions = Vector{Float32}(), Vector{Float32}(), Matrix{Float32}(undef, length(conf.action_space), 0), Vector{Int}()
	# @info "make_target" state_index conf.num_unroll_steps
    for current_index = state_index:state_index + conf.num_unroll_steps # returns vectors of length num_unroll_steps+1
        if current_index < length(history.root_values)
			value = compute_target_value(history, current_index)
            append!(target_values, value)
            append!(target_rewards, history.reward_history[current_index])
            target_policies = hcat(target_policies, history.child_visits[:, current_index])
            append!(actions, history.action_history[current_index])
        elseif current_index == length(history.root_values)
            append!(target_values, 0)
            append!(target_rewards, history.reward_history[current_index])
            # uniform policy
            target_policies = hcat(target_policies, fill(1.0f0 / length(conf.action_space), length(conf.action_space)))
            append!(actions, history.action_history[current_index])
        else
            # States past the end of games are treated as absorbing states
            append!(target_values, 0)
            append!(target_rewards, 0)
            target_policies = hcat(target_policies, fill(1.0f0 / length(conf.action_space), length(conf.action_space)))
            append!(actions, rand(rng, conf.action_space))
        end
    end
    return target_values, target_rewards, target_policies, actions
end

"""
Used in Reanalyze
"""
function update_history!(history::GameHistory, buffer::Dict{Int,GameHistory}, game_id::Int)::Nothing
    # The element could have been removed since its selection and training
	# Checks if the game_id still in buffer, where buffer keys are sorted Ints
	if first(keys(buffer)) ≤ game_id
		if conf.PER
			# Avoid read only array when loading replay buffer from disk
			history.priorities = copy(history.priorities)
		end
		buffer[game_id] = history
	end
	return nothing
end

"""
Sample position from game either uniformly or according to some priority.
See paper appendix Training.
"""
function sample_position(history::GameHistory, force_uniform=false)::Tuple{Int,Float32}
    position_prob = 0.0f0
    if conf.PER && !force_uniform
        position_probs = history.priorities ./ sum(history.priorities)
        position_index = rand(rng, Categorical(position_probs))
        position_prob = position_probs[position_index]
    else
        position_index = rand(1:length(history.root_values))
    end
	# @info "sample_position" length(history.root_values) position_index
    return position_index, position_prob
end

"""
Used in get_batch, either samples `n` games according to PER or randomly
"""
function sample_n_games(buffer::Dict{Int,GameHistory}, force_uniform=false)::Vector{Tuple{Int,GameHistory,Float32}}
    if conf.PER && !force_uniform
        game_id_list = Vector{Int}()
        game_probs = Vector{Float32}()
        for (game_id, history) in buffer
            append!(game_id_list, game_id)
            push!(game_probs, history.game_priority)
        end
        game_probs ./= sum(game_probs)
        game_prob_dict = Dict(game_id => prob for (game_id, prob) in zip(game_id_list, game_probs))
        selected_games = [game_id_list[i] for i in rand(rng, Categorical(game_probs), conf.batch_size)]
		n_games = [(game_id, buffer[game_id], game_prob_dict[game_id]) for game_id in selected_games]
    else
        selected_games = rand(collect(keys(buffer)), conf.batch_size)
        game_prob_dict = Dict()
		n_games  = [(game_id, buffer[game_id], 0.0f0) for game_id in selected_games]
    end
    return n_games
end

"""
Used in PER
"""
function sample_game(buffer::Dict{Int,GameHistory}, force_uniform=false)::Tuple{Int,GameHistory,Float32}
    game_prob = 0.0f0
    if conf.PER && !force_uniform
        game_probs = Vector{Float32}()
        for (_, history) in buffer
            append!(game_probs, history.game_priority)
        end
        game_probs ./= sum(game_probs)
        game_index = rand(rng, Categorical(game_probs))
        game_prob = game_probs[game_index]
    else
        game_index = rand(1:length(buffer))
    end
    game_id = fetch(num_played_games) - length(buffer) + game_index
    return game_id, buffer[game_id], game_prob
end

"""
Saves history in buffer
When using PER, it sets the initial priorities in history
"""
function save_game(history::GameHistory, remote_buffer::RemoteChannel{BufferChannel}, num_played_games,
num_played_steps,
total_samples,)
    if conf.PER
		# Initial priorities for the prioritized replay (See paper appendix Training)
		priorities = Vector{Float32}()
		for (i, root_value) in enumerate(history.root_values)
			priority =  abs(root_value - compute_target_value(history, i))^conf.PER_alpha
			append!(priorities, priority)
		end
		history.priorities = priorities
		history.game_priority = maximum(history.priorities)
	end

	len_hist= length(history.root_values)

	update_remote_counter!(num_played_games, 1)
	update_remote_counter!(num_played_steps, len_hist)
	update_remote_counter!(total_samples, len_hist)
	num_played_games_ = fetch(num_played_games)

	remote_buffer=put!(remote_buffer, num_played_games_, history)

	if conf.replay_buffer_size < num_played_games_
		del_id = num_played_games_ - conf.replay_buffer_size
		removed_hist= take!(remote_buffer, del_id)
		update_remote_counter!(total_samples, -(length(removed_hist.root_values)))
	end
end

"""
Used in Learning(only when PER) to Update game and position priorities with priorities calculated during the training.
See Distributed Prioritized Experience Replay (https://arxiv.org/abs/1803.00933)
"""
function update_priorities!(buffer::Dict{Int,GameHistory}, priorities::Matrix{Float32}, index_batch::Vector{Tuple})::Nothing
    for i = 1:length(index_batch)
        game_id, game_pos = index_batch[i]

        # checks if game_id still in buffer
		if first(keys(buffer)) ≤ game_id
			# Update position priorities
			priority = priorities[:,i]
			start_index = game_pos
			end_index = minimum(game_pos + length(priority), length(buffer[game_id].priorities))
			# Update game priorities
			buffer[game_id].priorities[start_index:end_index] = priority[1:end_index - start_index + 1]
			buffer[game_id].game_priority = maximum(buffer[game_id].priorities)
		end
    end
end

"""
Makes of batch of samples from the Buffer
"""
function get_batch(buffer::Dict{Int,GameHistory})::Tuple{Vector{Tuple{Int64,Float32}},Tuple{Array{Float32,4},Matrix{Float32},Matrix{Float32},Matrix{Float32},Array{Float32,3},Any,Vector{Float32}}}
    total_samples = sum([length(history.root_values) for history in values(buffer)])
    # @info "Replay buffer initialized with $(total_samples) samples"
    index_batch = Vector{Tuple{Int,Float32}}()
	observation_batch = Array{Float32}(undef, conf.observation_shape[1],conf.observation_shape[2], (conf.observation_shape[3]*(conf.stacked_observations+1)+conf.stacked_observations), 0)
	action_batch = Array{Float32}(undef, conf.num_unroll_steps+1, 0)
	reward_batch = Array{Float32}(undef, conf.num_unroll_steps+1, 0)
	value_batch = Array{Float32}(undef, conf.num_unroll_steps+1, 0)
	policy_batch = Array{Float32}(undef, length(conf.action_space), conf.num_unroll_steps+1, 0)
	gradient_scale_batch = Vector{Float32}()
    weight_batch = conf.PER ? Vector{Float32}() : nothing
	n_games= sample_n_games(buffer) # makes a batch of (batch_size) game unrolls
	# @info "Sampled N games" size(n_games)
    for (game_id, history, game_prob) in n_games
        game_pos, pos_prob = sample_position(history)
		# @info "get_batch"  size(history.observation_history) size(history.action_history) size(history.reward_history) size(history.to_play_history) size(history.child_visits) size(history.root_values)
        target_values, target_rewards, target_policies, actions = make_target(history, game_pos) # unrolls each game sample for num_unroll_steps
        push!(index_batch, (game_id, game_pos))
		# @info "Empty batch and Target sizes are:" size(observation_batch) size(action_batch)  size(reward_batch) size(value_batch ) size(policy_batch) size(gradient_scale_batch) size(target_values) size(target_rewards) size(target_policies) size(actions)
        observation_batch = cat(observation_batch, get_stacked_observations(history, game_pos, conf.stacked_observations), dims=ndims(observation_batch))
        action_batch = hcat(action_batch, actions)
        reward_batch = hcat(reward_batch, target_rewards)
        value_batch = hcat(value_batch, target_values)
        policy_batch = cat(policy_batch, target_policies, dims=ndims(target_policies) + 1)
        push!(gradient_scale_batch, minimum([conf.num_unroll_steps, length(history.action_history)+1 - game_pos]))
        conf.PER ? push!(weight_batch, 1 / (total_samples * game_prob * pos_prob)) : nothing
    end
    conf.PER ? weight_batch ./= maximum(weight_batch) : nothing
    return index_batch, (observation_batch, action_batch, value_batch, reward_batch, policy_batch, weight_batch, gradient_scale_batch)
end


text:
using Agents
using GLMakie
using InteractiveDynamics

model, flocking_agent_step!, flocking_model_step! = Models.flocking()
Bird = valtype(model.agents)

const bird_polygon = Polygon(Point2f0[(-0.5, -0.5), (1, 0), (-0.5, 0.5)])
function bird_marker(b::Bird)
    φ = atan(b.vel[2], b.vel[1]) #+ π/2 + π
    scale(rotate2D(bird_polygon, φ), 2)
end

abm_play(
    model, flocking_agent_step!, flocking_model_step!;
    am = bird_marker,
    title = "Flocking"
)


text:
# Tables interface
Tables.istable(::Type{ResultSet}) = true
Tables.rowaccess(::Type{ResultSet}) = true
Tables.rows(rs::ResultSet) = rs.rows
Tables.columnnames(rs::ResultSet) = [ Symbol(colname) for colname in names(rs) ]
Tables.columnnames(row::ResultSetRow) = [ Symbol(colname) for colname in names(row) ]
Tables.getcolumn(row::ResultSetRow, index::Symbol) = row[string(index)]
Tables.getcolumn(row::ResultSetRow, index::Integer) = row[index]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## async event notifications

"""
    AsyncCondition()

Create a async condition that wakes up tasks waiting for it
(by calling [`wait`](@ref) on the object)
when notified from C by a call to `uv_async_send`.
Waiting tasks are woken with an error when the object is closed (by [`close`](@ref).
Use [`isopen`](@ref) to check whether it is still active.
"""
mutable struct AsyncCondition
    handle::Ptr{Cvoid}
    cond::ThreadSynchronizer
    isopen::Bool

    function AsyncCondition()
        this = new(Libc.malloc(_sizeof_uv_async), ThreadSynchronizer(), true)
        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)
        err = ccall(:uv_async_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            eventloop(), this, uv_jl_asynccb::Ptr{Cvoid})
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_async_init", err))
        end
        return this
    end
end

"""
    AsyncCondition(callback::Function)

Create a async condition that calls the given `callback` function. The `callback` is passed one argument,
the async condition object itself.
"""
function AsyncCondition(cb::Function)
    async = AsyncCondition()
    waiter = Task(function()
        lock(async.cond)
        try
            while isopen(async)
                success = try
                    stream_wait(async, async.cond)
                    true
                catch exc # ignore possible exception on close()
                    isa(exc, EOFError) || rethrow()
                finally
                    unlock(async.cond)
                end
                success && cb(async)
                lock(async.cond)
            end
        finally
            unlock(async.cond)
        end
    end)
    # must start the task right away so that it can wait for the AsyncCondition before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return async
end

## timer-based notifications

"""
    Timer(delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object).

Waiting tasks are woken after an initial delay of `delay` seconds, and then repeating with the given
`interval` in seconds. If `interval` is equal to `0`, the timer is only triggered once. When
the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.
"""
mutable struct Timer
    handle::Ptr{Cvoid}
    cond::ThreadSynchronizer
    isopen::Bool

    function Timer(timeout::Real; interval::Real = 0.0)
        timeout ≥ 0 || throw(ArgumentError("timer cannot have negative timeout of $timeout seconds"))
        interval ≥ 0 || throw(ArgumentError("timer cannot have negative repeat interval of $interval seconds"))

        this = new(Libc.malloc(_sizeof_uv_timer), ThreadSynchronizer(), true)
        err = ccall(:uv_timer_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), eventloop(), this)
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_timer_init", err))
        end

        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)

        ccall(:jl_uv_update_time, Cvoid, (Ptr{Cvoid},), eventloop())
        ccall(:jl_uv_timer_start,  Cint,  (Ptr{Cvoid}, Ptr{Cvoid}, UInt64, UInt64),
              this, uv_jl_timercb::Ptr{Cvoid},
              UInt64(round(timeout * 1000)) + 1, UInt64(round(interval * 1000)))
        return this
    end
end

unsafe_convert(::Type{Ptr{Cvoid}}, t::Timer) = t.handle
unsafe_convert(::Type{Ptr{Cvoid}}, async::AsyncCondition) = async.handle

function wait(t::Union{Timer, AsyncCondition})
    lock(t.cond)
    try
        isopen(t) || throw(EOFError())
        stream_wait(t, t.cond)
    finally
        unlock(t.cond)
    end
end

isopen(t::Union{Timer, AsyncCondition}) = t.isopen

function close(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL && isopen(t)
        t.isopen = false
        isa(t, Timer) && ccall(:jl_uv_timer_stop, Cint, (Ptr{Cvoid},), t)
        ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
    end
    nothing
end

function uvfinalize(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL
        disassociate_julia_struct(t.handle) # not going to call the usual close hooks
        close(t)
        t.handle = C_NULL
    end
    t.isopen = false
    nothing
end

function _uv_hook_close(t::Union{Timer, AsyncCondition})
    lock(t.cond)
    try
        uvfinalize(t)
        notify_error(t.cond, EOFError())
    finally
        unlock(t.cond)
    end
    nothing
end

function uv_asynccb(handle::Ptr{Cvoid})
    async = @handle_as handle AsyncCondition
    lock(async.cond)
    try
        notify(async.cond)
    finally
        unlock(async.cond)
    end
    nothing
end

function uv_timercb(handle::Ptr{Cvoid})
    t = @handle_as handle Timer
    lock(t.cond)
    try
        if ccall(:uv_timer_get_repeat, UInt64, (Ptr{Cvoid},), t) == 0
            # timer is stopped now
            close(t)
        end
        notify(t.cond)
    finally
        unlock(t.cond)
    end
    nothing
end

"""
    sleep(seconds)

Block the current task for a specified number of seconds. The minimum sleep time is 1
millisecond or input of `0.001`.
"""
function sleep(sec::Real)
    sec ≥ 0 || throw(ArgumentError("cannot sleep for $sec seconds"))
    wait(Timer(sec))
    nothing
end

# timer with repeated callback
"""
    Timer(callback::Function, delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object) and
calls the function `callback`.

Waiting tasks are woken and the function `callback` is called after an initial delay of `delay` seconds,
and then repeating with the given `interval` in seconds. If `interval` is equal to `0`, the timer
is only triggered once. The function `callback` is called with a single argument, the timer itself.
When the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.

# Examples

Here the first number is printed after a delay of two seconds, then the following numbers are printed quickly.

```julia-repl
julia> begin
           i = 0
           cb(timer) = (global i += 1; println(i))
           t = Timer(cb, 2, interval = 0.2)
           wait(t)
           sleep(0.5)
           close(t)
       end
1
2
3
```
"""
function Timer(cb::Function, timeout::Real; interval::Real = 0.0)
    t = Timer(timeout, interval = interval)
    waiter = Task(function()
        while isopen(t)
            success = try
                wait(t)
                true
            catch exc # ignore possible exception on close()
                isa(exc, EOFError) || rethrow()
                false
            end
            success && cb(t)
        end
    end)
    # must start the task right away so that it can wait for the Timer before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return t
end

"""
    timedwait(testcb::Function, secs::Float64; pollint::Float64=0.1)

Waits until `testcb` returns `true` or for `secs` seconds, whichever is earlier.
`testcb` is polled every `pollint` seconds.
"""
function timedwait(testcb::Function, secs::Float64; pollint::Float64=0.1)
    pollint > 0 || throw(ArgumentError("cannot set pollint to $pollint seconds"))
    start = time()
    done = Channel(1)
    timercb(aw) = begin
        try
            if testcb()
                put!(done, :ok)
            elseif (time() - start) > secs
                put!(done, :timed_out)
            end
        catch e
            put!(done, :error)
        finally
            isready(done) && close(aw)
        end
    end

    if !testcb()
        t = Timer(timercb, pollint, interval = pollint)
        ret = fetch(done)
        close(t)
    else
        ret = :ok
    end
    ret
end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
## subsystem initialization

export version, Shutdown,
       ismultithreaded

function version()
    return libllvm_version
end

ismultithreaded() = convert(Core.Bool, API.LLVMIsMultithreaded())

for subsystem in [:Core, :TransformUtils, :ScalarOpts, :ObjCARCOpts, :Vectorization,
                  :InstCombine, :IPO, :Instrumentation, :Analysis, :IPA, :CodeGen, :Target]
    jl_fname = Symbol(:Initialize, subsystem)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname(R::PassRegistry) = API.$api_fname(ref(R))
    end
end

## target initialization

for component in [:TargetInfo, :Target, :TargetMC, :AsmPrinter, :AsmParser, :Disassembler]
    jl_fname = Symbol(:Initialize, :All, component, :s)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname() = API.$api_fname()
    end
end

for component in [:Target, :AsmPrinter, :AsmParser, :Disassembler]
    jl_fname = Symbol(:Initialize, :Native, component)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname() = convert(Core.Bool, API.$api_fname()) &&
                      throw(LLVMException($"Could not initialize native $component"))
    end
end

for target in libllvm_targets,
    component in [:Target, :AsmPrinter, :AsmParser, :Disassembler, :TargetInfo, :TargetMC]
    jl_fname = Symbol(:Initialize, target, component)
    api_fname = Symbol(:LLVM, jl_fname)

    @eval begin
        export $jl_fname
        $jl_fname() = API.@apicall($(QuoteNode(api_fname)), Cvoid, ())
    end
end


text: discrete_system(ω::AbstractFloat, source::AbstractSource, material::Material; kws...) = discrete_system(ω::AbstractFloat, source::AbstractSource, material::Material, setupsymmetry(source,material); kws...)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using SparseArrays, Arpack

mutable struct Modal <: AbstractSolver
    time :: Float64
    geometric_stiffness :: Bool
    eigvals :: Vector
    eigvecs :: Matrix
    nev :: Int
    which :: Symbol
    bc_invertible :: Bool
    P :: Vector{SparseMatrixCSC}
    symmetric :: Bool
    empty_assemblies_before_solution :: Bool
    dense :: Bool
    info_matrices :: Bool
    sigma :: Float64
end

function Modal(nev=10, which=:SM)
    solver = Modal(0.0, false, [], Matrix{Float64}(undef,0,0), nev, which,
                   false, [], true, true, false, false, 0.0)
end

"""
    calc_projection(problem, ndim)

A helper function to calculate P = D^-1*M
"""
function calc_projection(problem::T) where
        {T<:Union{Problem{Mortar}, Problem{Mortar2D}}}

    C1 = sparse(problem.assembly.C1)
    C2 = sparse(problem.assembly.C2)
    @assert C1 == C2
    @assert problem.properties.adjust == false

    s = get_nonzero_rows(C2)
    m = setdiff(get_nonzero_columns(C2), s)

    # Construct matrix P = D^-1*M
    D = C2[s,s]
    M = -C2[s,m]

    if !isdiag(D)
        @warn("Mortar matrix D is not diagonal. This might take a long time.")
        P = ldlt(1/2*(D + D')) \ M
    else
        P = D \ M
    end

    return s, m, P
end

function FEMBase.eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
    isempty(problem.assembly.C2) && return nothing
    C1 = sparse(problem.assembly.C1)
    C2 = sparse(problem.assembly.C2)
    C1 == C2 || error("Cannot eliminate boundary condition $P: C1 != C2.")
    isdiag(C1) || error("Cannot eliminate boundary condition $P: C is not diagonal")
    @info("Eliminating boundary condition $(problem.name) from global system.")
    fixed_dofs = get_nonzero_rows(C1)
    K[fixed_dofs,:] .= 0.0
    K[:,fixed_dofs] .= 0.0
    M[fixed_dofs,:] .= 0.0
    M[:,fixed_dofs] .= 0.0
    dropzeros!(K)
    dropzeros!(M)
    return nothing
end

"""
    eliminate_boundary_conditions!(problem, K, M, f)

Eliminate Mortar boundary condition from matrices K, M and force vector f.
"""
function FEMBase.eliminate_boundary_conditions!(problem::T, K, M, f) where
                            {T <: Union{Problem{Mortar}, Problem{Mortar2D}}}
    @info("Eliminating mesh tie constraint $(problem.name) using static condensation")
    s, m, P = calc_projection(problem)
    ndim = size(K, 1)
    Id = ones(ndim)
    Id[s] .= 0.0
    Q = sparse(Diagonal(Id))
    Q[s,m] += P
    K[:,:] .= Q'*K*Q
    M[:,:] .= Q'*M*Q
    return nothing
end

function FEMBase.run!(solver::Solver{Modal})
    time = solver.properties.time
    problems = get_problems(solver)
    properties = solver.properties

    @info("Starting natural frequency solver at time $time")

    @timeit "assemble matrices" begin
        assemble!(solver, time; with_mass_matrix=true)
        M, K, Kg, f = get_field_assembly(solver)
        if properties.geometric_stiffness
            K += Kg
        end
    end

    dim = size(K, 1)
    ndofs = size(K, 1)

    for P in properties.P
        @info("Using P to make transformation K_red = P'*K*P and M_red = P'*M*P")
        K[:,:] .= P'*K*P
        M[:,:] .= P'*M*P
    end

    for problem in get_problems(solver)
        eliminate_boundary_conditions!(problem, K, M, f)
    end

    # free up some memory before solution
    if properties.empty_assemblies_before_solution
        for problem in get_field_problems(solver)
            empty!(problem.assembly)
        end
    end

    SparseArrays.droptol!(K, 1.0e-9)
    SparseArrays.droptol!(M, 1.0e-9)
    nz = get_nonzero_rows(K)
    K = K[nz,nz]
    M = M[nz,nz]

    sigma = 0.0
    if properties.sigma != 0.0
        @info("Adding diagonal term $(properties.sigma) to stiffness matrix")
        sigma = properties.sigma
    end

    props = solver.properties

    @debug("Calculating $(props.nev) eigenvalues...")

    if properties.symmetric
        K = Symmetric(K)
        M = Symmetric(M)
    end

    if properties.info_matrices
        @info("is K symmetric? ", issymmetric(K))
        @info("is M symmetric? ", issymmetric(M))
        @info("is K positive definite? ", isposdef(K))
        @info("is M positive definite? ", isposdef(M))
    end

    if properties.dense
        K = Matrix(K)
        M = Matrix(M)
    end

    om2 = nothing
    X = nothing
    passed = false

    try
        @timeit "solve eigenvalue problem using `eigs`" begin
            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
        end
        passed = true
    catch
        @info("Failed to calculate eigenvalues for problem.",
              issymmetric(K), issymmetric(M), isposdef(K), isposdef(M))
        if !isapprox(properties.sigma, 0.0)
            @info("Stiffness matrix is not positive definite and Cholesky " *
                  "factorization is failing. Model is not supported enough " *
                  "with boundary conditions. To work around this problem, " *
                  "use `problem.properties.sigma = <some small value>` " *
                  "To add artificial stiffness to model. (Or add boundary " *
                  "conditions.)")
            rethrow()
        end
    end

    if !passed
        sigma = props.sigma = 1.0e-9
        @info("Calculation of eigenvalues failed. Stiffness matrix is not " *
              "positive definite and Cholesky factorization is failing. Trying " *
              "again by adjusting problem.properties.sigma to $sigma.")
        try
            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
            passed = true
        catch
            @info("Failed to calculate eigenvalues with sigma value $sigma. " *
                  "Manually set sigma to something larger and try again.")
            rethrow()
        end
    end

    @info("Squared eigenvalues: $om2.")

    props.eigvals = om2
    neigvals = length(om2)
    props.eigvecs = zeros(ndofs, neigvals)
    for i=1:neigvals
        props.eigvecs[nz,i] = X[:,i]
        for problem in get_boundary_problems(solver)
            isa(problem, Problem{Mortar}) || continue
            s, m, P = calc_projection(problem)
            # FIXME: store projection to boundary problem, i.e.
            # update!(problem, "master-slave projection", time => P)
            # us = P*um
            props.eigvecs[s,i] = P*props.eigvecs[m,i]
        end
    end

    @timeit "save results to Xdmf" update_xdmf!(solver)

    return nothing

end

function update_xdmf!(solver::Solver{Modal})

    results_writers = get_results_writers(solver)
    if length(results_writers) == 0
        @info("Xdmf is not attached to solver, not writing output to a file.")
        @info("To write results to Xdmf file, attach Xdmf to Solver, i.e.")
        @info("add_results_writer!(solver, Xdmf(\"results\"))")
        return
    end
    if maximum(abs.(imag(solver.properties.eigvals))) > 1.0e-9
        @info("Writing imaginary eigenvalues for Xdmf not supported.")
        return
    end

    xdmf = first(results_writers)

    @timeit "fetch geometry" X_ = solver("geometry", solver.properties.time)
    node_ids = keys(X_)
    @timeit "create node permutation" P = Dict(j=>i for (i, j) in enumerate(node_ids))
    nnodes = length(X_)
    ndofs = round(Int, size(solver.properties.eigvecs, 1)/nnodes)
    ndim = length(X_[first(node_ids)])
    @info("Number of nodes: $nnodes. ",
          "Number of dofs/node: $ndofs. ",
          "Dimension of geometry: $ndim.")
    @timeit "create ncoords array" begin
        X = zeros(ndim, nnodes)
        for j in node_ids
            X[:, P[j]] = X_[j]
        end
    end
    geom_type = (ndim == 2 ? "XY" : "XYZ")
    all_elements = get_all_elements(solver)
    element_types = unique(map(get_element_type, all_elements))
    nelements = length(all_elements)

    elcon_arrays  = Dict()
    @timeit "create topology arrays" for element_type in element_types
        elements = collect(filter_by_element_type(element_type, all_elements))
        nelements = length(elements)
        eldim = length(element_type)
        element_conn = zeros(Int, eldim, nelements)
        for (i, element) in enumerate(elements)
            for (j, conn) in enumerate(get_connectivity(element))
                element_conn[j,i] = P[conn]-1
            end
        end
        elcon_arrays[element_type] = element_conn
    end

    xdmf_element_mapping = Dict(
        "Poi1" => "Polyvertex",
        "Seg2" => "Polyline",
        "Tri3" => "Triangle",
        "Quad4" => "Quadrilateral",
        "Tet4" => "Tetrahedron",
        "Pyramid5" => "Pyramid",
        "Wedge6" => "Wedge",
        "Hex8" => "Hexahedron",
        "Seg3" => "Edge_3",
        "Tri6" => "Tri_6",
        "Quad8" => "Quad_8",
        "Tet10" => "Tet_10",
        "Pyramid13" => "Pyramid_13",
        "Wedge15" => "Wedge_15",
        "Hex20" => "Hex_20")

    # save modes

    temporal_collection = get_temporal_collection(xdmf)
    unknown_field_name = ucfirst(get_unknown_field_name(solver))
    frames = []

    @timeit "save modes" for (j, eigval) in enumerate(real(solver.properties.eigvals))
        if eigval < 0.0
            @warn("negative real eigenvalue found, om2=$eigval, setting to zero.")
            eigval = 0.0
        end
        freq = sqrt(eigval)/(2.0*pi)
        path = "/Results/Natural Frequency Analysis/$unknown_field_name/Mode $j"
        @info("Creating frequency frame f=$(round(freq; digits=3)), path=$path")

        frame = new_element("Grid")
        time = new_child(frame, "Time")
        set_attribute(time, "Value", freq)

        geometry = new_element("Geometry")
        set_attribute(geometry, "Type", geom_type)
        data_node_ids = new_dataitem(xdmf, "/Node IDs", collect(node_ids))
        data_geometry = new_dataitem(xdmf, "/Geometry", X)
        add_child(geometry, data_geometry)
        add_child(frame, geometry)

        for element_type in element_types
            timeit("save topology of element type $element_type") do
                elements = collect(filter_by_element_type(element_type, all_elements))
                nelements = length(elements)
                element_ids = map(get_element_id, elements)
                element_conn = elcon_arrays[element_type]
                #@timeit "creaet element_conn" element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements)
                #@timeit "hcat elcon" element_conn = hcat(element_conn...)
                element_code = split(string(element_type), ".")[end]
                dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
                dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
                topology = new_element("Topology")
                set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
                set_attribute(topology, "NumberOfElements", nelements)
                add_child(topology, dataitem)
                add_child(frame, topology)
            end
        end

        @timeit "store eigenmode" begin
            mode_ = reshape(solver.properties.eigvecs[:,j], ndofs, nnodes)
            @timeit "create mode array" begin
                mode = zeros(ndofs, nnodes)
                for nid in node_ids
                    mode[:,P[nid]] = mode_[:,nid]
                end
            end
            if ndofs == 1
                field_type = "Scalar"
            elseif ndofs == 2 # extend to 3d
                field_type = "Vector"
                mode = vcat(mode, zeros(1, nnodes))
            elseif ndofs == 3
                field_type = "Vector"
            elseif ndofs == 6 # has rotation dofs, drop them
                field_type = "Vector"
                mode = mode[1:3, :]
            else
                error("Number of dofs / node = $ndofs, I don't know how to store results to Xdmf!")
            end
            field_center = "Node"
            attribute = new_child(frame, "Attribute")
            set_attribute(attribute, "Name", unknown_field_name)
            set_attribute(attribute, "Center", field_center)
            set_attribute(attribute, "AttributeType", field_type)
            dataitem = new_dataitem(xdmf, path, mode)
            add_child(attribute, dataitem)
            add_child(frame, attribute)
            add_child(temporal_collection, frame)
        end

    end

    save!(xdmf)
end

function solve!(solver::Solver{Modal}, time::Float64)
    @info("solve!(analysis, time) is deprecated. Use run!(analysis) instead.")
    solver.properties.time = time
    run!(solver)
end


text:
"""
```julia
towel(u0 = [0.085, -0.121, 0.075])
```
```math
\\begin{aligned}
x_{n+1} &= a x_n (1-x_n) -0.05 (y_n +0.35) (1-2z_n) \\\\
y_{n+1} &= 0.1 \\left( \\left( y_n +0.35 \\right)\\left( 1+2z_n\\right) -1 \\right)
\\left( 1 -1.9 x_n \\right) \\\\
z_{n+1} &= 3.78 z_n (1-z_n) + b y_n
\\end{aligned}
```
The folded-towel map is a hyperchaotic mapping due to Rössler [1]. It is famous
for being a mapping that has the smallest possible dimensions necessary for hyperchaos,
having two positive and one negative Lyapunov exponent.
The name comes from the fact that when plotted looks like a folded towel, in every
projection.

Default values are the ones used in the original paper.

[1] : O. E. Rössler, Phys. Lett. **71A**, pp 155 (1979)
"""
function towel(u0=[0.085, -0.121, 0.075])
    return DDS(eom_towel, u0, nothing, jacob_towel)
end# should result in lyapunovs: [0.432207,0.378834,-3.74638]
function eom_towel(x, p, n)
    @inbounds x1, x2, x3 = x[1], x[2], x[3]
    SVector( 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
    0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
    3.78*x3*(1-x3)+0.2*x2 )
end
function jacob_towel(x, p, n)
    @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
    -0.19((x[2] + 0.35)*(1-2x[3]) - 1)  0.1*(1-2x[3])*(1-1.9x[1])  -0.2*(x[2] + 0.35)*(1-1.9x[1]);
    0.0  0.2  3.78(1-2x[3]) ]
end

function eom_towel_iip(dx, x, p, n)
    @inbounds begin
        x1, x2, x3 = x[1], x[2], x[3]
        dx[1] = 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3)
        dx[2] = 0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1)
        dx[3] = 3.78*x3*(1-x3)+0.2*x2
    end
end
function jacob_towel_iip(J, x, p, n)
    @inbounds begin
        J[1,1] = 3.8*(1 - 2x[1])
        J[2,1] = -0.19((x[2] + 0.35)*(1-2x[3]) - 1)
        J[3,1] = 0.0
        J[1,2] = -0.05*(1-2x[3])
        J[2,2] =  0.1*(1-2x[3])*(1-1.9x[1])
        J[3,2] = 0.2
        J[1,3] = 0.1*(x[2] + 0.35)
        J[2,3] = -0.2*(x[2] + 0.35)*(1-1.9x[1])
        J[3,3] = 3.78(1-2x[3])
    end
end

"""
```julia
standardmap(u0=0.001rand(2); k = 0.971635)
```
```math
\\begin{aligned}
\\theta_{n+1} &= \\theta_n + p_{n+1} \\\\
p_{n+1} &= p_n + k\\sin(\\theta_n)
\\end{aligned}
```
The standard map (also known as Chirikov standard map) is a two dimensional,
area-preserving chaotic mapping due to Chirikov [1]. It is one of the most studied
chaotic systems and by far the most studied Hamiltonian (area-preserving) mapping.

The map corresponds to the  Poincaré's surface of section of the kicked rotor system.
Changing the non-linearity parameter `k` transitions the system from completely periodic
motion, to quasi-periodic, to local chaos (mixed phase-space) and finally to global
chaos.

The default parameter `k` is the critical parameter where the golden-ratio torus is
destroyed, as was calculated by Greene [2]. The e.o.m. considers the angle variable
`θ` to be the first, and the angular momentum `p` to be the second, while
both variables
are always taken modulo 2π (the mapping is on the [0,2π)² torus).

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : B. V. Chirikov, Preprint N. **267**, Institute of
Nuclear Physics, Novosibirsk (1969)

[2] : J. M. Greene, J. Math. Phys. **20**, pp 1183 (1979)
"""
function standardmap(u0=0.001rand(2); k = 0.971635)
    return DDS(standardmap_eom, u0, [k], standardmap_jacob)
end
@inbounds function standardmap_eom(x, par, n)
    theta = x[1]; p = x[2]
    p += par[1]*sin(theta)
    theta += p
    while theta >= twopi; theta -= twopi; end
    while theta < 0; theta += twopi; end
    while p >= twopi; p -= twopi; end
    while p < 0; p += twopi; end
    return SVector(theta, p)
end
@inbounds standardmap_jacob(x, p, n) =
@SMatrix [1 + p[1]*cos(x[1])    1;
          p[1]*cos(x[1])        1]

"""
```julia
coupledstandardmaps(M::Int, u0 = 0.001rand(2M); ks = ones(M), Γ = 1.0)
```
```math
\\begin{aligned}
\\theta_{i}' &= \\theta_i + p_{i}' \\\\
p_{i}' &= p_i + k_i\\sin(\\theta_i) - \\Gamma \\left[
\\sin(\\theta_{i+1} - \\theta_{i}) + \\sin(\\theta_{i-1} - \\theta_{i})
\\right]
\\end{aligned}
```
A discrete system of `M` nonlinearly coupled standard maps, first
introduced in [1] to study diffusion and chaos thresholds.
The *total* dimension of the system
is `2M`. The maps are coupled through `Γ`
and the `i`-th map has a nonlinear parameter `ks[i]`.

The first `M` entries of the state are the angles, the last `M` are the momenta.

[1] : H. Kantz & P. Grassberger, J. Phys. A **21**, pp 127–133 (1988)
"""
function coupledstandardmaps end
using SparseArrays
function coupledstandardmaps(M::Int, u0 = 0.001rand(2M);
    ks = ones(M), Γ = 1.0)

    SV = SVector{M, Int}
    idxs = SV(1:M...) # indexes of thetas
    idxsm1 = SV(circshift(idxs, +1)...)  #indexes of thetas - 1
    idxsp1 = SV(circshift(idxs, -1)...)  #indexes of thetas + 1

    csm = CoupledStandardMaps{M}(idxs, idxsm1, idxsp1)
    J = zeros(eltype(u0), 2M, 2M)
    # Set ∂/∂p entries (they are eye(M,M))
    # And they dont change they are constants
    for i in idxs
        J[i, i+M] = 1
        J[i+M, i+M] = 1
    end
    sparseJ = sparse(J)
    p = (ks, Γ)
    csm(sparseJ, u0, p, 0)
    return DDS(csm, u0, p, csm, sparseJ)
end
struct CoupledStandardMaps{N}
    idxs::SVector{N, Int}
    idxsm1::SVector{N, Int}
    idxsp1::SVector{N, Int}
end
function (f::CoupledStandardMaps{N})(xnew::AbstractVector, x, p, n) where {N}
    ks, Γ = p
    @inbounds for i in f.idxs

        xnew[i+N] = mod2pi(
            x[i+N] + ks[i]*sin(x[i]) -
            Γ*(sin(x[f.idxsp1[i]] - x[i]) + sin(x[f.idxsm1[i]] - x[i]))
        )

        xnew[i] = mod2pi(x[i] + xnew[i+N])
    end
    return nothing
end
function (f::CoupledStandardMaps{M})(
    J::AbstractMatrix, x, p, n) where {M}

    ks, Γ = p
    # x[i] ≡ θᵢ
    # x[[idxsp1[i]]] ≡ θᵢ+₁
    # x[[idxsm1[i]]] ≡ θᵢ-₁
    @inbounds for i in f.idxs
        cosθ = cos(x[i])
        cosθp= cos(x[f.idxsp1[i]] - x[i])
        cosθm= cos(x[f.idxsm1[i]] - x[i])
        J[i+M, i] = ks[i]*cosθ + Γ*(cosθp + cosθm)
        J[i+M, f.idxsm1[i]] = - Γ*cosθm
        J[i+M, f.idxsp1[i]] = - Γ*cosθp
        J[i, i] = 1 + J[i+M, i]
        J[i, f.idxsm1[i]] = J[i+M, f.idxsm1[i]]
        J[i, f.idxsp1[i]] = J[i+M, f.idxsp1[i]]
    end
    return J
end

"""
```julia
henon(u0=zeros(2); a = 1.4, b = 0.3)
```
```math
\\begin{aligned}
x_{n+1} &= 1 - ax^2_n+y_n \\\\
y_{n+1} & = bx_n
\\end{aligned}
```
The Hénon map is a two-dimensional mapping due to Hénon [1] that can display a strange
attractor (at the default parameters). In addition, it also displays many other aspects
of chaos, like period doubling or intermittency, for other parameters.

According to the author, it is a system displaying all the properties of the
Lorentz system (1963) while being as simple as possible.
Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : M. Hénon, Commun.Math. Phys. **50**, pp 69 (1976)
"""
function henon(u0=zeros(2); a = 1.4, b = 0.3)
    return DDS(hoop, u0, [a,b], hoop_jac)
end # should give lyapunov exponents [0.4189, -1.6229]
hoop(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
hoop_jac(x, p, n) = @SMatrix [-2*p[1]*x[1] 1.0; p[2] 0.0]

function henon_iip(u0=zeros(2); a = 1.4, b = 0.3)
    return DDS(hiip, u0, [a, b], hiip_jac)
end
function hiip(dx, x, p, n)
    dx[1] = 1.0 - p[1]*x[1]^2 + x[2]
    dx[2] = p[2]*x[1]
    return
end
function hiip_jac(J, x, p, n)
    J[1,1] = -2*p[1]*x[1]
    J[1,2] = 1.0
    J[2,1] = p[2]
    J[2,2] = 0.0
    return
end

"""
```julia
logistic(x0 = rand(); r = 4.0)
```
```math
x_{n+1} = rx_n(1-x_n)
```
The logistic map is an one dimensional unimodal mapping due to May [1] and is used by
many as the archetypal example of how chaos can arise from very simple equations.

Originally intentend to be a discretized model of polulation dynamics, it is now famous
for its bifurcation diagram, an immensly complex graph that that was shown
be universal by Feigenbaum [2].

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : R. M. May, Nature **261**, pp 459 (1976)

[2] : M. J. Feigenbaum, J. Stat. Phys. **19**, pp 25 (1978)
"""
function logistic(x0=rand(); r = 4.0)
    return DDS(logistic_eom, x0, [r], logistic_jacob)
end
logistic_eom(x, p, n) = p[1]*x*(1-x)
logistic_jacob(x, p, n) = p[1]*(1-2x)

"""
    pomaeu_manneville(u0 = 0.2; z = 2.5)
The Pomeau-Manneville map is a one dimensional discrete map which is
characteristic for displaying intermittency [1]. Specifically, for
z > 2 the average time between chaotic bursts diverges, while
for z > 2.5, the map iterates are long range correlated [2].

Notice that here we are providing the "symmetric" version:
```math
x_{n+1} = \\begin{cases}
-4x_n + 3, & \\quad x_n \\in (0.5, 1] \\\\
x_n(1 + |2x_n|^{z-1}), & \\quad |x_n| \\le 0.5 \\\\
-4x_n - 3, & \\quad x_n \\in [-1, 0.5)
\\end{cases}
```

[1] : Manneville & Pomeau, Comm. Math. Phys. **74** (1980)

[2] : Meyer et al., New. J. Phys **20** (2019)
"""
function pomeau_manneville(u0 = 0.2, z = 2.5)
    return DDS(pm_eom, u0, [z], pm_jac)
end
function pm_eom(x, p, n)
    if x < -0.5
        -4x - 3
    elseif -0.5 ≤ x ≤ 0.5
        @inbounds x*(1 + abs(2x)^(p[1]-1))
    else
        -4x + 3
    end
end
function pm_jac(x, p, n)
    if x < -0.5
        -4.0
    elseif -0.5 ≤ x ≤ 0.5
        @inbounds z = p[1]
        0.5(x^2 * 2^z * (z-1)*abs(x)^(z-3) + 2^z * abs(x)^(z-1) + 2)
    else
        -4.0
    end
end