text:
using IntegratedOptics
######################
# Plot Options
ϵ_axes = Axes(title = "'ϵᵣ'",xlabel = "'x'", ylabel = "'y'", palette = :tempo, auto = "fix",size =" ratio -1")
Ez_axes = Axes(title = "'Ez'",xlabel = "'x'", ylabel = "'y'", palette = :balance, auto = "fix",size =" ratio -1")
######################
######################## MATERIALS ########################
Si =   Material(ϵᵣ = 3.64^2 )
SiO2 = Material(ϵᵣ = 1.45^2 )
Fg = Material(ϵᵣ = .5*(1.45^2 + 3.46^2) )
##########################################################
######################## SIMULATION #######################
λ₀ = 1.55
Lx , Ly =  λ₀ .* (4 , 4)
n_max = 3.64
n_min = 1.45
dx , dy = (λ₀/10/n_max) .* (1, 1)
sim = Simulation(λ₀ = 1.55 ;  grid = Grid(extent = (Lx,Ly) , spacing =  (dx,dy) ))
sim(setpml!,1.0)

##########################################################
###################### BUILD DEVICE ######################
#--------------------------------------------------------#
sim(setbackground!,SiO2)
#--------------------------------------------------------#
w_d = 1.5
h_d = 1.5
design_region = Cuboid(((.5(Lx - w_d) , .5(Ly - h_d)),(.5(Lx + w_d) , .5(Ly + h_d))))
sim(setmaterial!,Fg, design_region)
#--------------------------------------------------------#
h_wg_in = 0.5
input_waveguide  = Cuboid(((0 , .5(Ly - h_wg_in)),(.4(Lx) , .5(Ly + h_wg_in))))
sim(setmaterial!,Si,input_waveguide)
#--------------------------------------------------------#
h_wg_out_1 = 0.5
output_waveguide_1  = Cuboid(((.6Lx , .5(Ly - h_wg_out_1)),(Lx , .5(Ly + h_wg_out_1))))
sim(setmaterial!,Si,output_waveguide_1)
#--------------------------------------------------------#
#--------------------------------------------------------#
w_wg_out_2 = 0.5
output_waveguide_2  = Cuboid(((.5(Lx - w_wg_out_2),0),(.5(Lx + w_wg_out_2),.4Ly)))
sim(setmaterial!,Si,output_waveguide_2)
#--------------------------------------------------------#
##########################################################
sim(HeatMap,:ϵᵣ,ẑ, axes = ϵ_axes)
#sim(saveplot , output = "./device.png")
################### PREPARE TARGET FIELD #################
### INPUT
########################################################
x_cut = 1.2
sim_wg_in = sim(slice, x̂ , x_cut )
eig_vals , eig_vecs = sim_wg_in(solve_for_modes,solver = EigArpack())
mode_te_in = sim_wg_in(extractreshape, eig_vecs[:,1] ,x̂)
#mode_tm_in = sim_wg(Photon.extractreshape, eig_vecs[:,2] ,ŷ)
sim_wg_in(LinePlot,mode_te_in,lw = 2,title = "'mode profile'")
input_mode_sim = deepcopy(sim)
input_mode_sim(setbackground!,SiO2)
input_mode_sim(setmaterial!,Si,Cuboid(((0 , .5(Ly - h_wg_in)),(Lx , .5(Ly + h_wg_in)))))
input_mode_sim(HeatMap,:ϵᵣ,ẑ, axes = ϵ_axes)
input_mode_sim(attachmode!,mode_te_in, input_mode_sim(:J,ẑ), x̂ , x_cut )
input_mode_sim(HeatMap,:J,ẑ)
input_mode_sim(solve!)
input_mode_sim(HeatMap,:E,ẑ, axes = Ez_axes)

### OUTPUT1
########################################################
x_cut = 1.2
sim_wg_out_1 = sim(slice, x̂ , x_cut )
eig_vals , eig_vecs = sim_wg_out_1(solve_for_modes,solver = EigArpack())
mode_te_out_1 = sim_wg_out_1(extractreshape, eig_vecs[:,1] ,x̂)
#mode_tm_in = sim_wg(Photon.extractreshape, eig_vecs[:,2] ,ŷ)
sim_wg_out_1(LinePlot,mode_te_out_1,lw = 2,title = "'mode profile'")
output_1_mode_sim = deepcopy(sim)
output_1_mode_sim(setbackground!,SiO2)
output_1_mode_sim(setmaterial!,Si,Cuboid(((0 , .5(Ly - h_wg_out_1)),(Lx , .5(Ly + h_wg_out_1)))))
output_1_mode_sim(HeatMap,:ϵᵣ,ẑ, axes = ϵ_axes)
output_1_mode_sim(attachmode!,mode_te_out_1, output_1_mode_sim(:J,ẑ), x̂ , x_cut )
output_1_mode_sim(HeatMap,:J,ẑ)
output_1_mode_sim(solve!)
output_1_mode_sim(HeatMap,:E,ẑ, axes = Ez_axes)

## OUTPUT 2
##########################################################
y_cut = 1.9
sim_wg_out_2 = sim(slice, ŷ , y_cut )
eig_vals , eig_vecs = sim_wg_out_2(solve_for_modes)
mode_te_out_2 = sim_wg_out_2(extractreshape, eig_vecs[:,3] ,x̂)
sim_wg_out_2(LinePlot,mode_te_out_2,lw = 2,title = "'mode profile'")
output_2_mode_sim = deepcopy(sim)
output_2_mode_sim(setbackground!,SiO2)
output_2_mode_sim(setmaterial!,Si, Cuboid(((.5(Lx - w_wg_out_2),0),(.5(Lx + w_wg_out_2),Ly))))
output_2_mode_sim(HeatMap,:ϵᵣ,ẑ,axes =ϵ_axes)
output_2_mode_sim(attachmode!,mode_te_out_2, output_2_mode_sim(:J,ẑ), ŷ , Ly-y_cut )
output_2_mode_sim(HeatMap,:J,ẑ)
output_2_mode_sim(solve!)
output_2_mode_sim(HeatMap,:E,ẑ,axes = Ez_axes)
#########################################################
import LinearAlgebra: normalize
target_x = zeros(ComplexF64,size(sim.grid))
target_y = zeros(ComplexF64,size(sim.grid))
target_z = zeros(ComplexF64,size(sim.grid))

in_p = 1 / sqrt(abs(real.(input_mode_sim(:S,x̂,((x̂,2,2),(ŷ,3,3))))[1]))
input_mode_sim.E .*= in_p

out_1_p = 1 / sqrt(abs(real.(output_1_mode_sim(:S,x̂,((x̂,4.5,4.5),(ŷ,3,3))))[1]))
output_1_mode_sim.E .*= out_1_p

out_2_p = 1 / sqrt(abs(real.(output_2_mode_sim(:S,ŷ,((x̂,3,3),(ŷ,2,2))))[1]))
output_2_mode_sim.E .*= out_2_p

sim(add!,target_z , input_mode_sim(:E,ẑ),  x-> .4Lx > x[1] > x_cut)
sim(add!,target_z , (1/sqrt(4)) *  output_2_mode_sim(:E,ẑ), x-> 0 <= x[2] <= 2.3)
sim(add!,target_z , (sqrt(3)/sqrt(4)) *  output_1_mode_sim(:E,ẑ), x-> Lx-2.3 <= x[1] <= Lx)

sim(HeatMap,target_z, axes = Ez_axes)

###########################################################
sim(add!,sim(:J,ẑ) ,input_mode_sim(:E,ẑ), EVERYWHERE)
#sim(attachmode!,mode_te_in, sim(:J,ẑ), x̂ , x_cut )
#sim(HeatMap,  :J, ẑ)
sim(setTFSF!,x̂,LOW,x_cut + dx)
#sim.ϵᵣ_zz = input_mode_sim.ϵᵣ_zz
sim(solve!)
#############################################################
# Toplogy Optimization
#######################################
# Optimization Settings
obj =TargetObjective(target_x,target_y,target_z)
optimizer = Descent(1.0)
using ImageFiltering
using LocalFilters
kernel = Kernel.gaussian(2)
ϵ_limits = (1.45^2,3.64^2)
Nt = 100
binarization_rate(i) =  0.5* i * -(-(ϵ_limits...))/Nt
###########################################
# Optimization Loop
for i in 1:Nt
    f , ∇f  =  f_∇f(obj,sim,solver = LU())
    sim(update!,optimizer, sim(:ϵᵣ,ẑ), localmean(∇f[ẑ], 3),design_region)
    sim(update!,optimizer, sim(:ϵᵣ,ẑ), closing(∇f[ẑ], 5),design_region)

# if mod(i,10)==0   sim(set!, sim(:ϵᵣ,ẑ),localmean( sim(:ϵᵣ,ẑ),2),design_region) end
    sim(boxconstraint!,sim(:ϵᵣ,ẑ),ϵ_limits,design_region)
    sim(gradualbinarize!,sim(:ϵᵣ,ẑ),ϵ_limits,binarization_rate(i),design_region)
   ####
  #  sim(saveplot,term = "./plots/eps/eps_$(i).png")
    trans_2 = real.(sim(:S,ŷ,((x̂,3,3),(ŷ,2,2))))
    trans_1 = real.(sim(:S,x̂,((x̂,4.5,4.5),(ŷ,3,3))))

    println("iter:$(i)",": ",trans_1, " , ",trans_2 )
   # display(sim(HeatMap,  :ϵᵣ, ẑ,axes =ϵ_axes ))
  #  sim(saveplot,output = "./plots/devices/ninety_bend/eps/eps_$(i).png")
 #   display(sim(HeatMap,  :E, ẑ,axes =Ez_axes ))
#    sim(saveplot,output = "./plots/devices/ninety_bend/Ez/Ez_$(i).png")
       ####
end

max_stars_repo_path: examples/integrated optics/three ports/power dividers/TO_power_divider.jl


text:
# Test Sherman-Morrison inverse
A = rand(2, 2)
Ainv = inv(A)
u = rand(2, 1)
B = A + u * u'
@test NonlinearBandits.sherman_morrison_inv(Ainv, u, u) ≈ inv(B)
@test_throws(
    DimensionMismatch("u and v have mismatched dimensions"),
    NonlinearBandits.sherman_morrison_inv(Ainv, u, zeros(3, 1))
)
@test_throws(
    ArgumentError("u and v should have exactly one column"),
    NonlinearBandits.sherman_morrison_inv(Ainv, zeros(2, 2), zeros(2, 2))
)
@test_throws(
    DimensionMismatch("size of A does not match u/v"),
    NonlinearBandits.sherman_morrison_inv(Ainv, zeros(3, 1), zeros(3, 1))
)

## Test Bayesian linear model fitting procedure
d, n, σ = 2, 10000, 2
β = rand(Normal(0, 10), d, 1)
f(x) = β' * x
X, y = NonlinearBandits.gaussian_data(d, n, f; σ=σ)
shape0, scale0 = 0.001, 0.001

lm = BayesLM(d; shape0=shape0, scale0=scale0)
@test_throws(
    ArgumentError("normal-inverse-gamma has undefined mean for shape ≤ 1"), std(lm)
)
fit!(lm, X, y)

prediction_mae = mae(lm(X), mapslices(f, X; dims=1))
β_mae = mae(lm.β, β)
@test β_mae < 0.5
@test prediction_mae < 0.5
@test lm.shape == shape0 + n / 2
@test isapprox(std(lm), σ; atol=0.5)

# Check that sequential updates lead to the same results
sequential_lm = BayesLM(d; shape0=shape0, scale0=scale0)
for i in 1:n
    fit!(sequential_lm, X[:, i:i], y[:, i:i])
end

@test sequential_lm.β ≈ lm.β
@test sequential_lm.Σ ≈ lm.Σ
@test sequential_lm.shape == shape0 + (n / 2)
@test sequential_lm.scale ≈ lm.scale

# Check errors
@test_throws(
    ArgumentError("d, shape0 and scale0 must be strictly positive"),
    BayesLM(0; shape0=0.0, scale0=0.0)
)
@test_throws(ArgumentError("y should have exactly 1 row"), fit!(lm, X, y'))
@test_throws(
    DimensionMismatch("X and y don't match in dimension 2"), fit!(lm, X[:, 1:2], y)
)
@test_throws(
    DimensionMismatch("X does not have the same number of predictors as the model"),
    fit!(lm, X[1:1, :], y)
)
@test_throws(
    DimensionMismatch("X does not have the same number of predictors as the model"),
    lm(X[1:1, :])
)

ext: jl


text:
mutable struct AeroParameters
    Cl::Float64
    Cd::Float64
    A::Float64
    m::Float64
    function AeroParameters(α=deg2rad(12))

        # lift (Gallais, "Atmospheric Re-Entry Vehicle Mechanics," p.82-83)
        Cl1 = 1.42
        Cl =Cl1*α

        # drag (Gallais, "Atmospheric Re-Entry Vehicle Mechanics," p.82-83)
        Cd0 = 1.65
        Cd2 = -2.77
        Cd =Cd0 + Cd2*α^2

        # cross sectional area and mass
        A = π*(2.25)^2 # m²
        m = 2400.0     # kg

        return new(Cl, Cd, A, m)
    end
end

@inline function LD_mags(aero::AeroParameters,ρ::T,r::SVector{3,T},v::SVector{3,T}) where T

    dynamic_ρ = 0.5*ρ*aero.A*dot(v,v)/aero.m

    L = aero.Cl*dynamic_ρ
    D = aero.Cd*dynamic_ρ

    return L,D
end

@inline function e_frame(r::SVector{3, T},v::SVector{3, T}) where T
    e1 = cross(r,v)
    e1 = e1/norm(e1)
    e2 = cross(v,e1)
    e2 = e2/norm(e2)
    return e1,e2
end

# let
#
#
#     r = SA[1,4,.3]
#     v = SA[.3,.4,.2]
#
#
#     @btime e1,e2 = e_frame($r,$v)
#
# end


text:
function GeoMagneticElements(arg0::Vector3D)
    return GeoMagneticElements((Vector3D,), arg0)
end

function get_declination(obj::GeoMagneticElements)
    return jcall(obj, "getDeclination", jdouble, ())
end

function get_field_vector(obj::GeoMagneticElements)
    return jcall(obj, "getFieldVector", Vector3D, ())
end

function get_horizontal_intensity(obj::GeoMagneticElements)
    return jcall(obj, "getHorizontalIntensity", jdouble, ())
end

function get_inclination(obj::GeoMagneticElements)
    return jcall(obj, "getInclination", jdouble, ())
end

function get_total_intensity(obj::GeoMagneticElements)
    return jcall(obj, "getTotalIntensity", jdouble, ())
end

function to_string(obj::GeoMagneticElements)
    return jcall(obj, "toString", JString, ())
end


text:
using Documenter, JToric

DocMeta.setdocmeta!(JToric, :DocTestSetup, :(using JToric, Random; using JToric.Oscar); recursive = true)

doctest(JToric)


text:
content = readlines("input.txt")

drawingNums = split(popfirst!(content), ',')

# 2:6, 8:12, 14:18, 20:24

# Have to use a regex instead of just ' ' because there are multiple space characters
# Matrix(hcat(x...)') to convert to a matrix (transpose to get the correct orientation)

Base.adjoint(x::SubString) = x

boards = [Matrix(hcat(split.(content[x], r"\s+", keepempty=false)...)') for x in [range(i, stop=i+4) for i in 2:6:length(content) if i+4 <= length(content)]]

marker = "x"

markedString = repeat(marker, 5)

function checkBoard(board::Matrix)
    for row in eachrow(board)
        if join(row) == markedString
            return true
        end
    end

    for col in eachcol(board)
        if join(col) == markedString
            return true
        end
    end

    return false
end

boardScore(board, num) = sum([parse(Int, x) for x in board if x != marker]) * parse(Int, num)

function draw()

    for num in drawingNums
        for board in boards, i in 1:length(board)
            if board[i] == num
                board[i] = marker
                if checkBoard(board)
                    return boardScore(board, num)
                end
            end
        end
    end

end

draw()

#! Part 2

# Reset board
boards = [Matrix(hcat(split.(content[x], r"\s+", keepempty=false)...)') for x in [range(i, stop=i+4) for i in 2:6:length(content) if i+4 <= length(content)]]

function drawpt2()
    local trackedDrawingNum
    local winningBoardIndices::Vector{Int} = []
    for num in drawingNums
        for (boardIndex, board) in enumerate(boards), i in 1:length(board)
            if board[i] == num
                board[i] = marker
                if checkBoard(board) && !(boardIndex in winningBoardIndices)
                    push!(winningBoardIndices, boardIndex)
                    trackedDrawingNum = num
                    if length(winningBoardIndices) == length(boards)
                        return boardScore(boards[last(winningBoardIndices)], trackedDrawingNum)
                    end
                end
            end
        end
    end
end

drawpt2()

max_stars_repo_path: Day 4/main.jl

max_issues_repo_path: Day 4/main.jl

max_forks_repo_path: Day 4/main.jl


text:
abstract AbstractStateSpaceFilter

abstract AbstractKalmanFilter <: AbstractStateSpaceFilter
abstract LinearKalmanFilter <: AbstractKalmanFilter
abstract NonlinearKalmanFilter <: AbstractKalmanFilter
abstract NonlinearFilter <: AbstractStateSpaceFilter

# LinearGaussianSSM 	LinearKalmanFilter
#
# NonlinearGaussianSSM	NonlinearKalmanFilter
# 						NonlinearFilter
#
# NonlinearSSM 			NonlinearFilter

######################################################################
# Kalman filter
######################################################################

type KalmanFilter <: LinearKalmanFilter end
typealias KF KalmanFilter

## methods

# basic Kalman update, once we have the predicted state, an observation,
# and the observation matrix/covariance.
function update_kalman(pred, y, G, W)
    innovation = y - G * mean(pred)
    innovation_cov = G * cov(pred) * G' + W
    K = cov(pred) * G' * inv(innovation_cov)
    mean_update = mean(pred) + K * innovation
    cov_update = (eye(cov(pred)) - K * G) * cov(pred)
    return MvNormal(mean_update, cov_update)
end

function update(m::LinearGaussianSSM, pred::AbstractMvNormal, y::Vector,
		filter::KalmanFilter=KalmanFilter(), t::Real=0.0)
	return update_kalman(pred, y, m.G(t), m.W(t))
end

function update(m::LinearGaussianSSM, pred::AbstractMvNormal, y::Vector, t::Real=0.0)
    return update(m, pred, y, KalmanFilter(), t)
end

function update!(m::LinearGaussianSSM, fs::FilteredState, y::Vector;
		u::Vector=zeros(m.nu), filter::KalmanFilter=KalmanFilter(), t::Real=0.0)
	x_pred = predict(m, fs.state[end], u=u, t=t)
	x_filt = update_kalman(x_pred, y, m.G(t), m.W(t))
	push!(fs.state, x_filt)
	fs.observations = [fs.observations y]
end

######################################################################
# Extended Kalman filter
######################################################################
type ExtendedKalmanFilter <: NonlinearKalmanFilter end
typealias EKF ExtendedKalmanFilter

## methods

function update(m::NonlinearGaussianSSM, pred::AbstractMvNormal, y::Vector,
		filter::NonlinearKalmanFilter=EKF(), t::Real=0.0)
    G = observation_matrix(m, pred, t)
    W = m.W(t)
	return update_kalman(pred, y, G, W)
end

function update!(m::NonlinearGaussianSSM, fs::FilteredState, y::Vector;
		u::Vector=zeros(m.nu), filter::NonlinearKalmanFilter=EKF(), t::Real=0.0)
	x_pred = predict(m, fs.state[end], u=u, t=t-1)
    G = observation_matrix(m, x_pred, t)
    W = m.W(t)
	x_filt = update_kalman(x_pred, y, G, W)
	push!(fs.state, x_filt)
	fs.observations = [fs.observations y]
end

######################################################################
# Unscented Kalman filter
######################################################################

type UnscentedKalmanFilter{T<:Real} <: NonlinearKalmanFilter
	α::T
	β::T
	κ::T
end
typealias UKF UnscentedKalmanFilter
UnscentedKalmanFilter{T}(α::T=1e-3, β::T=2.0, κ::T=0.0) = UKF(α, β, κ)

"""
Data structure containing sigma points and their respective weights for the Unscented Kalman Filter

#### Fields
- `χ` : The matrix containing the 2L+1 sigma vectors, where L is the length of the system state vector.
- `wm` : Vector containing the weights for the corresponding sigma vectors required to reconsruct the predicted mean.
- `wc` : Vector containing the weights for the corresponding sigma vectors required to reconsruct the predicted covariance.
"""
type SigmaPoints{T<:Real}
    χ::Matrix{T}
    wm::Vector{T}
    wc::Vector{T}
end

"""
# calcSigmaPoints
Function to calculate the sigma points for the Unscented Transform

`calcSigmaPoints(state, α, β, κ)`
#### Parameters
- `state` : AbstractMvNormal type representing the current state estimate (mean vector with covariance matrix).
- `α` : Determines the spread of the sigma points around the mean of the system state. It is usually set to a small positive value. The default is set to 1e-3.
- `β` : Is used to incorporate prior knowledge of the distribution of the state. For Gaussian distributions a value of 2.0 is optimal. The default is set to 2.0.
- `κ` : Is a secondary scaling parameter which also determines the spread of sigma points around the mean of the system. This parameter allows for additional 'fine tuning'. It is usually set to a small value. The default is set to 0.0

`calcSigmaPoints(state, filter)`
#### Parameters
- `state` : AbstractMvNormal type representing the current state estimate (mean vector with covariance matrix).
- `filter` : UKFParameters type containing the α, β and κ parameters for the Unscented Kalman Filter.
"""
function calcSigmaPoints{T}(state::AbstractMvNormal, α::T, β::T, κ::T)
    x = mean(state)
    p = cov(state)
    L = length(x)
    χ = zeros(L, 2L+1)
    wm = zeros(2L+1)
    wc = zeros(2L+1)
    λ = α^2 * (L+κ) - L
    χ[:,1] = x
    wm[1] = λ/(L+λ)
    wc[1] = wm[1] + (1 - α^2 + β)
    wm[2:end] = 1/(2*(L+λ))
    wc[2:end] = wm[2:end]
    sp = sqrt(L+λ)*full(chol(p))
    for i = 2:L+1
        χ[:,i] = x + sp[:,i-1]
        χ[:,i+L] = x - sp[:,i-1]
    end
    return SigmaPoints(χ, wm, wc)
end
calcSigmaPoints(state::AbstractMvNormal, filter::UKF) = calcSigmaPoints(state, filter.α, filter.β, filter.κ)

"""
# timeUpdate
Function to calculate the predicted mean and the predicted covariance given a UKF model and the sigma points with their corresponding weights

`timeUpdate(m, sp)`
#### Parameters
- `m` : AbstractGaussianSSM type containing the parameters of the Unscented Kalman Filter model.
- `sp` : SigmaPoints type containing the matrix of sigma vectors and their corresponding weights
"""
function timeUpdate(m::AbstractGaussianSSM, sp::SigmaPoints, t::Real=0.0)
    L, M  = size(sp.χ)
    χ_x = zeros(L, M)
    x_pred = zeros(L)
    p_pred = zeros(L, L)
    for i in 1:2L+1
        χ_x[:,i] = m.f(sp.χ[:,i])
        x_pred += sp.wm[i] * χ_x[:,i]
    end
    for i in 1:2L+1
        p_pred += sp.wc[i] * (χ_x[:,i] - x_pred)*(χ_x[:,i] - x_pred)'
    end
    p_pred += m.V(t)
    return MvNormal(x_pred, p_pred), SigmaPoints(χ_x, sp.wm, sp.wc)
end

function predict(m::AbstractGaussianSSM, x::AbstractMvNormal, filter::UKF)
    sigPoints = calcSigmaPoints(x, filter)
    pred_state, new_sigPoints = timeUpdate(m, sigPoints)
    return pred_state, new_sigPoints
end

function observe(m::AbstractGaussianSSM, x::AbstractMvNormal, sp::SigmaPoints, y, t::Real=0.0)
    obsLength = length(y)
    L, M = size(sp.χ)
    y_trans = zeros(obsLength, M)
    y_pred = zeros(obsLength)
    for i in 1:2L+1
        y_trans[:,i] = m.g(sp.χ[:,i])
        y_pred += sp.wm[i] * y_trans[:,i]
    end

    P_xy = zeros(L, obsLength)
    P_yy = zeros(obsLength, obsLength)
    for i in 1:2L+1
        resy = (y_trans[:,i] - y_pred)
        P_xy += sp.wc[i] * (sp.χ[:,i] - mean(x)) * resy'
        P_yy += sp.wc[i] * resy * resy'
    end
    P_yy += m.W(t)
    return MvNormal(y_pred, P_yy), P_xy
end

function innovate(m::AbstractGaussianSSM, x::AbstractMvNormal, yPred::AbstractMvNormal, P_xy::Matrix, sp::SigmaPoints, y::Vector)
    kalmanGain = P_xy * inv(cov(yPred))
    new_x = mean(x) + kalmanGain * (y - mean(yPred))
    new_cov = cov(x) - kalmanGain * cov(yPred) * kalmanGain'
    return MvNormal(new_x, new_cov)
end

function update(m::AbstractGaussianSSM, x::AbstractMvNormal, sp::SigmaPoints, y::Vector)
    yPred, P_xy = observe(m, x, sp, y)
    return innovate(m, x, yPred, P_xy, sp, y)
end

function filter{T}(m::AbstractGaussianSSM, y::Array{T}, x0::AbstractMvNormal,
	filter::UKF=UKF())
    x_filtered = Array(AbstractMvNormal, size(y, 2))
    loglik = 0.0
	x_pred, sigma_points = predict(m, x0, filter)
	x_filtered[1] = update(m, x_pred, sigma_points, y[:, 1])
	for i in 2:size(y, 2)
		x_pred, sigma_points = predict(m, x_filtered[i-1], filter)
		# Check for missing values in observation
		if any(isnan(y[:, i]))
            x_filtered[i] = x_pred
        else
            x_filtered[i] = update(m, x_pred, sigma_points, y[:, i])
            loglik += logpdf(observe(m, x_filtered[i], calcSigmaPoints(x_filtered[i], filter.α, filter.β, filter.κ), y[:, 1])[1], y[:, 1])
        end
        loglik += logpdf(x_pred, mean(x_filtered[i]))
	end
	return FilteredState(y, x_filtered, loglik, false)
end

######################################################################
# Ensemble Kalman filter
######################################################################

type EnsembleKalmanFilter{I<:Integer} <: NonlinearFilter
    nparticles::I
end
typealias EnKF EnsembleKalmanFilter
EnsembleKalmanFilter() = EnsembleKalmanFilter(100)

function predict(m::AbstractGaussianSSM, ensemble::Matrix, filter::EnKF=EnKF();
        u::Vector=zeros(m.nu), t::Real=0.0)
    ensemble_new = similar(ensemble)
    CI = control_input(m, u, t)
    for i in 1:size(ensemble, 2)
        F = process_matrix(m, ensemble[:, i], t)
        ensemble_new[:, i] = F * ensemble[:, i] + CI
    end
    return ensemble_new
end

function update(m::AbstractGaussianSSM, ensemble::Matrix, y::Vector,
        filt::EnKF=EnKF(), t::Real=0.0)
    P = cov(ensemble')
    ensemble_updated = similar(ensemble)
    for i in 1:filt.nparticles
        G = observation_matrix(m, ensemble[:, i], t)
        W = observation_matrix(m, ensemble[:, i], t)
        innovation = y - G * ensemble[:, i]
        innovation_cov = G * P * G' + W
        K = P * G' * inv(innovation_cov)
        ensemble_updated[:, i] = ensemble[:, i] + K * innovation
        cov_update = (eye(P) - K * G) * P
    end
    return ensemble_updated
end

function filter{T}(m::AbstractGaussianSSM, y::Array{T}, x0::AbstractMvNormal,
        filt::EnKF=EnKF(); u::Matrix{T}=zeros(m.nu, size(y, 2)),
        times::Vector{T}=zeros(size(y, 2)), return_ensemble=false)
    nt = size(y, 2)
    ensemble = rand(x0, filt.nparticles)
    ensemble = predict(m, ensemble, filt, u=u[:,1], t=times[1])
    loglik = 0.0
    if return_ensemble
        x_filtered = zeros(m.nx, filt.nparticles, nt)
        x_filtered[:, :, 1] = update(m, ensemble, y[:, 1], filt, times[1])
        for i in 2:nt
            ensemble = predict(m, x_filtered[:, :, i-1], filt, u=u[:, i], t=times[i])
            x_filtered[:, :, i] = update(m, ensemble, y[:, i], filt, times[i])
        end
        return x_filtered
    else
        x_filtered = Array(AbstractMvNormal, size(y, 2))
        ensemble = update(m, ensemble, y[:, 1], filt, times[1])
        x_filtered[1] = MvNormal(vec(mean(ensemble, 2)), cov(ensemble'))
        for j in 1:filt.nparticles
            loglik += logpdf(x_filtered[1], ensemble[:, j])
        end
        if ! any(isnan(y[:, 1]))
            for j in 1:filt.nparticles
                loglik += logpdf(observe(m, x_filtered[1], times[1]), y[:,1])
            end
        end
        for i in 2:nt
            ensemble = predict(m, ensemble, filt, u=u[:, i], t=times[i])
            x_filtered[i] = MvNormal(vec(mean(ensemble, 2)), cov(ensemble'))
            for j in 1:filt.nparticles
                loglik += logpdf(x_filtered[i], ensemble[:, j])
            end
            if ! any(isnan(y[:, i]))
                for j in 1:filt.nparticles
                    loglik += logpdf(observe(m, x_filtered[1], times[1]), y[:,i])
                end
            end
        end
        return FilteredState(y, x_filtered, u, times, loglik, false)
    end
end


text:
module BioDataExpression

  using SQLite
  using Missings
  using StatsBase
  using MultivariateStats

  export SingleCellExpr, BioData, Normalize, Zscore, Log

  #=
  println("the comment has not worked!")
  basePath = dirname(dirname(Base.functionloc(BioDataExpression.eval, Tuple{Void})[1]))
  include(joinpath(basePath,"src","BioData/Constructors.jl"))
  include(joinpath(basePath,"src","BioData/Normalize.jl"))
  include(joinpath(basePath,"src","BioData/Zscore.jl"))
  include(joinpath(basePath,"src","BioData","Show.jl"))
  =#

  include("BioData/Constructors.jl")
  include("BioData/Normalize.jl")
  include("BioData/Zscore.jl")
  include("BioData/Log.jl")
  include("BioData/Show.jl")

end # module


text:
include("helper.jl")

my_slow_tests = [
  "problems/test_optimize_single_objective_problems.jl",
  "test_bboptimize.jl"
]

for t in my_slow_tests
  include(t)
end


text:
using Test
using LubyTransformCode

@testset "Encoder & Decoder" begin
    data = b"""A fountain (from the Latin "fons" (genitive "fontis"), a source or spring) is a piece
    of architecture which pours water into a basin or jets it into the air to supply drinking water
    and/or for a decorative or dramatic effect."""

    (drops, n_blocks) = encode(data)
    # `drops` is serialized and sent via network.
    decoded = decode(drops, n_blocks)
    if decoded != nothing
        @test decoded[1:length(data)] == data
    else
        println("Decoding failed.")
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
abstract type AbstractFOStateSpace end

mutable struct FOSS <: AbstractFOStateSpace
    a
    b
    c
    d
    α::Float64
    ioDelay
    E
    x0
end

foss(a, b, c, d, α, L, E, x0) =  FOSS(a, b, c, d, α, L, E, x0)

function fossdata(sys::FOSS)
    return [sys.a sys.b sys.c sys.d sys.E sys.ioDelay sys.x0 sys.α]
end

function minreal(G::FOSS)
    α = G.α
    G2 = foss2ss(G)
    G2 = minreal(G2)
end

# Convert fractional order state space to integer state space
function foss2ss(sys::FOSS)
    return ss(sys.a, sys.b, sys.c, sys.d)
end

# Show the created fractional order state space
Base.print(io::IO, sys::FOSS) = show(io, sys)

function Base.show(io::IO, sys::FOSS)
    println(typeof(sys))
    print("\n")
    println("A =")
    Base.print_matrix(io::IO, sys.a)
    print("\n")
    println("B =")
    Base.print_matrix(io::IO, sys.b)
    print("\n")
    println("C =")
    Base.print_matrix(io::IO, sys.c)
    print("\n")
    println("D =")
    Base.print_matrix(io::IO, sys.d)
    print("\n")
    println("\nDescriptor matrix:\n")
    print("E =\n")
    isa(sys.E, Matrix) ? Base.print_matrix(io::IO, sys.E) : print(sys.E)
    println("\n")
    println("Time delay is $(sys.ioDelay)")
    println("α = $(sys.α)")
    println("Initial state vector x₀ = $(sys.x0)")
end

"""
    ctrb(G::FOSS)

Controllability test matrix for FOSS
"""
ctrb(G::FOSS) = ControlSystems.ctrb(G.a, G.b)

"""
    obsv(G::FOSS)

Observability test matrix for FOSS
"""
obsv(G::FOSS) = ControlSystems.obsv(G.a, G.c)


text:
using GraphMatFun, LsqFit, LinearAlgebra, Random
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=6
target_n=200 # Measure error wrt this discr
n=50;        # Optimize wrt this discr
rho=2.7;  # SID

its=5;
droptol0=1e-10;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;
base.params[:kickit_mode]=1;

# First the optimization free
sid_org=init_state_mult!(deepcopy(base),:sid,m,showmeta=true)

#bbc_org=init_state_mult!(deepcopy(base),:bbc,m,showmeta=true)

#filename="simulations/newgraphs/exp_m4_sastre+GN_0_69.cgr"
sastre_org=init_state_mult!(deepcopy(base),:sastre,m,showmeta=true)

#sastre_org.cref=get_degopt_crefs(sastre_org.graph);

ps_org=init_state_mult!(deepcopy(base),:ps,m,showmeta=true)

mono_org=init_state_mult!(deepcopy(base),:mono,m,showmeta=true)

(mono,simlist,commandlist)=
     run_sequence(mono_org,"sknsssssssssddssddsssssssssssssdsddGGssssdddsGssdddsssdddssssddsssssssssssssddsssssddddsddddssssq");

(ps,simlist,commandlist)=
     run_sequence(ps_org, "kssnsKsssdddddsdddsGGsssGssdddssssssddssdddddsssddddsssddddsssssssssssddddddssGgsssssssq");

(sastre,simlist,commandlist)=
     run_sequence(sastre_org,
                  "sssnssNdddddddddddddssGGGsssddddddddsssdgggggssssssGGssssssssdsssGGssssddsddsddsssssssssq");

(sid,simlist,commandlist)=
     run_sequence(sid_org,
     "sssnssNdddddddddddddssGGGsssddsddsddsddsddddsdddsddddddssssddsssssddsssssddssssssq")

#(bbc,simlist,commandlist)=
#     run_sequence(bbc_org,
#     "sssssddssssddsssssssssssdddssdddssssdddssdddssssssssssdggdq")
#

include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module TestBsOptions

using Distributions
using Econometrics
using Base.Test

results = readcsv("/home/chris/research/julia/Econometrics/test/data/matlab_bsvals.csv")

include("/home/chris/research/julia/Econometrics/src/bsOptions.jl")

callPrices = [26.7, 824.9, 1166.7]
daxVals = [6175.05,
           6816.12,
           8455.83]
strikes = [6700, 6000, 7300]
r = [0.006,
     0.008,
     0.001]
T = [0.188,
     0.074,
     0.172]
sigmas = [0.4,
          0.2,
          0.2]

############
## prices ##
############

actOut = [bsCall(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
            for ii=1:3]

@test_approx_eq_eps actOut' results[1, :] 0.01

actOut = [bsPut(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
            for ii=1:3]

@test_approx_eq_eps actOut' results[2, :] 0.01

###########
## Delta ##
###########

actOut = [bsDeltaCall(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[3, :] 0.01

actOut = [bsDeltaPut(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[4, :] 0.01

###########
## Gamma ##
###########

actOut = [bsGamma(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[5, :] 0.01

##########
## Vega ##
##########

actOut = [bsVega(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[6, :] 0.01

###########
## Theta ##
###########

actOut = [bsThetaCall(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[7, :] 0.1

actOut = [bsThetaPut(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[8, :] 0.1

#########
## Rho ##
#########

actOut = [bsRhoCall(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[9, :] 0.1

actOut = [bsRhoPut(sigmas[ii], daxVals[ii], strikes[ii], r[ii], T[ii])
          for ii=1:3]

@test_approx_eq_eps actOut' results[10, :] 0.01

##################
## Implied vola ##
##################

sigma0 = 0.18
prec = 0.000001

volas = [implVolaCall(sigma0, callPrices[kk], daxVals[kk],
                     strikes[kk], r[kk], T[kk], prec)[1] for kk=1:3]

@test_approx_eq_eps volas' results[11, :] 0.01

end


text:
using Test
using LubyTransformCode

using LubyTransformCode: next_lehmer

@testset "Lehmer RNG" begin
    seeds = UInt32[
        0x00:0xFF;
        (0x7FFF_FFFF - 0xFF):(0x7FFF_FFFF + 0xFF);
        (0xFFFF_FFFF - 0xFF):0xFFFF_FFFF;
    ]

    for seed in seeds
        lehmer = Lehmer(seed)
        N = 5

        random_list = []
        for x in Iterators.take(lehmer, N)
            push!(random_list, x)
        end

        expected = []
        x = seed
        for i in 1:N
            x = next_lehmer(x)
            push!(expected, x)
        end

        @test random_list == expected
    end
end

@testset "Lehmer random sample: smoke test" begin
    seed = UInt32(42)
    lehmer = Lehmer(seed)
    n = UInt32(20)
    k = UInt32(10)
    samples = random_sample(lehmer, n, k)
    @test length(samples) == k && all(1 .<= samples .<= n)
end


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

using Test

@test true


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of headers and sub-headers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Headers and sub-headers" begin
    expected = """

    1
    2
    3
    4

    A
    B
    C
    D

    E
    F
    G
    H

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = Val(:html),
                          header = (["1", "2", "3", "4"],
                                    [:A,  :B,  :C,  :D],
                                    [:E,  :F,  :G,  :H]),
                          standalone = false)

    @test result == expected

    expected = """

    1
    2
    3
    4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = Val(:html),
                          header = (["1", "2", "3", "4"],
                                    [:A,  :B,  :C,  :D],
                                    [:E,  :F,  :G,  :H]),
                          nosubheader = true,
                          standalone = false)

    @test result == expected

    expected = """

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data;
                          backend = Val(:html),
                          header = (["1", "2", "3", "4"],
                                    [:A,  :B,  :C,  :D],
                                    [:E,  :F,  :G,  :H]),
                          noheader = true,
                          standalone = false)

    @test result == expected
end


text:
#
"""
 Functional based mesh adaptation

 FWD Pass: NonLin layer -> least sq. on linear layer
 BWD Pass: GD all the way through (takes adjoint on least sq.)

 Task: replace least sq. with PDE solve
       test problem: steady advection diffusion eqn at large
                     Peclet # (thin thin boundary layer)

"""

using Revise
using diffMesh
import Zygote
using NNlib, Flux, LinearAlgebra, Plots
#--------------------------------------#
np  = 2001 # data points
ng  = 31   # FEM nodes

function udata(x)
    u   = 2 .+ zero(x)
    u .+= @. sin(1*pi*x) * (-1 <= x <= 0)
    u .+= @. sin(4*pi*x) * ( 0 <  x <= 1)
    return u
end
#--------------------------------------#
function addEndPts(xe::AbstractArray
                  ,x0::Number,x1::Number)

    n = length(xe) + 2
    A,R = Zygote.ignore() do
        A = zeros(n,2)
        A[1,1] = 1
        A[n,2] = 1
        R = restrict(n)
        return A,R
    end

    xg  = A * [x0,x1] + R'*xe

    return xg
end
#--------------------------------------#
function model(xe,xp,up) # grid points, eval points, values to fit

    xg = addEndPts(xe,-1,1)
    NN = reluHatFunc(xg)
    ug = llsq(NN,xp,up)

    return reluHat(NN,ug,xp)
end
#--------------------------------------#
function llsq(NN,xp,up)

    A  = NN(xp)
    linOp(u) = reluHat(NN,u,xp)

#   u = udata(x)
    u = A \ up
    return u
end
#--------------------------------------#
function loss()
    fp   = model(xe,xp,up)
    loss = Flux.Losses.mse(fp,up)
    return loss
end
#--------------------------------------#
function makeplot()
    fp = model(xe,xp,up)

    xg = addEndPts(xe,-1,1)
    NN = reluHatFunc(xg)
    ug = llsq(NN,xp,up)

    plt = plot()
    plt = plot!(title="FEM AutoMesh",xlabel="x",ylabel="y")

    plt = plot!(xp,up,width=3,label="data")
    plt = plot!(xp,fp,width=3,label="prediction")
    plt = scatter!(xg,ug,width=4,color="red",label=false)

    # hat functions
    NNp = NN(xp')
    for i=1:ng
        scale = ug[i]
#       scale = 1.0
        global plt = plot!(xp,scale .* NNp[i,:]
                          ,width=2,color="black",label=false)
    end

    #savefig(plt,"diffMesh.png")
    display(plt)
    return plt
end
#--------------------------------------#
function callback(;doplot=true)

    global xe .= sort(xe)
    l2 = loss()
    fp = model(xe,xp,up)
    er = norm(up .- fp,Inf)
    println("Error ∞ norm: $er, loss: $l2")

    if(doplot) makeplot() end
    return false
end
#--------------------------------------#

# data
xp = Array(range(-1,stop=1,length=np))'
up = udata(xp)

# initialize
xg = Array(range(-1,stop=1,length=ng))
xg = sort(2*rand(ng).-1)

xe = xg[2:end-1]
fp = model(xe,xp,up)

callback(doplot=true)

ps = Flux.params(xe)
data = Iterators.repeated((),100)
opt = ADAM(1e-2)

gs = gradient(()->loss(),ps)

#anim = Animation()
for d in data
    global gs = gradient(()->loss(),ps)
    Flux.Optimise.update!(opt,ps,gs)
    callback(doplot=true)
    #frame(anim)
end
#gif(anim,"examples/AutoMesh.gif",fps=30)

#--------------------------------------#
nothing


text:
export simplify_coords, normalize_identified_coords

function simplify_coords(coords)
    new_coords = [ coords[1] ]
    for coord_ind in 2:length(coords)-1
        if (coords[coord_ind][2] != coords[coord_ind+1][2]) || (coords[coord_ind][2] != coords[coord_ind-1][2])
            push!(new_coords, coords[coord_ind])
        end
    end
    push!(new_coords, coords[end])
    return new_coords
end

function normalize_identified_coords(coords_nbman_id, opt_supp)
    supp_size = length(opt_supp)

    return normalized_coords = [ (t[1], 100 * t[2]/supp_size) for t in coords_nbman_id]
end

function runexpnums(problems, algorithms; FIGS_FOLDER = "./figs/expnums")

    ## Start from interesting point and optimize with ISTA, FISTA.
    problemclass_to_algstats = OrderedDict()
    for (problems_class_name, problem_class) in problems
        printstyled("\n---- Class $(problems_class_name)\n", color=:light_blue)

        algo_to_stats = OrderedDict()
        for (algo_id, algo) in enumerate(algorithms)
            printstyled("\n---- Algorithm $(algo.name)\n", color=:light_blue)

            algo_to_stats[algo] = Dict()
            for (pb_id, pb_data) in enumerate(problem_class)
                printstyled("\n----- Solving problem $pb_id, $(pb_data.name)\n", color=:light_blue)

                xopt, hist = solve_proxgrad(pb_data.pb, pb_data.xstart, algo.updatefunc; algo.params...)
                algo_to_stats[algo][pb_data] = (hist=hist, xopt=xopt)
            end
        end

        problemclass_to_algstats[problems_class_name] = algo_to_stats
    end

    ####################################################
    ## Display suboptimality, iterates position.
    basename(pwd()) == "src" && (FIGS_FOLDER = joinpath("..", FIGS_FOLDER))

    for (pbclass_name, algo_to_stats) in problemclass_to_algstats
        println("Generating graphs for class ", pbclass_name)

        ps = []

        colors_vec = distinguishable_colors(8) # transform=deuteranopic

        algo_id = 1
        for (algo, stats) in algo_to_stats
            @show algo.name

            ############################################
            ## Plot one for legend
            (pb, pbinstance_stats) = first(stats); delete!(stats, pb)

            reg = typeof(pb.pb).parameters[1]
            opt_supp = get_support(reg, pb.pb.x0)

            coords_nbman_y_identified = simplify_coords([ (k, length(intersect(v, opt_supp))) for (k, v) in pbinstance_stats.hist[:iter_support_y]])
            coords_nbman_y_identified = normalize_identified_coords(coords_nbman_y_identified, opt_supp)

            push!(ps, PlotInc(
                PGFPlotsX.Options(
                    "no marks" => nothing,
                    # "line width" => "0.3pt",
                    # "color" => colors_vec[algo_id],
                    get_iterates_algoparams(algo.name)...
                ),
                Coordinates(coords_nbman_y_identified)
            ))

            for (pb, pbinstance_stats) in stats
                reg = typeof(pb.pb).parameters[1]
                opt_supp = get_support(reg, pb.pb.x0)

                coords_nbman_y_identified = simplify_coords([ (k, length(intersect(v, opt_supp))) for (k, v) in pbinstance_stats.hist[:iter_support_y]])
                coords_nbman_y_identified = normalize_identified_coords(coords_nbman_y_identified, opt_supp)

                push!(ps, PlotInc(
                    PGFPlotsX.Options(
                        "no marks" => nothing,
                        "forget plot" => nothing,
                        # "line width" => "0.3pt",
                        # "color" => colors_vec[algo_id],
                        get_iterates_algoparams(algo.name)...
                    ),
                    Coordinates(coords_nbman_y_identified)
                ))

            end
            push!(ps, LegendEntry(get_legend(algo.name)))

            algo_id+=2
        end

        fig = @pgf Axis(
            {
                # height = "12cm",
                # width = "12cm",
                # xlabel = "iter",
                # ylabel = L"F(x_k)-F^\star",
                # ymode  = "log",
                xmode  = "log",
                # title = "Problem $(pb.name) -- Suboptimality",
                legend_pos = "north west",
                legend_cell_align = "left",
                legend_style = "font=\\footnotesize",
            },
            ps...
        )
        try
            pgfsave(joinpath(FIGS_FOLDER, "reg_$(pbclass_name)_identified.pdf"), fig)
        catch e
            println("Error while building pdf:")
            println(e)
        end
        pgfsave(joinpath(FIGS_FOLDER, "reg_$(pbclass_name)_identified.tikz"), fig)
        ##################################
    end

    return
end


text:
@testset "Memory: Static" begin

function memory_static_kernel(a,b)
    # Local
    ptr_local = alloc_special(Val(:local), Float32, Val(AS.Local), Val(1))
    unsafe_store!(ptr_local, a[1])
    b[1] = unsafe_load(ptr_local)

    # Region
    #= TODO: AMDGPU target cannot select
    ptr_region = alloc_special(Val(:region), Float32, Val(AS.Region), Val(1))
    unsafe_store!(ptr_region, a[2])
    b[2] = unsafe_load(ptr_region)
    =#

    # Private
    #= TODO
    ptr_private = alloc_special(Val(:private), Float32, Val(AS.Private), Val(1))
    unsafe_store!(ptr_private, a[3])
    b[3] = unsafe_load(ptr_private)
    =#

    nothing
end

A = ones(Float32, 1)
B = zeros(Float32, 1)

HA = HSAArray(A)
HB = HSAArray(B)

wait(@roc memory_static_kernel(HA, HB))

@test Array(HA) ≈ Array(HB)

end


text:
using Test
using CUDA
using CUDAKernels
using KernelAbstractions

# Function to use as a baseline for CPU metrics
function histogram(input)
    histogram_output = zeros(Int, maximum(input))
    for i = 1:length(input)
        histogram_output[input[i]] += 1
    end
    return histogram_output
end

# This is for a 1D histogram, for n-dim binning, we need to think a bit more
# The problem here is that the @synchronize calls are only synchronizing across
#     blocks, while we need a way to combine all the shmem histograms together
#     into a big global blob. This can be done by:
#     1. making sure we use lid instead of tid for the histogram at the end
#        and use atomic adds
#     2. lid doesn't make sense on the CPU with only 4 cores.
#     3. fix error with jl_apply_type for GPU
#     4. Think about svectors
#     5. dynamic shmem allocation is not allowed, so use threadsize or a const
#     6. We need to find a way to do a coalesced write to global memory
#         - I think this will be easier on the GPU when more threads are avail
@kernel function histogram_tag_kernel!(histogram_output, input)
    tid = @index(Global, Linear)
    #lid = @index(Local, Linear)

    #shared_histogram = @localmem Int32 length(histogram_output)
    shared_histogram = @localmem Int64 128

    # Note: this is 4 for CPU
    #@print("LID is: ", lid, '\n')
    @print("TID is: ", tid, '\n')

    # TODO: anything but this!
    if tid <= 128
        @inbounds shared_histogram[tid] = 0
    end

    @synchronize()

    #@print("values are: ", input[tid], "\n")

    if tid < length(histogram_output)
        @print("CHECK\n")
    end

    if tid <= length(input)
        @inbounds temp = shared_histogram[input[tid]]
        @inbounds shared_histogram[input[tid]] = temp + 1
        @print("input is: ", input[tid], '\n')
        @print("hist is: ", shared_histogram[input[tid]], "\n\n")
    end

    @synchronize()

    @print(length(histogram_output), '\n')

    if tid <= length(histogram_output)
        @print("histogram values: ", shared_histogram[tid], '\n')
        @inbounds histogram_output[input[tid]] += shared_histogram[input[tid]]
    end

    @synchronize()
end

function histogram_tag(input; numcores = 4, numthreads = 256)
    # I don't know how maximum is optimized on the GPU
    histogram_output = zeros(Int64, maximum(input))
    if isa(input, Array)
        event = histogram_tag!(histogram_output, input;
                               numcores = numcores, numthreads = numthreads)
        wait(event)
    else
        event = histogram_tag!(CuArray(histogram_output), input;
                               numcores = numcores, numthreads = numthreads)
        wait(event)
    end
    return histogram_output
end

function histogram_tag!(histogram_output, input;
                        numcores = 4, numthreads = 256)

    if isa(input, Array)
        kernel! = histogram_tag_kernel!(CPU(), numcores)
    else
        kernel! = histogram_tag_kernel!(CUDADevice(), numthreads)
    end

    kernel!(histogram_output, input, ndrange=size(input))
end

#=
@testset "histogram tests" begin

    a = [Int64(rand(1:128)) for i = 1:1000]
    final_histogram = histogram(a)

    CPU_tag_histogram = histogram_tag(a)

    @test isapprox(CPU_tag_histogram, final_histogram)

    if has_cuda_gpu()
        CUDA.allowscalar(false)
        d_a = CuArray(a)
        GPU_tag_histogram = histogram_tag(d_a)
        @test isapprox(Array(GPU_tag_histogram), final_histogram)
    end
end
=#

max_stars_repo_name: leios/simuleios

max_issues_repo_name: leios/simuleios

max_forks_repo_name: leios/simuleios


text:
module MouseEventTypes
    @enum MouseEventType begin
        out
        enter
        over
        leftdown
        rightdown
        middledown
        leftup
        rightup
        middleup
        leftdragstart
        rightdragstart
        middledragstart
        leftdrag
        rightdrag
        middledrag
        leftdragstop
        rightdragstop
        middledragstop
        leftclick
        rightclick
        middleclick
        leftdoubleclick
        rightdoubleclick
        middledoubleclick
        downoutside
    end
    export MouseEventType
end

using .MouseEventTypes

"""
    MouseEvent

Describes a mouse state change.
Fields:
- `type`: MouseEventType
- `t`: Time of the event
- `data`: Mouse position in data coordinates
- `px`: Mouse position in px relative to scene origin
- `prev_t`: Time of previous event
- `prev_data`: Previous mouse position in data coordinates
- `prev_px`: Previous mouse position in data coordinates
"""
struct MouseEvent
    type::MouseEventType
    t::Float64
    data::Point2f0
    px::Point2f0
    prev_t::Float64
    prev_data::Point2f0
    prev_px::Point2f0
end

struct MouseEventHandle
    obs::Makie.PriorityObservable{MouseEvent}
    observerfuncs::Vector{<:Observables.ObserverFunction}
end

"""
    clear!(handle::MouseEventHandle)

Cut observable connections to the scene and remove any listeners to the mouse events.
"""
function clear!(handle::MouseEventHandle)
    foreach(Observables.off, handle.observerfuncs)
    empty!(handle.observerfuncs)
    empty!(handle.obs.listeners)
    nothing
end

for eventtype in instances(MouseEventType)
    onfunctionname = Symbol("onmouse" * String(Symbol(eventtype)))
    @eval begin

        """
        Executes the function f whenever the `MouseEventHandle`'s observable is set to
        a MouseEvent with `event.type === $($eventtype)`.
        """
        function $onfunctionname(f, mev::MouseEventHandle; priority = Int8(0))
            sanitized_func = _sanitize_observer_function(f, (typeof(mev.obs.val),))
            on(mev.obs, priority = priority) do event
                if event.type === $eventtype
                    return sanitized_func(event)::Bool
                end
                return false
            end
        end
        export $onfunctionname
    end
end

"""
    addmouseevents!(scene, elements...)

Returns a `MouseEventHandle` with an observable inside which is triggered by all mouse
interactions with the `scene` and optionally restricted to all given
plot objects in `elements`.

To react to mouse events, use the onmouse... handlers.

Example:

```
mouseevents = addmouseevents!(scene, scatterplot)

onmouseleftclick(mouseevents) do event
    # do something with the mouseevent
end
```
"""
function addmouseevents!(scene, elements...; priority = Int8(1))
    is_mouse_over_relevant_area() = isempty(elements) ? Makie.is_mouseinside(scene) : mouseover(scene, elements...)
    _addmouseevents!(scene, is_mouse_over_relevant_area, priority)
end
function addmouseevents!(scene, bbox::Observables.AbstractObservable{<: Rect2D}; priority = Int8(1))
    is_mouse_over_relevant_area() = Makie.mouseposition_px(scene) in bbox[]
    _addmouseevents!(scene, is_mouse_over_relevant_area, priority)
end

function _addmouseevents!(scene, is_mouse_over_relevant_area, priority)
    Mouse = Makie.Mouse
    dblclick_max_interval = 0.2

    mouseevent = Makie.PriorityObservable{MouseEvent}(
        MouseEvent(MouseEventTypes.out, 0.0, Point2f0(0, 0), Point2f0(0, 0), 0.0, Point2f0(0, 0), Point2f0(0, 0))
    )

    # initialize state variables
    last_mouseevent = Ref{Mouse.Action}(Mouse.release)
    prev_data = Ref(mouseposition(scene))
    prev_px = Ref(Makie.mouseposition_px(scene))
    mouse_downed_inside = Ref(false)
    mouse_downed_button = Ref{Optional{Mouse.Button}}(nothing)
    drag_ongoing = Ref(false)
    mouse_was_inside = Ref(false)
    prev_t = Ref(0.0)
    t_last_click = Ref(0.0)
    b_last_click = Ref{Optional{Mouse.Button}}(nothing)
    last_click_was_double = Ref(false)

    # react to mouse position changes
    mousepos_observerfunc = on(events(scene).mouseposition, priority=priority) do mp
        consumed = false
        t = time()
        data = mouseposition(scene)
        px = Makie.mouseposition_px(scene)
        mouse_inside = is_mouse_over_relevant_area()

        # last_mouseevent can only be up or down

        # mouse moved while being pressed
        # this can mean a new drag started, or a drag continues if it is ongoing.
        # it can also mean that a drag that started outside and isn't related to this
        # object is going across it and should be ignored here
        if last_mouseevent[] == Mouse.press

            if drag_ongoing[]
                # continue the drag
                event = @match mouse_downed_button[] begin
                    Mouse.left => MouseEventTypes.leftdrag
                    Mouse.right => MouseEventTypes.rightdrag
                    Mouse.middle => MouseEventTypes.middledrag
                    x => error("No recognized mouse button $x")
                end
                x = setindex!(mouseevent,
                    MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                )
                consumed = consumed || x
            else
                # mouse was downed inside but no drag is ongoing
                # that means a drag started
                if mouse_downed_inside[]
                    drag_ongoing[] = true
                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseEventTypes.leftdragstart
                        Mouse.right => MouseEventTypes.rightdragstart
                        Mouse.middle => MouseEventTypes.middledragstart
                        x => error("No recognized mouse button $x")
                    end
                    x = setindex!(mouseevent,
                        MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x

                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseEventTypes.leftdrag
                        Mouse.right => MouseEventTypes.rightdrag
                        Mouse.middle => MouseEventTypes.middledrag
                        x => error("No recognized mouse button $x")
                    end
                    x = setindex!(mouseevent,
                        MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x
                end
            end
        else
            if mouse_inside
                x = if mouse_was_inside[]
                    setindex!(mouseevent,
                        MouseEvent(MouseEventTypes.over, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                else
                    setindex!(mouseevent,
                        MouseEvent(MouseEventTypes.enter, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                end
                consumed = consumed || x
            else
                if mouse_was_inside[]
                    x = setindex!(mouseevent,
                        MouseEvent(MouseEventTypes.out, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x
                end
            end
        end

        mouse_was_inside[] = mouse_inside
        prev_data[] = data
        prev_px[] = px
        prev_t[] = t
        return consumed::Bool
    end

    # react to mouse button changes
    mousedrag_observerfunc = on(events(scene).mousebutton, priority=priority) do event
        consumed = false
        t = time()
        data = prev_data[]
        px = prev_px[]

        # TODO: this could probably be simplified by just using event.button
        # though that would probably change the way this handles a bit
        pressed_buttons = events(scene).mousebuttonstate

        # mouse went down, this can either happen inside or outside the objects of interest
        # we also only react if one button is pressed, because otherwise things go crazy (pressed left button plus clicks from other buttons in between are not allowed, e.g.)
        if event.action == Mouse.press
            if length(pressed_buttons) == 1
                button = first(pressed_buttons)
                mouse_downed_button[] = button

                if mouse_was_inside[]
                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseEventTypes.leftdown
                        Mouse.right => MouseEventTypes.rightdown
                        Mouse.middle => MouseEventTypes.middledown
                        x => error("No recognized mouse button $x")
                    end
                    x = setindex!(mouseevent,
                        MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x
                    mouse_downed_inside[] = true
                else
                    mouse_downed_inside[] = false
                    x = setindex!(mouseevent,
                        MouseEvent(MouseEventTypes.downoutside, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x
                end
            end
            last_mouseevent[] = Mouse.press
        elseif event.action == Mouse.release
            # only register up events and clicks if the upped button matches
            # the recorded downed one
            # and it can't be nothing (if the first up event comes from outside)

            downed_button_missing_from_pressed = !(mouse_downed_button[] in pressed_buttons)
            some_mouse_button_had_been_downed = !isnothing(mouse_downed_button[])

            if downed_button_missing_from_pressed && some_mouse_button_had_been_downed

                if drag_ongoing[]
                    event = @match mouse_downed_button[] begin
                        Mouse.left => MouseEventTypes.leftdragstop
                        Mouse.right => MouseEventTypes.rightdragstop
                        Mouse.middle => MouseEventTypes.middledragstop
                        x => error("No recognized mouse button $x")
                    end
                    x = setindex!(mouseevent,
                        MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                    )
                    consumed = consumed || x
                    drag_ongoing[] = false

                    if mouse_was_inside[]
                        # up after drag done over element
                        event = @match mouse_downed_button[] begin
                            Mouse.left => MouseEventTypes.leftup
                            Mouse.right => MouseEventTypes.rightup
                            Mouse.middle => MouseEventTypes.middleup
                            x => error("No recognized mouse button $x")
                        end

                        x = setindex!(mouseevent,
                            MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                        )
                        consumed = consumed || x
                    else
                        # mouse could be not over elements after drag is over
                        x = setindex!(mouseevent,
                            MouseEvent(MouseEventTypes.out, t, data, px, prev_t[], prev_data[], prev_px[])
                        )
                        consumed = consumed || x
                    end
                else
                    if mouse_was_inside[]
                        dt_last_click = t - t_last_click[]
                        t_last_click[] = t

                        # guard against mouse coming in from outside, then mouse upping
                        if mouse_downed_inside[]
                            if dt_last_click < dblclick_max_interval && !last_click_was_double[] &&
                                    mouse_downed_button[] == b_last_click[]

                                event = @match mouse_downed_button[] begin
                                    Mouse.left => MouseEventTypes.leftdoubleclick
                                    Mouse.right => MouseEventTypes.rightdoubleclick
                                    Mouse.middle => MouseEventTypes.middledoubleclick
                                    x => error("No recognized mouse button $x")
                                end
                                x = setindex!(mouseevent,
                                    MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                                )
                                consumed = consumed || x
                                last_click_was_double[] = true
                            else
                                event = @match mouse_downed_button[] begin
                                    Mouse.left => MouseEventTypes.leftclick
                                    Mouse.right => MouseEventTypes.rightclick
                                    Mouse.middle => MouseEventTypes.middleclick
                                    x => error("No recognized mouse button $x")
                                end
                                x = setindex!(mouseevent,
                                    MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                                )
                                consumed = consumed || x
                                last_click_was_double[] = false
                            end
                            # save what type the last downed button was
                            b_last_click[] = mouse_downed_button[]
                        end
                        mouse_downed_inside[] = false

                        # up after click
                        event = @match mouse_downed_button[] begin
                            Mouse.left => MouseEventTypes.leftup
                            Mouse.right => MouseEventTypes.rightup
                            Mouse.middle => MouseEventTypes.middleup
                            x => error("No recognized mouse button $x")
                        end

                        x = setindex!(mouseevent,
                            MouseEvent(event, t, data, px, prev_t[], prev_data[], prev_px[])
                        )
                        consumed = consumed || x
                    end
                end
            end

            last_mouseevent[] = Mouse.release
        end

        prev_t[] = t
        return consumed::Bool
    end

    MouseEventHandle(mouseevent, [mousepos_observerfunc, mousedrag_observerfunc])
end


text:
using Oceananigans
using Oceananigans.Units
using GLMakie
using IcyAntics
using IcyAntics: time_step!

#####
##### Domain: 2D grid with 100km square extent, ~500m resolution
#####

arch = CPU() # change to GPU() to run coupled simulation on GPU.
Nx = Ny = 256
Lx = Ly = 100kilometers

grid = RectilinearGrid(arch,
                       size = (Nx, Ny),
                       halo = (3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       topology = (Periodic, Periodic, Flat))

#####
##### Barotropic turbulence ocean simulation
#####

# Ocean mesoscale closure: biharmonic diffusivity
Δh = grid.Δxᶜᵃᵃ
closure = HorizontalScalarBiharmonicDiffusivity(ν = Δh^4 / 4hour)

ocean_model = NonhydrostaticModel(; grid, closure,
                                  timestepper = :RungeKutta3,
                                  advection = WENO5(),
                                  buoyancy = nothing,
                                  tracers = nothing)

# Spin up the ocean model
ϵ(x, y, z) = randn()
set!(ocean_model, u=ϵ, v=ϵ) # random initial velocities
ocean_simulation = Simulation(ocean_model, Δt=2minutes, stop_time=12hours)
run!(ocean_simulation) # generate smooth ocean velocities
ocean_simulation.stop_time = Inf

#####
##### Simple viscoelastic ice model
#####

ocean_state = (; u=ocean_model.velocities.u, v=ocean_model.velocities.v, ρ=1024.0)
rheology = ViscoElasticRheology(modulus=1.0, viscosity=1e4)
@show ice_Δt = 1e-1 * min(Δh / sqrt(rheology.modulus), Δh^2 / (2 * rheology.viscosity))

# Spin up the ice model
ice_model = ContinuumIceModel(; grid, ocean=ocean_state, rheology)
ice_simulation = Simulation(ice_model, Δt=ice_Δt)

#####
##### Construct coupled simulation with substepping
#####

function coupled_time_step!(ocean_simulation, ice_simulation)
    # Floor ice_simulation Δt so `Nsubsteps` ice steps = one ocean step.
    Nsubsteps = ceil(Int, ocean_simulation.Δt / ice_simulation.Δt)
    ice_simulation.Δt = ocean_simulation.Δt / Nsubsteps

    # Substep ice model
    for substep = 1:Nsubsteps
        time_step!(ice_simulation)
    end

    time_step!(ocean_simulation)

    return nothing
end

#####
##### Visualization and coupled run
#####

n = Observable(1) # for visualization

# Ocean vorticity
uo, vo, wo = ocean_model.velocities
ocean_vorticity = Field(∂x(vo) - ∂y(uo))
ocean_vorticity_n = @lift ($n; interior(compute!(ocean_vorticity))[:, :, 1])

# Ice speed
ui, vi = ice_model.velocities
ice_speed = Field(sqrt(ui^2 + vi^2))
ice_speed_n = @lift ($n; interior(compute!(ice_speed))[:, :, 1])

# Make figure
x, y, z = nodes((Center, Center, Center), grid)
fig = Figure(resolution=(800, 400))
ax_i = Axis(fig[1, 1], aspect=1, xlabel="x (km)", ylabel="y (km)", title="Ice speed")
ax_o = Axis(fig[1, 2], aspect=1, xlabel="x (km)", ylabel="y (km)", title="Ocean vorticity")
hm_i = heatmap!(ax_i, 1e-3 * x, 1e-3 * y, ice_speed_n)
hm_o = heatmap!(ax_o, 1e-3 * x, 1e-3 * y, ocean_vorticity_n, colormap=:redblue)
title = @lift ($n; "Coupled ice-ocean simulation after " * prettytime(ice_model.clock.time))
Label(fig[0, :], title)
display(fig)

# Step forward and animate
N = 20
record(fig, "viscoelastic_ice.gif", 1:N, framerate=12) do i
    @info "Plotting frame $i of $N..."
    @time for _ = 1:10
        coupled_time_step!(ocean_simulation, ice_simulation)
    end
    n[] = 1
end


text:
using JuliaFormatter
using JuliaFormatter: DefaultStyle, YASStyle, Options, options, CONFIG_FILE_NAME
using CSTParser
using Test

fmt1(s; i = 4, m = 80, kwargs...) =
    JuliaFormatter.format_text(s; kwargs..., indent = i, margin = m)
fmt1(s, i, m; kwargs...) = fmt1(s; kwargs..., i = i, m = m)

# Verifies formatting the formatted text
# results in the same output
function fmt(s; i = 4, m = 80, kwargs...)
    kws = merge(options(DefaultStyle()), kwargs)
    s1 = fmt1(s; kws..., i = i, m = m)
    return fmt1(s1; kws..., i = i, m = m)
end
fmt(s, i, m; kwargs...) = fmt(s; kwargs..., i = i, m = m)

function run_pretty(text::String; style = DefaultStyle(), opts = Options())
    d = JuliaFormatter.Document(text)
    s = JuliaFormatter.State(d, opts)
    x = CSTParser.parse(text, true)
    t = JuliaFormatter.pretty(style, x, s)
    t
end
run_pretty(text::String, margin::Int) = run_pretty(text, opts = Options(margin = margin))

function run_nest(text::String; opts = Options(), style = DefaultStyle())
    d = JuliaFormatter.Document(text)
    s = JuliaFormatter.State(d, opts)
    x = CSTParser.parse(text, true)
    t = JuliaFormatter.pretty(style, x, s)
    JuliaFormatter.nest!(style, t, s)
    t, s
end
run_nest(text::String, margin::Int) = run_nest(text, opts = Options(margin = margin))

yasfmt1(str) = fmt1(str; style = YASStyle(), options(DefaultStyle())...)
yasfmt(str, i = 4, m = 80; kwargs...) = fmt(str, i, m; style = YASStyle(), kwargs...)

bluefmt1(str) = fmt1(str; style = BlueStyle(), options(DefaultStyle())...)
bluefmt(str, i = 4, m = 80; kwargs...) = fmt(str, i, m; style = BlueStyle(), kwargs...)

@testset "JuliaFormatter" begin
    include("default_style.jl")
    include("yas_style.jl")
    include("blue_style.jl")
    include("issues.jl")
    include("options.jl")
    include("config.jl")
    include("document.jl")
    include("interface.jl")
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/power.jl
# Contains definitions of inv and ^.
#############################################################################

# defines square operator
@inline sqr(x::Float64) = x*x
@inline cv_sqr_NS(x::Float64, xL::Float64, xU::Float64) = x*x
@inline dcv_sqr_NS(x::Float64, xL::Float64, xU::Float64) = 2.0*x
@inline function cc_sqr(x::Float64, xL::Float64, xU::Float64)
	if (xU > xL)
		cc = xL^2 + (xL + xU)*(x - xL)
	else
		cc = xU^2
	end
	return cc
end
@inline dcc_sqr(x::Float64, xL::Float64, xU::Float64) = (xU > xL) ? (xL + xU) : 0.0
@inline function cv_sqr(x::Float64, xL::Float64, xU::Float64)
    (0.0 <= xL || xU <= 0.0) && return x^2
	((xL < 0.0) && (0.0 <= xU) && (0.0 <= x)) && return (x^3)/xU
	return (x^3)/xL
end
@inline function dcv_sqr(x::Float64, xL::Float64, xU::Float64)
    (0.0 <= xL || xU <= 0.0) && return 2.0*x
	((xL < 0.0) && (0.0 <= xU) && (0.0 <= x)) && return (3.0*x^2)/xU
	return (3.0*x^2)/xL
end
@inline function sqr_kernel(x::MC{N,T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
	if (x.Intv.hi < 0.0)
       eps_min = x.Intv.hi
       eps_max = x.Intv.lo
    elseif (x.Intv.lo > 0.0)
       eps_min = x.Intv.lo
       eps_max = x.Intv.hi
    else
       eps_min = 0.0
       eps_max = (abs(x.Intv.lo) >= abs(x.Intv.hi)) ? x.Intv.lo : x.Intv.hi
    end
	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
	cc = cc_sqr(midcc, x.Intv.lo, x.Intv.hi)
	dcc = dcc_sqr(midcc, x.Intv.lo, x.Intv.hi)
	cv = cv_sqr_NS(midcv, x.Intv.lo, x.Intv.hi)
	dcv = dcv_sqr_NS(midcv, x.Intv.lo, x.Intv.hi)
	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function sqr_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N
	if (x.Intv.hi < 0.0)
       eps_min = x.Intv.hi
       eps_max = x.Intv.lo
    elseif (x.Intv.lo > 0.0)
       eps_min = x.Intv.lo
       eps_max = x.Intv.hi
    else
       eps_min = 0.0
       eps_max = (abs(x.Intv.lo) >= abs(x.Intv.hi)) ? x.Intv.lo : x.Intv.hi
    end
	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
	cc = cc_sqr(midcc, x.Intv.lo, x.Intv.hi)
    cv = cv_sqr(midcv, x.Intv.lo, x.Intv.hi)
    gdcv1 = dcv_sqr(x.cv, x.Intv.lo, x.Intv.hi)
    gdcv2 = dcv_sqr(x.cc, x.Intv.lo, x.Intv.hi)
    gdcc1 = dcc_sqr(x.cv, x.Intv.lo, x.Intv.hi)
    gdcc2 = dcc_sqr(x.cc, x.Intv.lo, x.Intv.hi)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline sqr(x::MC) = sqr_kernel(x, pow(x.Intv,2))

# convex/concave relaxation (Khan 3.1-3.2) of integer powers of 1/x for positive reals
@inline pow_deriv(x::Float64, c::T) where {T <: AbstractFloat} = c*x^(c-1)
@inline pow_deriv(x::Float64, n::Z) where {Z <: Integer} = n*x^(n-1)
@inline function cv_npp_or_pow4(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
	temp = x^(n-1)
	x*temp, n*temp
end
@inline function cc_npp_or_pow4(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
	dline_seg(^, pow_deriv, x, xL, xU, n)
end
# convex/concave relaxation of integer powers of 1/x for negative reals
@inline function cv_negpowneg(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    isodd(n) && (return dline_seg(^, pow_deriv, x, xL, xU, n))
    return x^n, n*x^(n-1)
end
@inline function cc_negpowneg(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    isodd(n) && (return x^n,n*x^(n-1))
    return dline_seg(^, pow_deriv, x, xL, xU, n)
end
# convex/concave relaxation of odd powers
@inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n))
    (0.0 <= xL) && (return x^n, n*x^(n - 1))
    val = (xL^n)*(xU - x)/(xU - xL) + (max(0.0, x))^n
    dval = -(xL^n)/(xU - xL) + n*(max(0.0, x))^(n-1)
    return val, dval
end
@inline function cc_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    (xU <= 0.0) && (return x^n, n*x^(n - 1))
    (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, n))
    val = (xU^n)*(x - xL)/(xU - xL) + (min(0.0, x))^n
    dval = (xU^n)/(xU - xL) + n*(min(0.0, x))^(n-1)
    return val, dval
end

@inline function npp_or_pow4(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}}
    if (x.Intv.hi < 0.0)
        eps_min = x.Intv.hi
        eps_max = x.Intv.lo
    elseif (x.Intv.lo > 0.0)
        eps_min = x.Intv.lo
        eps_max = x.Intv.hi
    else
        eps_min = 0.0
        eps_max = (abs(x.Intv.lo) >= abs(x.Intv.hi)) ? x.Intv.lo : x.Intv.hi
  	end
  	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
  	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
  	cc, dcc = cc_npp_or_pow4(midcc, x.Intv.lo, x.Intv.hi, c)
  	cv, dcv = cv_npp_or_pow4(midcv, x.Intv.lo, x.Intv.hi, c)
  	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
  	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
  	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function npp_or_pow4(x::MC{N,Diff}, c::Z, y::Interval{Float64}) where {N, Z<:Integer}
  	if (x.Intv.hi < 0.0)
    	eps_min = x.Intv.hi
    	eps_max = x.Intv.lo
  	elseif (x.Intv.lo > 0.0)
    	eps_min = x.Intv.lo
    	eps_max = x.Intv.hi
  	else
    	eps_min = 0.0
    	eps_max = (abs(x.Intv.lo) >= abs(x.Intv.hi)) ? x.Intv.lo : x.Intv.hi
  	end
  	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
  	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
  	cc, dcc = cc_npp_or_pow4(midcc, x.Intv.lo, x.Intv.hi, c)
  	cv, dcv = cv_npp_or_pow4(midcv, x.Intv.lo, x.Intv.hi, c)
  	gcc1, gdcc1 = cc_npp_or_pow4(x.cv, x.Intv.lo, x.Intv.hi, c)
  	gcv1, gdcv1 = cv_npp_or_pow4(x.cv, x.Intv.lo, x.Intv.hi, c)
  	gcc2, gdcc2 = cc_npp_or_pow4(x.cc, x.Intv.lo, x.Intv.hi, c)
  	gcv2, gdcv2 = cv_npp_or_pow4(x.cc, x.Intv.lo, x.Intv.hi, c)
  	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
  	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
  	return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end

@inline function pos_odd(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}}
    midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi)
    midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo)
    cc, dcc = cc_powodd_ns(midcc, x.Intv.lo, x.Intv.hi, c)
    cv, dcv = cv_powodd_ns(midcv, x.Intv.lo, x.Intv.hi, c)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function pos_odd(x::MC{N,Diff}, c::Z, y::Interval{Float64}) where {N, Z<:Integer}
    midcc, cc_id = mid3(x.cv, x.cc, x.Intv.hi)
    midcv, cv_id = mid3(x.cv, x.cc, x.Intv.lo)
    cc, dcc = cc_powodd(midcc, x.Intv.lo, x.Intv.hi, c)
    cv, dcv = cv_powodd(midcv, x.Intv.lo, x.Intv.hi, c)
    gcc1, gdcc1 = cc_powodd(x.cv, x.Intv.lo, x.Intv.hi, c)
    gcv1, gdcv1 = cv_powodd(x.cv, x.Intv.lo, x.Intv.hi, c)
    gcc2, gdcc2 = cc_powodd(x.cc, x.Intv.lo, x.Intv.hi, c)
    gcv2, gdcv2 = cv_powodd(x.cc, x.Intv.lo, x.Intv.hi, c)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function cv_powodd_ns(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n))
    (0.0 <= xL) && (return x^n, n*x^(n - 1))

	root_f(x::Float64, xL::Float64, xU::Float64) = (x^n - xL^n) - (x-xL)*(n)*(x^(n-1));
	root_df(x::Float64, xL::Float64, xU::Float64) = n*x^(n-1) + xL*n*(n-1)*x - n*x^(n-1) - n*(n-1)*x^(n-1);
	inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU)
	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
	if x <= inflection
		return dline_seg(^, pow_deriv, x, xL, inflection, n)
	else
		return x^n, n*x^(n-1)
	end
end
@inline function cc_powodd_ns(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer}
    (xU <= 0.0) && (return x^n, n*x^(n - 1))
    (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, n))

	root_f(x::Float64, xL::Float64, xU::Float64) = (xU^n-x^n) - (xU-x)*(n)*(x^(n-1));
	root_df(x::Float64, xL::Float64, xU::Float64) = -n*x^(n-1) - xU*n*(n-1)*x + n*x^(n-1) + n*(n-1)*x^(n-1);
	inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU)
	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
	if x <= inflection
		return x^n, n*x^(n-1)
	else
		return dline_seg(^, pow_deriv, x, inflection, xU, n)
	end
end

# neg_powneg_odd computes the McComrick relaxation of x^c where x < 0.0 and c is odd
@inline function cv_neg_powneg_odd(x::Float64, xL::Float64, xU::Float64, c::Z) where Z<:Integer
	if (xU == xL)
		dcv = 0.0
        cv = x^c
    else
		xUc = xU^c
        dcv = (xUc- xL^c)/(xU - xL)
        cv = xUc + dcv*(x - xL)
    end
	cv, dcv
end
@inline function cc_neg_powneg_odd(x::Float64, xL::Float64, xU::Float64, c::Z) where Z<:Integer
	x^c, c*x^(c-1)
end
@inline function neg_powneg_odd(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}}
  	xL = x.Intv.lo
  	xU = x.Intv.hi
  	midcc, cc_id = mid3(x.cc, x.cv, xU)
  	midcv, cv_id = mid3(x.cc, x.cv, xL)
	cc, dcc = cc_neg_powneg_odd(midcc, xL, xU, c)
	cv, dcv = cv_neg_powneg_odd(midcv, xL, xU, c)
  	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
  	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function neg_powneg_odd(x::MC{N,Diff}, c::Z, y::Interval{Float64}) where {N, Z<:Integer}
  	xL = x.Intv.lo
  	xU = x.Intv.hi
  	midcc = mid3v(x.cv, x.cc, xU)
  	midcv = mid3v(x.cv, x.cc, xL)
  	cc, dcc = cc_negpowneg(midcc, xL, xU, c)
  	cv, dcv = cv_negpowneg(midcv, xL, xU, c)
  	gcc1, gdcc1 = cc_negpowneg(x.cv, xL, xU, c)
  	gcv1, gdcv1 = cv_negpowneg(x.cv, xL, xU, c)
  	gcc2, gdcc2 = cc_negpowneg(x.cc, xL, xU, c)
  	gcv2, gdcv2 = cv_negpowneg(x.cc, xL, xU, c)
  	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
  	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
	return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end

# neg_powneg_odd computes the McComrick relaxation of x^c where x < 0.0 and c is even
@inline function neg_powneg_even(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}}
	xL = x.Intv.lo
	xU = x.Intv.hi
  	midcc, cc_id = mid3(x.cc, x.cv, xU)
  	midcv, cv_id = mid3(x.cc, x.cv, xL)
  	cc, dcc = cc_negpowneg(midcc, xL, xU, c)
  	cv, dcv = cv_negpowneg(midcv, xL, xU, c)
  	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
  	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
  	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function neg_powneg_even(x::MC{N,Diff}, c::Z, y::Interval{Float64}) where {N, Z<:Integer}
	xL = x.Intv.lo
	xU = x.Intv.hi
  	midcc = mid3v(x.cv, x.cc, xU)
  	midcv = mid3v(x.cv, x.cc, xL)
  	cc, dcc = cc_negpowneg(midcc, xL, xU, c)
  	cv, dcv = cv_negpowneg(midcv, xL, xU, c)
  	gcc1, gdcc1 = cc_negpowneg(x.cv, xL, xU, c)
  	gcv1, gdcv1 = cv_negpowneg(x.cv, xL, xU, c)
  	gcc2, gdcc2 = cc_negpowneg(x.cc, xL, xU, c)
  	gcv2, gdcv2 = cv_negpowneg(x.cc, xL, xU, c)
  	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
  	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
  	return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end

@inline function pow_kernel(x::MC{N,T}, c::Z, y::Interval{Float64}) where {Z<:Integer, N, T<:RelaxTag}
	c == 0 && (return one(MC{N,T}))
	c == 1 && (return x)
	if c > 0
        c == 2 && (return sqr_kernel(x, y))
		isodd(c) && (return pos_odd(x, c, y))
		return npp_or_pow4(x, c, y)
    else
        if lo(x) < 0.0
        	isodd(c) && (return neg_powneg_odd(x, c, y))
        	return neg_powneg_even(x, c, y)
		end
    end
	return npp_or_pow4(x, c, y)
end

if ~(VERSION < v"1.1-")
	@inline function pow(x::MC{N,T}, c::Z) where {Z<:Integer, N, T<:RelaxTag}
		if (x.Intv.lo <= 0.0 <= x.Intv.hi) && (c < 0)
			return nan(MC{N,T})
		end
		(c < 0) && pow_kernel(x, c, inv(pow(x.Intv,-c)))
		return pow_kernel(x, c, pow(x.Intv,c))
	end
else
	@inline function pow(x::MC{N,T}, c::Z) where {Z<:Integer, N, T<:RelaxTag}
		if (x.Intv.lo <= 0.0 <= x.Intv.hi) && (c < 0)
			return nan(MC{N,T})
		end
		(c < 0) && pow_kernel(x, c, inv(x.Intv^(-c)))
		return pow_kernel(x, c, x.Intv^c)
	end
end
@inline (^)(x::MC, c::Z) where {Z <: Integer} = pow(x,c)

# Power of MC to float
@inline cv_flt_pow_1(x::Float64, xL::Float64, xU::Float64, n::Float64) = dline_seg(^, pow_deriv, x, xL, xU, n)
@inline cc_flt_pow_1(x::Float64, xL::Float64, xU::Float64, n::Float64) = x^n, n*x^(n-1)
@inline function flt_pow_1(x::MC{N,T}, c::Float64, y::Interval{Float64}) where {N, T<:Union{NS,MV}}
	midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi)
	midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo)
	cc, dcc = cc_flt_pow_1(midcc, x.Intv.lo, x.Intv.hi, c)
	cv, dcv = cv_flt_pow_1(midcv, x.Intv.lo, x.Intv.hi, c)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function flt_pow_1(x::MC{N,Diff}, c::Float64, y::Interval{Float64}) where N
	midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi)
	midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo)
	cc, dcc = cc_flt_pow_1(midcc, x.Intv.lo, x.Intv.hi, c)
	cv, dcv = cv_flt_pow_1(midcv, x.Intv.lo, x.Intv.hi, c)
    gcc1, gdcc1 = cc_flt_pow_1(x.cv ,x.Intv.lo, x.Intv.hi, c)
    gcv1, gdcv1 = cv_flt_pow_1(x.cv, x.Intv.lo, x.Intv.hi, c)
    gcc2, gdcc2 = cc_flt_pow_1(x.cc, x.Intv.lo, x.Intv.hi, c)
    gcv2, gdcv2 = cv_flt_pow_1(x.cc, x.Intv.lo, x.Intv.hi, c)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end

@inline function (^)(x::MC{N,T}, c::Float64, y::Interval{Float64}) where {N, T<:Union{NS,MV}}
	if (x.Intv.lo <= 0.0 <= x.Intv.hi) && (c < 0)
		return nan(MC{N,T})
	end
    isinteger(c) && (return pow_kernel(x, Int(c), y))
    ((x.Intv.lo >= 0) && (0.0 < c < 1.0)) && (return flt_pow_1(x, c, y))
	z = exp(c*log(x))
    return MC{N,T}(z.cv, z.cc, y, z.cv_grad, z.cc_grad, x.cnst)
end
@inline function (^)(x::MC{N,Diff}, c::Float64, y::Interval{Float64}) where N
	if (x.Intv.lo <= 0.0 <= x.Intv.hi) && (c < 0)
		return nan(MC{N, Diff})
	end
    isinteger(c) && (return pow_kernel(x, Int(c), y))
    ((x.Intv.lo >= 0) && (0.0 < c < 1.0)) && (return flt_pow_1(x, c, y))
	z = exp(c*log(x))
    return MC{N,Diff}(z.cv, z.cc, y, z.cv_grad, z.cc_grad, x.cnst)
end

@inline (^)(x::MC, c::Float32, y::Interval{Float64}) = (^)(x, Float64(c), y)
@inline (^)(x::MC, c::Float16, y::Interval{Float64}) = (^)(x, Float64(c), y)
@inline (^)(x::MC, c::Float64) = (^)(x, c, x.Intv^c)
@inline (^)(x::MC, c::Float32) = x^Float64(c)
@inline (^)(x::MC, c::Float16) = x^Float64(c)
@inline (^)(x::MC, c::MC) = exp(c*log(x))
@inline pow(x::MC, c::F) where {F <: AbstractFloat} = x^c

# Define powers to MC of floating point number
@inline function pow(b::Float64, x::MC{N,T}) where {N,T<:RelaxTag}
	(b <= 0.0) && (return nan(MC{N,T}))
	exp(x*log(b))
end
@inline ^(b::Float64, x::MC) = pow(b, x) # DONE (no kernel)

########### Defines inverse
@inline function cc_inv1(x::Float64, xL::Float64, xU::Float64)
	t = (xL*xU)
	cc = (xU + xL - x)/t
	dcc = -1.0/t
	return cc, dcc
end
@inline function cv_inv1(x::Float64, xL::Float64, xU::Float64)
	cv = 1.0/x
	dcv = -1.0/(x*x)
	return cv, dcv
end
@inline function inv1(x::MC{N,T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
	eps_min = x.Intv.hi
	eps_max = x.Intv.lo
  	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
  	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
  	cc, dcc = cc_inv1(midcc, x.Intv.lo, x.Intv.hi)
  	cv, dcv = cv_inv1(midcv, x.Intv.lo, x.Intv.hi)
  	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
  	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
  	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function inv1(x::MC{N,Diff}, y::Interval{Float64}) where N
	xL = x.Intv.lo
	xU = x.Intv.hi
	midcc = mid3v(x.cv, x.cc, xL)
	midcv = mid3v(x.cv, x.cc, xU)
	cc, dcc = cc_inv1(midcc, xL, xU )
	cv, dcv = cv_inv1(midcv, xL, xU)
	gcc1, gdcc1 = cc_inv1(x.cv, xL, xU)
	gcv1, gdcv1 = cv_inv1(x.cv, xL, xU)
	gcc2, gdcc2 = cc_inv1(x.cc, xL, xU)
	gcv2, gdcv2 = cv_inv1(x.cc, xL, xU)
	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
  	return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function inv_kernel(x::MC{N,T}, y::Interval{Float64}) where {N,T<:RelaxTag}
	if x.Intv.lo <= 0.0 <= x.Intv.hi
		return nan(MC{N,T})
	end
	if x.Intv.hi < 0.0
		x = pos_odd(x, -1, y)
  	else x.Intv.lo > 0.0
		x = inv1(x, y)
	end
	return x
end
@inline inv(x::MC{N,T}) where {N,T<:RelaxTag} = inv_kernel(x, inv(x.Intv))


text:
using Distributions
using LinearAlgebra
using Random

# Define the DGP
function simulate_sample(N0, N1)
    dy1 = Weibull(2, 1);
    dy2 = Weibull(2, 1.3);
    dy = MixtureModel([Weibull(2, 0.5), Weibull(2, 1.5)], [0.5, 0.5])
    X0 = [ones(N0) rand([0, 1], N0, 4) rand(N0)]
    X1 = kron([ones(2, 4) [1, 0] [0, 0]], ones(N1))
    λ0 = @. exp(-1.0 + 1.0 * X0[:, end - 1] + 1.0 * X0[:, end])
    y0 = λ0 .* rand(dy, N0)
    y1 = LinRange(0, 6, N1) |> collect |> x -> repeat(x, 2)
    c0 = rand(Exponential(10), N0)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

# Not run
Random.seed!(1);
N0, N1 = 500, 0;
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
println(mean(z0 .> c0))
println(maximum(z0))
println(mean(z0))


text:
atreplinit() do REPL
  schedule(@task begin
    sleep(0.1)
    try
      @eval using Revise
    catch
      @warn "Could not load Revise."
    end

#    try
#      @eval using Rebugger
#      # Activate Rebugger's key bindings
#      Rebugger.keybindings[:stepin] = "\e[17~"      # Add the keybinding F6 to step into a function.
#      Rebugger.keybindings[:stacktrace] = "\e[18~"  # Add the keybinding F7 to capture a stacktrace.
#      Rebugger.repl_init(REPL)
#    catch
#      @warn "Could not turn on Rebugger key bindings."
#    end

    try
      @eval using OhMyREPL
    catch
      @warn "Could not load OhMyREPL."
    end
  end)
end

max_stars_repo_name: mdavezac/black-garlic

max_issues_repo_name: mdavezac/black-garlic

max_forks_repo_name: mdavezac/black-garlic


text:
### defining diff. rules for Dual numbers ###
using Distributions, LinearAlgebra

struct Dual{T} <: Number
    f::Union{T, NTuple{2, Dual{T}}}
    g::Union{T, AbstractArray}
end

DualRealArray = Union{Dual, Real, AbstractArray}
DualandReal = Union{Dual, Real}

### Differentiation rules via overloading ###
import Base: +,/,*,-,^,adjoint, convert, promote_rule
import Base: +,/,*,-,^, convert, promote_rule
+(x::Dual, y::Dual) = Dual(x.f + y.f, x.g + y.g)
+(x::Dual, y::Real) = Dual(x.f + y, x.g)
+(y::Real, x::Dual) = Dual(x.f + y, x.g)
-(x::Dual, y::Dual) = Dual(x.f - y.f, x.g - y.g)
-(x::Dual) = Dual(-x.f, -x.g)
-(x::Dual, y::Real) = Dual(x.f -y, x.g)
-(y::Real, x::Dual) = Dual(y-x.f, -x.g)
*(x::Dual, y::Dual) = Dual(x.f*y.f, x.f*y.g + y.f*x.g)
*(x::Dual, y::Real) = Dual(x.f*y, x.g*y)
*(y::Real, x::Dual) = Dual(x.f*y, x.g*y)
/(x::Dual, y::Dual) = Dual(x.f/y.f, (y.f*x.g - x.f*y.g)/y.f^2)
/(y::Real, x::Dual) = Dual(y/x.f, (-y*x.g) / x.f^2)
/(x::Dual, y::Real) = Dual(x.f/y, x.g/y)
^(x::Dual, k::Real) = Dual(x.f^(k), (x.g * k) * x.f ^ (k-1))
^(x::Dual, k::Int) = Dual(x.f^(k), (k * x.g) * x.f ^ (k-1))
Base.exp(x::Dual) = Dual(exp(x.f), x.g * exp(x.f))
Base.sqrt(x::Dual) = Dual(sqrt(x.f), x.g / (2 * sqrt(x.f)))
Base.log(x::Dual) = Dual(log(x.f), x.g/x.f)

Base.sin(x::Dual) = Dual(sin(x.f), cos(x.f)*x.g)
Base.cos(x::Dual) = Dual(cos(x.f), -sin(x.f)*x.g)

Base.abs(x::Dual) = sqrt(abs2(x))
Base.abs2(x::Dual) = real(conj(x)*x)
Base.conj(x::Dual) = Dual(conj(x.f) , conj(x.g))
Base.real(x::Dual) = Dual(real(x.f) , real(x.g))

Distributions.cdf(d, x::Dual) = Dual(cdf(d, x.f), pdf(d, x.f) * x.g)
Base.adjoint(x::Dual) = Dual(adjoint(x.f), adjoint(x.g))
LinearAlgebra.Adjoint(x::Dual) = Dual(LinearAlgebra.Adjoint(x.f), LinearAlgebra.Adjoint(x.g))
LinearAlgebra.dot(x::Dual, y::Dual) = Dual(dot(x.f,y.f), x.f * y.g + y.f * x.g)
# Base.zero(x::Dual) = Dual(zero(x.f), zero(x.g))

Base.isless(z::Dual,w::Dual) = z.f < w.f
Base.isless(z::Real,w::Dual) = z < w.f
Base.isless(z::Dual,w::Real) = z.f < w

Base.zero(x::Dual) = zero(x.f)

Base.one(x::Dual) = one(x.f)
### conversion, promotion rules ###
convert(::Type{Dual}, x::Real) = Dual(x, one(x))
convert(::Type{Dual}, x::AbstractArray) = DualArray(x)
convert(::Type{Array}, x::Real) = [x]
Dual(x) = convert(Dual, x)
Dual(x::Dual) = Dual(x, one(x.f))
promote_rule(::Type{Dual}, ::Type{<:Number}) = Dual

function DualArray(x::AbstractArray)
    l = length(x)
    eye = I(l)
    collect = []
    for i in 1:l
        push!(collect, Dual(x[i], view(eye,i,:,)))
    end
    return collect
end

### Recursive functions for getting derivatives ###
function chain(x::DualandReal, n::Int)
    dualone = one(x)
    if n == 1
        return Dual(x, dualone)
    else
        return chain(Dual(x, dualone), n-1)
    end
end

function chain(x::AbstractArray, n::Int)
    if n == 1
        return DualArray(x)
    else
        return chain(DualArray(x), n-1)
    end
end

function dechain(x::DualandReal)
    if x isa Real
        return x
    else
        return dechain(x.g)
    end
end

function dechain(x::DualRealArray, dim::Int)
    if x isa Dual
        return dechain(getfield(x, :g), dim)
    elseif eltype(x) <: Real
        l = size(x, 1)
        ar = Int.(ones(dim).*l)
        x = reshape(x, tuple(ar...)...)
        return x
    elseif eltype(x) <: Dual
        return dechain(getfield.(x, :g), dim)
    elseif eltype(x) <: Vector
        x = cat(x...;dims = dim)
        return dechain(x, dim)
    else
        return dechain(getfield.(x, :g), dim)
    end
end

derivative(f::Function, x::DualandReal, n::Int) = dechain(f.(chain([x],n)), 1)[1]
derivative(f::Function, x::DualandReal) = dechain(f.(chain([x],1)), 1)[1]
gradient(f::Function, x::AbstractArray) = dechain(f(chain(x,1)), 1)
gradient(f::Function, x::AbstractArray, n::Int) = dechain(f(chain(x, n)), n)
hessian(f::Function, x::AbstractArray) = dechain(f(chain(x, 2)), 2)

## to do, add support for multi dimensional arrays (change dual array function)

max_stars_repo_name: kgeoffrey/HigherOrderDerivatives

max_issues_repo_name: kgeoffrey/HigherOrderDerivatives

max_forks_repo_name: kgeoffrey/HigherOrderDerivatives


text:
import Base.println, Base.flush

const DEFAULT_MAX_FILE_SIZE = 5000 * 1028

"""
    FileRoller <: IO

Is responsible for managing a rolling log file.

# Fields
* `prefix::AbstractString`: filename prefix for the log.
* `folder::AbstractString`: directory where the log should be written.
* `filepath::AbstractString`: the full filepath for the log
* `file::IO`: the current file IO handle
* `byteswritten::Int64`: keeps track of how many bytes have been written to the current file.
* `max_sz::Int`: the maximum number of bytes written to a file before rolling over to another.
"""
mutable struct FileRoller <: IO
    prefix::AbstractString
    folder::AbstractString
    filepath::AbstractString
    file::IO
    byteswritten::Int64
    max_sz::Int
end

"""
    getsuffix(::Integer) -> String

Formats the nth file suffix.
"""
function getsuffix(n::Integer)
    str = string(n)
    for i = 1:(4 - length(str))
        str = "0"*str
    end
    str
end

"""
    getfile(folder::AbstractString, prefix::AbstractString) -> String, IO

Grabs the next log file.
"""
function getfile(folder::AbstractString, prefix::AbstractString)
    i = 1
    getpath(i) = joinpath(folder, "$(prefix).$(getsuffix(i))")
    p = getpath(i)
    while isfile(p)
        p = getpath(i)
        i += 1
    end
    p, open(p, "a")
end

"""
    FileRoller(prefix; max_size=DEFAULT_MAX_FILE_SIZE)

Creates a rolling log file in the current working directory
with the specified prefix.
"""
FileRoller(prefix; max_sz=DEFAULT_MAX_FILE_SIZE) = FileRoller(prefix, pwd(); max_sz=max_sz)

"""
    FileRoller(prefix, dir; max_size=DEFAULT_MAX_FILE_SIZE)

Creates a rolling log file in the specified directory with the given prefix.
"""
function FileRoller(prefix, dir; max_sz=DEFAULT_MAX_FILE_SIZE)
    FileRoller(prefix, dir, (getfile(dir, prefix))..., 0, max_sz)
end

"""
    println(::FileRoller, ::AbstractString)

Writes the string to a file and creates a new file if
we've reached the max file size.
"""
function println(f::FileRoller, s::AbstractString)
    if f.byteswritten > f.max_sz
        gf = getfile(f.folder, f.prefix)
        f.filepath = gf[1]
        f.file = gf[2]
        f.byteswritten = 0
    end
    f.byteswritten += write(f.file, "$s\n")
end

"""
    flush(::FileRoller)

Flushes the current open file.
"""
flush(f::FileRoller) = flush(f.file)


text:
##############################################################################
##
## Functions related to the interface with PrettyTables.jl.
##
##############################################################################

# Default DataFrames highlighter for text backend.
#
# This highlighter changes the text color to gray in cells with `nothing`,
# `missing`, `#undef`, and types related to DataFrames.jl.
function _pretty_tables_highlighter_func(data, i::Integer, j::Integer)
    try
        cell = data[i, j]
        return ismissing(cell) ||
            cell === nothing ||
            cell isa Union{AbstractDataset}
    catch e
        if isa(e, UndefRefError)
            return true
        else
            rethrow(e)
        end
    end
end

const _PRETTY_TABLES_HIGHLIGHTER = Highlighter(_pretty_tables_highlighter_func,
                                               Crayon(foreground = :dark_gray))

# Default DataFrames formatter for text backend.
#
# This formatter changes how the following types are presented when rendering
# the data frame:
#     - missing;
#     - nothing;
#     - Cells with types related to DataFrames.jl.

function _pretty_tables_general_formatter(v, i::Integer, j::Integer)
    if typeof(v) <: Union{AbstractDataset}

        # Here, we must not use `print` or `show`. Otherwise, we will call
        # `_pretty_table` to render the current table leading to a stack
        # overflow.
        return sprint(summary,v)
    # elseif ismissing(v)
    #     return "missing"
    # elseif v === nothing
    #     return ""
    else
        return v
    end
end


text:
@testset "LM algorithm + bound constraints" begin

    create_bc(l,u;n::Int) = BoundConstraints(l*ones(n),u*ones(n))

    @testset "Rosenbrock sol ∈ I" begin
        using NLS_Solver: LevenbergMarquardt_BC

        nls = Rosenbrock()
        θ=Float64[-10;-10]
        bc=create_bc(-1,2,n=parameter_size(nls))
        result=LevenbergMarquardt_BC(nls, θ, bc)

        @test converged(result)
        @test solution(result) ≈ Float64[1;1]
        @test iteration_count(result) == 14
    end

    @testset "Rosenbrock sol ∈ ∂I" begin
        using NLS_Solver: LevenbergMarquardt

        nls = Rosenbrock()
        θ=Float64[-10;-10]
        bc=create_bc(-1,1,n=parameter_size(nls))
        result=LevenbergMarquardt_BC(nls, θ, bc)

        @test converged(result)
        @test solution(result) ≈ Float64[1;1]
        @test iteration_count(result) == 14
    end

      @testset "Rosenbrock sol ∉ I" begin
        using NLS_Solver: LevenbergMarquardt

        nls = Rosenbrock()
        θ=Float64[-10;-10]
        bc=create_bc(-10,-2,n=parameter_size(nls))
        result=LevenbergMarquardt_BC(nls, θ, bc)

        @test converged(result)
        @test solution(result) ≈ Float64[-2;-2] # OK (MMA)
        @test iteration_count(result) == 4
    end

    # ================================================================

    @testset "Rosenbrock, use solve interface" begin
        nls = Rosenbrock()
        θ=Float64[-10;-10]
        bc=create_bc(-1,2,n=parameter_size(nls))

        conf = LevenbergMarquardt_BC_Conf()
        result=solve(nls, θ, bc, conf)

        @test converged(result)
        @test solution(result) ≈ Float64[1;1]
        @test iteration_count(result) == 14
    end

end


text:
module MatterPower
using OrdinaryDiffEq
using DoubleExponentialFormulas
using Roots
export t_nowiggle
export setup_growth
export sigma2, dsigma2dR
export sigma2gaus, dsigma2gausdR, d2sigma2gausdR2
export setup_halofit, halofit
include("eisensteinhu.jl")
include("growth.jl")
include("sigma2.jl")
include("halofit.jl")
end # module


text:
###############
###############
###############

T = 100                                 # Number of periods in the data
P = 2                                   # Number of lags in the factor equation
L = 1                                   # Number of AR lags in the observation equation
N = 2                                   # Total number of series
K = 1                                   # Total number of factors (only the global factor)

num_obs = T
γ1 = 1.0
γ2 = 1.0
ϕ1 = 0.5
ϕ2 = 0.1
ψ1 = 0.5
ψ2 = 0.25
σ2_1 = 0.1
σ2_2 = 0.1
H = [γ1 1.0 0.0 0.0; γ2 0.0 1.0 0.0]
A = zeros(2, 4)
F = [ϕ1 0.0 0.0 ϕ2; 0.0 ψ1 0.0 0.0; 0.0 0.0 ψ2 0.0; 1.0 0.0 0.0 0.0]
μ = [0.0, 0.0, 0.0, 0.0]
R = zeros(2, 2)
Q = zeros(4, 4)
Q[1, 1] = 1.0
Q[2, 2] = σ2_1
Q[3, 3] = σ2_2
Z = zeros(4, 4)

# Gather all SS parameters
ssmodel = SSModel(H, A, F, μ, R, Q, Z)

# Simulate common components model in state space form
data_y, data_z, data_β = simulateSSModel(num_obs, ssmodel)

priors = DFMStruct(P, L, 1000, 50)             # Set model priors

###############
###############
###############

results = OW1LevelEstimator(data_y, priors)

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, :]))
    push!(quant33, quantile(results.F[i, :], 0.33))
    push!(quant66, quantile(results.F[i, :], 0.66))
    push!(medians, median(results.F[i, :]))
end

plot(results.means.F[:, j])
plot(medians)
plot!(quant33)
plot!(quant66)
plot!(results.means.F[:, j] - stds)
plot!(results.means.F[:, j] + stds)


text:
using Test
using POMDPModels

let
    rng = MersenneTwister(41)
    p = LightDark1D()
    @test discount(p) == 0.9
    s0 = LightDark1DState(0,0)
    s0, _, r = gen(DDNOut(:sp, :o, :r), p, s0, +1, rng)
    @test s0.y == 1.0
    @test r == 0
    s1, _, r = gen(DDNOut(:sp, :o, :r), p, s0, 0, rng)
    @test s1.status != 0
    @test r == -10.0
    s2 = LightDark1DState(0, 5)
    obs = gen(DDNNode(:o), p, nothing, nothing, s2, rng)
    @test abs(obs-6.0) <= 1.1

    sv = convert_s(Array{Float64}, s2, p)
    @test sv == [0.0, 5.0]
    s = convert_s(LightDark1DState, sv, p)
    @test s == s2

    ov = convert_o(Array{Float64}, obs, p)
    @test ov == [obs]
    o = convert_o(Float64, ov, p)
    @test o == obs
end