text:
################################################################################
# GameObjective
################################################################################

mutable struct GameObjective
	probsize::ProblemSize  # Problem Size
	obj::Vector{Vector{Objective}} # objective of each player [number of players, number of obejctives]
	E::Vector{Vector{Objective}}   # quadratic expansion of the objective of each player [number of players, number of obejctives]
end

function GameObjective(Q::Vector{DQ}, R::Vector{DR}, xf::Vector{SVx}, uf::Vector{SVu},
	N::Int, model::AbstractGameModel) where {DQ,DR,SVx,SVu}
	T = eltype(xf[1])
	n = model.n
	m = model.m
	p = model.p
	probsize = ProblemSize(N,model)
	@assert length(Q) == length(R) == length(xf) == length(uf)

	obj = [Vector{Objective}(undef, 1) for i=1:p]
	E = [Vector{Objective}(undef, 1) for i=1:p]
	for i = 1:p
		Qi = Diagonal(SVector{n,T}(expand_vector(diag(Q[i]), model.pz[i], n)))
		Ri = Diagonal(SVector{m,T}(expand_vector(diag(R[i]), model.pu[i], m)))
		Rf = Diagonal(zeros(SVector{m,T}))
		xfi = SVector{n,T}(expand_vector(xf[i], model.pz[i], n))
		ufi = SVector{m,T}(expand_vector(uf[i], model.pu[i], m))
		cost = LQRCost(Qi,Ri,xfi,ufi,checks=false)
		cost_term = LQRCost(Qi,Rf,xfi,checks=false)
		obj[i][1] = Objective(cost, cost_term, N)
		E[i][1] = Objective([LQRCost(1e-10*ones(n,n)+I, 1e-10*zeros(m,m)+I, zeros(MVector{n,T})) for k=1:N])
	end
	return GameObjective(probsize,obj,E)
end

function expand_vector(v::AbstractVector{T}, inds::AbstractVector{Int}, n::Int) where {T}
	V = zeros(n)
	V[inds] = v
	return V
end

function cost_gradient!(game_obj::GameObjective, pdtraj::PrimalDualTraj)
	p = game_obj.probsize.p
	for i = 1:p
		# n_obj = length(game_obj.E[i])
		# for j = 1:n_obj
		# 	TrajectoryOptimization.cost_gradient!(game_obj.E[i][j], game_obj.obj[i][j], pdtraj.pr, true)
		# end
		cost_gradient!(game_obj, pdtraj, i)
	end
	return nothing
end

function cost_gradient!(game_obj::GameObjective, pdtraj::PrimalDualTraj, i::Int)
	n_obj = length(game_obj.E[i])
	for j = 1:n_obj
		TrajectoryOptimization.cost_gradient!(game_obj.E[i][j], game_obj.obj[i][j], pdtraj.pr, true)
	end
	return nothing
end

function cost_hessian!(game_obj::GameObjective, pdtraj::PrimalDualTraj)
	p = game_obj.probsize.p
	for i = 1:p
		# n_obj = length(game_obj.E[i])
		# for j = 1:n_obj
		# 	TrajectoryOptimization.cost_hessian!(game_obj.E[i][j], game_obj.obj[i][j], pdtraj.pr, true, true)
		# end
		cost_hessian!(game_obj, pdtraj, i)
	end
	return nothing
end

function cost_hessian!(game_obj::GameObjective, pdtraj::PrimalDualTraj, i::Int)
	n_obj = length(game_obj.E[i])
	for j = 1:n_obj
		TrajectoryOptimization.cost_hessian!(game_obj.E[i][j], game_obj.obj[i][j], pdtraj.pr, true, true)
	end
	return nothing
end

function add_collision_cost!(game_obj::GameObjective, radius::AbstractVector{T}, μ::AbstractVector{T}) where {T}
	N = game_obj.probsize.N
	n = game_obj.probsize.n
	m = game_obj.probsize.m
	p = game_obj.probsize.p
	px = game_obj.probsize.px
	@assert p == length(radius) == length(μ)
 	for i = 1:p
		for j ∈ setdiff(1:p,i)
			obj = Objective(CollisionCost{n,m,T,length(px[i])}(μ[i], radius[i], px[i], px[j]), N)
			E = Objective([LQRCost(1e-10*ones(n,n)+I, 1e-10*zeros(m,m)+I, zeros(MVector{n,T})) for k=1:N])
			push!(game_obj.obj[i], obj)
			push!(game_obj.E[i], E)
		end
	end
	return nothing
end

################################################################################
# Collision Cost
################################################################################

mutable struct CollisionCost{n,m,T,ni} <: TrajectoryOptimization.CostFunction
	μ::T
	r::T
	pxi::SVector{ni,Int}
	pxj::SVector{ni,Int}
    terminal::Bool

    function CollisionCost{n,m,T,ni}(μ::T, r::T, pxi::SVector{ni,Int},
		pxj::SVector{ni,Int}; terminal=false) where {n,m,T,ni}
        new{n,m,T,ni}(μ,r,pxi,pxj,terminal)
    end
end

function TrajectoryOptimization.stage_cost(cost::CollisionCost, x::AbstractVector, u::AbstractVector)
    J = TrajectoryOptimization.stage_cost(cost, x)
    return J
end

function TrajectoryOptimization.stage_cost(cost::CollisionCost, x::AbstractVector{T}) where T
	xi = x[cost.pxi]
	xj = x[cost.pxj]
	0.5 * cost.μ * max(0., cost.r - norm(xi-xj))^2
end

function TrajectoryOptimization.gradient!(E::TrajectoryOptimization.QuadraticCostFunction, cost::CollisionCost, x)
	n = length(x)
	ϵ = 1e-10
	ϵ_norm = ϵ*sqrt(n)
	xi = x[cost.pxi]
	xj = x[cost.pxj]
	Δ = (xi-xj)
	Δ_norm = norm(xi-xj)
	E.q .= 0.
	if max(0., cost.r - Δ_norm) > 0.
		g = cost.μ * (cost.r * (ϵ .+ Δ)/(ϵ_norm + Δ_norm) - Δ)
		E.q[cost.pxi] .= - g
	    E.q[cost.pxj] .=   g
	end
    return false
end

function TrajectoryOptimization.gradient!(E::TrajectoryOptimization.QuadraticCostFunction, cost::CollisionCost, x, u)
    TrajectoryOptimization.gradient!(E, cost, x)
    E.r .= 0.
    return false
end

function TrajectoryOptimization.hessian!(E::TrajectoryOptimization.QuadraticCostFunction, cost::CollisionCost, x)
	n = length(x)
	ϵ = 1e-10
	ϵ_norm = ϵ*sqrt(n)
	xi = x[cost.pxi]
	xj = x[cost.pxj]
	Δ = (xi-xj)
	Δ_norm = norm(xi-xj)
	E.Q .= 0.
	if max(0., cost.r - Δ_norm) > 0.
		E.Q[cost.pxi, cost.pxi] .=  cost.μ * (I - cost.r*I/Δ_norm + cost.r*(Δ*Δ')/(Δ_norm^3))
		E.Q[cost.pxi, cost.pxj] .= -cost.μ * (I - cost.r*I/Δ_norm + cost.r*(Δ*Δ')/(Δ_norm^3))
		E.Q[cost.pxj, cost.pxi] .= -cost.μ * (I - cost.r*I/Δ_norm + cost.r*(Δ*Δ')/(Δ_norm^3))
		E.Q[cost.pxj, cost.pxj] .=  cost.μ * (I - cost.r*I/Δ_norm + cost.r*(Δ*Δ')/(Δ_norm^3))
	end
    return true
end

function TrajectoryOptimization.hessian!(E::TrajectoryOptimization.QuadraticCostFunction, cost::CollisionCost, x, u)
    TrajectoryOptimization.hessian!(E, cost, x)
    E.R .= 0
    return true
end

import Base.copy
function Base.copy(c::CollisionCost{n,m,T,ni}) where {n,m,T,ni}
    CollisionCost{n,m,T,ni}(copy(c.μ), copy(c.r), copy(c.pxi), copy(c.pxj), terminal=c.terminal)
end

function TrajectoryOptimization.state_dim(cost::CollisionCost{n,m,T,ni}) where {n,m,T,ni}
	return n
end

function TrajectoryOptimization.control_dim(cost::CollisionCost{n,m,T,ni}) where {n,m,T,ni}
	return m
end


text:
export Performative, Message, GenericMessage, MessageClass, AbstractMessageClass, ParameterReq, ParameterRsp, set!

# global variables
const _messageclasses = Dict{String,DataType}()

"An action represented by a message."
module Performative
  const REQUEST = "REQUEST"
  const AGREE = "AGREE"
  const REFUSE = "REFUSE"
  const FAILURE = "FAILURE"
  const INFORM = "INFORM"
  const CONFIRM = "CONFIRM"
  const DISCONFIRM = "DISCONFIRM"
  const QUERY_IF = "QUERY_IF"
  const NOT_UNDERSTOOD = "NOT_UNDERSTOOD"
  const CFP = "CFP"
  const PROPOSE = "PROPOSE"
  const CANCEL = "CANCEL"
end

"Base class for messages transmitted by one agent to another."
abstract type Message <: AbstractDict{Symbol,Any} end

"""
    mtype = MessageClass(context, clazz[, superclass[, performative]])

Create a message class from a fully qualified class name. If a performative is not
specified, it is guessed based on the class name. For class names ending with "Req",
the performative is assumed to be REQUEST, and for all other messages, INFORM.

# Examples

```julia-repl
julia> MyShellExecReq = MessageClass(@__MODULE__, "org.arl.fjage.shell.ShellExecReq");
julia> req = MyShellExecReq(cmd="ps")
ShellExecReq: REQUEST [cmd:"ps"]
```
"""
function MessageClass(context, clazz::String, superclass=nothing, performative=nothing)
  sname = replace(string(clazz), "." => "_")
  tname = sname
  if performative === nothing
    performative = match(r"Req$",string(clazz))==nothing ? Performative.INFORM : Performative.REQUEST
  end
  if superclass === nothing
    superclass = "$(@__MODULE__).Message"
  else
    scname = string(superclass)
    ndx = findlast(isequal('.'), scname)
    if ndx !== nothing
      scname = scname[ndx+1:end]
    end
    if scname == tname
      tname = "$(scname)_"
    end
  end
  expr = Expr(:toplevel)
  expr.args = [Meta.parse("""
      struct $(tname) <: $(superclass)
        clazz::String
        data::Dict{String,Any}
        $(tname)(c::String, d::Dict{String,Any}) = new(c, d)
      end
    """),
    Meta.parse("""
      function $(sname)(; kwargs...)
        dict = Dict{String,Any}(
          "msgID" => string($(@__MODULE__).uuid4()),
          "perf" => "$(performative)"
        )
        for k in keys(kwargs)
          dict[string(k)] = kwargs[k]
        end
        return $(tname)("$(clazz)", dict)
      end
    """),
    Meta.parse("""
      $(@__MODULE__)._messageclasses["$(clazz)"] = $(tname)
    """)]
  if sname != tname
    push!(expr.args, Meta.parse("$(tname)(; kwargs...) = $(sname)(; kwargs...)"))
  end
  return context.eval(expr)
end

function AbstractMessageClass(context, clazz::String, performative=nothing)
  sname = replace(string(clazz), "." => "_")
  expr = Expr(:toplevel)
  expr.args = [Meta.parse("abstract type $sname <: $(@__MODULE__).Message end"), Meta.parse("$sname")]
  rv = context.eval(expr)
  MessageClass(context, clazz, rv, performative)
  return rv
end

"""
    registermessages()
    registermessages(messageclasses)

Register message classes with Fjage. Usually message classes are automatically registered on
creation with `MessageClass()`. However, when developing packages, if `MessageClass()` is used
at the module level, the types may be precompiled and the code to register the classes may not
get executed at runtime. In such cases, you may need to explicitly call `registermessages()`
in the `__init()__` funciton for the module.
"""
function registermessages(msg=subtypes(Fjage.Message))
  for t ∈ msg
    endswith(string(t), '_') && continue
    s = t().__clazz__
    Fjage._messageclasses[s] = t
    registermessages(subtypes(t))
  end
end

function _messageclass_lookup(clazz::String)
  haskey(_messageclasses, clazz) && return _messageclasses[clazz]
  Message
end

# helper function to see if a message matches a filter
function _matches(filt, msg)
  (msg === nothing || filt === nothing) && return true
  if typeof(filt) == DataType
    return typeof(msg) <: filt
  elseif typeof(filt) <: Message
    return msg.inReplyTo == filt.msgID
  elseif typeof(filt) <: Function
    return filt(msg)
  end
  false
end

# adds notation message.field
function Base.getproperty(s::Message, p::Symbol)
  if p == :__clazz__
    return getfield(s, :clazz)
  elseif p == :__data__
    return getfield(s, :data)
  else
    p1 = string(p)
    if p1 == "performative"
      p1 = "perf"
    elseif p1 == "messageID"
      p1 = "msgID"
    end
    v = getfield(s, :data)
    if !haskey(v, p1)
      return nothing
    end
    v = v[p1]
    return v
  end
end

# immutable dictionary interface for Messages

function Base.get(s::Message, p::Symbol, default)
  v = getproperty(s, p)
  v === nothing && return default
  v
end

Base.getindex(s::Message, p::Symbol) = getproperty(s, p)
Base.keys(s::Message) = Symbol.(keys(getfield(s, :data)))
Base.values(s::Message) = values(getfield(s, :data))
Base.eltype(s::Message) = Pair{Symbol,Any}
Base.length(s::Message) = length(getfield(s, :data))

function Base.iterate(s::Message)
  it = iterate(getfield(s, :data))
  it === nothing && return nothing
  (Symbol(it[1][1]) => it[1][2], it[2])
end

function Base.iterate(s::Message, state)
  it = iterate(getfield(s, :data), state)
  it === nothing && return nothing
  (Symbol(it[1][1]) => it[1][2], it[2])
end

# adds notation message.field
function Base.setproperty!(s::Message, p::Symbol, v)
  (p == :__clazz__ || p == :__data__) && throw(ArgumentError("read-only property cannot be set"))
  p1 = string(p)
  if p1 == "performative"
    p1 = "perf"
  elseif p1 == "messageID"
    p1 = "msgID"
  end
  getfield(s, :data)[p1] = v
  nothing
end

# pretty prints arrays without type names
function _repr(x)
  x = repr(x)
  m = match(r"[A-Za-z0-9]+(\[.+\])", x)
  m !== nothing && (x = m[1])
  x
end

# pretty printing of messages
function Base.show(io::IO, msg::Message)
  ndx = findlast(".", msg.__clazz__)
  s = ndx === nothing ? msg.__clazz__ : msg.__clazz__[ndx[1]+1:end]
  p = ""
  data_suffix = ""
  signal_suffix = ""
  suffix = ""
  data = msg.__data__
  for k in keys(data)
    x = data[k]
    if k == "perf"
      s *= ": " * x
    elseif k == "data"
      if typeof(x) <: Array
        data_suffix *= "($(length(x)) bytes)"
      else
        p *= " $k:" * _repr(data[k])
      end
    elseif k == "signal"
      if typeof(x) <: Array
        signal_suffix *= "($(length(x)) samples)"
      else
        p *= " $k:" * _repr(data[k])
      end
    elseif k != "sender" && k != "recipient" && k != "msgID" && k != "inReplyTo" && k != "sentAt"
      if typeof(x) <: Number || typeof(x) == String || typeof(x) <: Array || typeof(x) == Bool
        p *= " $k:" * _repr(x)
      else
        suffix = "..."
      end
    end
  end
  length(suffix) > 0 && (p *= " " * suffix)
  length(signal_suffix) > 0 && (p *= " " * signal_suffix)
  length(data_suffix) > 0 && (p *= " " * data_suffix)
  p = strip(p)
  length(p) > 0 && (s *= " [$p]")
  if msg.__clazz__ == "org.arl.fjage.GenericMessage"
    m = match(r"^GenericMessage: (.*)$", s)
    m === nothing || (s = m[1])
  end
  print(io, s)
end

"Generic message type that can carry arbitrary name-value pairs as data."
GenericMessage = MessageClass(@__MODULE__, "org.arl.fjage.GenericMessage")

"Parameter request message."
ParameterReq = MessageClass(@__MODULE__, "org.arl.fjage.param.ParameterReq")

"Parameter response message."
ParameterRsp = MessageClass(@__MODULE__, "org.arl.fjage.param.ParameterRsp")

"""
    msg = Message([perf])
    msg = Message(inreplyto[, perf])

Create a message with just a performative (`perf`) and no data. If the performative
is not specified, it defaults to INFORM. If the inreplyto is specified, the message
`inReplyTo` and `recipient` fields are set accordingly.
"""
Message(perf::String=Performative.INFORM) = _Message(perf=perf)
Message(inreplyto::Message, perf::String=Performative.INFORM) = _Message(perf=perf, inReplyTo=inreplyto.msgID, recipient=inreplyto.sender)

# message base class
_Message = MessageClass(@__MODULE__, "org.arl.fjage.Message")

# convenience methods and pretty printing for parameters

function org_arl_fjage_param_ParameterReq(vals...; index=-1)
  req = ParameterReq(index=index)
  qlist = Pair{String,Any}[]
  for v ∈ vals
    if v isa String
      push!(qlist, Pair{String,Any}(v, nothing))
    elseif v isa Symbol
      push!(qlist, Pair{String,Any}(string(v), nothing))
    elseif v isa Pair
      push!(qlist, Pair{String,Any}(string(v[1]), v[2]))
    end
  end
  if length(qlist) > 0
    q = popfirst!(qlist)
    req.param = q[1]
    req.value = q[2]
    if length(qlist) > 0
      req.requests = Dict{String,Any}[]
      for q ∈ qlist
        push!(req.requests, Dict{String,Any}("param" => q[1], "value" => q[2]))
      end
    end
  end
  req
end

"""
    get!(p::ParameterReq, param)

Request parameter `param` to be fetched.

# Examples

```julia-repl
julia> p = ParameterReq(index=1)
ParameterReq[index=1]
julia> get!(p, "modulation")
ParameterReq[index=1 modulation=?]
julia> get!(p, "fec")
ParameterReq[index=1 modulation=? ...]
```
"""
function Base.get!(p::ParameterReq, param)
  if p.param === nothing
    p.param = param
  else
    p.requests === nothing && (p.requests = Dict{String,Any}[])
    push!(p.requests, Dict{String,Any}("param" => param))
  end
  p
end

"""
    set!(p::ParameterReq, param, value)

Request parameter `param` to be set to `value`.

# Examples

```julia-repl
julia> p = ParameterReq(index=1)
ParameterReq[index=1]
julia> set!(p, "modulation", "ofdm")
ParameterReq[index=1 modulation=ofdm]
julia> set!(p, "fec", 1)
ParameterReq[index=1 modulation=ofdm ...]
```
"""
function set!(p::ParameterReq, param, value)
  if p.param === nothing
    p.param = param
    p.value = value
  else
    p.requests === nothing && (p.requests = Dict{String,Any}[])
    push!(p.requests, Dict{String,Any}("param" => param, "value" => value))
  end
  p
end

"""
    get(p::ParameterRsp, param)

Extract parameter `param` from a parameter response message.
"""
function Base.get(p::ParameterRsp, key)
  skey = string(key)
  dskey = "." * skey
  (p.param == skey || endswith(p.param, dskey)) && return p.value
  vals = p.values
  if vals !== nothing
    for (k, v) ∈ vals
      (k == skey || endswith(k, dskey)) && return v
    end
  end
  nothing
end

function Base.show(io::IO, p::ParameterReq)
  print(io, "ParameterReq[")
  p.index !== nothing && p.index ≥ 0 && print(io, "index=", p.index, ' ')
  p.param !== nothing && print(io, p.param, '=', (p.value === nothing ? "?" : string(p.value)))
  p.requests === nothing || print(io, " ...")
  print(io, ']')
end

function Base.show(io::IO, p::ParameterRsp)
  print(io, "ParameterRsp[")
  p.index !== nothing && p.index ≥ 0 && print(io, "index=", p.index, ' ')
  p.param !== nothing && print(io, p.param, '=', p.value)
  p.values === nothing || print(io, " ...")
  print(io, ']')
end

function Base.println(io::IO, p::ParameterRsp)
  plist = Pair{String,Any}[]
  if p.param !== nothing
    x = p.param
    occursin(".", x) || (x = "." * x)
    push!(plist, x => p.value)
    vs = p.values
    if vs !== nothing
      for v ∈ vs
        x = v[1]
        occursin(".", x) || (x = "." * x)
        push!(plist, x => v[2])
      end
    end
  end
  sort!(plist; by=(x -> x[1]))
  let n = findfirst(x -> x[1] == ".title", plist)
    n === nothing || println(io, "« ", plist[n][2], " »\n")
  end
  let n = findfirst(x -> x[1] == ".description", plist)
    n === nothing || plist[n][2] == "" || println(io, plist[n][2], "\n")
  end
  prefix = ""
  ro = p.readonly === nothing ? String[] : p.readonly
  for (k, v) ∈ plist
    k === ".type" && continue
    k === ".title" && continue
    k === ".description" && continue
    ks = split(k, '.')
    cprefix = join(ks[1:end-1], '.')
    if cprefix != prefix
      prefix != "" && println(io)
      prefix = cprefix
      println(io, '[', cprefix, ']')
    end
    println(io, "  ", ks[end], k ∈ ro ? " ⤇ " : " = ", v)
  end
end

ext: jl


text:
##############################################################################################################
##############################################################################################################
########################################         DATA          ###############################################
##############################################################################################################
##############################################################################################################
export Gene, Protein, GeneExpression, Methylation,
	RNASeq, RNASeqCall, ExomeSeq, RPPA, CNV,
	Drug, Pathway, DataView, CellLine, Outcome,
	Experiment, ViewType, KeyType, NAViewType, ValueViewType

""" An abstract type for samples, can be cell line or human tissue (latter not yet implemented) """
abstract type Sample end

""" Generic view type, i.e. a way to describe a cell line """
abstract type ViewType end

abstract type NAViewType <: ViewType end

abstract type ValueViewType <: ViewType end

abstract type SingleViewType <: ValueViewType end

abstract type VectorViewType <: ValueViewType end

"""
A gene. Mostly represented by its HGNC ID in the motivating Dreamchallenge data set.
Preferred ID is Entrez ID, Ensembl-ID can be stored as well.
As an additional data if available, can hold information on being a cancer gene.
"""
mutable struct Gene
	# is still mutable here, since we are continuously collecting information, maybe make immutable in the future
	entrez_id::Int64 # <- primary key preferred
	hgnc_id::String
	ensembl_id::String
	cancer_gene::Bool
	Gene() = new(-1,"","",false)
	function Gene(id::String;id_type::String="hgnc_id", is_cancer_gene::Bool=false)
		if id_type == "ensembl_id"
	 		return new(-1, "", id, is_cancer_gene)
	 	end
	 	# hgnc_id is standard type
	 	new(-1, id, "", is_cancer_gene)
	end
	Gene(entrez_id::Int64, hgnc_id::String, ensembl_id::String; is_cancer_gene::Bool=false) = new(entrez_id, hgnc_id, ensembl_id, is_cancer_gene)
	Gene(entrez_id::Int64, is_cancer_gene::Bool=false) = new(entrez_id, "", "", is_cancer_gene)
end

"""
A protein. If available, paralogs of the protein that are affected by the same antibody are encoded in the hgnc_id.
The field antibody_validated indicates whether cross-reactivity in the assay has been checked
and the antibodies behave linearly in the dilution series.
"""
struct Protein
	hgnc_id::String
	antibody_validated::Bool
end

import Base.isless
function isless(x::Gene, y::Gene)
	if x.entrez_id != -1 || y.entrez_id != -1
		return isless(x.entrez_id, y.entrez_id)
	end
	if x.ensembl_id != "" || y.ensembl_id != ""
		return isless(x.ensembl_id, y.ensembl_id)
	end
	if x.hgnc_id != "" || y.hgnc_id != ""
		return isless(x.hgnc_id, y.hgnc_id)
	end
	# if isless(x.entrez_id, y.entrez_id) return true
	# elseif x.entrez_id == y.entrez_id && (isless(x.ensembl_id, y.ensembl_id) || isless(x.hgnc_id, y.hgnc_id)) return true
	# else return false end
end
isless(x::Protein, y::Protein) = isless(x.hgnc_id, y.hgnc_id)

""" The types allowed as keys in data views """
const KeyType = Union{Gene, Protein}

""" A gene expression."""
mutable struct GeneExpression <: SingleViewType
	expression_value::Float64
	normalized_value::Float64
	GeneExpression(expression_value::Float64) = new(expression_value, expression_value)
end

"""
Illumina beta values indicating methylation.
CGct1 gives the number of CG dinucleotides in the probe sequence. The
U and M probes differ only at C's in CpG dinucleotides so there is a
good deal of cross-hyb between the two probes when the CG count is
low, and beta values shrink toward .5.  Probes are typically filtered
out with fewer than 3 CpGs as these values are unreliable.
Particularly, these probes will be less accurate for discrete M/U
calls.

Cct1 gives the number of off-CpG cytosines found in the probe. These
probes guarantee the specificity of the assay for successfully
bi-sulfite converted DNA, and probes with too little cytosine may
misinterpret unconverted DNA as methylated DNA, increasing beta
values. Probes are often filtered out that have fewer than 3
cytosines.
"""
mutable struct Methylation <: SingleViewType
	cgct1::Int64
	cct1::Int64
	illumina_beta_value::Float64
	is_methylated::Bool
	illumina_id::String
	normalized_value::Float64
	function Methylation(beta_value::Float64, cgct1::Int64, cct1::Int64, methylated_threshold::Float64=.2; illumina_id::String="")
		new(cgct1, cct1, beta_value, beta_value≥methylated_threshold, illumina_id, beta_value)
	end
end

"""
Whole transcriptome shotgun sequencing (RNA-seq) data.
Expression value is a fpkm value.
expression_status is 1 if the Ensembl gene model was detected above
the background noise level.
"""
mutable struct RNASeq <: SingleViewType
	expression_value::Float64
	normalized_value::Float64
	RNASeq(expression_value::Float64) = new(expression_value, expression_value)
end

"""
Whole transcriptome shotgun sequencing (RNA-seq) data.
expression_status is 1 if the Ensembl gene model was detected above
the background noise level.
"""
mutable struct RNASeqCall <: SingleViewType
	expression_status::Bool
	normalized_value::Bool
	RNASeqCall(expression_status::Bool) = new(expression_status, expression_status)
end

"""
Whole exome sequencing data. This is a vector view type because we may have multiple measurements for a single Gene key.

The different parts are:

num_cosmic - Number of samples in COSMIC that have mutation at this position.
variant_effect - Variant effect (silent, missense, nonsense, frame-shifiting Indel, etc.)
protein_change - Protein change
nucleotid_change - Change in DNAseq
variant_confidence - Confidence in variant (normalized in [0,1])
norm_zygosity - Zygosity (hom, het) of matched normal sample (blank)
norm_reference_count - # Reference alleles at position in matched normal (blank)
norm_variant_count - # Alternate alleles at position in matched normal (blank)
tumor_zygosity - Zygosity (hom, het) of Cell Line
tumor_reference_count - # Reference alleles at position in cell line
tumor_variant_count - # Alternate alleles at position in cell line
reference_mismatch_avg - METRIC: Average number of other mismatching bases in reads with reference base
variant_mismatch_avg - METRIC: Average number of other mismatching bases in reads with variant
reference_mismatch_sum - METRIC: Base quality sum of mismatching bases in reads with reference base
variant_mismatch_sum - METRIC: Base quality sum of mismatching bases in reads with variant
reference_dist3effective_avg - METRIC: Average normalized distance of reference bases from 3' end of their respective reads
variant_dist3effective_avg - METRIC: Average normalized distance of variant bases from 3' end of their respective reads
details - Various other information collected at this position.
"""
mutable struct ExomeSeq <: VectorViewType
	num_cosmic::Int64
	variant_effect::String
	protein_change::String
	nucleotid_change::String
	variant_confidence::Float64
	norm_zygosity::String
	norm_reference_count::Float64
	norm_variant_count::Float64
	tumor_zygosity::String
	tumor_reference_count::Float64
	tumor_variant_count::Float64
	reference_mismatch_avg::Float64
	variant_mismatch_avg::Float64
	reference_mismatch_sum::Float64
	variant_mismatch_sum::Float64
	reference_dist3effective_avg::Float64
	variant_dist3effective_avg::Float64
	details::String
	normalized_value::Float64
	function ExomeSeq(protein_change::AbstractString; reference_mismatch_sum::Float64=0., reference_mismatch_avg::Float64=0.,
		reference_dist3effective_avg::Float64=0., variant_mismatch_sum::Float64=0., variant_mismatch_avg::Float64=0.,
		variant_dist3effective_avg::Float64=0., num_cosmic::Int64=0, variant_effect::AbstractString="", nucleotid_change::AbstractString="",
		variant_confidence::Float64=1., norm_zygosity::AbstractString="", norm_reference_count::Float64=0., norm_variant_count::Float64=0.,
		tumor_zygosity::AbstractString="", tumor_reference_count::Float64=0., tumor_variant_count::Float64=0., details::AbstractString="")
		data_summary = reference_mismatch_avg + variant_mismatch_avg
						+ reference_mismatch_sum + variant_mismatch_sum
						+ reference_dist3effective_avg + variant_dist3effective_avg
		new(num_cosmic, variant_effect, protein_change, nucleotid_change, variant_confidence,
			norm_zygosity, norm_reference_count, norm_variant_count, tumor_zygosity, tumor_reference_count, tumor_variant_count,
			reference_mismatch_avg, variant_mismatch_avg, reference_mismatch_sum, variant_mismatch_sum,
			reference_dist3effective_avg, variant_dist3effective_avg, details,
			data_summary)
	end
end

""" Protein abundance value. """
mutable struct RPPA <: SingleViewType
	protein_abundance::Float64
	normalized_value::Float64
	RPPA(protein_abundance::Float64) = new(protein_abundance, protein_abundance)
end

""" The quantified value for DNA copy number analysis (aggregated across whole genes) """
mutable struct CNV <: SingleViewType
	gene_level_cnv::Float64
	normalized_value::Float64
	CNV(gene_level_cnv::Float64) = new(gene_level_cnv, gene_level_cnv)
end

"""
A drug, identified by a name.
Additional information may include knowledge about affected genes to
refine the number of features for predictions, this is useful for incorporating pathway information.
One useful approach could be to look at the generic pathways and select only those
that are in a known pathway.
More specific would be to select only genes of the pathway that is affected by a drug.

The chemical structure as a string can be used to compute chemical drug similarity
to additionally aid inclusion of correlations between drugs.
"""
mutable struct Drug
	id::String
	affected_genes::Vector{Gene}
	chemical_structure::String
	Drug(id::String; affected_genes::Vector{Gene}=Vector{Gene}(), chemical_structure::String="") = new(id, affected_genes, chemical_structure)
end

"""
Pathway information.
Pathway information is subsumed into a set of genes that affect the signaling pathway.
"""
struct Pathway{K<:KeyType}
	id::String
	name::String
	genes::Vector{K}
	Pathway{K}(genes::Vector{K}; name::String="", id::String="") where K = new(id, name, genes)
end

"""
A data view, typed with a identifier and data type.
The identifier is mostly of type Gene, for RPPA data it is Protein.
In general, similarity matrices for samples are computed on the basis of
these identifiers. In other words, for each data view, a cell line is represented
as a vector of genes (or proteins) with associated values. The identifiers determine
the ordering of this vector.
"""
struct DataView{K <: KeyType, V <: ViewType}
	cell_line_id::String
	key_type::Type{K}
	view_type::Type{V}
	used_keys::SortedSet{K}
	common_keys::Dict{String, Vector{K}}
	measurements::Union{Dict{K, V}, MultiDict{K, V}}
	DataView{K,V}(cell_line_id::String) where {K,V <: SingleViewType} = new(cell_line_id, K, V, SortedSet{K}(), Dict{String, Vector{K}}(), Dict{K,V}())
	DataView{K,V}(cell_line_id::String) where {K,V <: VectorViewType} = new(cell_line_id, K, V, SortedSet{K}(), Dict{String, Vector{K}}(), MultiDict{K,V}())
	DataView{K,V}(cell_line_id::String) where {K, V <: NAViewType} = new(cell_line_id, K, V, SortedSet{K}(), Dict{String, Vector{K}}(), Dict{K,V}())
end

"""
A cell line comprised with its ID, information on the cancer type and
dictionary of views, i.e. different data collected for this cell line.
"""
mutable struct CellLine <: Sample
	id::String
	cancer_type::String
	views::Dict{Type{<:ViewType}, DataView{<:KeyType, <:ViewType}}
	in_test_set::Bool
	CellLine(id::String, cancer_type::String="unknown"; in_test_set::Bool=false) = new(id, cancer_type, Dict{Type{<:ViewType}, DataView{<:KeyType, <:ViewType}}(), in_test_set)
end

"""
An experimental outcome.
An outcome should be associated with a drug.
The dictionary of outcome values associates the results to cell lines.
The reason for this is, that we believe the outcome across different
cell lines is correlated according to the molecular structure of the cell line tissue.
In other words, we model the outcome on available cell lines jointly.
"""
mutable struct Outcome
	drug_id::String
	outcome_values::OrderedDict{CellLine, Float64}
	outcome_type::String
	normalized_outcome_values::OrderedDict{CellLine, Float64}
	outcome_mean::Float64
	outcome_std::Float64
	Outcome(drug_id::String, outcome_type::String) = new(drug_id, OrderedDict{CellLine, Float64}(), outcome_type, OrderedDict{CellLine, Float64}(), .0, 1.)
end

# mutable struct Prediction
# 	predicted_outcome::Dict{Drug, Float64}
# 	predicted_rank::Dict{Drug, Int64}
# 	Prediction() = new(Dict{Drug, Float64}(), Dict{Drug,Int64}())
# end

"""
An experiment comprised of measurements (where each is associated
with a cell line) and results. This data structure holds all data that is provided
by the user, i.e. including pathway information, different data views, training and test responses etc.
"""
mutable struct Experiment
	results::OrderedDict{Drug, Outcome}
	# test_results::OrderedDict{Drug, Outcome}
	cell_lines::OrderedDict{String, CellLine}
	drugs::Dict{String, Drug}
	genes::Dict{Int64, Gene}
	genes_by_hgnc::Dict{String, Gene}
	genes_by_ensembl::Dict{String, Gene}
	proteins::Dict{String, Protein}
	views::OrderedSet{Type{<:ViewType}}
	statistics::Dict{Type{<:ViewType}, OrderedDict{KeyType, Tuple{Float64, Float64}}}
	pathway_information::OrderedDict{String, Pathway}
	is_normalized::Bool
	internal_id::String
	function Experiment(id::String)
		new(
			OrderedDict{Drug, Outcome}(), # results
			# OrderedDict{Drug, Outcome}(), # test results
			OrderedDict{String, CellLine}(), # cell_lines
			Dict{String, Drug}(), # drugs
			Dict{Int64, Gene}(), # genes
			Dict{String, Gene}(), # genes by hgnc
			Dict{String, Gene}(), # genes by ensembl
			Dict{String, Protein}(), # proteins
			OrderedSet{Type{<:ViewType}}(), # views
			Dict{Type{<:ViewType}, OrderedDict{KeyType, Tuple{Float64, Float64}}}(), # statistics
			OrderedDict{String, Pathway}(), #pathways
			false, # is normalized
			id
			)
	end
end


text:
#  Generates synthetic data.
module Synthetic
export gen_image!, gen_images!

using Celeste
using Celeste.Model
using Celeste.SensitiveFloats

import WCS
import Distributions
using ForwardDiff
using StaticArrays

function write_star_nmgy!(img::Image, ce::CatalogEntry)
    box = Model.box_around_point(img.wcs, ce.pos, 25)
    p = Model.ImagePatch(img, box)
    Model.write_star_nmgy!(ce.pos, ce.star_fluxes[img.b], p, img.pixels)
end

function write_galaxy_nmgy!(img::Image, ce::CatalogEntry)
    patches = Model.get_sky_patches([img], [ce], radius_override_pix=25.0)
    Model.write_galaxy_nmgy!(ce.pos, ce.gal_fluxes[img.b],
      ce.gal_frac_dev, ce.gal_axis_ratio, ce.gal_angle, ce.gal_radius_px,
      img.psf, patches, img.pixels)
end

function gen_image!(img::Image, n_bodies::Vector{CatalogEntry}; expectation=false)
    H, W = size(img.pixels)
    for w in 1:W, h in 1:H
        img.pixels[h, w] = img.sky[h, w]  # in nmgy
    end

    for body in n_bodies
        body.is_star ? write_star_nmgy!(img, body) : write_galaxy_nmgy!(img, body)
    end

    img.pixels .*= img.nelec_per_nmgy

    if !expectation
        for w in 1:W, h in 1:H
            img.pixels[h, w] = rand(Distributions.Poisson(img.pixels[h, w]))
        end
    end
end

"""
Generate a synthetic images based on a vector of catalog entries using
identity world coordinates.
"""
function gen_images!(images::Vector{<:Image}, n_bodies::Vector{CatalogEntry}; expectation=false)
    for img in images
        gen_image!(img, n_bodies; expectation=expectation)
    end
end

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

const pp = Model.load_prior()

function sample_fluxes(i::Int, r_s)
    k_s = rand(Distributions.Categorical(pp.k[i]))
    c_s = rand(Distributions.MvNormal(pp.color[i][:, k_s], pp.color[i][:, :, k_s]))

    l_s = Vector{Float64}(5)
    l_s[3] = r_s
    l_s[4] = l_s[3] * exp(c_s[3])
    l_s[5] = l_s[4] * exp(c_s[4])
    l_s[2] = l_s[3] / exp(c_s[2])
    l_s[1] = l_s[2] / exp(c_s[1])
    l_s
end

function synthetic_body(ce::CatalogEntry)
    ce2 = deepcopy(ce)
#    ce2.is_star = rand(Distributions.Bernoulli(pp.is_star[1]))
    ce2.star_fluxes[:] = sample_fluxes(1, ce.star_fluxes[3])
    ce2.gal_fluxes[:] = sample_fluxes(2, ce.gal_fluxes[3])
    ce2
end

function synthetic_bodies(n_bodies::Vector{CatalogEntry})
    CatalogEntry[synthetic_body(ce) for ce in n_bodies]
end

end


text:
using LineSearches
using Base.Test

lsfunctions = [hz_linesearch!, interpolating_linesearch!,
               mt_linesearch!, backtracking_linesearch!]
lsalphas = [0.5,0.5,0.49995,0.9]

f(x) = vecdot(x,x)

function g!(x, out)
    out[:] = 2x
end

df = LineSearches.DifferentiableFunction(f,g!)

for (i, linesearch!) in enumerate(lsfunctions)
    println("Testing $(string(linesearch!))")
    x = [-1., -1.]
    xtmp = copy(x)
    grtmp = similar(x)
    phi0 = df.fg!(x, grtmp)
    p = -grtmp
    dphi0 = dot(p, grtmp)

    lsr = LineSearchResults(eltype(x))
    push!(lsr, 0.0, phi0, dphi0)

    alpha = 1.0
    mayterminate = false

    alpha, f_update, g_update = linesearch!(df, x, p, xtmp, grtmp, lsr, alpha, mayterminate)
    #xnew = x + alpha*p

    @test_approx_eq alpha lsalphas[i]
end


text:
using Test, FinancialMonteCarlo;
@show "HestonModel"
S0 = 100.0;
K = 100.0;
r = [0.02; 0.02];
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
sigma_zero = 0.2;
kappa = 0.01;
theta = 0.03;
lambda = 0.01;
rho_ = 0.0;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8;

rfCurve = ZeroRate(r, T);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianData1 = AsianFloatingStrikeOption(T)
AsianData2 = AsianFixedStrikeOption(T, K)
Model = HestonProcess(sigma, sigma_zero, lambda, kappa, rho_, theta, Underlying(S0, d));

@show FwdPrice = pricer(Model, rfCurve, mc, FwdData);

@test abs(FwdPrice - 98.72567723404445) < toll

@show FwdPrice = pricer(Model, rfCurve, mc1, FwdData);
tollanti = 0.8
@test abs(FwdPrice - 98.72567723404445) < tollanti


text:
function generate(obj::ConvexHullGenerator, arg0::Collection)
    return jcall(obj, "generate", ConvexHull, (Collection,), arg0)
end


text:
using OptimizationDynamics
const iLQR = OptimizationDynamics.IterativeLQR
const Rotations = OptimizationDynamics.Rotations
using LinearAlgebra
using Random

# ## visualization
vis = Visualizer()
render(vis);

# ## mode
MODE = :projection
MODE = :nominal

# ## thrust max
u_max = 12.5

# ## state-space model
h = 0.05
T = 61
info = RocketInfo(rocket, u_max, h,
    eval(r_rocket_func), eval(rz_rocket_func), eval(rθ_rocket_func),
    eval(r_proj_func), eval(rz_proj_func), eval(rθ_proj_func))

nx = rocket.nq
nu = rocket.nu

# ## iLQR model
if MODE == :projection
    ilqr_dyn = iLQR.Dynamics((d, x, u, w) -> f_rocket_proj(d, info, x, u, w),
                        (dx, x, u, w) -> fx_rocket_proj(dx, info, x, u, w),
                        (du, x, u, w) -> fu_rocket_proj(du, info, x, u, w),
                        nx, nx, nu)
else
    ilqr_dyn = iLQR.Dynamics((d, x, u, w) -> f_rocket(d, info, x, u, w),
                    (dx, x, u, w) -> fx_rocket(dx, info, x, u, w),
                    (du, x, u, w) -> fu_rocket(du, info, x, u, w),
                    nx, nx, nu)
end

model = [ilqr_dyn for t = 1:T-1];

# ## initial conditions
x1 = zeros(rocket.nq)
x1[1] = 2.5
x1[2] = 2.5
x1[3] = 10.0
mrp = Rotations.MRP(Rotations.RotZ(0.25 * π) * Rotations.RotY(-0.5 * π))
x1[4:6] = [mrp.x; mrp.y; mrp.z]
x1[9] = -1.0

xT = zeros(rocket.nq)
xT[3] = rocket.length
mrpT = Rotations.MRP(Rotations.RotZ(0.25 * π) * Rotations.RotY(0.0))
xT[4:6] = [mrpT.x; mrpT.y; mrpT.z]

# ## objective
function objt(x, u, w)
	J = 0.0

	J += 0.5 * transpose(x - xT) * Diagonal(h * [1.0e-1 * ones(3); 1.0e-5 * ones(3); 1.0e-1 * ones(3); 1000.0 * ones(3)]) * (x - xT)
	J += 0.5 * transpose(u) * Diagonal(h * [1000.0; 1000.0; 100.0]) * u

	return J
end

function objT(x, u, w)
	J = 0.0

	J += 0.5 * transpose(x - xT) * Diagonal(h * 1000.0 * ones(rocket.nq)) * (x - xT)

	return J
end

ct = iLQR.Cost(objt, nx, nu)
cT = iLQR.Cost(objT, nx, 0)
obj = [[ct for t = 1:T-1]..., cT];

# ## constraints
x_con = [-0.5; 0.5]
y_con = [-0.75; 0.75]

function stage_con(x, u, w)
    if MODE == :projection
        [
         rocket.length - x[3];
        ]
    else
        [
         -1.0 - u[1];
         u[1] - 1.0;
         -1.0 - u[2];
         u[2] - 1.0;
         0.0 - u[3]; # control limit (lower)
         u[3] - u_max; # control limit (upper)
         rocket.length - x[3];
        ]
    end
end

function terminal_con(x, u, w)
    [
     x_con[1] - x[1];
     x[1] - x_con[2];
     y_con[1] - x[2];
     x[2] - y_con[2];
     (x - xT)[collect([3, 4, 5, 6, 7, 8, 9, 10, 11, 12])]
    ]
end

cont = iLQR.Constraint(stage_con, nx, nu, idx_ineq=collect(1:(1 + (MODE == :projection ? 0 : 6))))
conT = iLQR.Constraint(terminal_con, nx, 0, idx_ineq=collect(1:4))
cons = [[cont for t = 1:T-1]..., conT];

# ## rollout
Random.seed!(1)
ū = [1.0e-3 * randn(nu) for t = 1:T-1]
x̄ = iLQR.rollout(model, x1, ū)
visualize!(vis, rocket, x̄, Δt=h);

# ## solver
solver = iLQR.solver(model, obj, cons,
    opts=iLQR.Options(
        linesearch = :armijo,
        α_min=1.0e-5,
        obj_tol=1.0e-3,
        grad_tol=1.0e-3,
        max_iter=100,
        max_al_iter=15,
        con_tol=0.005,
        ρ_init=1.0,
        ρ_scale=10.0,
        verbose=false))
iLQR.initialize_controls!(solver, ū)
iLQR.initialize_states!(solver, x̄);

# ## solve
iLQR.reset!(solver.s_data)
@time iLQR.solve!(solver);

@show iLQR.eval_obj(solver.m_data.obj.costs, solver.m_data.x, solver.m_data.u, solver.m_data.w)
@show solver.s_data.iter[1]
@show norm(terminal_con(solver.m_data.x[T], zeros(0), zeros(0))[4 .+ (1:10)], Inf)
@show solver.s_data.obj[1] # augmented Lagrangian cost

# ## solution
x_sol, u_sol = iLQR.get_trajectory(solver)
visualize!(vis, rocket, x_sol, Δt=h);
open(vis)
# ## test thrust cone constraint
all([norm(u[1:2]) <= u[3] for u in u_sol])

# ## benchmark
solver.options.verbose = false
@benchmark iLQR.solve!($solver, x̄, ū) setup=(x̄=deepcopy(x̄), ū=deepcopy(ū));


text:
module Bases

import IterTools: groupby, imap
import LinearAlgebra: dot
import ..Utils: Interval

export Basis, Basis1D, domain, deriv, order,
       BSplineBasis,
       NURBSBasis

# Abstract Basis types
# ========================================================================

abstract type Basis end
abstract type Basis1D <: Basis end

mutable struct BasisFunction1D{B<:Basis1D}
    basis::B
    index::Int
    deriv::Int

    function BasisFunction1D{B}(basis, index, deriv) where {B<:Basis1D}
        if !(1 <= index <= length(basis)) throw(BoundsError()) end
        if !(0 <= deriv <= nderivs(basis)) throw(ArgumentError("Differentiation order not supported")) end
        new(basis, index, deriv)
    end

    BasisFunction1D(basis, index) = BasisFunction1D(basis, index, 0)
end

deriv(b::BasisFunction1D) = typeof(b)(b.basis, b.index, b.deriv+1)
deriv(b::BasisFunction1D, order) = typeof(b)(b.basis, b.index, b.deriv+order)

function (b::Basis1D)(pt::T) where {T<:Real}
    rng = supported(b, pt)
    (dropdims(evaluate_raw(b, [pt], b.deriv, rng), dims=2), rng)
end

function (b::BasisFunction1D)(pt::T) where {T<:Real}
    rng = supported(b.basis, pt)
    if b.index ∉ rng return 0.0 end
    evaluate_raw(b.basis, [pt], b.deriv, rng)[1 + b.index - rng.start, 1]
end

function (b::Basis1D)(pts::Vector{T}) where {T<:Real}
    tp = Tuple{Vector{T}, UnitRange{Int}}
    res = Array{tp}(undef, length(pts))

    j = 0
    for (subpts, rng) in supported(b, pts)
        out = evaluate_raw(b, subpts, b.deriv, rng)
        for i in 1:length(subpts)
            res[j+=1] = (out[:,i], rng)
        end
    end

    res
end

function (b::BasisFunction1D)(pts::Vector{T}) where {T<:Real}
    res = zeros(Float64, length(pts))

    i = 1
    for (subpts, rng) in supported(b.basis, pts)
        i += length(subpts)
        if b.index ∉ rng continue end

        out = evaluate_raw(b.basis, subpts, b.deriv, rng)
        res[i-length(subpts):i-1] = out[findall(in(b.index), rng), :]
    end

    res
end

function (b::Basis1D)(pt::S, coeffs::Vector{T}) where {S<:Real, T<:Real}
    (vals, idxs) = b(pt)
    dot(vals, coeffs[idxs])
end

function (b::Basis1D)(pt::S, coeffs::Array{T, 2}) where {S<:Real, T<:Real}
    (vals, idxs) = b(pt)
    vals' * coeffs[idxs,:]
end

function (b::Basis1D)(pt::S, coeffs::Array{T, 3}) where {S<:Real, T<:Real}
    (vals, idxs) = b(pt)
    #res = Array{Float64}(undef, size(coeffs, 2), size(coeffs, 3))
    #for i in 1:size(coeffs, 2)
    #    res[i, :] = vals' * coeffs[idxs, i, :]
    #end
    copy(reshape(
        vals' * reshape(
            coeffs[idxs, :, :], size(idxs, 1),
            size(coeffs, 2) * size(coeffs, 3)
        ),
        size(coeffs, 2),
        size(coeffs, 3)
    ))
end

(b::Basis1D)(pts::Vector{S}, coeffs::Vector{T}) where {S<:Real, T<:Real} =
    Float64[dot(vals, coeffs[idxs]) for (vals, idxs) in b(pts)]

function (b::Basis1D)(pts::Vector{S}, coeffs::Array{T, 2}) where {S<:Real, T<:Real}
    res = zeros(Float64, length(pts), size(coeffs, 2))

    for (i, (vals, idxs)) in enumerate(b(pts))
        res[i,:] = vals' * coeffs[idxs,:]
    end

    res
end

function (b::Basis1D)(pts::Vector{S}, coeffs::Array{T, 3}) where {S<:Real, T<:Real}
    res = zeros(Float64, length(pts), size(coeffs, 2), size(coeffs, 3))

    for (i, (vals, idxs)) in enumerate(b(pts)), j in 1:size(coeffs, 2)
        res[i, j, :] = vals' * coeffs[idxs, j, :]
    end

    res
end

# Include specific bases
# ========================================================================

include("BSplineBasis.jl")
include("NURBSBasis.jl")

end  # module Bases


text:
"""
    getcondvects_n_k(n::Int, k=1, val_vector=0:n-1)
    getcondvects_n_k(n::Vector{Int}, k::Int=length(n), val_vector::Union{Tuple, Vector}=missing)

Produces condition vectors. When n=2 and it produces a k-input truth table.

# Examples
julia> getcondvects_n_k(2, 3)
8×3 Matrix{Any}:
 0  0  0
 0  0  1
 0  1  0
 0  1  1
 1  0  0
 1  0  1
 1  1  0
 1  1  1
julia> getcondvects_n_k(3, 3, ['x', 'y', 'z'])
27×3 Matrix{Any}:
 'x'  'x'  'x'
 'x'  'x'  'y'
 'x'  'x'  'z'
 'x'  'y'  'x'
 'x'  'y'  'y'
 'x'  'y'  'z'
 'x'  'z'  'x'
 ⋮
 'z'  'y'  'x'
 'z'  'y'  'y'
 'z'  'y'  'z'
 'z'  'z'  'x'
 'z'  'z'  'y'
 'z'  'z'  'z'
julia> getcondvects_n_k([1,2,3], 3, [["x"], ['a','b'], [1, 2, 3]])
6×3 Matrix{Any}:
 "x"  'a'  1
 "x"  'a'  2
 "x"  'a'  3
 "x"  'b'  1
 "x"  'b'  2
 "x"  'b'  3

"""

function getcondvects_n_k(n::Int, k=1, val_vector=0:n-1)
	rows = n^k
	condvects=Matrix{Any}(missing, rows, k);
    for i = 1:k
        @views condvects[:, end+1-i] = reshape(repeat(permutedims(val_vector),
		Int(rows/n^(k-i+1)), Int(rows/n^i)), rows, 1);
    end
	return condvects
end

function getcondvects_n_k(n::Vector{Int}, k::Int=length(n),
	val_vector::Union{Tuple, Vector}=missing)
	@assert k == length(n) "n should have k entries"
	rows = prod(n)
	condvects=Matrix{Any}(missing, rows, k);
	if ismissing(val_vector)
		for i = 1:k
			@views condvects[:, i] = reshape(repeat(permutedims(1:n[i]),
			Int(rows/prod(n[1:i])), Int(rows/prod(n[i:end]))), rows, 1);
		end
	else
		@assert prod(length.(val_vector) .== n) "Length of each entry in val_vector should equal "
		for i = 1:k
	        @views condvects[:, i] = reshape(repeat(permutedims(val_vector[i]),
			Int(rows/prod(n[1:i])), Int(rows/prod(n[i:end]))), rows, 1);
		end
	end
	return condvects
end


text:
export run_spider,
       process_source!,
       process_html

process_source!(P,source::AbstractString,
                fileinfo::FileInfo;args...) = source

process_html(P,html::AbstractString,
             ::FileInfo; args...) = html

function run_spider(plugins...;
                    args...)
  sdir = get_arg(args,:source_dir)
  odir = get_arg(args,:output_dir)
  clear_output_dir = get_arg(args,:clear_output_dir,false)
  header_filename = get_arg(args,:header_file,"")
  footer_filename = get_arg(args,:footer_file,"")

  mdparser = CommonMark.Parser()

  if !isempty(header_filename)
    header_file = open(header_filename) do file read(file,String) end
  else
    header_file = Nothing
  end

  if !isempty(footer_filename)
    footer_file = open(footer_filename) do file read(file,String) end
  else
    footer_file = Nothing
  end

  if clear_output_dir
    run(`rm -fr $odir`)
  end
  run(`mkdir -p $odir`)

  for (root,dirs,files) in walkdir(sdir)
    offset = findfirst("/",root)
    folderstring = ""
    if !isnothing(offset)
      folderstring = root[offset[1]:end]
    end
    currs = sdir * folderstring
    curro = odir * folderstring

    for d in dirs
      run(`mkdir -p $curro/$d`)
    end

    for f in files
      (f[1]=='.') && continue
      if occursin(".",f)
        (basename,extension) = split(f,".")
      else
        basename = f
        extension = ""
      end
      sfname = currs*"/"*f
      if extension == "md"
        ofname = curro*"/"*basename*".html"
        mdstring = read(sfname,String)

        fileinfo = FileInfo()
        fileinfo["filename"] = f
        fileinfo["basename"] = basename
        fileinfo["extension"] = extension
        fileinfo["current_source_dir"] = currs
        fileinfo["current_output_dir"] = curro
        fileinfo["source_filename"] = sfname
        fileinfo["output_filename"] = ofname
        fileinfo["folderstring"] = folderstring
        fileinfo["source_dir"] = sdir
        fileinfo["output_dir"] = odir

        for P in plugins
          mdstring = process_source!(P,mdstring,fileinfo;args...)
        end

        #
        # Using custom parser
        #
        #open("_tmp_file.md","w") do tf
        #  print(tf,mdstring)
        #end
        #sp_md_command = split(md_parser)
        #html = read(`$sp_md_command _tmp_file.md`,String)

        #
        # Using CommonMark.jl
        #
        ast = mdparser(mdstring)
        htmlstr = CommonMark.html(ast)

        for P in plugins
          htmlstr = process_html(P,htmlstr,fileinfo;args...)
        end

        open(ofname,"w") do of
          if !isnothing(header_file)
            print(of,header_file)
          end

          print(of,htmlstr)

          if !isnothing(footer_file)
            print(of,footer_file)
          end
        end

        #
        # If using custom parser
        #
        #run(`rm -f _tmp_file.md`)

      else
        run(`cp $sfname $(curro*"/"*f)`)
      end
    end
  end
end


text:
using DSP
include("../notebook/DSPFilters.jl") #Defines DSPFilters.* & BodePlots.* utilities.
import InspectDR
using NumericIO
using Interact
using Blink

#=TODO
 - Make fl & fh sliders relative to fmax?
 - Disable unused parameters (depending on filter implementation)
 - Look for off by 1 errors (esp. for FFT vector lengths, etc).
 - Other details look not quite right (@ notches, phase response, ...)
=#

#Convenience aliases:
SI(x) = formatted(x, :SI, ndigits=3)
_throt(x) = throttle(0.1, x) #Default throttling for this example.

#When MIME"image/svg+xml" is disabled, Jupyter eventually requests PNG inline graphics:
#(Not certain if PNG is faster than SVG)
InspectDR.defaults.rendersvg = false

ltype_opt = OrderedDict("Series RLC" => :series, "Shunt RLC" => :shunt)

opt_ftype = ["Highpass", "Lowpass", "Bandstop", "Bandpass"]
opt_fimpl = ["Butterworth", "Elliptic", "ChebyshevPass", "ChebyshevStop"]

#Create base Observable-s that control plot output:
ftype = Observable(String("Lowpass")) #Filter type
fimpl = Observable(String("Elliptic")) #Filter implementation
fmax = Observable(Float64(500)) #Maximum plot frequency
forder = Observable(Int64(6)) #Filter order
fl = Observable(Float64(300)) #Low frequency point
fh = Observable(Float64(400)) #High frequency point
rpass = Observable(Float64(5)) #Pass-band ripple
rstop = Observable(Float64(40)) #Stop-band ripple

dispinline = Observable(Bool(true)) #Display inline?

#Construct widgets used to control plots:
txt_fmax = spinbox(value=fmax, label="max freq")
tgl_ftype = togglebuttons(opt_ftype, value=ftype)
tgl_fimpl = dropdown(opt_fimpl, value=fimpl)
sld_order = slider(2:1:30, value=forder, label="Filter order")
sld_rpass = slider(0.1:.1:20, value=rpass, label="Pass-band ripple (dB)")
sld_rstop = slider(5:1:200, value=rstop, label="Stop-band ripple (dB)")
sld_fl = slider(10:10:1e3, value=fl, label="Low frequency (Hz)")
sld_fh = slider(10:10:1e3, value=fh, label="High frequency (Hz)")

chk_inline = checkbox(dispinline, label="Display inline")
btn_creategtkplot = button("Create GTK plot")

#Build widget lists to control plot:
widgetlistA = [tgl_ftype, tgl_fimpl, sld_order, chk_inline, btn_creategtkplot]
widgetlistB = [sld_rpass, sld_rstop, sld_fl, sld_fh]

#Widget-controlled plot:
pobj = DSPFilters.newplot() #Create initial plot object
#Calculations are relatively expensive... throttle all inputs:
updater_plot = map(_throt(fmax), _throt(ftype), _throt(fimpl), _throt(forder), _throt(fl), _throt(fh), _throt(rpass), _throt(rstop)) do _fmax, _ftype, _fimpl, _forder, _fl, _fh, _rpass, _rstop
    _filttypeid = Symbol(lowercase(_ftype))
    _filtimplid = Symbol(lowercase(_fimpl))
    return DSPFilters.update(pobj, _fmax, _filttypeid, _filtimplid, _forder, _fl, _fh, _rpass, _rstop)
end

#Conditionally display inline plot with throttling to reduce traffic:
inline_plot = map(updater_plot, dispinline) do _plot, _inline
	if _inline
		return _plot
	else
		return nothing
	end
end

#Provide means to spawn new InspectDR-GTK plot windows:
creator_gtkplot = map(observe(btn_creategtkplot)) do _s
	#Display plot in Gtk GUI:
	gtkgui = display(InspectDR.GtkDisplay(), pobj) #Use same plot object

	#Refresh GUI when plot changes:
	updater_gtkgui = map(updater_plot) do _up
		InspectDR.refresh(gtkgui)
	end
end

#Display main window:
wnd = Window()
body!(wnd, vbox()) #first call to body! does not always work.
body!(wnd, vbox(hbox(vbox(widgetlistA...), vbox(widgetlistB...)), vbox(inline_plot)))

:DEMO_READY


text:
#####
##### A simple minmax player to be used as a baseline
#####

"""
A simple implementation of the minmax tree search algorithm, to be used as
a baseline against AlphaZero. Heuristic board values are provided by the
[`GameInterface.heuristic_value`](@ref) function.
"""
module MinMax

import ..GI, ..GameInterface, ..AbstractPlayer, ..think

amplify(r) = iszero(r) ? r : Inf * sign(r)

# Return the value of the current state for the player playing
function value(player, game, depth)
  if GI.game_terminated(game)
    return 0.
  elseif depth == 0
    return GI.heuristic_value(game)
  else
    qs = [qvalue(player, game, a, depth) for a in GI.available_actions(game)]
    return maximum(qs)
  end
end

function qvalue(player, game, action, depth)
  @assert !GI.game_terminated(game)
  next = copy(game)
  GI.play!(next, action)
  wr = GI.white_reward(next)
  r = GI.white_playing(game) ? wr : -wr
  if player.amplify_rewards
    r = amplify(r)
  end
  nextv = value(player, next, depth - 1)
  if GI.white_playing(game) != GI.white_playing(next)
    nextv = -nextv
  end
  return r + player.gamma * nextv
end

function minmax(player, game, actions, depth)
  return argmax([qvalue(player, game, a, player.depth) for a in actions])
end

"""
    MinMax.Player{Game} <: AbstractPlayer{Game}

A stochastic minmax player, to be used as a baseline.

    MinMax.Player{Game}(;depth, amplify_rewards, τ=0.)

The minmax player explores the game tree exhaustively at depth `depth`
to build an estimate of the Q-value of each available action. Then, it
chooses an action as follows:

- If there are winning moves (with value `Inf`), one of them is picked
  uniformly at random.
- If all moves are losing (with value `-Inf`), one of them is picked
  uniformly at random.

Otherwise,

- If the temperature `τ` is zero, a move is picked uniformly among those
  with maximal Q-value (there is usually only one choice).
- If the temperature `τ` is nonzero, the probability of choosing
  action ``a`` is proportional to ``e^{\\frac{q_a}{Cτ}}`` where ``q_a`` is the
  Q value of action ``a`` and ``C`` is the maximum absolute value of all
  finite Q values, making the decision invariant to rescaling of
  [`GameInterface.heuristic_value`](@ref).

If the `amplify_rewards` option is set to true, every received positive reward
is converted to ``∞`` and every negative reward is converted to ``-∞``.
"""
struct Player{G} <: AbstractPlayer{G}
  depth :: Int
  amplify_rewards :: Bool
  τ :: Float64
  gamma :: Float64
  function Player{G}(;depth, amplify_rewards, τ=0., γ=1.) where G
    return new{G}(depth, amplify_rewards, τ, γ)
  end
end

function think(p::Player, game)
  actions = GI.available_actions(game)
  n = length(actions)
  qs = [qvalue(p, game, a, p.depth) for a in actions]
  winning = findall(==(Inf), qs)
  if isempty(winning)
    notlosing = findall(>(-Inf), qs)
    best = argmax(qs)
    if isempty(notlosing)
      π = ones(n)
    elseif iszero(p.τ)
      π = zeros(n)
      all_best = findall(==(qs[best]), qs)
      π[all_best] .= 1.
    else
      qmax = qs[best]
      @assert qmax > -Inf
      C = maximum(abs(qs[a]) for a in notlosing) + eps()
      π = exp.((qs .- qmax) ./ C)
      π .^= (1 / p.τ)
    end
  else
    π = zeros(n)
    π[winning] .= 1.
  end
  π ./= sum(π)
  return actions, π, nothing
end

end

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
function DiffEqBase.__solve(prob::DiffEqBase.AbstractODEProblem,
                            alg::AbstractODEFilter, args...; kwargs...)
    @debug "Called solve with" args kwargs
    integrator = DiffEqBase.__init(prob, alg, args...; kwargs...)
    sol = DiffEqBase.solve!(integrator)
    return sol
end
DiffEqBase.__init(prob::DiffEqBase.AbstractODEProblem, alg::EKF0, args...; kwargs...) =
    DiffEqBase.__init(prob, ODEFilter(), args...; method=:ekf0, kwargs...)
DiffEqBase.__init(prob::DiffEqBase.AbstractODEProblem, alg::EKF1, args...; kwargs...) =
    DiffEqBase.__init(prob, ODEFilter(), args...; method=:ekf1, kwargs...)

function DiffEqBase.__init(prob::DiffEqBase.AbstractODEProblem, alg::ODEFilter;

                           method=:ekf1,
                           prior=:ibm,
                           q=1,
                           smooth=false,
                           initialize_derivatives=true,

                           bigfloat=false,

                           steprule=:standard,
                           dt=eltype(prob.tspan)(0),
                           abstol=1e-6, reltol=1e-3,
                           gamma=9//10,
                           qmin=2//10, qmax=10,
                           dtmin=DiffEqBase.prob2dtmin(prob; use_end_time=true),
                           dtmax=eltype(prob.tspan)((prob.tspan[end]-prob.tspan[1])),
                           beta1 = 7//(10(q+1)),
                           beta2 = 2//(5(q+1)),

                           qoldinit = 1//10^4,

                           sigmarule=:dynamic,

                           progressbar=false,
                           maxiters=1e5,
                           saveat=nothing,
                           save_everystep=true,
                           internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
                           unstable_check = DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK,

                           timeseries_errors=false,
                           dense_errors=false,

                           kwargs...)
    # Init
    bigfloat && (prob = remake(prob, u0=big.(prob.u0)))
    IIP = DiffEqBase.isinplace(prob)

    if method == :ekf1 && isnothing(prob.f.jac)
        error("""EKF1 requires the Jacobian. To automatically generate it with ModelingToolkit.jl
               use ProbNumoDE.remake_prob_with_jac(prob).""")
    end

    if length(prob.u0) == 1 && size(prob.u0) == ()
        @warn "prob.u0 is a scalar; In order to run, we remake the problem with u0 = [u0]."
        prob = remake(prob, u0=[prob.u0])
    end

    (timeseries_errors != false) && @warn("`timeseries_errors` currently not supported")
    (dense_errors != false) && @warn("`dense_errors` currently not supported")

    f = prob.f
    u0 = copy(prob.u0)
    t0, tmax = prob.tspan
    p = prob.p
    d = length(u0)

    # Model
    constants = GaussianODEFilterConstantCache(prob, q, prior, method)

    # Solver Options
    tType = eltype(prob.tspan)
    adaptive = steprule != :constant
    if !adaptive && dt == tType(0)
        error("Fixed timestep methods require a choice of dt")
    end
    steprules = Dict(
        :constant => ConstantSteps(),
        :standard => StandardSteps(),
        :PI => PISteps(),
    )
    steprule = steprules[steprule]

    sigmarules = Dict(
        :dynamic => DynamicSigma(),
        :dynamicMV => MVDynamicSigma(),
        :fixedMLE => MLESigma(),
        :fixedMLEMV => MVMLESigma(),
    )
    sigmarule = sigmarules[sigmarule]

    # Cache
    cache = GaussianODEFilterCache(d, q, prob, constants, initial_sigma(sigmarule, d, q), initialize_derivatives)

    xType = typeof(cache.x)
    sigmaType = typeof(cache.σ_sq)
    measType = typeof(cache.measurement)
    proposalType = NamedTuple{(:t, :prediction, :filter_estimate, :measurement, :H, :Q, :v, :σ², :accept, :dt),
                              Tuple{tType, xType, xType, measType, typeof(cache.H), typeof(cache.Qh),
                                    typeof(cache.h), sigmaType, Bool, tType}}
    empty_proposals = proposalType[]

    destats = DiffEqBase.DEStats(0)

    state_estimates = StructArray([copy(cache.x)])
    times = [t0]
    sigmas = []
    accept_step = false
    retcode = :Default

    isnothing(dtmin) && (dtmin = DiffEqBase.prob2dtmin(prob; use_end_time=true))
    dt_init = dt != 0 ? dt : 1e-3
    QT = tType
    opts = DEOptions{typeof(maxiters), typeof(abstol), typeof(reltol), QT,
    typeof(internalnorm), tType, typeof(unstable_check)}(
        maxiters, adaptive, abstol, reltol, QT(gamma), QT(qmin), QT(qmax),
        QT(beta1), QT(beta2), QT(qoldinit),
        internalnorm, unstable_check, dtmin, dtmax, false, true)

    return ODEFilterIntegrator{IIP, typeof(u0), typeof(t0), typeof(p), typeof(f), QT, typeof(opts), typeof(constants), typeof(cache),
                               typeof(sigmarule), typeof(steprule), typeof(empty_proposals),
                               xType, sigmaType, typeof(prob), typeof(alg)}(
        nothing, f, u0, t0, t0, t0, tmax, dt_init, p, one(QT), QT(qoldinit),
        constants, cache,
        # d, q, dm, mm, sigmarule, steprule,
        opts, sigmarule, steprule, smooth,
        #
        empty_proposals, state_estimates, times, sigmas,
        #
        0, 0, accept_step, retcode, prob, alg, destats,
    )
end

function DiffEqBase.solve!(integ::ODEFilterIntegrator)
    while integ.t < integ.tmax
        loopheader!(integ)
        if check_error!(integ) != :Success
            return integ.sol
        end
        perform_step!(integ)
        loopfooter!(integ)
    end
    postamble!(integ)
    if integ.sol.retcode == :Default
        integ.sol = solution_new_retcode(integ.sol, :Success)
    end
    return integ.sol
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
using Documenter, DataValues

makedocs(modules=[DataValues],
	sitename="DataValues.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ])

deploydocs(repo="github.com/queryverse/DataValues.jl.git")


text:
module AbstractGradOptim

export Adam, update!

include("Adam.jl")
# include("AbstractAPI.jl")
# include("SGD.jl")

end # module


text:
"""
    keyboarddriven()::TERMIOS.termios

Enter keyboard-driven mode.

Sets the terminal to non-canonical mode and turns off echo. Returns the current
settings so they may be restored later on.

"""
keyboarddriven() = keyboarddriven(false)

"""
    keyboarddriven(echo::Bool)::TERMIOS.termios

Enter keyboard-driven mode setting echo.

Sets the terminal to non-canonical mode and sets echo. Returns the current
settings so they may be restored later on.
"""
function keyboarddriven(echo::Bool)::TERMIOS.termios
    current = TERMIOS.termios()
    new = TERMIOS.termios()
    TERMIOS.tcgetattr(stdin, current)
    TERMIOS.tcgetattr(stdin, new)
    if !echo
        new.c_lflag &= ~TERMIOS.ECHO
    end
    new.c_lflag &= ~TERMIOS.ICANON
    setmode(new)

    current
end

"""
    normalmode()::TERMIOS.termios

Enter normal mode.

Sets the terminal to canonical mode and turns on echo. Returns the current
settings so the may be restored later on.
"""
function normalmode()::TERMIOS.termios
    current = TERMIOS.termios()
    new = TERMIOS.termios()
    TERMIOS.tcgetattr(stdin, current)
    TERMIOS.tcgetattr(stdin, new)
    new.c_lflag |= TERMIOS.ECHO
    new.c_lflag |= TERMIOS.ICANON
    setmode(new)

    current
end

"""
    setmode(mode::TERMIOS.termios)

Set a terminal mode.

# Examples

```
initial_mode = keyboarddriven()
# Do stuff in keyboard-driven mode
setmode(initial_mode)
```
"""
setmode(mode::TERMIOS.termios) = TERMIOS.tcsetattr(stdin, TERMIOS.TCSANOW, mode)


text:
using RegressionDiscontinuity
using Test

@testset "RegressionDiscontinuity.jl" begin
    include("test_runningvariable.jl")
    include("test_bandwidth.jl")
    include("test_density_test.jl")
    include("test_minmax.jl")
end


text:
using MultiFloats: Float64x3
using StaticArrays: SVector
using UUIDs: UUID

push!(LOAD_PATH, @__DIR__)
using PCREO
using PointGroups

const UUID_REGEX = Regex("[0-9a-fA-F]{8}-[0-9a-fA-F]{4}-" *
                         "[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12}")

function to_point_vector(points::AbstractMatrix{T}) where {T}
    dimension, num_points = size(points)
    @assert dimension == 3
    return [SVector{3,T}(view(points, :, i)) for i = 1 : num_points]
end

function add_to_database!(filepath::String)
    filename = basename(filepath)
    uuid = UUID(match(UUID_REGEX, filename).match)
    record = PCREORecord(filepath)
    dim_str = lpad(record.dimension, 2, '0')
    num_str = lpad(record.num_points, 8, '0')
    num_dir = joinpath(PCREO_DATABASE_DIRECTORY, num_str)
    if !isdir(num_dir)
        mkdir(num_dir)
    end
    for entry_name in readdir(num_dir)
        entry_dir = joinpath(num_dir, entry_name)
        reference_name = "PCREO-$dim_str-$num_str-$entry_name.csv"
        reference = PCREORecord(joinpath(entry_dir, reference_name))
        if isapprox(record.energy, reference.energy;
                    rtol=Float64x3(2^-80))
            if isometric(to_point_vector(record.points),
                         to_point_vector(reference.points), 2^-40)
                mv(filepath, joinpath(entry_dir, filename))
                return joinpath(entry_dir, filename)
            end
        end
    end
    new_dir = joinpath(num_dir, string(uuid))
    mkdir(new_dir)
    mv(filepath, joinpath(new_dir, filename))
    return joinpath(new_dir, filename)
end

function main(n_min::Int, n_max::Int)
    while true
        filenames = filter(startswith("PCREO-03-"),
                           readdir("D:\\Data\\PCREOCleaned"))
        if length(filenames) > 0
            filename = rand(filenames)
            num_points = parse(Int, filename[10:17])
            if n_min <= num_points <= n_max
                src = joinpath("D:\\Data\\PCREOCleaned", filename)
                try
                    dst = add_to_database!(src)
                    println(src, " => ", dst)
                catch e
                    if e isa AssertionError
                        println(filename, ": ", e)
                    else
                        rethrow(e)
                    end
                end
            end
        else
            sleep(1.0)
        end
    end
end

main(parse(Int, ARGS[1]), parse(Int, ARGS[2]))


text:
# Implements Decision Diagram data structures, constructors, basic manipulation
export
    DecisionDiagram,
    Terminal,
    Node,
    index,
    value,
    pliteral,
    nliteral,
    indicator

nextid = 1

abstract type DecisionDiagram{T} end

"""
  Implements a generic decision diagram whose leaves are of given type T
"""
struct Terminal{T} <: DecisionDiagram{T}
    "Unique identifier."
    id::Int
    "Constant."
    value::T
    # "Constructs a terminal."
    # Terminal{T}(v) where T = new(typemax(Int),v)
    function Terminal{T}(v)  where T
        global nextid
        id = nextid
        nextid = (nextid + 1) % typemax(Int)
        new{T}(id,v)
    end
end
Terminal(v::T) where T = Terminal{T}(v)
"Multiplicative identity for terminals of type T"
@inline Base.one(α::DecisionDiagram{T}) where {T <: Number} = Terminal(one(T))
"Additive identity for terminals of type T"
@inline Base.zero(α::DecisionDiagram{T}) where {T <: Number} = Terminal(zero(T))
"Returns the sum of the respective functions."
# GETTERS
@inline value(α::Terminal) = α.value
@inline index(α::Terminal) = typemax(Int)
# OUTPUT
function Base.show(io::IO,α::Terminal)
    print(io,"(",value(α),") [", id(α), "]")
    nothing
end
struct Node{T} <: DecisionDiagram{T}
    "Unique identifier."
    id::Int
    "Root vertex variable index."
    index::Int
    "Low child vertex of BDD."
    low::DecisionDiagram{T}
    "High child vertex of BDD."
    high::DecisionDiagram{T}
    "Constructs a variable."
    function Node{T}(i::Integer, low::DecisionDiagram{T}, high::DecisionDiagram{T}) where T
        global nextid
        id = nextid
        nextid = (nextid + 1) % typemax(Int)
        new{T}(id,i,low,high)
    end
    Node{T}(id,index,low,high) where T = new(id,index,low,high)
end
# CONSTRUCTORS
Node(index,low::DecisionDiagram{T},high::DecisionDiagram{T}) where T = Node{T}(index,low,high)
Node(index,low::T,high::T) where T = Node{T}(index,Terminal(low),Terminal(high))
Node(index,low::T,high::DecisionDiagram{T}) where T = Node{T}(index,Terminal(low),high)
Node(index,low::DecisionDiagram{T},high::T) where T = Node{T}(index,low,Terminal(high))
"""
   Creates Diagram representing indicator function on given variable.
   If complement is true, then returns complement indicator.
"""
indicator(T::Type,index::Int,complement::Bool=false) = complement ? Node(index,Terminal(one(T)),Terminal(zero(T))) : Node(index,Terminal(zero(T)),Terminal(one(T))) # Numeric terminals
""""
    pliteral(index::Int)

Creates a positive literal (an indicator over Boolean constants) for variable index.
"""
pliteral(index::Int) = indicator(Bool,index) # Positive Literal
""""
    nliteral(index::Int)

Creates a negated literal (an indicator over Boolean constants) for variable index.
"""
nliteral(index::Int) = indicator(Bool,index,true) # Negative Literal
# GETTER functions
id(α::DecisionDiagram) = α.id
index(α::Node) = α.index
low(α::Node) = α.low
high(α::Node) = α.high
# OUTPUT
Base.show(io::IO,α::Node) = print(io,"\n@"), Base.show(io,α," ")
function Base.show(io::IO,α::Node,prefix::String)
    print(io,"╮ V", index(α), " [", id(α), "] \n",prefix, "├─⊂−⊃─")
    show(io,low(α),prefix * "│     ")
    print(io,"\n",prefix, "╰─⊂+⊃─")
    show(io,high(α),prefix * "      ")
    nothing
end
function Base.show(io::IO,α::Terminal,prefix::String)
    show(io,α)
    nothing
end

# ITERATORS
mutable struct IteratorDiagram{T}
    frontier::Vector{DecisionDiagram{T}}
    visited::Set{Int}
    count::Int
end
"""
    Depth-first search traversal.
"""
function Base.iterate(α::DecisionDiagram{T}, state = IteratorDiagram(DecisionDiagram{T}[α], Set{Int}(), 0)) where T
    if isempty(state.frontier)
        return nothing
    end
    node = pop!(state.frontier) #  change to popfirst! to perform BFS instead
    # get next unvisited state
    while id(node) in state.visited
        if isempty(state.frontier) return nothing end
        node = pop!(state.frontier)
    end
    # mark state as visited
    push!(state.visited, id(node))
    state.count += 1
    if isa(node, Node) # not a terminal
        if !(id(low(node)) in state.visited)
            push!(state.frontier,low(node))
        end
        if !(id(high(node)) in state.visited)
            push!(state.frontier,high(node))
        end
    end
    (node, state)
end
"We can't know the number of nodes in advance (to traversing the diagram)."
Base.IteratorSize(::DecisionDiagram{T}) where T = Base.SizeUnknown()
"""
    A decision diagram is composed of other decision diagrams.
"""
Base.eltype(α::Type{DecisionDiagram{T}}) where T = DecisionDiagram{T}
"Returns the number of nodes in the diagram."
Base.length(α::DecisionDiagram{T}) where T = length(collect(α))


text: codegen(::Type{Val{:iterate}}, job::BasicMCJob) = codegen(Val{:iterate}, typeof(job.sampler), job)


text:
using MPI

testdir = @__DIR__

mpiexec() do cmd
    run(`$cmd -n 3 $(Base.julia_cmd()) $(joinpath(testdir, "mpi_queue.jl"))`)
end

for file in readdir()
    if endswith(file, "log")
        rm(file)
    end
end

mpiexec() do cmd
    run(`$cmd -n 4 $(Base.julia_cmd()) $(joinpath(testdir, "mpi_updates.jl"))`)
end


text:
"""
    Raibert hopper policy for Hopper2D
"""

@with_kw mutable struct RaibertOptions
    live_plotting::Bool=false # Use the live plotting tool to debug
end

mutable struct Raibert <: Policy
	kr_c::Real
	kr_p::Real
	kr_v_stance::Real
	kr_v_flight::Real
	kθ_c::Real
	kθ_p::Real
	kθ_v::Real
	u::AbstractVector
	contact::Bool
	v0::Real
	Tstance::Real
	Tflight::Real
	q0::AbstractVector
	q1::AbstractVector
	opts::RaibertOptions
end

function raibert_policy(model::Hopper2D;
		v0=0.5, Tstance=0.13, Tflight=0.62, h = 0.1,
		kr_c =  8e1,
		kr_p = -1e3,
		kr_v_stance = -1e-2,
		kr_v_flight = -1e1,
		kθ_c =  0.0,
		kθ_p = -3e1,
		kθ_v = -1e1,
		mpc_opts = RaibertOptions(),
		)

	u  = zeros(model.dim.u)
	q0 = zeros(model.dim.q)
	q1 = zeros(model.dim.q)
	contact = false

	Raibert(kr_c, kr_p, kr_v_stance, kr_v_flight, kθ_c, kθ_p, kθ_v,
		u, contact, v0, Tstance, Tflight, copy(q0), copy(q1), mpc_opts)
end

function policy(p::Raibert, x, traj, t)
	# Initialization
	h = traj.h
	p.q0 = copy(traj.q[t])
	p.q1 = copy(x)

	# Detect contact
	if any(traj.γ[max(1,t-1)] .> 1.5e-2)
		p.contact = true
	else
		p.contact = false
	end

	# Velocities
	qv = (p.q1 - p.q0)/h
	θv = qv[3]
	rv = qv[4]

	# References
	θ1 = p.q1[3]
	r1 = p.q1[4]
	rref = 0.5
	# td = touchdown
	θtd  = asin(p.v0*Tstance/(2*rref))/2

	# Gains
	kr_c = p.kr_c
	kr_p = p.kr_p
	kr_v_stance = p.kr_v_stance
	kr_v_flight = p.kr_v_flight

	kθ_c = p.kθ_c
	kθ_p = p.kθ_p
	kθ_v = p.kθ_v

	if p.contact
		# regulate around θtd using same gain kθ_p but taking into account
		# that stance time is much lower than flight time
		p.u[1] = kθ_c + kθ_p*(θ1 + θtd) * p.Tflight/p.Tstance
		p.u[2] = kr_c + kr_p*(r1 - rref) + kr_v_stance*rv
	else
		# regulate around θtd
		p.u[1] = kθ_p*(θ1 - θtd)  + kθ_v*θv
		p.u[2] = kr_p*(r1 - rref) + kr_v_flight*rv
	end

	return p.u .* h
end


text:
using Test
using TextSummarizer

docx = "The nativity of Jesus or birth of Jesus is described in the gospels of Luke and Matthew. The two accounts differ, but agree that Jesus was born in Bethlehem in Judea during the reign of King Herod the Great, his mother Mary was married to a man named Joseph, who was descended from King David and was not his biological father, and that his birth was caused by divine intervention.

       Luke's version says the birth took place during the Census of Quirinius, mentions an annunciation to the shepherds, the presentation of Jesus at the Temple, and gives the name of the angel who announces the coming birth to Mary (the archangel Gabriel). Matthew's version mentions the arrival of the Magi, the flight into Egypt by the family, and the Massacre of the Innocents by King Herod. The consensus of scholars is that both gospels were written about AD 75-85, and while it is possible that one account might be based on the other, or that the two share common source material, the majority conclusion is that the two nativity narratives are independent of each other.

       The nativity is the basis for the Christian holiday of Christmas and plays a major role in the Christian liturgical year. Christian congregations of the Western tradition (including the Catholic Church, the Western Rite Orthodox, the Anglican Communion, and many other Protestants, such as the Moravian Church) begin observing the season of Advent four Sundays before Christmas, the traditional feast-day of his birth, which falls on December 25."

@testset "Text Summarization " begin
    @test length(docx) > length(TextSummarizer.text_summarizer(docx))
    @test occursin(TextSummarizer.text_summarizer(docx),docx) == false
	@test TextSummarizer.text_summarizer(docx) !=  docx
end


text:
using DrWatson
@quickactivate "starter-kit-regions"
DrWatson.greet()

max_stars_repo_name: greimel/starter-kit-regions

max_issues_repo_name: greimel/starter-kit-regions

max_forks_repo_name: greimel/starter-kit-regions


text:
function initial_placement!(placement_struct)
    # Build graph
    graph, node_dict, component_dict = build_graph(placement_struct)
    # Run Bipartite Matching
    bipartite_match!(graph)
    @debug "Bipartite Matching Complete"
    do_assignment(placement_struct, graph, node_dict, component_dict)
    return graph
end

#=
Build light graph and translation tables to pass to bipartite_match!

        Taskgraph           Component
        Node                Indices
        Indices

        /---- a    ⋯    1 ----\
       /----- b    ⋯    2 -----\
source ------ c    ⋯    3 ------- sink
       \----- d    ⋯    4 -----/
        \---- e    ⋯    5 ----/
          ⋮                ⋮

An edge is created from a TaskgraphNode on the left to a component index
on the right if that Node can be mapped to the component represented by the
index.

Note: Bipartite matching is based on a maximum network flow algorithm, which
works with a single master source and master sink. Thus, we must observe the
following invariants:

* node index 1 is the ultimate source of the bipartite match.
* node index 2 is the ultimate sink of the bipartite match.
=#
function build_graph(sa::SAStruct)
    # Offset for the source and sink nodes.
    sourcesink_offset = 2

    # Create the node translation dictionary
    # Maps node_index to an index in the LightGraph
    node_dict = Dict(i => i + sourcesink_offset for i in 1:length(sa.nodes))

    # Maps Component to an index in the LightGraph
    component_dict = Dict{eltype(sa),Int64}()

    # Build the light graph. Initially populate with the master source, master
    # sink, and vertices representing the TaskgraphNodes
    graph = DiGraph(sourcesink_offset + length(node_dict))
    edges_added = 0

    for (index, node) in enumerate(sa.nodes)
        # Get the vertex number for this node.
        node_number = node_dict[index]

        # Get all locations that this node type can be mapped to.
        for key in getlocations(sa.maptable, getclass(node))
            # Add a vertex for the component if it hasn't been added yet.
            if !haskey(component_dict, key)
                add_vertex!(graph)
                component_dict[key] = nv(graph)
            end
            # Add an edge from node to component.
            add_edge!(graph, node_number, component_dict[key])
            edges_added += 1
        end
    end

    # Create the source and sink edges
    for i in values(node_dict)
        add_edge!(graph, 1, i)
    end
    for i in values(component_dict)
        add_edge!(graph, i, 2)
    end

    # Prints for debug.
    @debug """
        Translation Tables Complete.
        Number of Vertices: $(nv(graph)).
        Number of Tasks: $(length(node_dict)).
        Number of Components: $(length(component_dict)).
        Number of edges added: $(edges_added).
        """
    return graph, node_dict, component_dict
end

function do_assignment(placement_struct, graph, node_dict, component_dict)
    # Reverse the node and component dictionaries.
    node_dict_rev = rev_dict(node_dict)
    component_dict_rev = rev_dict(component_dict)
    # Iterate through all the indices belonging to nodes.
    for i in keys(node_dict_rev)
        # Iterate through all the "inneighbors" of this node.
        # Looking for an incoming edge that is not from the source.
        # This edge leads to the component this task will be mapped to.
        for neighbor in inneighbors(graph, i)
            neighbor == 1 && continue
            b = neighbor
            # Get the address and component number from the reversed
            # component dictionary.
            position = component_dict_rev[b]
            node = node_dict_rev[i]
            assign(placement_struct, node, position)
        end
    end
    return nothing
end

function bipartite_match!(g::AbstractGraph)
    ####################
    # Graph Definition #
    ####################
    # the graph g contains LHS vertices (referred to as "a"),
    #   RHS vertices (referred to as "b"), source vertex, and sink vertex
    # vertex #1 (source) is connected to all members of a (source -> v ∈ a)
    # vertex #2 (sink) is connected to all members of b (v ∈ b -> sink)

    ###################
    # Algorithm rules #
    ###################
    # the algorithm can go from left vertex to right vertex if they are
    #   connected only with "->" edge
    # once that path is being used, another edge "<-" is added
    # the algorithm can go from right vertex to left vertex if they are
    #   connected with both "->" and "<-" edges
    # once that path is being used, the edge "<-" is removed
    # the algorithm ends after exhausting all the "a" vertices
    # "a" set and "b" set will have a one-to-one relation when matching is complete
    # the matched pair will be connected with both "->" and "<-" edges

    # source and sink are defined as 1 and 2 respectively in LightGraph g
    source = 1
    sink = 2
    # group the vertices according to their number of preferences
    sort_dict = SortedDict{Int64,Array{Int64,1}}()
    for a in outneighbors(g, source)
        count = length(outneighbors(g, a))
        push_to_dict(sort_dict, count, a)
    end
    # unwrap vertex value array for each key in sort_dict
    a_set = Int64[] # all vertices are stored in a_set
    for v in values(sort_dict)
        append!(a_set, v)
    end
    # source -> a -> b -> sink links are made
    # if link requires a backward trace, it's skipped for now (taken care of later)
    for a in a_set
        for b in outneighbors(g, a)
            if has_edge(g, b => sink) && !has_edge(g, sink => b)
                add_edge!(g, a => source)
                add_edge!(g, b => a)
                add_edge!(g, sink => b)
                break
            elseif has_edge(g, b => sink) && has_edge(g, sink => b)
                continue
            end
        end
    end
    # links that require backward tracing are created below
    for a in outneighbors(g, source)
        has_edge(g, a => source) && continue
        two_found = false # initialize
        for b in outneighbors(g, a)
            two_found && break
            predecessor = Int64[] # create an array to keep track the path
            neighbor = b
            previous_neighbor = a
            exit = false # initialize
            while !two_found && !exit
                for new_neighbor in outneighbors(g, neighbor)
                    # check if neighbor can be connected to sink
                    if (new_neighbor == sink && !has_edge(g, sink => neighbor))
                        two_found = true
                        push!(predecessor, previous_neighbor)
                        push!(predecessor, neighbor)
                        push!(predecessor, new_neighbor)
                        break
                    end
                    # check if the vertex has no remaining moves
                    if (
                        (outneighbors(g, neighbor)) == [source, previous_neighbor] ||
                        (inneighbors(g, neighbor)) == [sink, previous_neighbor]
                    )
                        exit = true
                        break
                    end
                    # source, sink, and previous_neighbor are not valid vertices
                    # to move to
                    if (
                        new_neighbor == source ||
                        new_neighbor == sink ||
                        new_neighbor == previous_neighbor
                    )
                        continue
                    end
                    if (
                        has_edge(g, neighbor => sink) &&
                        has_edge(g, neighbor => new_neighbor) &&
                        has_edge(g, new_neighbor => neighbor)
                    )
                        push!(predecessor, previous_neighbor)
                        previous_neighbor = neighbor
                        neighbor = new_neighbor
                    elseif (
                        has_edge(g, source => neighbor) &&
                        has_edge(g, neighbor => new_neighbor) &&
                        !has_edge(g, new_neighbor => neighbor)
                    )
                        push!(predecessor, previous_neighbor)
                        previous_neighbor = neighbor
                        neighbor = new_neighbor
                    else
                        exit = true
                        break
                    end # end of if statement
                end # end of 3rd for loop
            end # end of while loop
            # if sink is found, trace back the predecessors and create appropriate
            # edges as mentioned in the "Algorithm Rules" above
            if two_found
                add_edge!(g, predecessor[1] => source)
                add_edge!(g, predecessor[2] => predecessor[1])
                add_edge!(g, sink => predecessor[end])
                for i in 2:(length(predecessor) - 1)
                    if i in outneighbors(g, source)
                        add_edge!(g, predecessor[i], predecessor[i + 1])
                    elseif i in inneighbors(g, sink)
                        rem_edge!(g, predecessor[i], predecessor[i + 1])
                    end # end if
                end # end for loop
                break
            end # if two_found
        end # end of 2nd for loop
    end # end of 1st for loop
    return g
end


text:
# Attempt to compute a reasonable default mu: at the starting
# position, the gradient of the input function should dominate the
# gradient of the barrier.
function initialize_mu{T}(gfunc::Array{T}, gbarrier::Array{T}, ops::Options)
    mu0::T
    mu0factor::T
    @defaults ops mu0=nan(T) mu0factor=0.001
    if isnan(mu0)
        gbarriernorm = sum(abs(gbarrier))
        if gbarriernorm > 0
            mu = mu0factor*sum(abs(gfunc))/gbarriernorm
        else
            # Presumably, there is no barrier function
            mu = zero(T)
        end
    else
        mu = mu0
    end
    @check_used ops
    return mu
end

function barrier_box{T}(g, x::Array{T}, l::Array{T}, u::Array{T})
    n = length(x)
    calc_grad = !(g === nothing)

    v = zero(T)
    for i = 1:n
        thisl = l[i]
        if isfinite(thisl)
            dx = x[i] - thisl
            if dx <= 0
                return inf(T)
            end
            v -= log(dx)
            if calc_grad
                g[i] = -one(T)/dx
            end
        else
            if calc_grad
                g[i] = zero(T)
            end
        end
        thisu = u[i]
        if isfinite(thisu)
            dx = thisu - x[i]
            if dx <= 0
                return inf(T)
            end
            v -= log(dx)
            if calc_grad
                g[i] += one(T)/dx
            end
        end
    end
    return v
end

function function_barrier{T}(gfunc, gbarrier, x::Array{T}, f::Function, fbarrier::Function)
    vbarrier = fbarrier(gbarrier, x)
    if isfinite(vbarrier)
        vfunc = f(gfunc, x)
    else
        vfunc = vbarrier
    end
    return vfunc, vbarrier
end

function barrier_combined{T}(g, gfunc, gbarrier, valeach::Vector{T}, x::Array{T}, fb::Function, mu::T)
    calc_grad = !(g === nothing)
    valfunc, valbarrier = fb(gfunc, gbarrier, x)
    valeach[1] = valfunc
    valeach[2] = valbarrier
    if calc_grad
        for i = 1:length(g)
            g[i] = gfunc[i] + mu*gbarrier[i]
        end
    end
    return convert(T, valfunc + mu*valbarrier) # FIXME make this unnecessary
end

function limits_box{T}(x::Array{T}, d::Array{T}, l::Array{T}, u::Array{T})
    alphamax = inf(T)
    for i = 1:length(x)
        if d[i] < 0
            alphamax = min(alphamax, ((l[i]-x[i])+eps(l[i]))/d[i])
        elseif d[i] > 0
            alphamax = min(alphamax, ((u[i]-x[i])-eps(u[i]))/d[i])
        end
    end
    epsilon = eps(max(alphamax, one(T)))
    if !isinf(alphamax) && alphamax > epsilon
        alphamax -= epsilon
    end
    return alphamax
end

# Default preconditioner for box-constrained optimization
# This creates the inverse Hessian of the barrier penalty
function precondprepbox(P, x, l, u, mu)
    for i = 1:length(x)
        xi = x[i]
        li = l[i]
        ui = u[i]
        P[i] = 1/(mu*(1/(xi-li)^2 + 1/(ui-xi)^2) + 1) # +1 like identity far from edges
    end
end

const PARAMETERS_MU = one64<<display_nextbit
display_nextbit += 1

function fminbox{T}(func::Function, x::Array{T}, l::Array{T}, u::Array{T}, ops::Options)
    tol::T
    mufactor::T
    @defaults ops tol=eps(T)^(2/3) mufactor=0.001 display=0 optimizer=(func, x, ops)->cgdescent(func, x, ops) P=Array(T,length(x)) precondprep=precondprepbox
    ops = copy(ops)  # to avoid passing back extended options
    x = copy(x)
    fbarrier = (gbarrier, x) -> barrier_box(gbarrier, x, l, u)
    fb = (gfunc, gbarrier, x) -> function_barrier(gfunc, gbarrier, x, func, fbarrier)
    gfunc = similar(x)
    gbarrier = similar(x)
    # Because we use the gradient to estimate the initial mu, we have
    # to be careful about one special case that might occur commonly
    # in practice: the initial guess x is exactly in the center of the
    # box. In that case, gbarrier is zero. But since the
    # initialization only makes use of the magnitude, we can fix this
    # by using the sum of the absolute values of the contributions
    # from each edge.
    for i = 1:length(gbarrier)
        thisx = x[i]
        thisl = l[i]
        thisu = u[i]
        if thisx < thisl || thisx > thisu
            error("Initial position must be inside the box")
        end
        gbarrier[i] = (isfinite(thisl) ? one(T)/(thisx-thisl) : zero(T)) + (isfinite(thisu) ? one(T)/(thisu-thisx) : zero(T))
    end
    valfunc = func(gfunc, x)
    mu = initialize_mu(gfunc, gbarrier, ops)
    if display > 0
        println("######## fminbox ########")
        println("Initial mu = ", mu)
    end

    g = similar(x)
    valboth = Array(T, 2)
    fval_all = Array(Vector{T}, 0)
    fcount_all = 0
    xold = similar(x)
    converged = false
    @set_options ops tol=10*tol alphamaxfunc=(x, d)->limits_box(x, d, l, u) P=P reportfunc=val->valboth[1]
    while true
        copy!(xold, x)
        # Optimize with current setting of mu
        funcc = (g, x) -> barrier_combined(g, gfunc, gbarrier, valboth, x, fb, mu)
        @set_options ops precondprep=(out, x)->precondprep(out, x, l, u, mu)
        if display > 0
            println("#### Calling optimizer with mu = ", mu, " ####")
        end
        x, fval, fcount, converged = optimizer(funcc, x, ops)
        if display & PARAMETERS_MU > 0
            println("x: ", x)
        end
        push!(fval_all, fval)
        fcount_all += fcount

        # Decrease mu
        mu *= mufactor

        # Test for convergence
        fcmp = abs(fval[1]) + abs(fval[end])
        tot = zero(T)
        for i = 1:length(x)
            tot += abs((x[i] - xold[i]) * (gfunc[i] + mu*gbarrier[i]))
        end
        if tot <= tol * (fcmp + eps(T))
            break
        end
    end
    @check_used ops
    return x, fval_all, fcount_all, converged
end
fminbox{T}(func::Function, x::Array{T}, l::Array{T}, u::Array{T}) = fminbox(func, x, l, u, Options())
export fminbox


text:
"""
    moments!(height, velx, vely, fout)

Computation of the hydrodynamic moments, `height` and `velocity`.

# Mathematics

The macroscopic quantities such as the height and the velocity are the moments of the distribution function,

`` h(\\mathbf{x},t) = \\sum_{i=0}^8 f_i , ``

and

`` \\mathbf{v}(\\mathbf{x},t) = \\frac{1}{h}\\sum_{i=0}^8 \\mathbf{c}_i f_i . ``

# Examples

```jldoctest
julia> using Swalbe, Test

julia> fout = zeros(5,5,9); fout[:,:,1] .= 1.0; fout[:,:,2] .= 0.1; # Dist with artifical velocity in x

julia> height = zeros(5,5); velx = zeros(5,5); vely = zeros(5,5);

julia> Swalbe.moments!(height,velx,vely,fout)

julia> @test all(height .== 1.1)
Test Passed

julia> @test all(velx .== 0.1/1.1)
Test Passed

julia> @test all(vely .== 0.0)
Test Passed
```

# References

- [Krüger](https://www.springer.com/gp/book/9783319446479)
- [Salmon](https://www.ingentaconnect.com/contentone/jmr/jmr/1999/00000057/00000003/art00005#)
- [Zitz, Scagliarini and Harting](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.033313)
"""
function moments!(height, velx, vely, fout)
    # Get views of the populations
    f0, f1, f2, f3, f4, f5, f6, f7, f8 = Swalbe.viewdists(fout)
    # Compute the height
    sum!(height, fout)
    # and the velocities (as simple as possible)
    velx .= (f1 .- f3 .+ f5 .- f6 .- f7 .+ f8) ./ height
    vely .= (f2 .- f4 .+ f5 .+ f6 .- f7 .- f8) ./ height
    return nothing
end

function moments!(height, vel, fout)
    # Get views of the populations
    f0, f1, f2 = Swalbe.viewdists_1D(fout)
    # Compute the height
    sum!(height, fout)
    # and the velocities (as simple as possible)
    vel .= (f1 .- f2) ./ height
    return nothing
end


text:
using Genser
using Documenter

push!(LOAD_PATH,"../src/")

DocMeta.setdocmeta!(Genser, :DocTestSetup, :(using Genser); recursive=true)

makedocs(;
    modules=[Genser],
    authors="Mark Wardle <mark@potient.com> and contributors",
    repo="https://github.com/wardlem/Genser.jl/blob/{commit}{path}#{line}",
    sitename="Genser.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wardlem.github.io/Genser.jl",
        assets=String[],
    ),
    pages=[
        "Genser" => "index.md",
        "Manual" => [
            "manual/overview.md",
            "manual/datamodel.md",
        ],
        "Library" => [
            "lib/primary.md",
            "lib/datamodel.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/wardlem/Genser.jl",
    devbranch="main",
)


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

name = "GR"
version = v"0.61.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sciapp/gr.git", "c1a183066bb7101d6805f2cb06ff950cbef6bab9"),
    FileSource("https://github.com/sciapp/gr/releases/download/v$version/gr-$version.js",
               "f0b0d7dc53077a07fa028125ad935caba0e1228ffd4c1503687af39681c8d7a0", "gr.js")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gr

if test -f "$prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake"; then
    sed -i 's/_qt5gui_find_extra_libs.*AGL.framework.*//' $prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake
fi

update_configure_scripts

make -C 3rdparty/qhull -j${nproc}

if [[ $target == *"mingw"* ]]; then
    winflags=-DCMAKE_C_FLAGS="-D_WIN32_WINNT=0x0f00"
    tifflags=-DTIFF_LIBRARY=${libdir}/libtiff-5.dll
else
    tifflags=-DTIFF_LIBRARY=${libdir}/libtiff.${dlext}
fi

if [[ "${target}" == *apple* ]]; then
    make -C 3rdparty/zeromq ZEROMQ_EXTRA_CONFIGURE_FLAGS="--host=${target}"
fi

if [[ "${target}" == arm-* ]]; then
    export CXXFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib"
fi

mkdir build
cd build
cmake $winflags -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DGR_USE_BUNDLED_LIBRARIES=ON $tifflags -DCMAKE_BUILD_TYPE=Release ..

VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
cp ../../gr.js ${libdir}/

install_license $WORKSPACE/srcdir/gr/LICENSE.md

if [[ $target == *"apple-darwin"* ]]; then
    cd ${bindir}
    ln -s ../Applications/gksqt.app/Contents/MacOS/gksqt ./
    ln -s ../Applications/GKSTerm.app/Contents/MacOS/GKSTerm ./
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("armv7l",  "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64",  "linux"; libc="glibc"),
    Platform("i686",  "linux"; libc="glibc"),
    Platform("powerpc64le",  "linux"; libc="glibc"),
    Platform("x86_64",  "windows"),
    Platform("i686",  "windows"),
    Platform("x86_64",  "macos"),
    Platform("aarch64", "macos"),
    Platform("x86_64",  "freebsd"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libGR", :libGR),
    LibraryProduct("libGR3", :libGR3),
    LibraryProduct("libGRM", :libGRM),
    LibraryProduct("libGKS", :libGKS),
    ExecutableProduct("gksqt", :gksqt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Bzip2_jll"; compat="1.0.8"),
    Dependency("Cairo_jll"; compat="1.16.1"),
    Dependency("FFMPEG_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("GLFW_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("libpng_jll"),
    Dependency("Libtiff_jll"; compat="4.3.0"),
    Dependency("Pixman_jll"),
#    Dependency("Qhull_jll"),
    Dependency("Qt5Base_jll"),
    BuildDependency("Xorg_libX11_jll"),
    BuildDependency("Xorg_xproto_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
# GCC version 7 because of ffmpeg, but building against Qt requires v8 on Windows.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version = v"8", julia_compat="1.6")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
abstract type AbstractIFParameter end
@with_kw struct IFParameter <: AbstractIFParameter
    τm::SNNFloat = 20ms
    τe::SNNFloat = 5ms
    τi::SNNFloat = 10ms
    Vt::SNNFloat = -50mV
    Vr::SNNFloat = -60mV
    El::SNNFloat = Vr
end

abstract type AbstractIF end
@with_kw mutable struct IF <: AbstractIF
    param::IFParameter = IFParameter()
    N::SNNInt = 100
    v::Vector{SNNFloat} = param.Vr .+ rand(N) .* (param.Vt - param.Vr)
    ge::Vector{SNNFloat} = zeros(N)
    gi::Vector{SNNFloat} = zeros(N)
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{SNNFloat} = zeros(N)
    records::Dict = Dict()
end

function integrate!(p::IF, param::IFParameter, dt::SNNFloat)
    @unpack N, v, ge, gi, fire, I = p
    @unpack τm, τe, τi, Vt, Vr, El = param
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
include(joinpath(@__DIR__, "../src/utils/data.jl"))
import .UtilsData: TestData

download(TestData)


text:
"""
    DerivativeUnitRange{m, n} <: AbstractDifferentialOp

Specifies a range of derivatives.

# Examples

Two ways of constructing derivative ranges:

```jldoctest
julia> Derivative(2):Derivative(4)
Derivative(2:4)

julia> Derivative(2:4)
Derivative(2:4)

julia> Tuple(Derivative(2:4))
(D{2}, D{3}, D{4})
```
"""
struct DerivativeUnitRange{m, n} <: AbstractDifferentialOp end

Base.show(io::IO, ::DerivativeUnitRange{m, n}) where {m, n} =
    print(io, "Derivative($m:$n)")

@inline Base.:(:)(::Derivative{m}, ::Derivative{n}) where {m, n} =
    DerivativeUnitRange{m, n}()

@inline Derivative(r::UnitRange{<:Integer}) =
    Derivative(r.start):Derivative(r.stop)

@generated function Base.Tuple(::DerivativeUnitRange{m, n}) where {m, n}
    ex = :(())
    for i ∈ m:n
        push!(ex.args, :(Derivative($i)))
    end
    ex
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines