text:
# Convert a 1-D timeseries into sliding window matrix for ML training
# using Plots

const gAggDict = Dict(
    :median => Statistics.median,
    :mean =>   Statistics.mean,
    :maximum => Statistics.maximum,
    :minimum => Statistics.minimum
)

mutable struct Matrifier <: Transformer
  model
  args

  function Matrifier(args=Dict())
    default_args = Dict{Symbol,Any}(
        :ahead => 1,
        :size => 7,
        :stride => 1,
    )
    new(nothing,mergedict(default_args,args))
  end
end

function fit!(mtr::Matrifier,xx::T,y::Vector=Vector()) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input is not a dataframe")
  x = deepcopy(xx[:Value])
  x isa Vector || error("data should be a vector")
  mtr.model = mtr.args
end

function transform!(mtr::Matrifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input is not a dataframe")
  x = deepcopy(xx[:Value])
  x isa Vector || error("data should be a vector")
  mtype = eltype(x)
  res=toMatrix(mtr,x)
  resarray=convert(Array{mtype},res) |> DataFrame
  rename!(resarray,names(resarray)[end] => :output)
end

function toMatrix(mtr::Transformer, x::Vector)
  stride=mtr.args[:stride];sz=mtr.args[:size];ahead=mtr.args[:ahead]
  @assert stride>0 && sz>0 && ahead > 0
  xlength = length(x)
  xlength > sz || error("data too short for the given size of sliding window")
  ndx=collect(xlength:-1:1)
  mtuples = slidingwindow(i->(i-ahead),ndx,sz,stride)
  height=size(mtuples)[1]
  mmatrix = Array{Union{DateTime,<:Real},2}(zeros(height,sz+1))
  ctr=1
  gap = xlength - mtuples[1][2][1]
  for (s,k) in mtuples
    v = [reverse(s);k] .+ gap
    mmatrix[ctr,:].=x[v]
    ctr+=1
  end
  mmatrix
end

### ====

# Convert a 1-D date series into sliding window matrix for ML training
mutable struct Dateifier <: Transformer
  model
  args

  function Dateifier(args=Dict())
    default_args = Dict{Symbol,Any}(
        :ahead => 1,
        :size => 7,
        :stride => 1
    )
    new(nothing,mergedict(default_args,args))
  end
end

function fit!(dtr::Dateifier,xx::T,y::Vector=[]) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input not a dataframe")
  x = deepcopy(xx[:Date])
  (eltype(x) <: DateTime || eltype(x) <: Date) || error("array element types are not dates")
  dtr.args[:lower] = minimum(x)
  dtr.args[:upper] = maximum(x)
  dtr.model = dtr.args
end

# transform to day of the month, day of the week, etc
function transform!(dtr::Dateifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input not a dataframe")
  x = deepcopy(xx[:Date])
  x isa Vector || error("data should be a vector")
  @assert eltype(x) <: DateTime || eltype(x) <: Date
  res=toMatrix(dtr,x)
  endpoints = convert(Array{DateTime},res)[:,end-1]
  dt = DataFrame()
  dt[:year]=Dates.year.(endpoints)
  dt[:month]=Dates.month.(endpoints)
  dt[:day]=Dates.day.(endpoints)
  dt[:hour]=Dates.hour.(endpoints)
  dt[:week]=Dates.week.(endpoints)
  dt[:dow]=Dates.dayofweek.(endpoints)
  dt[:doq]=Dates.dayofquarter.(endpoints)
  dt[:qoy]=Dates.quarterofyear.(endpoints)
  dtr.args[:header] = names(dt)
  #convert(Matrix{Int64},dt)
  return dt
end

### ====

# Date,Val time series
mutable struct DateValgator <: Transformer
  model
  args
  function DateValgator(args=Dict())
    default_args = Dict{Symbol,Any}(
        :dateinterval => Dates.Hour(1),
        :aggregator => :median
    )
    new(nothing,mergedict(default_args,args))
  end
end

function validdateval!(x::T) where {T<:DataFrame}
  size(x)[2] == 2 || error("Date Val timeseries need two columns")
  (eltype(x[:,1]) <: DateTime || eltype(x[:,1]) <: Date) || error("array element types are not dates")
  eltype(x[:,2]) <: Union{Missing,Real} || error("array element types are not values")
  cnames = names(x)
  rename!(x,Dict(cnames[1]=>:Date,cnames[2]=>:Value))
end

function fit!(dvmr::DateValgator,xx::T,y::Vector=[]) where {T<:Union{Matrix,DataFrame}}
  x = deepcopy(xx)
  validdateval!(x)
  aggr = dvmr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  dvmr.model=dvmr.args
end

function transform!(dvmr::DateValgator,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # make sure aggregator function exists
  aggr = dvmr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  # get the Statistics function
  aggfn = gAggDict[aggr]
  # pass the aggregator function to the generic aggregator function
  fn = aggregatorclskipmissing(aggfn)
  grpby = typeof(dvmr.args[:dateinterval])
  sym = Symbol(grpby)
  x[sym] = round.(x[:Date],grpby)
  aggr=by(x,sym,MeanValue = :Value=>fn)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = round(minimum(x[:Date]),grpby)
  upper = round(maximum(x[:Date]),grpby)
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dvmr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  joined
end

### ====

# Date,Val time series
# Normalize and clean date,val by replacing missings with medians
mutable struct DateValizer <: Transformer
  model
  args

  function DateValizer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :medians => DataFrame(),
        :dateinterval => Dates.Hour(1)
    )
    new(nothing,mergedict(default_args,args))
  end
end

function getMedian(t::Type{T},xx::DataFrame) where {T<:Union{TimePeriod,DatePeriod}}
  x = deepcopy(xx)
  sgp = Symbol(t)
  fn = Dict(Dates.Second=>Dates.second,
            Dates.Minute=>Dates.minute,
            Dates.Hour=>Dates.hour,
            Dates.Day=>Dates.day,
            Dates.Month=>Dates.month)
  try
    x[sgp]=fn[t].(x[:Date])
  catch
    error("unknown dateinterval")
  end
  gpmeans = by(x,sgp,Value = :Value => skipmedian)
  gpmeans
end

function fullaggregate!(dvzr::DateValizer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  x[sym] = round.(x[:Date],grpby)
  aggr = by(x,sym,MeanValue = :Value=>skipmedian)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = minimum(x[:Date])
  upper = maximum(x[:Date])
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dvzr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  joined
end

function fit!(dvzr::DateValizer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # get complete dates and aggregate
  joined = fullaggregate!(dvzr,x)
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  medians = getMedian(grpby,joined)
  dvzr.args[:medians] = medians
  dvzr.model=dvzr.args
end

function transform!(dvzr::DateValizer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # get complete dates, aggregate, and get medians
  joined = fullaggregate!(dvzr,x)
  # copy medians
  medians = dvzr.args[:medians]
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  fn = Dict(Dates.Hour=>Dates.hour,
            Dates.Minute=>Dates.minute,
            Dates.Second=>Dates.second,
            Dates.Day => Dates.day,
            Dates.Month=>Dates.month)
  try
    joined[sym]=fn[grpby].(joined[:Date])
  catch
    error("unknown dateinterval")
  end
  # find indices of missing
  missingndx = findall(ismissing.(joined[:Value]))
  jmndx=joined[missingndx,sym] .+ 1 # get time period index of missing, convert 0 index time to 1 index
  missingvals::SubArray = @view joined[missingndx,:Value]
  missingvals .= medians[jmndx,:Value] # replace missing with median value
  sum(ismissing.(joined[:,:Value])) == 0 || error("Aggregation by time period failed to replace missings")
  joined[:,[:Date,:Value]]
end

### ====

# fill-in missings with nearest-neighbors median
mutable struct DateValNNer <: Transformer
  model
  args

  function DateValNNer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :missdirection => :symmetric, #:reverse, # or :forward or :symmetric
        :dateinterval => Dates.Hour(1),
        :nnsize => 1,
        :strict => true,
        :aggregator => :median
    )
    new(nothing,mergedict(default_args,args))
  end
end

function fit!(dnnr::DateValNNer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  aggr = dnnr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  dnnr.model=dnnr.args
end

function transform!(dnnr::DateValNNer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # make sure aggregator function exists
  aggr = dnnr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function pass is not recognized: ",aggr)
  # get the Statistics function
  aggfn = gAggDict[aggr]
  # pass the aggregator function to the generic aggregator function
  fn = aggregatorclskipmissing(aggfn)
  grpby = typeof(dnnr.args[:dateinterval])
  sym = Symbol(grpby)
  # aggregate by time period
  x[sym] = round.(x[:Date],grpby)
  aggr = by(x,sym,MeanValue = :Value=>fn)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = round(minimum(x[:Date]),grpby)
  upper = round(maximum(x[:Date]),grpby)
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dnnr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  missingcount = sum(ismissing.(joined[:Value]))
  dnnr.args[:missingcount] = missingcount
  res = transform_worker!(dnnr,joined)
  count=1
  if dnnr.args[:missdirection] == :symmetric
    while sum(ismissing.(res[:Value])) > 0
      res = transform_worker!(dnnr,res)
      count += 1
    end
  end
  dnnr.args[:loopcount] = count
  res
end

function transform_worker!(dnnr::DateValNNer,joinc::T) where {T<:DataFrame}
  joined = deepcopy(joinc)
  maxrow = size(joined)[1]

  # to fill-in with nearest neighbors
  nnsize::Int64 = dnnr.args[:nnsize]
  themissing = findall(ismissing.(joined[:Value]))
  # ==== symmetric nearest neighbor
  missingndx = DataFrame()
  if dnnr.args[:missdirection] == :symmetric
    missed = themissing |> reverse
    missingndx[:Missed] = missed
    # get lower:upper range
    missingndx[:neighbors] = map(missingndx[:Missed]) do m
      lower = (m-nnsize >= 1) ? (m-nnsize) : 1
      upper = (m+nnsize <= maxrow) ? m+nnsize : maxrow
      lower:upper
    end
  else
    # ===== reverse and forward
    missed = (dnnr.args[:missdirection] == :reverse) ? (themissing |> reverse) : themissing
    missingndx[:Missed] = missed
    # dealing with boundary exceptions, default to range until the maxrow
    missingndx[:neighbors] = (m->((m+1>=maxrow) || (m+nnsize>=maxrow)) ? (m+1:maxrow) : (m+1:m+nnsize)).(missingndx[:Missed]) # NN ranges
  end
  #joined[missingndx[:Missed],:Value] = (r -> skipmedian(joined[r,:Value])).(missingndx[:neighbors]) # iterate to each range
  missingvals::SubArray = @view joined[missingndx[:Missed],:Value] # get view of only missings
  missingvals .=  (r -> skipmedian(joined[r,:Value])).(missingndx[:neighbors]) # replace with nn medians
  dnnr.args[:strict] && (sum(ismissing.(joined[:,:Value])) == 0 || error("Nearest Neigbour algo failed to replace missings"))
  joined
end

mutable struct CSVDateValReader <: Transformer
    model
    args
    function CSVDateValReader(args=Dict())
        default_args = Dict(
            :filename => "",
            :dateformat => ""
        )
        new(nothing,mergedict(default_args,args))
    end
end
function fit!(csvrdr::CSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvrdr.args[:filename]
    fmt = csvrdr.args[:dateformat]
    (fname != "" && fmt != "") || error("missing filename or date format")
    model = csvrdr.args
end

function transform!(csvrdr::CSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvrdr.args[:filename]
    fmt = csvrdr.args[:dateformat]
    df = CSV.read(fname)
    ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
    rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
    df[:Date] = DateTime.(df[:Date],fmt)
    df
end

mutable struct CSVDateValWriter <: Transformer
    model
    args
    function CSVDateValWriter(args=Dict())
        default_args = Dict(
            :filename => "",
            :dateformat => ""
        )
        new(nothing,mergedict(default_args,args))
    end
end

function fit!(csvwtr::CSVDateValWriter,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvwtr.args[:filename]
    fmt = csvwtr.args[:dateformat]
    (fname != "" && fmt != "") || error("missing filename or date format")
    model = csvwtr.args
end

function transform!(csvwtr::CSVDateValWriter,x::T) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvwtr.args[:filename]
    fmt = csvwtr.args[:dateformat]
    df = deepcopy(x)
    ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
    rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
    eltype(df[:Date]) <: DateTime || error("Date format error")
    df |> CSV.write(fname)
end


text:
function fib_producer(n)
    a, b = (0, 1)
    for i = 1:n
        produce(b)
        a, b = (b, a + b)
    end
end

tsk1 = Task( () -> fib_producer(10) ) # Task (runnable) @0x0000000005696180

consume(tsk1) #> 1
consume(tsk1) #> 1
consume(tsk1) #> 2
consume(tsk1) #> 3
consume(tsk1) #> 5
consume(tsk1) #> 8
consume(tsk1) #> 13
consume(tsk1) #> 21
consume(tsk1) #> 34
consume(tsk1) #> 55
consume(tsk1) #> nothing # Task (done) @0x0000000005696180

for n in tsk1
    println(n)
end
#> 1 1 1 2 3 5 8 13 21 34 55
tsk1 = @task fib_producer(10) # Task (runnable) @0x0000000005696d80

ext: jl

max_stars_repo_path: Module 1/Chapter04/tasks.jl

max_issues_repo_path: Module 1/Chapter04/tasks.jl

max_forks_repo_path: Module 1/Chapter04/tasks.jl


text:
module FootballForm

using Turing, CSV, Plots, StatsPlots, OnlineStats, DataFrames, DataFramesMeta, Dates, Pipe, ProgressMeter
Turing.setadbackend(:tracker)

include("data.jl")
include("game.jl")
include("model.jl")
include("util.jl")

export run_form_analysis, log_all_plots, load_games, load_league

end

max_stars_repo_name: shalinkpatel/football-form

max_issues_repo_name: shalinkpatel/football-form

max_forks_repo_name: shalinkpatel/football-form


text:
"""
$(TYPEDEF)

Returns a variable containing parameters to affect the `newton` algorithm when solving `F(x) = 0`.

# Arguments (with default values):
$(TYPEDFIELDS)

# Arguments for line search (Armijo)
- `linesearch = false`: use line search algorithm (i.e. Newton with Armijo's rule)
- `α = 1.0`: alpha (damping) parameter for line search algorithm
- `αmin  = 0.001 `: minimal vslue of the damping `alpha`

!!! tip "Mutating"
    For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package `Setfield.jl` to drastically simplify the mutation of different fields. See the tutorials for examples.
"""
@with_kw struct NewtonPar{T, L <: AbstractLinearSolver, E <: AbstractEigenSolver}
	"absolute tolerance for `F(x)`"
	tol::T           = 1e-10
	"number of Newton iterations"
	maxIter::Int64 	 = 50
	"display Newton iterations?"
	verbose::Bool    = false
	"linear solver, must be `<: AbstractLinearSolver`"
	linsolver::L 	 = DefaultLS()
	"eigen solver, must be `<: AbstractEigenSolver`"
	eigsolver::E 	 = DefaultEig()
	linesearch::Bool = false
	α::T             = convert(typeof(tol), 1.0)        # damping
	αmin::T          = convert(typeof(tol), 0.001)      # minimal damping
end

####################################################################################################
"""
		newton(F, J, x0, p0, options::NewtonPar; normN = norm, callback = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true, kwargs...)

This is the Newton-Krylov Solver for `F(x, p0) = 0` with Jacobian w.r.t. `x` written `J(x, p0)` and initial guess `x0`. The function `normN` allows to specify a norm for the convergence criteria. It is important to set the linear solver `options.linsolver` properly depending on your problem. This linear solver is used to solve ``J(x, p_0)u = -F(x, p_0)`` in the Newton step. You can for example use `linsolver = DefaultLS()` which is the operator backslash: it works well for Sparse / Dense matrices. See [Linear solvers (LS)](@ref) for more informations.

# Arguments:
- `F` is a function with input arguments `(x, p)` returning a vector `r` that represents the functional and for type stability, the types of `x` and `r` should match. In particular, it is not **inplace**.
- `J` is the jacobian of `F` at `(x, p)`. It can assume two forms. Either `J` is a function and `J(x, p)` returns a `::AbstractMatrix`. In this case, the default arguments of `NewtonPar` will make `newton` work. Or `J` is a function and `J(x, p)` returns a function taking one argument `dx` and returns `dr` of the same type of `dx`. In our notation, `dr = J * dx`. In this case, the default parameters of `NewtonPar` will not work and you have to use a Matrix Free linear solver, for example `GMRESIterativeSolvers`.
- `x0` initial guess
- `p0` set of parameters to be passed to `F` and `J`
- `options::NewtonPar` variable holding the internal parameters used by the `newton` method
- `callback` function passed by the user which is called at the end of each iteration. The default one is the following `cbDefault(x, f, J, res, it, itlinear, options; k...) = true`. Can be used to update a preconditionner for example. The arguments passed to the callback are as follows
    - `x` current solution
    - `f` current residual
    - `J` current jacobian
    - `res` current norm of the residual
    - `iteration` current newton iteration
    - `itlinear` number of iterations to solve the linear system
    - `optionsN` a copy of the argument `options` passed to `newton`
    - `kwargs` kwargs arguments, contain your initial guess `x0`
- `kwargs` arguments passed to the callback. Useful when `newton` is called from `continuation`

# Output:
- solution
- history of residuals
- flag of convergence
- number of newton iterations
- total number of linear iterations

# Simplified calls
When `J` is not passed, the jacobian **matrix** is then computed with finite differences (beware of large systems of equations!). The call is as follows:

	newton(F, x0, p0, options::NewtonPar; kwargs...)

You can also pass functions which do not have parameters `x -> F(x)`, `x -> J(x)` as follows

	newton(F, J, x0, options::NewtonPar;  kwargs...)

or

	newton(F, x0, options::NewtonPar;  kwargs...)

# Example

```
julia> F(x, p) = x.^3 .- 1
julia> Jac(x, p) = spdiagm(0 => 3 .* x.^2) # sparse jacobian
julia> x0 = rand(1_000)
julia> opts = NewtonPar()
julia> sol, hist, flag, _ = newton(F, Jac, x0, nothing, opts, normN = x -> norm(x, Inf))
```

!!! tip "Other formulation"
    If you don't have parameters, you can still use `newton` as follows `newton((x,p) -> F(x), (x,p)-> J(x), x0, nothing, options)`

!!! warning "Linear solver"
    Make sure that the linear solver (Matrix-Free...) corresponds to your jacobian (Matrix-Free vs. Matrix based).
"""
function newton(Fhandle, Jhandle, x0, p0, options::NewtonPar; normN = norm, callback = cbDefault, kwargs...)
	# Extract parameters
	@unpack tol, maxIter, verbose, α, αmin, linesearch = options

	# Initialize iterations
	x = _copy(x0)
	f = Fhandle(x, p0)
	d = _copy(f)

	res = normN(f)
	resHist = [res]

	# iterations count
	it = 0

	# total number of linear iterations
	itlineartot = 0

	# Displaying results
	verbose && displayIteration(it, res)

	# invoke callback before algo really starts
	compute = callback(x, f, nothing, res, it, 0, options; x0 = x0, resHist = resHist, fromNewton = true, kwargs...)

	# Main loop
	while (res > tol) & (it < maxIter) & compute
		J = Jhandle(x, p0)
		d, _, itlinear = options.linsolver(J, f)
		itlineartot += sum(itlinear)

		# Update solution: x .= x .- d
		minus!(x, d)

		f = Fhandle(x, p0)
		res = normN(f)

		push!(resHist, res)
		it += 1

		verbose && displayIteration(it, res, itlinear)

		compute = callback(x, f, J, res, it, itlinear, options; x0 = x0, resHist = resHist, fromNewton = true, kwargs...)
	end
	((resHist[end] > tol) && verbose) && @error("\n--> Newton algorithm failed to converge, residual = $(res[end])")
	flag = (resHist[end] < tol) & callback(x, f, nothing, res, it, nothing, options; x0 = x0, resHist = resHist, fromNewton = true, kwargs...)
	verbose && displayIteration(0, res, 0, true) # display last line of the table
	return x, resHist, flag, it, itlineartot
end

# simplified call to newton when no Jacobian is passed in which case we estimate it using finiteDifferences
function newton(Fhandle, x0, p0, options::NewtonPar; kwargs...)
	Jhandle = (u, p) -> finiteDifferences(z -> Fhandle(z, p), u)
	return newton(Fhandle, Jhandle, x0, p0, options; kwargs...)
end

# default callback
cbDefault(x, f, J, res, it, itlinear, options; k...) = true


text:
@testset "Compensator" begin

    bg = 0.5
    kappa = 0.5
    kernD = Distributions.Exponential(1/0.5)
    events = sort(rand(10))

    testval = HawkesProcesses.compensator(events[1], events, bg, kappa, kernD)

    @test length(testval) == 1

    inferTests && @inferred HawkesProcesses.compensator(events[1], events, bg, kappa, kernD)

end


text:
# ================================================================
#     AbstractVertexKernel
# ================================================================

abstract type AbstractVertexKernel <: KernelFunctions.Kernel end

function (kernel::AbstractVertexKernel)(gv1::Tuple{AbstractGraph, Integer}, gv2::Tuple{AbstractGraph, Integer})

    (g1, v1) = gv1
    (g2, v2) = gv2
    return kernel(g1, v1, g2, v2)
end

# ================================================================
#     ConstVertexKernel
# ================================================================

struct ConstVertexKernel <: AbstractVertexKernel
    c::Float64
end

(kernel::ConstVertexKernel)(g1, v1, g2, v2) = kernel.c

# ================================================================
#     DiracVertexKernel
# ================================================================

struct DiracVertexKernel <: AbstractVertexKernel end

(kernel::DiracVertexKernel)(g1::AbstractValGraph, v1, g2::AbstractValGraph, v2) =
    (get_vertexval(g1, v1, :) == get_vertexval(g2, v2, :)) ? 1.0 : 0.0

# ================================================================
#     DotVertexKernel
# ================================================================

struct DotVertexKernel <: AbstractVertexKernel end

(kernel::DotVertexKernel)(g1::AbstractValGraph, v1, g2::AbstractValGraph, v2) =
    Float64(dot(get_vertexval(g1, v1, :), get_vertexval(g2, v2, :)))


text:
#!/usr/bin/env julia

using EzXMLPatched

# Convert a Julia expression to an XML document.
function expr2xml(expr)
    doc = XMLDocument()
    setroot!(doc, expr2elem(expr))
    return doc
end

function expr2elem(expr)
    if isa(expr, Expr)
        # <expr/>
        expr_node = ElementNode("expr")
        # <head/>
        head_node = ElementNode("head")
        link!(head_node, TextNode(repr(expr.head)))
        link!(expr_node, head_node)
        # <args/>
        args_node = ElementNode("args")
        for arg in expr.args
            link!(args_node, expr2elem(arg))
        end
        link!(expr_node, args_node)
        return expr_node
    else
        # <literal/>
        literal_node = ElementNode("literal")
        literal_node["type"] = string(typeof(expr))
        link!(literal_node, TextNode(repr(expr)))
        return literal_node
    end
end

prettyprint(expr2xml(Meta.parse(String(read(stdin)))))


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Type shorthands #
# ========================== #
const _Variable = VariableRef
const _Decision = DecisionRef
const _Scalar = Union{_Variable, _Decision, _Constant}
const _DecisionAffOrQuadExpr{C} = Union{DecisionAffExpr{C}, DecisionQuadExpr{C}}

# Basic operations #
# ========================== #
MA.mutability(::Type{<:_DecisionAffOrQuadExpr}) = MA.IsMutable()
function MA.mutable_copy(expr::_DecisionAffOrQuadExpr)
    return map_coefficients(MA.copy_if_mutable, expr)
end

# Promote operation #
# ========================== #
# _Constant--_Decision
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               ::Type{<:_Constant}, S::Type{<:_Decision})
    return DAE
end
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               S::Type{<:_Decision}, ::Type{<:_Constant})
    return DAE
end
# _Variable--Decision
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               ::Type{<:_Variable}, S::Type{<:_Decision})
    return DAE
end
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               S::Type{<:_Decision}, ::Type{<:_Variable})
    return DAE
end
# _Constant--_DecisionAffOrQuadExpr
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               ::Type{<:_Constant}, S::Type{<:_DecisionAffOrQuadExpr})
    return S
end
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               S::Type{<:_DecisionAffOrQuadExpr}, ::Type{<:_Constant})
    return S
end
# _Variable/_Decision--_DecisionAffOrQuadExpr
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, ::Type{<:Union{_Variable, _Decision}},
    S::Type{<:_DecisionAffOrQuadExpr})
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, S::Type{<:_DecisionAffOrQuadExpr},
    ::Type{<:Union{_Variable, _Decision}})
    return S
end
# (_Variable/_Decision)AffExpr--(_Variable/_Decision)AffExpr
function MA.promote_operation(
    op::Union{typeof(+), typeof(-)}, ::Type{GenericAffExpr{C1, _Variable}},
    ::Type{GenericAffExpr{C2, _Decision}}) where {C1, C2}
    C = MA.promote_operation(op, C1, C2)
    return DecisionAffExpr{C}
end
function MA.promote_operation(
    op::Union{typeof(+), typeof(-)}, ::Type{GenericAffExpr{C1, _Decision}},
    ::Type{GenericAffExpr{C2, _Variable}}) where {C1, C2}
    C = MA.promote_operation(op, C1, C2)
    return DecisionAffExpr{C}
end
# (_Variable/_Decision)AffExpr--_DecisionAffOrQuadExpr
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, ::Type{GenericAffExpr{C,V}},
    S::Type{_DecisionAffOrQuadExpr{C}}) where {C, V <: Union{_Variable, _Decision}}
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, S::Type{_DecisionAffOrQuadExpr{C}},
    ::Type{GenericAffExpr{C,V}}) where {C, V <: Union{_Variable, _Decision}}
    return S
end
# _DecisionAffOrQuadExpr--_DecisionAffOrQuadExpr
function MA.promote_operation(::Union{typeof(+), typeof(-)}, ::Type{A},
                               ::Type{A}) where {A <: _DecisionAffOrQuadExpr}
    return A
end
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               ::Type{<:DecisionAffExpr},
                               S::Type{<:DecisionQuadExpr})
    return S
end
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               S::Type{<:DecisionQuadExpr},
                               ::Type{<:DecisionAffExpr})
    return S
end
# _Variable--DecisionRef
function MA.promote_operation(::typeof(*), ::Type{<:_Variable}, ::Type{<:_Decision})
    return DQE
end
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{<:_Variable})
    return DQE
end
# _Decision--_Decision
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{<:_Decision})
    return DQE
end
# _Variable/_Decision--(_Variable/_Decision)AffExpr
function MA.promote_operation(::typeof(*), ::Type{<:_Variable}, ::Type{GenericAffExpr{C, _Decision}}) where {C, D <: _Decision}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{C, D}}, ::Type{<:_Variable}) where {C, D <: _Decision}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{GenericAffExpr{C, V}}) where {C, V <: _Variable}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{C, V}}, ::Type{<:_Decision}) where {C, V <: _Variable}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{D}, ::Type{GenericAffExpr{C, D}}) where {C, D <: _Decision}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{C, D}}, ::Type{D}) where {C, D <: _Decision}
    return DecisionQuadExpr{C}
end
# _Variable/_Decision--DecisionAffExpr
function MA.promote_operation(::typeof(*), ::Type{V}, ::Type{DecisionAffExpr{C}}) where {C, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::typeof(*), ::Type{DecisionAffExpr{C}}, ::Type{V}) where {C, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{C}
end

function MA.scaling(aff::DecisionAffExpr{C}) where C
    if !isempty(aff.variables.terms) || !isempty(aff.decisions.terms)
        throw(InexactError("Cannot convert `$aff` to `$C`."))
    end
    return MA.scaling(aff.variables.constant)
end
function MA.scaling(quad::DecisionQuadExpr{C}) where C
    if !isempty(quad.variables.terms) || !isempty(quad.decisions.terms) || !isempty(quad.decisions.quad.terms) || !isempty(quad.cross_terms)
        throw(InexactError("Cannot convert `$quad` to `$C`."))
    end
    return MA.scaling(quad.variables.aff)
end

# Mutable operate #
# ========================== #
# zero/one #
function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, aff::DecisionAffExpr)
    MA.mutable_operate!(op, aff.variables)
    MA.mutable_operate!(op, aff.decisions)
    return aff
end
function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, quad::DecisionQuadExpr)
    MA.mutable_operate!(op, quad.variables)
    MA.mutable_operate!(op, quad.decisions)
    empty!(quad.cross_terms)
    return quad
end
# * #
function MA.mutable_operate!(::typeof(*), expr::_DecisionAffOrQuadExpr, α::_Constant)
    if iszero(α)
        return MA.mutable_operate!(zero, expr)
    else
        return map_coefficients_inplace!(x -> MA.mul!(x, α), expr)
    end
end
# +/- #
function MA.mutable_operate!(::typeof(+), expr::_DecisionAffOrQuadExpr, x)
    return JuMP.add_to_expression!(expr, x)
end
function MA.mutable_operate!(::typeof(-), expr::_DecisionAffOrQuadExpr, x)
    return JuMP.add_to_expression!(expr, -1.0, x)
end
# add/sub_mul #
function MA.mutable_operate!(::typeof(MA.add_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar)
    return JuMP.add_to_expression!(expr, x)
end
function MA.mutable_operate!(::typeof(MA.add_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar, y::_Scalar)
    return JuMP.add_to_expression!(expr, x, y)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar)
    return JuMP.add_to_expression!(expr, -1.0, x)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar, y::_Scalar)
    return JuMP.add_to_expression!(expr, -x, y)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::Union{_Variable, _Decision}, y::_Constant)
    return JuMP.add_to_expression!(expr, x, -y)
end
function MA.mutable_operate!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, x, y)
    return MA.mutable_operate!(op, expr, x * y)
end
@generated function _add_sub_mul_reorder!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, args::Vararg{Any, N}) where N
    n = length(args)
    @assert n ≥ 3
    varidx = findall(t -> (t <: Union{_Variable, _Decision}), collect(args))
    allscalar = all(t -> (t <: _Constant), args[setdiff(1:n, varidx)])
    idx = (allscalar && length(varidx) == 1) ? varidx[1] : n
    coef = Expr(:call, :*, [:(args[$i]) for i in setdiff(1:n, idx)]...)
    return :(MA.mutable_operate!(op, expr, $coef, args[$idx]))
end
function MA.mutable_operate!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, x, y, z, other_args::Vararg{Any, N}) where N
    return _add_sub_mul_reorder!(op, expr, x, y, z, other_args...)
end
function JuMP.add_to_expression!(expr::_DecisionAffOrQuadExpr, α::_Constant, β::_Constant)
    return JuMP.add_to_expression!(expr, *(α, β))
end


text:
module ChainRulesOverloadGeneration

using ChainRulesCore

export on_new_rule, refresh_rules

include("ruleset_loading.jl")
include("precompile.jl")

function __init__()
    # Need to refresh rules when a package is loaded
    push!(Base.package_callbacks, _package_hook)
end

end # module


text:
# Benchmarks for vectorized operations.

module LazyAlgebraVectorBenchmarks

using BenchmarkTools
using LazyAlgebra

const Vn = LazyAlgebra

#------------------------------------------------------------------------------
module V1
using ArrayTools
using LazyAlgebra:
    Floats, vscale!, vzero!,
    promote_multiplier, arguments_have_incompatible_axes

function vupdate!(y::AbstractArray{<:Floats,N},
                  α::Number,
                  x::AbstractArray{<:Floats,N}) where {N}
    I = all_indices(x, y)
    if α == 1
        @inbounds @simd for i in I
            y[i] += x[i]
        end
    elseif α == -1
        @inbounds @simd for i in I
            y[i] -= x[i]
        end
    elseif α != 0
        alpha = promote_multiplier(α, x)
        @inbounds @simd for i in I
            y[i] += alpha*x[i]
        end
    end
    return y
end

function vcombine!(dst::AbstractArray{<:Floats,N},
                   α::Number,
                   x::AbstractArray{<:Floats,N},
                   β::Number,
                   y::AbstractArray{<:Floats,N}) where {N}
    if α == 0
        axes(x) == axes(dst) || arguments_have_incompatible_axes()
        vscale!(dst, β, y)
    elseif β == 0
        axes(y) == axes(dst) || arguments_have_incompatible_axes()
        vscale!(dst, α, x)
    else
        I = all_indices(dst, x, y)
        if α == 1
            if β == 1
                @inbounds @simd for i in I
                    dst[i] = x[i] + y[i]
                end
            elseif β == -1
                @inbounds @simd for i in I
                    dst[i] = x[i] - y[i]
                end
            else
                beta = promote_multiplier(β, y)
                @inbounds @simd for i in I
                    dst[i] = x[i] + beta*y[i]
                end
            end
        elseif α == -1
            if β == 1
                @inbounds @simd for i in I
                    dst[i] = y[i] - x[i]
                end
            elseif β == -1
                @inbounds @simd for i in I
                    dst[i] = -x[i] - y[i]
                end
            else
                beta = promote_multiplier(β, y)
                @inbounds @simd for i in I
                    dst[i] = beta*y[i] - x[i]
                end
            end
        else
            alpha = promote_multiplier(α, x)
            if β == 1
                @inbounds @simd for i in I
                    dst[i] = alpha*x[i] + y[i]
                end
            elseif β == -1
                @inbounds @simd for i in I
                    dst[i] = alpha*x[i] - y[i]
                end
            else
                beta = promote_multiplier(β, y)
                @inbounds @simd for i in I
                    dst[i] = alpha*x[i] + beta*y[i]
                end
            end
        end
    end
    return dst
end

end # module V1

#------------------------------------------------------------------------------

n = 10_000
T = Float32
x = rand(T,n)
y = rand(T,n)
z1 = similar(x)
zn = similar(x)
alphas = (0,1,-1,-4.0)
betas = (0,1,-1,+2.0)
v1_vupdate!(dst, y, alpha, x) = V1.vupdate!(vcopy!(dst, y), alpha, x)
vn_vupdate!(dst, y, alpha, x) = LazyAlgebra.vupdate!(vcopy!(dst, y), alpha, x)
for a in alphas
    println("\nvupdate!(y, α=$a, x):")
    print("  v1:  ")
    @btime v1_vupdate!($z1,$y,$a,$x);
    print("  new: ")
    @btime vn_vupdate!($zn,$y,$a,$x);
    dz = vnorm2(z1 - zn)
    printstyled("  -> ‖z1 - zn‖ = $dz\n"; color=(dz==0 ? :green : :red))
end

for a in alphas, b in betas
    println("\nvcombine!(z, α=$a, x, β=$b, y):")
    print("  v1:  ")
    @btime V1.vcombine!($z1,$a,$x,$b,$y);
    print("  new: ")
    @btime LazyAlgebra.vcombine!($zn,$a,$x,$b,$y);
    dz = vnorm2(z1 - zn)
    printstyled("  -> ‖z1 - zn‖ = $dz\n"; color=(dz==0 ? :green : :red))
end

end # module


text:
"""script for running overnight to collect data"""

cd(dirname(@__FILE__))
include("../Percolation.jl")
using Plots, ColorSchemes, LaTeXStrings, DelimitedFiles, BenchmarkTools
figpath = "../../fig/percolation/"

for len in (10, 20, 40, 80, 160)
    samples, runs = 50, 1000

    probs = collect(range(0, 1, length=samples))
    gyration_avg, gyration_σ = zeros(samples), zeros(samples)

    for run in 1:runs
        randgrid = rand(len, len)
        for sample in 1:samples
            grid = randgrid .<= probs[sample]
            gyration = Percolation.clustgyration(grid)
            gyration_avg[sample] += gyration
            gyration_σ[sample] += gyration ^ 2
        end
    end

    gyration_avg ./= runs
    # absolute value added to ensure no negative numbers resulting from floating-point error at small σ
    gyration_σ = .√(abs.(gyration_σ ./ runs - gyration_avg .^ 2))

    # save data
    open("correlation-length-data/gyration-full-$len.txt", "w") do io
        writedlm(io, [probs gyration_avg gyration_σ])
    end

    p = scatter(probs, gyration_avg, yerr=gyration_σ, markerstrokecolor=:red,
    color=:mediumblue, lw=2, legend=false, fontfamily="Computer Modern",
    title=L"L=%$len,\ \textrm{averaged\ over\ %$runs \ runs}", xlabel=L"p", ylabel=L"\xi",
    titlefontsize=18, tickfontsize=10, labelfontsize=18)
    plot!(p, probs, gyration_avg, color=:mediumblue, lw=2)

    savefig(p, figpath * "gyration-full-$len.pdf")

    # use full range data to find narrow range for finding p_c and xi_max with better precision
    maxgyr = argmax(gyration_avg)
    gyrstep = probs[2] - probs[1]
    minprob, maxprob = probs[maxgyr - 1] - gyrstep / 2, probs[maxgyr + 1] + gyrstep / 2

    samples, runs = 50, 10000

    probs = collect(range(minprob, maxprob, length=samples))
    gyration_avg, gyration_σ = zeros(samples), zeros(samples)

    for run in 1:runs
        randgrid = rand(len, len)
        for sample in 1:samples
            grid = randgrid .<= probs[sample]
            gyration = Percolation.clustgyration(grid)
            gyration_avg[sample] += gyration
            gyration_σ[sample] += gyration ^ 2
        end
    end

    gyration_avg ./= runs
    # absolute value added to ensure no negative numbers resulting from floating-point error at small σ
    gyration_σ = .√(abs.(gyration_σ ./ runs - gyration_avg .^ 2))

    # save data
    open("correlation-length-data/gyration-zoom-$len.txt", "w") do io
        writedlm(io, [probs gyration_avg gyration_σ])
    end

    # find peak
    argpc = argmax(gyration_avg)
    pc = probs[argpc]
    ximax = gyration_avg[argpc]
    ximaxerr = gyration_σ[argpc]

    # save peak data
    open("correlation-length-data/gyration-peak.txt", "a") do io
        writedlm(io, [len pc ximax ximaxerr])
    end

    p = plot(probs, gyration_avg, ribbon=gyration_σ, fillcolor=:lightskyblue,
        color=:mediumblue, marker=:circle, lw=3, legend=false, fontfamily="Computer Modern",
        title=L"L=%$len,\ \textrm{averaged\ over\ %$runs \ runs}", xlabel=L"p", ylabel=L"\xi",
        titlefontsize=18, tickfontsize=10, labelfontsize=18)

    vline!(p, [pc], color=:red, linestyle=:dash)
    hline!(p, [ximax], color=:red, linestyle=:dash)

    annotetey = ylims(p)[1] + (ylims(p)[2] - ylims(p)[1]) / 15
    annotatex = pc + (xlims(p)[2] - xlims(p)[1]) / 100
    annotate!(p, annotatex, annotetey, text(L"p_c=%$(round(pc, digits=3))",
        :blueviolet, :left, :dash, 12))

    savefig(p, figpath * "gyration-zoom-$len.pdf")
end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
module ErrorHandler end

import Base.Test

custom_handler(r::Test.Success) = println("Teste '$(r.expr)' executado com sucesso")

custom_handler(r::Test.Failure) = error("Teste '$(r.expr)' falhou. $(r)")

custom_handler(r::Test.Error) = rethrow(r);

export custom_handler

max_stars_repo_name: raphaelpanta/julia-lang-exemplos

max_issues_repo_name: raphaelpanta/julia-lang-exemplos

max_forks_repo_name: raphaelpanta/julia-lang-exemplos


text:
using Test

using Algorithm101
const ALGO = Algorithm101

using DataStructure101
const DS = DataStructure101

@testset "TopologicalSort" begin

    # dict = Dict(
    #     1 => "Underwear",
    #     2 => "T-Shirt",
    #     3 => "Pents",
    #     4 => "Jasket",
    #     5 => "Socks",
    #     6 => "Shoes"
    # )

    g = DS.Graph(6, :di_matrix)

    DS.relate!(g, 1, 2)
    DS.relate!(g, 1, 3)

    DS.relate!(g, 2, 4)
    DS.relate!(g, 2, 6)

    DS.relate!(g, 3, 6)

    DS.relate!(g, 5, 6)

    @test ALGO.find_entrances(g) == [1, 5]
    @test ALGO.topological_sort(g) == [1, 5, 2, 3, 6, 4]

end


text:
using HTTP
using HTTP: hasheader
using MbedTLS

@static if is_apple()
    launch(x) = run(`open $x`)
elseif is_linux()
    launch(x) = run(`xdg-open $x`)
elseif is_windows()
    launch(x) = run(`cmd /C start $x`)
end

@async begin
    sleep(2)
    launch("https://127.0.0.1:8000/examples/mwe")
end

HTTP.listen("127.0.0.1", 8000;
            ssl = true,
            sslconfig = MbedTLS.SSLConfig(joinpath(dirname(@__FILE__), "cert.pem"),
                                          joinpath(dirname(@__FILE__), "key.pem"))) do http
    if HTTP.WebSockets.is_websocket_upgrade(http.message)

        HTTP.WebSockets.upgrade(http) do client
            count = 1
            while !eof(client);
                msg = String(readavailable(client))
                println(msg)
                write(client, "Hello JavaScript! From Julia $count")
                count += 1
            end
        end
    else
        HTTP.Servers.handle_request(http) do req::HTTP.Request
            HTTP.Response(200,read(joinpath(dirname(@__FILE__),"mwe.html"), String))
        end
    end
end


text:
module TestAbstractFuncs

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running abstract function tests")

## create instance
vals = rand(30, 4)
vals[3, 4] = 0.5
dats = Date[Date(2013, 7, ii) for ii=1:30]
nams = [:A, :B, :C, :D]

allTypes = (:Timedata, :Timenum, :Timematr)

###########################
## names, idx, get, core ##
###########################

## test information retrieval functions
for t in allTypes
    eval(quote
        td = $(t)(vals, nams, dats)
        @test isequal(names(td), nams)
        @test isequal(TimeData.idx(td), DataArray(dats))
        @test get(td, 3, 4) == 0.5
        @test get(td)[3, 4] == 0.5
        @test core(td)[3, 4] == 0.5
        @test size(td) == (30, 4)
        @test ndims(td) == 2
        @test td[:, 1:2] == hcat(td[1], td[2])
        #        @test td == hcat(td[1], td[2], td[3], td[4])
    end)
end

## get with :
##-----------

tm = TimeData.testcase(TimeData.Timematr, 1)
@test TimeData.get(tm[1, 1], :) == 100
@test_throws Exception TimeData.get(tm[1, :], :)

df = DataFrame(a = [3, 5])
@test get(df[1, :], :) == 3
@test_throws Exception TimeData.get(df, :)

#############
## isnaElw ##
#############

## set up Timenum with NAs
df = DataFrame()
df[:a] = @data([4, 5, 6, NA, 8])
df[:b] = @data([3, 8, NA, NA, 2])
dats = [Date(2014,1,1):Date(2014,1,5)]
tn = Timenum(df, dats)

## test outcome
naVals = [false false; false false; false true; true true; false false]
isNaTd = Timedata(naVals, names(tn), idx(tn))
@test isnaElw(tn) == isNaTd

###################
## isequal tests ##
###################

df = DataFrame()
df[:a] = @data([3, NA])
df[:b] = @data([4, NA])

df2 = DataFrame()
df2[:a] = @data([3, 4])
df2[:b] = @data([4, NA])

df3 = DataFrame()
df3[:a] = @data([10, 4])
df3[:b] = @data([4, NA])

tn = Timenum(df)
tn2 = Timenum(df2)
tn3 = Timenum(df3)
td = Timedata(df)
tn4 = Timenum(df, [Date(2010, 1, 1):Date(2010, 1, 2)])

## testing isequal
@test isequal(tn, tn)
@test !isequal(tn, tn2)
@test !isequal(tn, tn3)

## testing ==
@test isna(tn == tn)
@test isna(tn == tn2)
@test !(tn == tn3)

## differences in names, indices or types
@test !isequal(tn, td)
@test !(tn == td)
@test !isequal(tn, tn4)
@test !(tn == tn4)

## elementwise equal objects
@test isequalElw(tn, tn) == Timedata([true true; true true],
                                          names(td), idx(td))

## elementwise with unequal entries

## differences in names, indices or types throw errors
@test_throws ErrorException !isequalElw(tn, td)
@test_throws ErrorException !isequalElw(tn, tn4)

@test isequalElw(tn, tn2) == Timedata([true true; false true],
                                           names(td), idx(td))

##############
## isapprox ##
##############

tm = testcase(Timenum, 1)
tm2 = tm .+ 0.00000001
@test !isequal(tm, tm2)
@test isapprox(tm, tm2)

tm = testcase(Timenum, 2)
tm2 = tm .+ 0.00000001
@test !isequal(tm, tm2)
@test isapprox(tm, tm2)

td = Timedata(rand(2, 2))
td2 = deepcopy(td)
td[1, 1] = get(td, 1, 1) .+ 0.00000001
@test isapprox(td, td2)

####################
## complete_cases ##
####################

df = DataFrame()
df[:a] = @data([3, NA])
df[:b] = @data([4, NA])
td = Timedata(df)
tn = Timenum(df)
tm = Timematr(rand(2, 3))

@test complete_cases(td) == [true, false]
@test complete_cases(tn) == [true, false]
@test complete_cases(tm) == [true, true]

###############
## equMeta ##
###############

df = DataFrame()
df[:a] = @data([4, 5, 6, NA, 8])
df[:b] = @data([3, 8, NA, NA, 2])
dats = [Date(2014,1,1):Date(2014,1,5)]
tn = TimeData.Timenum(df, dats)

## similar
@test TimeData.equMeta(tn, tn)
@test TimeData.similarMeta(tn, tn)

## different indices
tn2 = TimeData.Timenum(df)
@test !TimeData.equMeta(tn, tn2)
@test TimeData.similarMeta(tn, tn2)

## different names
df = DataFrame()
df[:a] = @data([4, 5, 6, NA, 8])
df[:b] = @data([3, 8, NA, NA, 2])
dats = [Date(2014,1,1):Date(2014,1,5)]
tn2 = TimeData.Timenum(df, dats)
names!(tn2.vals, [:c, :d])
@test !TimeData.equMeta(tn, tn2)
@test TimeData.similarMeta(tn, tn2)

## different types
tm = TimeData.Timedata(df, dats)
@test !TimeData.equMeta(tn, tm)
@test !TimeData.similarMeta(tn, tm)

end


text:
#push!(LOAD_PATH,"../src/")
using Documenter, BridgeSDEInference

makedocs(
    modules = [BridgeSDEInference],
    format = Documenter.HTML(
        #prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename="BridgeSDEInference.jl",
    pages = [
        "Home" => "index.md",
        "Tutorial" => Any[
            "Workflow" => joinpath("man", "overview.md"),
            "Defining a diffusion" => joinpath("man", "model_definition.md"),
            "Data generation" => joinpath("man", "generate_data.md"),
            "Setup" => joinpath("man", "setup.md"),
            "Running the sampler" => joinpath("man", "run.md"),
            "Querying the results" => joinpath("man", "query.md"),
            ],
        "Examples" => Any[
            "FitzHugh-Nagumo model" => joinpath("examples", "fitzhugh_nagumo.md"),
            "Jansen-Rit model" => joinpath("examples", "jansen_rit.md"),
            "Lorenz63 system" => joinpath("examples", "lorenz63.md"),
            "Lorenz96 system" => joinpath("examples", "lorenz96.md"),
            "Prokaryotic autoregulatory gene network" => joinpath("examples", "prokaryote.md"),
            "Sine diffusion" => joinpath("examples", "sine.md")
        ],
        "Visualisation tools" => joinpath("vis", "visualisation.md")
    ],
)

deploydocs(
    repo = "github.com/mmider/BridgeSDEInference.jl.git",
)


text:
# Put all the graphical tests here.

@testset "periodize" begin
    c = Cylinder([0,0,0], 1/4, 1)  # circular base
    # c = PolygonalPrism([0,0,0], [1/4 1/4; -1/4 1/4; -1/4 -1/4; 1/4 -1/4], 1)  # square base
    # c = Prism([0,0,0], regpoly(5, 1/4), 1)  # regular pentagonal base
    ∆range = Cuboid([0,0,0], [8,2*√3,1])
    A = [1 0 0; 0.5 √3/2 0; 0 0 1]'
    c_array = periodize(c, A, ∆range)

    b = Cuboid([0,0,0], [12,4*√3,1])
    bnd = bounds(b)
    x = range(bnd[1][1], stop=bnd[2][1], length=120*5)
    y = range(bnd[1][2], stop=bnd[2][2], length=70*5)
    Nx, Ny = length(x), length(y)
    X = repeat(x, outer=(1,Ny))
    Y = repeat(y', outer=(Nx,1))
    V = zeros(Bool, Nx, Ny)

    kd = KDTree(c_array)
    for j = 1:Ny, i = 1:Nx
        p = [x[i], y[j], 0.0]
        s = findfirst(p, kd)
        V[i,j] = s≠nothing
    end

    using PyPlot
    clf()
    axis("equal")
    pcolor(X, Y, V, cmap="gray_r")

    using PyCall
    @pyimport matplotlib.patches as patch
    b∆range = bounds(∆range)
    cnr = b∆range[1][1:2]
    w, h = -((-)(b∆range...)[1:2])
    gca()[:add_artist](patch.Rectangle(cnr, w, h, fill=false, edgecolor="r", linestyle="--"))
end  # @testset "periodize"


text:
# This file includes the testset for pins

@testset "PinTestSet" begin
    @info "Running PinTestSet ..."

    # Construction of Outpin
    op = Outpin()
    @test isa(op.links, Missing)
    @test !isbound(op)

    # Construction of Inpin
    ip = Inpin()
    @test !isbound(ip)

    # Connection of pins
    op, ip = Outpin(), Inpin()
    l = connect!(op, ip)
    @test isa(l, Link)
    @test isbound(op)
    @test isbound(ip)
    @test l.masterid == op.id
    @test l.slaveid == ip.id
    @test isconnected(op, ip)
    op2, ip2 = Outpin(), Inpin()
    @test_throws Exception connect!(op, op2)  # Outpin cannot be connected to Inpin
    @test_throws Exception connect!(ip, ip2)  # Inpin cannot be connected to Inpin
    @test_throws Exception connect!(ip, op)   # Inpin cannot be connected to Outpin

    op1 = Outpin()
    op2 = Outpin()
    ip = Inpin()
    @test_throws MethodError connect!(ip, op1)  # Inpin cannot drive and Outpin
    @test_throws MethodError connect!(op1, op2) # Outpoin cannot drive and Outpin
    @test !isbound(op1)
    @test !isbound(op2)
    @test !isbound(ip)
    l = connect!(op1, ip)   # Outpin can drive Inpin
    @test isbound(op1)
    @test isbound(ip)
    @test_throws ErrorException connect!(op1, ip)   # Reconnection is not possible
    @test_throws ErrorException connect!(op2, ip)   # `ip` is bound. No new connections are allowed.
    @test isconnected(op1, ip)
    disconnect!(op1, ip)
    @test !isconnected(op1, ip)
    l = connect!(op2, ip)
    @test isconnected(op2, ip)

    @test isbound(ip)
    # Connection of multiple Inpins to an Outpin
    op = Outpin()
    ips = [Inpin() for i in 1 : 5]
    ls = map(ip -> connect!(op, ip), ips)
    for (l, ip) in zip(ls, ips)
        @test l.masterid == op.id
        @test l.slaveid == ip.id
    end

    # Data transfer through pins
    op, ip = Outpin(), Inpin()
    @test_throws MethodError take!(op)  # Data cannot be read from Outpin
    @test_throws MethodError put!(ip, 1.)  # Data cannot be written into Inpin
    l = connect!(op, ip)
    t = @async while true
        take!(ip) === NaN && break  # Take from ip
    end
    for val in 1 : 5
        put!(op, val)   # Write into op.
        @test !istaskdone(t)
        @test read(l.buffer) == val
    end
    put!(op, NaN)
    @test istaskdone(t)

    # Disconnection of Outpin and Inpin
    op, ip = Outpin(), Inpin()
    l = connect!(op, ip)
    @test isconnected(op, ip)
    disconnect!(op, ip)
    @test !isconnected(op, ip)

    @info "Done PinTestSet."
end  # testset


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayAvailableSslPredefinedPolicies <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{ApplicationGatewaySslPredefinedPolicy} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function ApplicationGatewayAvailableSslPredefinedPolicies(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(ApplicationGatewayAvailableSslPredefinedPolicies, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(ApplicationGatewayAvailableSslPredefinedPolicies, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type ApplicationGatewayAvailableSslPredefinedPolicies

const _property_map_ApplicationGatewayAvailableSslPredefinedPolicies = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_ApplicationGatewayAvailableSslPredefinedPolicies = Dict{Symbol,String}(Symbol("value")=>"Vector{ApplicationGatewaySslPredefinedPolicy}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ ApplicationGatewayAvailableSslPredefinedPolicies }) = collect(keys(_property_map_ApplicationGatewayAvailableSslPredefinedPolicies))
Swagger.property_type(::Type{ ApplicationGatewayAvailableSslPredefinedPolicies }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayAvailableSslPredefinedPolicies[name]))}
Swagger.field_name(::Type{ ApplicationGatewayAvailableSslPredefinedPolicies }, property_name::Symbol) =  _property_map_ApplicationGatewayAvailableSslPredefinedPolicies[property_name]

function check_required(o::ApplicationGatewayAvailableSslPredefinedPolicies)
    true
end

function validate_property(::Type{ ApplicationGatewayAvailableSslPredefinedPolicies }, name::Symbol, val)
end


text:
SupportedEstimation = Union{StatsModels.TableRegressionModel,FixedEffectModel,GLFixedEffectModel}

function parse(est::SupportedEstimation; level::Real=0.95, drop_cons::Bool=true)
    if est isa FixedEffectModel
        coefplot = parse(StatsModels.coef(est), StatsModels.confint(est;level=level), StatsModels.coefnames(est))
    else
        coefplot = parse(StatsModels.coef(est), StatsModels.confint(est,level), StatsModels.coefnames(est))
    end
    if drop_cons
        drop_cons!(coefplot::Coefplot)
        equidist!(coefplot)
    end
    return coefplot
end

function parse(coefvec::Vector{T} where T<:Real, confint::Matrix{T} where T<:Real ,coefnames::Vector{T} where T<:Any)
    # create Coefplot from combining regmodel output and reglabel, regnote, xtitle, ytitle
    @assert size(coefvec,1) == size(confint,1) == size(coefnames,1) "Dimension doesn't match"
    dict = OrderedDict([Symbol(coefnames[i])=>SinglecoefPlot(coefvec[i], confint[i,:]..., coefnames[i], i) for i in 1:size(coefvec,1)])
    Coefplot(dict)
end

"""
    drop_cons!(coefplot[, cons])

Drop constant term from the coefplot.

# Examples
```julia-repl
julia> coefplot = parse(ols; drop_cons=false);
julia> drop_cons!(coefplot);
julia> plot(coefplot)
```
"""
drop_cons!(coefplot::Coefplot,cons::Symbol=Symbol("(Intercept)")) = delete!(coefplot,cons)

"""
    equidist!(coefplot)

Rearrange the location of each coefficients according to their index in the subplot vector.

# Examples
```julia-repl
julia> coefplot = parse(coefplot; drop_cons=false);
julia> equidist!(coefplot);
julia> plot(coefplot)
```
"""
function equidist!(coefplot::Coefplot;interval::Real=1,preserve_locorder::Bool=false)
    if preserve_locorder
        loc = [singlecoefplot.thiscoef_loc for singlecoefplot in values(coefplot.dict)]
        @assert ~all(ismissing.(loc)) "there are missing singlecoefplot location"
        newloc = invperm(sortperm(loc))
        for (thiscoef_newloc, singlecoefplot) in zip(newloc,values(coefplot.dict))
            singlecoefplot.thiscoef_loc = thiscoef_newloc * interval
        end
    else
        for (index, singlecoefplot) in enumerate(values(coefplot.dict))
            singlecoefplot.thiscoef_loc = index * interval
        end
    end
    return coefplot
end


text:
#nb # %% slideshow [markdown] {"slideshow": {"slide_type": "slide"}}
#md # # [Lecture 6 - Puzzle time! (Regular Expressions)](@id lecture6)
#!md # # Lecture 6 - Puzzle time! (Regular Expressions)
#
# **Essential Skills For Computational Biology**
#
# *2020-06-30*
#
# Kevin Bonham, PhD

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Strings contain patterns
#
# Here's a typical list of files from a DNA sequencing run:

sequences = [
    "C0005-3F_L001_R1_001.fastq.gz",
    "C0005-3F_L001_R2_001.fastq.gz",
    "C0005-3F_L002_R1_001.fastq.gz",
    "C0005-3F_L002_R2_001.fastq.gz",
    "C0016-3E_L001_R1_001.fastq.gz",
    "C0016-3E_L001_R2_001.fastq.gz",
    "C0016-3E_L002_R1_001.fastq.gz",
    "C0016-3E_L002_R2_001.fastq.gz",
]

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ##Simple Patterns:
#
# - `abc` matches literal "abc"
# - `[abc]` matches any one of `a`, `b`, or `c`
# - `[^abc]` matches any character _except_ `a`, `b`, or `c`
# - `[a-z0-9]` matches any lowercase letter or digit

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# Puzzle: https://regexr.com/57ijc

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ##Repeating Patterns:
#
# - `+` matches 1 or more of a pattern
#   - eg. `[0-9]+` matches one or more digits
# - `*` matches 0 or more of a pattern
# - `?` matches 0 or 1 of a pattern
# - `{N}` matches `N` of a pattern
#   - eg. `[abcdr]{4}` matches `abra` and `cada`, but not `abracadabra`
#   - Can also to ranges; `{2,4}` matches 2 to 4, `{3,}` is 3 or more

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# Puzzle: https://regexr.com/57ik4

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Special Patterns:
#
# - `.` matches any character, except a new line
# - `\w` matches any word character (same as `[A-Za-z0-9_]`)
# - `\d` matches any digit (same as `[0-9]`)
# - `\s` matches any whitespace (space, tab, newline, etc)
# - If you need to match a character that is otherwise special,
#   "escape" it with `\`
#   - eg. to match a literal period, `\.`

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# Puzzles:
# - https://regexr.com/57ikd
# - https://regexr.com/57ikg

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## More fun with regular Expressions
#
# - Use the regexr cheatsheet (or find another one - they're all over)
#   - There are _a lot_ more patterns that can be used
# - https://regexcrossword.com/
# - My most complicated regex to date:

# ```
# ^([`~]{3,})(?:(?:(?:\\{|\\{\\.|)(julia)(?:;|))|(@(docs|autodocs|ref|meta|index|content|example|repl|eval|setup|raw)))\\s*(?:.*?)(\\}|)\\s*$
# ```

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# > Some people, when confronted with a problem, think "I know, I'll use regular expressions." Now they have two problems.

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ![xkcd regex](https://imgs.xkcd.com/comics/regular_expressions.png)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Regular expressions in julia

for seq in sequences
    m = match(r"^([MC]\d+\-[0-9][EF])_L00.+fastq\.gz", seq)
    println(m.captures[1])
end


text:
"""
Calculate total energy components.
Potentials are NOT updated.
psi is assumed to be already orthogonalized properly
rho is calculated outside this function.
"""
function calc_Energies( LF, ∇2, Potentials::PotentialsT, rho, Focc, psi::Array{Float64,2} )
  #
  Nstates = size(psi)[2]
  ΔV  = LF.LFx.h * LF.LFy.h * LF.LFz.h
  #
  E_Hartree = 0.5*sum( rho .* Potentials.Hartree ) * ΔV
  #
  E_xc = sum(excVWN( rho ) .* rho) * ΔV
  #
  E_ionic = sum( rho .* Potentials.Ionic ) * ΔV

  #
  E_kin = 0.0
  for ist = 1:Nstates
    nabla2_v = ∇2 * psi[:,ist]
    E_kin = E_kin + -0.5*dot( psi[:,ist], nabla2_v ) * ΔV * Focc[ist]
  end

  Energies = EnergiesT(0.0, 0.0, 0.0, 0.0, 0.0)
  Energies.Kinetic = E_kin
  Energies.Ionic   = E_ionic
  Energies.Hartree = E_Hartree
  Energies.XC      = E_xc
  Energies.Total   = E_kin + E_ionic + E_Hartree + E_xc
  #
  return Energies
end

"""
Calculate total energy components AND update Potentials.
psi is assumed to be already orthogonalized properly
This function modifies `Energies` and `Potentials`.
"""
function calc_Energies!( LF, Gv, Energies, Potentials, Focc, psi::Array{Float64,2} )

  Npoints = size(psi)[1]
  Nstates = size(psi)[2]
  ΔV  = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # Electron density
  rho = calc_rho( Focc, psi )
  #
  # solve Poisson equation to get Hartree potential
  V_Hartree = solve_poisson_FFT( Gv, rho )
  # Calculate Hartree energy
  E_Hartree = 0.5*sum( rho .* V_Hartree ) * ΔV
  #
  V_xc = excVWN( rho ) + rho .* excpVWN( rho )
  E_xc = sum(excVWN( rho ) .* rho) * ΔV
  #
  E_ionic = sum( rho .* Potentials.Ionic ) * ΔV

  #
  E_kin = 0.0
  for ist = 1:Nstates
    nabla2_v = op_nabla2( LF, psi[:,ist] )
    E_kin = E_kin + -0.5*dot( psi[:,ist], nabla2_v ) * ΔV * Focc[ist]
  end

  #
  Energies.Kinetic = E_kin
  Energies.Ionic   = E_ionic
  Energies.Hartree = E_Hartree
  Energies.XC      = E_xc
  Energies.Total   = E_kin + E_ionic + E_Hartree + E_xc
  #
  Potentials.Hartree = V_Hartree[:]
  Potentials.XC      = V_xc[:]

  return
end

function calc_Energies!( LF, Gv, ∇2, Energies, Potentials, Focc, psi::Array{Float64,2} )

  Npoints = size(psi)[1]
  Nstates = size(psi)[2]
  ΔV  = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # Electron density
  rho = calc_rho( Focc, psi )
  #
  # solve Poisson equation to get Hartree potential
  V_Hartree = solve_poisson_FFT( Gv, rho )
  # Calculate Hartree energy
  E_Hartree = 0.5*sum( rho .* V_Hartree ) * ΔV
  #
  V_xc = excVWN( rho ) + rho .* excpVWN( rho )
  E_xc = sum(excVWN( rho ) .* rho) * ΔV
  #
  E_ionic = sum( rho .* Potentials.Ionic ) * ΔV

  #
  E_kin = 0.0
  for ist = 1:Nstates
    nabla2_v = ∇2 * psi[:,ist]
    E_kin = E_kin + -0.5*dot( psi[:,ist], nabla2_v ) * ΔV * Focc[ist]
  end

  #
  Energies.Kinetic = E_kin
  Energies.Ionic   = E_ionic
  Energies.Hartree = E_Hartree
  Energies.XC      = E_xc
  Energies.Total   = E_kin + E_ionic + E_Hartree + E_xc
  #
  Potentials.Hartree = V_Hartree[:]
  Potentials.XC      = V_xc[:]

  return
end


text:
using MultiStartLH
using Dates, Test

mdl = MultiStartLH;

isdir(test_dir())  ||  mkpath(test_dir());

@testset "MultiStartLH.jl" begin
    include("history_test.jl");
    include("update_rule_test.jl");
    include("optim_test.jl");
end


text:
### auxiliary methods
###

## return index which are not NAN
function  isnotnan(arr)
    index = map(!isnan, arr[:])

    return(index)
end

###return blacklist if any
function read_blacklist(blackname)
    if isfile(blackname)
        df= CSV.read(blackname, delim=";")
        blacklist= df.votname
    else
        blacklist= [""]
    end

    return(blacklist)
end

## transform struct to DF
function convertStruct2Df(s)::DataFrame
    T= typeof(s)
    f= fieldnames(T)
    n=length(f)

     field1= fieldname(T,1) ; val1= getfield(s, field1)
     df= DataFrame(field1 => val1)

     for i in 2:n
         field= fieldname(T,i)
         val= getfield(s, field)
         insertcols!(df,i,field=>val)
     end

     return(df)
end

## special printing code
## return string
function specialstr(str, CODE)
    d= Dict( "PURPLE" => "\033[95m" ,
             "CYAN" => "\033[96m" ,
             "DARKCYAN" => "\033[36m" ,
             "BLUE" => "\033[94m" ,
             "GREEN" => "\033[92m" ,
             "YELLOW" => "\033[93m" ,
             "RED" => "\033[91m" ,
             "BOLD" => "\033[1m" ,
             "UNDERLINE" => "\033[4m" ,
             "END" => "\033[0m")

    if haskey(d, CODE)
        res= @sprintf("%s%s%s", d[CODE], str, d["END"])
        return(res)
    else
        return(str)
    end
end

function bold(s) GaiaClustering.specialstr(s,"BOLD") end
function yellow(s) GaiaClustering.specialstr(s,"YELLOW") end
function purple(s) GaiaClustering.specialstr(s,"PURPLE") end
function cyan(s) GaiaClustering.specialstr(s,"CYAN") end
function blue(s) GaiaClustering.specialstr(s,"BLUE") end
function red(s) GaiaClustering.specialstr(s,"RED") end

function header_extract()
    println("#")
    println("# >>>>>>>>>> EXTRActing Stellar Clusters from GAIA data <<<<<<<<<<")
    println("# Version: $VERSION  (https://github.com/bosscha/gaia-shock)")
    println("#\n#")
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
using Test
using GradientBoost
using GradientBoost.GB
using GradientBoost.LossFunctions

struct DummyGradientBoost <: GBAlgorithm; end

struct StubGradientBoost <: GBAlgorithm
  loss_function::LossFunction
  sampling_rate::Float64
  learning_rate::Float64
  num_iterations::Int
end

function GradientBoost.GB.build_base_func(
  gb::StubGradientBoost,
  instances,
  labels,
  prev_func_pred,
  psuedo)

  function pred(instances)
    num_instances = size(instances, 1)
    predictions = Array{Float64}(undef, num_instances)
    for i = 1:num_instances
      predictions[i] = sum(instances[i,:])
    end
    predictions
  end

  model_const = 0.5
  return (instances) -> model_const .* pred(instances)
end

sgb_instances = [
  2 2;
  2 4
]
sgb_labels = [
  1.0;
  3.0
]

@testset "Gradient Boost" begin
  @testset "not implemented functions throw an error" begin
    dgb = DummyGradientBoost()
    emp_mat = fill(1,1,1)
    emp_vec = Array[]
    @test_throws ErrorException build_base_func(
      dgb,
      emp_mat,
      emp_vec,
      emp_vec,
      emp_vec
    )
  end

  @testset "stochastic_gradient_boost works" begin
    # Sanity check
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    @test 1 == 1
  end

  @testset "fit returns model" begin
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    @test model isa GBModel
  end

  @testset "predict works" begin
    expected = [
      3.0
      3.5
    ]
    sgb = StubGradientBoost(LeastSquares(), 1.0, 0.5, 1)
    model = stochastic_gradient_boost(sgb, sgb_instances, sgb_labels)
    predictions = predict(model, sgb_instances)
    @test predictions == expected
  end

  @testset "create_sample_indices works" begin
    instances = [1:5 6:10]
    labels = [1:5]

    sgb = StubGradientBoost(LeastSquares(), 1, 1, 1)
    indices = create_sample_indices(sgb, instances, labels)
    @test length(indices) == 5
    @test length(unique(indices)) == 5
    @test minimum(indices) >= 1
    @test maximum(indices) <= 5

    sgb = StubGradientBoost(LeastSquares(), 0.5, 1, 1)
    indices = create_sample_indices(sgb, instances, labels)
    @test length(indices) == 3
    @test length(unique(indices)) == 3
    @test minimum(indices) >= 1
    @test maximum(indices) <= 5
  end
end


text:
# nufft/z-list.jl

include("dtft.jl")
include("nufft.jl")


text:
using Formatting, FFTW, Jets, JetPackWaveFD, LinearAlgebra, SpecialFunctions, Statistics, Test, WaveFD

function make_op(modeltype, interpmethod, fs; v₀=1500, ϵ₀=0.2, η₀=0.4, comptype = Float32, st = 0.0, wavelet = WaveletCausalRicker(f=5.0))
    nsponge = 10
    pad = 20
    nx = 50+2*pad
    ny = 40+2*pad
    nz = 30+2*pad
    nt = 101
    dx,dy,dz,dt = 25.0,25.0,25.0,0.002
    xmin,ymin,zmin,tmin = -pad*dx,-pad*dy,-pad*dz,0.0
    xmax,ymax,zmax,tmax = xmin+dx*(nx-pad-1),ymin+dy*(ny-pad-1),zmin+dz*(nz-pad-1),tmin+dt*(nt-1)
    sx = dx*div(nx,2)
    sy = dy*div(ny,2)
    sz = dz
    rz = [dz for iy = 1+pad:ny-pad, ix = 1+pad:nx-pad][:]
    ry = [(iy-1)*dy for iy = 1+pad:ny-pad, ix=1+pad:nx-pad][:]
    rx = [(ix-1)*dx for iy = 1+pad:ny-pad, ix=1+pad:nx-pad][:]

    local ϵ,η
    if modeltype  == "v"
        ϵ = zeros(Float32,nz,ny,nx)
        η = zeros(Float32,nz,ny,nx)
    else
        ϵ = Array{Float32}(undef,0,0,0)
        η = Array{Float32}(undef,0,0,0)
    end

    b = ones(Float32,nz,ny,nx)

    F = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(b = b, f = 0.85, ϵ = ϵ, η = η,
        z0 = zmin, y0 = ymin, x0 = xmin, dz = dz, dy = dy, dx = dx,
        sz = sz, sy = sy, sx = sx, st = 0.0, rz = rz, ry = ry, rx = rx,
        dtrec = dt, dtmod = dt, ntrec = nt, nsponge = nsponge, comptype = comptype, compscale = 1e-4,
        freqQ = 5.0, qMin = 0.1, qInterior = 100.0, wavelet = wavelet, freesurface = fs,
        reportinterval = 0, interpmethod = interpmethod, isinterior = false)

    m₀ = zeros(domain(F))
    m₀[:,:,:,modelindex(F,"v")] .= v₀
    if modeltype == "vϵη"
        m₀[:,:,:,modelindex(F,"ϵ")] .= ϵ₀
        m₀[:,:,:,modelindex(F,"η")] .= η₀
    end

    m₀, F
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- close" begin
    m₀, F = make_op("vϵη", :hicks, false)
    d = F * m₀

    @test isfile(state(F).srcfieldfile*"-pspace")
    @test isfile(state(F).srcfieldfile*"-mspace")
    close(F)
    @test !(isfile(state(F).srcfieldfile*"-pspace"))
    @test !(isfile(state(F).srcfieldfile*"-mspace"))
    d = F * m₀
    @test isfile(state(F).srcfieldfile*"-pspace")
    @test isfile(state(F).srcfieldfile*"-mspace")
    close(F)
    @test !(isfile(state(F).srcfieldfile*"-pspace"))
    @test !(isfile(state(F).srcfieldfile*"-mspace"))

    close(F)
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- perfstat" begin
    m₀, F = make_op("vϵη", :hicks, false)
    d = F * m₀

    s = Jets.perfstat(F)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test s["%imaging"] ≈ 0.0

    J = jacobian!(F, m₀)
    d = J * m₀
    s = Jets.perfstat(J)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test isa(s["%imaging"], Float64)

    m = J' * d
    s = Jets.perfstat(J)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test isa(s["%imaging"], Float64)

    close(F)
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- non-causal source injection" begin
    wavelet = WaveletRicker(;f=5.0)
    mₒ, F₁ = make_op(:hicks, "v", false; wavelet, st=-1.0)
    mₒ, F₂ = make_op(:hicks, "v", false; wavelet, st=-2.0)

    d₁ = F₁ * mₒ
    d₂ = F₂ * mₒ
    @test d₁ ≈ d₂
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- linearization, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    d = F * m₀

    mmask = zeros(domain(F))
    mmask[:,:,:,modelindex(F,"v")] .= 1
    if modeltype == "vϵη"
        mmask[:,:,:,modelindex(F,"ϵ")] .= 0.001
        mmask[:,:,:,modelindex(F,"η")] .= 0.001
    end

    μobs, μexp = linearization_test(F, m₀,
        mmask = mmask,
        μ = 100*sqrt.([1.0,1.0/2,1.0/4,1.0/8,1.0/16,1.0/32,1.0/64,1.0/128,1.0/256,1.0/512]))

    @show μobs, μexp
    @show minimum(abs, μobs - μexp)
    @test minimum(abs, μobs - μexp) < 0.01

    close(F)
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- dot product, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear, :hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    J = jacobian!(F, m₀)
    m = zeros(domain(J))
    m[:,:,:,modelindex(F,"v")] .= 100 * (-1 .+ 2*rand(Float32,size(m)[1:3]))
    if modeltype == "vϵη"
        m[:,:,:,modelindex(F,"ϵ")] .= rand(Float32,size(m)[1:3]...)
        m[:,:,:,modelindex(F,"η")] .= rand(Float32,size(m)[1:3]...)
    end
    d = -1 .+ 2*rand(range(J))

    lhs, rhs = dot_product_test(J, m, d)
    err = abs(lhs - rhs) / abs(lhs + rhs)
    @show lhs,rhs,(lhs-rhs)/(lhs+rhs)
    @test lhs ≈ rhs rtol=1e-4
    @test err < 1e-4

    close(F)
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- linearity, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    J = jacobian!(F, m₀)
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs

    close(F)
end

# note the compression is exercised on the second pass of F * m₀
@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- serialization, modeltype=$modeltype, C=$C" for modeltype in ("vϵη", "v"), C in (Float32,UInt32)
    m₀, F = make_op(modeltype, :hicks, false, comptype = C)
    d₁ = F * m₀
    d₂ = F * m₀
    @test d₁ ≈ d₂

    close(F)
end

@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- srcillum, C=$C, modeltype=$modeltype" for C in (Float32, UInt32), modeltype in ("vϵη", "v")
    m₀, F = make_op(:hicks, modeltype, false, comptype = C)
    d₁ = F * m₀;
    J = jacobian(F, m₀);
    s1 = srcillum(F, m₀);
    s2 = srcillum(J);
    close(F)
    @test s1 ≈ s2
    @test maximum(s1) > 0
end

#=
@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks)
    z,y,x,dz,dy,dx,dtrec,dtmod,tmax,sz,sy,sx,c = 2000.0,500.0,2500.0,40.0,40.0,40.0,0.016,0.004,2.0,1000.0,250.0,20.0,1500.0
    nz,ny,nx,nt = round(Int,z/dz)+1,round(Int,y/dy)+1,round(Int,x/dx)+1,round(Int,tmax/dtrec)+1

    fmax = c/6/max(dx,dy,dz)
    fdef = WaveletCausalRicker(f=fmax/2)

    function analytic_result()
        # non-negative frequency axis
        nt_pad = nextprod([2,3,5,7],4*nt)
        nw = div(nt_pad,2)+1
        w = range(0, stop=pi/dtrec, length=nw)

        # wavelet
        f = WaveFD.get(fdef, dtrec*collect(0:nt_pad-1))
        F = rfft(f)

        # Morse and Feshbach, 1953, p. 891
        u_ana = zeros(Float32,nt,nz,ny,nx)
        U = zeros(Complex{Float64},nw)
        for ix = 1:nx, iy = 1:ny, iz = 1:nz
            z = (iz-1)*dz
            y = (iy-1)*dy
            x = (ix-1)*dx
            r = sqrt((z-sz)^2+(y-sy)^2+(x-sx)^2)
            for iw = 2:nw
                if r*w[iw] > eps(Float32)
                    U[iw] = exp(-im*w[iw]/c*r)*F[iw]
                end
            end
            u_ana[:,iz,iy,ix] = irfft(U,nt_pad)[1:nt]
            printfmt("{:.2f}", ((ix-1)*ny*nz+(iy-1)*nz+iz-1)/(nz*ny*nx)*100))
        end
        u_ana
    end
    @info "computing analytic result (be patient, can take several minutes)..."
    u_ana = analytic_result()
    @info "...done."

    local ϵ,η
    if modeltype == "v"
        ϵ = zeros(Float32,nz,ny,nx)
        η = zeros(Float32,nz,ny,nx)
    else
        ϵ = Float32[]
        η = Float32[]
    end

    rx = (dx*[0:nx-1;] * ones(ny)')[:]
    ry = (ones(nx) * (dy*[0:ny-1;])')[:]
    rz = 2*dz*ones(length(rx))

    # Finite difference
    M = JopNlProp3DAcoVTIDenQ_DEO2_FDTD(
        comptype = UInt32,
        compscale = 1e-3,
        srcfieldfile = "$_scratch/srcfield-$(randstring()).bin",
        b = (1.0f0/2.0f0)*ones(Float32,nz,ny,nx),
        f = 0.75,
        ϵ = ϵ,
        η = η,
        sz = sz,
        sy = sy,
        sx = sx,
        rz = rz,
        ry = ry,
        rx = rx,
        interpmethod = interpmethod,
        dz = dz,
        dy = dy,
        dx = dx,
        dtrec = dtrec,
        dtmod = dtmod,
        ntrec = nt,
        wavelet = fdef,
        nsponge = 200,
        freqQ = 5.0,
        qMin = 0.1,
        qInterior = 100.0,
        freesurface = false,
        reportinterval = 250)

    m = zeros(domain(M))
    m[:,:,:,modelindex(M,"v")] .= c
    @info "computing numerical result (be patient)"
    d_tmp = M*m
    @info "...done."

    # Read FD field back from disk
    u_fd = zeros(Float32,state(M).ntrec,nz,ny,nx)
    u_snap = zeros(Float32,nz,ny,nx)
    u_snap_ginsu = zeros(Float32,size(state(M).ginsu, interior=false)...)
    io = open("$(state(M).srcfieldfile)-pold")

    open(state(M).compressor["pold"])
    for i = 1:state(M).ntrec
        WaveFD.compressedread!(io, state(M).compressor["pold"], i, u_snap_ginsu)
        super!(u_snap,state(M).ginsu,u_snap_ginsu,interior=false)
        u_fd[i,:,:,:] = u_snap[:,:,:]
    end
    close(state(M).compressor["pold"])

    # compute Pearson correlation (https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient)
    xx = ones(nz,ny,nx)
    for ix = 1:nx, iy = 1:ny, iz = 1:nz
        xx[iz,iy,ix] = cor(u_fd[:,iz,iy,ix],u_ana[:,iz,iy,ix])
    end
    xx_mean = mean(xx)
    @test xx_mean > .9

    close(M)
    @info "be patient, the garbage collection is probably deleting large scratch-files"
end
=#

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
function aux_sample!(Ω, lik::AbstractLikelihood, y, f)
    return aux_sample!(GLOBAL_RNG, Ω, lik, y, f)
end

function aux_sample!(
    rng::AbstractRNG, Ω, lik::AbstractLikelihood, y::AbstractVector, f::AbstractVector
)
    map!(Ω, y, f) do yᵢ, fᵢ
        ntrand(rng, aux_full_conditional(lik, yᵢ, fᵢ))
    end
    return Ω
end

function aux_sample(lik::AbstractLikelihood, y, f)
    return aux_sample(GLOBAL_RNG, lik, y, f)
end

function aux_sample(rng::AbstractRNG, lik::AbstractLikelihood, y, f)
    return aux_sample!(rng, init_aux_variables(lik, length(y)), lik, y, f)
end

function aux_posterior(lik::AbstractLikelihood, y, f)
    return aux_posterior!(init_aux_posterior(lik, length(y)), lik, y, f)
end

function aux_full_conditional(lik::AbstractLikelihood, y, f)
    return For(1:length(y)) do i
        aux_full_conditional(lik, y[i], f[i])
    end
end

function init_aux_variables(lik::AbstractLikelihood, n::Int)
    return init_aux_variables(GLOBAL_RNG, lik, n)
end

function init_aux_posterior(lik::AbstractLikelihood, n::Int)
    return init_aux_posterior(Float64, lik, n)
end

function logtilt(
    lik::AbstractLikelihood, Ω::TupleVector, y::AbstractVector, f::AbstractVector
)
    return mapreduce(+, aux_field(lik, Ω), y, f) do Ωᵢ, yᵢ, fᵢ
        logtilt(lik, Ωᵢ, yᵢ, fᵢ)
    end
end

function aug_loglik(lik::AbstractLikelihood, Ω, y, f)
    return logtilt(lik, Ω, y, f) + logdensity_def(aux_prior(lik, y), Ω)
end

function aux_kldivergence(lik::AbstractLikelihood, qΩ::For, y)
    return aux_kldivergence(lik, qΩ, aux_prior(lik, y))
end

function aux_kldivergence(::AbstractLikelihood, qΩ::For, pΩ::For)
    return mapreduce(Distributions.kldivergence, +, marginals(qΩ), marginals(pΩ))
end

function auglik_potential_and_precision(lik::AbstractLikelihood, Ω, y)
    return (auglik_potential(lik, Ω, y), auglik_precision(lik, Ω, y))
end

function expected_auglik_potential_and_precision(lik::AbstractLikelihood, qΩ, y)
    return (expected_auglik_potential(lik, qΩ, y), expected_auglik_precision(lik, qΩ, y))
end

# Generic wrapper for prior not taking any argument
aux_prior(lik::AbstractLikelihood, ::Real) = aux_prior(lik)

nlatent(::AbstractLikelihood) = 1 # Default number of latent for each likelihood
# This should potentially move to GPLikelihoods.jl


text:
"""
Linearly interpolates two vector x and y on a linear grid consisting of `nsteps`.
"""
function interp(x::Vector, y::Vector, nsteps::Int)
    # Interpolate
    itp = interpolate((x,), y, Gridded(Linear()))
    # Interpolate at the given resolution
    x_fills = LinRange(minimum(x), maximum(x), nsteps)
    y_fills = itp(x_fills)

    return collect(x_fills), y_fills
end

export interp


text:
# Todo:
#
# Abstract integrands and allow for other f-Divergences to work with the implemented methods
# Create more elaborate method for default domain of integration: Check that both means are inside cube and then choose widths of HyperCube proportional to Σ?
# Add other distributions with known analytical expressions for KL-divergence: Wishart, Beta, Gompertz, generalized Gamma

"""
    KullbackLeibler(p::Function,q::Function,Domain::HyperCube=HyperCube([-15,15]); tol=2e-15, N::Int=Int(3e7), Carlo::Bool=(length(Domain)!=1))
Computes the Kullback-Leibler divergence between two probability distributions `p` and `q` over the `Domain`.
If `Carlo=true`, this is done using a Monte Carlo Simulation with `N` samples.
If the `Domain` is one-dimensional, the calculation is performed without Monte Carlo to a tolerance of ≈ `tol`.
```math
D_{\\text{KL}}[p,q] \\coloneqq \\int \\mathrm{d}^m y \\, p(y) \\, \\mathrm{ln} \\bigg( \\frac{p(y)}{q(y)} \\bigg)
```
"""
function KullbackLeibler(p::Function, q::Function, Domain::HyperCube; tol::Real=1e-9, N::Int=Int(3e7), Carlo::Bool=false)
    function Integrand(x)
        P = p(x)[1];   Q = q(x)[1];   Rat = P / Q
        (Rat ≤ 0. || !isfinite(Rat)) && throw(ArgumentError("Ratio p(x)/q(x) = $Rat in log(p/q) for x=$x."))
        P * log(Rat)
    end
    if !Carlo
        return IntegrateND(Integrand, Domain; tol=tol)
    else
        return length(Domain) == 1 ? MonteCarloArea(x->Integrand(x[1]),Domain,N)[1] : MonteCarloArea(Integrand,Domain,N)
    end
end

function KullbackLeibler(p::Product{Continuous}, q::DiagNormal, Domain::HyperCube=HyperCube([[-20,20] for i in 1:length(p)]); tol::Real=1e-12, kwargs...)
    !(length(p) == length(q) == length(Domain)) && throw("KL: Sampling dimension mismatch: dim(p) = $(length(p)), dim(q) = $(length(q)), Domain = $(length(Domain))")
    sum(KullbackLeibler(p.v[i], Normal(q.μ[i],sqrt(q.Σ.diag[i])), HyperCube([Domain.L[i], Domain.U[i]]); tol=tol) for i in 1:length(p))
end
function KullbackLeibler(p::DiagNormal, q::Product{Continuous}, Domain::HyperCube=HyperCube([[-20,20] for i in 1:length(p)]); tol::Real=1e-12, kwargs...)
    !(length(p) == length(q) == length(Domain)) && throw("KL: Sampling dimension mismatch: dim(p) = $(length(p)), dim(q) = $(length(q)), Domain = $(length(Domain))")
    sum(KullbackLeibler(Normal(p.μ[i],sqrt(p.Σ.diag[i])), q.v[i], HyperCube([Domain.L[i], Domain.U[i]]); tol=tol) for i in 1:length(p))
end
function KullbackLeibler(p::Product{Continuous}, q::Product{Continuous}, Domain::HyperCube=HyperCube([[-20,20] for i in 1:length(p)]); tol::Real=1e-12, kwargs...)
    !(length(p) == length(q) == length(Domain)) && throw("KL: Sampling dimension mismatch: dim(p) = $(length(p)), dim(q) = $(length(q)), Domain = $(length(Domain))")
    sum(KullbackLeibler(p.v[i], q.v[i], HyperCube([Domain.L[i], Domain.U[i]]); tol=tol) for i in 1:length(p))
end

function KullbackLeibler(p::Distribution, q::Distribution, Domain::HyperCube=HyperCube([[-20,20] for i in 1:length(p)]); tol::Real=1e-9, N::Int=Int(3e7), Carlo::Bool=false)
    !(length(p) == length(q) == length(Domain)) && throw("KL: Sampling dimension mismatch: dim(p) = $(length(p)), dim(q) = $(length(q)), Domain = $(length(Domain))")
    function Integrand(x)
        P = logpdf(p,x);   Q = logpdf(q,x)
        exp(P) * (P - Q)
    end
    if !Carlo
        return IntegrateND(Integrand, Domain; tol=tol)
    else
        return length(p) == 1 ? MonteCarloArea(x->Integrand(x[1]),Domain,N)[1] : MonteCarloArea(Integrand,Domain,N)
    end
end

# Known analytic expressions
function KullbackLeibler(P::MvNormal, Q::MvNormal, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...)
    (1/2) * (logdet(Q.Σ) - logdet(P.Σ) - length(P.μ) + tr(inv(Q.Σ) * P.Σ) + InnerProduct(inv(Q.Σ), Q.μ-P.μ))
end
KullbackLeibler(P::Normal, Q::Normal, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...) = log(Q.σ / P.σ) + (1/2) * ((P.σ / Q.σ)^2 + (P.μ - Q.μ)^2 * Q.σ^(-2) -1.)
KullbackLeibler(P::Cauchy, Q::Cauchy, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...) = log(((P.σ + Q.σ)^2 + (P.μ - Q.μ)^2) / (4P.σ * Q.σ))

# Note the sign difference between the conventions (1/θ)*exp(-x/θ) and λ*exp(-λx). Distributions.jl uses the former.
KullbackLeibler(P::Exponential, Q::Exponential, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...) = log(Q.θ / P.θ) + P.θ / Q.θ - 1.
function KullbackLeibler(P::Weibull, Q::Weibull, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...)
    log(P.α / (P.θ^P.α)) - log(Q.α / (Q.θ^Q.α)) + (P.α - Q.α)*(log(P.θ) - Base.MathConstants.γ / P.α) + (P.θ / Q.θ)^Q.α * SpecialFunctions.gamma(1 + Q.α / P.α) - 1.
end
function KullbackLeibler(P::Distributions.Gamma, Q::Distributions.Gamma, Domain::HyperCube=HyperCube([-Inf,Inf]); kwargs...)
    (P.α - Q.α) * digamma(P.α) - loggamma(P.α) + loggamma(Q.α) + Q.α*log(Q.θ / P.θ) + P.α*(P.θ / Q.θ - 1.)
end

# """
#     NormalDist(DM::DataModel,p::AbstractVector) -> Distribution
# Constructs either `Normal` or `MvNormal` type from `Distributions.jl` using data and a parameter configuration.
# This makes the assumption, that the errors associated with the data are normal.
# """
# function NormalDist(DM::DataModel,p::AbstractVector)::Distribution
#     if length(ydata(DM)[1]) == 1
#         length(ydata(DM)) == 1 && return Normal(ydata(DM)[1] .- DM.model(xdata(DM)[1],p),ysigma(DM)[1])
#         return MvNormal(ydata(DM) .- map(x->DM.model(x,p),xdata(DM)),diagm(floatify(ysigma(DM).^2)))
#     else
#         throw("Not programmed yet.")
#     end
# end

# """
#     KullbackLeibler(DM::DataModel,p::AbstractVector,q::AbstractVector)
# Calculates Kullback-Leibler divergence under the assumption of a normal likelihood.
# """
# KullbackLeibler(DM::DataModel,p::AbstractVector,q::AbstractVector) = KullbackLeibler(NormalDist(DM,p),NormalDist(DM,q))
#
# KullbackLeibler(DM::DataModel,p::AbstractVector) = KullbackLeibler(MvNormal(zeros(length(ydata(DM))),inv(yInvCov(DM))),NormalDist(DM,p))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Matrix factorizations and decompositions

abstract type Factorization{T} end

eltype(::Type{<:Factorization{T}}) where {T} = T
size(F::Adjoint{<:Any,<:Factorization}) = reverse(size(parent(F)))
size(F::Transpose{<:Any,<:Factorization}) = reverse(size(parent(F)))

macro assertposdef(A, info)
   :($(esc(info)) == 0 ? $(esc(A)) : throw(PosDefException($(esc(info)))))
end

macro assertnonsingular(A, info)
   :($(esc(info)) == 0 ? $(esc(A)) : throw(SingularException($(esc(info)))))
end

"""
    issuccess(F::Factorization)

Test that a factorization of a matrix succeeded.

```jldoctest
julia> F = cholesky([1 0; 0 1]);

julia> LinearAlgebra.issuccess(F)
true

julia> F = lu([1 0; 0 0]);

julia> LinearAlgebra.issuccess(F)
false
```
"""
issuccess(F::Factorization)

function logdet(F::Factorization)
    d, s = logabsdet(F)
    return d + log(s)
end

function det(F::Factorization)
    d, s = logabsdet(F)
    return exp(d)*s
end

convert(::Type{T}, f::T) where {T<:Factorization} = f
convert(::Type{T}, f::Factorization) where {T<:Factorization} = T(f)

convert(::Type{T}, f::Factorization) where {T<:AbstractArray} = T(f)

### General promotion rules
Factorization{T}(F::Factorization{T}) where {T} = F
inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n)))

Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash, h, 1:nfields(F))
Base.:(==)(  F::T, G::T) where {T<:Factorization} = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F))
Base.isequal(F::T, G::T) where {T<:Factorization} = all(f -> isequal(getfield(F, f), getfield(G, f)), 1:nfields(F))::Bool

# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
    c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2))
    x = ldiv!(F, c2r)
    return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))), _ret_size(F, B))
end

function \(F::Factorization, B::AbstractVecOrMat)
    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    BB = similar(B, TFB, size(B))
    copyto!(BB, B)
    ldiv!(F, BB)
end
function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
    F = adjF.parent
    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    BB = similar(B, TFB, size(B))
    copyto!(BB, B)
    ldiv!(adjoint(F), BB)
end

# support the same 3-arg idiom as in our other in-place A_*_B functions:
function ldiv!(Y::AbstractVecOrMat, A::Factorization, B::AbstractVecOrMat)
    m, n = size(A, 1), size(A, 2)
    if m > n
        ldiv!(A, B)
        return copyto!(Y, view(B, 1:n, :))
    else
        return ldiv!(A, copyto!(Y, view(B, 1:m, :)))
    end
end
function ldiv!(Y::AbstractVecOrMat, adjA::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
    checksquare(adjA)
    return ldiv!(adjA, copyto!(Y, B))
end
function ldiv!(Y::AbstractVecOrMat, transA::Transpose{<:Any,<:Factorization}, B::AbstractVecOrMat)
    checksquare(transA)
    return ldiv!(transA, copyto!(Y, B))
end

# fallback methods for transposed solves
\(F::Transpose{<:Any,<:Factorization{<:Real}}, B::AbstractVecOrMat) = adjoint(F.parent) \ B
\(F::Transpose{<:Any,<:Factorization}, B::AbstractVecOrMat) = conj.(adjoint(F.parent) \ conj.(B))

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
push!(LOAD_PATH, ".")

using ISCLib

using HDF5, JLD, NIfTI, Logging
Logging.configure(level=DEBUG)

function process!(v::Array{Float64,1})
  map!(x-> isnan(x)?0:x, v ,v) # no nan`s
  map!(x-> x<0.1? 0: x, v, v) # no insignificant values-
  v
end

"""
rb - how many rows should be dropped from the back
"""
function run_isc(stim::Array{String,1}=["S1","S2"],pref::String="S",rb::Int=1)
  const means = [ get_mean(subj) for subj in 1:NSUBJ ]
  debug("means done")
  covs = zeros(Float32, SIZE,div(NSUBJ*(NSUBJ-1),2))
  disps = zeros(Float32, (SIZE,NSUBJ))
  indexes = [get_indexes(i,stim,rb) for i in 1:NSUBJ]
  for j in 1:size(indexes[1],1)
    local seq = [get_niis(i,indexes[i][j,1]:indexes[i][j,2]) for i in 1:NSUBJ]
    debug("Got niis")
    cov_n!(covs,disps,seq, means, 16) # 16 - 6 sec before 2 sec after => 9sec = 16 "frames"
    debug("Done with : $j")
  end
  debug("covs_n done")
  h5open("out/$pref|covs.h5","w") do file
    write(file, "data", covs)
  end
  h5open("out/$pref|disps.h5","w") do file
    write(file, "data", disps/(16*length(indexes)))
  end
  niscores = ISC_res(covs,disps)
  info("Scores done")

  res = process!(niscores)
  ni = get_nii(1,1)
  res = NIfTI.NIVolume(ni.header, ni.extensions, reshape(res, SHAPE)) # Cheap, as reshape returns a view
  niwrite("out/$pref|segmentation_result.nii",res)
  niscores
end

run_isc()
debug("S done!")
run_isc(["V1","V2","V3","V4"],"V")
debug("V done!")
run_isc(["S1","S2","V1","V2","V3","V4"],"S_V",2)
debug("S and V done!")


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

@testset ExtendedTestSet "`ProjectiveTransforms`" begin
    @testset ExtendedTestSet "image" begin
        transform = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)

        ## We run `ProjectiveTransforms` in the different [`Context`]s:
        imagetrain = FastAI.run(transform, Training(), image)
        @test size(imagetrain) == (32, 32)

        imagevalid = FastAI.run(transform, Validation(), image)
        @test size(imagevalid) == (32, 32)

        imageinference = FastAI.run(transform, Inference(), image)
        @test size(imageinference) == (32, 48)

        ## During inference, the aspect ratio should stay the same
        @test size(image, 1) / size(image, 2) == size(imageinference, 1) / size(imageinference, 2)
    end

    @testset ExtendedTestSet "keypoints" begin
        transform = ProjectiveTransforms((32, 32))
        ks = [SVector(0., 0), SVector(64, 96)]
        keypoints = Keypoints(ks, (64, 96))
        kstrain = FastAI.run(transform, Training(), keypoints)
        ksvalid = FastAI.run(transform, Validation(), keypoints)
        ksinference = FastAI.run(transform, Inference(), keypoints)

        @test ksvalid[1][1] == 0
        @test ksvalid[2][1] == 32
        @test ksinference[2] == ks[2] ./ 2
    end

    @testset ExtendedTestSet "image and keypoints" begin
        transform = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)
        ks = [SVector(0., 0), SVector(64, 96)]

        @test_nowarn FastAI.run(transform, Training(), (image, ks))
        @test_nowarn FastAI.run(transform, Validation(), (image, ks))
        @test_nowarn FastAI.run(transform, Inference(), (image, ks))
    end

    @testset ExtendedTestSet "`run!`" begin
        transform = ProjectiveTransforms((32, 32))
        image1 = rand(RGB, 64, 96)
        image2 = rand(RGB, 64, 96)
        buf = FastAI.run(transform, Validation(), image1)
        cbuf = copy(buf)
        FastAI.run!(buf, transform, Validation(), image1)
        # run and run! should return the same result
        @test buf ≈ cbuf

        FastAI.run!(buf, transform, Validation(), image2)
        # buf should be modified on different input
        @test !(buf ≈ cbuf)

    end
end


text:
using ForwardDiff

function anomaloushall(kx,ky,T,μ,Ham!,H1,H2,H3)
    ϵ = -1:0.1:1
    kb = 8.6*10^(-5) #eV / K
    β = 1/(kb*T) # ≈ 40 at 300 K
    ħ = 6.582*10^(-16) # eV*s
    Ham!(kx,ky,H1) ; Hx = similar(H1) ; Hy = similar(H1) ;
    val , vec = eigen!(H1)

    curv = similar(val)
    Dx(kx,ky,Ham!,H2,H3,Hx)
    vx = vec'*Hx*vec

    denominator = [ ifelse( n-m ≈ 0, 0, 1/(n-m)^2 ) for n in val , m in val]

    Dy(kx,ky,Ham!,H2,H3,Hy)
    vy = vec'*Hy*vec

    @. vy *= denominator

    @inbounds for c in eachindex(curv)
        curv[c] = -2*imag(sum(vx[c,m]*vy[m,c] for m in 1:length(val)))
    end

    @. curv /= T*ħ    # make curv larger, for curv is small

    nernst = 0
    @inbounds for n in 1:length(val)
        f = exp(β*(val[n] - μ)) |> x -> ifelse(x≈0, 0, x)
        nernst += (val[n]-μ)/( 1+f )+log( 1+1.0/f)/β |>  x -> x*curv[n]
    end
    -nernst
end