text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA
import ParallelStencil.ParallelKernel: @require, longnameof
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. allocator macros" begin
            @require !@is_initialized()
            @init_parallel_kernel($package, Float16)
            @testset "mapping to package" begin
                @require @is_initialized()
                @test @zeros(2,3) == parentmodule($package).zeros(Float16,2,3)
                @test @ones(2,3) == parentmodule($package).ones(Float16,2,3)
                @static if $package == CUDA
                    @test typeof(@rand(2,3)) == typeof(CUDA.CuArray(rand(Float16,2,3)))
                else
                    @test typeof(@rand(2,3)) == typeof(parentmodule($package).rand(Float16,2,3))
                end
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;


text:
#
# Linear Algebra of order 2 Tensors
#
# Even though CuDenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray

function Base.:*(
  T1::Tensor{ElT1,2,StoreT1,IndsT1}, T2::Tensor{ElT2,2,StoreT2,IndsT2}
) where {ElT1,StoreT1<:CuDense,IndsT1,ElT2,StoreT2<:CuDense,IndsT2}
  RM = matrix(T1) * matrix(T2)
  indsR = IndsT1(ind(T1, 1), ind(T2, 2))
  pT = promote_type(ElT1, ElT2)
  return Tensor(Dense(vec(RM)), indsR)
end
#= FIX ME
function LinearAlgebra.exp(T::CuDenseTensor{ElT,2}) where {ElT,IndsT}
  expTM = exp(matrix(T))
  return Tensor(Dense(vec(expTM)),inds(T))
end

function expHermitian(T::CuDenseTensor{ElT,2}) where {ElT,IndsT}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(Hermitian(matrix(T))))
  return Tensor(Dense(vec(expTM)),inds(T))
end
=#

# svd of an order-2 tensor
function LinearAlgebra.svd(T::CuDenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
  mindim::Int = get(kwargs, :mindim, 1)
  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
  absoluteCutoff::Bool = get(kwargs, :absoluteCutoff, false)
  doRelCutoff::Bool = get(kwargs, :doRelCutoff, true)
  fastSVD::Bool = get(kwargs, :fastSVD, false)
  aT = array(T)
  @timeit "CUSOLVER svd" begin
    MU, MS, MV = CUSOLVER.svd!(aT)
  end
  # for consistency with cpu version,
  # ITensors.jl/NDTensors/src/linearalgebra.jl/svd
  # need conj!(MV)
  conj!(MV)
  P = MS .^ 2
  truncerr, docut, P = truncate!(
    P;
    mindim=mindim,
    maxdim=maxdim,
    cutoff=cutoff,
    absoluteCutoff=absoluteCutoff,
    doRelCutoff=doRelCutoff,
  )
  spec = Spectrum(P, truncerr)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:, 1:dS]
    MS = MS[1:dS]
    MV = MV[:, 1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T, 1), u))
  Sinds = IndsT((u, v))
  Vinds = IndsT((ind(T, 2), v))
  U = Tensor(Dense(vec(MU)), Uinds)
  Sdata = CUDA.zeros(ElT, dS * dS)
  dsi = diagind(reshape(Sdata, dS, dS), 0)
  Sdata[dsi] = MS
  S = Tensor(Dense(Sdata), Sinds)
  V = Tensor(Dense(vec(MV)), Vinds)
  return U, S, V, spec
end

function LinearAlgebra.eigen(
  T::Hermitian{ElT,<:CuDenseTensor{ElT,2,IndsT}}; kwargs...
) where {ElT<:Union{Real,Complex},IndsT}
  ispossemidef::Bool = get(kwargs, :ispossemidef, false)
  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
  mindim::Int = get(kwargs, :mindim, 1)
  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
  absoluteCutoff::Bool = get(kwargs, :absoluteCutoff, false)
  doRelCutoff::Bool = get(kwargs, :doRelCutoff, true)
  @timeit "CUSOLVER eigen" begin
    local DM, UM
    if ElT <: Complex
      DM, UM = CUSOLVER.heevd!('V', 'U', matrix(parent(T)))
    else
      DM, UM = CUSOLVER.syevd!('V', 'U', matrix(parent(T)))
    end
  end
  DM_ = reverse(DM)
  @timeit "truncate" begin
    truncerr, docut, DM = truncate!(
      DM_;
      maxdim=maxdim,
      cutoff=cutoff,
      absoluteCutoff=absoluteCutoff,
      doRelCutoff=doRelCutoff,
    )
  end
  spec = Spectrum(DM, truncerr)
  dD = length(DM)
  dV = reverse(UM; dims=2)
  if dD < size(dV, 2)
    dV = CuMatrix(dV[:, 1:dD])
  end
  # Make the new indices to go onto U and V
  l = eltype(IndsT)(dD)
  r = eltype(IndsT)(dD)
  Vinds = IndsT((dag(ind(T, 2)), dag(r)))
  Dinds = IndsT((l, dag(r)))
  U = Tensor(Dense(vec(dV)), Vinds)
  D = Tensor(Diag(real.(DM)), Dinds)
  return D, U, spec
end

function LinearAlgebra.qr(T::CuDenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  QM, RM = qr(matrix(T))
  # Make the new indices to go onto Q and R
  q, r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T, 1), q))
  Rinds = IndsT((q, ind(T, 2)))
  QM = CuMatrix(QM)
  Q = Tensor(Dense(vec(QM)), Qinds)
  R = Tensor(Dense(vec(RM)), Rinds)
  return Q, R
end

function polar(T::CuDenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM, RM = polar(matrix(T))
  dim = size(QM, 2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T, 1), q))
  Rinds = IndsT((q, ind(T, 2)))
  Q = Tensor(Dense(vec(QM)), Qinds)
  R = Tensor(Dense(vec(RM)), Rinds)
  return Q, R
end

ext: jl


text:
######################################################
# functions related to computation of reduced hessian
######################################################

# Q xx' skron yy' Qt
function skron(
  x::AbstractVector,
  y::AbstractVector,
  indRow::AbstractVector,
  indCol::AbstractVector)

  nRow = length(indRow)
  nCol = length(indCol)
  out = zeros( nRow, nCol )
  for ri=1:nRow
    for ci=1:nCol
      out[ri, ci] = skron(x, y, indRow[ri], indCol[ci])
    end
  end
  out
end

function skron(
  x::AbstractVector,
  y::AbstractVector,
  indRow::Int,
  indCol::Int
  )

  p = length(x)
  row1, col1 = ind2subLowerTriangular(p, indRow)
  row2, col2 = ind2subLowerTriangular(p, indCol)

  vrow = (row1 == col1) ? x[row1] * y[row1] : (x[row1] * y[col1] + y[row1] * x[col1]) / sqrt(2.)
  vcol = (row2 == col2) ? x[row2] * y[row2] : (x[row2] * y[col2] + y[row2] * x[col2]) / sqrt(2.)

  vrow * vcol
end

#
# Q Sx skron Sy Qt
function skron(
  Sx::Symmetric,
  Sy::Symmetric)

  p = size(Sx, 1)
  return skron(Sx, Sy, 1:div(p*(p+1),2))
end

function skron(
  Sx::Symmetric,
  Sy::Symmetric,
  indRow::AbstractVector,
  indCol::AbstractVector)

  nRow = length(indRow)
  nCol = length(indCol)
  out = zeros( nRow, nCol )
  for ri=1:nRow
    for ci=1:nCol
      out[ri, ci] = skron(Sx, Sy, indRow[ri], indCol[ci])
    end
  end
  out
end

function skron(
  Sx::Symmetric,
  Sy::Symmetric,
  I::AbstractVector)

  m = length(I)

  out = zeros( m, m )

  for ci = 1:m
    for ri = ci:m
      out[ri, ci] = skron(Sx, Sy, I[ri], I[ci])
      if ri != ci
        out[ci, ri] = out[ri, ci]
      end
    end
  end
  out
end

function skron(
  Sx::Symmetric,
  Sy::Symmetric,
  ind_row::Int,
  ind_col::Int
  )

  p = size(Sx, 1)
  row1, col1 = ind2subLowerTriangular(p, ind_row)
  row2, col2 = ind2subLowerTriangular(p, ind_col)

  if row1 == col1
    if row2 == col2
      return Sx[row1, row2] * Sy[row1, row2]
    else
      return ( Sx[row1, row2] * Sy[row1, col2] + Sx[row1, col2] * Sy[row1, row2] ) / sqrt(2.)
    end
  else
    if row2 == col2
      return ( Sx[row1, row2] * Sy[col1, row2] + Sx[col1, row2] * Sy[row1, row2] ) / sqrt(2.)
    else
      return (Sx[row1, row2] * Sy[col1, col2] +
              Sx[row1, col2] * Sy[col1, row2] +
              Sx[col1, row2] * Sy[row1, col2] +
              Sx[col1, col2] * Sy[row1, row2] ) / 2.
    end
  end
end

#
# Q Sx skron yy' Qt
function skron(
  Sx::Symmetric,
  y::AbstractVector,
  indRow::AbstractVector,
  indCol::AbstractVector)

  nRow = length(indRow)
  nCol = length(indCol)
  out = zeros( nRow, nCol )
  for ri=1:nRow
    for ci=1:nCol
      out[ri, ci] = skron(Sx, y, indRow[ri], indCol[ci])
    end
  end
  out
end

function skron(
  Sx::Symmetric,
  y::AbstractVector,
  ind_row::Int,
  ind_col::Int
  )

  p = size(Sx, 1)
  row1, col1 = ind2subLowerTriangular(p, ind_row)
  row2, col2 = ind2subLowerTriangular(p, ind_col)

  if row1 == col1
    if row2 == col2
      return Sx[row1, row2] * y[row1] * y[row2]
    else
      return ( Sx[row1, row2] * y[row1] * y[col2] + Sx[row1, col2] * y[row1] * y[row2] ) / sqrt(2.)
    end
  else
    if row2 == col2
      return ( Sx[row1, row2] * y[col1] * y[row2] + Sx[col1, row2] * y[row1] * y[row2] ) / sqrt(2.)
    else
      return (Sx[row1, row2] * y[col1] * y[col2] +
              Sx[row1, col2] * y[col1] * y[row2] +
              Sx[col1, row2] * y[row1] * y[col2] +
              Sx[col1, col2] * y[row1] * y[row2] ) / 2.
    end
  end
end

######################################################
# functions related to application of Q to Symmetric matrix
######################################################

# Q(S) == svec(S)

function svec(S::Symmetric)
  p = size(S, 1)
  return svec(S, 1:div(p*(p+1),2))
end

function svec(
  S::Symmetric,
  I::AbstractVector)

  m = length(I)
  out = zeros( m )
  for ri = 1:m
    out[ri] = svec(S, I[ri])
  end
  out
end

function svec(
  S::Symmetric,
  indRow::Int)

  row, col = ind2subLowerTriangular(size(S, 1), indRow)
  (row == col) ? S[row, row] : S[row, col] * sqrt(2.)
end

function svec(
  x::AbstractVector,
  y::AbstractVector,
  I::AbstractVector)

  m = length(I)
  out = zeros( m )
  for ri = 1:m
    out[ri] = svec(x, y, I[ri])
  end
  out
end

function svec(
  x::AbstractVector,
  y::AbstractVector,
  indRow::Int
  )

  row, col = ind2subLowerTriangular(length(x), indRow)

  if row == col
    return x[row] * y[row]
  else
    return (x[row] * y[col] + y[row] * x[col]) / sqrt(2.)
  end
end


text:
import Base: size, getindex, setindex!, length, IndexStyle, +, -, getfield

# Parted Array Type
struct PartedArray{T,N,M<:AbstractArray{T,N},P} <: AbstractArray{T,N}
    A::M
    parts::P
    function PartedArray(A::AbstractArray{T,N}, partition::P) where {T,N,I,P<:Dict{Symbol,I}}
        new{T,N,typeof(A),P}(A, partition)
    end
    function PartedArray(A::AbstractArray{T,N}, partition::P) where {T,N,P<:NamedTuple}
        new{T,N,typeof(A),P}(A, partition)
    end
end

PartedVector{T,P} = PartedArray{T,1,P}
PartedMatrix{T,P} = PartedArray{T,2,P}

# Constructors
function PartedVector(A::AbstractVector, lengths::NTuple{N,Int}, names::NTuple{N,Symbol}) where N
    part = Dict(zip(names, create_partition(lengths)))
    PartedArray(A, part)
end
PartedVector(A::AbstractVector, lengths::Vector{Int}, names::Vector{Symbol}) = PartedVector(A, Tuple(lengths), Tuple(names))
function PartedVector(A::AbstractVector, lengths::NTuple{N,Int}, ::Val{names}) where {N,names}
    part = NamedTuple{names}(create_partition(lengths))
    PartedArray(A, part)
end
function PartedVector(A::AbstractVector, part::NamedTuple)
    PartedArray(A, part)
end

function PartedMatrix(A::AbstractMatrix, lengths::NTuple{N,Int}, names::NTuple{N,Symbol}) where N
    part = create_partition2(lengths, lengths)
    names = combine_names(names, names)
    partition = Dict(zip(names, part))
    PartedArray(A, partition)
end
PartedMatrix(A::AbstractMatrix, lengths::Vector{Int}, names::Vector{Symbol}) = PartedVector(A, Tuple(lengths), Tuple(names))
function PartedMatrix(A::AbstractMatrix, lengths::NTuple{N,Int}, ::Val{names}) where {N,names}
    part = create_partition2(lengths, lengths)
    partition = NamedTuple{names}(part)
    PartedArray(A, partition)
end
function PartedMatrix(A::AbstractMatrix, part::NamedTuple)
    PartedArray(A, part)
end

# Basic Methods
size(A::PartedArray) = size(A.A)
getindex(A::PartedArray, i::Int) = getindex(A.A, i)
getindex(A::PartedArray, I::Vararg{Int,2}) = A.A[I[1], I[2]]
getindex(A::PartedArray, I...) = getindex(A.A, I...)
setindex(A::PartedArray, v, I...) = setindex(A.A, v, I...)
setindex!(A::PartedArray, v, i::Int) = setindex!(A.A, v, i)
setindex!(A::PartedArray, v, I::Vararg{Int, 2}) = A.A[I[1], I[2]] = v
IndexStyle(A::PartedArray) = IndexCartesian()
length(A::PartedArray) = length(A.A)
Base.keys(A::PartedArray) = keys(A.part)
Base.copy(A::PartedArray) = PartedArray(copy(A.A),A.parts)
Base.zero(A::PartedArray) = PartedArray(zero(A.A), A.parts)
+(A::PartedArray, B::AbstractArray) = PartedArray(A.A + B, A.parts)
+(B::AbstractArray, A::PartedArray) = PartedArray(B + A.A, A.parts)
+(A::PartedArray, B::PartedArray) = PartedArray(A.A + B.A, merge(A.parts, B.parts))
+(A::PartedArray, b::Real) = PartedArray(A.A .+ b, A.parts)
getindex(A::PartedArray, p::Symbol) = view(getfield(A,:A), getfield(A,:parts)[p]...)
getindex(A::PartedVector, p::Symbol) = view(getfield(A,:A), getfield(A,:parts)[p])
function Base.getproperty(A::PartedArray{T,N,I}, p::Symbol) where {T,N,I}
    if p == :A || p == :parts
        getfield(A,p)
    else
        A[p]
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -2.8031108847289734 <= q <= 2.6473335630382056)

                     add_NL_constraint(m, :(sigmoid(-0.8475658981369358 + 0.9991334286088005*$(x[1]) + 0.6614487188155533*$(x[2])) + sigmoid(0.7696772372915519 + 0.2658705903136487*$(x[1]) + 0.798769486145587*$(x[2])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ReferringExpressions
using Base.Test

@test try
    refexp = RefExpData("refclef", "unc")
    true
catch err
    println("ERROR: $(err.msg)")
    catch_stacktrace()
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

if !isdefined(@__MODULE__(), :MyGame) @__MODULE__().eval(:(module MyGame __precompile__(false); import FlatBuffers end)) end
if !isdefined(MyGame, :OtherNameSpace) MyGame.eval(:(module OtherNameSpace __precompile__(false); import FlatBuffers end)) end
if !isdefined(MyGame, :Example2) MyGame.eval(:(module Example2 __precompile__(false); import FlatBuffers end)) end
if !isdefined(MyGame, :Example) MyGame.eval(:(module Example __precompile__(false); import FlatBuffers end)) end
include("MyGame/InParentNamespace.jl")
include("MyGame/Example2/Monster.jl")
include("MyGame/Example/Color.jl")
include("MyGame/Example/Race.jl")
include("MyGame/Example/Test.jl")
include("MyGame/Example/Vec3.jl")
include("MyGame/Example/Stat.jl")
include("MyGame/Example/Ability.jl")
include("MyGame/Example/Referrable.jl")
include("MyGame/Example/TestSimpleTableWithEnum.jl")
include("MyGame/Example/AnyUniqueAliases.jl")
include("MyGame/Example/AnyAmbiguousAliases.jl")
include("MyGame/Example/Monster.jl")
include("MyGame/Example/Any_.jl")
include("MyGame/Example/TypeAliases.jl")

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
function is_equilateral(sides)
end

function is_isosceles(sides)
end

function is_scalene(sides)
end


text:
## Fractional Calculus utilities

function df(f, h)
    # Simple derivative of function f
    function _f(x)
        (f(x + h) - f(x - h)) / (2 * h)
    end;
    _f
end;

max_stars_repo_name: pedromxavier/fraccalc

max_issues_repo_name: pedromxavier/fraccalc

max_forks_repo_name: pedromxavier/fraccalc


text:
function scale!{N}(c::Composable, val)
    c[Scale] = val
    return
end
function translate!{N}(c::Composable, val)
    c[Position] = val
    return
end

# TODO allow also Vec2, Vec3, Axis + Float.
# This will be done by implementing symmetric transforms on the Rotation field type
function rotate!(c::Composable, val)
    c[Rotation] = val
    return
end


text:
import Base: show

# Commands
export AbstractCommand, LeafCommand, NodeCommand, MainCommand

const DocString = Union{String, Markdown.MD}

"""
Abstract type for Commands

## Interface
- `name(cmd)`: returns the name of this cmd.
- `docstring`: returns the docstring
- all subtype of `AbstractCommand` should be callable (like a function).
- overload `show` to support `help`
"""
abstract type AbstractCommand end

"""
    name(cmd) -> String

get a command's name
"""
name(cmd::AbstractCommand) = string(cmd.name)

"""
    docstring(cmd) -> Union{String, Markdown.MD}

get command's docstring
"""
docstring(cmd::AbstractCommand) = cmd.doc

"""
    findcmd(s, cmds)

find a command in `cmds` named `s`, return `nothing` if it does not
in `cmds`.
"""
function findcmd(s::AbstractString, cmds::Vector{AbstractCommand})
    for cmd in cmds
        if s == name(cmd)
            return cmd
        end
    end
    return nothing
end

"""
default command call. accepts an option `:help` and exits when receives
option `:help`.
"""
function eval_default_command(cmd::AbstractCommand;kwargs...)
    for (key, val) in kwargs
        if key == :help
            print(cmd)
            exit()
        end
    end
end

"""
The last sub-command in a CLI
"""
mutable struct LeafCommand{D<:DocString} <: AbstractCommand
    name::Symbol
    callable::Function
    sig::Signature
    doc::D
end

LeafCommand(name::Symbol, callable::Function, sig::Signature) =
    LeafCommand(name, callable, sig, "No Document")

function (cmd::LeafCommand)(args...;kwargs...)
    eval_default_command(cmd;kwargs...)

    if length(args) != length(cmd.sig.args)
        throw(
            CLIError(
                "expect $(length(cmd.sig.args)) arguments, got $(length(args))"
            )
        )
    end

    cmd.callable(args...;kwargs...)
end

# TODO: use Strong for SYNOPSIS
function show(io::IO, cmd::LeafCommand)
    print(io, "[SYNOPSIS]\n")
    print(io, "\t$(cmd.name) ")
    for op in cmd.sig.options
        print(io, "[$(optionstring(op))] ")
    end

    if length(cmd.sig.args) > 0
        print(io, "<args>\n")
    end

    print(io, "[DESCRIPTION]\n")
    print(io, "\t$(cmd.doc)")
end

mutable struct NodeCommand{C<:AbstractCommand, D<:DocString} <: AbstractCommand
    name::String
    cmds::Vector{C}
    doc::D
end

NodeCommand(name::String) = NodeCommand(name, AbstractCommand[])
NodeCommand(name::String, cmds::Vector) = NodeCommand(name, cmds, "No Document")

function (cmd::NodeCommand)(;kwargs...)
    for (key, val) in kwargs
        if key == :help
            print(cmd)
            exit()
        else
            throw(CLIError("Unknown flag $(key)"))
        end
    end
end

function show(io::IO, cmd::NodeCommand)
    println(io, "[SYNOPSIS]")
    println(io, "\t$(cmd.name) <command>")

    println(io, "[DESCRIPTION]")
    println(io, "\t$(cmd.doc)")

    println(io, "[COMMANDS]")
    # TODO: print doc in one block
    for subcmd in cmd.cmds
        println(io, "\t$(name(subcmd))")
        println(io, "\t  $(docstring(subcmd))")
    end
end

"""
```shell
> julia script.jl -a -b<arg> --option --option=arg cmd1 cmd2 arg1 arg2 -c -d<arg> --option
```

`-a` is equivalent to `-atrue`
`--option` is equivalent to `--option=true`

each short option (start with `-` followed with a letter) must have
corresponding option (start with `--` followed with ascii letter and
`-`, e.g `--long-option-xxx`)

a command line interface has three part

<main command options> <command node> <command node> ... <leaf command> <leaf command args and options>

main command only accepts options
<main> is followed by some <command node> or <leaf command>, e.g

```shell
> git remote add <url>
```

here `git` is the main command name, determined by your script name
`remote` is a <command node>, and `add` is the <leaf node>

## Invalid Syntax

`-c arg` is not allowed

N: number of arguments
"""
mutable struct MainCommand{O<:AbstractOption, C<:AbstractCommand, D<:DocString} <: AbstractCommand
    name::String
    options::Vector{O}
    cmds::Vector{C}
    doc::D
    version::VersionNumber
    entry::Function
end

__EMPTY__(;) = nothing

MainCommand() = MainCommand(AbstractOption[], AbstractCommand[])
MainCommand(name::String, options::Vector, cmds::Vector) =
    MainCommand(name, [options..., help, version], cmds, "", v"0.0.1", __EMPTY__)
MainCommand(options::Vector, cmds::Vector) = MainCommand("main", options, cmds)

function eval_default_command(cmd::MainCommand;kwargs...)
    for (key, val) in kwargs
        if key == :help
            print(cmd)
            exit()
        end

        if key == :version
            print("$(cmd.name) $(cmd.version)")
            exit()
        end
    end
end

function (cmd::MainCommand)(;kwargs...)
    eval_default_command(cmd;kwargs...)
    return cmd.entry(;kwargs...)
end

function show(io::IO, cmd::MainCommand)
    println(io, "[SYNOPSIS]")

    print(io, "\t$(cmd.name) ")
    for op in cmd.options
        print(io, "[$(optionstring(op))] ")
    end

    println(io, "<command>")

    println(io, "[DESCRIPTION]")
    println(io, "\t$(cmd.doc)")

    println(io, "[COMMANDS]")
    # TODO: print doc in one block
    for subcmd in cmd.cmds
        println(io, "\t$(name(subcmd))")
        println(io, "\t  $(docstring(subcmd))")
    end
end


text:
include("../src/utils/readfiles.jl")
include("../src/setting/constants.jl")

"""
generate perturbation file according to the reference bin file
"""
function generate_perturbation(target_basedir::String, reference_basedir::String, output_basedir::String, nproc::Int64, nspec::Int64)
    tags = ["vph","vpv","vsh","vsv","eta","qmu","rho"]
    model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    for iproc in 0:nproc - 1
        if(iproc % 10 == 0)
            @info "step" iproc
        end
        for tag in tags
            sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
            sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

            model_gll_output = (model_gll_target .- model_gll_reference) ./ model_gll_reference
            @assert count(isnan.(model_gll_output)) == 0
            sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
        end
    end
end

"""
generate the real vaue bin file based on reference model and perturbation bin file
"""
function generate_real(target_basedir::String, reference_basedir::String, output_basedir::String, nproc::Int64, nspec::Int64)
    tags = ["vph","vpv","vsh","vsv","eta","qmu","rho"]
    model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    for iproc in 0:nproc - 1
        if(iproc % 10 == 0)
            @info "step" iproc
        end
        for tag in tags
            sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
            sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

            model_gll_output = model_gll_reference .* (1 .+ model_gll_target)
            @assert count(isnan.(model_gll_output)) == 0
            sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
        end
    end
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
using ApproxFun
using SpecialExponentialFamilies
using Test

@testset "SpecialExponentialFamilies.jl" begin
    include("test_twintower.jl")
end


text:
"""
    F = JopLog(spc)

where `F` is the log operator with domain and range given by `spc::JetSpace`.
"""
JopLog(spc::JetSpace) = JopNl(dom = spc, rng = spc, f! = JopLog_f!, df! = JopLog_df!, df′! = JopLog_df′!)
export JopLog

# TODO: assert denominator in linearization is > 0

JopLog_f!(d::AbstractArray, m::AbstractArray) = d .= log.(m)
JopLog_df!(δd::AbstractArray, δm::AbstractArray; mₒ) = δd .= δm ./ mₒ
JopLog_df′!(δm::AbstractArray, δd::AbstractArray; mₒ) = δm .= δd ./ conj(mₒ)


text:
export fctoresp

"""
    fctoresp(f)
    fctoresp(f, c)

Create PZResp or PZResp64 instrument response from lower corner frequency `f` and damping constant `c`. If no damping constant is supplies, assumes `c = 1/sqrt(2)`.

### See Also
PZResp, PZResp64
"""
function fctoresp(f::AbstractFloat, c::AbstractFloat=1.0f0/sqrt(2.0f0))
  T = typeof(f)
  fxp = T(2.0*pi)*f
  r = sqrt(Complex(c^2 - one(T)))
  z = zeros(Complex{T}, 1)
  p = Complex{T}[r-c, -r-c].*fxp
  if T == Float32
    return PZResp(f0 = Float32(f), p = p, z = z)
  else
    return PZResp64(f0 = Float64(f), p = p, z = z)
  end
end


text:
import Base:show
import JGraph: Edge, EdgeIter

abstract type AbstractGraph end

mutable struct Graph{T <: Integer} <: AbstractGraph
    ne::T
    adjlist::Vector{Vector{T}}

    function Graph()
        new{Int64}(0, [])
    end
end

"""
    Graph(n, m)

generate a graph with n vertices and m edges. Return a graph G.
"""
function Graph(n, m)
    g = Graph();
    add_vertex!(g, n);
    num_vertice = 0

    while num_vertice < m
        s = rand(1:n)
        t = rand(1:n)
        if s != t && !has_edge(g, s, t)
            add_edge!(g, s, t);
            num_vertice += 1
        end
    end
    return g
end

nv(G::AbstractGraph) = length(G.adjlist)

ne(G::AbstractGraph) = G.ne

"""
    has_edge(G::AbstractGraph, u::T, v::T)

Test if egde(u,v) exist in the graph G.
"""
function has_edge(G::AbstractGraph, u::T, v::T) where T <: Integer
    if max(u,v) > nv(G)
        return false
    elseif v in G.adjlist[u]
        return true
    else
        return false
    end
end

has_edge(G::AbstractGraph, e::Edge) = has_edge(G, src(e), dst(e))

"""
    add_edge(G::AbstractGraph, u::T, v::T)

Add the egde(u,v) in the graph G. Return False is the edge already exit, true otherwise.
"""
function add_edge!(G::AbstractGraph, u::T, v::T) where T <: Integer
    if has_edge(G, u, v)
        return false
    end
    if nv(G) < max(u,v)
        add_vertex!(G, max(u,v));
    end
    push!(G.adjlist[u], v);
    push!(G.adjlist[v], u);
    G.ne += 1;
    return true
end

add_edge!(G::AbstractGraph, e::Edge) = add_edge(G, src(e), dst(e))

"""
    add_edge(G::AbstractGraph, num_vertex)

Add "num_vertex" vertices to the graph G.
"""
function add_vertex!(G::AbstractGraph, num_vertex::T) where T <: Integer
    while nv(G) < num_vertex
        add_vertex!(G);
    end
    return true
end

function add_vertex!(G::AbstractGraph)
    return push!(G.adjlist, Vector{Int}());
end

function remove_edge!(G::AbstractGraph, u, v)
    if has_edge(G, u, v)
        G.adjlist[u] = deleteat!(G.adjlist[u], findall(x-> x == v, G.adjlist[u]));
        G.adjlist[v] = deleteat!(G.adjlist[v], findall(x-> x == u, G.adjlist[v]));
        return true
    else
        return false
    end
end

remove_edge!(G::AbstractGraph, e::Edge) = remove_edge!(G, src(e), dst(e))

function remove_vertex!(G::AbstractGraph, s)
    for t in neighbors(G, s)
        remove_edge!(G, s, t)
    end
    return true
end

"""
    degree(G::Graph)

Return the list of nodes degree.
"""
function degree(G::Graph)
    degree_list = []
    for u in G.adjlist
        push!(degree_list, length(u))
    end
    degree_list
end

"""
    neighbors(G::Graph, u)

return the neighbors list adjacent to the node u
"""
function neighbors(G::Graph, u)
    return copy(G.adjlist[u])
end

"""
    vertices(G::Graph)

vertices itertor
"""
vertices(G::Graph) = Base.OneTo(nv(G))

edges(g::AbstractGraph) = EdgeIter(g)

"""
    fadj(G::Graph)

return the adjacency list of the graph
"""
fadj(G::Graph) = G.adjlist
fadj(G::Graph, v::Integer) = G.adjlist[v]

show(io::IO, G::Graph) = print(io, "Undirected Graph G(V,E) with |V|=$(nv(G)), |E|=$(ne(G))");


text:
using LinearAlgebra
using Random, Distributions, Distances, StatsFuns
using Zygote

include("./kernel.jl")
include("./transform.jl")

struct Model{K<:AbstractKernel, G<:AbstractTransform}
    X::Matrix{Float64}
    Y::Vector{Float64}
    k::K
    g::G
    dims::Int
end

function make_fg(f, unpack)
    return function(F, G, x)
        v, back = pullback(x) do x
            f(unpack(x)...)
        end
        if G !== nothing
            G .= back(1.0)[1]
        end
        if F !== nothing
            return v
        end
        return nothing
    end
end

function logprob(mdl, ℓ, θ, ϵ, λ)
    n = length(mdl.Y)
    L = mdl.X / ℓ
    D = pairwise(SqEuclidean(), L; dims=mdl.dims)
    K = (τ -> radial(mdl.k, θ, τ)).(D) + max(ϵ, 1e-8) * I
    cholK = cholesky(Symmetric(K))
    Z = (y -> transform(mdl.g, λ, y)).(mdl.Y)
    lnJ = logjacobian(mdl.g, λ, mdl.Y)
    lnDet = logdet(cholK)
    return lnJ - lnDet / 2 - dot(Z, cholK \ Z) / 2 - n * log2π / 2
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    convert(ListRecord, qrec::QueueRecord{E})

Converts a QueueRecord into the corresponding ListRecord.
"""
function Base.convert{S,D,I}(::Type{ListRecord{S,D,I}}, qrec::QueueRecord{Entity{S,D,I}})

    frames = Array(RecordFrame, length(qrec))
    states = Array(RecordState{S,I}, nstates(qrec))
    defs = Dict{I, D}()

    lo = 1
    for (i,pastframe) in enumerate(1-length(qrec) : 0)
        frame = qrec[pastframe]

        hi = lo
        for entity in frame
            defs[entity.id] = entity.def
            states[hi] = RecordState{S,I}(entity.state, entity.id)
            hi += 1
        end

        frames[i] = RecordFrame(lo, hi)
        lo = hi
    end

    return ListRecord{S,D,I}(get_timestep(qrec), frames, states, defs)
end
Base.convert{S,D,I}(::Type{ListRecord}, qrec::QueueRecord{Entity{S,D,I}}) = convert(ListRecord{S,D,I}, qrec)

"""
    convert(QueueRecord, lrec::ListRecord)

Converts a ListRecord into the corresponding QueueRecord{Entity{S,D,I}}.
Note that the timesteps for a ListRecord are not necessarily constant timesteps.
"""
function Base.convert{S,D,I}(::Type{QueueRecord{Entity{S,D,I}}}, lrec::ListRecord{S,D,I})

    N = nframes(lrec)
    M = maximum(n_objects_in_frame(lrec, i) for i in 1 : N)
    retval = QueueRecord(Entity{S,D,I}, N, get_timestep(lrec), M)

    frame = Frame(Entity{S,D,I}, M)
    for i in 1 : N
        get!(frame, lrec, i)
        update!(retval, frame)
    end

    return retval
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using Test
Random.seed!(233)

function tri_solve(a,b,c,d)
    tri_solve_ = load_op_and_grad("./build/libTriSolve","tri_solve")
    a,b,c,d = convert_to_tensor(Any[a,b,c,d], [Float64,Float64,Float64,Float64])
    tri_solve_(a,b,c,d)
end

n = 10
a = rand(n-1)
b = rand(n).+10
c = rand(n-1)
d = rand(n)

A = diagm(0=>b, -1=>a, 1=>c)
x = A\d

# TODO: specify your input parameters
u = tri_solve(a,b,c,d)
sess = Session(); init(sess)
@test run(sess, u)≈x

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(tri_solve(a,b,x,d)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(n-1))
v_ = rand(n-1)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
## RESAMPLING STRATEGIES

abstract type ResamplingStrategy <: MLJType end
show_as_constructed(::Type{<:ResamplingStrategy}) = true

# resampling strategies are `==` if they have the same type and their
# field values are `==`:
function ==(s1::S, s2::S) where S <: ResamplingStrategy
    return all(getfield(s1, fld) == getfield(s2, fld) for fld in fieldnames(S))
end

# fallbacks:
train_test_pairs(s::ResamplingStrategy, rows, X, y, w) =
    train_test_pairs(s, rows, X, y)
train_test_pairs(s::ResamplingStrategy, rows, X, y) =
    train_test_pairs(s, rows, y)
train_test_pairs(s::ResamplingStrategy, rows, y) =
    train_test_pairs(s, rows)

# Helper to interpret rng, shuffle in case either is `nothing` or if
# `rng` is an integer:
function shuffle_and_rng(shuffle, rng)
    if rng isa Integer
        rng = MersenneTwister(rng)
    end

    if shuffle === nothing
        shuffle = ifelse(rng===nothing, false, true)
    end

    if rng === nothing
        rng = Random.GLOBAL_RNG
    end

    return shuffle, rng
end

"""
    holdout = Holdout(; fraction_train=0.7,
                         shuffle=nothing,
                         rng=nothing)

Holdout resampling strategy, for use in `evaluate!`, `evaluate` and in tuning.

    train_test_pairs(holdout, rows)

Returns the pair `[(train, test)]`, where `train` and `test` are
vectors such that `rows=vcat(train, test)` and
`length(train)/length(rows)` is approximatey equal to fraction_train`.

Pre-shuffling of `rows` is controlled by `rng` and `shuffle`. If `rng`
is an integer, then the `Holdout` keyword constructor resets it to
`MersenneTwister(rng)`. Otherwise some `AbstractRNG` object is
expected.

If `rng` is left unspecified, `rng` is reset to `Random.GLOBAL_RNG`,
in which case rows are only pre-shuffled if `shuffle=true` is
specified.

"""
struct Holdout <: ResamplingStrategy
    fraction_train::Float64
    shuffle::Bool
    rng::Union{Int,AbstractRNG}

    function Holdout(fraction_train, shuffle, rng)
        0 < fraction_train < 1 ||
            error("`fraction_train` must be between 0 and 1.")
        return new(fraction_train, shuffle, rng)
    end
end

# Keyword Constructor
Holdout(; fraction_train::Float64=0.7, shuffle=nothing, rng=nothing) =
    Holdout(fraction_train, shuffle_and_rng(shuffle, rng)...)

function train_test_pairs(holdout::Holdout, rows)

    train, test = partition(rows, holdout.fraction_train,
                          shuffle=holdout.shuffle, rng=holdout.rng)
    return [(train, test),]

end

"""
    cv = CV(; nfolds=6,  shuffle=nothing, rng=nothing)

Cross-validation resampling strategy, for use in `evaluate!`,
`evaluate` and tuning.

    train_test_pairs(cv, rows)

Returns an `nfolds`-length iterator of `(train, test)` pairs of
vectors (row indices), where each `train` and `test` is a sub-vector
of `rows`. The `test` vectors are mutually exclusive and exhaust
`rows`. Each `train` vector is the complement of the corresponding
`test` vector. With no row pre-shuffling, the order of `rows` is
preserved, in the sense that `rows` coincides precisely with the
concatenation of the `test` vectors, in the order they are
generated. All but the last `test` vector have equal length.

Pre-shuffling of `rows` is controlled by `rng` and `shuffle`. If `rng`
is an integer, then the `CV` keyword constructor resets it to
`MersenneTwister(rng)`. Otherwise some `AbstractRNG` object is
expected.

If `rng` is left unspecified, `rng` is reset to `Random.GLOBAL_RNG`,
in which case rows are only pre-shuffled if `shuffle=true` is
explicitly specified.

"""
struct CV <: ResamplingStrategy
    nfolds::Int
    shuffle::Bool
    rng::Union{Int,AbstractRNG}
    function CV(nfolds, shuffle, rng)
        nfolds > 1 || error("Must have nfolds > 1. ")
        return new(nfolds, shuffle, rng)
    end
end

# Constructor with keywords
CV(; nfolds::Int=6,  shuffle=nothing, rng=nothing) =
    CV(nfolds, shuffle_and_rng(shuffle, rng)...)

function train_test_pairs(cv::CV, rows)

    n_observations = length(rows)
    nfolds = cv.nfolds

    if cv.shuffle
        rows=shuffle!(cv.rng, collect(rows))
    end

    # number of observations per fold
    k = floor(Int, n_observations/nfolds)
    k > 0 || error("Inusufficient data for $nfolds-fold cross-validation.\n"*
                   "Try reducing nfolds. ")

    # define the (trainrows, testrows) pairs:
    firsts = 1:k:((nfolds - 1)*k + 1) # itr of first `test` rows index
    seconds = k:k:(nfolds*k)          # itr of last  `test` rows index

    ret = map(1:nfolds) do k
        f = firsts[k]
        s = seconds[k]
        k < nfolds || (s = n_observations)
        return (vcat(rows[1:(f - 1)], rows[(s + 1):end]), # trainrows
                rows[f:s])                                # testrows
    end

    return ret
end

"""
    stratified_cv = StratifiedCV(; nfolds=6,
                                   shuffle=false,
                                   rng=Random.GLOBAL_RNG)

Stratified cross-validation resampling strategy, for use in
`evaluate!`, `evaluate` and in tuning. Applies only to classification
problems (`OrderedFactor` or `Multiclass` targets).

    train_test_pairs(stratified_cv, rows, y)

Returns an `nfolds`-length iterator of `(train, test)` pairs of
vectors (row indices) where each `train` and `test` is a sub-vector of
`rows`. The `test` vectors are mutually exclusive and exhaust
`rows`. Each `train` vector is the complement of the corresponding
`test` vector.

Unlike regular cross-validation, the distribution of the levels of the
target `y` corresponding to each `train` and `test` is constrained, as
far as possible, to replicate that of `y[rows]` as a whole.

Specifically, the data is split into a number of groups on which `y`
is constant, and each individual group is resampled according to the
ordinary cross-validation strategy `CV(nfolds=nfolds)`. To obtain the
final `(train, test)` pairs of row indices, the per-group pairs are
collated in such a way that each collated `train` and `test` respects
the original order of `rows` (after shuffling, if `shuffle=true`).

Pre-shuffling of `rows` is controlled by `rng` and `shuffle`. If `rng`
is an integer, then the `StratifedCV` keyword constructor resets it to
`MersenneTwister(rng)`. Otherwise some `AbstractRNG` object is
expected.

If `rng` is left unspecified, `rng` is reset to `Random.GLOBAL_RNG`,
in which case rows are only pre-shuffled if `shuffle=true` is
explicitly specified.

"""
struct StratifiedCV <: ResamplingStrategy
    nfolds::Int
    shuffle::Bool
    rng::Union{Int,AbstractRNG}
    function StratifiedCV(nfolds, shuffle, rng)
        nfolds > 1 || error("Must have nfolds > 1. ")
        return new(nfolds, shuffle, rng)
    end
end

# Constructor with keywords
StratifiedCV(; nfolds::Int=6,  shuffle=nothing, rng=nothing) =
       StratifiedCV(nfolds, shuffle_and_rng(shuffle, rng)...)

function train_test_pairs(stratified_cv::StratifiedCV, rows, X, y)

    n_observations = length(rows)
    nfolds = stratified_cv.nfolds

    if stratified_cv.shuffle
        rows=shuffle!(stratified_cv.rng, collect(rows))
    end

    st = scitype(y)
    st <: AbstractArray{<:Finite} ||
        error("Supplied target has scitpye $st but stratified "*
              "cross-validation applies only to classification problems. ")

    freq_given_level = countmap(y[rows])
    minimum(values(freq_given_level)) >= nfolds ||
        error("The number of observations for which the target takes on a "*
              "given class must, for each class, exceed `nfolds`. Try "*
              "reducing `nfolds`. ")

    levels_seen = keys(freq_given_level) |> collect

    cv = CV(nfolds=nfolds)

    # the target is constant on each stratum, a subset of `rows`:
    class_rows = [rows[y[rows] .== c] for c in levels_seen]

    # get the cv train/test pairs for each level:
    train_test_pairs_per_level = (train_test_pairs(cv, class_rows[m])
                              for m in eachindex(levels_seen))

    # just the train rows in each level:
    trains_per_level = map(x -> first.(x),
                           train_test_pairs_per_level)

    # just the test rows in each level:
    tests_per_level  = map(x -> last.(x),
                                train_test_pairs_per_level)

    # for each fold, concatenate the train rows over levels:
    trains_per_fold = map(x->vcat(x...), zip(trains_per_level...))

    # for each fold, concatenate the test rows over levels:
    tests_per_fold = map(x->vcat(x...), zip(tests_per_level...))

    # restore ordering specified by rows:
    trains_per_fold = map(trains_per_fold) do train
        filter(in(train), rows)
    end
    tests_per_fold = map(tests_per_fold) do test
        filter(in(test), rows)
    end

    # re-assemble:
    return zip(trains_per_fold, tests_per_fold) |> collect

end

## EVALUATION TYPE

const PerformanceEvaluation = NamedTuple{(:measure, :measurement,
                               :per_fold, :per_observation)}
# pretty printing:
round3(x) = round(x, sigdigits=3)
_short(v::Vector{<:Real}) = MLJBase.short_string(v)
_short(v::Vector) = string("[", join(_short.(v), ", "), "]")
_short(::Missing) = missing

function Base.show(io::IO, ::MIME"text/plain", e::PerformanceEvaluation)
    data = hcat(e.measure, round3.(e.measurement),
                [round3.(v) for v in e.per_fold])
    header = ["_.measure", "_.measurement", "_.per_fold"]
    PrettyTables.pretty_table(io, data, header;
                              header_crayon=Crayon(bold=false),
                              alignment=:l)
    println(io, "_.per_observation = $(_short(e.per_observation))")
end

function Base.show(io::IO, e::PerformanceEvaluation)
    summary = Tuple(round3.(e.measurement))
    print(io, "PerformanceEvaluation$summary")
end

## EVALUATION METHODS

function _check_measure(model, measure, y, operation, override)

    override && (return nothing)

    T = scitype(y)

    T == Unknown && (return nothing)
    target_scitype(measure) == Unknown && (return nothing)
    prediction_type(measure) == :unknown && (return nothing)

    avoid = "\nTo override measure checks, set check_measure=false. "

    T <: target_scitype(measure) ||
        throw(ArgumentError(
            "\nscitype of target = $T but target_scitype($measure) = "*
            "$(target_scitype(measure))."*avoid))

    if model isa Probabilistic
        if operation == predict
            if prediction_type(measure) != :probabilistic
                suggestion = ""
                if target_scitype(measure) <: Finite
                    suggestion = "\nPerhaps you want to set operation="*
                    "predict_mode. "
                elseif target_scitype(measure) <: Continuous
                    suggestion = "\nPerhaps you want to set operation="*
                    "predict_mean or operation=predict_median. "
                else
                    suggestion = ""
                end
                throw(ArgumentError(
                   "\n$model <: Probabilistic but prediction_type($measure) = "*
                      ":$(prediction_type(measure)). "*suggestion*avoid))
            end
        end
    end

    model isa Deterministic && prediction_type(measure) != :deterministic &&
        throw(ArgumentError("$model <: Deterministic but "*
                            "prediction_type($measure) ="*
              ":$(prediction_type(measure))."*avoid))

    return nothing

end

function _process_weights_measures(weights, measures, mach,
                                   operation, verbosity, check_measure)

    if measures === nothing
        candidate = default_measure(mach.model)
        candidate ===  nothing && error("You need to specify measure=... ")
        _measures = [candidate, ]
    elseif !(measures isa AbstractVector)
        _measures = [measures, ]
    else
        _measures = measures
    end

    y = mach.args[2]

    [ _check_measure(mach.model, m, y, operation, !check_measure)
     for m in _measures ]

    if weights != nothing
        weights isa AbstractVector{<:Real} ||
            throw(ArgumentError("`weights` must be a `Real` vector."))
        length(weights) == nrows(y) ||
            throw(DimensionMismatch("`weights` and target "*
                                    "have different lengths. "))
        _weights = weights
    elseif  length(mach.args) == 3
        verbosity < 1 ||
            @info "Passing machine sample weights to any supported measures. "
        _weights = mach.args[3]
    else
        _weights = weights
    end

    return _weights, _measures

end

"""
    evaluate!(mach,
              resampling=CV(),
              measure=nothing,
              weights=nothing,
              operation=predict,
              n = 1,
              acceleration=default_resource(),
              force=false,
              verbosity=1)

Estimate the performance of a machine `mach` wrapping a supervised
model in data, using the specified `resampling` strategy (defaulting
to 6-fold cross-validation) and `measure`, which can be a single
measure or vector.

Do `subtypes(MLJ.ResamplingStrategy)` to obtain a list of available
resampling strategies. If `resampling` is not an object of type
`MLJ.ResamplingStrategy`, then a vector of pairs (of the form
`(train_rows, test_rows)` is expected. For example, setting

    resampling = [(1:100), (101:200)),
                   (101:200), (1:100)]

gives two-fold cross-validation using the first 200 rows of data.

The resampling strategy is applied repeatedly if `repeats > 1`. For
`resampling = CV(nfolds=5)`, for example, this generates a total of
`5n` test folds for evaluation and subsequent aggregation.

If `resampling isa MLJ.ResamplingStrategy` then one may optionally
restrict the data used in evaluation by specifying `rows`.

An optional `weights` vector may be passed for measures that support
sample weights (`MLJ.supports_weights(measure) == true`), which is
ignored by those that don't.

*Important:* If `mach` already wraps sample weights `w` (as in `mach =
machine(model, X, y, w)`) then these weights, which are used for
*training*, are automatically passed to the measures for
evaluation. However, for evaluation purposes, any `weights` specified
as a keyword argument will take precedence over `w`.

User-defined measures are supported; see the manual for details.

If no measure is specified, then `default_measure(mach.model)` is
used, unless this default is `nothing` and an error is thrown.

The `acceleration` keyword argument is used to specify the compute resource (a
subtype of `ComputationalResources.AbstractResource`) that will be used to
accelerate/parallelize the resampling operation.

Although evaluate! is mutating, `mach.model` and `mach.args` are
untouched.

### Return value

A property-accessible object of type `PerformanceEvaluation` with
these properties:

- `measure`: the vector of specified measures

- `measurements`: the corresponding measurements, aggregated across the
  test folds using the aggregation method defined for each measure (do
  `aggregation(measure)` to inspect)

- `per_fold`: a vector of vectors of individual test fold evaluations
  (one vector per measure)

- `per_observation`: a vector of vectors of individual observation
  evaluations of those measures for which
  `reports_each_observation(measure)` is true, which is otherwise
  reported `missing`.

See also [`evaluate`](@ref)

"""
function evaluate!(mach::Machine{<:Supervised};
                   resampling=CV(),
                   measures=nothing, measure=measures, weights=nothing,
                   operation=predict, acceleration=default_resource(),
                   rows=nothing, repeats=1, force=false,
                   check_measure=true, verbosity=1)

    # this method just checks validity of options, preprocess the
    # weights and measures, and dispatches a strategy-specific
    # `evaluate!`

    repeats > 0 || error("Need n > 0. ")

    if resampling isa TrainTestPairs
        if rows !== nothing
            error("You cannot specify `rows` unless `resampling "*
                  "isa MLJ.ResamplingStrategy` is true. ")
        end
        if repeats != 1 && verbosity > 0
            @warn "repeats > 1 not supported unless "*
            "`resampling<:ResamplingStrategy. "
        end
    end

    _weights, _measures =
        _process_weights_measures(weights, measure, mach,
                                  operation, verbosity, check_measure)

    if verbosity >= 0 && weights !== nothing
        unsupported = filter(_measures) do m
            !supports_weights(m)
        end
        if !isempty(unsupported)
            unsupported_as_string = string(unsupported[1])
            unsupported_as_string *=
                reduce(*, [string(", ", m) for m in unsupported[2:end]])
                @warn "Sample weights ignored in evaluations of the following"*
            " measures, as unsupported: \n$unsupported_as_string "
        end
    end

    evaluate!(mach, resampling, _weights, rows, verbosity, repeats,
                   _measures, operation, acceleration, force)

end

"""
    evaluate(model, X, y; measure=nothing, options...)
    evaluate(model, X, y, w; measure=nothing, options...)

Evaluate the performance of a supervised model `model` on input data
`X` and target `y`, optionally specifying sample weights `w` for
training, where supported. The same weights are passed to measures
that support sample weights, unless this behaviour is overridden by
explicitly specifying the option `weights=...`.

See the machine version `evaluate!` for the complete list of options.

"""
evaluate(model::Supervised, args...; kwargs...) =
    evaluate!(machine(model, args...); kwargs...)

const AbstractRow = Union{AbstractVector{<:Integer}, Colon}
const TrainTestPair = Tuple{AbstractRow,AbstractRow}
const TrainTestPairs = AbstractVector{<:TrainTestPair}

function _evaluate!(func::Function, res::CPU1, nfolds, verbosity)
    p = Progress(nfolds + 1, dt=0, desc="Evaluating over $nfolds folds: ",
                 barglyphs=BarGlyphs("[=> ]"), barlen=25, color=:yellow)
    verbosity > 0 && next!(p)
    return reduce(vcat, (func(k, p, verbosity) for k in 1:nfolds))
end
function _evaluate!(func::Function, res::CPUProcesses, nfolds, verbosity)
    # TODO: use pmap here ?:
    return @distributed vcat for k in 1:nfolds
        func(k)
    end
end
@static if VERSION >= v"1.3.0-DEV.573"
    function _evaluate!(func::Function, res::CPUThreads, nfolds, verbosity)
        task_vec = [Threads.@spawn func(k) for k in 1:nfolds]
        return reduce(vcat, fetch.(task_vec))
    end
end

# Evaluation when resampling is a TrainTestPairs (core evaluator):
function evaluate!(mach::Machine, resampling, weights,
                   rows, verbosity, repeats,
                   measures, operation, acceleration, force)

    # Note: `rows` and `n` are ignored here

    resampling isa TrainTestPairs ||
        error("`resampling` must be an "*
              "`MLJ.ResamplingStrategy` or tuple of pairs "*
              "of the form `(train_rows, test_rows)`")

    X = mach.args[1]
    y = mach.args[2]

    nfolds = length(resampling)

    nmeasures = length(measures)

    function get_measurements(k)
        train, test = resampling[k]
        fit!(mach; rows=train, verbosity=verbosity-1, force=force)
        Xtest = selectrows(X, test)
        ytest = selectrows(y, test)
        if weights == nothing
            wtest = nothing
        else
            wtest = weights[test]
        end
        yhat = operation(mach, Xtest)
        return [value(m, yhat, Xtest, ytest, wtest)
                for m in measures]
    end
    function get_measurements(k, p, verbosity) # p = progress meter
        ret = get_measurements(k)
        verbosity > 0 && next!(p)
        return ret
    end

    measurements_flat = if acceleration isa CPUProcesses
        ## TODO: progress meter for distributed case
        if verbosity > 0
            @info "Distributing cross-validation computation " *
                  "among $(nworkers()) workers."
        end
    end
    measurements_flat =
        _evaluate!(get_measurements, acceleration, nfolds, verbosity)

    # in the following rows=folds, columns=measures:
    measurements_matrix = permutedims(
        reshape(measurements_flat, (nmeasures, nfolds)))

    # measurements for each observation:
    per_observation = map(1:nmeasures) do k
        m = measures[k]
        if reports_each_observation(m)
            [measurements_matrix[:,k]...]
        else
            missing
        end
    end

    # measurements for each fold:
    per_fold = map(1:nmeasures) do k
        m = measures[k]
        if reports_each_observation(m)
            broadcast(MLJBase.aggregate, per_observation[k], [m,])
        else
            [measurements_matrix[:,k]...]
        end
    end

    # overall aggregates:
    per_measure = map(1:nmeasures) do k
        m = measures[k]
        MLJBase.aggregate(per_fold[k], m)
    end

    ret = (measure=measures,
           measurement=per_measure,
           per_fold=per_fold,
           per_observation=per_observation)

    return ret

end

function actual_rows(rows, N, verbosity)
    unspecified_rows = (rows === nothing)
    _rows = unspecified_rows ? (1:N) : rows
    unspecified_rows || @info "Creating subsamples from a subset of all rows. "
    return _rows
end

# Evaluation when resampling is a ResamplingStrategy:
function evaluate!(mach::Machine, resampling::ResamplingStrategy,
                   weights, rows, verbosity, repeats, args...)

    y = mach.args[2]
    _rows = actual_rows(rows, length(y), verbosity)

    repeated_train_test_pairs =
        vcat([train_test_pairs(resampling, _rows, mach.args...)
              for i in 1:repeats]...)

    return evaluate!(mach::Machine,
                     repeated_train_test_pairs,
                     weights, nothing, verbosity, repeats, args...)

end

## RESAMPLER - A MODEL WRAPPER WITH `evaluate` OPERATION

"""
    resampler = Resampler(model=ConstantRegressor(),
                          resampling=CV(),
                          measure=nothing,
                          weights=nothing,
                          operation=predict,
                          repeats = 1,
                          acceleration=default_resource(),
                          check_measure=true)

Resampling model wrapper, used internally by the `fit` method of
`TunedModel` instances. See [`evaluate!](@ref) for options. Not
intended for general use.

Given a machine `mach = machine(resampler, args...)` one obtains a
performance evaluation of the specified `model`, performed according
to the prescribed `resampling` strategy and other parameters, using
data `args...`, by calling `fit!(mach)` followed by
`evaluate(mach)`. The advantage over using `evaluate(model, X, y)` is
that the latter call always calls `fit` on the `model` but
`fit!(mach)` only calls `update` after the first call.

The sample `weights` are passed to the specified performance
measures that support weights for evaluation.

*Important:* If `weights` are left unspecified, then any weight vector
`w` used in constructing the resampler machine, as in
`resampler_machine = machine(resampler, X, y, w)` (which is then used
in *training* the model) will also be used in evaluation.

"""
mutable struct Resampler{S,M<:Supervised} <: Supervised
    model::M
    resampling::S # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    operation
    acceleration::AbstractResource
    check_measure::Bool
    repeats::Int
end

MLJBase.is_wrapper(::Type{<:Resampler}) = true
MLJBase.supports_weights(::Type{<:Resampler{<:Any,M}}) where M =
    supports_weights(M)
MLJBase.is_pure_julia(::Type{<:Resampler}) = true

function MLJBase.clean!(resampler::Resampler)
    warning = ""
    if resampler.measure === nothing
        measure = default_measure(resampler.model)
        if measure === nothing
            error("No default measure known for $(resampler.model). "*
                  "You must specify measure=... ")
        else
            warning *= "No `measure` specified. "*
            "Setting `measure=$measure`. "
        end
    end
    return warning
end

function Resampler(; model=ConstantRegressor(), resampling=CV(),
            measure=nothing, weights=nothing, operation=predict,
            acceleration=default_resource(), check_measure=true, repeats=1)

    resampler = Resampler(model, resampling, measure, weights, operation,
                          acceleration, check_measure, repeats)
    message = MLJBase.clean!(resampler)
    isempty(message) || @warn message

    return resampler

end

function MLJBase.fit(resampler::Resampler, verbosity::Int, args...)

    mach = machine(resampler.model, args...)

    weights, measures =
        _process_weights_measures(resampler.weights, resampler.measure,
                                  mach, resampler.operation,
                                  verbosity, resampler.check_measure)

    fitresult = evaluate!(mach, resampler.resampling,
                          weights, nothing, verbosity - 1, resampler.repeats,
                          measures, resampler.operation,
                          resampler.acceleration, false)

    cache = (mach, deepcopy(resampler.resampling))
    report = NamedTuple()

    return fitresult, cache, report

end

# in special case of holdout, we can reuse the underlying model's
# machine, provided the training_fraction has not changed:
function MLJBase.update(resampler::Resampler{Holdout},
                        verbosity::Int, fitresult, cache, args...)

    old_mach, old_resampling = cache

    if old_resampling.fraction_train == resampler.resampling.fraction_train
        mach = old_mach
    else
        mach = machine(resampler.model, args...)
        cache = (mach, deepcopy(resampler.resampling))
    end

    weights, measures =
        _process_weights_measures(resampler.weights, resampler.measure,
                                  mach, resampler.operation,
                                  verbosity, resampler.check_measure)
    mach.model = resampler.model
    fitresult = evaluate!(mach, resampler.resampling,
                          weights, nothing, verbosity - 1, resampler.repeats,
                          measures, resampler.operation,
                          resampler.acceleration, false)

    report = NamedTuple

    return fitresult, cache, report

end

MLJBase.input_scitype(::Type{<:Resampler{S,M}}) where {S,M} =
    MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:Resampler{S,M}}) where {S,M} =
    MLJBase.target_scitype(M)
MLJBase.package_name(::Type{<:Resampler}) = "MLJBase"

MLJBase.load_path(::Type{<:Resampler}) = "MLJBase.Resampler"

evaluate(resampler::Resampler, fitresult) = fitresult

function evaluate(machine::AbstractMachine{<:Resampler})
    if isdefined(machine, :fitresult)
        return evaluate(machine.model, machine.fitresult)
    else
        throw(error("$machine has not been trained."))
    end
end


text:
import SparseGrids:NGrid,Quadratic,Linear
function initS(slist::Vector{Variable}, B = Quadratic, G1 = nothing)
    if G1 == nothing
        G = NGrid(Int[SparseGrids.cc_iM(length(x)) - 1 for x in slist], hcat(Vector{Float64}[x.bounds for x in slist]...), B = B)
    else
        G = G1
    end
    S = values(G)
    for s in slist
        s.x = sort(unique(S[:,find(slist, s.name)]))
    end

    exogenousfuture = map(x -> repmat(x[:]', length(G), 1), ndgrid([s.x for s in Stochastic(slist)]...))
    nP = prod(map(length, Stochastic(slist)))
    ProbWeights =  Float64[prod([Stochastic(slist)[N].T[findfirst(Stochastic(slist)[N].x, S[i,findfirst(names(State(slist)), names(Stochastic(slist))[N])]), findfirst(Stochastic(slist)[N].x, exogenousfuture[N][i,j])] for N = 1:length(Stochastic(slist))]) for i = 1:length(G), j = 1:nP]
    SP = zeros(length(G) * nP, size(S, 2))
    for i = 1:length(Stochastic(slist))
        SP[:,findfirst(names(slist), Stochastic(slist)[i].name)] = exogenousfuture[i][:]
    end
    return G, S, SP, ProbWeights
end

function initX(M::Model)
    for i = 1:length(M.variables)
        v = M.variables[i]
        if isa(v, Policy) || isa(v, Exogenous) || isa(v, Aggregate)
            if isa(v.init, Number)
                M.X[:,i] = v.init
            elseif isa(v.init, Expr)
                init = addindex!(v.init)
                initv = getv(init)
                @assert all(map(x -> in(x.args[1], names(State(M.variables))), initv)) "Initialisation of $(v.name) includes a non state variable"
                for j = 1:length(M)
                    initD = Dict(x => M.X[j,find(State(M.variables), x.args[1])] for x in initv)
                    M.X[j,i] = eval(current_module(), subs(init, initD))
                end
            end
        end
    end
end

function initD(M::Model)
    hloc = hardloc(M.variables)
    for v in Dependant(M.variables)
        F = quote
                function $(gensym(v.name))(M)
                    for i = 1:length(M)
                        $(M.X)[i,$(find(M.variables, v.name))] = $(subs(v.e, hloc))
                    end
                end
            end
        v.update = eval(current_module(), F)
    end
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
"""
    sparsity(A)

Calculates the sparsity of the given array.
The sparsity is defined as number of zero-valued elements divided by
total number of elements.
"""
function sparsity(A::AbstractArray{T}) where T<:Number
    (length(A)-countnz(A))/length(A)
end

"""
    reldiff(A, B)

Relative difference of absolute values of `A` and `B` defined as

``
\\operatorname{reldiff} = 2 \\dfrac{\\operatorname{abs}(A - B)}{\\operatorname{abs}(A+B)}.
``
"""
function reldiff(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
    return 2*abs.(A-B)./abs.(A+B)
end

"""
    effreldiff(A, B, threshold=1e-14)

Same as `reldiff(A,B)` but with all elements set to zero where corresponding element of
`absdiff(A,B)` is smaller than `threshold`. This is useful in avoiding artificially large
relative errors.
"""
function effreldiff(A::AbstractArray{T}, B::AbstractArray{S}, threshold::Float64=1e-14) where T<:Number where S<:Number
    r = reldiff(A,B)
    r[findall(x -> abs.(x)<threshold, absdiff(A,B))] .= 0.
    return r
end

"""
    absdiff(A, B)

Difference of absolute values of `A` and `B`.
"""
function absdiff(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
    return abs.(A-B)
end

"""
    compare(A, B)

Compares two matrices `A` and `B`, prints out the maximal absolute and relative differences
and returns a boolean indicating wether `isapprox(A,B)`.
"""
function compare(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
    @printf("max absdiff: %.1e\n", maximum(absdiff(A,B)))
    @printf("mean absdiff: %.1e\n", mean(absdiff(A,B)))
    @printf("max reldiff: %.1e\n", maximum(reldiff(A,B)))
    @printf("mean reldiff: %.1e\n", mean(reldiff(A,B)))

    r = effreldiff(A,B)
    @printf("effective max reldiff: %.1e\n", maximum(r))
    @printf("effective mean reldiff: %.1e\n", mean(r))

    return isapprox(A,B)
end

# NOTE
# Currenlty julia/sparearrays does not implement this function (type signature)
# once it does this can be removed/depracted in favor of mul!
# see also: test/slice_matrices.jl
occursin(
    "SparseArrays",
    string(which(mul!, (Matrix, Matrix, SparseMatrixCSC)).file)
) && @warn(
    "A Method `mul!(::Matrix, ::Matrix, ::SparseMatrixCSC)` now exists in " *
    "`SparseArrays`. The method defined in `helpers.jl` is likely to be  " *
    "unnecessary now."
)

function SparseArrays.mul!(C::StridedMatrix, X::StridedMatrix, A::SparseMatrixCSC)
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch())
    fill!(C, zero(eltype(C)))
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for  col = 1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        ki=rowval[k]
        kv=nzval[k]
        for multivec_row=1:mX
            C[multivec_row, col] += X[multivec_row, ki] * kv
        end
    end
    C
end


text:
@doc "All the abstract domains" ->
all_domains() = Set(Interval,AbstractBool,HyperBox)

@doc "Return a type Union of all the abstract domains which can abstract a type T" ->
abstractdomains(CT::DataType) = Union{filter(AT->isabstract(CT,AT),all_domains())...}

## Random Number Generation
## =========================
@doc "Random number between `a` and `b`" ->
rand_interval(a::Float64, b::Float64) = a + (b - a) * rand()


text:
struct PermuteFactor{S<:Unsigned, T<:Unsigned} <: AbstractFactor{T}
    name::String
    levels::AbstractVector{<:AbstractString}
    basefactor::AbstractFactor{S}
    newindex::Vector{T}
end

Base.length(factor::PermuteFactor{T}) where {T<:Unsigned} = length(factor.basefactor)

# function PermuteFactor(name::String, basefactor::AbstractFactor{T}, islessfun::Function) where {T<:Unsigned}
#     baselevels = getlevels(basefactor)
#     levelcount = length(baselevels)
#     perm = sortperm(baselevels, lt = islessfun)
#     sortlevels = baselevels[perm]
#     newindex = Vector{T}(levelcount + 1)
#     newindex[1] = 0
#     for i in 1:levelcount
#         newindex[i + 1] = perm[i]
#     end
#     PermuteFactor{T, T}(name, sortlevels, basefactor, newindex)
# end

# function PermuteFactor(name::String, basefactor::AbstractFactor{T}) where {T<:Unsigned}
#     PermuteFactor{T}(name, basefactor, isless)
# end

function slice(factor::PermuteFactor{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned} where {S<:Unsigned}
    newindex = factor.newindex
    f = (i::S -> newindex[i + 1])
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(factor.basefactor, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, T)
end

function PermuteFactor(basefactor::AbstractFactor{T}, factortopermute::AbstractFactor{S}) where {T<:Unsigned} where {S<:Unsigned}
    baselevels = getlevels(basefactor)
    baselen = length(baselevels)
    permlevels = getlevels(factortopermute)
    unionlevels = union(Set(baselevels), Set(permlevels))
    difflevels = [x for x in setdiff(unionlevels, Set(baselevels))]
    levels = [i <= baselen ? baselevels[i] : difflevels[i - baselen] for i in 1:length(unionlevels)]
    newindex = [i == 0 ? 0 : findfirst(levels, permlevels[i]) for i in 0:length(permlevels)]
    if length(levels) <= typemax(UInt8)
        PermuteFactor{S, UInt8}(getname(factortopermute), levels, factortopermute, newindex)
    elseif length(levels) <= typemax(UInt16)
        PermuteFactor{S, UInt16}(getname(factortopermute), levels, factortopermute, newindex)
    else
        PermuteFactor{S, UInt32}(getname(factortopermute), levels, factortopermute, newindex)
    end
end

function Base.map(factor::PermuteFactor, dataframe::AbstractDataFrame)
    islessfun = (x, y) -> isless(findfirst(factor.levels, x), findfirst(factor.levels, y))
    OrdinalFactor(factor.name, map(factor.basefactor, dataframe), islessfun)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
include("run_base.jl")
include("run_utils.jl")

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA