text:
"""
    get_cluster(g::AbstractGraph, node::Integer)

Determines the cluster in `g` which `node` is a member of

Arguments
* `g`   : An instance of type AbstractGraph
* `node`: Node in `keys(g.cluster_ids)`
Returns
* Set of nodes representing the cluster which `node` is a member of
"""
function get_cluster(g::AbstractGraph, node::Integer)
    return g.clusters[g.cluster_ids[node]]
end

"""
    get_largest_cluster_size(g::AbstractGraph)

Determines the size of the largest cluster in `g`

Arguments
* `g`: An instance of type AbstractGraph
Returns
* Integer representing the number of nodes in the largest cluster in `g`
"""
function get_largest_cluster_size(g::AbstractGraph)
    return maximum(keys(g.cluster_sizes))
end

"""
    get_avg_cluster_size(g::AbstractGraph)

Determines the average cluster size in `g`

Arguments
* `g`: An instance of type AbstractGraph
Returns
* Float64 representing the average cluster size in `g`
"""
function get_avg_cluster_size(g::AbstractGraph)
    return sum([key * value for (key, value) in g.cluster_sizes]) / length(g.clusters)
end

"""
    get_largest_clusters(g::AbstractGraph, n_clusters::Integer)

Determines the `n_clusters` largest clusters in 'g'

Arguments
* `g`         : An instance of type AbstractGraph
* `n_clusters`: The number of largest clusters to return
Returns
* Sorted (descending) array of the `n_clusters` largest clusters in `g`
"""
function get_largest_clusters(g::AbstractGraph, n_clusters::Integer)
    return sort!(collect(values(g.clusters)), by = length, rev = true)[1:n_clusters]
end

"""
    update_clusters!(g::AbstractGraph, edge::Tuple)

Updates `g` with the newly merged cluster and the largest cluster size
Arguments
* `g`   : An instance of type AbstractGraph
* `edge`: Edge added to `g` at step `g.t`
Returns
* None, updates `g` in-place
"""
function update_clusters!(g::AbstractGraph, edge::NTuple{2,Integer})
    if length(g.clusters[g.cluster_ids[edge[1]]]) >
       length(g.clusters[g.cluster_ids[edge[2]]])
        larger_cluster_id = g.cluster_ids[edge[1]]
        smaller_cluster_id = g.cluster_ids[edge[2]]
    else
        larger_cluster_id = g.cluster_ids[edge[2]]
        smaller_cluster_id = g.cluster_ids[edge[1]]
    end

    if g.cluster_ids[edge[1]] ≠ g.cluster_ids[edge[2]]
        update_cluster_sizes!(g, larger_cluster_id, smaller_cluster_id)
        update_cluster_ids!(g, larger_cluster_id, smaller_cluster_id)
        merge_clusters!(g, larger_cluster_id, smaller_cluster_id)
        update_observables!(g, larger_cluster_id, smaller_cluster_id)
    else
        update_observables!(g, larger_cluster_id, smaller_cluster_id)
    end
end

"""
    update_cluster_sizes!(g::AbstractGraph, larger_cluster_id::Integer, smaller_cluster_id::Integer)

Updates the cluster size distribution dictionary

Arguments
* `g`                 : An instance of type AbstractGraph
* `larger_cluster_id` : Cluster ID of the larger cluster
* `smaller_cluster_id`: Cluster ID of the smaller cluster
Returns
* None, updates `g` in-place
"""
function update_cluster_sizes!(
    g::AbstractGraph,
    larger_cluster_id::Integer,
    smaller_cluster_id::Integer,
)
    if g.cluster_sizes[length(g.clusters[smaller_cluster_id])] ≠ 1
        g.cluster_sizes[length(g.clusters[smaller_cluster_id])] -= 1
    else
        delete!(g.cluster_sizes, length(g.clusters[smaller_cluster_id]))
    end

    if g.cluster_sizes[length(g.clusters[larger_cluster_id])] ≠ 1
        g.cluster_sizes[length(g.clusters[larger_cluster_id])] -= 1
    else
        delete!(g.cluster_sizes, length(g.clusters[larger_cluster_id]))
    end

    if haskey(
        g.cluster_sizes,
        (length(g.clusters[larger_cluster_id]) + length(g.clusters[smaller_cluster_id])),
    )
        g.cluster_sizes[(length(g.clusters[larger_cluster_id])+length(g.clusters[smaller_cluster_id]))] +=
            1
    else
        g.cluster_sizes[(length(g.clusters[larger_cluster_id])+length(g.clusters[smaller_cluster_id]))] =
            1
    end
end

"""
    update_cluster_ids!(g::AbstractGraph, larger_cluster_id::Integer, smaller_cluster_id::Integer)

Updates the cluster IDs of the nodes in the smaller cluster to that of the larger cluster it is being merged into

Arguments
* `g`                 : An instance of type AbstractGraph
* `larger_cluster_id` : Cluster ID of the larger cluster
* `smaller_cluster_id`: Cluster ID of the smaller cluster
Returns
* None, updates `g` in-place
"""
function update_cluster_ids!(
    g::AbstractGraph,
    larger_cluster_id::Integer,
    smaller_cluster_id::Integer,
)
    for node in g.clusters[smaller_cluster_id]
        g.cluster_ids[node] = larger_cluster_id
    end
end

"""
    merge_clusters!(g::AbstractGraph, larger_cluster_id::Integer, smaller_cluster_id::Integer)

Mergers the smaller cluster into the larger cluster in-place

Arguments
* `g`                 : An instance of type AbstractGraph
* `larger_cluster_id` : Cluster ID of the larger cluster
* `smaller_cluster_id`: Cluster ID of the smaller cluster
Returns
* None, updates `g` in-place
"""
function merge_clusters!(
    g::AbstractGraph,
    larger_cluster_id::Integer,
    smaller_cluster_id::Integer,
)
    union!(g.clusters[larger_cluster_id], g.clusters[smaller_cluster_id])
    delete!(g.clusters, smaller_cluster_id)
end

"""
    update_observables!(g::AbstractGraph, larger_cluster_id::Integer, smaller_cluster_id::Integer)

Updates the largest cluster size, average cluster size, and cluster heterogeneity

Arguments
* `g`                 : An instance of type AbstractGraph
* `larger_cluster_id` : Cluster ID of the larger cluster
* `smaller_cluster_id`: Cluster ID of the smaller cluster
Returns
* None, updates `g` in-place
"""
function update_observables!(
    g::AbstractGraph,
    larger_cluster_id::Integer,
    smaller_cluster_id::Integer,
)
    push!(
        g.observables.largest_cluster_size,
        maximum((
            g.observables.largest_cluster_size[g.t],
            length(g.clusters[larger_cluster_id]),
        )),
    )
    push!(g.observables.avg_cluster_size, get_avg_cluster_size(g))
    push!(g.observables.heterogeneity, length(g.cluster_sizes))
end


text:
## The standard univariate linear Gaussian state-space model example

module LinearGaussian

using SequentialMonteCarlo
using RNGPool
import SMCExamples.Particles.Float64Particle

struct LGTheta
  A::Float64
  Q::Float64
  C::Float64
  R::Float64
  x0::Float64
  v0::Float64
end

function makeLGModel(theta::LGTheta, ys::Vector{Float64})
  n::Int64 = length(ys)
  sqrtQ::Float64 = sqrt(theta.Q)
  invRover2::Float64 = 0.5/theta.R
  sqrtv0::Float64 = sqrt(theta.v0)
  logncG::Float64 = -0.5 * log(2 * π * theta.R)
  @inline function lG(p::Int64, particle::Float64Particle, ::Nothing)
    @inbounds v::Float64 = theta.C*particle.x - ys[p]
    return logncG - v * invRover2 * v
  end
  @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
    particle::Float64Particle, ::Nothing)
    if p == 1
      newParticle.x = theta.x0 + sqrtv0*randn(rng)
    else
      newParticle.x = theta.A*particle.x + sqrtQ*randn(rng)
    end
  end
  return SMCModel(M!, lG, n, Float64Particle, Nothing)
end

function makelM(theta::LGTheta)
  tQ::Float64 = theta.Q
  negInvQover2::Float64 = -0.5/tQ
  tA::Float64 = theta.A
  @inline function lM(::Int64, particle::Float64Particle,
    newParticle::Float64Particle, ::Nothing)
    x::Float64 = particle.x
    y::Float64 = newParticle.x
    t::Float64 = tA*x - y
    return negInvQover2 * t * t
  end
  return lM
end

function simulateLGModel(theta::LGTheta, n::Int64)
  model = makeLGModel(theta, Vector{Float64}(undef, 0))
  ys = Vector{Float64}(undef, n)
  xParticle = Float64Particle()
  rng = getRNG()
  for p in 1:n
    model.M!(xParticle, rng, p, xParticle, nothing)
    ys[p] = theta.C*xParticle.x + sqrt(theta.R)*randn(rng)
  end
  return ys
end

function defaultLGModel(n::Int64)
  theta = LGTheta(0.9, 0.6, 1.2, 1.5, 1.0, 2.0)
  ys = simulateLGModel(theta, n)

  ko = kalman(theta, ys)

  lgModel = makeLGModel(theta, ys)

  return lgModel, theta, ys, ko
end

include("lgKalman.jl")
include("lgLOPModel.jl")
include("lgAPFModel.jl")

end

ext: jl


text:
using Stopping

import Stopping: _init_max_counters, fill_in!, _resources_check!, _unbounded_problem_check!, _optimality_check
import Stopping: start!, stop!, update_and_start!, update_and_stop!, reinit!, status

"""
Type: MPCCStopping (specialization of GenericStopping)
Methods: start!, stop!, update_and_start!, update_and_stop!, fill_in!, reinit!, status

Stopping structure for non-linear programming problems using NLPModels.
    Input :
       - pb         : an AbstractMPCCModel
       - state      : The information relative to the problem, see GenericState
       - (opt) meta : Metadata relative to stopping criterion, see *StoppingMeta*.
       - (opt) main_stp : Stopping of the main loop in case we consider a Stopping
                          of a subproblem.
                          If not a subproblem, then nothing.
       - (opt) listofstates : ListStates designed to store the history of States.
       - (opt) user_specific_struct : Contains any structure designed by the user.

`MPCCStopping(:: AbstractMPCCModel, :: AbstractState; meta :: AbstractStoppingMeta = StoppingMeta(), main_stp :: Union{AbstractStopping, Nothing} = nothing, list :: Union{ListStates, Nothing} = nothing, user_specific_struct :: Any = nothing, kwargs...)`

 Note:
 * optimality_check : takes two inputs (AbstractMPCCModel, MPCCAtX)
 and returns a Float64 to be compared at 0.
 * designed for MPCCAtX State. Constructor checks that the State has the
 required entries.

 Warning:
 * optimality_check does not necessarily fill in the State.
 """
mutable struct MPCCStopping <: AbstractStopping

    # problem
    pb :: AbstractMPCCModel

    # Common parameters
    meta      :: AbstractStoppingMeta

    # current state of the problem
    current_state :: AbstractState

    # Stopping of the main problem, or nothing
    main_stp :: Union{AbstractStopping, Nothing}

    listofstates :: Union{ListStates, Nothing}

    # User-specific structure
    user_specific_struct :: Any

    function MPCCStopping(pb            :: AbstractMPCCModel,
                         current_state  :: AbstractState;
                         meta           :: AbstractStoppingMeta = StoppingMeta(;max_cntrs = _init_max_counters_mpcc(), optimality_check = MStat),
                         main_stp       :: Union{AbstractStopping, Nothing} = nothing,
                         listofstates   :: Union{ListStates, Nothing} = nothing,
                         user_specific_struct :: Any = nothing,
                         kwargs...)

        if !(isempty(kwargs))
           meta = StoppingMeta(;max_cntrs = _init_max_counters_mpcc(), optimality_check = MStat, kwargs...)
        end

        #current_state is an AbstractState with requirements
        try
            current_state.evals
            current_state.fx, current_state.gx, current_state.Hx
            #if there are bounds:
            current_state.mu
            if pb.meta.ncon > 0 #if there are constraints
               current_state.Jx, current_state.cx, current_state.lambda
            end
        catch
            throw("error: missing entries in the given current_state")
        end

        return new(pb, meta, current_state, main_stp, listofstates, user_specific_struct)
    end

end

"""
MPCCStopping(pb): additional default constructor
The function creates a Stopping where the State is by default and the
optimality function is the function KKT().

key arguments are forwarded to the classical constructor.
"""
function MPCCStopping(pb :: AbstractMPCCModel; kwargs...)

 #Create a default MPCCAtX
 nlp_at_x = MPCCAtX(pb.meta.x0, pb.meta.y0)
 admissible = (x,y; kwargs...) -> MStat(x,y; kwargs...)

 return MPCCStopping(pb, nlp_at_x, optimality_check = admissible; kwargs...)
end

"""
_init_max_counters_mpcc(): initialize the maximum number of evaluations on each of
                        the functions present in the MPCCCounters.
"""
function _init_max_counters_mpcc(; obj     :: Int64 = 40000,
                                   grad    :: Int64 = 40000,
                                   cons    :: Int64 = 40000,
                                   consG   :: Int64 = 40000,
                                   consH   :: Int64 = 40000,
                                   jcon    :: Int64 = 40000,
                                   jgrad   :: Int64 = 40000,
                                   jac     :: Int64 = 40000,
                                   jacG    :: Int64 = 40000,
                                   jacH    :: Int64 = 40000,
                                   jprod   :: Int64 = 40000,
                                   jGprod  :: Int64 = 40000,
                                   jHprod  :: Int64 = 40000,
                                   jtprod  :: Int64 = 40000,
                                   jGtprod :: Int64 = 40000,
                                   jHtprod :: Int64 = 40000,
                                   hess    :: Int64 = 40000,
                                   hessG   :: Int64 = 40000,
                                   hessH   :: Int64 = 40000,
                                   hprod   :: Int64 = 40000,
                                   hGprod  :: Int64 = 40000,
                                   hHprod  :: Int64 = 40000,
                                   jhprod  :: Int64 = 40000,
                                   sum     :: Int64 = 40000*11)

  cntrs = Dict([(:neval_obj,       obj), (:neval_grad,   grad),
                (:neval_cons,     cons), (:neval_jcon,   jcon),
                (:neval_consG,   consG), (:neval_consH,   consH),
                (:neval_jgrad,   jgrad), (:neval_jac,    jac),
                (:neval_jacG,     jacG), (:neval_jacH,     jacH),
                (:neval_jprod,   jprod), (:neval_jtprod, jtprod),
                (:neval_jGprod, jGprod), (:neval_jGtprod, jGtprod),
                (:neval_jHprod, jHprod), (:neval_jHtprod, jHtprod),
                (:neval_hess,     hess), (:neval_hprod,  hprod),
                (:neval_hessG,   hessG), (:neval_hessH,   hessH),
                (:neval_hGprod,  hGprod),(:neval_hHprod,  hHprod),
                (:neval_jhprod, jhprod), (:neval_sum,    sum)])

 return cntrs
end

"""
fill_in!: a function that fill in the required values in the State

kwargs are forwarded to _compute_mutliplier function
"""
function fill_in!(stp     :: MPCCStopping,
                  x       :: Iterate;
                  fx      :: Iterate    = nothing,
                  gx      :: Iterate    = nothing,
                  Hx      :: Iterate    = nothing,
                  cx      :: Iterate    = nothing,
                  Jx      :: Iterate    = nothing,
                  lambda  :: Iterate    = nothing,
                  mu      :: Iterate    = nothing,
                  lambdaG :: Iterate    = nothing,
                  lambdaH :: Iterate    = nothing,
                  cGx     :: Iterate    = nothing,
                  JGx     :: MatrixType = nothing,
                  cHx     :: Iterate    = nothing,
                  JHx     :: MatrixType = nothing,
                  matrix_info :: Bool   = true,
                  kwargs...)

 gfx = fx == nothing  ? obj(stp.pb, x)   : fx
 ggx = gx == nothing  ? grad(stp.pb, x)  : gx

 if Hx == nothing && matrix_info
   gHx = hess(stp.pb, x)
 else
   gHx = Hx
 end

 if stp.pb.meta.ncon > 0
     gJx = Jx == nothing ? jac_nl(stp.pb, x)  : Jx
     gcx = cx == nothing ? cons_nl(stp.pb, x) : cx
 else
     gJx = Jx
     gcx = cx
 end

 gcGx, gcHx, gJGx, gJHx = cGx, cHx, JGx, JHx
 if stp.pb.meta.ncc > 0
     gcGx = cGx == nothing ? consG(stp.pb, x) : cGx
     gcHx = cHx == nothing ? consH(stp.pb, x) : cHx
     gJGx = JGx == nothing ? jacG(stp.pb, x)  : JGx
     gJHx = JHx == nothing ? jacH(stp.pb, x)  : JHx
 end

 #update the Lagrange multiplier if one of the 2 is asked
 #if (stp.pb.meta.ncc > 0 && stp.pb.meta.ncon > 0 && has_bounds(stp.pb)) && (lambdaG == nothing || lambdaH == nothing ||lambda == nothing || mu == nothing)
 if (stp.pb.meta.ncc > 0) && (lambdaG == nothing || lambdaH == nothing ||lambda == nothing || mu == nothing)
  lb, lc, lG, lH = _compute_mutliplier(stp.pb, x, ggx, gcx, gJx, gcGx, gJGx, gcHx, gJHx, kwargs...)
 #elseif (stp.pb.meta.ncon > 0 && has_bounds(stp.pb)) && (lambda == nothing || mu == nothing)
 elseif (stp.pb.meta.ncon > 0 || has_bounds(stp.pb)) && (lambda == nothing || mu == nothing)
  lb, lc = _compute_mutliplier(stp.pb.mp, x, ggx, gcx, gJx, kwargs...)
  lG, lH = lambdaG, lambdaH
 elseif  stp.pb.meta.ncon == 0 && !has_bounds(stp.pb) && lambda == nothing
  lb, lc = mu, stp.current_state.lambda
  lG, lH = lambdaG, lambdaH
 else
  lb, lc = mu, lambda
  lG, lH = lambdaG, lambdaH
 end

 return update!(stp.current_state, x=x, fx = gfx,    gx = ggx, Hx = gHx,
                                        cx = gcx,    Jx = gJx, mu = lb,
                                        lambda = lc, lambdaG = lG, lambdaH = lH,
                                        cGx = gcGx, cHx = gcHx,
                                        JGx = gJGx, JHx = gJHx)
end

"""
_resources_check!: check if the optimization algorithm has exhausted the resources.
                   This is the MPCC specialized version that takes into account
                   the evaluation of the functions following the sum_counters
                   structure from MPCC.

Note: function uses counters in stp.pb, and update the counters in the state
"""
function _resources_check!(stp    :: MPCCStopping,
                           x      :: Iterate)

  cntrs = stp.pb.counters
  update!(stp.current_state, evals = cntrs)

  max_cntrs = stp.meta.max_cntrs

  # check all the entries in the counter
  max_f = false
  for f in fieldnames(MPCCCounters)
      max_f = max_f || (getfield(cntrs, f) > max_cntrs[f])
  end

 # Maximum number of function and derivative(s) computation
 max_evals = sum_counters(stp.pb) > max_cntrs[:neval_sum]

 # global user limit diagnostic
 stp.meta.resources = max_evals || max_f

 return stp
end

"""
_unbounded_problem_check!: This is the NLP specialized version that takes into account
                   that the problem might be unbounded if the objective or the
                   constraint function are unbounded.

Note: * evaluate the objective function if state.fx is void.
      * evaluate the constraint function if state.cx is void.
"""
function _unbounded_problem_check!(stp  :: MPCCStopping,
                                   x    :: Iterate)

 if stp.current_state.fx == nothing
	 stp.current_state.fx = obj(stp.pb, x)
 end
 f_too_large = norm(stp.current_state.fx) >= stp.meta.unbounded_threshold

 c_too_large = false
 if stp.pb.meta.ncon != 0 #if the problems has constraints, check |c(x)|
  if stp.current_state.cx == nothing
   stp.current_state.cx = cons(stp.pb, x)
  end
  c_too_large = norm(stp.current_state.cx) >= abs(stp.meta.unbounded_threshold)
 end

 stp.meta.unbounded_pb = f_too_large || c_too_large

 return stp
end

################################################################################
# Nonlinear problems admissibility functions
# Available: unconstrained_check(...), optim_check_bounded(...), KKT
################################################################################
include("mpcc_admissible_functions.jl")

################################################################################
# Functions computing Lagrange multipliers of a nonlinear problem
# Available: _compute_mutliplier(...)
################################################################################
include("nlp_compute_multiplier.jl")


text:
using Random
export customop,
xavier_init,
load_op_and_grad,
load_op,
use_gpu,
install,
load_system_op,
register,
debug,
doctor,
nnuq,
compile,
list_physical_devices,
MCMCSimple,
simulate,
diagnose,
get_placement,
timestamp,
load_library,
sleep_for,
get_library_symbols

"""
    xavier_init(size, dtype=Float64)

Returns a matrix of size `size` and its values are from Xavier initialization.
"""
function xavier_init(size, dtype=Float64)
    in_dim = size[1]
    xavier_stddev = 1. / sqrt(in_dim / 2.)
    return randn(dtype, size...)*xavier_stddev
end

############### custom operators ##################
"""
    cmake(DIR::String=".."; CMAKE_ARGS::Union{Array{String}, String} = "")

The built-in Cmake command for building C/C++ libraries. If extra Cmake arguments are needed, please specify it through `CMAKE_ARGS`.

# Example
```
ADCME.cmake(CMAKE_ARGS=["SHARED=YES", "STAITC=NO"])
```

The executed command might be:
```
/home/darve/kailaix/.julia/adcme/bin/cmake -G Ninja -DCMAKE_MAKE_PROGRAM=/home/darve/kailaix/.julia/adcme/bin/ninja -DJULIA=/home/darve/kailaix/julia-1.3.1/bin/julia -DCMAKE_C_COMPILER=/home/darve/kailaix/.julia/adcme/bin/x86_64-conda_cos6-linux-gnu-gcc -DCMAKE_CXX_COMPILER=/home/darve/kailaix/.julia/adcme/bin/x86_64-conda_cos6-linux-gnu-g++ SHARED=YES STATIC=NO ..
```
"""
function cmake(DIR::String=".."; CMAKE_ARGS::Union{Array{String}, String} = "")
    DIR = abspath(DIR)
    ENV_ = copy(ENV)
    LD_PATH = Sys.iswindows() ? "PATH" : "LD_LIBRARY_PATH"
    if haskey(ENV_, LD_PATH)
        ENV_[LD_PATH] = ENV[LD_PATH]*":$LIBDIR"
    else
        ENV_[LD_PATH] = LIBDIR
    end

    if has_mpi(false)
        ENV_["MPI_INCLUDE_PATH"], ENV_["MPI_C_LIBRARIES"] = get_mpi()
    end

    if Sys.iswindows()
        if !haskey(ENV_, "VS150COMNTOOLS")
            # @warn "VS150COMNTOOLS is not set, default to /c/Program Files (x86)/Microsoft Visual Studio/2017/Community/Common7/Tools" maxlog=1
            ENV_["VS150COMNTOOLS"] = "/c/Program Files (x86)/Microsoft Visual Studio/2017/Community/Common7/Tools"
        end
        # @info "Do remember to add ADD_DEFINITIONS(-DNOMINMAX) to your CMakeLists.txt" maxlog=1
        run(setenv(`$CMAKE -G"Visual Studio 15" -DJULIA="$(joinpath(Sys.BINDIR, "julia"))" -A x64 $CMAKE_ARGS $DIR`, ENV_)) # very important, x64
    else
        run(setenv(`$CMAKE -G Ninja -DCMAKE_MAKE_PROGRAM=$NINJA -DJULIA="$(joinpath(Sys.BINDIR, "julia"))" -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX $CMAKE_ARGS $DIR`, ENV_))
    end
end

function make()
    ENV_ = copy(ENV)
    LD_PATH = Sys.iswindows() ? "PATH" : "LD_LIBRARY_PATH"
    if haskey(ENV_, LD_PATH)
        ENV_[LD_PATH] = ENV[LD_PATH]*":$LIBDIR"
    else
        ENV_[LD_PATH] = LIBDIR
    end
    if Sys.iswindows()
        sln_file = filter(x->endswith(x, ".sln"), readdir())
        if length(sln_file)==0
            error("No .sln file found. Did you run `ADCME.cmake()`?")
        elseif length(sln_file)>1
            error("More than 1 .sln file found. Check your program.")
        else
            sln_file = sln_file[1]
        end
        run(`cmd /c $CMAKE  --build . -j --target ALL_BUILD --config Release`)
    else
        if haskey(ENV, "TRAVIS_BRANCH")
            run(setenv(`$NINJA -j1`, ENV_))
        else
            run(setenv(`$NINJA -j20`, ENV_))
        end
    end
end

"""
    make_library(Libdir::String)

Make shared library in `Libdir`. The structure of the source codes files are

```
- Libdir
  - *.cpp
  - *.h
  - CMakeLists
  - build (Optional)
```
"""
function make_library(Libdir::String)
   if !isdir(Libdir)
    error("$Libdir is not a valid directory.")
   end
   PWD = pwd()
   cd(Libdir)
   if !isdir("build")
    mkdir("build")
   end
   cd("build")
   if !isfile("Makefile")
    ADCME.cmake()
   end
   ADCME.make()
   cd(PWD)
end

load_op_dict = Dict{Tuple{String, String}, PyObject}()
load_op_grad_dict = Dict{Tuple{String, String}, PyObject}()

@doc """
    load_op(oplibpath::Union{PyObject, String}, opname::String; verbose::Union{Missing, Bool} = missing)

Loads the operator `opname` from library `oplibpath`.
"""
function load_op(oplibpath::Union{PyObject, String}, opname::String; verbose::Union{Missing, Bool} = missing)
    verbose = coalesce(verbose, options.customop.verbose)
    oplibpath = get_library(oplibpath)
    if haskey(load_op_dict, (oplibpath,opname))
        return load_op_dict[(oplibpath,opname)]
    end
    s = getproperty(load_library(oplibpath), opname)
    load_op_dict[(oplibpath,opname)] = s
    verbose && printstyled("Load library operator: $oplibpath ==> $opname\n", color=:green)
    return s
end

@doc """
    load_op_and_grad(oplibpath::Union{PyObject, String}, opname::String; multiple::Bool=false)

Loads the operator `opname` from library `oplibpath`; gradients are also imported.
If `multiple` is true, the operator is assumed to have multiple outputs.
"""
function load_op_and_grad(oplibpath::Union{PyObject, String}, opname::String;
        multiple::Bool=false, verbose::Union{Missing, Bool} = missing)
    verbose = coalesce(verbose, options.customop.verbose)
    is_system_op = oplibpath == LIBADCME

    if isa(oplibpath, String)
        if Sys.iswindows()
            a, b = splitdir(oplibpath)
            if length(b) >=3 && b[1:3]=="lib"
                b = b[4:end]
            end
            oplibpath = joinpath(a, b)
        end
        if splitext(oplibpath)[2]==""
            oplibpath = oplibpath * (Sys.islinux() ?
                            ".so" : Sys.isapple() ? ".dylib" : ".dll")
        end
        oplibpath = abspath(oplibpath)
        if haskey(load_op_grad_dict, (oplibpath,opname))
            return load_op_grad_dict[(oplibpath,opname)]
        end
        if !isfile(oplibpath)
            error("Library $oplibpath does not exist.")
        end
    end

    opname_grad = opname*"_grad"
    fn_name = opname*randstring(8)
    try
        if is_system_op
py"""
import tensorflow as tf
lib$$fn_name = $libadcme
"""
        elseif isa(oplibpath, String)
py"""
import tensorflow as tf
lib$$fn_name = tf.load_op_library($oplibpath)
"""
        else
py"""
import tensorflow as tf
lib$$fn_name = $oplibpath
"""
        end
if !multiple
py"""
@tf.custom_gradient
def $$fn_name(*args):
    u = lib$$fn_name.$$opname(*args)
    def grad(dy):
        return lib$$fn_name.$$opname_grad(dy, u, *args)
    return u, grad
"""
else
py"""
@tf.custom_gradient
def $$fn_name(*args):
    u = lib$$fn_name.$$opname(*args)
    def grad(*dy):
        dy = [y for y in dy if y is not None and y.dtype in [tf.float64, tf.float32]] # only float64 and float32 can backpropagate gradients
        return lib$$fn_name.$$opname_grad(*dy, *u, *args)
    return u, grad
"""
end
catch(e)
    printstyled("Failed load $oplibpath or its symbols. Error Message from the TensorFlow backend\n$(string(e))\n", color=:red)
    Libdl.dlopen(oplibpath)
end
        s = py"$$fn_name"
        if isa(oplibpath, String)
            load_op_grad_dict[(oplibpath,opname)] = s
            verbose && printstyled("Load library operator (with gradient, multiple outputs = $multiple): $oplibpath ==> $opname\n", color=:green)
        end
        return s
end

"""
    load_system_op(opname::String, grad::Bool=true; multiple::Bool=false)

Loads custom operator from CustomOps directory (shipped with ADCME instead of TensorFlow)
For example
```
s = "SparseOperator"
oplib = "libSO"
grad = true
```
this will direct Julia to find library `CustomOps/SparseOperator/libSO.dylib` on MACOSX
"""
function load_system_op(opname::String, grad::Bool=true; multiple::Bool=false)
    if !isfile(LIBADCME)
        @info "$LIBADCME does not exist. Precompiling..."
        ADCME.precompile()
    end
    if grad
        load_op_and_grad(LIBADCME, opname; multiple=multiple, verbose=false)
    else
        load_op(LIBADCME, opname, verbose=false)
    end
end

"""
    compile(s::String; force::Bool=false)

Compiles the library given by path `deps/s`. If `force` is false, `compile` first check whether
the binary product exists. If the binary product exists, return 2. Otherwise, `compile` tries to
compile the binary product, and returns 0 if successful; it return 1 otherwise.
"""
function compile(s::String; force::Bool=false, customdir::Bool = false)
    PWD = pwd()
    dir = s
    if !customdir
        dir = joinpath(joinpath("$(@__DIR__)", "../deps/CustomOps"), s)
    end
    if !isdir(dir)
        @warn("Folder for the operator $s does not exist: $dir")
        return 1
    end
    cd(dir)

    local surfix
    if Sys.isapple()
        surfix = ".dylib"
    elseif Sys.islinux()
        surfix = ".so"
    elseif Sys.iswindows()
        surfix = ".dll"
    end
    if !force && isdir("build") # check if product exists
        files = readdir("build")
        if any([endswith(x, surfix) for x in files])
            @warn("The binary product exists.")
            cd(PWD)
            return 2
        end
    end
    rm("build",force=true,recursive=true)
    mkdir("build")
    cd("build")
    try
        cmake()
        make()
        cd(PWD)
        return 0
    catch e
        error("Compilation error: $e")
        cd(PWD)
        return 1
    end
end

"""
    precompile(force::Bool=false)

Precompile the built-in custom operators.
"""
function Base.:precompile(force::Bool=false)
    PWD = pwd()
    cd("$(@__DIR__)/../deps/CustomOps")
    if force
        try
            rm("build", force=true, recursive=true)
        catch
            error("""Failed to remove build directory. Follow the following steps and try again:
1. Quit ALL julia processes that use ADCME;
2. Remove $(joinpath(pwd(), "build")) manually.""")
        end
    end
    change_directory("build")
    require_cmakecache() do
        if Sys.isapple()
            try
                ADCME.cmake()
            catch
                @info "Use system clang..."
                mv(joinpath(BINDIR, "clang"), joinpath(BINDIR, "clang_original"))
                mv(joinpath(BINDIR, "clang++"), joinpath(BINDIR, "clang++_original"))
                symlink("/usr/bin/clang",joinpath(BINDIR, "clang"))
                symlink("/usr/bin/clang++",joinpath(BINDIR, "clang++"))
                ADCME.cmake()
            end
        else
            ADCME.cmake()
        end
    end
    require_library("adcme") do
        ADCME.make()
    end
    cd(PWD)
    global libadcme = tf.load_op_library(LIBADCME)
end

"""
    compile()

Compile a custom operator in the current directory. A `CMakeLists.txt` must be present.
"""
function compile()
    PWD = pwd()
    if !isfile("CMakeLists.txt")
        error(SystemError("No CMakeLists.txt in the current directory found."))
    end
    if !isdir("build")
        mkdir("build")
    end
    cd("build")
    try
        cmake()
        make()
    catch e
        @warn "Compiling failed: $e"
    finally
        cd(PWD)
    end
end

"""
    customop(;with_mpi::Bool = false)

Create a new custom operator. Typically users call `customop` twice: the first call generates a `customop.txt`,
users edit the content in the file; the second all generates C++ source code, CMakeLists.txt, and gradtest.jl from `customop.txt`.

# Example
```julia-repl
julia> customop() # create an editable `customop.txt` file
[ Info: Edit custom_op.txt for custom operators
julia> customop() # after editing `customop.txt`, call it again to generate interface files.
```

# Options
- `with_mpi`: Whether the custom operator uses MPI
"""
function customop(;with_mpi::Bool = false)
    # install_custom_op_dependency()
    py_dir = "$(@__DIR__)/../deps/CustomOpsTemplate"
    if !("custom_op.txt" in readdir("."))
        cp("$(py_dir)/custom_op.example", "custom_op.txt")
        @info "Edit custom_op.txt for custom operators"
        return
    else
        python = PyCall.python
        with_mpi = with_mpi ? 1 : 0
        run(`$python $(py_dir)/customop.py custom_op.txt $py_dir $with_mpi`)
    end
    nothing
end

function use_gpu(i::Union{Nothing,Int64}=nothing)
    if length(CUDA_INC)==0
        error("""ADCME is not built against GPU. Set ENV["GPU"]=1 and rebuild GPU.""")
    end
    dl = pyimport("tensorflow.python.client.device_lib")
    if !isnothing(i) && i>=1
        i = join(collect(0:i-1),',')
        ENV["CUDA_VISIBLE_DEVICES"] = i
    elseif !isnothing(i) && i==0
        ENV["CUDA_VISIBLE_DEVICES"] = ""
    end
    local_device_protos = dl.list_local_devices()
    return [x.name for x in local_device_protos if x.device_type == "GPU"]
end

function list_physical_devices(cpu_or_gpu::String = "all")
    dl = pyimport("tensorflow.python.client.device_lib")
    local_device_protos = dl.list_local_devices()
    CPU = [x.name for x in local_device_protos if x.device_type == "CPU"]
    GPU = [x.name for x in local_device_protos if x.device_type == "GPU"]
    if cpu_or_gpu == "all"
        return [CPU;GPU]
    elseif cpu_or_gpu == "GPU"
        return GPU
    elseif cpu_or_gpu == "CPU"
        return CPU
    else
        error(ArgumentError("$cpu_or_gpu is not a valid input. Expected: all, CPU, or GPU"))
    end
end

"""
    install(s::String; force::Bool = false, islocal::Bool = false)

Install a custom operator from a URL, a directory (when `islocal` is true), or a string. In any of the three case,
`install` copy the folder to $(abspath(joinpath(LIBADCME, "../../Plugin"))).
When `s` is a string, `s` is converted to

https://github.com/ADCMEMarket/<s>
"""
function install(s::String; force::Bool = false, islocal::Bool = false)
    if !islocal && !startswith(s, "http")
        s = "https://github.com/ADCMEMarket/"*s
    end
    _, name = splitdir(s)
    if force
        rm(joinpath(LIBPLUGIN, name), force=true, recursive=true)
    elseif isdir(joinpath(LIBPLUGIN, name)) && ("build" in readdir(joinpath(LIBPLUGIN, name)))
        return _plugin_lib(joinpath(LIBPLUGIN, name))
    end
    try
        LibGit2.clone(s, joinpath(LIBPLUGIN, name))
    catch
        LibGit2.clone("git://$(s[9:end]).git", joinpath(LIBPLUGIN, name))
    end
    PWD = pwd()
    cd(joinpath(LIBPLUGIN, name))
    if isfile(joinpath(LIBPLUGIN, name, "build.jl"))
        include(joinpath(LIBPLUGIN, name, "build.jl"))
    end
    cmakelists = String(read("CMakeLists.txt"))
    if !occursin("cmake_minimum_required", cmakelists)
        cmakelists = replace(String(read(joinpath(LIBPLUGIN, "CMakeLists.txt"))), "[INSTRUCTION]"=>cmakelists)
    end
    open("CMakeLists.txt", "w") do io
        write(io, cmakelists)
    end
    mkdir("build")
    cd("build")
    ADCME.cmake()
    ADCME.make()
    cd(PWD)
    return _plugin_lib(joinpath(LIBPLUGIN, name))
end

function _plugin_lib(D)
    files = readdir(joinpath(D, "build"))
    dylib = filter(x->endswith(x, ".$dlext"), files)
    if length(dylib)==0
        error(SystemError("No dynamic library found."))
    elseif length(dylib)>1
        error(SystemError("More then one dynamic libraries found."))
    else
        return joinpath(D, "build", dylib[1])
    end
end

@doc raw"""
    register(forward::Function, backward::Function; multiple::Bool=false)

Register a function `forward` with back-propagated gradients rule `backward` to the backward.
∘ `forward`: it takes $n$ inputs and outputs $m$ tensors. When $m>1$, the keyword `multiple` must be true.
∘ `backward`: it takes $\tilde m$ top gradients from float/double output tensors of `forward`, $m$ outputs of the `forward`,
   and $n$ inputs of the `forward`. `backward` outputs $n$ gradients for each input of `forward`. When input $i$ of
   `forward` is not float/double, `backward` should return `nothing` for the corresponding gradients.

# Example
```julia
forward = x->log(1+exp(x))
backward = (dy, y, x)->dy*(1-1/(1+y))
f = register(forward, backward)
```
"""
function register(forward::Function, backward::Function; multiple::Bool=false)
    fn_name = "customgrad_"*randstring(8)
    if !multiple
py"""
import tensorflow as tf
@tf.custom_gradient
def $$fn_name(*args):
    u = $forward(*args)
    def grad(dy):
        return $backward(dy, u, *args)
    return u, grad
"""
    else
py"""
import tensorflow as tf
@tf.custom_gradient
def $$fn_name(*args):
    u = forward_$$fn_name(*args)
    def grad(*dy):
        dy = [y for y in dy if y is not None and y.dtype in [tf.float64, tf.float32]] # only float64 and float32 can backpropagate gradients
        return backward_$$fn_name(*dy, *u, *args)
    return u, grad
"""
    end
    return py"$$fn_name"
end

"""
    debug(sess::PyObject, o::PyObject)

In the case a session run yields an error from the TensorFlow backend, this function can help print the exact error.
For example, you might encounter  `InvalidArgumentError()` with no detailed error information, and this function can be useful for debugging.
"""
function debug(sess::PyObject, o::PyObject)
py"""
import tensorflow as tf
import traceback
try:
    $sess.run($o)
except Exception:
    print(traceback.format_exc())
"""
end

"""
    debug(libfile::String = "")

Loading custom operator shared library. If the loading fails, detailed error message is printed.
"""
function debug(libfile::String = "")
    if libfile==""
        libfile = LIBADCME
    end
py"""
import tensorflow as tf
import traceback
try:
    tf.load_op_library($libfile)
except Exception:
    print(traceback.format_exc())
"""
end

"""
    doctor()

Reports health of the current installed ADCME package. If some components are broken, possible fix is proposed.
"""
function doctor()
    function yes(name)
        printstyled("[✔] $name\n", color=:green, bold=true)
    end
    function no(name, diagnose, instruction)
        printstyled("[✘] $name\n", color=:red, bold=true)
        printstyled("\n[Reason]\n", color=:magenta)
        printstyled("$diagnose\n\n", color=:blue)
        printstyled("\n[Instruction]\n", color=:magenta)
        printstyled("$instruction\n\n", color=:blue)
    end

    c = true
    if VERSION>=v"1.4" && VERSION<v"1.5" && Sys.isapple()
        c = false
    end

    if c
        yes("Julia version")
    else
        no("Julia version",
"""Your Julia version is $VERSION, and your system is MACOSX. This combination may have a compatibility issue.""",
"""If `using ADCME` freezes, consider install another version of Julia.""")
    end

    c = (length(tf.__version__)>=6) && (tf.__version__[1:6]=="1.15.0")
    if c
        yes("TensorFlow version")
    else
        no("TensorFlow version",
"""Your TensorFlow version is $(tf.__version__). The TensorFlow version shipped with ADCME is 1.15.0.""",
"""Set ENV["FORCE_REINSTALL_ADCME"] = 1 and rebuild ADCME
julia> ENV["FORCE_REINSTALL_ADCME"] = 1
julia> ]
pkg> build ADCME""")
    end

    c = (PyCall.python==ADCME.PYTHON)
    if c
        yes("Python executable file")
    else
        no("Python executable file",
"""PyCall Python path $(PyCall.python) does not match the ADCME-compatible Python $(ADCME.PYTHON)""",
"""Rebuild PyCall with a compatible Python version:

using Pkg
ENV["PYTHON"] = "$(ADCME.PYTHON)"
Pkg.build("PyCall")
""")
    end

    c = true
    try
        if Sys.iswindows()
            run(`cmd /c where julia`)
        else
            read(`which julia`)
        end
    catch
        c = false
    end

    if c
        yes("Julia path")
    else
        no("Julia path (Optional)",
"""`julia` outputs nothing. This will break custom operator compilation.""",
"""Add your julia binary path to your environment path, e.g. (Unix systems)

export PATH=$(Sys.BINDIR):\$PATH

For convenience, you can add the above line to your `~/.bashrc` (Linux) or `~/.bash_profile` (Apple).
For Windows, you need to add it to system environment.""")
    end

    c = Sys.iswindows() ?
        haskey(ENV, "PATH") && occursin(ADCME.LIBDIR, ENV["PATH"]) :
        haskey(ENV, "LD_LIBRARY_PATH") && occursin(ADCME.LIBDIR, ENV["LD_LIBRARY_PATH"])
    if c
        yes("Dynamic library path")
    else
        no("Dynamic library path (Optional)",
"""$(ADCME.LIBDIR) is not in LD_LIBRARY_PATH. This MAY break custom operator compilation. However, in most cases, ADCME automatic fixes this problem for you.""",
"""Add your dynamic library path path to your environment path, e.g. (Unix systems)

export LD_LIBRARY_PATH=$(ADCME.LIBDIR):\$LD_LIBRARY_PATH

For convenience, you can add the above line to your `~/.bashrc` (Linux or Apple).
For Windows, you need to add it to PATH instead of LD_LIBRARY_PATH.""")
    end

    c  = Sys.WORD_SIZE==64
    if c
        yes("Memory Address Length =  64")
    else
        no("Memory Address Length",
"""Your memory address length is $(Sys.WORD_SIZE). ADCME is only tested against 64-bit machine.""",
"""If you do not need custom operators, then it's fine. Otherwise you need to switch to a 64-bit machine""")
    end

    if Sys.iswindows()
        c = isfile(ADCME.MAKE*".exe") && occursin("15", (ADCME.MAKE)) && occursin("2017", ADCME.MAKE)
        if c
            yes("C Compiler")
        else
            no("C Compiler",
"""You specified that the C compiler for custom operators is
$(ADCME.MAKE)
However, one of the following requirements is not met:
1*. The file you specified $(ADCME.MAKE*".exe") does not exist.
2**. (Optional) For compatibility, we suggest you use Microsoft Visual Studio 2017 (Version number: 15).

* The path is actually not needed in compilation, but we raise such an issue here in case you obtain some compilation errors in the future.

* We check the version by looking for "15" and "2017" in the path specification. If you are sure your compiler is correct, you can ignore this message. """,
"""Manually edit $(abspath(joinpath(splitdir(pathof(ADCME))[1], "../deps/deps.jl"))) and modify `MAKE` to be the correct compiler.""")
        end

    end

    c = haskey(ENV, "PATH") && occursin(ADCME.BINDIR, ENV["PATH"])
    if c
        yes("Binaries path")
    else
        no("Binaries path",
"""$(ADCME.BINDIR) is not in PATH. This path contains compatible tools such as a GCC compiler, `cmake`, `make`, or any other tools you want to use directly from terminal.
However, setting the path is NOT a requirement, and ADCME works totally fine without any action.""",
"""(Optional) Add your binary path to your environment path, e.g. (Unix systems)

export PATH=$(ADCME.BINDIR):\$PATH

For convenience, you can add the above line to your `~/.bashrc` (Linux) or `~/.bash_profile` (Apple).
For Windows, you need to add it to system environment.""")
    end

    if length(ADCME.CUDA_INC)>0
        c = Sys.iswindows() ?
                haskey(ENV, "PATH") && occursin(ADCME.LIBCUDA, ENV["PATH"]) :
                haskey(ENV, "LD_LIBRARY_PATH") && occursin(ADCME.LIBCUDA, ENV["LD_LIBRARY_PATH"])
        if c
            yes("CUDA LD_LIBRARY_PATH")
        else
            no("CUDA LD_LIBRARY_PATH",
    """$(ADCME.LIBCUDA) is not in LD_LIBRARY_PATH. This path contains compatible tools such as a GCC compiler, `cmake`, `make`, etc.""",
    """The fix is OPTIONAL.
    Add your binary path to your environment path, e.g. (Unix systems)

    export LD_LIBRARY_PATH=$(ADCME.LIBCUDA):\$LD_LIBRARY_PATH

    For convenience, you can add the above line to your `~/.bashrc` (Linux or Apple).
    For Windows, you need to add it to PATH instead of LD_LIBRARY_PATH.""")
        end

        try
            if !Sys.iswindows()
                Libdl.dlpath("libcuda")
                Libdl.dlpath("libcudnn")
                Libdl.dlpath("libcublas")
            else
                Libdl.dlpath("cudart64_100")
                Libdl.dlpath("cudnn64_7")
                Libdl.dlpath("cublas64_100")
            end
            yes("CUDA Shared Library")
        catch
            no("CUDA Shared Library",
    """libcuda, libcudnn, and (or) libcublas can not be loaded.""",
    """If you intend to use GPU, this fix is mandatory. Make sure cudatoolkit and cudnn libraries can be found in
$(ADCME.LIBCUDA)
and `nvcc` is in your path.""")
        end

        c = isdir(ADCME.CUDA_INC) && "cuda.h" in readdir(ADCME.CUDA_INC)
        if c
            yes("CUDA Header Files")

            if !isfile(joinpath(ADCME.TF_INC, "third_party/gpus/cuda/include/cuda_fp16.h"))
                println("Fixing third_party/gpus/cuda/include...")
                if !ispath(joinpath(ADCME.TF_INC, "third_party/gpus/cuda/"))
                    mkpath(joinpath(ADCME.TF_INC, "third_party/gpus/cuda/"))
                end
                rm(joinpath(ADCME.TF_INC, "third_party/gpus/cuda/include/"), force=true, recursive=true)
                symlink(ADCME.CUDA_INC, joinpath(ADCME.TF_INC, "third_party/gpus/cuda/include"))
            end

        else
            no("CUDA Header Files",
            """Cuda include library does not exist or `cuda.h` is missing.""",
        """It might be possible that your cuda include library is located somewhere else other than $(ADCME.CUDA_INC). Fix the dependency file.""")
        end
    else
        no("GPU Support (Optional)",
    """ADCME is not compiled against GPU.""",
    """If you intend to use GPU, set ENV["GPU"] = 1 and then rebuild ADCME.""")
    end

    depsfile = abspath(joinpath(@__DIR__, "../deps/deps.jl"))
    println("============================================================\n$(depsfile)\n============================================================")
    println(read(depsfile, String))

    @assert isdir(ADCME.BINDIR)
    @assert isdir(ADCME.LIBDIR)
    @assert isdir(ADCME.TF_INC)
    @assert isdir(ADCME.PREFIXDIR)
    if !Sys.iswindows()
        @assert isfile(ADCME.CC) || islink(ADCME.CC)
        @assert isfile(ADCME.CXX) || islink(ADCME.CXX)
        @assert isfile(ADCME.NINJA)
    else
        @assert isfile(ADCME.MAKE)
    end
    @assert isfile(ADCME.CMAKE)
    @assert isfile(ADCME.TF_LIB_FILE) || islink(ADCME.TF_LIB_FILE)
    @assert isdir(ADCME.INCDIR)
end

"""
    test_gpu()

Tests the GPU ultilities
"""
function test_gpu()
    PWD = pwd()
    run(`which nvcc`)
    cd("$(@__DIR__)/../examples/gpu_custom_op")
    mkdir("build")
    cd("build")
    ADCME.cmake()
    ADCME.make()
    cd("..")
    include("gputest.jl")
    cd(PWD)
end

@doc raw"""
    nnuq(H::Array{Float64,2}, invR::Union{Float64, Array{Float64,2}}, invQ::Union{Float64, Array{Float64,2}})

Returns the variance matrix for the Baysian inversion.

The negative log likelihood function is

$$l(s) =\frac{1}{2} (y-h(s))^T R^{-1} (y-h(s)) + \frac{1}{2} s^T Q^{-1} s$$

The covariance matrix is computed by first linearizing $h(s)$

$$h(s)\approx h(s_0) + \nabla h(s_0) (s-s_0)$$

and then computing the second order derivative

$$V = \left(\frac{\partial^2 l}{\partial s^T\partial s}\right)^{-1} = (H^T R^{-1} H + Q^{-1})^{-1}$$

Note the result is independent of $s_0$, $y_0$, and only depends on $\nabla h(s_0)$
"""
function nnuq(H::Array{Float64,2}, invR::Union{Float64, Array{Float64,2}}, invQ::Union{Float64, Array{Float64,2}})
    if isa(invQ, Float64)
        invQ = invQ * I
    end
    Σ = inv(H' * invR * H + invQ)
    (Σ + Σ')/2
end

mutable struct MCMCSimple
    logf::Function
    proposal::Function
    θ0::Array{Float64, 1}
    ub::Float64
    lb::Float64
    L::Array{Float64, 1}
    AC::Array{Float64, 1}
end

@doc raw"""
    MCMCSimple(obs::Array{Float64, 1}, h::Function,
    σ::Float64, θ0::Array{Float64,1}, lb::Float64, ub::Float64)

A very simple yet useful interface for MCMC simulation in many scientific computing problems.

- `obs`: Observations
- `h`: Forward computation function
- `σ`: Noise standard deviation for the observed data
- `ub`, `lb`: upper and lower bound
- `θ0`: Initial guess

The mathematical model is

```math
y_{obs} = h(\theta)
```

and we have a hard constraint `lb\leq \theta \leq ub`.
"""
function MCMCSimple(obs::Array{Float64, 1}, h::Function,
        σ::Float64, θ0::Array{Float64,1}, lb::Float64, ub::Float64,
        δ::Union{Missing, Float64}=missing)
    τ = (ub-lb)/6
    δ = coalesce(δ, (ub-lb)/100)
    function logf(x)
        -sum((obs - h(x)).^2/2σ^2) - sum((x-θ0).^2)/2τ^2
    end
    function proposal(x)
        x + (rand(length(x)) .- 0.5)*2 * δ
    end
    MCMCSimple(logf, proposal, θ0, ub, lb, zeros(0), zeros(0))
end

function simulate(ms::MCMCSimple, N::Int64, burnin::Union{Int64, Missing} = missing)
    burnin = coalesce(burnin, Int64(round(N*0.2)))
    sim = zeros(N, length(ms.θ0))
    sim[1,:] = ms.θ0
    L = zeros(N)
    AC = ones(N)
    AC[1] = NaN
    L[1] = ms.logf(ms.θ0)
    k = 1
    for i = 2:N
        sim[i,:], L[i], k_ = _MCMCSimple_simulate(ms, sim[i-1,:])
        k += k_
        AC[i] = k/i
    end
    ms.AC = AC
    ms.L = L
    return sim
end

function diagnose(ms::MCMCSimple)
    if !isdefined(Main, :PyPlot)
        error("Package PyPlot.jl must be imported in the main module using `import PyPlot` or `using PyPlot`")
    end
    Main.PyPlot.figure(figsize = (10, 4))
    Main.PyPlot.subplot(121)
    Main.PyPlot.title("Acceptance Rate")
    Main.PyPlot.plot(ms.AC)
    Main.PyPlot.ylim(0,1.05)
    Main.PyPlot.subplot(122)
    Main.PyPlot.title("Log Likelihood")
    Main.PyPlot.plot(ms.L)
end

function _MCMCSimple_simulate(ms::MCMCSimple, x::Array{Float64})
    local x_star
    while true
        x_star = ms.proposal(x)
        if all(x_star.<=ms.ub) && all(x_star.>=ms.lb)
            break
        end
    end
    Δ =  ms.logf(x_star) - ms.logf(x)
    if log(rand())<Δ
        return x_star, ms.logf(x_star), 1
    else
        return x, ms.logf(x), 0
    end
end

"""
    get_placement()

Returns the operation placements.
"""
function get_placement()
    sess = Session(config=tf.ConfigProto(log_device_placement=true))
    originalSTDOUT = stdout
    (outRead, outWrite) = redirect_stdout()
    init(sess)
    close(outWrite)
    data = readavailable(outRead)
    close(outRead)
    redirect_stdout(originalSTDOUT)
    lines = split(String(data), '\n')[1:end-1]
end

"""
    sleep_for(t::Union{PyObject, <:Real})

Sleeps for `t` seconds.
"""
function sleep_for(t::Union{PyObject, <:Real})
    sleep_for_ = load_system_op("sleep_for", false)
    sleep_for_(convert_to_tensor(t, dtype=Float64))
end

"""
    timestamp(deps::Union{PyObject, <:Real, Missing}=missing)

These functions are usually used with [`bind`](@ref) for profiling.
Note the timing is not very accurate in a multithreaded environment.

- `deps`: `deps` is always executed before returning the timestamp.
# Example
```julia
a = constant(3.0)
t0 = timestamp(a)
sleep_time = sleep_for(a)
t1 = timestamp(sleep_time)
sess = Session(); init(sess)
t0_, t1_ = run(sess, [t0, t1])
time = t1_ - t0_
```
"""
function timestamp(deps::Union{PyObject, <:Real, Missing}=missing)
    deps = coalesce(deps, 0.0)
    deps = convert_to_tensor(deps, dtype = Float64)
    timer_ = load_system_op("timer", false)
    timer_(deps)
end

"""
    load_library(filename::String)

Load custom operator libraries. If used with
"""
function load_library(filename::String)
    filename = get_library(filename)
    keyname = "custom_op_library_"*filename
    if haskey(STORAGE, keyname)
        return STORAGE[keyname]
    end
    if !isfile(filename)
        error("File $filename not found. If you intend to load a library by absolute path, try `tf.load_op_library(filename)`.")
    end
    try
        STORAGE[keyname] =  tf.load_op_library(filename)
    catch e
        error("Failed to load library: $filename. Original error message:\n$e")
    end
    return STORAGE[keyname]
end

"""
    get_library_symbols(file::Union{String, PyObject})

Returns the symbols in the custom op library `file`.
"""
function get_library_symbols(files::Union{String, PyObject})
    if isa(files, String)
        files = load_library(files)
    end
    files = keys(files)
    filter(x->islowercase(String(x)[1]) && String(x)!="tf_export" && !(occursin("fallback", String(x))) && String(x)!="deprecated_endpoints", files)
end

function Base.:NamedTuple(df::PyObject)
    names = tuple(Symbol.(df.index.tolist())...)
    vals = tuple(df.values...)
    NamedTuple{names}(vals)
end


text:
#
# types.jl -
#
# Type definitions and (some) constructors for linear algebra.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

struct SingularSystem <: Exception
    msg::String
end
showerror(io::IO, err::SingularSystem) =
    print(io, "singular linear system (", err.msg, ")")

struct NonPositiveDefinite <: Exception
    msg::String
end
showerror(io::IO, err::NonPositiveDefinite) =
    print(io, "non-positive definite operator (", err.msg, ")")

struct UnimplementedOperation <: Exception
    msg::String
end
showerror(io::IO, err::UnimplementedOperation) =
    print(io, err.msg)

struct UnimplementedMethod <: Exception
    msg::String
end
showerror(io::IO, err::UnimplementedMethod) =
    print(io, err.msg)

"""

Type `Reals` is the set of the floating point types.  It is the numerical
approximation of reals in the mathematical sense.

This definition closely follows the semantic used in the BLAS module that
`BlasReal` are all the real types supported by the BLAS library.

"""
const Reals = AbstractFloat

"""

Type `Complexes` is the set of the complexes whose real and imaginary parts are
floating point.  It is the numerical approximation of complexes in the
mathematical sense.

This definition closely follows the semantic used in the BLAS module that
`BlasComplex` are all the complex types supported by the BLAS library.

"""
const Complexes = Complex{<:Reals}

"""

Type `Floats` is the union of all floating-point types (reals and complexes).

This definition closely follows the semantic used in the BLAS module that
`BlasFloat` are all floating-point types supported by the BLAS library.

"""
const Floats = Union{Reals,Complexes}

"""

A `Mapping` is any function between two variables spaces.  Assuming upper case
Latin letters denote mappings, lower case Latin letters denote variables, and
Greek letters denote scalars, then:

* `A*x` or `A⋅x` yields the result of applying the mapping `A` to `x`;

* `A\\x` yields the result of applying the inverse of `A` to `x`;

Simple constructions are allowed for any kind of mappings and can be used to
create new instances of mappings which behave correctly.  For instance:

* `B = α*A` (where `α` is a real) is a mapping which behaves as `A` times `α`;
  that is `B⋅x` yields the same result as `α*(A⋅x)`.

* `C = A + B + ...` is a mapping which behaves as the sum of the mappings `A`,
  `B`, ...; that is `C⋅x` yields the same result as `A⋅x + B⋅x + ...`.

* `C = A*B` or `C = A⋅B` is a mapping which behaves as the composition of the
  mappings `A` and `B`; that is `C⋅x` yields the same result as `A⋅(B.x)`.  As
  for the sum of mappings, there may be an arbitrary number of mappings in a
  composition; for example, if `D = A*B*C` then `D⋅x` yields the same result as
  `A⋅(B⋅(C⋅x))`.

* `C = A\\B` is a mapping such that `C⋅x` yields the same result as `A\\(B⋅x)`.

* `C = A/B` is a mapping such that `C⋅x` yields the same result as `A⋅(B\\x)`.

These constructions can be combined to build up more complex mappings.  For
example:

* `D = A*(B + C)` is a mapping such that `C⋅x` yields the same result as
  `A⋅(B⋅x + C⋅x)`.

A `LinearMapping` is any linear mapping between two spaces.  This abstract
subtype of `Mapping` is introduced to extend the notion of *matrices* and
*vectors*.  Assuming the type of `A` inherits from `LinearMapping`, then:

* `A'⋅x` and `A'*x` yields the result of applying the adjoint of the mapping
  `A` to `x`;

* `A'\\x` yields the result of applying the adjoint of the inverse of mapping
  `A` to `x`.

* `B = A'` is a mapping such that `B⋅x` yields the same result as `A'⋅x`.

The following methods should be implemented for a mapping `A` of specific type
`M <: Mapping`:

```julia
vcreate(::Type{P}, A::M, x, scratch::Bool) -> y
apply!(α::Number, ::Type{P}, A::M, x, , scratch::Bool, β::Number, y) -> y
```

for any supported operation `P ∈ Operations` (`Direct`, `Adjoint`, `Inverse`
and/or `InverseAdjoint`).  See the documentation of these methods for
explanations.  Optionally, methods `P(A)` may be extended, *e.g.* to throw
exceptions if operation `P` is forbidden (or not implemented).  By default, all
these operations are assumed possible (except `Adjoint` and `InverseAdjoint`
for a nonlinear mapping).

See also: [`apply`](@ref), [`apply!`](@ref), [`vcreate`](@ref),
          [`LinearType`](@ref), [`Scalar`](@ref), [`Direct`](@ref),
          [`Adjoint`](@ref), [`Inverse`](@ref), [`InverseAdjoint`](@ref).

"""
abstract type Mapping <: Function end

abstract type LinearMapping <: Mapping end
@doc @doc(Mapping) LinearMapping

"""

    Identity()

yields the identity linear mapping.  The purpose of this mapping is to be as
efficient as possible, hence the result of applying this mapping may be the
same as the input argument.

The identity is a singleton and is also available as:

    Id

The `LinearAlgebra` module of the standard library exports a constant `I` which
also corresponds to the identity (but in the sense of a matrix).  When `I` is
combined with any LazyAlgebra mapping, it is recognized as an alias of `Id`.
So that, for instance, `I/A`, `A\\I`, `Id/A` and `A\\Id` all yield `inv(A)` for
any LazyAlgebra mappings `A`.

"""
struct Identity <: LinearMapping; end

"""

The abstract type `Trait` is inherited by types indicating specific traits.

See also: [`LinearType`](@ref), [`SelfAdjointType`](@ref),
          [`DiagonalType`](@ref), [`MorphismType`](@ref).

"""
abstract type Trait end

# Trait indicating whether a mapping is certainly linear.
abstract type LinearType <: Trait end
struct NonLinear <: LinearType end
struct Linear <: LinearType end

# Trait indicating whether a mapping is certainly a self-adjoint linear map.
abstract type SelfAdjointType <: Trait end
struct NonSelfAdjoint <: SelfAdjointType end
struct SelfAdjoint <: SelfAdjointType end

# Trait indicating whether a mapping is certainly an endomorphism.
abstract type MorphismType <: Trait end
struct Morphism <: MorphismType end
struct Endomorphism <: MorphismType end

# Trait indicating whether a mapping is certainly a diagonal linear mapping.
abstract type DiagonalType <: Trait end
struct NonDiagonalMapping <: DiagonalType end
struct DiagonalMapping <: DiagonalType end

"""

Type `Direct` is a singleton type to indicate that a linear mapping should
be directly applied.  This type is part of the union `Operations`.

See also: [`LinearMapping`](@ref), [`apply`](@ref), [`Operations`](@ref).

"""
struct Direct; end

"""

    Adjoint(A) -> obj

yields an object instance `obj` representing `A'`, the adjoint of the linear
mapping `A`.

Directly calling this constructor is discouraged, use an expression like `A'`
instead and benefit from automatic simplification rules.

Call [`unveil(obj)`](@ref) to reveal the linear mapping `A` embedded in `obj`.

See also [`DecoratedMapping`](@ref).

"""
struct Adjoint{T<:Mapping} <: LinearMapping
    op::T

    # The outer constructors prevent most illegal calls to `Adjoint(A)` we
    # just have to check that the argument is a simple linear mapping.
    function Adjoint{T}(A::T) where {T<:Mapping}
        is_linear(A) || throw_forbidden_adjoint_of_non_linear_mapping()
        return new{T}(A)
    end
end

"""
    Inverse(A) -> obj

yields an object instance `obj` representing the inverse of the mapping `A`.

Directly calling this constructor is discouraged, call `inv(A)` or use an
expression like `Id/A` instead and benefit from automatic simplification rules.

Call [`unveil(obj)`](@ref) to reveal the mapping `A` embedded in `obj`.

See also [`DecoratedMapping`](@ref).

"""
struct Inverse{T<:Mapping} <: Mapping
    op::T

    # The outer constructors prevent all illegal calls to `Inverse(A)` so there
    # is nothing more to check.
    Inverse{T}(A::T) where {T<:Mapping} = new{T}(A)
end

"""
    InverseAdjoint(A) -> obj

yields an object instance `obj` representing the inverse of the adjoint of the
linear mapping `A`.

Directly calling this constructor is discouraged, use expressions like
`inv(A')`, `inv(A')` or `Id/A'` instead and benefit from automatic
simplification rules.

Call [`unveil(obj)`](@ref) to reveal the mapping `A` embedded in `obj`.

`AdjointInverse` is an alias for `InverseAdjoint`.

See also [`DecoratedMapping`](@ref).

"""
struct InverseAdjoint{T<:Mapping} <: LinearMapping
    op::T

    # The outer constructors prevent most illegal calls to `InverseAdjoint(A)`
    # we just have to check that the argument is a simple linear mapping.
    function InverseAdjoint{T}(A::T) where {T<:Mapping}
        is_linear(A) ||
            bad_argument("taking the inverse adjoint of non-linear mappings is not allowed")
        return new{T}(A)
    end
end

const AdjointInverse{T} = InverseAdjoint{T}
@doc @doc(InverseAdjoint) AdjointInverse

"""

    Gram(A) -> obj

yields an object instance `obj` representing the composition `A'*A` for the
linear mapping `A`.

Directly calling this constructor is discouraged, call [`gram(A)`](@ref) or use
expression `A'*A` instead and benefit from automatic simplification rules.

Call [`unveil(obj)`](@ref) to reveal the linear mapping `A` embedded in `obj`.

See also [`gram`](@ref), [`unveil`](@ref) and [`DecoratedMapping`](@ref).

"""
struct Gram{T<:Mapping} <: LinearMapping
    op::T

    # The outer constructors prevent most illegal calls to `Gram(A)` we
    # just have to check that the argument is a simple linear mapping.
    function Gram{T}(A::T) where {T<:Mapping}
        is_linear(A) || throw_forbidden_Gram_of_non_linear_mapping()
        return new{T}(A)
    end
end

"""

`DecoratedMapping` is the union of the *decorated* mapping types:
[`Adjoint`](@ref), [`Inverse`](@ref), [`InverseAdjoint`](@ref), and
[`Gram`](@ref).

The method [`unveil(A)`](@ref) can be called to reveal the mapping embedded in
a decorated mapping `A`.

"""
const DecoratedMapping = Union{Adjoint,Inverse,InverseAdjoint,Gram}

"""

    Jacobian(A,x) -> obj

yields an object instance `obj` representing the Jacobian `∇(A,x)` of the non-linear
mapping `A` for the variables `x`.

Directly calling this constructor is discouraged, call [`jacobian(A,x)`](@ref) or
[`∇(A,x)`](@ref) instead and benefit from automatic simplification rules.

"""
struct Jacobian{M<:Mapping,T} <: Mapping
    A::M
    x::T

    # The outer constructors prevent most illegal calls to `Jacobian(A)` we
    # just have to check that the argument is not a simple linear mapping.
    function Jacobian{M,T}(A::M, x::T) where {M<:Mapping,T}
        is_linear(A) &&
            bad_argument("the Jacobian of a linear mapping of type `",
                         M, "` should be the mapping itself")
        return new{M,T}(A, x)
    end
end

"""

`Operations` is the union of the possible variants to apply a mapping:
[`Direct`](@ref), [`Adjoint`](@ref), [`Inverse`](@ref) and
[`InverseAdjoint`](@ref) (or its alias [`AdjointInverse`](@ref)).

See also: [`apply`](@ref) and [`apply!`](@ref).

"""
const Operations = Union{Direct,Adjoint,Inverse,InverseAdjoint}

"""

    Scaled(λ, M) -> obj

yields an object instance `obj` representing `λ*M`, the mapping `M` multiplied
by a scalar `λ`.

Directly calling this constructor is discouraged, use expressions like `λ*M`
instead and benefit from automatic simplification rules.

Call [`multiplier(obj)`](@ref) and [`unscaled(obj)`](@ref) with a scaled
mapping `obj = λ*M` to retrieve `λ` and `M` respectively.

"""
struct Scaled{T<:Mapping,S<:Number} <: Mapping
    λ::S
    M::T
    Scaled{T,S}(λ::S, M::Mapping) where {S<:Number,T<:Mapping} =
        new{T,S}(λ, M)
end

"""

    Sum(A, B, ...) -> obj

yields an object instance `obj` representing the sum `A + B + ...` of the
mappings `A`, `B`, ...

Directly calling this constructor is discouraged, use expressions like `A + B +
...` instead and benefit from automatic simplification rules.

Call [`terms(obj)`](@ref) retrieve the tuple `(A,B,...)` of the terms of the
sum stored in `obj`.

"""
struct Sum{N,T<:NTuple{N,Mapping}} <: Mapping
    ops::T

    # The inner constructor ensures that the number of arguments is at least 2.
    function Sum{N,T}(ops::T) where {N,T<:NTuple{N,Mapping}}
        N ≥ 2 ||
            throw(ArgumentError("a sum of mappings has at least 2 components"))
        new{N,T}(ops)
    end
end

"""

    Composition(A, B, ...) -> obj

yields an object instance `obj` representing the composition `A*B*...` of the
mappings `A`, `B`, ...

Directly calling this constructor is discouraged, use expressions like
`A*B*...` `A∘B∘...` or `A⋅B⋅...` instead and benefit from automatic
simplification rules.

Call [`terms(obj)`](@ref) retrieve the tuple `(A,B,...)` of the terms of the
composition stored in `obj`.

"""
struct Composition{N,T<:NTuple{N,Mapping}} <: Mapping
    ops::T

    # The inner constructor ensures that the number of arguments is at least 2.
    function Composition{N,T}(ops::T) where {N,T<:NTuple{N,Mapping}}
        N ≥ 2 ||
            throw(ArgumentError("a composition of mappings has at least 2 components"))
        new{N,T}(ops)
    end
end


text:
# The prime factors of 13195 are 5, 7, 13 and 29.
# What is the largest prime factor of the number 600851475143 ?

num = BigInt(600851475143)
#num = BigInt(13195)

is_prime(n) = !any(n%i == 0 for i in 2:BigInt(floor(sqrt(n))))

function largest_factor(num)
    #step_ = (num % 2 == 0) ? -1 : -2
    num_leftover = num
    for p in (BigInt(p) for p=2:num_leftover if is_prime(p))
        if p == num_leftover
            println(p)
            return p
        elseif num_leftover % p == 0
            num_leftover = BigInt(num_leftover / p)
        end
    end
end
largest_factor(num)
println()

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
# advent of code 2020 day 2

function load_data(path:: String):: AbstractArray
    println(path)
    file = open(path)
    data_array = readlines(file)
    close(file)
    return data_array
end

function split_first(qstring:: String, qchar:: String)
    split_qstring = split(qstring, qchar)
    first_match = split_qstring[1]
    rest = join(split_qstring[2: length(split_qstring)])
    return first_match, rest
end

function check_pw_policy1(pw, qchar, low_bound, up_bound)
    low_bound, up_bound = parse(Int, low_bound), parse(Int, up_bound)
    num_match = length(findall(qchar, pw))
    if low_bound <= num_match <= up_bound
        return 1
    else
        return 0
    end

end

function check_pw_policy2(pw, qchar, low_bound, up_bound)
    low_bound, up_bound = parse(Int, low_bound), parse(Int, up_bound)
    char1, char2 = pw[low_bound:low_bound], pw[up_bound:up_bound]
    if (char1 == qchar) ⊻ (char2 == qchar)
        return 1
    else
        return 0
    end

end

data_array = load_data("./data/day2.txt")
use_policy = 3

let num_valid_pw = 0
    if use_policy == 1
        check_pw = check_pw_policy1
    elseif use_policy ==  2
        check_pw = check_pw_policy2
    else
        println("Policy: $use_policy Not implemented, defaulting to policy 1")
        global use_policy = 1
        check_pw = check_pw_policy1
    end

    for rule_pw in data_array
        low_bound, rule_pw = split_first(rule_pw, "-")
        up_bound, rule_pw = split_first(rule_pw, " ")
        char, pw = split_first(rule_pw, ":")
        num_valid_pw += check_pw(pw, char, low_bound, up_bound)
    end
    println("number of valid pw: $num_valid_pw, using policy $use_policy")
end


text:
using CoDa
using Tables
using DataDeps
using Distances
using ScientificTypes
using RData
using CSV
using Test
using LinearAlgebra

# accept downloads without interaction
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
isappveyor = "APPVEYOR" ∈ keys(ENV)
isCI = istravis || isappveyor
datadir = joinpath(@__DIR__,"data")

# download and setup data dependencies
register(DataDep("juraset",
  "A geochemical dataset from the Swiss Jura",
  "https://github.com/cran/compositions/raw/master/data/juraset.rda",
  "0fcf23fbca1d6fcb58ae0de6f11365f39fa3df02828128708185ebf45b002382"))
rda  = joinpath(datadep"juraset", "juraset.rda")
jura = RData.load(rda)["juraset"]
CSV.write(joinpath(datadir,"jura.csv"), jura)

# list of tests
testfiles = [
  "compositions.jl",
  "codaarrays.jl",
  "transforms.jl",
  "covariances.jl",
  "distances.jl",
  "matrices.jl",
  "scitypes.jl"
]

@testset "CoDa.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


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

module ModelReduction

using Reexport
@reexport using FEMBase

include("guyan_reduction.jl")
include("craig_bampton.jl")
include("model_reduction_craig_bampton.jl")
export CraigBampton
end


text:
using DebuggableCompiler
using Test

@testset "DebuggableCompiler.jl" begin
    # Write your tests here.
end


text:
using PredictiveControl
using Test

H = [ 10.0 0.0;
       0.0 2.0 ];

b = [ 2.0;
      4.0 ];

# Run for only 40 iterations
x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, stopconditions = nothing, maxiter = 40, dispinterval = 0 )
@test iter == 40

# Run until converged
x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, stopconditions = [FGM.Conjugate(1e-8)], dispinterval = 0 )
@test isapprox( x, [-0.2, -2.0], atol = 1e-4 )

x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, stopconditions = [FGM.Gradient(1e-8)], dispinterval = 0 )
@test isapprox( x, [-0.2, -2.0], atol = 1e-4 )

x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, stopconditions = [FGM.Best(1e-8)], dispinterval = 0 )
@test isapprox( x, [-0.2, -2.0], atol = 1e-4 )

# Start on the actual point - should converge very fast
x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, x₀ = [-0.2; -2.0], stopconditions = [FGM.Best(1e-8)], dispinterval = 0 )
@test iter <= 2
@test isapprox( x, [-0.2, -2.0], atol = 1e-4 )

# Test a clipped region (project onto constraints)
x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> max.(-1.0, min.(1.0, x)), stopconditions = [FGM.Best(1e-8)], dispinterval = 0 )
@test isapprox( x, [-0.2, -1.0], atol = 1e-4 )

# Test the variable steping
x, iter = PredictiveControl.FastGradientMethod.fastgradientmethod( H, b, x -> x, step=FGM.VariableStep(), stopconditions = [FGM.Best(1e-8)], dispinterval = 0 )
@test isapprox( x, [-0.2, -2.0], atol = 1e-4 )


text:
export AbstractMetricReader, CompositMetricReader, MetricReader, PeriodicMetricReader

"""
All metric readers should implement `shut_down!(::AbstractMetricReader)` and `(r::AbstractMetricReader)()`.

Builtin readers:

  - [`CompositMetricReader`](@ref)
  - [`MetricReader`](@ref)
  - [`PeriodicMetricReader`](@ref)
"""
abstract type AbstractMetricReader end

struct CompositMetricReader <: AbstractMetricReader
    readers::Vector{Any}
    function CompositMetricReader(readers...)
        new([readers...])
    end
end

function (r::CompositMetricReader)()
    for x in r.readers
        x()
    end
end

function shut_down!(r::CompositMetricReader)
    for x in r.readers
        shut_down!(x)
    end
end

#####

"""
    MetricReader([global_meter_provider()], [ConsoleExporter()])

Note that all metrics will be read on initialization.
"""
struct MetricReader{P,E} <: AbstractMetricReader
    provider::P
    exporter::E
    function MetricReader(
        provider::P,
        exporter::E,
    ) where {P<:AbstractMeterProvider,E<:AbstractExporter}
        r = new{P,E}(provider, exporter)
        r()
        r
    end
end

MetricReader(p::AbstractMeterProvider) = MetricReader(p, ConsoleExporter())
MetricReader(e::AbstractExporter) = MetricReader(global_meter_provider(), e)
MetricReader() = MetricReader(global_meter_provider(), ConsoleExporter())

"""
    (r::MetricReader)()

For async instruments in `r`, their callbacks will be executed first before reading all the metrics.
"""
function (r::MetricReader)()
    for ins in keys(r.provider.async_instruments)
        ins()
    end
    export!(r.exporter, metrics(r.provider))
end

# ??? shut_down! provider?
shut_down!(r::MetricReader) = shut_down!(r.exporter)

#####

"""
    PeriodicMetricReader(reader; export_interval_seconds = 60, export_timeout_seconds = 30)

Periodically call the `reader` in the background.
"""
struct PeriodicMetricReader{R<:AbstractMetricReader} <: AbstractMetricReader
    reader::R
    export_interval_seconds::Int
    export_timeout_seconds::Int
    timer::Timer
    function PeriodicMetricReader(
        reader;
        export_interval_seconds = 60,
        export_timeout_seconds = 30,
    )
        timer = Timer(0; interval = export_interval_seconds) do t
            res = timedwait(export_timeout_seconds; pollint = 1) do
                reader()
                true
            end
            if res == :timed_out
                @warn "timed out when exporting metrics"
                # how to stop exporting?
            end
        end

        new{typeof(reader)}(reader, export_interval_seconds, export_timeout_seconds, timer)
    end
end

"""
For `PeriodicMetricReader`, there's no need to call this method since the read
operation will be done periodically in the background.
"""
function (r::PeriodicMetricReader)() end

function shut_down!(r::PeriodicMetricReader)
    close(r.timer)
    shut_down!(r.reader)
end


text:
@doc """
Calculates the F test as is commonly implemented in SSR research.
TODO: Add references to MASTER and Luts et al

#### Parameters

* Sweep measurements. Samples x Sweeps x Channels
* Frequency(ies) of interest (Hz)
* Sampling rate (Hz)
* The amount of data to use on each side of frequency of interest to estimate noise (Hz)
* Filter used on the sweep data. If provided then is compensated for
* The number of bins to ignore on each side of the frequency of interest

#### Returns

* Signal to noise ratio in dB
* Signal phase at frequency of interest
* Signal power at frequency of interest
* Noise power estimated of side frequencies
* F statistic

""" ->
function ftest(sweeps::Union{Array{Float64, 3}, Array{Float32, 3}}, freq_of_interest::Real,
                fs::Real, side_freq::Real, used_filter::Union{Filter, Void}, spill_bins::Int)

    spectrum    = EEG._ftest_spectrum(sweeps)
    # No compensation is made here for prior filtering. See SSR.jl for an example of how to compensate filtering.
    frequencies = linspace(0, 1, Int(size(spectrum, 1)))*float(fs)/2

    ftest(spectrum, frequencies, freq_of_interest, side_freq, spill_bins)
end

function ftest{T <: AbstractFloat}(spectrum::Array{Complex{T},2}, frequencies::AbstractArray,
            freq_of_interest::Real, side_freq::Real, spill_bins::Int)

    Logging.info("Calculating F statistic on $(size(spectrum)[end]) channels at $freq_of_interest Hz +-$(side_freq) Hz")

    idx      = _find_closest_number_idx(frequencies, freq_of_interest)
    idx_Low  = _find_closest_number_idx(frequencies, freq_of_interest - side_freq)
    idx_High = _find_closest_number_idx(frequencies, freq_of_interest + side_freq)

    # Determine signal phase
    signal_phase = angle(spectrum[idx, :])                             # Biased response phase

    # Determine signal power
    signal_power = vec(abs(spectrum[idx, :]).^2)                       # Biased response power

    # Determine noise power
    noise_idxs  = [idx_Low - div(spill_bins, 2) : idx - spill_bins; idx + spill_bins : idx_High + div(spill_bins, 2)]
    noise_bins  = spectrum[noise_idxs,:]
    noise_bins  = abs(noise_bins)
    noise_power = vec(sum(noise_bins .^2, 1) ./ size(noise_bins,1))     # Recording noise power

    # Calculate SNR
    snr = (signal_power ./ noise_power)                                # Biased recording SNR
    snrDb = 10 * log10(snr)

    # Calculate statistic
    continuous_distribution = FDist(2, 2*size(noise_bins,1))
    statistic = ccdf(continuous_distribution, snr)

    # Debugging information
    debug("Frequencies = [$(freq_of_interest - side_freq), $(freq_of_interest), $(freq_of_interest + side_freq)]")
    debug("Indicies    = [$(minimum(noise_idxs)), $(idx), $(maximum(noise_idxs))]")
    debug("Noise bins  = $(size(noise_bins,1))")
    debug("Signal      = $(signal_power)")
    debug("Noise       = $(noise_power)")
    debug("SNR         = $(snr)")
    debug("SNR dB      = $(snrDb)")
    debug("Stat        = $(statistic)")

    return snrDb, signal_phase, signal_power, noise_power, statistic
end

# Calculates the spectrum for ftest and plotting
function _ftest_spectrum(sweep::Union{Array{Float64,1}, Array{Float64,2}}; ref::Int=0)
    # First dimension is samples, second dimension if existing is channels

    sweepLen      = size(sweep)[1]

    # Calculate amplitude sweepe at each frequency along first dimension
    fftSweep    = 2 / sweepLen * fft(sweep, 1)
    spectrum    = fftSweep[1:round(Int, sweepLen / 2 + 1), :]

    if ref > 0
        refspec = spectrum[:,ref]
        for i = 1:size(spectrum)[2]
            spectrum[:,i] = spectrum[:,i] - refspec
        end
    end

    return spectrum
end

function _ftest_spectrum(sweeps::Array{Float64,3};ref=0); _ftest_spectrum(squeeze(mean(sweeps,2),2),ref=ref); end
#=function _ftest_spectrum(s::Array{Float32}; ref=0); _ftest_spectrum(convert(Array{AbstractFloat}, s), ref=ref); end=#
#=function _ftest_spectrum(s::Array{Float64}; ref=0); _ftest_spectrum(convert(Array{AbstractFloat}, s), ref=ref); end=#


text:
@testset "basic tests" begin

Base.@kwdef mutable struct TestActor
    state::Union{Nothing, Int} = nothing
end

struct CurrentStateMsg <: AbstractMessage
    state
end

Base.@kwdef struct ReadStateMsg <: AbstractMessage
    from = self()
end

struct IncMsg <: AbstractMessage end
struct DecMsg <: AbstractMessage end

(x::TestActor)(msg::StartMsg{Tuple{Int}}) = x.state = msg.args[1]
(x::TestActor)(msg::StopMsg) = x.state = nothing
(x::TestActor)(::IncMsg) = x.state += 1
(x::TestActor)(::DecMsg) = x.state -= 1
(x::TestActor)(msg::ReadStateMsg) = put!(msg.from, CurrentStateMsg(x.state))

x = actor(TestActor())
put!(x, StartMsg(0))

put!(x, ReadStateMsg())
@test take!(self()).state == 0

@sync begin
    for _ in 1:100
        Threads.@spawn put!(x, IncMsg())
        Threads.@spawn put!(x, DecMsg())
    end
    for _ in 1:10
        for _ in 1:10
            Threads.@spawn put!(x, IncMsg())
        end
        for _ in 1:10
            Threads.@spawn put!(x, DecMsg())
        end
    end
end

put!(x, ReadStateMsg())
@test take!(self()).state == 0

y = actor(TestActor())
put!(x, ProxyMsg(;to=y,msg=StartMsg(0)))
put!(x, ProxyMsg(;to=y,msg=ReadStateMsg()))
@test take!(self()).state == 0

end


text:
using Documenter, AMA

makedocs(modules=[AMA],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/es335mathwiz/AMA.jl.git",
    julia  = "0.6.0",

    osname = "linux")


text:
using LightGraphs, MetaGraphs, Compose
using DataFrames, GraphPlot, Colors

"""
    read_graph

Reads a comma delimitted text file where each line represents: `source, destination, weight` for each edge.
"""
function read_graph(filename)
    f = open(filename)
    s = readlines(f)
    srcs = Vector{Int}(undef, length(s))
    dsts = Vector{Int}(undef, length(s))
    wgts = Vector{Float64}(undef, length(s))
    for (i, line) in enumerate(s)
        spl = split(line, ",")
        srcs[i] = parse(Int64, spl[1])
        dsts[i] = parse(Int64, spl[2])
        wgts[i] = parse(Float64, spl[3])
    end
    return construct_graph(srcs, dsts, wgts)
end

"""
    construct_graph(src, dst, wgt)

Takes in three vectors (sources, destinations, weights) and constructs a MetaGraph out of them.
"""
function construct_graph(srcs, dsts, wgts)
    # the maximal value in sources and destinations is the number of nodes
    n = maximum([srcs; dsts])
    G = MetaGraph(SimpleGraph(n))
    defaultweight!(G, 0.0)
    # For each src=>dst, wgt, create a weighted edge for it
    for (s,d,w) in zip(srcs, dsts, wgts)
        # ignore zero weight edges, since the default is 0 anyway
        if w > 0
            add_edge!(G, s, d)
            set_prop!(G, s, d, :weight, w)
        end
    end
    return G
end

"""
    plot_graph(G, fn)

Save the graph `G` as an SVG in the local directory with the filename `fn`.
Uses the default `spring_layout` layout.
"""
function plot_graph(G, fn)
    mn, mx = extrema(weights(G))
    line_colors = map(edges(G)) do ε
        src, dst = Tuple(ε)
        weighted_color_mean(weights(G)[src, dst], colorant"white", colorant"black")
    end

    nodelabel = 1:nv(G)
    p = gplot(G, edgestrokec = line_colors, nodelabel = nodelabel)
    draw(SVG(fn, 10cm, 10cm), p)
end

"""
    main

I don't believe in the concept of "main" functions, but whatever.
"""
function main()
    G = read_graph("weightsSmall.csv")
    plot_graph(G, "test-graph.svg")
end


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

@testset ExtendedTestSet "crop_foreground" begin
    @testset ExtendedTestSet "crop_foreground" begin
    x = [1 1 0 0; 1 1 0 0]
    x_crop = crop_foreground(x, x -> x > 0)
    @test x_crop == [1 1;1 1]
    end

    @testset ExtendedTestSet "crop_foreground" begin
    x = [1 1 0 0; 1 1 0 0]
    x_crop = crop_foreground(x, x -> x < 1)
    @test x_crop == [0 0; 0 0]
    end
end


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

@testset "Default objective sense" begin
    MOIT.default_objective_test(MOIU.MockOptimizer(MOIU.Model{Float64}()))
end

@testset "Default statuses" begin
    model = MOIU.MockOptimizer(MOIU.Model{Float64}())
    MOIT.default_status_test(model)
    MOI.empty!(model)
    MOIT.default_status_test(model)
end

@testset "Name test" begin
    MOIT.nametest(MOIU.MockOptimizer(MOIU.Model{Float64}()))
end

struct NoFreeModel <: MOI.ModelLike end
MOI.supports_add_constrained_variables(::NoFreeModel, ::Type{MOI.Reals}) = false

@testset "supports_add_constrained_variable" begin
    optimizer = MOIU.MockOptimizer(MOIU.Model{Float64}())
    @test MOI.supports_add_constrained_variable(
        optimizer,
        MOI.GreaterThan{Float64},
    )
    @test !MOI.supports_add_constrained_variable(
        optimizer,
        MOIT.UnknownScalarSet{Float64},
    )
    @test MOI.supports_add_constrained_variables(optimizer, MOI.Nonnegatives)
    @test !MOI.supports_add_constrained_variables(
        optimizer,
        MOIT.UnknownVectorSet,
    )

    nofree_optimizer = MOIU.MockOptimizer(NoFreeModel())
    @test !MOI.supports_add_constrained_variable(
        nofree_optimizer,
        MOI.GreaterThan{Float64},
    )
    @test !MOI.supports_add_constrained_variables(
        nofree_optimizer,
        MOI.Nonnegatives,
    )
    @test !MOI.supports_add_constrained_variables(nofree_optimizer, MOI.Reals)
end

@testset "Optimizer attributes" begin
    optimizer = MOIU.MockOptimizer(MOIU.Model{Float64}())
    @test MOI.supports(optimizer, MOIU.MockModelAttribute())
    MOI.set(optimizer, MOIU.MockModelAttribute(), 10)
    @test MOI.get(optimizer, MOIU.MockModelAttribute()) == 10

    v1 = MOI.add_variable(optimizer)
    @test MOI.supports(optimizer, MOIU.MockVariableAttribute(), typeof(v1))
    MOI.set(optimizer, MOIU.MockVariableAttribute(), v1, 11)
    @test MOI.get(optimizer, MOIU.MockVariableAttribute(), v1) == 11
    MOI.set(optimizer, MOIU.MockVariableAttribute(), [v1], [-11])
    @test MOI.get(optimizer, MOIU.MockVariableAttribute(), [v1]) == [-11]

    @test MOI.supports_constraint(
        optimizer,
        MOI.SingleVariable,
        MOI.GreaterThan{Float64},
    )
    c1 = MOI.add_constraint(
        optimizer,
        MOI.SingleVariable(v1),
        MOI.GreaterThan(1.0),
    )
    @test MOI.supports(optimizer, MOIU.MockConstraintAttribute(), typeof(c1))
    MOI.set(optimizer, MOIU.MockConstraintAttribute(), c1, 12)
    @test MOI.get(optimizer, MOIU.MockConstraintAttribute(), c1) == 12
    MOI.set(optimizer, MOIU.MockConstraintAttribute(), [c1], [-12])
    @test MOI.get(optimizer, MOIU.MockConstraintAttribute(), [c1]) == [-12]
end

@testset "Optimizer solve no result" begin
    optimizer = MOIU.MockOptimizer(MOIU.Model{Float64}())

    v1 = MOI.add_variable(optimizer)

    # Load fake solution
    MOI.set(optimizer, MOI.TerminationStatus(), MOI.INFEASIBLE)
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE

    MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.compute_conflict!(optimizer)
    @test MOI.get(optimizer, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
end

@testset "Optimizer solve with result" begin
    optimizer = MOIU.MockOptimizer(
        MOIU.Model{Float64}(),
        eval_objective_value = false,
        eval_variable_constraint_dual = false,
    )

    v = MOI.add_variables(optimizer, 2)
    c1 = MOI.add_constraint(
        optimizer,
        MOI.SingleVariable(v[1]),
        MOI.GreaterThan(1.0),
    )
    soc = MOI.add_constraint(
        optimizer,
        MOI.VectorOfVariables(v),
        MOI.SecondOrderCone(2),
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(v[1]),
    )
    MOI.set(optimizer, MOI.ResultCount(), 1)
    @test_throws(
        ErrorException("No mock primal is set for variable `$(v[1])`."),
        MOI.get(optimizer, MOI.VariablePrimal(), v[1])
    )
    @test_throws(
        MOI.InvalidIndex(MOI.VariableIndex(-1)),
        MOI.get(optimizer, MOI.VariablePrimal(), MOI.VariableIndex(-1))
    )
    MOI.set(optimizer, MOI.ResultCount(), 0)
    err = MOI.ResultIndexBoundsError(MOI.VariablePrimal(1), 0)
    @test_throws err MOI.get(optimizer, MOI.VariablePrimal(), v[1])
    err = MOI.ResultIndexBoundsError(MOI.ConstraintDual(1), 0)
    @test_throws err MOI.get(optimizer, MOI.ConstraintDual(), c1)

    # Load fake solution
    # TODO: Provide a more compact API for this.
    MOI.set(optimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(optimizer, MOI.ObjectiveValue(), 1.0)
    MOI.set(optimizer, MOI.DualObjectiveValue(2), 2.0)
    MOI.set(optimizer, MOI.ResultCount(), 2)
    MOI.set(optimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(optimizer, MOI.DualStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(optimizer, MOI.VariablePrimal(), v, [1.0, 2.0])
    MOI.set(optimizer, MOI.VariablePrimal(), v[1], 3.0)
    MOI.set(optimizer, MOI.ConstraintDual(2), c1, 5.9)
    MOI.set(optimizer, MOI.ConstraintDual(2), soc, [1.0, 2.0])

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.ResultCount()) == 2
    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 1.0
    @test isnan(MOI.get(optimizer, MOI.ObjectiveValue(2)))
    optimizer.eval_objective_value = true
    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 3.0
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.ObjectiveValue(2))
    )
    @test_throws(
        ErrorException(
            "No mock dual is set for constraint `$c1` at result index `1`.",
        ),
        MOI.get(optimizer, MOI.DualObjectiveValue())
    )
    @test MOI.get(optimizer, MOI.DualObjectiveValue(2)) == 5.9
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.DualStatus(2)) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.VariablePrimal(), v) == [3.0, 2.0]
    @test MOI.get(optimizer, MOI.VariablePrimal(), v[1]) == 3.0
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.VariablePrimal(2), v[1])
    )
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c1) == 3.0
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), soc) == [3.0, 2.0]
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        @show MOI.get(optimizer, MOI.ConstraintPrimal(2), c1)
    )
    @test_throws(
        ErrorException(
            "No mock primal is set for variable `$(v[1])` at result index `2`.",
        ),
        MOI.get(optimizer, MOI.ConstraintPrimal(2), soc)
    )
    @test MOI.get(optimizer, MOI.ConstraintDual(2), c1) == 5.9
    @test MOI.get(optimizer, MOI.ConstraintDual(2), soc) == [1.0, 2.0]
    @test_throws(
        ErrorException(
            "No mock dual is set for constraint `$c1` at result index `1`.",
        ),
        MOI.get(optimizer, MOI.ConstraintDual(1), c1)
    )
end

@testset "Delete" begin
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
    MOIT.delete_test(mock)
end

@testset "CanonicalConstraintFunction" begin
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    fx, fy = MOI.SingleVariable.(MOI.add_variables(mock, 2))
    cx = MOI.add_constraint(mock, fx, MOI.LessThan(0))
    c = MOI.add_constraint(mock, 1fx + fy, MOI.LessThan(1))
    @test MOIU.is_canonical(MOI.get(mock, MOI.ConstraintFunction(), cx))
    @test MOIU.is_canonical(
        MOI.get(mock, MOI.CanonicalConstraintFunction(), cx),
    )
    @test !MOIU.is_canonical(MOI.get(mock, MOI.ConstraintFunction(), c))
    @test MOIU.is_canonical(MOI.get(mock, MOI.CanonicalConstraintFunction(), c))
end

@testset "Conflict access" begin
    mock = MOIU.MockOptimizer(MOIU.Model{Int}())
    fx, fy = MOI.SingleVariable.(MOI.add_variables(mock, 2))
    cx = MOI.add_constraint(mock, fx, MOI.LessThan(0))
    c = MOI.add_constraint(mock, 1fx + fy, MOI.LessThan(1))
    MOI.set(mock, MOI.ConstraintConflictStatus(), cx, MOI.NOT_IN_CONFLICT)
    MOI.set(mock, MOI.ConstraintConflictStatus(), c, MOI.IN_CONFLICT)
    MOI.compute_conflict!(mock)

    @test MOI.get(mock, MOI.ConstraintConflictStatus(), cx) ==
          MOI.NOT_IN_CONFLICT
    @test MOI.get(mock, MOI.ConstraintConflictStatus(), c) == MOI.IN_CONFLICT
end