text:
# Included from build.jl, ../test/test_build.jl and ../src/PyCall.jl

import Libdl

hassym(lib, sym) = Libdl.dlsym_e(lib, sym) != C_NULL

# call dlsym_e on a sequence of symbols and return the symbol that gives
# the first non-null result
function findsym(lib, syms...)
    for sym in syms
        if hassym(lib, sym)
            return sym
        end
    end
    error("no symbol found from: ", syms)
end

# Static buffer to make sure the string passed to libpython persists
# for the lifetime of the program, as CPython API requires:
const __buf_programname = Vector{UInt8}(undef, 1024)
const __buf_pythonhome = Vector{UInt8}(undef, 1024)

# Need to set PythonHome before calling GetVersion to avoid warning (#299).
# Unfortunately, this poses something of a chicken-and-egg problem because
# we need to know the Python version to set PythonHome via the API.  Note
# that the string (or array) passed to Py_SetPythonHome needs to be a
# constant that lasts for the lifetime of the program, which is why we
# prepare static buffer __buf_pythonhome, copy the string to it, and then
# pass the pointer to the buffer to the CPython API.
function Py_SetPythonHome(libpy, pyversion, PYTHONHOME::AbstractString)
    isempty(PYTHONHOME) && return
    if pyversion.major < 3
        ccall(Libdl.dlsym(libpy, :Py_SetPythonHome), Cvoid, (Cstring,),
              _preserveas!(__buf_pythonhome, Cstring, PYTHONHOME))
    else
        ccall(Libdl.dlsym(libpy, :Py_SetPythonHome), Cvoid, (Ptr{Cwchar_t},),
              _preserveas!(__buf_pythonhome, Cwstring, PYTHONHOME))
    end
end

function Py_SetProgramName(libpy, pyversion, programname::AbstractString)
    isempty(programname) && return
    if pyversion.major < 3
        ccall(Libdl.dlsym(libpy, :Py_SetProgramName), Cvoid, (Cstring,),
              _preserveas!(__buf_programname, Cstring, programname))
    else
        ccall(Libdl.dlsym(libpy, :Py_SetProgramName), Cvoid, (Ptr{Cwchar_t},),
              _preserveas!(__buf_programname, Cwstring, programname))
    end
end

"""
    _preserveas!(dest::Vector{UInt8}, (Cstring|Cwstring), x::String) :: Ptr

Copy `x` as `Cstring` or `Cwstring` to `dest` and return a pointer to
`dest`.  Thus, this pointer is safe to use as long as `dest` is
protected from GC.
"""
function _preserveas!(dest::Vector{UInt8}, ::Type{Cstring}, x::AbstractString)
    s = transcode(UInt8, String(x))
    copyto!(dest, s)
    dest[length(s) + 1] = 0
    return pointer(dest)
end

function _preserveas!(dest::Vector{UInt8}, ::Type{Cwstring}, x::AbstractString)
    s = Base.cconvert(Cwstring, x)
    copyto!(dest, reinterpret(UInt8, s))
    return pointer(dest)
end

# need to be able to get the version before Python is initialized
Py_GetVersion(libpy) = unsafe_string(ccall(Libdl.dlsym(libpy, :Py_GetVersion), Ptr{UInt8}, ()))

# Fix the environment for running `python`, and setts IO encoding to UTF-8.
# If cmd is the Conda python, then additionally removes all PYTHON* and
# CONDA* environment variables.
function pythonenv(cmd::Cmd)
    env = copy(ENV)
    if dirname(cmd.exec[1]) == abspath(Conda.PYTHONDIR)
        pythonvars = String[]
        for var in keys(env)
            if startswith(var, "CONDA") || startswith(var, "PYTHON")
                push!(pythonvars, var)
            end
        end
        for var in pythonvars
            pop!(env, var)
        end
    end
    # set PYTHONIOENCODING when running python executable, so that
    # we get UTF-8 encoded text as output (this is not the default on Windows).
    env["PYTHONIOENCODING"] = "UTF-8"
    setenv(cmd, env)
end

function pythonhome_of(pyprogramname::AbstractString)
    if Sys.iswindows()
        # PYTHONHOME tells python where to look for both pure python
        # and binary modules.  When it is set, it replaces both
        # `prefix` and `exec_prefix` and we thus need to set it to
        # both in case they differ. This is also what the
        # documentation recommends.  However, they are documented
        # to always be the same on Windows, where it causes
        # problems if we try to include both.
        script = """
        import sys
        if hasattr(sys, "base_exec_prefix"):
            sys.stdout.write(sys.base_exec_prefix)
        else:
            sys.stdout.write(sys.exec_prefix)
        """
    else
        script = """
        import sys
        if hasattr(sys, "base_exec_prefix"):
            sys.stdout.write(sys.base_prefix)
            sys.stdout.write(":")
            sys.stdout.write(sys.base_exec_prefix)
        else:
            sys.stdout.write(sys.prefix)
            sys.stdout.write(":")
            sys.stdout.write(sys.exec_prefix)
        """
        # https://docs.python.org/3/using/cmdline.html#envvar-PYTHONHOME
    end
    return read(pythonenv(`$pyprogramname -c $script`), String)
end
# To support `venv` standard library (as well as `virtualenv`), we
# need to use `sys.base_prefix` and `sys.base_exec_prefix` here.
# Otherwise, initializing Python in `__init__` below fails with
# unrecoverable error:
#
#   Fatal Python error: initfsencoding: unable to load the file system codec
#   ModuleNotFoundError: No module named 'encodings'
#
# This is because `venv` does not symlink standard libraries like
# `virtualenv`.  For example, `lib/python3.X/encodings` does not
# exist.  Rather, `venv` relies on the behavior of Python runtime:
#
#   If a file named "pyvenv.cfg" exists one directory above
#   sys.executable, sys.prefix and sys.exec_prefix are set to that
#   directory and it is also checked for site-packages
#   --- https://docs.python.org/3/library/venv.html
#
# Thus, we need point `PYTHONHOME` to `sys.base_prefix` and
# `sys.base_exec_prefix`.  If the virtual environment is created by
# `virtualenv`, those `sys.base_*` paths point to the virtual
# environment.  Thus, above code supports both use cases.
#
# See also:
# * https://docs.python.org/3/library/venv.html
# * https://docs.python.org/3/library/site.html
# * https://docs.python.org/3/library/sys.html#sys.base_exec_prefix
# * https://github.com/JuliaPy/PyCall.jl/issues/410


text:
#
using Flux, NNlib, DiffEqFlux
using Plots
import SEM
#--------------------------------------#
d   = 2
Lb2 = 1e1
np  = 100
nl  = 10
#--------------------------------------#
xx = SEM.linspace(-Lb2,Lb2,np)

if(d==1)
    x = SEM.ndgrid(xx)
    xy = Array{Float64}(undef,np,d)
    xy .= x
elseif(d==2)
    x,y = SEM.ndgrid(xx,xx)
    xy = Array{Float64}(undef,np,np,d)
    xy[:,:,1] .= x
    xy[:,:,2] .= y
end

xy = flatten(xy)'

NN = Chain(Dense(rand(nl, d),rand(nl),relu)
          ,Dense(rand(nl,nl),rand(nl),relu)
          ,Dense(rand(nl,nl),rand(nl),relu)
          ,Dense(rand(nl,nl),rand(nl),relu)
          ,Dense(rand(nl,nl),rand(nl),relu)
          ,Dense(rand( 1,nl),rand( 1)))

f = NN(xy)
f = reshape(f,size(x))
#--------------------------------------#
# partitions
if(d==1)
    # first layer discontinuitites
    xpart = - NN[1].b ./ NN[1].W
    fpart = NN(xpart')'
elseif(d==2)
    # first layer discontinuities
    #NN[1].W * [x,y] + NN[1].b == 0
    x1 = -Lb2 * ones(size(NN[1].b))
    x2 =  Lb2 * ones(size(NN[1].b))
    y1 = @. (-NN[1].b - x1*NN[1].W[:,1]) / NN[1].W[:,2]
    y2 = @. (-NN[1].b - x2*NN[1].W[:,1]) / NN[1].W[:,2]

    xpart = vcat(x1',x2')
    fpart = vcat(y1',y2')
end
#--------------------------------------#
plt = plot()
if(d==1)
    plt = plot!(x,f)
    plt = scatter!(xpart,fpart,color="red")
    plt = plot!(title="f=NN(x)",xlabel="x",ylabel="y=NN(x)")
    plt = plot!(xlims=(-Lb2,Lb2))
elseif(d==2)
    plt = plot!(xpart,fpart,legend=false,color="black")
    plt = plot!(title="Shallow ReLU DNN Partitions",xlabel="x",ylabel="y")
#   plt = plot!(xlims=(-Lb2,Lb2),ylims=(-Lb2,Lb2))
end

plt
#savefig(plt,"reluPartition.png")
#display(plt)

#plt1 = plot()
#if(d==2)
#    plt1 = plot()
#    plt1 = heatmap!(f')
#end
#
#savefig(plt1,"reluPartitionHmp.png")
#display(plt1)
#

ext: jl


text:
function center(U, mass)
    d = length(U[1].x)

    c = zeros(Float64, d)

    for i = 1:length(mass)
        c += mass[i] .* U[i].x
    end

    return c / sum(mass)
end

function center(U::Array)
    n = length(U)

    mass = getMass(U)

    return center(U, mass),
    argworst(U), # worst
    argbest(U)  # best
end

function getU(P::Array, K::Int, I::Vector{Int}, i::Int, N::Int)
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    return P[U_ids]
end

function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int)
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    return U_ids
end

function getU_ids(K::Int, I::Vector{Int}, i::Int, N::Int, feasible_solutions)
    # at least half of the population is feasible to generate random centers
    if length(feasible_solutions) >= 0.5N || length(feasible_solutions) == 0
        return getU_ids(K, I, i, N)
    end

    # center of mass is generated with at least one feasible solution
    K -= 1
    if i <= N - K
        U_ids = I[i:K+i]
    else
        j = (i:K+i) .% N
        U_ids = I[j.+1]
    end

    push!(U_ids, rand(feasible_solutions))

    return U_ids
end

function crossover(
    x::Vector{Float64},
    y::Vector{Float64},
    p_cr::Vector{Float64},
)
    D = length(x)
    tmp2 = zeros(D)
    for j = 1:D
        if rand() < p_cr[j]
            y[j] = x[j]
            tmp2[j] += 1
        end
    end

    return y, tmp2
end


text:
# Transient CellField
abstract type TransientCellField <: CellField end

get_data(f::TransientCellField) = @abstractmethod
get_triangulation(f::TransientCellField) = @abstractmethod
DomainStyle(::Type{TransientCellField}) =  @abstractmethod
gradient(f::TransientCellField) =  @abstractmethod
∇∇(f::TransientCellField) =  @abstractmethod
function change_domain(f::TransientCellField,trian::Triangulation,target_domain::DomainStyle)
  @abstractmethod
end

struct TransientSingleFieldCellField{A} <: TransientCellField
  cellfield::A
  derivatives::Tuple#{Vararg{A,B} where B}
end

SingleFieldTypes = Union{GenericCellField,SingleFieldFEFunction}

function TransientCellField(single_field::SingleFieldTypes,derivatives::Tuple)
  TransientSingleFieldCellField(single_field,derivatives)
end

# CellField methods
get_data(f::TransientSingleFieldCellField) = get_data(f.cellfield)
get_triangulation(f::TransientSingleFieldCellField) = get_triangulation(f.cellfield)
DomainStyle(::Type{<:TransientSingleFieldCellField{A}}) where A = DomainStyle(A)
gradient(f::TransientSingleFieldCellField) = gradient(f.cellfield)
∇∇(f::TransientSingleFieldCellField) = ∇∇(f.cellfield)
change_domain(f::TransientSingleFieldCellField,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.cellfield,trian,target_domain)

# Skeleton related Operations
function Base.getproperty(f::TransientSingleFieldCellField, sym::Symbol)
  if sym in (:⁺,:plus,:⁻, :minus)
    derivatives = ()
    if sym in (:⁺,:plus)
      cellfield = CellFieldAt{:plus}(f.cellfield)
      for iderivative in f.derivatives
        derivatives = (derivatives...,CellFieldAt{:plus}(iderivative))
      end
    elseif sym in (:⁻, :minus)
      cellfield = CellFieldAt{:minus}(f.cellfield)
      for iderivative in f.derivatives
        derivatives = (derivatives...,CellFieldAt{:plus}(iderivative))
      end
    end
    return TransientSingleFieldCellField(cellfield,derivatives)
  else
    return getfield(f, sym)
  end
end

# Transient FEBasis
struct TransientFEBasis{A} <: FEBasis
  febasis::A
  derivatives::Tuple{Vararg{A}}
end

# FEBasis methods
get_data(f::TransientFEBasis) = get_data(f.febasis)
get_triangulation(f::TransientFEBasis) = get_triangulation(f.febasis)
DomainStyle(::Type{<:TransientFEBasis{A}}) where A = DomainStyle(A)
BasisStyle(::Type{<:TransientFEBasis{A}}) where A = BasisStyle(A)
gradient(f::TransientFEBasis) = gradient(f.febasis)
∇∇(f::TransientFEBasis) = ∇∇(f.febasis)
change_domain(f::TransientFEBasis,trian::Triangulation,target_domain::DomainStyle) = change_domain(f.febasis,trian,target_domain)

# Time derivative
function ∂t(f::Union{TransientCellField,TransientFEBasis})
  cellfield, derivatives = first_and_tail(f.derivatives)
  TransientCellField(cellfield,derivatives)
end

∂tt(f::Union{TransientCellField,TransientFEBasis}) = ∂t(∂t(f::Union{TransientCellField,TransientFEBasis}))


text:
@testset "basic" begin
    let A=rand(10)
        s1 = VectorShape(Float64, 10)
        s2 = VectorShape(Float64, 11)
        @test_throws DimensionMismatch ShapedView(A, -1, s1)
        @test_throws DimensionMismatch ShapedView(A, 11, s1)
        @test_throws DimensionMismatch ShapedView(A, 0, s2)
    end

    ms = MultiShape(
        x = MatrixShape(Int,5,5),
        y = ScalarShape(Float64),
        z = VectorShape(Float64, 2),
    )
    lx = length(ms.x)
    ly = length(ms.y)
    lz = length(ms.z)
    d = Float64.(collect(1:(lx + ly + lz)))
    sv = ms(d)

    x(sv) = sv.x
    y(sv) = sv.y
    z(sv) = sv.z
    @test @inferred(x(sv)) == reshape(view(d, 1:lx), size(ms.x))
    @test @inferred(y(sv)) === d[1+lx]
    @test @inferred(z(sv)) == view(d, (1 + lx + ly):(lx + ly + lz))

    d1 = rand(1:100, 5, 5)
    d2 = rand()
    d3 = rand(2)

    sv.x .= d1
    sv.y = d2
    sv.z .= d3

    @test sv.x == d1
    @test sv.y == d2
    @test sv.z == d3
end

(s1=rand(5,10), s2=rand(5), )

@testset "nested" begin
    s1 = MatrixShape(Real,5,10)
    s2 = VectorShape(Float32,5)
    s3 = ScalarShape(Int)
    ms1 = MultiShape(s1 = s1, s2 = s2, s3=s3)

    s4 = VectorShape(Float32,3)
    s5 = ScalarShape(Int)
    ms2 = MultiShape(s4=s4, s5=s5)

    s6 = ScalarShape(Float64)
    ms3 = MultiShape(s6=s6, ms2=ms2)

    ms = MultiShape(ms1=ms1, ms2=ms2, ms3=ms3)
    x = collect(1:length(ms))
    sv = ms(x)

    @test sv.ms1 == 1:56
    @test vec(sv.ms1.s1) == 1:50
    @test vec(sv.ms1.s2) == 51:55
    @test sv.ms1.s3 == 56

    @test vec(sv.ms2.s4) == 57:59
    @test sv.ms2.s5 == 60

    @test sv.ms3.s6 == 61
    @test vec(sv.ms3.ms2.s4) == 62:64
    @test sv.ms3.ms2.s5 == 65
end


text:
"""
    half_saturation(max, half, x)
Half satration curve.
"""
half_saturation(max, half, x) = max/(oneunit(half/x) + half/x)


text:
"""
@MStructType name fieldnames...

This macro gives a convenient syntax for declaring mutable `StructType`s for reading specific variables from a JSONLines file. Also defines `row[:col]` access for rows of the resulting type.

* `name`: Name of the `StructType`
* `fieldnames...`: Names of the variables to be read (must be the same as in the file)
"""
macro MStructType(name, fieldnames...)
quote
    mutable struct $name
        $(fieldnames...)
        $(esc(name))() = new(fill(missing, $(esc(length(fieldnames))))...)
    end
    StructTypes.StructType(::Type{$(esc(name))}) = StructTypes.Mutable()
    StructTypes.names(::Type{$(esc(name))}) = tuple()
    Base.getindex(x::$(esc(name)), col::Symbol) = getproperty(x, col)
end
end

function MStructType(name, vars...)
eval(:(@MStructType $name $(vars...)))
end

macro readcols(path::String, nworkers, cols...)
    name = gensym(basename(path))
    MStructType(name, cols...)
    quote
       LineIndex($path, structtype = $name, nworkers = $nworkers)
    end
end

"""
    readcols(path::String, cols...; nworkers = 1) => LineIndex

* `path`: Path to JSONLines file
* `cols...`: Columnnames to be selected
* Keyword Argument:
    * `nworkers=1`: Number of threads to use for operations on the resulting LineIndex
"""
function readcols(path::String, cols...; nworkers = 1)
    eval(:(@readcols $path $nworkers $(cols...)))
end


text:
using DFTK
using Test
include("testcases.jl")

@testset "Symmetrization and not symmetrization yields the same density and energy" begin
    args = ((kgrid=[2, 2, 2], kshift=[1/2, 0, 0]),
            (kgrid=[2, 2, 2], kshift=[1/2, 1/2, 0]),
            (kgrid=[2, 2, 2], kshift=[0, 0, 0]),
            (kgrid=[3, 2, 3], kshift=[0, 0, 0]),
            (kgrid=[3, 2, 3], kshift=[0, 1/2, 1/2]))
    testcase = silicon
    for case in args
        model_nosym = model_LDA(testcase.lattice, testcase.atoms, testcase.positions;
                                symmetries=false)
        basis = PlaneWaveBasis(model_nosym; Ecut=5, case...)
        DFTK.check_group(basis.symmetries)

        scfres = self_consistent_field(basis; is_converged=DFTK.ScfConvergenceDensity(1e-10))
        ρ1 = scfres.ρ
        E1 = scfres.energies.total

        model_sym = model_LDA(testcase.lattice, testcase.atoms, testcase.positions)
        basis = PlaneWaveBasis(model_sym; Ecut=5, case...)
        scfres = self_consistent_field(basis; is_converged=DFTK.ScfConvergenceDensity(1e-10))
        ρ2 = scfres.ρ
        E2 = scfres.energies.total

        @test abs(E1 - E2) < 1e-10
        @test norm(ρ1 - ρ2) .* sqrt(basis.dvol) < 1e-8
    end
end


text:
# Shark: Chess Utilities
@inline color(p::UInt8)::UInt8 = p & -p;
@inline class(p::UInt8)::UInt8 = p & OFF;
@inline opposite(p::UInt8)::UInt8 = p ⊻ 3;

# Piece "Constructor"
function piece(x::Char = '.')::UInt8
    color = islowercase(x) ? BLACK : WHITE
    class = lowercase(x)
    if class == 'p'
        return color | PAWN;
    elseif class == 'n'
        return color | KNIGHT;
    elseif class == 'b'
        return color | BISHOP;
    elseif class == 'r'
        return color | ROOK;
    elseif class == 'q'
        return color | QUEEN;
    elseif class == 'k'
        return color | KING;
    elseif class == '.'
        return NONE;
    else # class == ' '
        return OFF;
    end
end

max_stars_repo_name: Joshua-Robison/Shark

max_issues_repo_name: Joshua-Robison/Shark

max_forks_repo_name: Joshua-Robison/Shark


text:
module LaplaceCoefficients

VERSION < v"0.4-dev" && using Docile

export laplace_coefficients_initialize
export laplace_wisdom
include("laplace_wisdom.jl")

"""
# This computes the Laplace coefficients via recursion.
"""
function initialize(jmax::Integer,alpha::Number)
const nmax=7
b0=Array(eltype(alpha),nmax,jmax+1) # Array to hold the coefficients
# Compute the highest two Laplace coefficients using Wisdom's series approach:
for j=0:jmax
  for i=0:nmax-1
    b0[i+1,j+1]=laplace_wisdom(1//2,j,i,alpha)/alpha^i
  end
end
return b0
end

end # module

max_stars_repo_name: ericagol/TTVFaster

max_issues_repo_name: ericagol/TTVFaster

max_forks_repo_name: ericagol/TTVFaster


text:
using Plots
using Reproduce
using JLD2
using FileIO

function plot_single_is_ratio(results, log_interval)
    previous_GKSwstype = get(ENV, "GKSwstype", "")
    ENV["GKSwstype"] = "100"

    is_ratio = results[:is_ratio]

    println(size(is_ratio))
    # println(size(emphasis))
    # println(emphasis[:, 1:4])
    # println(maximum(emphasis))
    # println(minimum(emphasis))
    # println(mean(emphasis))
    # println(median(emphasis))

    # println(is_ratio[:,])

    p = plot()

    for gvf_i in 1:4
        plot!(p, 1:log_interval:size(is_ratio)[2] * log_interval, is_ratio[gvf_i, :])
    end
    save_file = "plotting/chunlok/generated_plots/single_is_ratio.svg"
    savefig(save_file)

    ylim = (0, 2.0)

    p = plot()
    plot!(p, 1:log_interval:size(is_ratio)[2] * log_interval, is_ratio[1, :], ylim=ylim, title="gvf 1")
    save_file = "plotting/chunlok/generated_plots/single_is_ratio_gvf_1.svg"
    savefig(save_file)

    p = plot()
    plot!(p, 1:log_interval:size(is_ratio)[2] * log_interval, is_ratio[2, :], ylim=ylim, title="gvf 2")
    save_file = "plotting/chunlok/generated_plots/single_is_ratio_gvf_2.svg"
    savefig(save_file)
    p = plot()
    plot!(p, 1:log_interval:size(is_ratio)[2] * log_interval, is_ratio[3, :], ylim=ylim, title="gvf 3")
    save_file = "plotting/chunlok/generated_plots/single_is_ratio_gvf_3.svg"
    savefig(save_file)
    p = plot()
    plot!(p, 1:log_interval:size(is_ratio)[2] * log_interval, is_ratio[4, :], ylim=ylim, title="gvf 4")
    save_file = "plotting/chunlok/generated_plots/single_is_ratio_gvf_4.svg"
    savefig(save_file)
    # asdasds

    # println(findall(x->x!=1.0, step_sizes))
    # dasd
    # println(step_sizes[:, :, 1])
    # println(step_sizes[:, :, 1] - step_sizes[:, :, 100])

    # println(minimum(step_sizes))
    # println(maximum(step_sizes))

    # function plot_gvf_step_sizes(range)
    #     gvf_stepsizes = minimum(step_sizes[range, :, :], dims=(1, 2))[1, 1, :]
    #     plot!(p, 1:log_interval:size(gvf_stepsizes)[1] * log_interval, gvf_stepsizes)
    # end

    # plot_gvf_step_sizes(1:4)
    # plot_gvf_step_sizes(5:8)
    # plot_gvf_step_sizes(9:12)
    # plot_gvf_step_sizes(13:16)

    # for i in 1:4
    #     # println(size(error[i, :]))
    #     plot!(p, 1:size(error)[2], error[i, :])
    # end

    # # print(error[3,10])

    println("Single autostep stepsize plot saved to $(save_file)")

    ENV["GKSwstype"] = previous_GKSwstype
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module PlottingUtils

using Gadfly
using ColorSchemes
using Measures
using DataFrames
using Statistics: mean, quantile
import Cairo,
      Fontconfig

const PLOT_FONT = "Arial"
const VTR = 0.125

"""
Adds noise to a Vector of Vector
Used for adding simulated noise to
the states of a system.
Assumes that the values can be negative.
"""
function add_noise(data::Vector{Vector{T}}, percent::T,
                   positive::Bool=false)::Vector{Vector{T}} where T
    if percent == 0.0
        return data
    end

    noise_data = deepcopy(data)
    epsilon_arr = [0.01*mean(getindex.(data,i)) for i in 1:length(data[1])]
    sigma = zero(T)
    @inbounds for i in 1:length(noise_data)
        for j in 1:length(noise_data[1])
            sigma = abs(percent*noise_data[i][j] + epsilon_arr[j])
            d = Normal(0,sigma)
            noise_data[i][j] += rand(d)
            noise_data[i][j] = positive ? abs(noise_data[i][j]) : noise_data[i][j]
        end
    end
    return noise_data
end

function step_success_rate(x::T)::Int64 where T
    if x < VTR
       return 1
    else
        return 0
    end
end

function box_data(arr::Array{<:AbstractFloat})::Array{<:AbstractFloat}
    q = quantile(arr,[0.25,0.5,0.75])
    iqr = q[3] - q[1]
    min_val = q[1] - 1.5*iqr
    max_val = q[3] + 1.5*iqr
    return [maximum([min_val,minimum(arr)]),
            q[1],q[2],q[3],
            minimum([max_val,maximum(arr)])]
end

function box_scatter_plot(arr::Array{<:AbstractFloat})::Array{<:AbstractFloat}
    qarr = box_data(arr)
    return [qarr[1],qarr[3],qarr[5]]
end

function get_plot_data(results::Dict,
                        noise_level::AbstractArray{Symbol},
                        methods::Vector{Symbol})::Dict
    plot_data = Dict()
    for m in methods
        plot_data[m] = Dict()
        plot_data[m]["error"] = []
        plot_data[m]["time"] = []
        plot_data[m]["est"] = []
    end
    for v in noise_level
        for m in methods
            error = [e[1][1] for e in results[v][m]]
            if length(error) > 0
                push!(plot_data[m]["error"], error)
            else
                push!(plot_data[m]["error"], [NaN])
            end
            time = [e[2][1] for e in results[v][m]]
            if length(time) > 0
                push!(plot_data[m]["time"], time)
            else
                push!(plot_data[m]["time"], [NaN])
            end
            est = [e[3] for e in results[v][m]]
            if length(time) > 0
                push!(plot_data[m]["est"], est)
            else
                push!(plot_data[m]["est"], [NaN])
            end
        end
    end
    return plot_data
end

"""
Error Plots
xaxis: Noise Percentage
yaxis: Mean Error
"""
function error_plots(plot_data::Dict,
                    noise_level::AbstractArray{<:AbstractFloat},
                    methods::Array,
                    method_label::Dict,
                    method_color::Dict,
                    sam::Int,
                    path::String)::String

    nrmse_log = "NRMSE\n"
    p = plot(x=noise_level, y=[],Geom.line,
            Theme(background_color=colorant"white",
                panel_fill=colorant"white",
                major_label_font=PLOT_FONT,
                minor_label_font=PLOT_FONT))

    p = plot(x=noise_level, xlabel="Noise Percentage", ylabel="NRMSE", legend=:outertopright)
    ylim_arr = []
    for m in methods
        nrmse_log *= "$(m), "
        p2 = plot(x=noise_level, xlabel="Noise Percentage", ylabel="NRMSE", legend=:outertopright)
        error = plot_data[m]["error"]
        # Proceed only if there are no NaN
        if !any(isnan.(vcat(error...)))
            qerror = hcat(box_scatter_plot.(error)...)
            # Don't append to plot p if data is too far
            # from the expected values (i.e. plot only if
            # data < 5)
            if !any(qerror[3,:] .> 5.0)

                append!(p.layers,layer(x=x1, y=y1,
                                       size=map(x->x*mm, repeat([1.25], length(x1))),
                                       Theme(default_color=colorant"#e6612f"), #Vivid Vermilion
                                       Geom.PointGeometry))

#                plot!(p, noise_level, qerror[2,:],
#                            grid=true,
#                            ribbon=(qerror[1,:],
#                                    qerror[3,:]),
#                            fillalpha=.5, label=method_label[m], color=method_color[m])
                push!(ylim_arr, ylims(p))
            end

            plot!(p2, noise_level, qerror[2,:],
                        grid=true,
                        ribbon=(qerror[1,:],
                                qerror[3,:]),
                        fillalpha=.5, label=method_label[m], color=method_color[m])
            #display(p2)
            nrmse_log *= "$(qerror)\n"

            savefig(p,path*"/error_inter_$(m)_$(sam).pdf")
            savefig(p2,path*"/error_$(m)_$(sam).pdf")
        end
    end
    #=>
    if length(ylim_arr) > 1
        ylims!(p, (minimum(abs.([yl[1] for yl in ylim_arr])),
                    minimum(abs.([yl[2] for yl in ylim_arr]))))
    end
    <=#
    #display(p)

    savefig(p,path*"/error_all_$(sam).pdf")
    return nrmse_log
end

"""
Box Error Plots
xaxis: Method
yaxis: Error Distribution
"""
function box_error_plots(plot_data::Dict,
                    var::Symbol,
                    methods::Vector{Symbol},
                    method_label::Dict,
                    method_color::Dict,
                    sam::Symbol,
                    path::String)::String

    log_str = "NRMSE\n"
    df = DataFrame()
    for m in methods
        # Substitute infinite for missing
        # so that boxplot can still work
        # with the data
        log_str *= "$(m), "
        data = [ifelse(isinf(x),missing,x) for x in plot_data[m]["error"][1]]
        aux_df = DataFrame(method=method_label[m],
                    data=data, color=method_color[m])
        log_str *= "$(round(box_scatter_plot(data)[2], digits=4))\n "
        try
            append!(df, aux_df)
        catch e
            @show e
        end
    end
    p = plot(df, x="method", y="data", color="color", Geom.boxplot,
                Guide.xlabel("Method"), Guide.ylabel("NRMSE"),
                Theme(background_color=colorant"white",
                    panel_fill=colorant"white",
                    major_label_font=PLOT_FONT,
                    minor_label_font=PLOT_FONT))
    try
        p |> PDF(joinpath(path,replace("box_$(sam)_$(var)","."=>"")*".pdf"))
    catch e
        @show e
        println(df)
        println()
    end
    return log_str
end

"""
Parameter Distribution Plots
xaxis: Parameter
yaxis: Value Distribution
"""
function parameter_plots(plot_data::Dict,
                    var::Symbol,
                    methods::Array,
                    method_label::Dict,
                    method_color::Dict,
                    sam::Symbol,
                    path::String)::Nothing
    num_pars = length(plot_data[methods[1]]["est"][1][1])
    p_arr = []
    for m in methods
        p = plot(legend=:outertopright, ylabel="Value", xlabel="Parameter")
        for i in 1:num_pars
            data = getindex.(plot_data[m]["est"][1],i)
            if i == 1
                boxplot!(p, [string(i)], log10.(data), color=method_color[m], label=m)
            else
                boxplot!(p, [string(i)], log10.(data), color=method_color[m], label="")
            end
        end
        savefig(p,joinpath(path,replace("par_$(m)_$(sam)_$(var)","."=>"")*".pdf"))
        push!(p_arr,p)
    end

    p = plot(p_arr...,layout=(length(methods),1))
    savefig(p,joinpath(path,replace("par_all_$(sam)_$(var)","."=>"")*".pdf"))
end

"""
Success Rate vs Time Plots
xaxis: Mean Computation Time
yaxis: 1 / Success Rate
"""
function sr_plots(plot_data::Dict,
                    noise_level::AbstractArray{Symbol},
                    methods::Vector{Symbol},
                    method_label::Dict,
                    method_color::Dict,
                    sam::Symbol,
                    path::String)::String

    log_str = "Method, SR, Cost, ISR\n"

    p = plot(x=[], y=[], Geom.point,
                Guide.xlabel("Time"), Guide.ylabel("1 / Success Rate"),
                Theme(background_color=colorant"white",
                    panel_fill=colorant"white",
                    major_label_font=PLOT_FONT,
                    minor_label_font=PLOT_FONT),
                    Guide.manual_color_key("Method",
                    [l[end] for l in method_label],
                    [c[end] for c in method_color]))

    p2 = plot(x=[], y=[], Geom.point,
                Guide.xlabel("Time"), Guide.ylabel("1 / Success Rate"),
                Theme(background_color=colorant"white",
                    panel_fill=colorant"white",
                    major_label_font=PLOT_FONT,
                    minor_label_font=PLOT_FONT),
                    Guide.manual_color_key("Method",
                    [l[end] for l in method_label],
                    [c[end] for c in method_color]))

    ylim_arr = []
    for m in methods
        log_str *= "$(m), "
        p3 = plot(x=[], y=[],Geom.point,
                Guide.xlabel("Time"), Guide.ylabel("1 / Success Rate"),
                Theme(background_color=colorant"white",
                    panel_fill=colorant"white",
                    major_label_font=PLOT_FONT,
                    minor_label_font=PLOT_FONT))

        sr = mean.([step_success_rate.(e) for e in plot_data[m]["error"]])
        qsr = box_scatter_plot(sr)
        log_str *= "$(round(qsr[2], digits=4)), "
        isr = 1.0./sr
        timed = plot_data[m]["time"]
        if !any(isnan.(vcat(timed...)))
            qtime = hcat(box_scatter_plot.(timed)...)
            qqtime = box_scatter_plot(qtime[2,:])
            qisr = box_scatter_plot(isr)
            log_str *= "$(round(qqtime[2], digits=4)), "
            log_str *= "$(round(qisr[2], digits=4))\n"
            if qisr[2] <= 5.0
                append!(p.layers,
                    layer(x=[qqtime[2]], y=[qisr[2]],
                            xmin=[qqtime[1]],
                            xmax=[qqtime[3]],
                            ymin=[qisr[1]],
                            ymax=[qisr[3]],
                            color=[method_color[m]],
                            Geom.point, Geom.errorbar))
            end

            append!(p2.layers, layer(x=qtime[2,:], y=isr,
                                     size=map(x->x*mm,
                                              repeat([1.25], length(qtime[2,:]))),
                                     color=[method_color[m]], Geom.point))

            append!(p3.layers, layer(x=qtime[2,:], y=isr,
                                     color=[method_color[m]],
                                     label=String.(noise_level),
                                     Geom.point, Geom.label))

            try
                p3 |> PDF(joinpath(path,"sr_$(m)_$(sam).pdf"))
            catch e
                @show e
            end
        end
    end
    # display(p2)
    # display(p)

    try
        p2 |> PDF(joinpath(path,"sr_all_$(sam).pdf"))
    catch e
        @show e
    end

    try
        p |> PDF(joinpath(path,"sr_all_medians_$(sam).pdf"))
    catch e
        @show e
    end

    log_str *= "\n"
    return log_str
end

function sr_plots(plot_data::Dict,
                    var::Symbol,
                    methods::Vector{Symbol},
                    method_label::Dict,
                    method_color::Dict,
                    sam::Symbol,
                    path::String)::String

    log_str = "Method, SR, Cost, ISR\n"

    p = plot(x=[], y=[], Geom.point,
                Guide.xlabel("Time"), Guide.ylabel("1 / Success Rate"),
                Theme(background_color=colorant"white",
                    panel_fill=colorant"white",
                    major_label_font=PLOT_FONT,
                    minor_label_font=PLOT_FONT),
                    Guide.manual_color_key("Method",
                    [l[end] for l in method_label],
                    [c[end] for c in method_color]))

    ylim_arr = []
    for m in methods
        log_str *= "$(m), "
        sr = mean([step_success_rate(e) for e in plot_data[m]["error"][1]])
        log_str *= "$(round(sr, digits=4)), "
        isr = [1.0./sr]
        timed = plot_data[m]["time"][1]

        if !any(isnan.(vcat(timed...)))
            qtime = box_scatter_plot(timed)
            qisr = box_scatter_plot(isr)
            log_str *= "$(round(qtime[2], digits=4)), "
            log_str *= "$(round(qisr[2], digits=4))\n"
            if qisr[2] <= 5.0
                append!(p.layers,
                    layer(x=[qtime[2]], y=[qisr[2]],
                            xmin=[qtime[1]],
                            xmax=[qtime[3]],
                            ymin=[qisr[1]],
                            ymax=[qisr[3]],
                            color=[method_color[m]],
                            Geom.point, Geom.errorbar))
            end
        end
    end

    # display(p)
    try
        p |> PDF(joinpath(path,replace("sr_all_medians_$(sam)_$(var)","."=>"")*".pdf"))
    catch  e
        @show e
    end

    log_str *= "\n"
    return log_str
end

"""
Normalized Mean Square Error
"""
function nmse(data::Vector{T}, data_est::Vector{T})::T where T
    normalizer = abs2(maximum(data) - minimum(data))
    res = mean(abs2.(data_est-data))/normalizer
    return res
end

function plot_compare(data::Vector, data_est::Vector)
    alphabet='A':'Z'
    label_n=reshape(["$i (Nominal)" for i in alphabet[1:length(data[1])]],(1,length(data[1])))
    label_e=reshape(["$i (Estimated)" for i in alphabet[1:length(data[1])]],(1,length(data[1])))
    err = sqrt(nmse(reduce(vcat,data), reduce(vcat,data_est)))
    p = plot(title = "RMSE = $(err)")
    plot!(p, reduce(hcat,data)', label=label_n, markershape=:circle, linestyle=:solid)
    plot!(p, reduce(hcat,data_est)', label=label_e, markershape=:cross, linestyle=:dash)
    display(p)
    return p
end

"""
Overall Efficiency (OE) Plots
xaxis: Method
yaxis: OE Score
"""
function oe_plots(plot_data::Dict,
                    noise_level::AbstractArray{Symbol},
                    methods::Vector{Symbol},
                    method_label::Dict,
                    method_color::Dict,
                    sam::Symbol,
                    path::String)::String

    log_str = "Overall Efficiency\n"
    t_succ = zeros(length(methods))
    for (i,m) in enumerate(methods)
        error = plot_data[m]["error"]
        sr = mean.([step_success_rate.(e) for e in error])
        time = plot_data[m]["time"]
        mean_time = mean.(time)
        t_succ_arr = mean_time./sr
        t_succ[i] = mean(t_succ_arr)
        log_str *= "$(m), "
    end

    oe = minimum(t_succ)./t_succ
    log_str *= "\n$(round.(oe, digits=4))\n"

    df = DataFrame()
    for (i,m) in enumerate(methods)
        aux_df = DataFrame(method=method_label[m],
                    data=oe[i], color=method_color[m])
        try
            append!(df, aux_df)
        catch e
            @show e
        end
    end

    p = plot(df, x="method", y="data", color="color", Geom.bar,
            Guide.xlabel("Method"), Guide.ylabel("Overall Efficiency"),
            Theme(background_color=colorant"white",
                panel_fill=colorant"white",
                major_label_font=PLOT_FONT,
                minor_label_font=PLOT_FONT))

    methods_str = reduce((acc,x) -> acc*"_"*string(x), methods; init="")
    try
        if length(noise_level) > 1
            p |> PDF(joinpath(path,"oe_$(sam)$(methods_str)_all.pdf"))
        else
            var = noise_level[1]
            p |> PDF(joinpath(path,replace("oe_$(sam)$(methods_str)_$(var)","."=>"")*".pdf"))
        end
    catch e
        @show e
    end

    return log_str
end

function heatmap(x::AbstractArray,
                y::AbstractArray,
                z::AbstractArray;
                xlabel::String="",
                ylabel::String="",
                zlabel::String="",
                title::String="",
                fillcolor::ColorScheme=ColorSchemes.vik)::Gadfly.Plot
    is = repeat(x, length(x))
    js = vec(repeat(y', length(y)))
    values = z
    p = Gadfly.plot(x=is, y=js, color=values,
        Gadfly.Guide.xlabel(xlabel),
        Gadfly.Guide.ylabel(ylabel),
        Gadfly.Guide.colorkey(zlabel),
        Gadfly.Guide.title(title),
        Gadfly.Theme(minor_label_font=PLOT_FONT,
                    major_label_font=PLOT_FONT,
                    key_title_font=PLOT_FONT,
                    key_label_font=PLOT_FONT),
        # Gadfly.Guide.manual_color_key("Legend",
        #                                 ["True k1",
        #                                 "True k2"],
        #                                 ["magenta",
        #                                 "cyan"]),
        Gadfly.Coord.cartesian(yflip=false,
                        fixed=false,
                        xmin=minimum(x),
                        xmax=maximum(x),
                        ymin=minimum(y),
                        ymax=maximum(y)),
        Gadfly.Scale.ContinuousColorScale(p ->
                ColorSchemes.get(fillcolor, p)),
        Gadfly.Geom.rectbin,
        Gadfly.Scale.x_continuous,
        Gadfly.Scale.y_continuous)
    return p
end

function problem_plot(data::Matrix,
                        t::Vector,
                        type::String)::Gadfly.Plot

    m = size(data)[1]
    alphabet = 'A':'Z'
    label = ["$i" for i in alphabet[1:m]]

    df = DataFrame()
    for i in 1:m
        aux_df = DataFrame(x=data[i,:], State=label[i], t=t)
        append!(df, aux_df)
    end

    if type == "scatter"
        p = plot(df, x=:t, y=:x, color=:State,
                    Guide.xlabel("Time"), Guide.ylabel("State"),
                    Geom.point,
                    Gadfly.Coord.cartesian(yflip=false,
                                            fixed=false,
                                            xmin=minimum(df[:t]),
                                            xmax=maximum(df[:t]),
                                            ymin=minimum(df[:x]),
                                            ymax=maximum(df[:x])),
                    Gadfly.Theme(minor_label_font=PLOT_FONT,
                                    major_label_font=PLOT_FONT,
                                    key_title_font=PLOT_FONT,
                                    key_label_font=PLOT_FONT)
                    )

    elseif type == "line"
        p = plot(df, x=:t, y=:x, color=:State,
                    Guide.xlabel("Time"), Guide.ylabel("State"),
                    Geom.line,
                    Gadfly.Coord.cartesian(yflip=false,
                                            fixed=false,
                                            xmin=minimum(df[:t]),
                                            xmax=maximum(df[:t]),
                                            ymin=minimum(df[:x]),
                                            ymax=maximum(df[:x])),
                    Gadfly.Theme(minor_label_font=PLOT_FONT,
                                    major_label_font=PLOT_FONT,
                                    key_title_font=PLOT_FONT,
                                    key_label_font=PLOT_FONT)
                    )

    elseif type == "scatter_line"
        p = plot(df, x=:t, y=:x, color=:State,
                    Guide.xlabel("Time"), Guide.ylabel("State"),
                    Geom.point, Geom.line,
                    Gadfly.Coord.cartesian(yflip=false,
                                            fixed=false,
                                            xmin=minimum(df[:t]),
                                            xmax=maximum(df[:t]),
                                            ymin=minimum(df[:x]),
                                            ymax=maximum(df[:x])),
                    Gadfly.Theme(minor_label_font=PLOT_FONT,
                                    major_label_font=PLOT_FONT,
                                    key_title_font=PLOT_FONT,
                                    key_label_font=PLOT_FONT)
                    )
    else
        p = plot()
    end

    return p
end

end

max_stars_repo_name: atgmello/DEParamEst

max_issues_repo_name: atgmello/DEParamEst

max_forks_repo_name: atgmello/DEParamEst


text:
using DelimitedFiles
using Dates
using StatsBase

df = DateFormat("yyyy-mm-dd HH:MM")

function solve(input)
    raw = readdlm(input,'\n')
    shifts = []
    guardIDs = Array{Int16,1}()
    mapreduce(parseShift, fillShifts, raw, init=(shifts,guardIDs))
    sort!(shifts, by = x -> x[1])
    unique!(guardIDs)
    numGuards = length(guardIDs)
    guardIndices = Dict(zip(guardIDs,1:numGuards))
    reverseGuardIndices = Dict(zip(1:numGuards,guardIDs))
    #part 1
    sleepGrid = zeros(Int16,12,31,numGuards,60)
    lastGuardIndex = nothing
    lastGuardTimeStamp = nothing
    reduce(fillGrid, shifts, init=(lastGuardIndex,lastGuardTimeStamp,sleepGrid,guardIndices))
    @show sum(sleepGrid[:,:,:,:],dims=[1,2,4])
    sleepyGuardIndex = argmax(sum(sleepGrid[:,:,:,:],dims=[1,2,4]))[3]
    sleepyGuard = reverseGuardIndices[sleepyGuardIndex]
    println(sleepyGuard)
    sleepyMinute = argmax(sum(sleepGrid[:,:,sleepyGuardIndex,:],dims=[1,2]))[3] - 1
    println(sleepyMinute)
    part1=sleepyGuard*sleepyMinute
    println("result part 1: $part1")
    #part 2
    marmotte=argmax(sum(sleepGrid[:,:,:,:],dims=[1,2]))
    dormeur=reverseGuardIndices[marmotte[3]]
    minutemax=marmotte[4]-1
    part2=dormeur*minutemax
    println("result part 2: $part2")
end

function parseShift(s)
    m = match(r"\[(?<timestamp>.+)\] (Guard #(?<guardID>.+) begins shift|(?<sleeping>falls asleep$)|(?<awake>wakes up$))", s)
    t = DateTime(m[:timestamp],df)
    guardID = m[:guardID]
    if guardID!= nothing
        guardID = tryparse(Int16,guardID)
    end
    sleeping = (m[:sleeping]!=nothing)
    awake = (m[:awake]!=nothing)
    return (t,guardID,sleeping,awake)
end

function fillShifts((shifts,guardIDs),(t,guardID,sleeping,awake))
    if guardID!=nothing
        push!(guardIDs,guardID)
    end
    push!(shifts,(t,guardID,sleeping,awake))
    return (shifts,guardIDs)
end

function fillGrid((lastGuardIndex,lastGuardTimeStamp,grid,guardIndices),(timestamp,guardID,sleeping,awake))
    if guardID!=nothing
        lastGuardIndex=guardIndices[guardID]
    elseif sleeping
        grid[Dates.month(timestamp),Dates.day(timestamp),lastGuardIndex,Dates.minute(timestamp)+1]=1
    elseif awake
        sleepMinute = lastGuardTimeStamp
        while timestamp > sleepMinute
            sleepMinute = sleepMinute + Dates.Minute(1)
            grid[Dates.month(sleepMinute),Dates.day(sleepMinute),lastGuardIndex,Dates.minute(sleepMinute)+1]=1
        end
        grid[Dates.month(timestamp),Dates.day(timestamp),lastGuardIndex,Dates.minute(timestamp)+1]=0
    end
    lastGuardTimeStamp=timestamp
    return (lastGuardIndex,lastGuardTimeStamp,grid,guardIndices)

end

solve("input.txt")


text:
using MultivariateStats
using Base.Test
import StatsBase

srand(15678)

## icagfun

f = icagfun(:tanh)
u, v = evaluate(f, 1.5)
@test u ≈ 0.905148253644866438242
@test v ≈ 0.180706638923648530597

f = icagfun(:tanh, Float32)
u, v = evaluate(f, 1.5f0)
@test u ≈ 0.90514827f0
@test v ≈ 0.18070662f0

f = icagfun(:tanh, 1.5)
u, v = evaluate(f, 1.2)
@test u ≈ 0.946806012846268289646
@test v ≈ 0.155337561057228069719

f = icagfun(:tanh, 1.5f0)
u, v = evaluate(f, 1.2f0)
@test u ≈ 0.94680610f0
@test v ≈ 0.15533754f0

f = icagfun(:gaus)
u, v = evaluate(f, 1.5)
@test u ≈ 0.486978701037524594696
@test v ≈ -0.405815584197937162246

f = icagfun(:gaus, Float32)
u, v = evaluate(f, 1.5f0)
@test u ≈ 0.4869787f0
@test v ≈ -0.40581557f0

## data

# sources
n = 1000
k = 3
m = 8

t = linspace(0.0, 10.0, n)
s1 = sin.(t * 2)
s2 = s2 = 1.0 - 2.0 * Bool[isodd(floor(Int, x / 3)) for x in t]
s3 = Float64[mod(x, 5.0) for x in t]

s1 += 0.1 * randn(n)
s2 += 0.1 * randn(n)
s3 += 0.1 * randn(n)

S = hcat(s1, s2, s3)'
@assert size(S) == (k, n)

A = randn(m, k)

X = A * S
mv = vec(mean(X,2))
@assert size(X) == (m, n)
C = cov(X, 2)

# FastICA

M = fit(ICA, X, k; do_whiten=false, tol=Inf)
@test isa(M, ICA)
@test indim(M) == m
@test outdim(M) == k
@test mean(M) == mv
W = M.W
@test transform(M, X) ≈ W' * (X .- mv)
@test W'W ≈ eye(k)

M = fit(ICA, X, k; do_whiten=true, tol=Inf)
@test isa(M, ICA)
@test indim(M) == m
@test outdim(M) == k
@test mean(M) == mv
W = M.W
@test W'C * W ≈ eye(k)

@test_throws StatsBase.ConvergenceException fit(ICA, X, k; do_whiten=true, tol=1)

# Use data of different type

XX = convert(Matrix{Float32}, X)

M = fit(ICA, XX, k; do_whiten=true, tol=Inf)
@test eltype(mean(M)) == Float32
@test eltype(M.W) == Float32

M = fit(ICA, XX, k; do_whiten=false, tol=Inf)
@test isa(M, ICA)
@test eltype(mean(M)) == Float32
@test eltype(M.W) == Float32
W = M.W
@test transform(M, X) ≈ W' * convert(Matrix{Float32}, (X .- mv))
@test W'W ≈ eye(k)


text:
"""
    AbstractMatrix(k::Kernel)

Convert `k` into an `AbstractMatrix`, if such a representation exists.
"""
function AbstractMatrix(k::Kernel)
    @assert isfinite(size(k, 1))
    return pairwise(k, eachindex(k, 1))
end
AbstractMatrix(k::EmpiricalKernel) = k.Σ

"""
    AbstractMatrix(k::CrossKernel)

Convert `k` into an `AbstractMatrix`, if such a representation exists.
"""
function AbstractMatrix(k::CrossKernel)
    @assert isfinite(size(k, 1))
    @assert isfinite(size(k, 2))
    return pairwise(k, eachindex(k, 1), eachindex(k, 2))
end


text:
"""
    Patch{T <: Any}

Data structure which contains the `diff` between two `Vector`. The second
`Vector` can then be rebuilt by patching the first one with the `Patch`.
"""
struct Patch{T <: Any}
    range::Vector{UnitRange{Int32}}
    data::Vector{T}
end

"""
    diff(old::Vector{T}, new::Vector{T}, sizeof_T::Int)

Generate a `Patch` by delta compressing the difference between `old` and `new`.
"""
function diff(old::Vector{T}, new::Vector{T}, sizeof_T::Int)::Patch{T} where T <: Any
    min_range_dist = (2 * sizeof(Int32)) / sizeof_T - 1
    range = Vector{UnitRange{Int32}}()
    data = Vector{T}()

    start = 1
    for i in 1:length(old)
        # current elements do not match
        if old[i] != new[i]
            # is the last serie of identical elements big enough?
            if i - start >= min_range_dist
                push!(range, start:i-1)
                append!(data, new[i])
            else # not enough elements...
                append!(data, new[start:i])
            end
            start = i + 1
        end
    end

    # we may have reached the end with identical elements...
    if start <= length(old)
        if length(old) - start >= min_range_dist
            push!(range, start:length(old))
        else
            append!(data, new[start:length(old)])
        end
    end

    return Patch{T}(range, data)
end

"""
    patch(old::Vector{T}, patch::Patch{T})

Apply `patch` to `old`.
"""
function patch(old::Vector{T}, patch::Patch{T})::Vector{T} where T <: Any
    patched = similar(old)

    curr = curr_data = 1
    for range in patch.range

        if range.start != curr
            num_val_to_copy = range.start - curr
            patched[curr:curr+num_val_to_copy-1] = patch.data[curr_data:curr_data+num_val_to_copy-1]
            curr += num_val_to_copy
            curr_data += num_val_to_copy
        end

        patched[range] = old[range]
        curr += length(range)
    end

    patched[curr:end] = patch.data[curr_data:end]

    return patched
end


text:
using GCMAES
using Zygote
using Test

try @mpirun using Elemental catch end

rastrigin(x) = 10length(x) + sum(x.^2 .- 10 .* cos.(2π .* x))

D = 5000
x0 = fill(0.3f0, D)
σ0 = 0.2f0
lo = fill(-5.12f0, D)
hi = fill(5.12f0, D)

@mpirun GCMAES.minimize(rastrigin, x0, σ0, lo, hi, λ = 200, maxiter = 200)

@mpirun GCMAES.minimize(rastrigin, x0, σ0, lo, hi, λ = 200, maxiter = 200, autodiff = true)

@mpirun GCMAES.minimize(x0, σ0, lo, hi, λ = 3, maxiter = 200) do x
    GCMAES.allmean(rastrigin(x), GCMAES.localcomm())
end

rm("CMAES.bson", force = true)


text:
# aerodyne_parsetime.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 13.12.2016
# Last Edit: 11.04.2019

"""# aerodyne_parsetime

Parse the STR and STC file name timestamps into DateTime and sort

`time,Files = aerodyne_parsetime(F)`\n
* **time**::DateTime = Sorted timestamps
* **Files**::Array{String} = Temporally sorted file list
* **F**::Array{String} = List of files to parse and sort
"""
function aerodyne_parsetime(F::Array{T,1}) where T <: String
	## Parse the names of files into a readable time ##
	t = Array{DateTime}(undef,Int(length(F))) # Preallocate time column
	df = Dates.DateFormat("yyyymmdd_HHMMSS") # Date format
	for i = 1:1:length(F)
		t[i] = DateTime("20" * F[i][end-16:end-4],df)
	end
	I = sortperm(t) # Sort index (by time)
	t = t[I] # Sort time by time
	F = F[I] # Sort the files by time

	return t,F
end # End of aerodyne_parsetime(F::Array{String,1})

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
@json struct Comment
    url::String
    pull_request_review_id::Int
    id::Int
    node_id::String
    diff_hunk::String
    path::String
    position::Int
    original_position::Int
    commit_id::String
    original_commit_id::String
    in_reply_to_id::Int
    user::User
    body::String
    created_at::DateTime
    updated_at::DateTime
    html_url::String
    pull_request_url::String
    author_association::String
    _links => links::Links
    start_line::Int
    original_start_line::Int
    start_side::String
    line::Int
    original_line::Int
    side::String
end

## Pull Request (Issue) Comments

# List All Pull Request (Issue) Comments
# https://docs.github.com/en/rest/reference/issues#list-issue-comments
function endpoint(
    ::GitHubAPI,
    ::typeof(list_pull_request_comments),
    owner::AStr,
    repo::AStr,
    pull_request_id::Integer,
)
    return Endpoint(:GET, "/repos/$owner/$repo/issues/$pull_request_id/comments")
end
into(::GitHubAPI, ::typeof(list_pull_request_comments)) = Vector{Comment}

# Get Single Pull Request (Issue) Comment
# https://docs.github.com/en/rest/reference/issues#get-an-issue-comment
function endpoint(
    ::GitHubAPI,
    ::typeof(get_pull_request_comment),
    owner::AStr,
    repo::AStr,
    comment_id::Integer,
)
    return Endpoint(
        :GET, "/repos/$owner/$repo/issues/comments/$comment_id")
end
info(::GitHubAPI, ::typeof(get_pull_request_comment)) = Comment

# Create Pull Request (Issue) Comment
# https://docs.github.com/en/rest/reference/issues#create-an-issue-comment
function endpoint(
    ::GitHubAPI,
    ::typeof(create_pull_request_comment),
    owner::AStr,
    repo::AStr,
    pull_request_id::Integer,
)
    return Endpoint(:POST, "/repos/$owner/$repo/issues/$pull_request_id/comments")
end
into(::GitHubAPI, ::typeof(create_pull_request_comment)) = Comment

# Update Pull Request (Issue) Comment
# https://docs.github.com/en/rest/reference/issues#update-an-issue-comment
function endpoint(
    ::GitHubAPI,
    ::typeof(update_pull_request_comment),
    owner::AStr,
    repo::AStr,
    comment_id::Integer,
)
    return Endpoint(:PATCH, "/repos/$owner/$repo/issues/comments/$comment_id")
end
into(::GitHubAPI, ::typeof(update_pull_request_comment)) = Comment

# Delete Pull Request (Issue) Comment
# https://docs.github.com/en/rest/reference/issues#delete-an-issue-comment
function endpoint(
    ::GitHubAPI,
    ::typeof(delete_pull_request_comment),
    owner::AStr,
    repo::AStr,
    comment_id::Integer,
)
    return Endpoint(:DELETE, "/repos/$owner/$repo/issues/comments/$comment_id")
end
into(::GitHubAPI, ::typeof(delete_pull_request_comment)) = Comment


text:
# An observer pattern example
# See a similar example at https://easymock.org/getting-started.html

using Pretend
Pretend.activate()

using BinaryTraits
using BinaryTraits.Prefix: Can

using Test

# Formal interface
@trait HandleNewDocument
@implement Can{HandleNewDocument} by document_added(_, document)

# A ProofReader can handle any new documents
struct ProofReader end
document_added(reader::ProofReader, document) = println("notified document_added => $document")

# We will claim that it satisfies the HandleNewDocument trait
@assign ProofReader with Can{HandleNewDocument}
@check ProofReader

# A DocumentManager deals with all documentation needs
struct DocumentManager
    listeners::Vector
end

# Add a new document to the system
# 1. persist the document
# 2. notify any collaborator that might be interested about this event
function add(dm::DocumentManager, title, content)
    println("added doc '$title' with content '$content'")
    foreach(x -> document_added(x, title), dm.listeners)
end

reader = ProofReader()
dm = DocumentManager([reader])

# just a check
add(dm, "Life is boring", "Really? I don't think so.")

# Annotate the mockable function
@mockable document_added(reader::ProofReader, document) = println("notified document_added => $document")

# Create a mock(patch)
apply(document_added => (reader::ProofReader, document) -> println("mocked! ha!")) do
    add(dm, "Life is boring", "Really? I don't think so.")
    @test called_exactly_once(document_added, reader, "Life is boring")
end


text:
using CSV,DataFrames,Statistics,StatsBase

include( "matching/matching_funcs.jl")

frs_all_years_scot_he = create_frs( 2015:2019 )
shs_all_years = create_shs( 2016:2019 )

donor,recip = create_donor_and_recip_datasets( frs_all_years_scot_he, shs_all_years )

#
# save everything
#
#CSV.write( "data/merging/shs_donor_data.tab", donor )
#CSV.write( "data/merging/frs_recip_data.tab", recip )
CSV.write( "data/merging/shs_all_years.tab", shs_all_years )
CSV.write( "data/merging/frs_all_years_scot_he.tab", frs_all_years_scot_he )

make_initial_match!( recip, donor )

CSV.write( "data/merging/shs_donor_data.tab", donor; quotestrings=true, delim='\t' )
CSV.write( "data/merging/frs_shs_merging_indexes.tab", recip; quotestrings=true, delim='\t' )

# donor = CSV.File( "data/merging/shs_donor_data.tab"; types=Dict(:uniqidnew => String))|>DataFrame
# recip = CSV.File( "data/merging/frs_recip_data.tab" ) |> DataFrame

unmatched = recip[(recip.shs_datayear_1 .== 0),:]
unmatched[!,critmatche1]

# match shelter on random shelter
unmatched = Vector(((recip.shs_datayear_1 .== 0).&(recip.shelter_1 .== 1)))
n = match_up_unmatched!( recip, donor, [:shelter], unmatched )

CSV.write( "data/merging/shs_donor_data.tab", donor; quotestrings=true, delim='\t' )
CSV.write( "data/merging/frs_shs_merging_indexes.tab", recip; quotestrings=true, delim='\t' )

"""

final missing 7 cases after matching any sheltered home:

 Row │ sernum  datayear  shelter_1  singlepar_1  numadults_2  numkids_2  empstathigh_1  agehigh_1
     │ Int64   Int64     Int64      Int64        Int64        Int64      Int64          Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────
   1 │   6171        16          0            1            1          2              1         62
   2 │  16510        16          0            0            3          1              1         19
   3 │  10973        17          0            0            3          4              1         42
   4 │  18803        17          0            0            3          2              5         80
   5 │   4861        18          0            0            3          2              1         19
   6 │   5234        18          0            0            3          0              6         18
   7 │  14536        18          0            0            2          2              5         80

"""

# reload so we can just start the script here
# donor = CSV.File( "data/merging/shs_donor_data.tab"; delim='\t' ) |> DataFrame
#
# no idea whatsoever why the 1st uniq needs cast, but it does seem to ..
#
# recip = CSV.File( "data/merging/frs_shs_merging_indexes.tab"; delim='\t', types=Dict("shs_uniqidnew_1"=>String) ) |> DataFrame

# 1) so drop age & emp status for last 7

final_unmatched = recip[(recip.shs_datayear_1 .== 0),:]

final_targets = [:singlepar,:numadults,:numkids]

# get rid of this cast
final_unmatched = Vector((recip.shs_datayear_1 .== 0))

#
# matching the remaining 7 on sp, ads, kids
#
n = match_up_unmatched!( recip, donor, final_targets, final_unmatched )
#
# & we're down to 1 unmatched ..
#
"""
 Row │ sernum  datayear  shelter_1  singlepar_1  numadults_2  numkids_2  empstathigh_1  agehigh_1
     │ Int64   Int64     Int64      Int64        Int64        Int64      Int64          Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────
   1 │  10973        17          0            0            3          4              1         42

"""

# try with just kids, no adults
final_targets_2 = [:tenure,:singlepar,:numkids,:datayear]
final_unmatched_2 = Vector(recip.shs_datayear_1 .== 0)
n = match_up_unmatched!( recip, donor, final_targets_2, final_unmatched_2 )

# no more than 4 children in SHS
# maximum( donor.numkids_1 ) => 4
# maximum( donor.

CSV.write( "data/merging/shs_donor_data.tab", donor; quotestrings=true, delim='\t' )
CSV.write( "data/merging/frs_shs_merging_indexes.tab", recip; quotestrings=true, delim='\t' )

mhh = CSV.File( "data/model_households_scotland.tab"; delim='\t') |> DataFrame
shs_councils = CSV.File( "data/merging/la_mappings.csv"; delim=',') |> DataFrame
target_pops = CSV.File( "data/merging/hhlds_and_people_2019_nrs_estimates.csv" ) |> DataFrame
shs_hhn = count_councils(shs_all_years, shs_councils )
#
# actual sampling frequencies by council, pooled over all shs years
# we use these as probability weights
#
inv_freqs = Dict()

for p in eachrow( target_pops )
    inv_freqs[p.code] = p.hhlds_2019/shs_hhn[p.code]
end
#
# idiot check that we sum back up
#
s = 0.0
for p in eachrow( target_pops )
    global s
    s += inv_freqs[p.code]*shs_hhn[p.code]
end

add_in_las_to_recip!( recip, shs_all_years, shs_councils )

@assert s  ≈ 2495622

tot_hhlds_19 = 2_495_622

n = add_council_to_frs!(
    # ;
    mhh   = mhh,
    recip = recip,
    shs_councils = shs_councils,
    inv_freqs = inv_freqs )

aw = tot_hhlds_19/sum(values(n))

println( "|code|name|shs sample freq| n |" )
println( "|----|----|-----------|------|" )
for p in eachrow( shs_councils )
    v = inv_freqs[p.Code]
    c = n[p.Code]*aw
    println( "| $(p.Code) | $(p.name) | $v | $c")
end

for r in eachrow(mhh)
    n[r.council] += 1
end

CSV.write( "data/model_households_scotland.tab", mhh; delim='\t')

#
# todo : add bedrooms bedroom6 frs capped at 6 hc4 shs
#

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
mutable struct BidLimitOrder <: LimitOrder
    dealer::Union{AbstractInvestor, String}
    price::Int64
    quantity::Int64

    from::AbstractAsset # This is the asset that is traded away (ie. currency in buy order)
    to::AbstractAsset # This is the asset that is traded to (ie. stock in buy order)

    timestamp::Int64
    function BidLimitOrder(dealer::Union{AbstractInvestor, String}; price, quantity, from, to)
        new(
            dealer,
            price,
            quantity,
            from,
            to,
            get_time()
        )
    end
    function BidLimitOrder(dealer::Union{AbstractInvestor, String}; price, quantity, market)
        new(
            dealer,
            price,
            quantity,
            market.currency,
            market.traded_asset,
            get_time()
        )
    end
end

mutable struct AskLimitOrder <: LimitOrder
    dealer::Union{AbstractInvestor, String}
    price::Int64
    quantity::Int64

    from::AbstractAsset # This is the asset that is traded away (ie. stock in sell order)
    to::AbstractAsset # This is the asset that is traded to (ie. currency in sell order)

    timestamp::Int64
    function AskLimitOrder(dealer::Union{AbstractInvestor, String}; price::Int64, quantity::Int64, from::AbstractAsset, to::AbstractAsset)
        new(
            dealer,
            price,
            quantity,
            from,
            to,
            get_time()
        )
    end
    function AskLimitOrder(dealer::Union{AbstractInvestor, String}; price::Int64, quantity::Int64, market::AbstractMarket)
        new(
            dealer,
            price, quantity,
            market.traded_asset,
            market.currency,
            get_time()
        )
    end
end

"Get amount of asset commited in the order (commits currency)"
function get_amount_commited(order::BidLimitOrder)
    return order.price * order.quantity
end

"Get amount of asset commited in the order (commits the asset)"
function get_amount_commited(order::AskLimitOrder)
    return order.quantity
end

"Push order to the orderbook of the market"
function Base.push!(market::AbstractMarket, order::BidLimitOrder)
    push!(market.bid_limit_orders, order)
end

"Push order to the orderbook of the market"
function Base.push!(market::AbstractMarket, order::AskLimitOrder)
    push!(market.ask_limit_orders, order)
end

get_clearable_quantity(buy::BidLimitOrder, sell::AskLimitOrder) = min(buy.quantity, sell.quantity)

CURRENT_TIME = 0
# Let defines new name space
let state = 0
    global get_time() = (state += 1)
 end


text:
const MOI  = Gurobi.MOI
const MOIT = MOI.Test
const MOIB = MOI.Bridges

const GUROBI_ENV = Gurobi.Env()

@testset "Unit Tests" begin
    config = MOIT.TestConfig()
    solver = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
    MOIT.basic_constraint_tests(solver, config)
    MOIT.unittest(solver, config,
        ["solve_affine_interval", "solve_qcp_edge_cases"])
    @testset "solve_affine_interval" begin
        MOIT.solve_affine_interval(
            MOIB.SplitInterval{Float64}(Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)),
            config
        )
    end
    @testset "solve_qcp_edge_cases" begin
        MOIT.solve_qcp_edge_cases(solver,
            MOIT.TestConfig(atol=1e-3)
        )
    end
    MOIT.modificationtest(solver, config, [
        "solve_func_scalaraffine_lessthan"
    ])
end

@testset "Linear tests" begin
    @testset "Default Solver"  begin
        solver = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        MOIT.contlineartest(solver, MOIT.TestConfig(), [
            # This requires interval constraint.
            "linear10", "linear10b",
            # This requires an infeasiblity certificate for a variable bound.
            "linear12"
        ])
    end
    @testset "linear10" begin
        MOIT.linear10test(
            MOIB.SplitInterval{Float64}(Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)),
            MOIT.TestConfig()
        )
        MOIT.linear10btest(
            MOIB.SplitInterval{Float64}(Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)),
            MOIT.TestConfig()
        )
    end
    @testset "No certificate" begin
        MOIT.linear12test(
            Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0, InfUnbdInfo=0),
            MOIT.TestConfig(infeas_certificates=false)
        )
    end
end

@testset "Quadratic tests" begin
    MOIT.contquadratictest(
        Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0),
        MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals=false, query=false)
    )
end

@testset "Linear Conic tests" begin
    MOIT.lintest(
        Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0),
        MOIT.TestConfig()
    )
end

@testset "Integer Linear tests" begin
    MOIT.intlineartest(
        Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0),
        MOIT.TestConfig(),
        ["int3"]  # int3 has interval constriants
    )
    @testset "int3" begin
        MOIT.int3test(
            MOIB.SplitInterval{Float64}(Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)),
            MOIT.TestConfig()
        )
    end
end
@testset "ModelLike tests" begin
    solver = Gurobi.Optimizer(GUROBI_ENV)
    @test MOI.get(solver, MOI.SolverName()) == "Gurobi"
    @testset "default_objective_test" begin
         MOIT.default_objective_test(solver)
     end
     @testset "default_status_test" begin
         MOIT.default_status_test(solver)
     end
    @testset "nametest" begin
        MOIT.nametest(solver)
    end
    @testset "validtest" begin
        MOIT.validtest(solver)
    end
    @testset "emptytest" begin
        MOIT.emptytest(solver)
    end
    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(solver)
    end
    @testset "copytest" begin
        MOIT.copytest(solver, Gurobi.Optimizer(GUROBI_ENV))
    end
end

@testset "Gurobi Callback" begin
    @testset "Generic callback" begin
        m = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0)
        x = MOI.add_variable(m)
        MOI.add_constraint(m, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
        MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction{Float64}(
                [MOI.ScalarAffineTerm{Float64}(1.0, x)],
                0.0
            )
        )

        cb_calls = Int32[]
        function callback_function(cb_data::Gurobi.CallbackData, cb_where::Int32)
            push!(cb_calls, cb_where)
            nothing
        end

        MOI.set(m, Gurobi.CallbackFunction(), callback_function)
        MOI.optimize!(m)

        @test length(cb_calls) > 0
        @test Gurobi.CB_MESSAGE in cb_calls
        @test Gurobi.CB_PRESOLVE in cb_calls
        @test !(Gurobi.CB_MIPSOL in cb_calls)
    end

    @testset "Lazy cut" begin
        m = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0, Cuts=0, Presolve=0, Heuristics=0, LazyConstraints=1)
        MOI.Utilities.loadfromstring!(m,"""
            variables: x, y
            maxobjective: y
            c1: x in Integer()
            c2: y in Integer()
            c3: x in Interval(0.0, 2.0)
            c4: y in Interval(0.0, 2.0)
        """)
        x = MOI.get(m, MOI.VariableIndex, "x")
        y = MOI.get(m, MOI.VariableIndex, "y")

        # We now define our callback function that takes two arguments:
        #   (1) the callback handle; and
        #   (2) the location from where the callback was called.
        # Note that we can access m, x, and y because this function is defined
        # inside the same scope
        cb_calls = Int32[]
        function callback_function(cb_data::Gurobi.CallbackData, cb_where::Int32)
            push!(cb_calls, cb_where)
            if cb_where == Gurobi.CB_MIPSOL
                Gurobi.loadcbsolution!(m, cb_data, cb_where)
                x_val = MOI.get(m, MOI.VariablePrimal(), x)
                y_val = MOI.get(m, MOI.VariablePrimal(), y)
                # We have two constraints, one cutting off the top
                # left corner and one cutting off the top right corner, e.g.
                # (0,2) +---+---+ (2,2)
                #       |xx/ \xx|
                #       |x/   \x|
                #       |/     \|
                # (0,1) +       + (2,1)
                #       |       |
                # (0,0) +---+---+ (2,0)
                TOL = 1e-6  # Allow for some impreciseness in the solution
                if y_val - x_val > 1 + TOL
                    Gurobi.cblazy!(cb_data, m,
                        MOI.ScalarAffineFunction{Float64}(
                            MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                            0.0
                        ),
                        MOI.LessThan{Float64}(1.0)
                    )
                elseif y_val + x_val > 3 + TOL
                    Gurobi.cblazy!(cb_data, m,
                        MOI.ScalarAffineFunction{Float64}(
                            MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                            0.0
                        ),
                        MOI.LessThan{Float64}(3.0)
                    )
                end
            end
        end

        MOI.set(m, Gurobi.CallbackFunction(), callback_function)
        MOI.optimize!(m)

        @test MOI.get(m, MOI.VariablePrimal(), x) == 1
        @test MOI.get(m, MOI.VariablePrimal(), y) == 2

        @test length(cb_calls) > 0
        @test Gurobi.CB_MESSAGE in cb_calls
        @test Gurobi.CB_PRESOLVE in cb_calls
        @test Gurobi.CB_MIPSOL in cb_calls
    end
end

@testset "LQOI Issue #38" begin
    # https://github.com/JuliaOpt/LinQuadOptInterface.jl/issues/38#issuecomment-407625187
    _getinner(opt::Gurobi.Optimizer) = opt.inner
    @inferred _getinner(Gurobi.Optimizer(GUROBI_ENV))
end

@testset "User limit handling (issue #140)" begin
    # Verify that we return the correct status codes when a mixed-integer
    # problem has been solved to a *feasible* but not necessarily optimal
    # solution. To do that, we will set up an intentionally dumbed-down
    # Gurobi Gurobi.Optimizer (with all heuristics and pre-solve turned off) and
    # ask it to solve a classic knapsack problem. Setting SolutionLimit=1
    # forces the solver to return after its first feasible MIP solution,
    # which tests the right part of the code without relying on potentially
    # flaky or system-dependent time limits.
    m = Gurobi.Optimizer(GUROBI_ENV, OutputFlag=0, SolutionLimit=1,
                         Heuristics=0.0, Presolve=0)
    N = 100
    x = MOI.add_variables(m, N)
    for xi in x
        MOI.add_constraint(m, MOI.SingleVariable(xi), MOI.ZeroOne())
        MOI.set(m, MOI.VariablePrimalStart(), xi, 0.0)
    end
    # Given a collection of items with individual weights and values,
    # maximize the total value carried subject to the constraint that
    # the total weight carried is less than 10.
    if VERSION >= v"0.7-"
        Random.seed!(1)
    else
        srand(1)
    end
    item_weights = rand(N)
    item_values = rand(N)
    MOI.add_constraint(m,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0))
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(.-item_values, x), 0.0))
    MOI.optimize!(m)

    @test MOI.get(m, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    # We should have a primal feasible solution:
    @test MOI.get(m, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    # But we have no dual status:
    @test MOI.get(m, MOI.DualStatus()) == MOI.NO_SOLUTION
end

@testset "Constant objective (issue #111)" begin
    m = Gurobi.Optimizer(GUROBI_ENV)
    x = MOI.add_variable(m)
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 2.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 2.0
    @test Gurobi.get_dblattr(m.inner, "ObjCon") == 2.0

    MOI.modify(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarConstantChange(3.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 3.0
    @test Gurobi.get_dblattr(m.inner, "ObjCon") == 3.0
end

@testset "Env" begin
    @testset "User-provided" begin
        env = Gurobi.Env()
        model_1 = Gurobi.Optimizer(env)
        @test model_1.inner.env === env
        model_2 = Gurobi.Optimizer(env)
        @test model_2.inner.env === env
        # Check that finalizer doesn't touch env when manually provided.
        finalize(model_1.inner)
        @test Gurobi.is_valid(env)
    end
    @testset "Automatic" begin
        model_1 = Gurobi.Optimizer()
        model_2 = Gurobi.Optimizer()
        @test model_1.inner.env !== model_2.inner.env
        # Check that env is finalized with model when not supplied manually.
        finalize(model_1.inner)
        @test !Gurobi.is_valid(model_1.inner.env)
    end
    @testset "Env when emptied" begin
        @testset "User-provided" begin
            env = Gurobi.Env()
            model = Gurobi.Optimizer(env)
            @test model.inner.env === env
            @test Gurobi.is_valid(env)
            MOI.empty!(model)
            @test model.inner.env === env
            @test Gurobi.is_valid(env)
        end
        @testset "Automatic" begin
            model = Gurobi.Optimizer()
            env = model.inner.env
            MOI.empty!(model)
            @test model.inner.env !== env
            @test Gurobi.is_valid(model.inner.env)
        end
    end
end


text:
### Define the target distribution
using Distributions: logpdf, MvNormal

D = 10
ℓπ(θ) = logpdf(MvNormal(zeros(D), ones(D)), θ)

### Build up a HMC sampler to draw samples
using AdvancedHMC, ForwardDiff

# Parameter settings
n_samples, n_adapts = 12_000, 2_000
θ₀ = rand(D)    # draw a random starting points

# Define metric space, Hamiltonian, sampling method and adaptor
metric = DiagEuclideanMetric(D)
hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)  # or, Hamiltonian(metric, ℓπ, ∂ℓπ∂θ) for hand-coded gradient ∂ℓπ∂θ
integrator = Leapfrog(find_good_eps(hamiltonian, θ₀))
proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(Preconditioner(metric), NesterovDualAveraging(0.8, integrator))

# Draw samples via simulating Hamiltonian dynamics
# - `samples` will store the samples
# - `stats` will store statistics for each sample
samples, stats = sample(hamiltonian, proposal, θ₀, n_samples, adaptor, n_adapts; progress=true)


text:
"""
    LRhoKLocalOpProblem <: AbstractProblem

Problem or finding the steady state of a ℒdagℒ matrix by computing
𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ |ρ⟩⟩|² using the sparse Liouvillian matrix.

DO NOT USE WITH COMPLEX-WEIGHT NETWORKS, AS IT DOES NOT WORK
"""
struct LRhoKLocalOpProblem{B, SM1, SM} <: LRhoSquaredProblem where {B<:Basis,
                                                 SM<:SparseMatrixCSC}
    HilbSpace::B            # 0
    HnH::SM1
    L_ops::Vector{SM}       # 4
    ρss
end

LRhoKLocalOpProblem(gl::GraphLindbladian) = LRhoKLocalOpProblem(STD_REAL_PREC, gl)
function LRhoKLocalOpProblem(T, gl::GraphLindbladian)
    HnH, c_ops, c_ops_t = to_linear_operator(gl, Complex{real(T)})
    return LRhoKLocalOpProblem(basis(gl), HnH, c_ops, 0.0)
end

basis(prob::LRhoKLocalOpProblem) = prob.HilbSpace

# Standard method dispatched when the state is generic (non lut).
# will work only if 𝝝 and 𝝝p are the same type (and non lut!)
function compute_Cloc!(LLO_i, ∇lnψ, prob::LRhoKLocalOpProblem,
                       net::MatrixNet, 𝝝::S,
                       lnψ=net(𝝝), 𝝝p::S=deepcopy(𝝝)) where {S}
    # hey
    HnH = prob.HnH
    L_ops = prob.L_ops

    set_index!(𝝝p, index(𝝝))
    𝝝p_row = row(𝝝p)
    𝝝p_col = col(𝝝p)

    for el=LLO_i
      el .= 0.0
    end

    C_loc = zero(Complex{real(out_type(net))})

    # ⟨σ|Hρ|σt⟩ (using hermitianity of HdH)
    # diffs_hnh = row_valdiff(HnH, row(𝝝))
    set_index!(𝝝p_col, index(col(𝝝)))
    for op=operators(HnH)
        r=local_index(row(𝝝), sites(op))
        for (mel, changes)=op.op_conns[r] #diffs_hnh
            set_index!(𝝝p_row, index(row(𝝝)))
            for (site,val)=changes
                setat!(𝝝p_row, site, val)
            end

            lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
            C_loc_i  =  -1.0im * mel * exp(lnψ_i - lnψ)
            for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
              LLOave .+= C_loc_i .* _∇lnψ
            end
            C_loc  += C_loc_i
        end
    end #operators

    # ⟨σ|ρHᴴ|σt⟩
    # diffs_hnh = row_valdiff(HnH, col(𝝝))
    set_index!(𝝝p_row, index(row(𝝝)))
    for op=operators(HnH)
        r=local_index(col(𝝝), sites(op))
        for (mel, changes)=op.op_conns[r]
            set_index!(𝝝p_col, index(col(𝝝)))
            for (site,val)=changes
                setat!(𝝝p_col, site, val)
            end

            lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
            C_loc_i  =  1.0im * conj(mel) * exp(lnψ_i - lnψ)
            for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
              LLOave .+= C_loc_i .* _∇lnψ
            end
            C_loc  += C_loc_i
        end
    end

    # L rho Ldag H #ok
    # -im ⟨σ|L ρ Lᴴ|σt⟩
    for L=L_ops
        #diffs_r = row_valdiff(L, row(𝝝)) # TODO Not allocate!
        #diffs_c = row_valdiff(L, col(𝝝))
        for op_r=operators(L)
            r_r=local_index(row(𝝝), sites(op_r))
            for op_c=operators(L)
                r_c=local_index(col(𝝝), sites(op_c))

                for (mel_r, changes_r)=op_r.op_conns[r_r]
                    set_index!(𝝝p_row, index(row(𝝝)))
                    for (site,val)=changes_r
                        setat!(𝝝p_row, site, val)
                    end

                    for (mel_c, changes_c)=op_c.op_conns[r_c]
                        set_index!(𝝝p_col, index(col(𝝝)))
                        for (site,val)=changes_c
                            setat!(𝝝p_col, site, val)
                        end

                        lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
                        C_loc_i  =  (mel_r) * conj(mel_c) *  exp(lnψ_i - lnψ)

                        for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
                          LLOave .+= C_loc_i .* _∇lnψ
                        end
                        C_loc  += C_loc_i
                    end
                end
            end
        end
    end

    return C_loc
end

# pretty printing
Base.show(io::IO, p::LRhoKLocalOpProblem) = print(io,
    "LRhoKLocalOpProblem on space $(basis(p)) computing the variance of Lrho using the sparse liouvillian")

# Variant for when the state has a LookUpTable and resorts to computing
# only lut updates.
function compute_Cloc!(LLO_i, ∇lnψ, prob::LRhoKLocalOpProblem,
                       net::MatrixNet, 𝝝::S,
                       _lnψ=nothing, _𝝝p::NS=nothing) where {S<:LUState, NS<:Union{Nothing, S}}
    # hey
    HnH = prob.HnH
    L_ops = prob.L_ops

    for el=LLO_i
      el .= 0.0
    end
    𝝝s = state(𝝝)
    no_changes = changes(row(𝝝s))

    C_loc = zero(Complex{real(out_type(net))})

    # ⟨σ|Hρ|σt⟩ (using hermitianity of HdH)
    # TODO should be non allocating!
    diffs_hnh = row_valdiff(HnH, raw_state(row(𝝝s)))
    for (mel, changes)=diffs_hnh
        Δ_lnψ, ∇lnψ_i = Δ_logψ_and_∇logψ!(∇lnψ, net, 𝝝, changes, no_changes)

        C_loc_i  =  -1.0im * mel * exp(Δ_lnψ)
        for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
          LLOave .+= C_loc_i .* _∇lnψ
        end
        C_loc  += C_loc_i
    end

    # ⟨σ|ρHᴴ|σt⟩
    resize!(diffs_hnh, 0)
    row_valdiff!(diffs_hnh, HnH, raw_state(col(𝝝s)))
    for (mel, changes)=diffs_hnh
        Δ_lnψ, ∇lnψ_i = Δ_logψ_and_∇logψ!(∇lnψ, net, 𝝝, no_changes, changes)

        C_loc_i  =  1.0im * conj(mel) * exp(Δ_lnψ)
        for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
          LLOave .+= C_loc_i .* _∇lnψ
        end
        C_loc  += C_loc_i
    end

    # L rho Ldag H #ok
    # -im ⟨σ|L ρ Lᴴ|σt⟩
    for L=L_ops
        diffs_r = row_valdiff(L, raw_state(row(𝝝s)))
        diffs_c = row_valdiff(L, raw_state(col(𝝝s)))

        for (mel_r, changes_r)=diffs_r

            for (mel_c, changes_c)=diffs_c
                Δ_lnψ, ∇lnψ_i = Δ_logψ_and_∇logψ!(∇lnψ, net, 𝝝, changes_r, changes_c)

                C_loc_i  =  (mel_r) * conj(mel_c) *  exp(Δ_lnψ)
                for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
                  LLOave .+= C_loc_i .* _∇lnψ
                end
                C_loc  += C_loc_i
            end
        end
    end

    return C_loc
end


text:
export detect, ignore, fmap, take, keep, drop, cut, distinct, span, merge, zip, sum, average, max, min, count

using DataStructures
using Base.Threads: Atomic, atomic_sub!

# / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /
#
# Various types of operators for creating new `Observable`s
#
# / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /

"""
    detect(fn::Function)

Apply a filter function such that only values for which the function returns
true will be passed onto [`Observer`](@ref)s. See [`ignore`](@ref).
"""
function detect(fn::Function)
  react() do observers, value
    if fn(value)
      evt = ValueEvent(value)
      notify!(observers, evt)
    end
  end
end

"""
    ignore(fn::Function)::Reactor

Apply a filter function such that only values for which the function returns
false will be passed onto [`Observer`](@ref)s. See [`detect`](@ref).
"""
function ignore(fn::Function)::Reactor
  react() do observers, value
    if !fn(value)
      evt = ValueEvent(value)
      notify!(observers, evt)
    end
  end
end

"""
    fmap(fn::Function)

Apply function to all observed values, and emit the result. Functional equivalent of `map`.

"""
function fmap(fn::Function)
  react() do observers, value
    evt = ValueEvent(fn(value))
    notify!(observers, evt)
  end
end

"""
    take(n)

Take only the first n values, discarding the rest. If less than n values observed,
emit only values observed. See [`keep`](@ref).
"""
function take(n)
  let counter = n
    react() do observers, value
      counter -= 1
      if counter >= 0
        notify!(observers, ValueEvent(value))
      end
    end
  end
end

"""
    keep(n)

Keep only the last n values, discarding the rest. If less than n values observed,
emit only values observed. See [`take`](@ref).
"""
function keep(n)
  let backlog = Queue{Any}()
    dispatch() do observers, evt
      if isa(evt, ValueEvent)
        enqueue!(backlog, evt)
        if length(backlog) > n
          dequeue!(backlog)
        end
      elseif isa(evt, CompletedEvent)
        for evt in backlog
          notify!(observers, evt)
        end
        complete!(observers)
      else
        notify!(observers, evt)
      end
    end
  end
end

"""
    drop(n)

Drop the first n events observed, emitting all others that precede them. If less than n
events observed, emit nothing. See [`cut`](@ref).
"""
function drop(n)
  let counter = n
    react() do observers, value
      if counter > 0
        counter -= 1
      else
        notify!(observers, ValueEvent(value))
      end
    end
  end
end

"""
    cut(n)

Drop the last n events observed, emitting all others that precede them. If less than n
events observed, emit nothing. See [`drop`](@ref)
"""
function cut(n)
  let backlog = Queue{Any}()
    counter = n
    react() do observers, value
      enqueue!(backlog, value)
      if counter > 0
        counter -= 1
      else
        val = dequeue!(backlog)
        notify!(observers, ValueEvent(val))
      end
    end
  end
end

"""
    distinct()

Only emit unique values observed
"""
function distinct()
  let values = Set([])
    react() do observers, value
      if !(value in values)
        union!(values, Set([value]))
        notify!(observers, ValueEvent(value))
      end
    end
  end
end

"""
    span(m, n)

Emit only the integers in the range of m to n, inclusive. Effectively an explicit alternative
to using range literal `m:n`.
"""
function span(m, n)
  collect(range(m, length = (n - m + 1)))
end

struct Merge
  closeCount::Atomic{Int}
  collector::Collector
end

Merge(observables) = Merge(Atomic{Int}(length(observables) - 1), collector())

function onEvent(merger::Merge, event::Event)
  if isa(event, CompletedEvent)
    count = atomic_sub!(merger.closeCount, 1)
    if count <= 0
      put!(merger.collector.events, event)
      close(merger.collector.events)
    end
  else
    put!(merger.collector.events, event)
    if isa(event, ErrorEvent)
      close(merger.collector.events)
    end
  end
end

function subscribe!(merger::Merge, observer)
  subscribe!(merger.collector, observer)
end

"""
    merge(observables...)

Produce an [`Observable`](@ref) that emits from the supplied [`Observable`](@ref)s,
until either all have emitted a [`CompletedEvent`](@ref) or one of them emits
an [`ErrorEvent`](@ref)
"""
function merge(observables...)
  let merger = Merge(observables)
    for observable in observables
      @async begin
        subscribe!(observable, merger)
      end
    end
    merger
  end
end

"""
    zip(observables...)

Given one or more [`Observable`](@ref)s, read one value at a time from each
and collect into a `Tuple`, stopping when either there is an [`ErrorEvent`](@ref)
or any [`CompletedEvent`](@ref)
"""
function zip(observables...)
  collectors = map(observables) do observable
    @rx() do
      observable
    end
  end
  Base.Iterators.zip(collectors...)
end

"""
    sum()

Compute a sum of all values observed.
"""
function sum()
  let total = 0
    dispatch() do observers::Observers, event::Event
      if isa(event, ValueEvent)
        total += event.value
      elseif isa(event, CompletedEvent)
        notify!(observers, ValueEvent(total))
        complete!(observers)
      else
        notify!(observers, event)
      end
    end
  end
end

"""
    average()

Compute an average of all values observed.
"""
function average()
  let total = 0
    count = 0
    dispatch() do observers::Observers, event::Event
      if isa(event, ValueEvent)
        total += event.value
        count += 1
      elseif isa(event, CompletedEvent)
        notify!(observers, ValueEvent(total / count))
        complete!(observers)
      else
        notify!(observers, event)
      end
    end
  end
end

"""
    max()

Compute the max of all values observed.
"""
function max()
  let maxSoFar = -Inf
    dispatch() do observers::Observers, event::Event
      if isa(event, ValueEvent)
        maxSoFar = Base.max(maxSoFar, event.value)
      elseif isa(event, CompletedEvent)
        notify!(observers, ValueEvent(maxSoFar))
        complete!(observers)
      else
        notify!(observers, event)
      end
    end
  end
end

"""
    min()

Compute min of all values observed.
"""
function min()
  let minSoFar = Inf
    dispatch() do observers::Observers, event::Event
      if isa(event, ValueEvent)
        minSoFar = Base.min(minSoFar, event.value)
      elseif isa(event, CompletedEvent)
        notify!(observers, ValueEvent(minSoFar))
        complete!(observers)
      else
        notify!(observers, event)
      end
    end
  end
end

"""
    count()

Count the number of values observed
"""
function count()
  let count = 0
    dispatch() do observers::Observers, event::Event
      if isa(event, ValueEvent)
        count += 1
      elseif isa(event, CompletedEvent)
        notify!(observers, ValueEvent(count))
        complete!(observers)
      else
        notify!(observers, event)
      end
    end
  end
end


text:
#####
##### AlphaZero.jl
##### Jonathan Laurent, Carnegie Mellon University (2019-2020)
#####

module AlphaZero

# Submodules
export MCTS, MinMax, GameInterface, GI, Report, Network, Benchmark
# AlphaZero parameters
export Params, SelfPlayParams, LearningParams, ArenaParams
export MctsParams, MemAnalysisParams
export SamplesWeighingPolicy, CONSTANT_WEIGHT, LOG_WEIGHT, LINEAR_WEIGHT
export AbstractSchedule, ConstSchedule, PLSchedule, StepSchedule
# Players and games
export AbstractGame, AbstractPlayer
export think, select_move, reset_player!, player_temperature, apply_temperature
export play_game, interactive!
export MctsPlayer, RandomPlayer, EpsilonGreedyPlayer, NetworkPlayer, Human
export ColorPolicy, ALTERNATE_COLORS, BASELINE_WHITE, CONTENDER_WHITE
# Networks
export AbstractNetwork, OptimiserSpec, Nesterov, CyclicNesterov, Adam
export SimpleNet, SimpleNetHP, ResNet, ResNetHP
# Training environments
export Env, train!, get_experience
# User interface
export UserInterface
export Session, resume!, save, play_interactive_game
export Explorer, start_explorer

include("util.jl")
import .Util
using .Util: Option, @unimplemented, apply_temperature

include("game.jl")
using .GameInterface
const GI = GameInterface

include("mcts.jl")
import .MCTS

include("networks/network.jl")
using .Network

using Base: @kwdef
using DataStructures: CircularBuffer
using Distributions: Categorical, Dirichlet
using Statistics: mean

include("schedule.jl")
include("params.jl")
include("report.jl")
include("memory.jl")
include("learning.jl")
include("play.jl")
include("training.jl")

include("minmax.jl")
import .MinMax

include("benchmark.jl")
using .Benchmark

# We provide a library of predefined network architectures for convenience.
# Right now, it is included in the main AlphaZero.jl package. In the future,
# we may want to separate it so as to drop the Knet and Flux dependencies.

const USE_KNET_FOR_NETLIB = true # The Flux netlib is currently broken

if USE_KNET_FOR_NETLIB
  @eval begin
    include("networks/knet.jl")
    using .KNets
  end
else
  @eval begin
    include("networks/flux.jl")
    using .FluxNets
  end
end

# The default user interface is included here for convenience but it could be
# replaced or separated from the main AlphaZero.jl package (which would also
# enable dropping some dependencies such as Crayons or JSON3).

include("ui/ui.jl")
using .UserInterface

end


text:
# Declaration order is important here, some types depend on eachother

type OutlookAuthentication
    accesstoken::AbstractString
    refreshtoken::AbstractString

    OutlookAuthentication(accesstoken) = new(accesstoken)
    OutlookAuthentication(accesstoken, refreshtoken) = new(accesstoken, refreshtoken)
end

type OutlookMessage
    id::AbstractString
    subject::AbstractString
    draft::Bool
    hasAttachments::Bool
    fromAddress::AbstractString
    fromName::AbstractString
    body::AbstractString
    isread::Bool
    importance::AbstractString
end

type OutlookFolder
    id::AbstractString
    name::AbstractString
    parentFolder::AbstractString
    itemCount::Integer
    unreadItemCount::Integer
    childFolderCount::Integer
    messages::Array{OutlookMessage, 1}
end

type OutlookClient
    email::AbstractString
    clientid::AbstractString
    clientsecret::AbstractString
    authentication::OutlookAuthentication
    folders::Array{OutlookFolder, 1}
end

immutable OutlookException <: Exception
    msg::AbstractString
end


text:
mutable struct Particle{M,T}

    x::SVector{M,T} # Current Position
    v::SVector{M,T} # Current Velocity
    p::SVector{M,T} # Best Position

    fx::T # Current Objective Function Value
    fp::T # Best Objective Function Value

end

function Particle(x::SVector{M,T}, v::SVector{M,T}) where {M,T}
    Particle(x, v, NaN.*x, Inf*x[1], Inf*x[1])
end

function Particle(x::SVector{M,T}) where {M,T}
    Particle(x, NaN.*x, NaN.*x, Inf*x[1], Inf*x[1])
end

function Particle(M)
    x = SVector{M}(NaN.*zeros(M))
    Particle(x, NaN.*x, NaN.*x, Inf*x[1], Inf*x[1])
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
function write_asdf( hdf_out::String, S::GphysData, chan_numbers::Array{Int64,1} ;
  add           ::Bool      = false,            # add traces
  evid          ::String    = "",               # event ID
  ovr           ::Bool      = false,            # overwrite trace data
  len           ::Period    = Day(1),           # length of added traces
  tag           ::String    = "",               # trace tag
  v             ::Int64     = KW.v              # verbosity
  )

  # "add" implies "ovr"
  if add == true
    ovr = true
  end

  # Precheck for degenerate time structs, time shift
  if ovr
    for i in chan_numbers
      if (size(S.t[i],1) < 2) && (S.fs[i] != 0.0)
        error(string(S.id[i] *
          ": malformed :t; can't write with ovr=true."))
      end

      # ensure each channel starts on an exact sample
      t0 = S.t[i][1,2]
      Δ = round(Int64, 1.0e6/S.fs[i])
      t1 = div(t0,Δ)*Δ
      δ = t0-t1
      S.t[i][1,2] = t1
      S.misc[i]["tc"] = δ
    end
  end

  xml_buf = IOBuffer()
  netsta, cha, nsid = mk_netsta(S)

  if isfile(hdf_out)
    io = h5open(hdf_out, "r+")
    fmt = read(attrs(io)["file_format"])
    (fmt == "ASDF") || (close(io); error("invalid ASDF file!"))
    if has(io, "Waveforms")
      wav = io["Waveforms"]
    else
      wav = g_create(io, "Waveforms")
    end
  else
    io = h5open(hdf_out, "cw")
    attrs(io)["file_format"] = "ASDF"
    attrs(io)["file_format_version"] = "1.0.2"
    wav = g_create(io, "Waveforms")
  end

  # =======================================================================
  # Add "empty" traces of all NaNs
  if add
    nc = length(chan_numbers)
    ts = zeros(Int64, nc)
    te = zeros(Int64, nc)
    for (i,j) in enumerate(chan_numbers)
      ts[i] = S.t[j][1,2]*1000
      te[i] = endtime(S.t[j], S.fs[j])*1000
    end
    p = convert(Nanosecond, len).value
    asdf_mktrace(S, xml_buf, chan_numbers, wav, ts, te, p, v, tag)
  end

  # write channels to net.sta waveform groups
  for j in 1:length(nsid)
    (v > 0) && println("writing ", nsid)
    id = nsid[j]
    chans = Int64[]
    for i in chan_numbers
      if (netsta[i] == id) && (S.fs[i] != 0.0)
        push!(chans, i)
      end
    end
    nc = length(chans)

    # does net.sta exist?
    if has(wav, id)
      sta = wav[id]

      if ovr
        # ====================================================================
        # overwrite StationXML
        (v > 1) && println("merging XML")

        # read old XML
        SX = SeisData()
        sxml = String(UInt8.(read(sta["StationXML"])))
        read_station_xml!(SX, sxml, msr=true)

        # merge S headers into SX, overwriting SX
        SM = SeisData(length(chan_numbers))
        for f in (:id, :name, :loc, :fs, :gain, :resp, :units)
          setfield!(SM, f, deepcopy(getindex(getfield(S, f), chan_numbers)))
        end
        sxml_mergehdr!(SX, SM, nofs=true, app=false)

        # remake channel list; SX ordering differs from S
        cc = Int64[]
        for i in 1:SX.n
          idx = split_id(SX.id[i])
          ns = idx[1]*"."*idx[2]
          if ns == id
            push!(cc, i)
          end
        end
        o_delete(sta, "StationXML")
        asdf_sxml(xml_buf, SX, cc, sta)

        # ==================================================================
        # overwrite trace data
        trace_ids = S.id[chans]
        t = Array{Array{Int64,2},1}(undef,nc)
        for (i,j) in enumerate(chans)
          t[i] = t_win(S.t[j], S.fs[j])
          broadcast!(*, t[i], t[i], 1000)
        end
        (v > 2) && println("t = ", t)

        # loop over trace waveforms using id, start time, end time
        for n in names(sta)
          (n == "StationXML") && continue
          cha_id = String(split(n, "_", limit=2, keepempty=true)[1])
          x = sta[n]
          (v > 1) && println("checking ", n)
          nx = length(x)
          t0 = read(x["starttime"])
          fs = read(x["sampling_rate"])
          (v > 2) && println("t0 = ", t0, "; fs = ", fs, "; nx = ", nx)

          # convert fs to sampling interval in ns
          Δ = round(Int64, 1.0e9/fs)
          t1 = t0 + (nx-1)*Δ

          for i in 1:length(chans)
            (trace_ids[i] == cha_id) || continue
            nk = size(t[i], 1)

            for k in 1:nk
              (v > 2) && println("segment k = ", k, "/", nk)
              # overlap
              ts = t[i][k,1]
              te = t[i][k,2]
              if (ts ≤ t1) && (te ≥ t0)
                # set channel index j in S
                j = chans[i]
                (v > 2) && println("j = ", j)

                # check for fs mismatch
                trace_fs = S.fs[j]
                (trace_fs == fs) || (@warn(string("Can't write ", S.id[j], "; fs mismatch!")); continue)
                lx = div(te-ts, Δ)+1

                # determine start, end indices in x that are overwritten
                (v > 2) && println("ts = ", ts, "; te = ", te, "; t0 = ", t0, "; t1 = ", t1, "; nx = ", nx, "; lx = ", lx)
                i0, i1, t2 = get_trace_bounds(ts, te, t0, t1, Δ, nx)

                # determine start, end indices in X to copy
                si = get_trace_bound(t0, ts, Δ, lx)
                ei = si + min(i1-i0, lx-1)

                # overwrite
                if v > 2
                  println("writing ", si, ":", ei, " to ", i0, ":", i1)
                end
                save_data!(S.x[j], x, i0:i1, si:ei)
              end
            end
          end

        end
      else
        for i in chans
          asdf_write_chan(S, sta, i, cha[i], evid, v)
        end
      end

    elseif ovr == false
      # Create Waveforms/net.sta
      sta = g_create(wav, id)

      # Create Waveforms/net.sta/chan_str
      for i in chans
        asdf_write_chan(S, sta, i, isempty(tag) ? cha[i] : tag, evid, v)
      end

      # Write StationXML to sta
      asdf_sxml(xml_buf, S, chans, sta)
    end
  end

  close(xml_buf)
  close(io)

  # Correct :t
  if ovr
    for i in chan_numbers
      δ = get(S.misc[i], "tc", 0)
      t0 = S.t[i][1,2]
      S.t[i][1,2] = t0+δ
      delete!(S.misc[i], "tc")
    end
  end
  return nothing
end