text:
using RandomNumbers.Xorshifts
using POMDPModelTools

@testset "PMDPe.jl" begin

    # initialize
    pp = simple_pp()
    rng = Xorshift128Plus(1)
    me = PMDPs.PMDPe(pp)

    @test typeof(me) == PMDPs.PMDPe

    # Test methods
    @testset "states and indices" begin
        tiny_pp = PMDPs.PMDPProblem(
            [PMDPs.Product([true], 3),],
            SA[1],
            BernoulliScheme(3, [0.4]),
            [DiscreteNonParametric([10.], [1.]),],
            [0., 5., 10., 15.],
            :revenue
        )

        tiny_me = PMDPs.PMDPe(tiny_pp)

        states = [
            PMDPs.State(SA[1], 1, 2),
            PMDPs.State(SA[1], 2, 2),
            PMDPs.State(SA[1], 3, 2),
            PMDPs.State(SA[1], 1, 1),
            PMDPs.State(SA[1], 2, 1),
            PMDPs.State(SA[1], 3, 1),
            PMDPs.State(SA[0], 1, 2),
            PMDPs.State(SA[0], 2, 2),
            PMDPs.State(SA[0], 3, 2),
            PMDPs.State(SA[0], 1, 1),
            PMDPs.State(SA[0], 2, 1),
            PMDPs.State(SA[0], 3, 1)
            ]

        @test Set(PMDPs.generate_states(tiny_pp)) == Set(states)

        state_indices = PMDPs.stateindices(tiny_pp)
        @test Set(map(s->POMDPs.stateindex(tiny_me, s), states)) == Set(1:12)
    end

    # transition
    @testset "POMDPs.transition" begin
        s₀ = PMDPs.State(pp.c₀, 1, 1)
        a = 10.
        td = POMDPs.transition(me, s₀, a)
        @test isa(td, SparseCat)
        @test length(td.vals) == 8 # requested product is either sold or not and pp has 4 different products (including empty) -> 8 outcomes
        @test sum(td.probs) == 1.
    end
end


text:
"""
   ApproximateBayesianComputing
Module providing
   types (abc_pmc_plan_type, abc_population_type) and methods (generate_theta, init_abc, update_abc_pop, run_abc) for using Approximate Bayesian Computing
"""
module ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing

# package code goes here

#using Distributions
using PDMats
using DistributedArrays
using Compat
if VERSION >= v"0.7"
  using Statistics
  using Distributed
  import Base: length
  import Statistics: mean, median, maximum, minimum, quantile, std, var, cov, cor
  import LinearAlgebra.isposdef
else
  #import Base: sum, mean, std
  import Compat: sum, mean, std
  import Base: length, median, maximum, minimum, quantile, var, cov, cor
  using Compat.Statistics
  using Compat.Distributed
end

export
  # types
  abc_plan_type,
  abc_pmc_plan_type,
  abc_population_type,

  # util methods

  # abc alg methods
  generate_theta,
  init_abc,
  update_abc_pop,
  run_abc

import StatsBase.sample
import Base: rand
import Compat: @compat

# until added to distributions, use our own
include("GaussianMixtureModelCommonCovar.jl")

include("types.jl")
include("util.jl")
include("alg.jl")
include("alg_serial.jl")
include("alg_parallel.jl")
#include("alg_parallel_custom.jl")
include("make_proposal.jl")
include("log.jl")
#include("emulator.jl")

# Bonus code!
include("composite.jl")
export CompositeDistributions
include("beta_linear_transformed.jl")
export TransformedBetaDistributions

end # module

ext: jl


text:
import Base.string
import Base.print
import Base.show
import Base.convert
import Base.length
import Base.==
import Base.hash

import Dates, Intervals

export DbId, SQLType, AbstractModel
export SQLInput, SQLColumn, SQLColumns, SQLLogicOperator
export SQLWhere, SQLWhereExpression, SQLWhereEntity, SQLLimit, SQLOrder, SQLQuery, SQLRaw
export SQLJoin, SQLOn, SQLJoinType, SQLHaving

export @sql_str

abstract type SearchLightAbstractType end
abstract type SQLType <: SearchLightAbstractType end

abstract type AbstractModel <: SearchLightAbstractType end

function hash(a::T) where {T<:AbstractModel}
  Base.hash(string(typeof(a)) * string(getfield(a, pk(a) |> Symbol)))
end

function ==(a::A, b::B) where {A<:AbstractModel,B<:AbstractModel}
  hash(a) == hash(b)
end

function Base.print(io::IO, t::T) where {T<:SearchLightAbstractType}
  props = []
  for (k,v) in to_string_dict(t)
    push!(props, "$k=$v")
  end
  print(io, string("$(typeof(t))(", join(props, ','), ")"))
end

Base.show(io::IO, t::T) where {T<:SearchLightAbstractType} = print(io, searchlightabstracttype_to_print(t))

"""
    searchlightabstracttype_to_print{T<:SearchLightAbstractType}(m::T) :: String

Pretty printing of SearchLight types.
"""
function searchlightabstracttype_to_print(m::T) :: String where {T<:SearchLightAbstractType}
  string(typeof(m), "\n", Millboard.table(to_string_dict(m)), "\n")
end

mutable struct DbId
  value::Union{Nothing,Int,String}
end
DbId() = DbId(nothing)
DbId(id::Number) = DbId(Int(id))
DbId(id::AbstractString) = DbId(id)

function hash(a::DbId)
  Base.hash(a.value)
end

function ==(a::DbId, b::DbId)
  hash(a) == hash(b)
end

Base.convert(::Type{DbId}, v::Union{Number,String}) = DbId(v)
Base.convert(::Type{DbId}, v::Nothing) = DbId(nothing)

Base.convert(::Type{String}, id::DbId) = string(id.value)
Base.convert(::Type{Union{Int,String}}, v::DbId) = v.value

function Base.convert(::Type{Int}, v::DbId)
  isa(v.value, Int) && return v.value

  x = if isa(v.value, Number) || isa(v.value, String)
    tryparse(Int, v.value)
  else
    nothing
  end

  x === nothing && return 0
end

Base.show(io::IO, dbid::DbId) = print(io, (dbid.value === nothing ? "NULL" : string(dbid.value)))

#
# SQLRaw
#

"""
Wrapper around a raw SQL query part.
"""
struct SQLRaw <: SQLType
  value::String
end

macro sql_str(q)
  SQLRaw(q)
end

string(io::IO, t::SQLRaw) = t.value

#
# SQLInput
#

"""
Provides safe input into SQL queries and operations related to that.
"""
mutable struct SQLInput <: SQLType
  value::Union{AbstractString,Real}
  escaped::Bool
  raw::Bool
  SQLInput(v::Union{AbstractString,Real}; escaped = false, raw = false) = new(v, escaped, raw)
end
SQLInput(a::Dates.Date) = string(a) |> SQLInput
SQLInput(a::Vector{T}) where {T} = SQLInput.(a)
SQLInput(s::SubString{T}) where {T} = convert(String, s) |> SQLInput
SQLInput(i::SQLInput) = i
SQLInput(s::Symbol) = string(s) |> SQLInput
SQLInput(r::SQLRaw) = SQLInput(r.value, raw = true)
SQLInput(n::Nothing) = SQLInput("NULL", escaped = true, raw = true)
SQLInput(i::Intervals.Interval) = SQLInput(replace(string(i),".."=>","))
SQLInput(a::Any) = string(a) |> SQLInput

==(a::SQLInput, b::SQLInput) = a.value == b.value

string(s::SQLInput) = "$(safe(s).value)"
string(a::Vector{SQLInput}) = join(map(x -> string(x), a), ",")
endof(s::SQLInput) = endof(s.value)
length(s::SQLInput) = length(s.value)
next(s::SQLInput, i::Int) = next(s.value, i)
safe(s::SQLInput) = escape_value(s)

print(io::IO, s::SQLInput) = print(io, string(s))
show(io::IO, s::SQLInput) = print(io, string(s))

convert(::Type{SQLInput}, r::Real) = SQLInput(parse(r))
convert(::Type{SQLInput}, s::Symbol) = SQLInput(string(s))
convert(::Type{SQLInput}, d::Dates.DateTime) = SQLInput(string(d))
convert(::Type{SQLInput}, d::Dates.Date) = SQLInput(string(d))
convert(::Type{SQLInput}, d::Dates.Time) = SQLInput(string(d))
convert(::Type{SQLInput}, id::DbId) = SQLInput(id.value)

"""
    escape_value(i::SQLInput)

Sanitizes input to be used as values in SQL queries.
"""
function escape_value(i::SQLInput) :: SQLInput
  (i.value == "NULL" || i.value == "NOT NULL") && return i

  if ! i.escaped && ! i.raw
    i.value = SearchLight.escape_value(i.value)
    i.escaped = true
  end

  return i
end

#
# SQLColumn
#

"""
Represents a SQL column when building SQL queries.
"""
mutable struct SQLColumn <: SQLType
  value::String
  escaped::Bool
  raw::Bool
  table_name::String
  column_name::String
end
SQLColumn(v::Union{String,Symbol}; escaped = false, raw = false, table_name = "", column_name = "") = begin
  v = string(v)

  if ! raw
    v == "*" && (raw = true)
    is_fully_qualified(v) && ((table_name, v) = from_fully_qualified(v))
  end

  SQLColumn(string(v), escaped, raw, string(table_name), string(v))
end
SQLColumn(a::Array) = map(x -> SQLColumn(x), a)
SQLColumn(t::Tuple) = SQLColumn([t...])
SQLColumn(c::SQLColumn) = c
SQLColumn(r::SQLRaw) = SQLColumn(r.value, raw = true)
SQLColumn(a::Any) = SQLColumn(string(a))

==(a::SQLColumn, b::SQLColumn) = a.value == b.value

string(a::Vector{SQLColumn}) = join(map(x -> string(x), a), ", ")
string(s::SQLColumn) = safe(s).value
safe(s::SQLColumn) = escape_column_name(s)

convert(::Type{SQLColumn}, s::Symbol) = SQLColumn(string(s))
convert(::Type{SQLColumn}, s::String) = SQLColumn(s)
convert(::Type{SQLColumn}, v::String, e::Bool, r::Bool) = SQLColumn(v, escaped = e, raw = r)
convert(::Type{SQLColumn}, v::String, e::Bool, r::Bool, t::Any) = SQLColumn(v, escaped = e, raw = r, table_name = string(t))

print(io::IO, a::Vector{SQLColumn}) = print(io, string(a))
show(io::IO, a::Vector{SQLColumn}) = print(io, string(a))
print(io::IO, s::SQLColumn) = print(io, string(s))
show(io::IO, s::SQLColumn) = print(io, string(s))

function SQLColumns(m::Type{T}, substitutes = NamedTuple())::Vector{SQLColumn} where {T<:AbstractModel}
  isempty(substitutes) && return SQLColumn(fieldnames(m))

  cols = Any[fieldnames(m)...]
  setdiff!(cols, substitutes |> keys |> collect)

  push!(cols, (substitutes |> values |> collect)...) |> SQLColumn
end

"""
    escape_column_name(c::SQLColumn) :: SQLColumn
    escape_column_name(s::String)

Sanitizes input to be use as column names in SQL queries.
"""
function escape_column_name(c::SQLColumn) :: SQLColumn
  if ! c.escaped && ! c.raw
    c.column_name = isempty(c.column_name) && is_fully_qualified(c.value) ? from_fully_qualified(c.value)[2] : c.value
    val = ! isempty(c.table_name) && ! startswith(c.value, (c.table_name * ".")) && ! is_fully_qualified(c.value) ? c.table_name * "." * c.value : c.value
    c.value = escape_column_name(val)
    c.escaped = true

  end

  c
end

#
# SQLLogicOperator
#

"""
Represents the logic operators (OR, AND) as part of SQL queries.
"""
struct SQLLogicOperator <: SQLType
  value::String
  SQLLogicOperator(v::String) = new( v == "OR" ? "OR" : "AND" )
end
SQLLogicOperator(v::Any) = SQLLogicOperator(string(v))
SQLLogicOperator() = SQLLogicOperator("AND")

string(s::SQLLogicOperator) = s.value

#
# SQLWhere
#

"""
Provides functionality for building and manipulating SQL `WHERE` conditions.
"""
struct SQLWhere <: SQLType
  column::SQLColumn
  value::SQLInput
  condition::SQLLogicOperator
  operator::String

  SQLWhere(column::SQLColumn, value::SQLInput, condition::SQLLogicOperator, operator::String) =
    new(column, value, condition, operator)
end
SQLWhere(column::Any, value::Any, condition::Any, operator::String) = SQLWhere(SQLColumn(column), SQLInput(value), SQLLogicOperator(condition), operator)
SQLWhere(column::SQLColumn, value::SQLInput, operator::String) = SQLWhere(column, value, SQLLogicOperator("AND"), operator)
SQLWhere(column::SQLColumn, value::SQLInput, condition::SQLLogicOperator) = SQLWhere(column, value, condition, "=")
SQLWhere(column::Any, value::Any, operator::Any) = SQLWhere(SQLColumn(column), SQLInput(value), SQLLogicOperator("AND"), operator)
SQLWhere(column::SQLColumn, value::SQLInput) = SQLWhere(column, value, SQLLogicOperator("AND"))
SQLWhere(column::Any, value::Any) = SQLWhere(SQLColumn(column), SQLInput(value))

string(w::SQLWhere) = "$(w.condition.value) ($(w.column) $(w.operator) $(enclosure(w.value, w.operator)))"
function string(w::SQLWhere, m::Type{T}) where {T <: AbstractModel}
  w.column = SQLColumn(w.column.value, escaped = w.column.escaped, raw = w.column.raw, table_name = table(m))
  "$(w.condition.value) ($(w.column) $(w.operator) $(enclosure(w.value, w.operator)))"
end
print(io::IO, w::T) where {T<:SQLWhere} = print(io, searchlightabstracttype_to_print(w))
show(io::IO, w::T) where {T<:SQLWhere} = print(io, searchlightabstracttype_to_print(w))

convert(::Type{Vector{SQLWhere}}, w::SQLWhere) = [w]

#
# SQLWhereExpression
#

"""
    SQLWhereExpression(sql_expression::String, values::T)
    SQLWhereExpression(sql_expression::String[, values::Vector{T}])

Constructs an instance of SQLWhereExpression, replacing the `?` placeholders inside `sql_expression` with
properly quoted values.

# Examples:
```julia
julia> SQLWhereExpression("slug LIKE ?", "%julia%")

SearchLight.SQLWhereExpression
+================+=============+
|            key |       value |
+================+=============+
|      condition |         AND |
+----------------+-------------+
| sql_expression | slug LIKE ? |
+----------------+-------------+
|         values |   '%julia%' |
+----------------+-------------+

julia> SQLWhereExpression("id BETWEEN ? AND ?", [10, 20])

SearchLight.SQLWhereExpression
+================+====================+
|            key |              value |
+================+====================+
|      condition |                AND |
+----------------+--------------------+
| sql_expression | id BETWEEN ? AND ? |
+----------------+--------------------+
|         values |              10,20 |
+----------------+--------------------+

julia> SQLWhereExpression("question LIKE 'what is the question\\?'")

SearchLight.SQLWhereExpression
+================+========================================+
|            key |                                  value |
+================+========================================+
|      condition |                                    AND |
+----------------+----------------------------------------+
| sql_expression | question LIKE 'what is the question?'  |
+----------------+----------------------------------------+
|         values |                                        |
+----------------+----------------------------------------+
```
"""
struct SQLWhereExpression <: SQLType
  sql_expression::String
  values::Vector{SQLInput}
  condition::String
end

const QUESTION_MARK_REPLACEMENT = "&QM;"

function SQLWhereExpression(sql_expression::String, values::Vector{SQLInput})
  condition = "AND"

  parts = split(sql_expression, " ")
  if in(parts[1], ["AND", "OR", "and", "or"])
    condition = parts |> first |> uppercase
    sql_expression = parts[2:end] |> strip
  end

  SQLWhereExpression(sql_expression, [SQLInput(isa(v.value, AbstractString) ? replace(v.value, "?" => QUESTION_MARK_REPLACEMENT) : v.value) for v in values], condition)
end

SQLWhereExpression(sql_expression::String, values...) = SQLWhereExpression(sql_expression, [values...])
SQLWhereExpression(sql_expression::String) = SQLWhereExpression(sql_expression, SQLInput[])
SQLWhereExpression(sql_expression::String, values::Dates.Date) = SQLWhereExpression(sql_expression, [SQLInput(string(values))])
SQLWhereExpression(sql_expression::String, values::Vector{T}) where {T} = SQLWhereExpression(sql_expression, SQLInput(values))
SQLWhereExpression(sql_expression::String, values::T) where {T} = SQLWhereExpression(sql_expression, [SQLInput(values)])

function string(we::SQLWhereExpression)
  string_value = we.sql_expression

  for pl in (m.match for m = eachmatch(r":[a-zA-Z0-9_-]*", string_value))
    string_value = replace(string_value, pl => SQLColumn(string(pl[2:end])))
  end

  counter = 1
  result = string_value
  pos = 0
  while ( pos = something(findnext(isequal('?'), string_value, pos+1), 0) ) > 0
    result = replace(result, '?'=>string(we.values[counter]), count = 1)

    counter += 1
  end

  result = replace(result, QUESTION_MARK_REPLACEMENT => '?')

  string(we.condition, " ", result)
end

const SQLWhereEntity = Union{SQLWhere,SQLWhereExpression}
const SQLHaving = Union{SQLWhere,SQLWhereExpression}

convert(::Type{Vector{SQLWhereEntity}}, s::String) = SQLWhereEntity[SQLWhereExpression(s)]
convert(::Type{SQLWhereEntity}, s::String) = SQLWhereExpression(s);

#
# SQLLimit
#

const SQLLimit_ALL = "ALL"
export SQLLimit_ALL

struct UnparsableSQLLimitException <: Exception
  value
end
Base.showerror(io::IO, e::UnparsableSQLLimitException) = print(io, "Can't parse SQLLimit value $(e.value)")

"""
Wrapper around SQL `limit` operator.
"""
struct SQLLimit <: SQLType
  value::Union{Int,String}
  SQLLimit(v::Int) = new(v)

  function SQLLimit(v::String)
    v = strip(uppercase(v))
    if v == SQLLimit_ALL
      return new(SQLLimit_ALL)
    else
      i = tryparse(Int, v)
      i === nothing ? throw(UnparsableSQLLimitException(v)) : i
    end
  end
end
SQLLimit() = SQLLimit(SQLLimit_ALL)

string(l::SQLLimit) = string(l.value)

convert(::Type{SQLLimit}, v::Int) = SQLLimit(v)

#
# SQLOrder
#

"""
Wrapper around SQL `order` operator.
"""
struct SQLOrder <: SQLType
  column::SQLColumn
  direction::String
  SQLOrder(column::SQLColumn, direction::String) =
    new(column, uppercase(string(direction)) == "DESC" ? "DESC" : "ASC")
end
SQLOrder(column::Union{String,Symbol}, direction::Any; raw::Bool = false) = SQLOrder(SQLColumn(column, raw = raw), string(direction))
function SQLOrder(s::Union{String,Symbol}; raw::Bool = false)
  s = String(s)

  if endswith(rstrip(uppercase(s)), " ASC") || endswith(rstrip(uppercase(s)), " DESC")
    parts = split(s, " ")
    SQLOrder(String(parts[1]), String(parts[2]), raw = raw)
  else
    SQLOrder(s, "ASC", raw = raw)
  end
end
SQLOrder(r::SQLRaw, direction::Any = "ASC") = SQLOrder(r.value, direction, raw = true)

string(o::SQLOrder) = "($(o.column) $(o.direction))"

convert(::Type{SQLOrder}, s::String) = SQLOrder(s)
convert(::Type{Vector{SQLOrder}}, o::SQLOrder) = [o]
convert(::Type{Vector{SQLOrder}}, s::Symbol) = [SQLOrder(s)]
convert(::Type{Vector{SQLOrder}}, s::String) = [SQLOrder(s)]
convert(::Type{Vector{SQLOrder}}, t::Tuple{Symbol,Symbol}) = [SQLOrder(t[1], t[2])]

#
# SQLJoin
#

#
# SQLJoin - SQLOn
#

"""
Represents the `ON` operator used in SQL `JOIN`
"""
struct SQLOn <: SQLType
  column_1::SQLColumn
  column_2::SQLColumn
  conditions::Vector{SQLWhereEntity}

  SQLOn(column_1, column_2; conditions = SQLWhereEntity[]) = new(column_1, column_2, conditions)
end
function string(o::SQLOn)
  on = " ON $(o.column_1) = $(o.column_2) "
  if ! isempty(o.conditions)
    on *= " AND " * join( map(x -> string(x), o.conditions), " AND " )
  end

  on
end

convert(::Type{Vector{SQLOn}}, j::SQLOn) = [j]

#
# SQLJoin - SQLJoinType
#

struct InvalidJoinTypeException <: Exception
  jointype::String
end

Base.showerror(io::IO, e::InvalidJoinTypeException) = print(io, "Invalid join type $(e.jointype)")

"""
Wrapper around the various types of SQL `join` (`left`, `right`, `inner`, etc).
"""
struct SQLJoinType <: SQLType
  join_type::String

  function SQLJoinType(t::Union{String,Symbol})
    t = string(t)
    accepted_values = ["inner", "INNER", "left", "LEFT", "right", "RIGHT", "full", "FULL"]
    if in(t, accepted_values)
      new(uppercase(t))
    else
      @error  """Accepted JOIN types are $(join(accepted_values, ", "))"""
      throw(InvalidJoinTypeException(t))
    end
  end
end

convert(::Type{SQLJoinType}, s::Union{String,Symbol}) = SQLJoinType(s)

string(jt::SQLJoinType) = jt.join_type

#
# SQLJoin
#

"""
Builds and manipulates SQL `join` expressions.
"""
struct SQLJoin{T<:AbstractModel} <: SQLType
  model_name::Type{T}
  on::Vector{SQLOn}
  join_type::SQLJoinType
  outer::Bool
  where::Vector{SQLWhereEntity}
  natural::Bool
  columns::Vector{SQLColumn}
end

SQLJoin(model_name::Type{T},
        on::Vector{SQLOn};
        join_type = SQLJoinType("INNER"),
        outer = false,
        where = SQLWhereEntity[],
        natural = false,
        columns = SQLColumn[]) where {T<:AbstractModel} = SQLJoin{T}(model_name, on, join_type, outer, where, natural, columns)

SQLJoin(model_name::Type{T},
        on_column_1::Union{String,SQLColumn},
        on_column_2::Union{String,SQLColumn};
        join_type = SQLJoinType("INNER"),
        outer = false,
        where = SQLWhereEntity[],
        natural = false,
        columns = SQLColumn[]) where {T<:AbstractModel} = SQLJoin(model_name, SQLOn(on_column_1, on_column_2), join_type = join_type, outer = outer, where = where, natural = natural, columns = columns)

function string(j::SQLJoin)
  sql = """ $(j.natural ? "NATURAL " : "") $(string(j.join_type)) $(j.outer ? "OUTER " : "") JOIN $( escape_column_name(table(j.model_name), SearchLight.connection())) $(join(string.(j.on), " AND ")) """
  sql *=  if ! isempty(j.where)
          SearchLight.to_where_part(j.where)
        else
          ""
        end

  sql = replace(sql, "  " => " ")

  replace(sql, " AND ON " => " AND ")
end

convert(::Type{Vector{SQLJoin}}, j::SQLJoin) = [j]

#
# SQLQuery
#

"""
    SQLQuery( columns = SQLColumn[],
              where   = SQLWhereEntity[],
              limit   = SQLLimit("ALL"),
              offset  = 0,
              order   = SQLOrder[],
              group   = SQLColumn[],
              having  = SQLWhereEntity[])

Returns a new instance of SQLQuery.

# Examples
```julia
julia> SQLQuery(where = [SQLWhereExpression("id BETWEEN ? AND ?", [10, 20])], offset = 5, limit = 5, order = :title)

SearchLight.SQLQuery
+=========+==============================================================+
|     key |                                                        value |
+=========+==============================================================+
| columns |                                                              |
+---------+--------------------------------------------------------------+
|   group |                                                              |
+---------+--------------------------------------------------------------+
|  having | Union{SearchLight.SQLWhere,SearchLight.SQLWhereExpression}[] |
+---------+--------------------------------------------------------------+
|   limit |                                                            5 |
+---------+--------------------------------------------------------------+
|  offset |                                                            5 |
+---------+--------------------------------------------------------------+
|         |                                        SearchLight.SQLOrder[ |
|         |                                         SearchLight.SQLOrder |
|         |                                      +===========+=========+ |
|         |                                      |       key |   value | |
|   order |                                                 +========... |
+---------+--------------------------------------------------------------+
|         |  Union{SearchLight.SQLWhere,SearchLight.SQLWhereExpression}[ |
|         |                               SearchLight.SQLWhereExpression |
|   where |                                                 +========... |
+---------+--------------------------------------------------------------+
```
"""
mutable struct SQLQuery <: SQLType
  columns::Vector{SQLColumn}
  where::Vector{SQLWhereEntity}
  limit::SQLLimit
  offset::Int
  order::Vector{SQLOrder}
  group::Vector{SQLColumn}
  having::Vector{SQLWhereEntity}

  SQLQuery(;  columns = SQLColumn[], where = SQLWhereEntity[], limit = SQLLimit("ALL"), offset = 0,
              order = SQLOrder[], group = SQLColumn[], having = SQLWhereEntity[]) =
    new(columns, where, limit, offset, order, group, having)
end

string(q::SQLQuery, m::Type{T}) where {T<:AbstractModel} = to_fetch_sql(m, q)


text:
module SemiseparableMatrices
using LinearAlgebra: BlasFloat
using ArrayLayouts, BandedMatrices, LazyArrays, LinearAlgebra, MatrixFactorizations, LazyBandedMatrices, Base

import Base: size, getindex, setindex!, convert, copyto!
import MatrixFactorizations: QR, QRPackedQ, getQ, getR, QRPackedQLayout, AdjQRPackedQLayout
import LinearAlgebra: qr, qr!, lmul!, ldiv!, rmul!, triu!, factorize, rank
import BandedMatrices: _banded_qr!, bandeddata
import LazyArrays: arguments, applylayout, _cache, CachedArray, CachedMatrix, ApplyLayout, resizedata!
import ArrayLayouts: MemoryLayout, sublayout, sub_materialize, MatLdivVec, materialize!, triangularlayout,
                        triangulardata, zero!, _copyto!, colsupport, rowsupport,
                        _qr, _qr!, _factorize
import LazyBandedMatrices: resize, ScalarOrBandedLayouts, ApplyBandedLayout

export SemiseparableMatrix, AlmostBandedMatrix, LowRankMatrix, ApplyMatrix, ApplyArray, almostbandwidths, almostbandedrank

const LowRankMatrix{T,A,B} = MulMatrix{T,Tuple{A,B}}
LowRankMatrix(A::AbstractArray, B::AbstractArray) = ApplyMatrix(*, A, B)
LowRankMatrix(S::SubArray) = LowRankMatrix(map(Array,arguments(S))...)
LowRankMatrix{T}(::UndefInitializer, (m,n)::NTuple{2,Integer}, r::Integer) where T =
    ApplyMatrix(*, Array{T}(undef, m, r), Array{T}(undef, r, n))
LowRankMatrix(Z::Zeros{T}, r::Integer) where T = ApplyMatrix(*, zeros(T, size(Z,1), r), zeros(T, r, size(Z,2)))
LowRankMatrix{T}(Z::Zeros, r::Integer) where T = ApplyMatrix(*, zeros(T, size(Z,1), r), zeros(T, r, size(Z,2)))

separablerank(A) = size(arguments(ApplyLayout{typeof(*)}(),A)[1],2)

include("SemiseparableMatrix.jl")
include("AlmostBandedMatrix.jl")
include("invbanded.jl")

end # module


text:
# Largely copied from Ti's example in ApproximateGPs.

using ApproximateGPs
using LinearAlgebra
using Distributions
using FillArrays
using Optim
using ParameterHandling
using Zygote

using Plots
default(; legend=:outertopright, size=(700, 400))

using Random
Random.seed!(1);

# Specify parameters.
N_tr = 50;
x_tr = range(-10.0, 10.0; length=N_tr);
θ_init = (
    gp = (
        σ² = positive(1.0),
        l = positive(3.0),
    ),
    α = positive(1.0),
    qu = (
        z = fixed(x_tr),
        m = zeros(N_tr),
        C = positive_definite(Matrix(Eye(N_tr))),
    ),
    jitter = fixed(1e-6),
);

# Specify functions to build LatentGP.
build_gp(θ) = GP(θ.σ² * with_lengthscale(SEKernel(), θ.l))

build_latent_gp(θ) = LatentGP(build_gp(θ.gp), GammaLikelihood(θ.α), θ.jitter)

# Generate synthetic data.

(f_tr, y_tr) = rand(build_latent_gp(ParameterHandling.value(θ_init))(x_tr));

function build_svgp(θ)
    lf = build_latent_gp(θ)
    q = MvNormal(θ.qu.m, θ.qu.C + 1e-9I)
    fz = lf(θ.qu.z).fx
    return SVGP(fz, q), lf
end

function loss(θ)
    svgp, f = build_svgp(θ)
    fx = f(x_tr)
    return -elbo(svgp, fx, y_tr)
end

# Optimise everything.
θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init);

# L-BFGS parameters chosen because they seems to work well empirically.
# You could also try with the defaults.
optimisation_result = optimize(
    loss ∘ unflatten,
    θ -> only(Zygote.gradient(loss ∘ unflatten, θ)),
    θ_flat_init + randn(length(θ_flat_init)),
    LBFGS(;
        alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
        linesearch=Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(; iterations=1_000);
    inplace=false,
);
θ_opt = unflatten(optimisation_result.minimizer);

#
x_pr = range(-15.0, 15.0; length=250);
svgp_opt, lf_opt = build_svgp(θ_opt);

# Abstract these lines? Would be nice to write
# lf_post_opt = posterior(lf_opt, svgp_opt)
# or something.
post_opt = posterior(svgp_opt);
lf_post_opt = LatentGP(post_opt, lf_opt.lik, lf_opt.Σy);

f = [rand(lf_post_opt(x_pr)).f for _ in 1:20];
p = map(_f -> lf_opt.lik.invlink.(_f), f);

let
    plt = plot()
    plot!(plt, x_pr, p; color=:blue, linealpha=0.2, label="")
    plot!(plt, x_pr, mean(p); color=:blue, linalpha=0.5, linewidth=2, label="")
    scatter!(plt, x_tr, y_tr; color=:red, label="")
end


text:
module CLHMM
    import Distributions:Distribution,VariateForm,Univariate,Categorical,logpdf,isprobvec
    import Distributed: RemoteChannel, myid
    import StatsFuns: logsumexp, logaddexp

    export AbstractHMM, HMM, linear_step, linear_hmm_converger!, lin_obs_set_lh, lps

    include("hmm.jl")
    include("linear.jl")
end # module


text:
export primepi
export legendrephi
export nthprime
export primeLi
export primeLiinv
export primepi_xmax
export primepi_num_threads
export primepi_num_threads
export prime_set_print_status
export primepi_test

const libccountname = "libcprimecount.so"

for (f,c) in ( # (:primepi, :(:pi_int64)), use function with keyword
               (:pi_deleglise_rivat, :(:pi_deleglise_rivat)),
              (:pi_legendre, :(:pi_legendre)), (:pi_lehmer, :(:pi_lehmer)),
              (:pi_meissel, :(:pi_meissel)), (:pi_lmo, :(:pi_lmo)),
              (:piprimesieve, :(:pi_primesieve)), (:nthprimecount, :(:nth_prime)),
              (:primeLi, :(:prime_Li)), (:primeLiinv, :(:prime_Li_inverse)))
    @eval begin
        # are c++ exceptions the culprit ? Don't see it in the code
        # function ($f){T<:Real}(n::T)   # try-catch not preventing segfaults
        #     res = try
        #         ccall(($c, libccountname), Int64, (Int64,), convert(Int64,n))
        #     catch
        #         throw(InterruptException())
        #     end
        #     return convert(T,res)
        # end
        ($f){T<:Real}(n::T) = ccall(($c, libccountname), Int64, (Int64,), convert(Int64,n))
        ($f){T<:AbstractString}(n::T) = ($f)(conv128(n))
        # Base.@vectorize_1arg Real $f
        # Base.@vectorize_1arg AbstractString $f
    end
end

function legendrephi(x,a)
    ccall((:prime_phi, libccountname), Int64, (Int64, Int64), convert(Int64,x), convert(Int64,a))
end
#Base.@vectorize_2arg Integer legendrephi

function nthprime(x; alg::Symbol = :count)
    if alg == :count
        return nthprimecount(x)
    elseif alg == :sieve
        return nthprimea(x)
    else error("algorithm must be one of :count, :sieve")
    end
end
#Base.@vectorize_1arg Integer nthprime

# libprimecount has a member function converts a string to Int128, but we probably handle more cases this way
function primepi{T<:AbstractString}(s::T)
    n1 = conv128(s)
    s1 = string(n1)
    parse(Int128, bytestring(ccall((:pi_string,libccountname),Ptr{UInt8},(Ptr{UInt8},),s1)))
end

#Base.@vectorize_1arg AbstractString primepi

# Can't get access to Int128 routine, so we convert back and forth many times.
pi_deleglise_rivat(x::Int128) = primepi(string(x))

# :auto is not perfect, fails often
# Tables or interpolation, or a crude fit for both methods is a better way.
# The two parameters are x and rem, the distance to the previous table value.
function primepi(x; alg::Symbol = :auto)
    x < 2 && return zero(x)
    if alg == :auto
        if (x < 7*10^11)  # This is a fairly sharp crossover here
            return countprimes(x)
        else
            rem = piandrem(x)[3]
            rem == 0 && return countprimes(x)
            if rem <= 10^9 || x/rem >= 10^7  # this gets a lot of cases correctly.
                return countprimes(x)
            else
                return pi_deleglise_rivat(x)
            end
        end
    elseif alg == :deleglise_rivat || alg == :dr
        pi_deleglise_rivat(x)
    elseif alg == :tabsieve
        countprimes(x)
    elseif alg == :lehmer   # The remaining are of academic interest (... well all are, really)
        pi_lehmer(x)
    elseif alg == :meissel
        pi_meissel(x)
    elseif alg == :lmo
        pi_lmo(x)
    elseif alg == :legendre
        pi_legendre(x)
    elseif alg == :sieve
        piprimesieve(x)
    else
        error("Algorithm must be one of :auto, :deleglise_rivat (:dr), :tabsieve, :legendre, " *
               ":lehmer, :meissel, :lmo, :sieve.")
    end
end

# using macro will work, too
#Base.@vectorize_1arg Integer primepi
function primepi(arr::AbstractArray; alg::Symbol = :auto)
    arrout = similar(arr)
    for i in 1:length(arr) arrout[i] = primepi(arr[i]; alg=alg) end
    arrout
end

#register_sigint() = ccall((:cprimecount_register_sigint, libccountname), Void, ())

# Sun Apr  3 18:21:07 CEST 2016:  These two are broken. This syntax worked at one point. Julia has changed
#primepi_xmax() = parse(Int128, bytestring(ccall((:pi_xmax, libccountname), Ptr{UInt8}, ())))
#const PRIMEPI_XMAX = primepi_xmax()

primepi_num_threads(n) = ccall((:prime_set_num_threads,libccountname),Void,(Int,), convert(Int,n))
primepi_num_threads() = ccall((:prime_get_num_threads,libccountname),Int,())
prime_set_print_status(stat::Bool) = ccall((:prime_set_num_threads,libccountname),Void,(Int,), stat ? 1 : 0)
primepi_test() = ccall((:prime_test,libccountname),Int,())

# libprimecount does not really set the number of threads until you ask.
# In primecount, the initial value returned by init_num_threads == num cores.
# In primesieve it was a huge integer. I am paranoid, so I check.
function init_num_threads()
    n = primepi_num_threads()
    n < 10000 && primepi_num_threads(n)
end

# Set multithreading now
init_num_threads()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
import Pkg: PackageSpec

name = "libcgal_julia"
rversion = v"0.17.0"
version = VersionNumber(rversion.major,
                        rversion.minor,
                        100rversion.patch + julia_version.minor)

isyggdrasil = get(ENV, "YGGDRASIL", "") == "true"
rname = "libcgal-julia"

# Collection of sources required to build CGAL
sources = [
    isyggdrasil ?
        GitSource("https://github.com/rgcv/$rname.git",
                  "91ab24f45b689a5fcd760db068f0c8fd06539744") :
        DirectorySource(joinpath(ENV["HOME"], "src/github/rgcv/$rname"))
]

# Bash recipe for building across all platforms
jlcgaldir = ifelse(isyggdrasil, rname, ".")
script = raw"""
## pre-build setup
# exit on error
set -eu

macosflags=
case $target in
  *apple-darwin*)
    macosflags="-DCMAKE_CXX_COMPILER_ID=AppleClang"
    macosflags="$macosflags -DCMAKE_CXX_COMPILER_VERSION=10.0.0"
    macosflags="$macosflags -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
    ;;
esac
""" * """
## configure build
cmake $jlcgaldir """ * raw"""\
  -B /tmp/build \
  `# cmake specific` \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_FIND_ROOT_PATH=$prefix \
  -DCMAKE_INSTALL_PREFIX=$prefix \
  -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN \
  `# tell jlcxx where julia is` \
  -DJulia_PREFIX=$prefix \
  $macosflags

## and away we go..
VERBOSE=ON cmake --build /tmp/build --config Release --target install -- -j$nproc
""" * """
install_license $jlcgaldir/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
# generates an abundance of linker errors and notes about using older versions
# of GCC.  Among many things, this could be related to boost as well.  However,
# requiring newer versions would, much like libsingular_julia, require the
# dropping of older julia versions <1.6.  CGAL also contains several deprecation
# notices when configured against newer versions of boost, so it's probably best
# to avoid it for now as well.
filter!(p -> arch(p) ≠ "armv7l", platforms)
# filter experimental platforms
filter!(p -> arch(p) ≠ "armv6l", platforms)
filter!(p -> !(Sys.isapple(p) && arch(p) == "aarch64"), platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcgal_julia_exact", :libcgal_julia_exact),
    LibraryProduct("libcgal_julia_inexact", :libcgal_julia_inexact),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    BuildDependency(PackageSpec(name="GMP_jll", version=v"6.1.2")),
    BuildDependency(PackageSpec(name="MPFR_jll", version=v"4.0.2")),

    Dependency("CGAL_jll", compat="~5.3"),
    Dependency("libcxxwrap_julia_jll", VersionNumber(0, 8, julia_version.minor)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=gcc_version,
               julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 46f59f50-7d41-11eb-28e1-8fc4afcc6837
using Optim

# ╔═╡ c60c577a-7d41-11eb-2390-79a7c6acb0d0
f(x) = (1.0 - x[1]) ^2 + 100.0  * (x[2] - x[1]^2)^2

# ╔═╡ 028dbf9a-7d42-11eb-218b-657847801c0f
x0 = [0.0 0.0]

# ╔═╡ 0b11dd5e-7d42-11eb-15e9-bde4ac3d7662
optimize(f, x0)

# ╔═╡ 24734878-7d42-11eb-10e7-35421889faaf

# ╔═╡ 167d66f2-7d42-11eb-1967-3fff5fa301dd

# ╔═╡ 11e7141e-7d42-11eb-2d63-7155e958f07d

# ╔═╡ Cell order:
# ╠═46f59f50-7d41-11eb-28e1-8fc4afcc6837
# ╠═c60c577a-7d41-11eb-2390-79a7c6acb0d0
# ╠═028dbf9a-7d42-11eb-218b-657847801c0f
# ╠═0b11dd5e-7d42-11eb-15e9-bde4ac3d7662
# ╠═24734878-7d42-11eb-10e7-35421889faaf
# ╠═167d66f2-7d42-11eb-1967-3fff5fa301dd
# ╠═11e7141e-7d42-11eb-2d63-7155e958f07d


text:
immutable UTURUniChannelCSK <: DigitalUTURUniChannelSystem
  coherent::Bool
  carrier::Carrier
  noise::Distribution
  decoder::Decoder

  function UTURUniChannelCSK(coh::Bool, c::Carrier, n::Distribution, d::Decoder)
    if isa(n, Normal)
      @assert n.μ == 0. "Normal noise with μ!=0 not supported in UTURUniChannelCSK."

      if coh
        @assert isa(d, CorDecoder) "$d decoder not supported in coherent UTURUniChannelCSK."
      else
        @assert isa(d, CorDecoder) || isa(d, MCMLDecoder) "$d decoder not supported in non-coherent UTURUniChannelCSK."
      end
    else
      throw(ArgumentError("$n channel noise not supported in UTURUniChannelCSK."))
    end

    new(coh, c, n, d)
  end
end

UTURUniChannelCSK(coh::Bool, c::Carrier; noise::Distribution=Normal(), decoder::Decoder=CorDecoder()) =
  UTURUniChannelCSK(coh, c, noise, decoder)

UTURUniChannelCSK(; coherent::Bool=false, carrier::Carrier=LogisticCarrier(5), noise::Distribution=Normal(),
  decoder::Decoder=CorDecoder()) = UTURUniChannelCSK(coherent, carrier, noise, decoder)

uturunichannelcsk_snrdb2var(s::Float64; carrier::Carrier=LogisticCarrier(5)) = var(carrier)/(10.0^(s/10.0))

uturunichannelcsk_var2snrdb(v::Float64; carrier::Carrier=LogisticCarrier(5)) = 10.0*log10(var(carrier)/v)

uturunichannelcsk_ebn0db2var(eb::Float64; carrier::Carrier=LogisticCarrier(5)) =
  carrier.len*var(carrier)/(10.0^(eb/10.0))

uturunichannelcsk_var2ebn0db(v::Float64; carrier::Carrier=LogisticCarrier(5)) = 10.0*log10(carrier.len*var(carrier)/v)

uturunichannelcsk_ebn0db2snrdb(eb::Float64; carrier::Carrier=LogisticCarrier(5)) = eb-10.0*log10(carrier.len)

uturunichannelcsk_snrdb2ebn0db(s::Float64; carrier::Carrier=LogisticCarrier(5)) = s+10.0*log10(carrier.len)

function logmcml(d::UTURUniChannelCSKMCMLDecoder, v::Float64, t::Vector{Float64})
  logsumexp(-0.5*t/v)-d.carrier.len*(log(v)+log2π)-log(d.nmc)
end

function gradlogmcml(d::UTURUniChannelCSKMCMLDecoder, v::Float64, t::Vector{Float64})
  expterms = exp(-0.5*t/v)
  (0.5*dot(t, expterms)/(sum(expterms)*v)-d.carrier.len)/v
end

function decode(d::UTURUniChannelCSKMCMLDecoder, r1::Vector{Float64}, r2::Vector{Float64})
  tp1, tm1 = Array(Float64, d.nmc), Array(Float64, d.nmc)
  x = Array(Float64, d.carrier.len)
  local maxfp1, maxxp1, retp1, maxfm1, maxxm1, retm1
  opt = Array(Opt, 2)

  for i = 1:d.nmc
    generate!(d.carrier, x)
    subtract!(x, mean(d.carrier))
    dotr2 = sumsqdiff(r2, x)
    tp1[i] = sumsqdiff(r1, x)+dotr2
    tm1[i] = sumsq(r1+x)+dotr2
  end

  function objective_fp1(v::Vector{Float64}, grad::Vector{Float64})
    if length(grad) > 0
      grad[1] = gradlogmcml(d, v[1], tp1)
    end
    logmcml(d, v[1], tp1)
  end

  opt[1] = convert(Opt, d.opt[1])
  max_objective!(opt[1], objective_fp1)

  try
    (maxfp1, maxxp1, retp1) = optimize(opt[1], [d.init[1]])
  catch
    return -3
  end

  if in(retp1, opt_failure)
    return -4
  end

  function objective_fm1(v::Vector{Float64}, grad::Vector{Float64})
    if length(grad) > 0
      grad[1] = gradlogmcml(d, v[1], tm1)
    end
    logmcml(d, v[1], tm1)
  end

  opt[2] = convert(Opt, d.opt[2])
  max_objective!(opt[2], objective_fm1)

  try
    (maxfm1, maxxm1, retm1) = optimize(opt[2], [d.init[2]])
  catch
    return -5
  end

  if in(retm1, opt_failure)
    return -6
  end

  dec = maxfp1-maxfm1

  if dec > 0
    return 1
  elseif dec < 0
    return -1
  elseif dec == 0
    return -2
  end
end

function sim_sys(s::UTURUniChannelCSK, bit::Int)
  @assert in(bit, [-1, 1]) "Bit must be either -1 or 1."

  x = Array(Float64, s.carrier.len)
  generate!(s.carrier, x)
  subtract!(x, mean(s.carrier))
  r1 = bit*x+rand(s.noise, s.carrier.len)
  r2 = s.coherent ? x : x+rand(s.noise, s.carrier.len)
  decode(s.decoder, r1, r2)
end

function uturunichannelcsk_gen_sys(ebn0db::Ranges{Float64}, sprlen::Ranges{Int}; coherent::Bool=true,
  carrier::Carrier=LogisticCarrier(), decoder::Decoder=CorDecoder())
  elen, slen = length(ebn0db), length(sprlen)
  systems = Array(UTURUniChannelCSK, elen, slen)
  local d

  for i = 1:elen
    for j = 1:slen
      ckwargs = Dict()
      for k in typeof(carrier).names
        ckwargs[k] = k != :len ? carrier.(k) : sprlen[j]
      end
      c = typeof(carrier)(; ckwargs...)

      v = ebn0db2var(ebn0db[i]; system=:UTURUniChannelCSK, carrier=c)

      if isa(decoder, CorDecoder)
        d = decoder
      elseif isa(decoder, MCMLDecoder)
        dkwargs = Dict()
        for k in typeof(decoder).names
          dkwargs[k] = k != :carrier ? decoder.(k) : c
        end
        d = mcml_decoder(system=:UTURUniChannelCSK; dkwargs...)
      else
        error("$decoder is not a supported type of decoder for the uturunichannelcsk_gen_sys function.")
      end

      systems[i, j] = UTURUniChannelCSK(coherent, c; noise=Normal(0., sqrt(v)), decoder=d)
    end
  end

  return systems
end

function sim_ber(s::Vector{UTURUniChannelCSK}, bit::Int, n::Int64; verbose::Bool=false)
  slen = length(s)
  nbiterrors, ndecfails, bers = Array(Int64, slen), Array(Int64, slen), Array(Float64, slen)
  local meter

  if verbose
    println("Execution started at ", now(), "\n\nSimulating BER...")
    meter = ProgressMeter.Progress(slen, "  ")
  end

  for i = 1:slen
    try
      nbiterrors[i], ndecfails[i], bers[i] = sim_ber(s[i], bit, n)
    catch
      nbiterrors[i] = ndecfails[i] = bers[i] = NaN
    end

    if verbose; ProgressMeter.next!(meter); end
  end

  if verbose; println("\nExecution completed at ", now()); end

  return nbiterrors, ndecfails, bers
end

function sim_ber(s::Matrix{UTURUniChannelCSK}, bit::Int, n::Int64; verbose::Bool=false)
  nrows, ncols = size(s)
  nbiterrors, ndecfails, bers = Array(Int64, nrows, ncols), Array(Int64, nrows, ncols), Array(Float64, nrows, ncols)
  local meter

  if verbose
    println("Execution started at ", now(), "\n\nSimulating BER...")
    meter = ProgressMeter.Progress(nrows*ncols, "  ")
  end

  for i = 1:nrows
    for j = 1:ncols
      try
        nbiterrors[i, j], ndecfails[i, j], bers[i, j] = sim_ber(s[i, j], bit, n)
      catch
        nbiterrors[i, j] = ndecfails[i, j] = bers[i, j] = NaN
      end

      if verbose; ProgressMeter.next!(meter); end
    end
  end

  if verbose; println("\nExecution completed at ", now()); end

  return nbiterrors, ndecfails, bers
end

function psim_ber(s::Vector{UTURUniChannelCSK}, bit::Int, n::Int64; ptype::Symbol=:pfor, verbose::Bool=false)
  slen = length(s)
  nbiterrors, ndecfails, bers = Array(Int64, slen), Array(Int64, slen), Array(Float64, slen)
  local meter

  if verbose
    println("Execution started at ", now(), "\n\nSimulating BER...")
    meter = ProgressMeter.Progress(slen, "  ")
  end

  for i = 1:slen
    try
      nbiterrors[i], ndecfails[i], bers[i] = psim_ber(s[i], bit, n; ptype=ptype)
    catch
      nbiterrors[i] = ndecfails[i] = bers[i] = NaN
    end

    if verbose; ProgressMeter.next!(meter); end
  end

  if verbose; println("\nExecution completed at ", now()); end

  return nbiterrors, ndecfails, bers
end

function psim_ber(s::Matrix{UTURUniChannelCSK}, bit::Int, n::Int64; ptype::Symbol=:pfor, verbose::Bool=false)
  nrows, ncols = size(s)
  nbiterrors, ndecfails, bers = Array(Int64, nrows, ncols), Array(Int64, nrows, ncols), Array(Float64, nrows, ncols)
  local meter

  if verbose
    println("Execution started at ", now(), "\n\nSimulating BER...")
    meter = ProgressMeter.Progress(nrows*ncols, "  ")
  end

  for i = 1:nrows
    for j = 1:ncols
      try
        nbiterrors[i, j], ndecfails[i, j], bers[i, j] = psim_ber(s[i, j], bit, n; ptype=ptype)
      catch
        nbiterrors[i, j] = ndecfails[i, j] = bers[i, j] = NaN
      end

      if verbose; ProgressMeter.next!(meter); end
    end
  end

  if verbose; println("\nExecution completed at ", now()); end

  return nbiterrors, ndecfails, bers
end

function ber_lb(s::UTURUniChannelCSK; btype::Symbol=:jensen)
  if btype == :jensen
    if isa(s.noise, Normal) && s.noise.μ == 0.0
      cdf(Normal(), -sqrt(s.carrier.len*var(s.carrier))/s.noise.σ)
    else
      error("Jensen's BER lower bound defined for UTURUniChannelCSK only for Normal noise with zero mean.")
    end
  else
    error("$btype type of BER lower bound not known for UTURUniChannelCSK.")
  end
end


text:
using DashBootstrapComponents

modal = html_div([
    dbc_button("Open modal", id = "open-fs"),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Fullscreen modal")),
            dbc_modalbody("Wow this thing takes up a lot of space..."),
        ],
        id = "modal-fs",
        fullscreen = true,
    ),
])

callback!(
    app,
    Output("modal-fs", "is_open"),
    Input("open-fs", "n_clicks"),
    State("modal-fs", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function getDiagnostics(filename)
    lines = readlines(filename)
    ndigits = length(lines[1])
    nums = parse.(Int32, lines, base=2)
    counts = zeros(Int32, ndigits)
    for n in nums
        for i in 1:ndigits
            pow = 2^(i-1)
            d = n & pow
            counts[i] += (d/pow)
        end
    end
    gamma = counts .> length(lines)/2
    g = 0
    for i in 1:length(gamma)
        g += gamma[i] * 2^(i-1)
    end
    epsilon = counts .< length(lines)/2
    e = 0
    for i in 1:length(epsilon)
        e += epsilon[i] * 2^(i-1)
    end
    return g * e
end

function nthDigit(num, digit, ndigits)
    pow = 2^(ndigits-digit)
    return convert(Int32, (num & pow) / pow)
end

function getRatings(filename)
    lines = readlines(filename)
    ndigits = length(lines[1])
    nums = parse.(Int32, lines, base=2)
    nums2 = copy(nums)
    round = 1
    while length(nums) > 1
        count = 0
        for n in nums
            count += nthDigit(n, round, ndigits)
        end
        oneorzero = convert(Int32, count >= length(nums)/2)
        filter!( (x) -> nthDigit(x, round, ndigits) == oneorzero, nums)
        @show nums
        round += 1
    end
    oxygen = nums[1]
    nums = nums2
    # shameless copy-paste here
    round = 1
    while length(nums) > 1
        count = 0
        for n in nums
            count += nthDigit(n, round, ndigits)
        end
        oneorzero = convert(Int32, count >= length(nums)/2)
        filter!( (x) -> nthDigit(x, round, ndigits) != oneorzero, nums)
        @show nums
        round += 1
    end
    co2 = nums[1]
    return(oxygen * co2)
end


text:
using Test
using SimplePartitions
using Permutations

@testset "Set Partitions" begin
n = 20
p = RandomPermutation(n)
P = Partition(p)
@test P >= Partition(p*p')
@test length(cycles(p)) == num_parts(P)
@test num_elements(P) == n
@test ground_set(P) == Set(1:n)

Plist = collect(parts(P))
for A in Plist
    @test in(A,P)
end

@test P == PartitionBuilder(parts(P))

P = Partition(n)
for k=1:n-1
    merge_parts!(P,k,k+1)
end
@test num_parts(P) == 1

@test length(all_partitions(Set(1:4))) == 15

P = Partition(RandomPermutation(10))
Q = Partition(RandomPermutation(10))

@test P ∧ Q <= P
@test P ∨ Q >= Q

end  # end testset

@testset "Integer Partitions" begin

p = IntegerPartition(1,2,3)
@test sum(p)==6
@test p' == p
@test sum(p+p) == 2*sum(p)

p = IntegerPartition(5)
@test string(p) == "(5)"
@test parts(p') == ones(Int,5)

end # end testset


text:
function forward{T}(::typeof(batchnorm), x::CuArray{T}, scale::CuArray, bias,
    mean, var; epsilon=0.001)

    h = CUDNN.handle(x)
    y = similar(x)
    xdesc = CUDNN.TensorDesc(x)
    sbmvdesc = CUDNN.TensorDesc(scale)
    if inference
        cudnnBatchNormalizationForwardInference(h, mode, T[1], T[0], xdesc, x,
            xdesc, y, sbmvdesc, scale, bias, estimated_mean, estimated_var, Cdouble(epsilon))
    else # training
        running_mean = zeros(scale)
        running_var = zeros(scale)
        save_mean = similar(scale)
        save_invvar = similar(scale)
        cudnnBatchNormalizationForwardTraining(h, mode, T[1], T[0], xdesc, x, ydesc, y,
            sbmvdesc, scale, bias, Cdouble(averagefactor), running_mean, running_var,
    	    Cdouble(epsilon), save_mean, save_invvar)
        function backward!(gy, gx)
            if !isvoid(gx)
                cudnnBatchNormalizationBackward(h, mode, T[1], T[0], T[1], T[0],
                    xdesc, x, dydesc, dy, dxdesc, dx, dsbdesc, scale, dscale, dbias,
                    Cdouble(epsilon), saved_mean, saved_invvar)
            end
        end
    end
    y, backward!
end

function batchnorm_training(mode, x, scale, bias, averagefactor::Float64, epsilon::Float64)
    T = eltype(x)
    h = handle(x)
    y = similar(x)
    xdesc = TensorDesc(x)
    ydesc = TensorDesc(y)
    sbmvdesc = TensorDesc(scale)
    running_mean = zeros(scale)
    running_var = zeros(scale)
    save_mean = similar(scale)
    save_invvar = similar(scale)

    cudnnBatchNormalizationForwardTraining(h, mode, T[1], T[0], xdesc, x, ydesc, y,
        sbmvdesc, scale, bias, Cdouble(averagefactor), running_mean, running_var,
	    Cdouble(epsilon), save_mean, save_invvar)

    y, running_mean, running_var, save_mean, save_invvar
end

function ∇batchnorm!(mode, x, dy, dx, scale, dscale, dbias, epsilon::Float64, saved_mean, saved_invvar)
    T = eltype(x)
    h = handle(x)
    xdesc = tensor_desc(x)
    dydesc = tensor_desc(dy)
    dxdesc = tensor_desc(dx)
    dsbdesc = tensor_desc(scale)

    cudnnBatchNormalizationBackward(h, mode, T[1], T[0], T[1], T[0],
    xdesc, x, dydesc, dy, dxdesc, dx, dsbdesc, scale, dscale, dbias,
        Cdouble(epsilon), saved_mean, saved_invvar)
end


text:
# Henon fractal zoom
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

fig, axes = subplots(1,3, figsize = (figx,figx/4))
botrow = axes

he = Systems.henon()
N = 10000
tr = trajectory(he, N, Ttr = 100)
integ = integrator(he)
data1 = Dataset([integ.u])
data2 = Dataset([integ.u])
n, m = 1, 1

zbox1 = ((0.93, 0.05), (1.17, 0.13))
zbox2 = ((1.1275, 0.078), (1.1425, 0.085))
while n < N
    step!(integ)
    if zbox1[1][1] < integ.u[1] < zbox1[2][1] && zbox1[1][2] < integ.u[2] < zbox1[2][2] && m < N
        push!(data1, integ.u)
        global m+=1
    end
    if zbox2[1][1] < integ.u[1] < zbox2[2][1] && zbox2[1][2] < integ.u[2] < zbox2[2][2]
        push!(data2, integ.u)
        global n+=1
    end
end

for ax in botrow[1:2]
    ax.clear()
    ax.plot(tr[:, 1], tr[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
end
botrow[1].set_xticks([])
botrow[1].set_yticks([])
botrow[2].clear()
botrow[2].plot(data1[:, 1], data1[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
botrow[2].axis("off")
botrow[3].clear()
botrow[3].plot(data2[:, 1], data2[:, 2], ls = "None", marker = ".", color = COLORS[1], ms = 1)
botrow[3].axis("off")

for (i, zbox) in enumerate((zbox1, zbox2))
    axis_zoomin!(botrow[i+1], botrow[i], zbox, zbox, COLORS[[4,3][i]])
    botrow[i+1].set_xlim(zbox[1][1], zbox[2][1])
    botrow[i+1].set_ylim(zbox[1][2], zbox[2][2])
end

fig.tight_layout(pad = 0.3)
fig.subplots_adjust(top = 0.98, bottom = 0.05, right = 0.95, left = 0.05, wspace=0.1, hspace = 0.1)
wsave(plotsdir("5", "henon_zoom"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using StationXML, Test
import EzXML
using EzXML: ElementNode

@testset "Utils" begin
    @testset "XML escaping" begin
        let s = "a&b<c>d"e'f", s′ = "a&b<c>d\"e'f"
            @test StationXML.xml_unescape(s) == s′
            @test StationXML.xml_escape(s′) == s
        end
    end

    @testset "Name transform" begin
        let f = StationXML.transform_name
            @test f("Module") == :module_name
            @test f("ModuleURI") == :module_uri
            @test f("Email") == :email
            @test f("SelectedNumberChannels") == :selected_number_channels
            @test f("FrequencyDBVariation") == :frequency_db_variation
            @test f("URI") == :uri
            @test f("sourceID") == :source_id
            @test f("measurementMethod") == :measurement_method
            @test f("end") == :end_
            @test f("maximumTimeTear") == :maximum_time_tear
            @test f("type") == :type
        end
        let f = StationXML.xml_element_name
            @test f(:module_name) == "Module"
            @test f(:module_uri) == "ModuleURI"
            @test f(:email) == "Email"
            @test f(:selected_number_channels) == "SelectedNumberChannels"
            @test f(:fir) == "FIR"
            @test f(:input_units) == "InputUnits"
            @test f(:frequency_db_variation) == "FrequencyDBVariation"
            @test f(:uri) == "URI"
        end
        let f = StationXML.xml_attribute_name
            @test f(:resource_id) == "resourceId"
            @test f(:i) == "i"
            @test f(:start_date) == "startDate"
            @test f(:code) == "code"
            @test f(:restricted_status) == "restrictedStatus"
            @test f(:number) == "number"
            @test f(:source_id) == "sourceID"
            @test f(:measurement_method) == "measurementMethod"
            @test f(:start) == "start"
            @test f(:end_) == "end"
            @test f(:maximum_time_tear) == "maximumTimeTear"
            @test f(:type) == "type"
            @test_throws ArgumentError f(:made_up_attribute)
        end
        # Round-trip all names in the spec
        @testset "Specification v$(version)" for version in ("1.0", "1.1")
            let file = joinpath(@__DIR__, "data", "fdsn-station-$(version).xsd.xml"),
                    xml = EzXML.root(EzXML.readxml(file))
                # Elements
                for node in EzXML.findall("//xs:element", xml,
                        ["xs"=>"http://www.w3.org/2001/XMLSchema",
                        "fsx"=>"http://www.fdsn.org/xml/station/1"])
                    for att in EzXML.eachattribute(node)
                        if att.name == "name" && att.content != "FDSNStationXML"
                            @test StationXML.xml_element_name(StationXML.transform_name(att.content)) == att.content
                        end
                    end
                end
                # Attributes
                for node in EzXML.findall("//xs:attribute", xml,
                        ["xs"=>"http://www.w3.org/2001/XMLSchema",
                        "fsx"=>"http://www.fdsn.org/xml/station/1"])
                    for att in EzXML.eachattribute(node)
                        if att.name == "name"
                            @test StationXML.xml_attribute_name(StationXML.transform_name(att.content)) == att.content
                        end
                    end
                end
            end
        end
    end

    @testset "Enumeration macro" begin
        @eval StationXML.@enumerated_struct(ExampleStruct, ("A", "B", "D"))
        # FIXME: Work out correct escaping so permitted_values is evaluated
        #        in the correct scope (i.e., always StationXML)
        @test Main.permitted_values(ExampleStruct) == ("A", "B", "D")
        @test_throws LoadError @eval StationXML.@enumerated_struct(ExampleStruct2, ())
        @test_throws LoadError @eval StationXML.@enumerated_struct(ExampleStruct3, 1.0)
        @test isdefined(Main, :ExampleStruct)
        @test ExampleStruct("A") == ExampleStruct(value="A") isa ExampleStruct
        @test_throws ArgumentError ExampleStruct("C")
        @test convert(String, ExampleStruct("A")) == "A"
        @test convert(SubString{String}, ExampleStruct("B")) == SubString{String}("B")
        @test convert(ExampleStruct, "A") == ExampleStruct("A")
        @test convert(ExampleStruct, SubString{String}("B")) == ExampleStruct("B")
        @test parse(ExampleStruct, "A") == ExampleStruct("A")
        @test StationXML.local_parse(ExampleStruct, "B") == ExampleStruct("B")
        node = ElementNode("nodeName")
        node.content = "B"
        @test StationXML.parse_node(ExampleStruct, node, true) == ExampleStruct("B")
        @test StationXML.attribute_fields(ExampleStruct) == ()
        @test StationXML.element_fields(ExampleStruct) == ()
        @test StationXML.has_text_field(ExampleStruct) == true
        @test StationXML.text_field(ExampleStruct) == :value
    end
end


text:
# Following the approach in JuMP
# Must support interactions between the Following
# 1. Numbers
# 2. Quantities
# 3. UnitVariableRef
# 4. VariableRef
# 5. UnitExpression

# Unitful.Quantity/Number -- VariableRef/UnitVariableRef

# Use the Unitful.unit of the lhs when adding or subtracting different quantities
function Base.:+(lhs::Unitful.Quantity, rhs::UnitVariableRef)
    factor = Unitful.ustrip(
        Unitful.uconvert(Unitful.unit(lhs), Unitful.Quantity(1.0, rhs.unit)),
    )
    return UnitExpression(
        Unitful.ustrip(lhs) + factor * rhs.variable,
        Unitful.unit(lhs),
    )
end

function Base.:+(lhs::UnitVariableRef, rhs::Unitful.Quantity)
    rhsval = Unitful.ustrip(Unitful.uconvert(lhs.unit, rhs))
    return UnitExpression(lhs.variable + rhsval, lhs.unit)
end

function Base.:*(lhs::Unitful.Quantity, rhs::VariableRef)
    return UnitExpression(Unitful.ustrip(lhs) * rhs, Unitful.unit(lhs))
end

function Base.:*(lhs::Number, rhs::UnitVariableRef)
    return UnitExpression(lhs * rhs.variable, rhs.unit)
end

function Base.:*(lhs::Unitful.Quantity, rhs::UnitVariableRef)
    val = lhs * Unitful.Quantity(1.0, rhs.unit)
    return UnitExpression(Unitful.ustrip(val) * rhs.variable, Unitful.unit(val))
end

Base.:-(lhs::UnitVariableRef) = UnitExpression(-lhs.variable, lhs.unit)
Base.:-(lhs::Unitful.Quantity, rhs::UnitVariableRef) = (+)(lhs, -rhs)
Base.:-(lhs::UnitVariableRef, rhs::Unitful.Quantity) = (+)(lhs, -rhs)

Base.:*(lhs::UnitVariableRef, rhs::Number) = (*)(rhs, lhs)
Base.:*(lhs::VariableRef, rhs::Unitful.Quantity) = (*)(rhs, lhs)
Base.:*(lhs::UnitVariableRef, rhs::Unitful.Quantity) = (*)(rhs, lhs)

Base.:/(lhs::UnitVariableRef, rhs::Number) = (*)(1.0 / rhs, lhs)
Base.:/(lhs::VariableRef, rhs::Unitful.Quantity) = (*)(1.0 / rhs, lhs)
Base.:/(lhs::UnitVariableRef, rhs::Unitful.Quantity) = (*)(1.0 / rhs, lhs)

# UnitVariableRef -- UnitVariableRef
# Use the Unitful.unit of the lhs term
function Base.:+(lhs::UnitVariableRef, rhs::UnitVariableRef)
    factor = Unitful.ustrip(
        Unitful.uconvert(lhs.unit, Unitful.Quantity(1.0, rhs.unit)),
    )
    return UnitExpression(lhs.variable + factor * rhs.variable, lhs.unit)
end

Base.:-(lhs::UnitVariableRef, rhs::UnitVariableRef) = (+)(lhs, -rhs)

# Unitful.Quantity/Number -- UnitExpression
# Use the Unitful.unit of the lhs term
function Base.:+(lhs::Unitful.Quantity, rhs::UnitExpression)
    factor = Unitful.ustrip(
        Unitful.uconvert(Unitful.unit(lhs), Unitful.Quantity(1.0, rhs.unit)),
    )
    return UnitExpression(
        Unitful.ustrip(lhs) + factor * rhs.expr,
        Unitful.unit(lhs),
    )
end

function Base.:+(lhs::UnitExpression, rhs::Unitful.Quantity)
    rhsval = Unitful.ustrip(Unitful.uconvert(lhs.unit, rhs))
    return UnitExpression(lhs.expr + rhsval, lhs.unit)
end

function Base.:*(lhs::Unitful.Quantity, rhs::UnitExpression)
    factor = lhs * Unitful.Quantity(1.0, rhs.unit)
    return UnitExpression(
        Unitful.ustrip(factor) * rhs.expr,
        Unitful.unit(factor),
    )
end

function Base.:*(lhs::Number, rhs::UnitExpression)
    return UnitExpression(lhs * rhs.expr, rhs.unit)
end

Base.:-(lhs::UnitExpression) = UnitExpression(-lhs.expr, lhs.unit)
Base.:-(lhs::Unitful.Quantity, rhs::UnitExpression) = (+)(lhs, -rhs)
Base.:-(lhs::UnitExpression, rhs::Unitful.Quantity) = (+)(lhs, -rhs)

Base.:*(lhs::UnitExpression, rhs::Unitful.Quantity) = (*)(rhs, lhs)
Base.:/(lhs::UnitExpression, rhs::Unitful.Quantity) = (*)(1.0 / rhs, lhs)
Base.:*(lhs::UnitExpression, rhs::Number) = (*)(rhs, lhs)
Base.:/(lhs::UnitExpression, rhs::Number) = (*)(1.0 / rhs, lhs)

# UnitVariableRef -- UnitExpression
# Use the Unitful.unit of the lhs term
function Base.:+(lhs::UnitVariableRef, rhs::UnitExpression)
    factor = Unitful.ustrip(
        Unitful.uconvert(lhs.unit, Unitful.Quantity(1.0, rhs.unit)),
    )
    return UnitExpression(lhs.variable + factor * rhs.expr, lhs.unit)
end

function Base.:+(lhs::UnitExpression, rhs::UnitVariableRef)
    factor = Unitful.ustrip(
        Unitful.uconvert(lhs.unit, Unitful.Quantity(1.0, rhs.unit)),
    )
    return UnitExpression(lhs.expr + factor * rhs.variable, lhs.unit)
end

Base.:-(lhs::UnitVariableRef, rhs::UnitExpression) = (+)(lhs, -rhs)
Base.:-(lhs::UnitExpression, rhs::UnitVariableRef) = (+)(lhs, -rhs)

# UnitExpression -- UnitExpression
# Use the Unitful.unit of the lhs term
function Base.:+(lhs::UnitExpression, rhs::UnitExpression)
    factor = Unitful.ustrip(
        Unitful.uconvert(lhs.unit, Unitful.Quantity(1, rhs.unit)),
    )
    aff = JuMP.add_to_expression!(copy(lhs.expr), factor, rhs.expr)
    return UnitExpression(aff, lhs.unit)
end

Base.:-(lhs::UnitExpression, rhs::UnitExpression) = (+)(lhs, -rhs)

# Do not allow combination of terms with or without units
function Base.:+(lhs::_UnitAffOrVar, rhs::_AffOrVar)
    return error("Can not combine variables with or without units")
end
function Base.:+(lhs::_UnitAffOrVar, rhs::Real)
    return error("Can not combine unit expression with parameter without unit")
end
Base.:+(lhs::_AffOrVar, rhs::_UnitAffOrVar) = (+)(rhs, lhs)
Base.:+(lhs::Real, rhs::_UnitAffOrVar) = (+)(rhs, lhs)
Base.:-(lhs::_UnitAffOrVar, rhs::_AffOrVar) = (+)(lhs, -rhs)
Base.:-(lhs::_AffOrVar, rhs::_UnitAffOrVar) = (+)(lhs, -rhs)
Base.:-(lhs::_UnitAffOrVar, rhs::Real) = (+)(lhs, -rhs)
Base.:-(lhs::Real, rhs::_UnitAffOrVar) = (+)(lhs, -rhs)


text:
export jd
export jd_betcke
export jd_effenberger

using IterativeSolvers
using LinearAlgebra
using Random
using ..NEPTypes:DeflatedNEP
using ..NEPTypes:deflated_nep_compute_Q

"""
 For Jacobi-Davidson type solvers see
[`jd_betcke`](@ref) or
[`jd_effenberger`](@ref).
"""
function jd()
    error("The function `jd` is not implemented. See `jd_betcke` or jd_effenberger`.")
end

"""
    jd_betcke([eltype]], nep::ProjectableNEP; [neigs=1], [tol=eps(real(T))*100], [maxit=100], [λ=zero(T)], [orthmethod=DGKS],  [errmeasure], [linsolvercreator=DefaultLinSolverCreator()], [v = randn(size(nep,1))], [logger=0], [inner_logger=0], [inner_solver_method=DefaultInnerSolver()], [projtype=:PetrovGalerkin], [target=zero(T)])
The function computes eigenvalues using Jacobi-Davidson method, which is a projection method.
The projected problems are solved using a solver spcified through the type `inner_solver_method`.
The logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.
For numerical stability the basis is kept orthogonal, and the method for orthogonalization is specified by `orthmethod`, see the package `IterativeSolvers.jl`.
The function tries to compute `neigs` number of eigenvalues, and throws a `NoConvergenceException` if it cannot.
The value `λ` and the vector `v` are initial guesses for an eigenpair. `linsolvercreator` is a function which specifies how the linear system is created and solved.
The `target` is the center around which eiganvlues are computed.
By default the method uses a Petrov-Galerkin framework, with a trial (left) and test (right) space, hence ``W^H T(λ) V`` is the projection considered. By specifying  `projtype` to be `:Galerkin` then `W=V`.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("dep0",50);
julia> λ,v=jd_betcke(nep,tol=1e-8,maxit=20);
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
1.9570222439914163e-10
```

# References
* T. Betcke and H. Voss, A Jacobi-Davidson-type projection method for nonlinear eigenvalue problems. Future Gener. Comput. Syst. 20, 3 (2004), pp. 363-372.
* H. Voss, A Jacobi–Davidson method for nonlinear eigenproblems. In: International Conference on Computational Science. Springer, Berlin, Heidelberg, 2004. pp. 34-41.
See also
* C. Effenberger, Robust successive computation of eigenpairs for nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl. 34, 3 (2013), pp. 1231-1256.
"""
jd_betcke(nep::NEP;kwargs...) = jd_betcke(ComplexF64,nep;kwargs...)
function jd_betcke(::Type{T},
                   nep::ProjectableNEP;
                   maxit::Int = 100,
                   neigs::Int = 1,
                   projtype::Symbol = :PetrovGalerkin,
                   inner_solver_method = DefaultInnerSolver(),
                   orthmethod::Type{T_orth} = IterativeSolvers.DGKS,
                   errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                   linsolvercreator=DefaultLinSolverCreator(),
                   tol::Number = eps(real(T))*100,
                   λ::Number = zero(T),
                   v::Vector = randn(size(nep,1)),
                   target::Number = zero(T),
                   logger = 0,
                   inner_logger = 0) where {T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    @parse_logger_param!(logger)
    @parse_logger_param!(inner_logger)

    # Initial logical checks
    n = size(nep,1)
    if (maxit > n)
        error("maxit = ", maxit, " is larger than size of NEP = ", n,".")
    end
    if (projtype != :Galerkin) && projtype != :PetrovGalerkin
        error("Only accepted values of 'projtype' are :Galerkin and :PetrovGalerkin.")
    end
    if (projtype != :Galerkin) && isa(inner_solver_method,SGIterInnerSolver)
        error("Need to use 'projtype' :Galerkin in order to use SGITER as inner solver.")
    end

    # Allocations and preparations
    λ::T = T(λ)
    target::T = T(target)
    tol::real(T) = real(T)(tol)
    λ_vec::Vector{T} = Vector{T}(undef,neigs)
    u_vec::Matrix{T} = zeros(T,n,neigs)
    u::Vector{T} = Vector{T}(v)
    normalize!(u)
    conveig = 0

    # Initial check for convergence
    err = estimate_error(errmeasure,λ,u)
    if (err < tol) #Frist check, no other eiganvalues can be converged
        conveig += 1
        λ_vec[conveig] = λ
        u_vec[:,conveig] = u
    end
    if (conveig == neigs)
        return (λ_vec,u_vec)
    end

    # More allocations and preparations
    pk::Vector{T} = zeros(T,n)
    s_memory::Vector{T} = zeros(T,n)
    proj_nep = create_proj_NEP(nep,maxit,T)
    dummy_vector::Vector{T} = zeros(T,maxit+1)

    V_memory::Matrix{T} = zeros(T, size(nep,1), maxit+1)
    V_memory[:,1] = u
    if( projtype == :PetrovGalerkin ) # Petrov-Galerkin uses a left (test) and a right (trial) space
        W_memory = zeros(T, size(nep,1), maxit+1)
        W_memory[:,1] = compute_Mlincomb(nep,λ,u)
        normalize!(view(W_memory,:,1))
    else # Galerkin uses the same trial and test space
        W_memory = view(V_memory, :, :);
    end

    # Main loop
    for k = 1:maxit
        # Projection matrices
        V = view(V_memory, :, 1:k); W = view(W_memory, :, 1:k); # extact subarrays, memory-CPU efficient
        v = view(V_memory, :, k+1); w = view(W_memory, :, k+1); # next vector position
        s = view(s_memory,1:k)

        # Project and find the eigenvalue projected NEP
        expand_projectmatrices!(proj_nep, W, V)
        λv,sv = inner_solve(inner_solver_method, T, proj_nep,
                            j = conveig+1, # For SG-iter
                            λv = λ .* ones(T,conveig+1),
                            σ = target,
                            neigs = conveig+1,
                            inner_logger = inner_logger)
        λ,s = jd_eig_sorter(λv, sv, conveig+1, target)
        normalize!(s)

        # the approximate eigenvector
        u[:] = V*s

        # Check for convergence
        err = estimate_error(errmeasure,λ,u)
        push_iteration_info!(logger,k,λ=λ,err=err,continues=true)
        push_info!(logger," conveig=$conveig");

        if (err < tol) && (conveig == 0 ||
                all( abs.(λ .- λ_vec[1:conveig])./abs.(λ_vec[1:conveig]) .> sqrt(sqrt(eps(real(T)))) ) )
            conveig += 1
            λ_vec[conveig] = λ
            u_vec[:,conveig] = u
        end
        if (conveig == neigs)
            return (λ_vec,u_vec)
        end

        # Solve for basis extension using comment on top of page 367 of Betcke
        # and Voss, to avoid matrix access. Orthogonalization to u comes anyway
        # since u in V. OBS: Non-standard in JD-literature
        pk[:] = compute_Mlincomb(nep,λ,u,[one(T)],1)
        linsolver::LinSolver=create_linsolver(linsolvercreator,nep,λ)
        v[:] = lin_solve(linsolver, pk, tol=tol) # M(λ)\pk
        orthogonalize_and_normalize!(V, v, view(dummy_vector, 1:k), orthmethod)

        if( projtype == :PetrovGalerkin )
            w[:] = compute_Mlincomb(nep,λ,u);
            orthogonalize_and_normalize!(W, w, view(dummy_vector, 1:k), orthmethod)
        end
    end #End main loop

    msg="Number of iterations exceeded. maxit=$(maxit) and only $(conveig) eigenvalues converged out of $(neigs)."
    throw(NoConvergenceException(cat(λ_vec[1:conveig], λ, dims = 1), cat(u_vec[:,1:conveig], u, dims = 2), err, msg))
end

function jd_eig_sorter(λv::Vector{T}, V, N, target::T) where T <: Number
    NN = min(N, length(λv))
    c = sortperm(abs.(λv.-target))
    λ::T = λv[c[NN]]
    s::Vector{T} = V[:,c[NN]]
    return λ, s
end

"""
    jd_effenberger([eltype]], nep::ProjectableNEP; [maxit=100], [neigs=1], [inner_solver_method=DefaultInnerSolver()], [orthmethod=DGKS], [linsolvercreator=DefaultLinSolverCreator()], [tol=eps(real(T))*100], [λ=zero(T)], [v = rand(T,size(nep,1))], [target=zero(T)],  [logger=0], [inner_logger=0])
The function computes eigenvalues using the Jacobi-Davidson method, which is a projection method.
Repreated eigenvalues are avoided by using deflation, as presented in the reference by Effenberger.
The projected problems are solved using a solver spcified through the type `inner_solver_method`.
The logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.
For numerical stability the basis is kept orthogonal, and the method for orthogonalization is specified by `orthmethod`, see the package `IterativeSolvers.jl`.
The function tries to compute `neigs` number of eigenvalues, and throws a `NoConvergenceException` if it cannot.
The value `λ` and the vector `v` are initial guesses for an eigenpair. `linsolvercreator` is a function which specifies how the linear system is created and solved.
The `target` is the center around which eiganvalues are computed. For further specifications on the `deflation_mode`, see the function `deflate_eigpair`.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("dep0",100);
julia> λ,v=jd_effenberger(nep,maxit=30,v=ones(size(nep,1)),λ=0);
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
9.577305772525487e-15
```

# References
* C. Effenberger, Robust successive computation of eigenpairs for nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl. 34, 3 (2013), pp. 1231-1256.
See also
* T. Betcke and H. Voss, A Jacobi-Davidson-type projection method for nonlinear eigenvalue problems. Future Gener. Comput. Syst. 20, 3 (2004), pp. 363-372.
* H. Voss, A Jacobi–Davidson method for nonlinear eigenproblems. In: International Conference on Computational Science. Springer, Berlin, Heidelberg, 2004. pp. 34-41.
"""
jd_effenberger(nep::NEP;kwargs...) = jd_effenberger(ComplexF64,nep;kwargs...)
function jd_effenberger(::Type{T},
                        nep::ProjectableNEP;
                        maxit::Int = 100,
                        neigs::Int = 1,
                        inner_solver_method = DefaultInnerSolver(),
                        orthmethod::Type{T_orth} = IterativeSolvers.DGKS,
                        linsolvercreator=DefaultLinSolverCreator(),
                        tol::Number = eps(real(T))*100,
                        λ::Number = rand(T),
                        v::Vector = rand(T,size(nep,1)),
                        target::Number = zero(T),
                        deflation_mode = :Auto,
                        logger = 0,
                        inner_logger = 0) where {T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    @parse_logger_param!(logger)
    @parse_logger_param!(inner_logger)

    # Initial logical checks
    n = size(nep,1)
    if (maxit > n)
        error("maxit = ", maxit, " is larger than size of NEP = ", n,".")
    end
    if isa(inner_solver_method,SGIterInnerSolver)
        error("Inner solver 'SGIterInnerSolver' not accepted since deflated problem not min-max.")
    end

    # Allocations and preparations
    λ::T = T(λ)
    u::Vector{T} = Vector{T}(v)
    normalize!(u)
    λ_init::T = λ
    u_init::Vector{T} = u

    target::T = T(target)
    tol::real(T) = real(T)(tol)
    conveig = 0
    tot_nrof_its = 0
    V_memory_base::Matrix{T} = zeros(T, n+neigs, maxit+1)
    W_memory_base::Matrix{T} = zeros(T, n+neigs, maxit+1)

    # Initial check for convergence
    err = norm(compute_Mlincomb(nep,λ,u))
    if (err < tol) # Check if initial guess is good enough, otherwise compute initial invariant pair
        λ_init = T(rand())
        u_init = rand(n+1)
    else
        V_memory = view(V_memory_base, 1:n, tot_nrof_its+1:(maxit+1))
        W_memory = view(W_memory_base, 1:n, tot_nrof_its+1:(maxit+1))
        λ, u, tot_nrof_its, u_init, λ_init = dispatch_inner!(T, V_memory, W_memory,
                                                          nep, maxit, tot_nrof_its,
                                                          conveig, inner_solver_method, orthmethod,
                                                          linsolvercreator, tol, target, logger,
                                                          neigs, u_init, λ_init, inner_logger)
    end

    conveig += 1
    deflated_nep = deflate_eigpair(nep, λ, u, mode=deflation_mode)

    while true # Can only escape the loop on convergence (return) or too many iterations (error)
        # Check for fulfillment. If so, compute the eigenpairs and return
        if (conveig == neigs)
            (λ_vec, u_vec,) = get_deflated_eigpairs(deflated_nep)
            return λ_vec, u_vec
        end

        V_memory = view(V_memory_base, 1:(n+conveig), tot_nrof_its+1:(maxit+1))
        W_memory = view(W_memory_base, 1:(n+conveig), tot_nrof_its+1:(maxit+1))
        λ, u, tot_nrof_its, u_init, λ_init = dispatch_inner!(T, V_memory, W_memory,
                                                          deflated_nep, maxit, tot_nrof_its,
                                                          conveig, inner_solver_method, orthmethod,
                                                          linsolvercreator, tol, target, logger,
                                                          neigs, u_init, λ_init, inner_logger)
        conveig += 1 #OBS: minimality index = 1, hence only exapnd by one

        deflated_nep = deflate_eigpair(deflated_nep, λ, u)
    end
end

# Two dispatch functions depending on if it is a DeflatedNEP or a ProjectableNEP
# A bit hacky, but avoids corner cases with "empty deflation"
function dispatch_inner!(::Type{T},
                              V_memory::SubArray,
                              W_memory::SubArray,
                              target_nep::DeflatedNEP,
                              args...) where {T<:Number}
    return  jd_effenberger_inner!(T, V_memory, W_memory,
                                  target_nep, target_nep.V0, target_nep.S0, target_nep.orgnep,
                                  args...)
end
function dispatch_inner!(::Type{T},
                              V_memory::SubArray,
                              W_memory::SubArray,
                              target_nep::ProjectableNEP,
                              args...) where {T<:Number}
    return  jd_effenberger_inner!(T, V_memory, W_memory,
                                  target_nep, zeros(T,size(target_nep,1),0), zeros(T,0,0), target_nep,
                                  args...)
end

function jd_effenberger_inner!(::Type{T},
                              V_memory::SubArray,
                              W_memory::SubArray,
                              target_nep::Union{DeflatedNEP,ProjectableNEP}, # OBS: target_nep can either be a DeflatedNEP, or another NEP if we have no initial invariant pair
                              X::Matrix{T}, # This is explicitly allocated by the dispatch_inner!(...)
                              Λ::Matrix{T}, # This is explicitly allocated by the dispatch_inner!(...)
                              orgnep::NEP, # This is explicitly allocated by the dispatch_inner!(...)
                              maxit::Int,
                              nrof_its::Int,
                              conveig::Int,
                              inner_solver_method,
                              orthmethod::Type,
                              linsolvercreator,
                              tol::Number,
                              target::Number,
                              logger::Logger,
                              neigs::Int,
                              u::Vector{T},
                              λ::T,
                              inner_logger::Logger) where {T<:Number}
    # Allocations and preparations
    n = size(orgnep,1) # Size of original problem
    m = size(Λ,1) # Size of deflated subspace

    normalize!(u)

    λ_temp::T = λ
    newton_step::Vector{T} = Vector{T}(rand(n+m))
    pk::Vector{T} = zeros(T,n+m)
    s_memory::Vector{T} = zeros(T,maxit+1-nrof_its)
    proj_nep = create_proj_NEP(target_nep,maxit+1-nrof_its,T)

    dummy_vector::Vector{T} = zeros(T,maxit+1-nrof_its)

    V_memory[:,1] = u
    W_memory[:,1] = compute_Mlincomb(target_nep, λ, u)
    normalize!(view(W_memory,:,1))

    err = Inf
    for loop_counter = (nrof_its+1):maxit
        k = loop_counter - nrof_its # Which index are we doing on THIS specific level of deflation
        V = view(V_memory, :, 1:k); W = view(W_memory, :, 1:k); # extact subarrays, memory-CPU efficient
        v = view(V_memory, :, k+1); w = view(W_memory, :, k+1); # next vector position
        s = view(s_memory,1:k)

        # Project and solve the projected NEP
        expand_projectmatrices!(proj_nep, W, V)
        λv,sv = inner_solve(inner_solver_method, T, proj_nep,
                            tol = tol/10,
                            λv = λ .* ones(T,2),
                            σ = target,
                            neigs = 2,
                            inner_logger = inner_logger)
        λ_temp,s = jd_eig_sorter(λv, sv, 1, target) #Always closest to target, since deflated
        normalize!(s)

        # If inner solver converged to a solution to the projected problem use that.
        # Otherwise take the "newton step". Not implemented exactly as in Effenbergerm but similar
        if !isnan(λ_temp) && !any(isnan.(s[1:k])) && (norm(compute_Mlincomb(proj_nep, λ_temp, s[1:k])) < tol*50) # Has converged to an eigenvalue of the projected problem
            u[:] = V*s # The approximate eigenvector
            λ = λ_temp
        else
            u[:] = u[:] + newton_step
            normalize!(u)
        end

        # Compute residual and check for convergence
        rk::AbstractVector = compute_Mlincomb(target_nep, λ, u)
        err = norm(rk) #Error measure, see (9)/(3.2) in Effenberger # TODO: Can we use another error measure?
        push_iteration_info!(logger,loop_counter,λ=λ,err=err,continues=true);
        push_info!(logger," conveig=$conveig, subspace dim=$k");

        if (err < tol) #Frist check, no other eiganvalues can be converged
            info_str = "One eigenvalue converged."
            if conveig+1 < neigs
                info_str = string(info_str, " Deflating and restarting.")
            end
            push_info!(logger, info_str)

            # TODO: Here one can implement a continuation with the same basis as in Effenberger section 4.2.5
            # What is here is only a light kind adapted to an "unknown" inner solver
            λ2,s2 = jd_eig_sorter(λv, sv, 2, target)
            if (size(sv,2) > 1) && (abs(λ-λ2)/abs(λ) > sqrt(eps(real(T))))
                normalize!(s2)
                u2 = vcat(V*s2, zero(T))
            else
                λ2 = T(rand())
                u2 = Vector{T}(rand(n+m+1))
            end
            return (λ, u, loop_counter, u2, λ2)

        end

        # Extend the basis
        # The following analogy can be made if w is chosen as [y_k; v_k] in (26)/(4.6) in Effenberger
        # Solve for basis extension using comment on top of page 367 of Betcke
        # and Voss, to avoid matrix access. Orthogonalization to u comes anyway
        # since u in V. OBS: Non-standard in JD-literature
        pk[:] = compute_Mlincomb(target_nep, λ, u, [one(T)], 1)
        linsolver::LinSolver=create_linsolver(linsolvercreator,orgnep,λ)
        jd_inner_effenberger_linear_solver!(v, target_nep, λ, linsolver, pk, tol)
        newton_step = copy(v)
        orthogonalize_and_normalize!(V, v, view(dummy_vector, 1:k), orthmethod)
        w[:] = rk
        orthogonalize_and_normalize!(W, w, view(dummy_vector, 1:k), orthmethod)

    end

    msg="Number of iterations exceeded. maxit=$(maxit) and only $(conveig) eigenvalues converged out of $(neigs)."
    # Compute the eigenvalues we have and throw these.
    F = eigen(Λ)
    λ_vec = F.values
    uv = F.vectors
    u_vec = X * uv
    throw(NoConvergenceException([λ_vec; λ], [u_vec u[1:n]], err, msg))

end

function jd_inner_effenberger_linear_solver!(v, deflated_nep::DeflatedNEP, λ::T, linsolver, pk::Vector{T}, tol) where{T}
    # If it is a deflated NEP we solve with a Schur complement strategy such that
    # the user specified solve of M can be used.
    # (M, U; X^T, 0)(v1;v2) = (pk1;pk2)
    # OBS: Assume minimality index = 1
    # OBS: Forms the Schur complement. Assume that only a few eigenvalues are deflated

    X = deflated_nep.V0
    Λ = deflated_nep.S0
    orgnep = deflated_nep.orgnep
    n = size(orgnep,1)
    m = size(Λ,1)

    v1 = view(v, 1:n)
    v2 = view(v, (n+1):(n+m))
    pk1 = pk[1:n]
    pk2 = pk[(n+1):(n+m)]
    U = deflated_nep_compute_Q(deflated_nep, λ, 0)

    # Precompute some reused entities
    pk1tilde = lin_solve(linsolver, pk1, tol=tol) # pk1tilde = M^{-1}pk1
    Z::Matrix{T} = zeros(T,n,m)
    for i = 1:m
        Z[:,i] = lin_solve(linsolver, vec(U[:,i]), tol=tol) # Z = M^{-1}U
    end
    S = -X'*Z #Schur complement
    v2[:] = S\(pk2 - X'*pk1tilde)
    v1[:] = pk1tilde - Z*v2

    return v
end

function jd_inner_effenberger_linear_solver!(v, nep::NEP, λ::T, linsolver, pk::Vector{T}, tol) where{T}
    v[:] = lin_solve(linsolver, pk, tol=tol)
    return v
end


text:
# parser

function read_instance(MyFileName::String)
  path = "./Instances_ECMA/" * MyFileName
	# Si le fichier path existe
	if isfile(path) # exemple : "./Instances_ECMA/20_USA-road-d.BAY.gr"
		# L’ouvrir
		myFile = open(path)

    # Lire les premieres lignes
    n = parse(Int64, split(readline(myFile), " ")[3])
    s = parse(Int64, split(readline(myFile), " ")[3])
    t = parse(Int64, split(readline(myFile), " ")[3])
    S = parse(Int64, split(readline(myFile), " ")[3])
    d1 = parse(Int64, split(readline(myFile), " ")[3])
    d2 = parse(Int64, split(readline(myFile), " ")[3])

    p = Vector{Int64}(undef,0)
    line = split(readline(myFile), " ")
    number = line[3]
    number = chop(number, head=1,tail=1)
    number = parse(Int64, number)
    append!(p, number)
    popfirst!(line)
    popfirst!(line)
    popfirst!(line)
    for chunk in line
      chunk = chop(chunk)
      append!(p, parse(Int64, chunk))
    end

    ph = Vector{Int64}(undef,0)
    line = split(readline(myFile), " ")
    number = line[3]
    number = chop(number, head=1,tail=1)
    number = parse(Int64, number)
    append!(ph, number)
    popfirst!(line)
    popfirst!(line)
    popfirst!(line)
    for chunk in line
      chunk = chop(chunk)
      append!(ph, parse(Int64, chunk))
    end

    readline(myFile)
    # valeur 0 si il n'y a pas d'arc
    Mat_d = Array{Int64,2}(zeros(n,n))
    Mat_D = Array{Float64,2}(zeros(n,n))
		# Lire toutes les lignes d’un fichier
		data = readlines(myFile)
		for datum in data
      line = split(datum, " ")
      line[3] = chop(line[3])
      line[4] = chop(line[4])
      Mat_d[parse(Int64,line[1]),parse(Int64,line[2])] = parse(Int64,line[3])
      Mat_D[parse(Int64,line[1]),parse(Int64,line[2])] = parse(Float64,line[4])
		end

		return n, s, t, S, d1, d2, p, ph, Mat_d, Mat_D
	end
end

max_stars_repo_name: boumpteryx/ECMA

max_issues_repo_name: boumpteryx/ECMA

max_forks_repo_name: boumpteryx/ECMA


text:
"""
    grdview(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; kwargs...)

Reads a 2-D grid and produces a 3-D perspective plot by drawing a mesh, painting a
colored/grayshaded surface made up of polygons, or by scanline conversion of these polygons
to a raster image.

Full option list at [`grdview`]($(GMTdoc)grdview.html)

- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** | **drapefile** :: [Type => Str | GMTgrid | a Tuple with 3 GMTgrid types]

    Drape the image in drapefile on top of the relief provided by relief_file.
    ($(GMTdoc)grdview.html#g)
- **I** | **shade** | **intensity** | **intensfileintens** :: [Type => Str | GMTgrid]

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    ($(GMTdoc)grdview.html#i)
- **N** | **plane** :: [Type => Str | Int]

    Draws a plane at this z-level.
    ($(GMTdoc)grdview.html#n)
- $(GMT.opt_P)
- **Q** | **surftype** | **surf_type** :: [Type => Str | Int]

    Specify **m** for mesh plot, **s* for surface, **i** for image.
    ($(GMTdoc)grdview.html#q)
- **S** | **smooth** :: [Type => Number]

    Smooth the contours before plotting.
    ($(GMTdoc)grdview.html#s)
- **T** | **no_interp** :: [Type => Str]

    Plot image without any interpolation.
    ($(GMTdoc)grdview.html#t)
- **W** | **pens** :: [Type => Str]

    Draw contour, mesh or facade. Append pen attributes.
    ($(GMTdoc)grdview.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdview(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdview", cmd0, arg1)
	arg2 = nothing;	arg3 = nothing;	arg4 = nothing;	arg5 = nothing;

	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, = parse_BJR(d, "", "grdview", O, " -JX12c/0")
	cmd = parse_common_opts(d, cmd, [:UVXY :JZ :c :f :n :p :t :params], first)
	cmd = add_opt(cmd, 'S', d, [:S :smooth])
	if ((val = find_in_dict(d, [:N :plane])[1]) !== nothing)
		cmd *= " -N" * parse_arg_and_pen(val, "+g", false)
	end
	cmd = add_opt(cmd, 'Q', d, [:Q :surf :surftype],
				  (mesh=("m", add_opt_fill), waterfall=("mx", add_opt_fill), surface="_s",
				   surf="_s", img="i", image="i", nan_alpha="_c", monochrome="_+m"))
	cmd = add_opt(cmd, 'W', d, [:W :pens :pen], (contour=("c", add_opt_pen),
	              mesh=("m", add_opt_pen), facade=("f", add_opt_pen)) )
	cmd = add_opt(cmd, 'T', d, [:T :no_interp], (skip="_+s", outlines=("+o", add_opt_pen)) )

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted

	if (isa(arg1, Array{<:Number}))  arg1 = mat2grid(arg1)  end

	cmd, N_used, arg1, arg2, arg3 = get_cpt_set_R(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, arg3, "grdview")
	cmd, arg1, arg2, arg3, arg4 = common_shade(d, cmd, arg1, arg2, arg3, arg4, "grdview")

	if ((val = find_in_dict(d, [:G :drapefile])[1]) !== nothing)
		if (isa(val, String))				# Uff, simple. Either a file name or a -A type modifier
			cmd *= " -G" * val
		elseif (isa(val, GMTgrid))			# A single drape grid (arg1-3 may be used already)
			cmd, N_used = put_in_slot(cmd, val, 'G', [arg1, arg2, arg3, arg4])
			if     (N_used == 1)  arg1 = val
			elseif (N_used == 2)  arg2 = val
			elseif (N_used == 3)  arg3 = val
			elseif (N_used == 4)  arg4 = val
			end
		elseif (isa(val, Tuple) && length(val) == 3)
			cmd, N_used = put_in_slot(cmd, val[1], 'G', [arg1, arg2, arg3, arg4, arg5])
			cmd *= " -G -G"		# Because the above only set one -G and we need 3
			if     (N_used == 1)  arg1 = val[1];	arg2 = val[2];		arg3 = val[3]
			elseif (N_used == 2)  arg2 = val[1];	arg3 = val[2];		arg4 = val[3]
			elseif (N_used == 3)  arg3 = val[1];	arg4 = val[2];		arg5 = val[3]
			end
		else
			error("Wrong way of setting the drape (G) option.")
		end
	end

    return finish_PS_module(d, "grdview " * cmd, "", K, O, true, arg1, arg2, arg3, arg4, arg5)
end

# ---------------------------------------------------------------------------------------------------
grdview!(cmd0::String="", arg1=nothing; first=false, kw...) = grdview(cmd0, arg1; first=first, kw...)
grdview(arg1; first=true, kw...) = grdview("", arg1; first=first, kw...)
grdview!(arg1; first=false, kw...) = grdview("", arg1; first=first, kw...)


text:
using Gurobi
using Mosek

# train adversarial method using constraint generation
function train_adv_cg(X::Matrix, y::Vector, C::Real=1.0;
  perturb::Real=0.0, tol::Real=1e-6, psdtol::Real=1e-6, obj_reltol::Real=0.0,
  log::Real=0, n_thread::Int=0, solver::Symbol=:gurobi, verbose::Bool=true)

  n = length(y)
  # add one
  X1 = [ones(n) X]'   # transpose
  m = size(X1, 1)

  # number of class
  n_c = maximum(y)
  n_f = n_c * m   # number of features

  # parameters. init with zero
  w = zeros(n_f)

  alpha = zeros(n)

  # array of tuple
  constraints = Tuple{Integer, Vector}[]

  # prepare saved vars
  idmi = map(i -> idi(m, i), collect(1:n_c))

  # precompute xi dot xj
  K = [( i >= j ? dot(view(X1,:,i), view(X1,:,j)) : 0.0 )::Float64 for i=1:n, j=1:n]
  K = [( i >= j ? K[i,j] : K[j,i] )::Float64 for i=1:n, j=1:n]

  if solver == :gurobi
    # gurobi solver
    # gurobi environtment
    env = Gurobi.Env()
    # Method : 0=primal simplex, 1=dual simplex, 2=barrier  ; default for QP: barrier
    # Threads : default = 0 (use all threads)
    setparams!(env, PSDTol=psdtol, LogToConsole=log, Method=2, Threads=n_thread)
  elseif solver == :mosek
    # mosek environtment
    env = makeenv()
  end

  # params for Gurobi
  Q = zeros(0,0)
  nu = zeros(0)
  A = spzeros(n,0)        # sparse matrix
  b = ones(n) * C

  Q_prev = zeros(0,0)
  nu_prev = zeros(0)
  A_prev = spzeros(n,0)   # sparse matrix

  # additional for w
  L = zeros(n_f,0)
  L_prev = zeros(n_f,0)

  iter = 0
  dual_obj = 0.0
  dual_obj_prev = -Inf

  while true
    iter += 1

    if verbose
      println("Iteration : ", iter)
      tic();
    end

    # previous constraints
    const_prev = copy(constraints)

    ## add to constraints
    const_added = Tuple{Integer, Vector}[]

    # find constraint for each sample
    for i=1:n
      psis = psi_list(w, X1, y, i, n_c, idmi)
      psis_id, val = best_psis(psis)      # most violated constraints

      # current xi_i
      id_i = find(x -> x[1] == i, constraints)
      xi_i_list =  map(x -> x[2], constraints[id_i])
      max_xi_i = 0
      for j = 1:length(xi_i_list)
        a = calc_const(psis::Vector, xi_i_list[j])
        if a > max_xi_i
          max_xi_i = a
        end
      end

      if val > max_xi_i
        cs = (i, sort!(psis_id))
        if findfirst(constraints .== cs) == 0
          push!(constraints, cs)
          push!(const_added, cs)
        end
      end
    end

    # if no constraints added
    if length(const_added) == 0
      break
    end

    n_const = length(constraints)

    #### Start QP ###

    if verbose
      println(">> Start QP")
      toc();
      tic();
    end

    n_prev = length(const_prev)
    n_added = length(const_added)

    Q_aug = [ ( calc_dot((const_prev[i][1], const_prev[i][2], const_added[j][1], const_added[j][2]), K, y) )::Float64
          for i=1:n_prev, j=1:n_added]

    Q_aug_diag = [
          ( i >= j ? calc_dot((const_added[i][1], const_added[i][2], const_added[j][1], const_added[j][2]), K, y) : 0.0 )::Float64
          for i=1:n_added, j=1:n_added]
    Q_aug_diag = [ (i >= j ? Q_aug_diag[i,j] : Q_aug_diag[j,i] )::Float64 for i=1:n_added, j=1:n_added ]

    Q = [
          ( (i <= n_prev && j <= n_prev) ? Q_prev[i, j] : ( i <= n_prev ? Q_aug[i, j-n_prev] :
          ( j <= n_prev ? Q_aug[j, i-n_prev] : Q_aug_diag[i-n_prev, j-n_prev] ) ) )::Float64
          for i=1:n_const, j=1:n_const]

    nu_aug = [ ( calc_cconst(const_added[i][2]) )::Float64 for i=1:n_added ]
    nu = [ ( i <= n_prev ? nu_prev[i] : nu_aug[i-n_prev] )::Float64 for i=1:n_const]

    A_aug = spzeros(n, n_added)
    for j = 1:n_added
      A_aug[const_added[j][1], j] = 1.0
    end
    A = [A_prev A_aug]

    ## add perturbation
    for i=1:n_const
      Q[i,i] = Q[i,i] + perturb
    end

    if verbose
      toc();
      tic();
    end

    if solver == :gurobi

      if verbose println(">> Optim :: Gurobi") end

      ## init model
      model = gurobi_model(env,
                 sense = :minimize,
                 H = Q,
                 f = -nu,
                 A = A,
                 b = b,
                 lb = zeros(n_const)
                 )
      # Print the model to check correctness
      # print(model)

      # Solve with Gurobi
      Gurobi.optimize(model)

      if verbose
        toc();
        println("<< End QP")
      end

      dual_obj = -get_objval(model)
      # Solution
      if verbose println("Objective value: ", dual_obj) end

      # get alpha
      alpha = get_solution(model)

      if verbose println("n constraints = ", length(constraints)) end

      ### end QP ###

    elseif solver == :mosek

      if verbose println(">> Optim :: Mosek") end

      task = maketask(env)

      # set params
      putintparam(task, Mosek.MSK_IPAR_LOG, 1)
      putintparam(task, Mosek.MSK_IPAR_LOG_CHECK_CONVEXITY, 1)
      putdouparam(task, Mosek.MSK_DPAR_CHECK_CONVEXITY_REL_TOL, psdtol)

      # variables
      appendvars(task, n_const)
      # bound on var
      for i::Int32 = 1:n_const
        putbound(task, Mosek.MSK_ACC_VAR, i, Mosek.MSK_BK_RA, 0.0, C)
      end

      # objective
      putobjsense(task, Mosek.MSK_OBJECTIVE_SENSE_MINIMIZE)
      qi = zeros(Int32, (n_const * (n_const+1)) ÷ 2 )
      qj = zeros(Int32, (n_const * (n_const+1)) ÷ 2 )
      qv = zeros(Float64, (n_const * (n_const+1)) ÷ 2 )
      ix = 1
      for j::Int32 = 1:n_const
        for i::Int32 = j:n_const
          qi[ix] = i
          qj[ix] = j
          qv[ix] = Q[i,j]
          ix += 1
        end
      end
      putqobj(task, qi, qj, qv)

      putclist(task, collect(1:n_const), -nu)

      # constraints
      ## sparse array
      appendcons(task, n)
      for i::Int32 = 1:n_const
        id_nz = A[:,i].nzind
        putacol(task, i, id_nz, ones(length(id_nz)))
      end

      for i::Int32 = 1:n
        putbound(task, Mosek.MSK_ACC_CON, i, Mosek.MSK_BK_RA, 0.0, C)
      end

      if verbose
        toc(); tic();
      end

      Mosek.optimize(task)

      if verbose
        toc();
        println("<< End QP")
      end

      # Solution
      dual_obj, _ = getsolutioninfo(task, Mosek.MSK_SOL_ITR)
      if verbose println("Objective value: ", -dual_obj) end

      # get alpha
      alpha = getxx(task, Mosek.MSK_SOL_ITR)

      if verbose println("n constraints = ", length(constraints)) end
    end

    L_aug = zeros(n_f, n_added)
    for i=1:n_added
      L_aug[:, i] = calc_dconst(const_added[i], X1, y, n_c, idmi)
    end
    L = [L_prev L_aug]

    # recover w
    w = zeros(n_f)
    for i = 1:n_const
      w -= alpha[i] * L[:, i]
    end

    if obj_reltol > 0.0
      if (dual_obj - dual_obj_prev) / dual_obj_prev < obj_reltol && iter > 1
        if verbose println(">> Iteration STOPPED | Objective relative tolerance : ", obj_reltol) end
        break
      end
    end

    if verbose println() end

    Q_prev = Q
    nu_prev = nu
    A_prev = A
    L_prev = L
    dual_obj_prev = dual_obj
  end

  # finalizing losses
  gv_aug = zeros(n)
  gv_01 = zeros(n)
  l_adv = zeros(n)
  l_01 = zeros(n)
  for i=1:n
    psis = psi_list(w, X1, y, i, n_c, idmi)
    psis_id, val = best_psis(psis)      # most violated constraints
    n_ps = length(psis_id)

    gv_aug[i] = val

    # compute probs
    p_hat = zeros(n_c)
    p_check = zeros(n_c)
    for j=1:n_c
      if j in psis_id
        p_hat[j] = ( (n_ps-1.0)*psis[j] - sum(psis[psis_id[psis_id .!= j]]) + 1.0 ) / n_ps
        p_check[j] = 1.0 / n_ps
      else
        p_hat[j] = 0.0
        p_check[j] = 0.0
      end
    end

    C01 = 1 - eye(n_c)    # 01 loss matrix
    v = p_hat' * C01 * p_check   # the result is vector size 1 not a number
    gv_01[i] = v[1]     #

    ## training loss
    l_adv[i] = 1.0 - p_hat[y[i]]
    l_01[i] = 1.0 - round(Int, indmax(p_hat) == y[i])
  end

  game_value_01 = mean(gv_01)
  game_value_augmented = mean(gv_aug)

  # create model
  adv_model = MultiAdversarialModel(w, alpha, constraints, n_c, game_value_01, game_value_augmented, mean(l_adv), mean(l_01))

  return adv_model::MultiAdversarialModel
end

function predict_adv(model::MultiAdversarialModel, X_test::Matrix)

  w = model.w
  n_c = model.n_class
  n = size(X_test, 1)

  X1 = [ones(n) X_test]'   # transpose
  m = size(X1, 1)

  # prepare saved vars
  idmi = map(i -> idi(m, i), collect(1:n_c))

  prob = zeros(n, n_c)
  pred = zeros(n)
  for i=1:n
    psis = psi_list(w, X1, i, n_c, idmi)
    psis_id, val = best_psis(psis)      # most violated constraints
    n_ps = length(psis_id)

    for j=1:n_c
      if j in psis_id
        prob[i,j] = ( (n_ps-1.0)*psis[j] - sum(psis[psis_id[psis_id .!= j]]) + 1.0 ) / n_ps
      else
        prob[i,j] = 0.0
      end
    end

    pred[i] = indmax(psis)

  end

  return prob::Matrix{Float64}, pred::Vector{Float64}
end

function test_adv(model::MultiAdversarialModel, X_test::Matrix, y_test::Vector)
  n = size(X_test, 1)

  y_prob, y_pred = predict_adv(model, X_test)

  # calculate testing loss
  losses = zeros(n)
  losses01 = zeros(n)
  for i=1:n
    losses[i] = 1.0 - y_prob[i, y_test[i]]
    losses01[i] = 1.0 - round(Int, y_pred[i] == y_test[i])
  end

  loss = sum(losses) / n
  loss01 = sum(losses01) / n

  return loss::Float64, losses::Vector{Float64}, loss01::Float64, losses01::Vector{Float64},
   y_prob::Matrix{Float64}, y_pred::Vector{Float64}

end

max_stars_repo_name: rizalzaf/adversarial-multiclass

max_issues_repo_name: rizalzaf/adversarial-multiclass

max_forks_repo_name: rizalzaf/adversarial-multiclass


text:
module SchedulersTests

include("utils.jl")
include("test_basics.jl")
include("test_multi_scheduler.jl")
include("test_global_scheduler.jl")
include("test_benchmarks.jl")

end  # module SchedulersTests