text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(MatLang, Symbol("##freqspaceM#28")) && precompile(Tuple{getfield(MatLang, Symbol("##freqspaceM#28")), Int64, typeof(identity), Int64, Symbol})
    isdefined(MatLang, Symbol("##onesM#11")) && precompile(Tuple{getfield(MatLang, Symbol("##onesM#11")), Array{Float64, 1}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("##onesM#11")) && precompile(Tuple{getfield(MatLang, Symbol("##onesM#11")), Array{Float64, 2}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("##randM#17")) && precompile(Tuple{getfield(MatLang, Symbol("##randM#17")), Array{Float64, 1}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("##randM#17")) && precompile(Tuple{getfield(MatLang, Symbol("##randM#17")), Array{Float64, 2}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("##sortM#29")) && precompile(Tuple{getfield(MatLang, Symbol("##sortM#29")), Bool, typeof(identity), Array{Int64, 1}, Symbol})
    isdefined(MatLang, Symbol("##sortM#30")) && precompile(Tuple{getfield(MatLang, Symbol("##sortM#30")), Bool, typeof(identity), Array{Int64, 2}, Int64, Symbol})
    isdefined(MatLang, Symbol("##sortM#30")) && precompile(Tuple{getfield(MatLang, Symbol("##sortM#30")), Bool, typeof(identity), Array{Integer, 3}, Int64, Symbol})
    isdefined(MatLang, Symbol("##sortperm#31")) && precompile(Tuple{getfield(MatLang, Symbol("##sortperm#31")), Int64, Bool, typeof(identity), Array{Int64, 2}})
    isdefined(MatLang, Symbol("##zerosM#5")) && precompile(Tuple{getfield(MatLang, Symbol("##zerosM#5")), Array{Float64, 1}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("##zerosM#5")) && precompile(Tuple{getfield(MatLang, Symbol("##zerosM#5")), Array{Float64, 2}, typeof(identity), Array{Int64, 1}})
    isdefined(MatLang, Symbol("#kw##eyeM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##eyeM")), NamedTuple{(:like,), Tuple{Array{Int8, 2}}}, typeof(MatLang.eyeM), Int64, Int64})
    isdefined(MatLang, Symbol("#kw##freqspaceM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##freqspaceM")), NamedTuple{(:dim,), Tuple{Int64}}, typeof(MatLang.freqspaceM), Int64, Symbol})
    isdefined(MatLang, Symbol("#kw##freqspaceM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##freqspaceM")), NamedTuple{(:dim,), Tuple{Int64}}, typeof(MatLang.freqspaceM), Int64})
    isdefined(MatLang, Symbol("#kw##ndgridM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##ndgridM")), NamedTuple{(:dim,), Tuple{Int64}}, typeof(MatLang.ndgridM), Base.StepRange{Int64, Int64}})
    isdefined(MatLang, Symbol("#kw##onesM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##onesM")), NamedTuple{(:like,), Tuple{Array{Float64, 2}}}, typeof(MatLang.onesM), Array{Int64, 1}})
    isdefined(MatLang, Symbol("#kw##onesM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##onesM")), NamedTuple{(:like,), Tuple{Array{Int32, 2}}}, typeof(MatLang.onesM), Symbol, Int64})
    isdefined(MatLang, Symbol("#kw##onesM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##onesM")), NamedTuple{(:like,), Tuple{Array{Int64, 2}}}, typeof(MatLang.onesM), Int64})
    isdefined(MatLang, Symbol("#kw##randM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##randM")), NamedTuple{(:like,), Tuple{Array{Float64, 2}}}, typeof(MatLang.randM), Array{Int64, 1}})
    isdefined(MatLang, Symbol("#kw##randM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##randM")), NamedTuple{(:like,), Tuple{Array{Int32, 2}}}, typeof(MatLang.randM), Symbol, Int64})
    isdefined(MatLang, Symbol("#kw##randM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##randM")), NamedTuple{(:like,), Tuple{Array{Int64, 2}}}, typeof(MatLang.randM), Int64})
    isdefined(MatLang, Symbol("#kw##sortM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##sortM")), NamedTuple{(:I,), Tuple{Bool}}, typeof(MatLang.sortM), Array{Int64, 1}})
    isdefined(MatLang, Symbol("#kw##sortM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##sortM")), NamedTuple{(:I,), Tuple{Bool}}, typeof(MatLang.sortM), Array{Int64, 2}, Int64, Symbol})
    isdefined(MatLang, Symbol("#kw##zerosM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##zerosM")), NamedTuple{(:like,), Tuple{Array{Float64, 2}}}, typeof(MatLang.zerosM), Array{Int64, 1}})
    isdefined(MatLang, Symbol("#kw##zerosM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##zerosM")), NamedTuple{(:like,), Tuple{Array{Int32, 2}}}, typeof(MatLang.zerosM), Symbol, Int64})
    isdefined(MatLang, Symbol("#kw##zerosM")) && precompile(Tuple{getfield(MatLang, Symbol("#kw##zerosM")), NamedTuple{(:like,), Tuple{Array{Int64, 2}}}, typeof(MatLang.zerosM), Int64})
    precompile(Tuple{typeof(MatLang.__init__)})
    precompile(Tuple{typeof(MatLang.blkdiagM), Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.blkdiagM), Symbol, Array{Float64, 2}, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.catM), Int64, Array{Float64, 2}, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.catM), Int64, Array{Float64, 3}, Array{Float64, 3}})
    precompile(Tuple{typeof(MatLang.clcM)})
    precompile(Tuple{typeof(MatLang.diagM), Array{Int64, 1}, Int64})
    precompile(Tuple{typeof(MatLang.diagM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.diagM), Array{Int64, 2}, Int64})
    precompile(Tuple{typeof(MatLang.diagM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.doubleM), Float64})
    precompile(Tuple{typeof(MatLang.eyeM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.eyeM), Int64, Int64})
    precompile(Tuple{typeof(MatLang.eyeM), Int64})
    precompile(Tuple{typeof(MatLang.eyeM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.eyeM), Type{Int}, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.eyeM), Type{Int}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.eyeM), Type{Int}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.falseM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.falseM), Int64})
    precompile(Tuple{typeof(MatLang.falseM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.falseM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.flipM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.flipM), Array{Int64, 2}, Int64})
    precompile(Tuple{typeof(MatLang.flipM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.flipM), String})
    precompile(Tuple{typeof(MatLang.int16M), Int64})
    precompile(Tuple{typeof(MatLang.int32M), Int64})
    precompile(Tuple{typeof(MatLang.int64M), Int64})
    precompile(Tuple{typeof(MatLang.int8M), Int64})
    precompile(Tuple{typeof(MatLang.iscolumnM), Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Array{String, 2}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Array{Float64, 1}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Array{String, 1}})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.iscolumnM), Symbol, String})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, Array{Float64, 1}})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, Array{String, 1}})
    precompile(Tuple{typeof(MatLang.ismatrixM), Symbol, String})
    precompile(Tuple{typeof(MatLang.isrowM), Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.isrowM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.isrowM), Array{String, 2}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Array{Float64, 1}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Array{String, 1}})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.isrowM), Symbol, String})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Array{Float64, 3}})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Array{String, 1}})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Array{String, 2}})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Bool})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Char})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.isscalarM), Symbol, String})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Array{Float64, 1}})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Array{Float64, 2}})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Array{String, 1}})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.isvectorM), Symbol, String})
    precompile(Tuple{typeof(MatLang.lengthM), Array{Any, 2}})
    precompile(Tuple{typeof(MatLang.lengthM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.linspaceM), Int64, Int64, Int64})
    precompile(Tuple{typeof(MatLang.linspaceM), Int64, Int64})
    precompile(Tuple{typeof(MatLang.linspaceM), Symbol, Int64, Int64, Int64})
    precompile(Tuple{typeof(MatLang.linspaceM), Symbol, Int64, Int64})
    precompile(Tuple{typeof(MatLang.logspaceM), Int64, Int64, Int64})
    precompile(Tuple{typeof(MatLang.logspaceM), Int64, Int64, Symbol})
    precompile(Tuple{typeof(MatLang.logspaceM), Int64, Int64})
    precompile(Tuple{typeof(MatLang.meshgridM), Base.StepRange{Int64, Int64}, Base.StepRange{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.meshgridM), Base.StepRange{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.meshgridM), Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}})
    precompile(Tuple{typeof(MatLang.ndgridM), Base.StepRange{Int64, Int64}, Base.StepRange{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.ndgridM), Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}})
    precompile(Tuple{typeof(MatLang.ndgrid_fill), Array{Int64, 3}, Base.UnitRange{Int64}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.onesM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.onesM), Int64})
    precompile(Tuple{typeof(MatLang.onesM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.onesM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.onesM), Type{Int32}, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.onesM), Type{Int}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.onesM), Type{Int}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.randM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.randM), Int64})
    precompile(Tuple{typeof(MatLang.randM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.randM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.randM), Type{Int32}, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.randM), Type{Int}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.randM), Type{Int}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.repelemM), Array{Int64, 2}, Int64})
    precompile(Tuple{typeof(MatLang.repmatM), Array{Int64, 1}, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.repmatM), Array{Int64, 1}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.repmatM), Array{Int64, 1}, Int64})
    precompile(Tuple{typeof(MatLang.repmatM), Array{Int64, 1}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.repmatM), Int64, Int64, Int64})
    precompile(Tuple{typeof(MatLang.singleM), Float64})
    precompile(Tuple{typeof(MatLang.sizeM), Array{Float64, 3}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.sizeM), Array{Int64, 2}, Int64})
    precompile(Tuple{typeof(MatLang.sizeM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.sizeM), Symbol, Array{Float64, 3}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.sizeM), Symbol, Array{Int64, 2}, Int64})
    precompile(Tuple{typeof(MatLang.sizeM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.sortM), Array{Int64, 1}, Symbol})
    precompile(Tuple{typeof(MatLang.sortM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.sortM), Array{Int64, 2}, Int64, Symbol})
    precompile(Tuple{typeof(MatLang.sortM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.sortM), Array{Integer, 3}, Int64})
    precompile(Tuple{typeof(MatLang.squeezeM), Array{Float64, 3}})
    precompile(Tuple{typeof(MatLang.transposeM), Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.transposeM), Symbol, Array{Int64, 2}})
    precompile(Tuple{typeof(MatLang.trueM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.trueM), Int64})
    precompile(Tuple{typeof(MatLang.trueM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.trueM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.uint16M), Int64})
    precompile(Tuple{typeof(MatLang.uint32M), Int64})
    precompile(Tuple{typeof(MatLang.uint64M), Int64})
    precompile(Tuple{typeof(MatLang.uint8M), Int64})
    precompile(Tuple{typeof(MatLang.zerosM), Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.zerosM), Int64})
    precompile(Tuple{typeof(MatLang.zerosM), Symbol, Int64})
    precompile(Tuple{typeof(MatLang.zerosM), Tuple{Int64, Int64}})
    precompile(Tuple{typeof(MatLang.zerosM), Type{Int32}, Array{Int64, 1}})
    precompile(Tuple{typeof(MatLang.zerosM), Type{Int}, Int64, Int64, Int64})
    precompile(Tuple{typeof(MatLang.zerosM), Type{Int}, Int64, Int64})
    precompile(Tuple{typeof(MatLang.zerosM), Type{Int}, Tuple{Int64, Int64}})
end

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
using Plots

include("rw_1d_stat.jl")

num = length(ARGS)>0 ? parse(ARGS[1]) : 1000
nsteps = length(ARGS)>1 ? parse(ARGS[2]) : 1000
prob = length(ARGS)>2 ? parse(ARGS[3]) : 0.5

means, vars = randomwalk(num, nsteps, prob)

pl = plot(means, xaxis=("t"), yaxis=("E[x]"), color=[:black], legend=false)
savefig(pl, "mean.pdf")
pl = plot(vars, xaxis=("t"), yaxis=("V[x]"), color=[:black], legend=false)
savefig(pl, "var.pdf")

ext: jl

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
"""
  Merge{T}(s::T, offset::Int, nunits::Int, extra::Union{T, Nothing})

Similar to `SubString` but different.
"""
struct Merge{T}
  string::T
  offset::Int
  nunits::Int
  extra::Union{T, Nothing}
end

Merge(a::Merge, e) = Merge(a.string, a.offset, a.nunits, e)
Merge(s::SubString) = Merge(s, nothing)
Merge(s::SubString, e) = Merge(s.string, s.offset, s.ncodeunits, e)
Merge(s::String, e) = Merge(SubString(s), e)

SubString(a::Merge{String}) = SubString(a.string, a.offset+1, prevind(a.string, a.offset+a.nunits+1))

function String(a::Merge{String})
  s = SubString(a)
  if a.extra === nothing
    return String(s)
  else
    return string(s, a.extra)
  end
end

change_extra(a::Merge, s, o, e) = a.extra == o ? Merge(a, e) : Merge(a, s)

function Base.:(*)(a::Merge{T}, b::Merge{T}) where T
  if a.string === b.string
    if a.offset < b.offset
      offset = a.offset
      @assert offset + a.nunits == b.offset "merge un-adjacent Merge"
    elseif a.offset > b.offset
      offset = b.offset
      @assert offset + b.nunits == a.offset "merge un-adjacent Merge"
    else
      error("merge two Merge at same offset: partial string?")
    end
    nunits = a.nunits + b.nunits
    return Merge(a.string, offset, nunits, b.extra)
  else
    error("merge different Merge")
  end
end

function Base.show(io::IO, a::Merge)
  show(io, SubString(a))
  if a.extra !== nothing
    print(io, " + ")
    show(io, a.extra)
  end
  return io
end


text:
using LinearAlgebra, Test
#JSO packages
using NLPModels, NLPModelsIpopt, NLPModelsKnitro, SolverTools
#This package
using Stopping, StoppingInterface

@test true

@testset "Rosenbrock with ∑x = 1" begin
    nlp = ADNLPModel(x->(x[1] - 1.0)^2 + 100 * (x[2] - x[1]^2)^2, [-1.2; 1.0],
                     x->[sum(x)-1], [0.0], [0.0])
    sol = [-1.612771347383541; 2.612771347383541]
    stp = NLPStopping(nlp)
    @test status_stopping_to_stats(stp) == :unknown
    fill_in!(stp, sol)
    @test stop!(stp)
    status_stopping_to_stats(stp) == :first_order
    reinit!(stp, rstate = true, x = zeros(2))
    stp.meta.nb_of_stop = stp.meta.max_iter + 1
    fill_in!(stp, stp.current_state.x)
    @test stop!(stp)
    @test status_stopping_to_stats(stp) == :max_iter

    reinit!(stp)
    stp = knitro(stp)
    @show status(stp), stp.current_state.x
    reinit!(stp)
    stp = ipopt(stp)
    @show status(stp), stp.current_state.x
end

@testset "GSE status to meta" begin
  nlp = ADNLPModel(x->(x[1] - 1.0)^2 + 100 * (x[2] - x[1]^2)^2, [-1.2; 1.0],
                     x->[sum(x)-1], [0.0], [0.0])
  stp = NLPStopping(nlp)
  StoppingInterface.stats_status_to_meta!(stp, :first_order)
  @test stp.meta.optimal
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :acceptable)
  @test stp.meta.suboptimal
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :infeasible)
  @test stp.meta.infeasible
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_step)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_eval)
  @test stp.meta.resources
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_iter)
  @test stp.meta.iteration_limit
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :max_time)
  @test stp.meta.tired
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :neg_pred)
  @test stp.meta.fail_sub_pb
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :not_desc)
  @test stp.meta.fail_sub_pb
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :unbounded)
  @test stp.meta.unbounded
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :user)
  @test stp.meta.stopbyuser
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :stalled)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_residual)
  @test stp.meta.stalled
  reinit!(stp)
  StoppingInterface.stats_status_to_meta!(stp, :small_step)
  @test stp.meta.stalled
end


text:
#################################
### Meteorological variables ####
#################################

#' Air Density
#'
#' Air density of moist air from air temperature and pressure.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - constants Kelvin - conversion degC to Kelvin
#'                  Rd - gas constant of dry air (J kg-1 K-1)
#'                  kPa2Pa - conversion kilopascal (kPa) to pascal (Pa)
#'
#' # Details
 Air density (``\\rho``) is calculated as:
#'
#'   ``\\rho = pressure / (Rd * Tair)``
#'
#' # Value
 - ``\\rho``: air density (kg m-3)
#'
#' ```@example; output = false
#' ```
#' # air density at 25degC and standard pressure (101.325kPa)
#' air_density(25,101.325)
#'
#' #References
#' Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
#'
"""
"""
function air_density(Tair,pressure,constants=BigleafConstants())

  Tair     = Tair + constants.Kelvin
  pressure = pressure * constants.kPa2Pa

  rho = pressure / (constants.Rd * Tair)

  return(rho)
end

#' Atmospheric Pressure from Hypsometric Equation
#'
#' An estimate of mean pressure at a given elevation as predicted by the
#'              hypsometric equation.
#'
#' - elev      Elevation a_s_l. (m)
#' - Tair      Air temperature (deg C)
#' - VPD       Vapor pressure deficit (kPa); optional
#' - constants Kelvin- conversion degC to Kelvin
#'                  pressure0 - reference atmospheric pressure at sea level (Pa)
#'                  Rd - gas constant of dry air (J kg-1 K-1)
#'                  g -  gravitational acceleration (m s-2)
#'                  Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 Atmospheric pressure is approximated by the hypsometric equation:
#'
#'          ``pressure = pressure_0 / (exp(g * elevation / (Rd Temp)))``
#'
#' # Note
#' The hypsometric equation gives an estimate of the standard pressure
#'       at a given altitude.
#'       If VPD is provided, humidity correction is applied and the
#'       virtual temperature instead of air temperature is used. VPD is
#'       internally converted to specific humidity.
#'
#' # Value
 - pressure -: Atmospheric pressure (kPa)
#'
#' #References
#' Stull B., 1988: An Introduction to Boundary Layer Meteorology.
#'             Kluwer Academic Publishers, Dordrecht, Netherlands.
#'
#' ```@example; output = false
#' ```
#' # mean pressure at 500m altitude at 25 deg C and VPD of 1 kPa
#' pressure_from_elevation(500,Tair=25,VPD=1)
#'
#' @export
function pressure_from_elevation(elev,Tair,VPD=nothing,constants=BigleafConstants())

  Tair     = Tair + constants.Kelvin

  if(isnothing(VPD))

    pressure = constants.pressure0] / exp(constants[:g] * elev / (constants[:Rd*Tair))

else

    pressure1   = constants.pressure0] / exp(constants[:g] * elev / (constants[:Rd*Tair))
    Tv          = virtual_temp(Tair - constants.Kelvin],pressure1 * constants[:Pa2kPa,
                                VPD,Esat_formula="Sonntag_1990",constants) + constants.Kelvin

    pressure    = constants.pressure0] / exp(constants[:g] * elev / (constants[:Rd*Tv))
end

  pressure = pressure * constants.Pa2kPa
  return(pressure)
end

#' Saturation Vapor Pressure (Esat) and Slope of the Esat Curve
#'
#' Calculates saturation vapor pressure (Esat) over water and the
#'              corresponding slope of the saturation vapor pressure curve.
#'
#' - Tair      Air temperature (deg C)
#' - Esat_formula   Esat_formula to be used. Either `"Sonntag_1990"` (Default),
#'                  `"Alduchov_1996"`, or `"Allen_1998"`.
#' - constants Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 Esat (kPa) is calculated using the Magnus equation:
#'
#'    ``Esat = a * exp((b * Tair) / (c + Tair)) / 1000``
#'
#'  where the coefficients a, b, c take different values depending on the formula used.
#'  The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This version
#'  of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternatively,
#'  parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can be
#'  used (see references).
#'  The slope of the Esat curve (``\\Delta``) is calculated as the first derivative of the function:
#'
#'    ``\\Delta = dEsat / dTair``
#'
#'  which is solved using `\link[stats]{D`}.
#'
#' # Value
 A dataframe with the following columns:
#'  - Esat: Saturation vapor pressure (kPa)
#'  - Delta: Slope of the saturation vapor pressure curve (kPa K-1)
#'
#' #References
#' Sonntag D. 1990: Important new values of the physical constants of 1986, vapor
#'             pressure formulations based on the ITS-90 and psychrometric formulae.
#'             Zeitschrift fuer Meteorologie 70, 340-344.
#'
#'             World Meteorological Organization 2008: Guide to Meteorological Instruments
#'             and Methods of Observation (WMO-No.8). World Meteorological Organization,
#'             Geneva. 7th Edition.
#'
#'             Alduchov, O. A. & Eskridge, R. E., 1996: Improved Magnus form approximation of
#'             saturation vapor pressure. Journal of Applied Meteorology, 35, 601-609
#'
#'             Allen, R_G., Pereira, L_S., Raes, D., Smith, M., 1998: Crop evapotranspiration -
#'             Guidelines for computing crop water requirements - FAO irrigation and drainage
#'             paper 56, FAO, Rome.
#'
#' ```@example; output = false
#' ```
#' Esat_slope(seq(0,45,5))[,"Esat"]  # Esat in kPa
#' Esat_slope(seq(0,45,5))[,"Delta"] # the corresponding slope of the Esat curve (Delta) in kPa K-1
#'
#' @importFrom stats D
"""
"""
function Esat_slope(Tair,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                       constants=BigleafConstants())

  Esat_formula = match_arg(Esat_formula)

  if (Esat_formula == "Sonntag_1990")
    a = 611.2
    b = 17.62
    c = 243.12
elseif (Esat_formula == "Alduchov_1996")
    a = 610.94
    b = 17.625
    c = 243.04
elseif (Esat_formula == "Allen_1998")
    a = 610.8
    b = 17.27
    c = 237.3
end

  # saturation vapor pressure
  Esat = a * exp((b * Tair) / (c + Tair))
  Esat = Esat * constants.Pa2kPa

  # slope of the saturation vapor pressure curve
  Delta = eval(D(expression(a * exp((b * Tair) / (c + Tair))),name="Tair"))
  Delta = Delta * constants.Pa2kPa

  return(DataFrame(Esat,Delta))
end

#' Psychrometric Constant
#'
#' Calculates the psychrometric 'constant'.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - constants cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  eps - ratio of the molecular weight of water vapor to dry air (-)
#'
#' # Details
 The psychrometric constant (``\\gamma``) is given as:
#'
#'    ``\\gamma = cp * pressure / (eps * \\lambda)``
#'
#'  where ``\\lambda`` is the latent heat of vaporization (J kg-1),
#'  as calculated from [`latent_heat_vaporization`](@ref).
#'
#' # Value
 - ``\\gamma`` -: the psychrometric constant (kPa K-1)
#'
#' #References
#' Monteith J_L., Unsworth M_H., 2008: Principles of Environmental Physics.
#'             3rd Edition. Academic Press, London.
#'
#' ```@example; output = false
#' ```
#' psychrometric_constant(seq(5,45,5),100)
#'
"""
"""
function psychrometric_constant(Tair,pressure,constants=BigleafConstants())

  lambda = latent_heat_vaporization(Tair)
  gamma  = (constants.cp] * pressure) / (constants[:eps * lambda)

  return(gamma)
end

#' Latent Heat of Vaporization
#'
#' Latent heat of vaporization as a function of air temperature.
#'
#' - Tair  Air temperature (deg C)
#'
#' # Details
 The following Esat_formula is used:
#'
#'   ``\\lambda = (2.501 - 0.00237*Tair)10^6``
#'
#' # Value
 - ``\\lambda`` -: Latent heat of vaporization (J kg-1)
#'
#' #References
#' Stull, B., 1988: An Introduction to Boundary Layer Meteorology (p.641)
#'             Kluwer Academic Publishers, Dordrecht, Netherlands
#'
#'             Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
#'
#' ```@example; output = false
#' ```
#' latent_heat_vaporization(seq(5,45,5))
#'
"""
"""
function latent_heat_vaporization(Tair)

  k1 = 2.501
  k2 = 0.00237
  lambda = ( k1 - k2 * Tair ) * 1e+06

  return(lambda)
end

#' Solver Function for Wet-Bulb Temperature
#'
#' Solver function used in wetbulb_temp()
#'
#' - ea           Air vapor pressure (kPa)
#' - Tair         Air temperature (degC)
#' - gamma        Psychrometric constant (kPa K-1)
#' - accuracy     Accuracy of the result (degC)
#' - Esat_formula Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                     One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                     See [`Esat_slope`](@ref).
#' - constants    Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Note
#' Arguments `accuracy` and `Esat_formula` are passed to this function by wetbulb_temp().
#'
#' @importFrom stats optimize
#'
#' @keywords internal
function wetbulb_solver(ea,Tair,gamma,accuracy,Esat_formula,constants=BigleafConstants())
  wetbulb_optim = optimize(function(Tw){abs(ea - c((Esat_slope(Tw,Esat_formula,constants)[,"Esat"] - 0.93*gamma*(Tair - Tw))))},
                            interval=c(-100,100),tol=accuracy)
  return(wetbulb_optim)
end

#' Wet-Bulb Temperature
#'
#' calculates the wet bulb temperature, i.e. the temperature
#'              that the air would have if it was saturated.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - VPD       Vapor pressure deficit (kPa)
#' - accuracy  Accuracy of the result (deg C)
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  eps - ratio of the molecular weight of water vapor to dry air (-)
#'                  Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 Wet-bulb temperature (Tw) is calculated from the following expression:
#'
#'            ``e = Esat(Tw) - gamma* (Tair - Tw)``
#'
#'          The equation is solved for Tw using `\link[stats]{optimize`}.
#'          Actual vapor pressure e (kPa) is calculated from VPD using the function [`VPD_to_e`](@ref).
#'          The psychrometric constant gamma (kPa K-1) is calculated from [`psychrometric_constant`](@ref).
#'
#' # Value
 - Tw -: wet-bulb temperature (degC)
#'
#' #References
#' Monteith J_L., Unsworth M_H., 2008: Principles of Environmental Physics.
#'             3rd edition. Academic Press, London.
#'
#' ```@example; output = false
#' ```
#' wetbulb_temp(Tair=c(20,25),pressure=100,VPD=c(1,1.6))
#'
#' @importFrom stats optimize
"""
"""
function wetbulb_temp(Tair,pressure,VPD,accuracy=1e-03,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                         constants=BigleafConstants())

  if (!is_numeric(accuracy))
    stop("'accuracy' must be numeric!")
end

  if (accuracy > 1)
    print("'accuracy' is set to 1 degC")
    accuracy = 1
end

  # determine number of digits to print
  ndigits = as_numeric(strsplit(format(accuracy,scientific = true),"-")[[1]][2])
  ndigits = ifelse(ismissing(ndigits),0,ndigits)

  gamma  = psychrometric_constant(Tair,pressure)
  ea     = VPD_to_e(VPD,Tair,Esat_formula)

  Tw = sapply(seq_along(ea),function(i) round(wetbulb_solver(ea[i],Tair[i],gamma[i],
                                                              accuracy=accuracy,Esat_formula,constants)$minimum,ndigits))

  return(Tw)

end

#' Solver function for dew point temperature
#'
#' Solver function used in dew_point()
#'
#' - ea           Air vapor pressure (kPa)
#' - accuracy     Accuracy of the result (degC)
#' - Esat_formula Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                     One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                     See [`Esat_slope`](@ref).
#' - constants    Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Note
#' Arguments `accuracy` and `Esat_formula` are passed to this function by dew_point().
#'
#' @importFrom stats optimize
#'
#' @keywords internal
function dew_point_solver(ea,accuracy,Esat_formula,constants=BigleafConstants())

  Td_optim = optimize(function(Td){abs(ea - Esat_slope(Td,Esat_formula,constants)[,"Esat"])},
                       interval=c(-100,100),tol=accuracy)
  return(Td_optim)
end

#' Dew Point
#'
#' calculates the dew point, the temperature to which air must be
#'              cooled to become saturated (i.e. e = Esat(Td))
#'
#' - Tair     Air temperature (degC)
#' - VPD      Vapor pressure deficit (kPa)
#' - accuracy Accuracy of the result (deg C)
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 Dew point temperature (Td) is defined by:
#'
#'           ``e = Esat(Td)``
#'
#'          where e is vapor pressure of the air and Esat is the vapor pressure deficit.
#'          This equation is solved for Td using `\link[stats]{optimize`}.
#'
#' # Value
 - Td -: dew point temperature (degC)
#'
#' #References
#' Monteith J_L., Unsworth M_H., 2008: Principles of Environmental Physics.
#'             3rd edition. Academic Press, London.
#'
#' ```@example; output = false
#' ```
#' dew_point(c(25,30),1.5)
#'
#' @importFrom stats optimize
#' @export
function dew_point(Tair,VPD,accuracy=1e-03,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                      constants=BigleafConstants())

  if (!is_numeric(accuracy))
    stop("'accuracy' must be numeric!")
end

  if (accuracy > 1)
    print("'accuracy' is set to 1 degC")
    accuracy = 1
end

  # determine number of digits to print
  ndigits = as_numeric(strsplit(format(accuracy,scientific = true),"-")[[1]][2])
  ndigits = ifelse(ismissing(ndigits),0,ndigits)

  ea = VPD_to_e(VPD,Tair,Esat_formula)
  Td = sapply(seq_along(ea),function(i) round(dew_point_solver(ea[i],accuracy=accuracy,
                                                                Esat_formula,constants)$minimum,ndigits))

  return(Td)
end

#' Virtual Temperature
#'
#' Virtual temperature, defined as the temperature at which dry air would have the same
#'              density as moist air at its actual temperature.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - VPD       Vapor pressure deficit (kPa)
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  eps - ratio of the molecular weight of water vapor to dry air (-)
#'
#' # Details
 the virtual temperature is given by:
#'
#'    ``Tv = Tair / (1 - (1 - eps) e/pressure)``
#'
#'  where Tair is in Kelvin (converted internally). Likewise, VPD is converted
#'  to actual vapor pressure (e in kPa) with [`VPD_to_e`](@ref) internally.
#'
#' # Value
 - Tv -: virtual temperature (deg C)
#'
#' #References
#' Monteith J_L., Unsworth M_H., 2008: Principles of Environmental Physics.
#'             3rd edition. Academic Press, London.
#'
#' ```@example; output = false
#' ```
#' virtual_temp(25,100,1.5)
#'
"""
"""
function virtual_temp(Tair,pressure,VPD,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                         constants=BigleafConstants())

  e    = VPD_to_e(VPD,Tair,Esat_formula)
  Tair = Tair + constants.Kelvin

  Tv = Tair / (1 - (1 - constants.eps) * e/pressure)
  Tv = Tv - constants.Kelvin

  return(Tv)
end

#' Kinematic Viscosity of Air
#'
#' calculates the kinematic viscosity of air.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  pressure0 - reference atmospheric pressure at sea level (Pa)
#'                  Tair0 - reference air temperature (K)
#'                  kPa2Pa - conversion kilopascal (kPa) to pascal (Pa)
#'
#' # Details
 where v is the kinematic viscosity of the air (m2 s-1),
#'          given by (Massman 1999b):
#'
#'            ``v = 1.327 * 10^-5(pressure0/pressure)(Tair/Tair0)^1.81``
#'
#' # Value
 - v -: kinematic viscosity of air (m2 s-1)
#'
#' #References
#' Massman, W_J., 1999b: Molecular diffusivities of Hg vapor in air,
#'             O2 and N2 near STP and the kinematic viscosity and thermal diffusivity
#'             of air near STP. Atmospheric Environment 33, 453-457.
#'
#' ```@example; output = false
#' ```
#' kinematic_viscosity(25,100)
#'
#' @export
function kinematic_viscosity(Tair,pressure,constants=BigleafConstants())

  Tair     = Tair + constants.Kelvin
  pressure = pressure * constants.kPa2Pa

  v  = 1.327e-05*(constants.pressure0]/pressure)*(Tair/constants[:Tair0)^1.81
  return(v)
end


text:
# let
    net = Network()

    # add a success node
    success = add_node(net, DSL_CPT, "Success")
    somenames = IdArray()
    add(somenames, "Success")
    add(somenames, "Failure")
    set_number_of_outcomes(definition(get_node(net, success)), somenames)

    # add a forecast node
    forecast = add_node(net, DSL_CPT, "Forecast")
    flush(somenames)
    add(somenames, "Good")
    add(somenames, "Moderate")
    add(somenames, "Poor")
    set_number_of_outcomes(definition(get_node(net, forecast)), somenames)

    # add an arc from "Success" to "Forecast" to represent conditional independence
    add_arc(net, success, forecast)

    # add the distributions
    theprobs = DoubleArray()
    set_size(theprobs, 2)
    set_at(theprobs, 0, 0.2)
    set_at(theprobs, 1, 0.8)
    set_definition(definition(get_node(net, success)), theprobs)

    thecoordinates = SysCoordinates(definition(get_node(net, forecast)))
    set_unchecked_value(thecoordinates, 0.4); next(thecoordinates)
    set_unchecked_value(thecoordinates, 0.4); next(thecoordinates)
    set_unchecked_value(thecoordinates, 0.2); next(thecoordinates)
    set_unchecked_value(thecoordinates, 0.1); next(thecoordinates)
    set_unchecked_value(thecoordinates, 0.3); next(thecoordinates)
    set_unchecked_value(thecoordinates, 0.6)

    test_vec_approx_equal(get_cpt_probability_vec(net, success, Dict{Cint, Cint}()), [0.2, 0.8], 1e-7)
    test_vec_approx_equal(get_cpt_probability_vec(net, forecast, Dict(success=>Int32(0))), [0.4, 0.4, 0.2], 1e-7)
    test_vec_approx_equal(get_cpt_probability_vec(net, forecast, Dict(success=>Int32(1))), [0.1, 0.3, 0.6], 1e-7)

    set_default_BN_algorithm(net, DSL_ALG_BN_LAURITZEN)
    clear_all_evidence(net)
    rand(net)

    evidence = Dict(forecast=>Int32(0)) # forecast is good
    @test( does_assignment_match_evidence(Dict(success=>Int32(0), forecast=>Int32(0)), evidence))
    @test( does_assignment_match_evidence(Dict(success=>Int32(1), forecast=>Int32(0)), evidence))
    @test(!does_assignment_match_evidence(Dict(success=>Int32(0), forecast=>Int32(1)), evidence))

    @test(abs(monte_carlo_probability_estimate(net, evidence, 1000) - 0.16) < 0.1)
    @test(abs(monte_carlo_probability_estimate(net, Dict(success=>Int32(1)), evidence, 100000) - 0.5) < 0.1)

    # println(rand!(net, [forecast=>Int32(2)]))
    # println(rand!(net, [success =>Int32(1)]))

    @test(isapprox(marginal_probability(net, Dict(success=>Int32(0)), Dict(forecast=>Int32(0))), 0.5))
    @test(isapprox(marginal_probability(net, Dict(success=>Int32(1)), Dict(forecast=>Int32(0))), 0.5))
    @test(isapprox(marginal_probability(net, Dict(success=>Int32(0)), Dict(forecast=>Int32(1))), 0.25))
    @test(isapprox(marginal_probability(net, Dict(success=>Int32(1)), Dict(forecast=>Int32(1))), 0.75))
    @test(isapprox(marginal_probability(net, Dict(success=>Int32(0)), Dict(forecast=>Int32(2))), 0.0769, atol=0.001))
    @test(isapprox(marginal_probability(net, Dict(success=>Int32(1)), Dict(forecast=>Int32(2))), 0.9231, atol=0.001))

    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(0)), Dict(success=>Int32(0))), 0.4))
    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(1)), Dict(success=>Int32(0))), 0.4))
    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(2)), Dict(success=>Int32(0))), 0.2))
    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(0)), Dict(success=>Int32(1))), 0.1))
    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(1)), Dict(success=>Int32(1))), 0.3))
    @test(isapprox(marginal_probability(net, Dict(forecast=>Int32(2)), Dict(success=>Int32(1))), 0.6))

    @test(isapprox(probability(net, Dict(success=>Int32(0), forecast=>Int32(0))), 0.2*0.4))
    @test(isapprox(probability(net, Dict(success=>Int32(0), forecast=>Int32(1))), 0.2*0.4))
    @test(isapprox(probability(net, Dict(success=>Int32(0), forecast=>Int32(2))), 0.2*0.2))
    @test(isapprox(probability(net, Dict(success=>Int32(1), forecast=>Int32(0))), 0.8*0.1))
    @test(isapprox(probability(net, Dict(success=>Int32(1), forecast=>Int32(1))), 0.8*0.3))
    @test(isapprox(probability(net, Dict(success=>Int32(1), forecast=>Int32(2))), 0.8*0.6))

    @test(isapprox(probability(net, Dict(success=>Int32(0))), 0.2))
    @test(isapprox(probability(net, Dict(success=>Int32(1))), 0.8))

    @test(isapprox(probability(net, Dict(forecast=>Int32(0))), 0.16))
    @test(isapprox(probability(net, Dict(forecast=>Int32(1))), 0.32))
    @test(isapprox(probability(net, Dict(forecast=>Int32(2))), 0.52))

    @test(isapprox(logprob(net, Dict(success=>Int32(0), forecast=>Int32(0))), log(0.2*0.4)))
    @test(isapprox(logprob(net, Dict(success=>Int32(0), forecast=>Int32(1))), log(0.2*0.4)))
    @test(isapprox(logprob(net, Dict(success=>Int32(0), forecast=>Int32(2))), log(0.2*0.2)))
    @test(isapprox(logprob(net, Dict(success=>Int32(1), forecast=>Int32(0))), log(0.8*0.1)))
    @test(isapprox(logprob(net, Dict(success=>Int32(1), forecast=>Int32(1))), log(0.8*0.3)))
    @test(isapprox(logprob(net, Dict(success=>Int32(1), forecast=>Int32(2))), log(0.8*0.6)))

    @test(isapprox(logprob(net, Dict(success=>Int32(0))), log(0.2)))
    @test(isapprox(logprob(net, Dict(success=>Int32(1))), log(0.8)))

    @test(isapprox(logprob(net, Dict(forecast=>Int32(0))), log(0.16)))
    @test(isapprox(logprob(net, Dict(forecast=>Int32(1))), log(0.32)))
    @test(isapprox(logprob(net, Dict(forecast=>Int32(2))), log(0.52)))
# end


text:
const difffn = Dict(
    :Central => :CenteredDifference,
    :Center => :CenteredDifference,
    :Centered => :CenteredDifference,
    :central => :CenteredDifference,
    :center => :CenteredDifference,
    :centered => :CenteredDifference,
    :CenteredDifference => :CenteredDifference,
    :Upwind => :UpwindDifference,
    :upwind => :UpwindDifference,
    :UpwindDifference => :UpwindDifference,
)

const bcfn = Dict(
    :Period => :PeriodicBC,
    :Periodic => :PeriodicBC,
    :period => :PeriodicBC,
    :Period => :PeriodicBC,
    :PeriodBC => :PeriodicBC,
    :Dirichlet => :DirichletBC,
    :dirichlet => :DirichletBC,
    :DirichletBC => :DirichletBC,
    :none => nothing,
)

*(Δ::AbstractDerivativeOperator, B::Nothing) = Δ

"""
$(SIGNATURES)

Finite difference operation

Note that for upwind difference, the coeff = 1/-1 can't take the default value as nothing.
"""
function finite_difference(
    y::AA,
    x::AA,
    coeff = 1,
    uL = first(y),
    uR = last(y);
    method = :central,
    dimension = 1,
    derivative = 1,
    truncation = 2,
    bc = :dirichlet,
)

    oprtfn = eval(difffn[method])

    i0 = firstindex(x)
    i1 = lastindex(x)

    if bc in (nothing, :none)
        dx = x[i0+1:i1] .- x[i0:i1-1]
        Δ = oprtfn{dimension}(derivative, truncation, dx, length(y) - 2, coeff)
        B = nothing
    else
        dx = [x[i0+1] - x[i0]; x[i0+1:i1] .- x[i0:i1-1]; x[i1] - x[i1-1]]
        Δ = oprtfn{dimension}(derivative, truncation, dx, length(y), coeff)
    end

    bcnm = bcfn[bc]
    bcrtfn = eval(bcnm)

    if bcnm == :PeriodicBC
        B = bcrtfn(eltype(y))
    elseif bcnm == :DirichletBC
        B = bcrtfn(uL, uR)
    end

    return Δ * B * y

end

function finite_difference(y, dx::Real, args...; kwargs...)
    x = linspace(0, length(y) - 1, length(y)) .* dx
    return finite_difference(y, x, args...; kwargs...)
end

"""
$(SIGNATURES)

Central difference
"""
function central_diff(y::T, x::T) where {T<:AV{<:Number}}
    dy = zeros(eltype(y), axes(y))

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    dy[i0] = (y[i0+1] - y[i0]) / (x[i0+1] - x[i0] + 1e-7)
    dy[i1] = (y[i1] - y[i1-1]) / (x[i1] - x[i1-1] + 1e-7)
    for i = i0+1:i1-1
        dy[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1] + 1e-7)
    end

    return dy
end

"""
$(SIGNATURES)
"""
function central_diff(y::AV{T}, dx) where {T}
    x = linspace(0, length(y) - 1, length(y)) .* dx
    dy = central_diff(y, x)

    return dy
end

"""
$(SIGNATURES)

Central difference for second-order derivatives
"""
function central_diff2(y::T, x::T) where {T<:AV{<:Number}}
    dy = zeros(eltype(y), axes(y))

    i0 = eachindex(y) |> first
    i1 = eachindex(y) |> last

    for i = i0+1:i1-1
        dy[i] = (y[i+1] - 2.0 * y[i] + y[i-1]) / (x[i+1] - x[i-1])^2
    end

    return dy
end

"""
$(SIGNATURES)
"""
function central_diff2(y::AV{T}, dx) where {T}
    x = linspace(0, length(y) - 1, length(y)) .* dx
    dy = central_diff2(y, x)

    return dy
end

"""
$(SIGNATURES)

Central difference
"""
function central_diff!(dy::AV{T1}, y::T2, x::T2) where {T1,T2<:AV{<:Number}}
    @assert axes(dy) == axes(y) == axes(x)

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    dy[i0] = (y[i0+1] - y[i0]) / (x[i0+1] - x[i0] + 1e-7)
    dy[i1] = (y[i1] - y[i1-1]) / (x[i1] - x[i1-1] + 1e-7)
    for i = i0+1:i1-1
        dy[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1] + 1e-7)
    end

    return nothing
end

"""
$(SIGNATURES)
"""
function central_diff!(dy::AV{T1}, y::AV{T2}, dx) where {T1,T2}
    x = linspace(0, length(y) - 1, length(y)) .* dx
    central_diff!(dy, y, x)

    return nothing
end

"""
$(SIGNATURES)

Central difference
"""
function central_diff2!(dy::AV{T1}, y::T2, x::T2) where {T1,T2<:AV{<:Number}}
    @assert axes(dy) == axes(y) == axes(x)

    i0 = eachindex(y) |> first
    i1 = eachindex(y) |> last

    dy[i0] = 0.0
    dy[i1] = 0.0
    for i = i0+1:i1-1
        dy[i] = (y[i+1] - 2.0 * y[i] + y[i-1]) / (x[i+1] - x[i-1])^2
    end

    return nothing
end

"""
$(SIGNATURES)
"""
function central_diff2!(dy::AV{T1}, y::AV{T2}, dx) where {T1,T2}
    x = ones(eltype(y), axes(y)) .* dx
    central_diff2!(dy, y, x)

    return nothing
end

"""
$(SIGNATURES)

Upwind difference
"""
function upwind_diff(y::AV{T1}, x::AV{T2}; stream = :right::Symbol) where {T1,T2}
    dy = zeros(eltype(y), axes(y))

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    if stream == :right
        dy[i0] = 0.0
        for i = i0+1:i1
            dy[i] = (y[i] - y[i-1]) / (x[i] - x[i-1] + 1e-7)
        end
    elseif stream == :left
        dy[i1] = 0.0
        for i = i0:i1-1
            dy[i] = (y[i+1] - y[i]) / (x[i+1] - x[i] + 1e-7)
        end
    else
        throw("streaming direction should be :left or :right")
    end

    return dy
end

"""
$(SIGNATURES)
"""
function upwind_diff(y::AV{T}, dx; stream = :right::Symbol) where {T}
    x = linspace(0, length(y) - 1, length(y)) .* dx
    dy = upwind_diff(y, x, stream = stream)

    return dy
end

"""
$(SIGNATURES)

Upwind difference
"""
function upwind_diff!(
    dy::AV{T1},
    y::AV{T2},
    x::AV{T3};
    stream = :right::Symbol,
) where {T1,T2,T3}

    @assert axes(dy) == axes(y) == axes(x)

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    if stream == :right
        dy[i0] = 0.0
        for i = i0+1:i1
            dy[i] = (y[i] - y[i-1]) / (x[i] - x[i-1] + 1e-7)
        end
    elseif stream == :left
        dy[i1] = 0.0
        for i = i0:i1-1
            dy[i] = (y[i+1] - y[i]) / (x[i+1] - x[i] + 1e-7)
        end
    else
        throw("streaming direction should be :left or :right")
    end

    return nothing

end

"""
$(SIGNATURES)
"""
function upwind_diff!(dy::AV{T1}, y::AV{T2}, dx; stream = :right::Symbol) where {T1,T2}
    x = linspace(0, length(y) - 1, length(y)) .* dx
    upwind_diff!(dy, y, x, stream = stream)

    return nothing
end

"""
$(SIGNATURES)

Finite difference for pseudo-unstructured mesh
"""
function unstruct_diff(
    u::AV{T1},
    x::AV{T2},
    nx::Integer;
    mode = :central::Symbol,
) where {T1,T2}

    uu = reshape(u, (nx, :))
    xx = reshape(x, (nx, :))

    dux = similar(xx)
    for i = 1:nx
        if mode == :central
            dux[i, :] .= central_diff(uu[i, :], xx[i, :])
        elseif mode == :upwind
            dux[i, :] .= upwind_diff(uu[i, :], xx[i, :])
        else
            throw("difference mode should be central or upwind")
        end
    end

    return reshape(dux, (1, :))

end

"""
$(SIGNATURES)
"""
function unstruct_diff(
    u::Function,
    x::AV{T},
    nx::Integer,
    dim::Integer;
    mode = :central::Symbol,
) where {T}

    uu = reshape(u.(x), (nx, :))
    xx = reshape(x, (nx, :))
    dux = zeros(eltype(x), axes(xx))

    if dim == 1
        for j in axes(xx, 2)
            if mode == :central
                dux[:, j] .= central_diff(uu[:, j], xx[:, j])
            elseif mode == :upwind
                dux[:, j] .= upwind_diff(uu[:, j], xx[:, j])
            end
        end
    elseif dim == 2
        for i in axes(xx, 1)
            if mode == :central
                dux[i, :] .= central_diff(uu[i, :], xx[i, :])
            elseif mode == :upwind
                dux[i, :] .= upwind_diff(uu[i, :], xx[i, :])
            end
        end
    end

    return dux[:]

end


text:
module ThinPlateSplines

using LinearAlgebra

export ThinPlateSpline
export tps_solve, tps_energy, tps_deform

struct ThinPlateSpline
	λ  # Stiffness.
	x1 # control points
	Y  # Homogeneous control point coordinates
	Φ  # TPS kernel
	d  # Affine component
	c  # Non-affine component
end

# Thin-plate splines.
# based on http://en.wikipedia.org/wiki/Thin_plate_spline

is_zero(r::AbstractFloat) = abs(r)<eps(r)
is_zero(r) = false

tps_basis(r::T) where T<:Any  = is_zero(r) ? zero(T) : r*r*log(r)

my_norm(a) = sqrt(sum(a.^2))

# x: matrix of size KxD
tps_kernel(x) = [tps_basis(my_norm(x[i,:] - x[j,:])) for i=1:size(x,1),j=1:size(x,1)]

# find solution of tps transformation
# compute_affine: compute affine component
# (required for some operations but takes additional time.)
function tps_solve(x,y,λ,compute_affine=true)
	# D: number of dimensions
	# K: number of data points
	K,D = size(x)

	# homogeneous coordinates
	X=cat(dims=2,ones(K,1),x)
	Y=cat(dims=2,ones(K,1),y)

	# compute TPS kernel
	Φ = tps_kernel(x)

	# full QR decomposition
	Q,r = qr(X)
	q1 = Q[:,1:(D+1)]
	q2 = Q[:,(D+2):end]

	# warping coefficients
	c = q2*inv(UniformScaling(λ) + q2'*Φ*q2)*q2'*Y

	# affine component
	d = compute_affine ?  r\(q1'*(Y - Φ*c)) : []
	return ThinPlateSpline(λ,x,Y,Φ,d,c)
end

# Thin-plate spline bending energy at minimum
tps_energy(tps::ThinPlateSpline) = tps.λ*tr(tps.c*tps.Y')

# Deform points according to thin-plate spline.
# tps is the thin-plate spline.
# x2 are coordinates of points to be deformed.
function tps_deform(x2::AbstractArray,tps::ThinPlateSpline) where T<:Any
	x1,d,c=tps.x1,tps.d,tps.c
	d==[] && throw(ArgumentError("Affine component not available; run tps_solve with compute_affine=true."))

	# deform
	y2 = zeros(eltype(x2),size(x2,1),size(x2,2)+1)
	for i = 1 : size(x2,1)
		z = x2[i:i,:]
		defc = zeros(eltype(x2),1,size(x2,2)+1)
		for j = 1 : size(x1,1)
			n = my_norm(z - x1[j:j,:])
			defc += tps_basis(n)*c[j:j,:]
		end
		y2[i:i,:] = cat(dims=2, 1.0, z)*d + defc
	end

	y2[:,2:end]
end

# Deform points according to thin-plate spline,
# solving to find c and d.
# x1 are coordinates of control points.
# y are coordinates of deformed control points;
# x2 are coordinates of points to be deformed.
function tps_deform(x1,x2,y,λ)
	tps = tps_solve(x1,y,λ,compute_affine=true)

	tps_deform_with_param(x1,x2,tps)
end

end


text:
### jacvec functions
### should be upstreamed to ForwardDiff

struct JacVecTag end

# J(f(x))*v
function auto_jacvec!(du, f, x, v,
                 cache1 = ForwardDiff.Dual{JacVecTag}.(x, v),
                 cache2 = ForwardDiff.Dual{JacVecTag}.(x, v))
    cache1 .= ForwardDiff.Dual{JacVecTag}.(x, v)
    f(cache2,cache1)
    du .= ForwardDiff.partials.(cache2, 1)
end
function auto_jacvec(f, x, v)
    ForwardDiff.partials.(f(ForwardDiff.Dual{JacVecTag}.(x, v)), 1)
end

function num_jacvec!(du,f,x,v,cache1 = similar(v),
                 cache2 = similar(v);
                 compute_f0 = true)
    compute_f0 && (f(cache1,x))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    f(cache2,x)
    @. x -= ϵ*v
    @. du = (cache2 - cache1)/ϵ
end
function num_jacvec(f,x,v,f0=nothing)
    f0 === nothing ? _f0 = f(x) : _f0 = f0
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    (f(x.+ϵ.*v) .- f(x))./ϵ
end

#=
### vecjec
### Not useful right now

"""
    vecjac(f, x, v) -> u
``v'J(f(x))``
"""
function vecjac(f, x, v)
    tp = ReverseDiff.InstructionTape()
    tx = ReverseDiff.track(x, tp)
    ty = f(tx)
    ReverseDiff.increment_deriv!(ty, v)
    ReverseDiff.reverse_pass!(tp)
    return ReverseDiff.deriv(tx)'
end
=#

### Operator Implementation

mutable struct JacVecOperator{T,F,T1,T2,uType,P,tType,O} <: DiffEqBase.AbstractDiffEqLinearOperator{T}
    f::F
    cache1::T1
    cache2::T2
    u::uType
    p::P
    t::tType
    autodiff::Bool
    ishermitian::Bool
    opnorm::O

    function JacVecOperator{T}(f,p=nothing,t::Union{Nothing,Number}=nothing;autodiff=true,ishermitian=false,opnorm=true) where T
        p===nothing ? P = Any : P = typeof(p)
        t===nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),Nothing,Nothing,Any,P,tType,typeof(opnorm)}(f,nothing,nothing,nothing,nothing,nothing,autodiff,ishermitian)
    end
    function JacVecOperator{T}(f,u::AbstractArray,p=nothing,t::Union{Nothing,Number}=nothing;autodiff=true,ishermitian=false,opnorm=true) where T
        if autodiff
            cache1 = ForwardDiff.Dual{JacVecTag}.(u, u)
            cache2 = ForwardDiff.Dual{JacVecTag}.(u, u)
        else
            cache1 = similar(u)
            cache2 = similar(u)
        end
        p===nothing ? P = Any : P = typeof(p)
        t===nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),typeof(cache1),typeof(cache2),typeof(u),P,tType,typeof(opnorm)}(f,cache1,cache2,u,p,t,autodiff,ishermitian,opnorm)
    end
    function JacVecOperator(f,u,args...;kwargs...)
        JacVecOperator{eltype(u)}(f,u,args...;kwargs...)
    end
end

LinearAlgebra.opnorm(L::JacVecOperator, p::Real=2) = L.opnorm
LinearAlgebra.ishermitian(L::JacVecOperator) = L.ishermitian

Base.size(L::JacVecOperator) = (length(L.cache1),length(L.cache1))
Base.size(L::JacVecOperator,i::Int) = length(L.cache1)
function update_coefficients!(L::JacVecOperator,u,p,t)
    L.u = u
    L.p = p
    L.t = t
    !L.autodiff && L.cache1 !== nothing && L.f(L.cache1,L.u,L.p,L.t)
end

# Interpret the call as df/du * u

function (L::JacVecOperator)(u,p,t::Number)
    update_coefficients!(L,u,p,t)
    L*u
end

function (L::JacVecOperator)(du,u,p,t::Number)
    update_coefficients!(L,u,p,t)
    mul!(du,L,u)
end

Base.:*(L::JacVecOperator,x::AbstractVector) = L.autodiff ? auto_jacvec(_u->L.f(_u,L.p,L.t),L.u,x) : num_jacvec(_u->L.f(_u,L.p,L.t),L.u,x)

function LinearAlgebra.mul!(du::AbstractVector,L::JacVecOperator,x::AbstractVector)
    let p=L.p,t=L.t
        if L.cache1 === nothing
            if L.autodiff
                auto_jacvec!(du,(_du,_u)->L.f(_du,_u,p,t),L.u,x)
            else
                num_jacvec!(du,(_du,_u)->L.f(_du,_u,p,t),L.u,x;compute_f0=false)
            end
        else
            if L.autodiff
                auto_jacvec!(du,(_du,_u)->L.f(_du,_u,p,t),L.u,x,L.cache1,L.cache2)
            else
                num_jacvec!(du,(_du,_u)->L.f(_du,_u,p,t),L.u,x,L.cache1,L.cache2;compute_f0=false)
            end
        end
    end
end

### AnalyticalOperator Implementation

mutable struct AnalyticalJacVecOperator{T,F,uType,P,tType,O} <: DiffEqBase.AbstractDiffEqLinearOperator{T}
    f::F
    u::uType
    p::P
    t::tType
    ishermitian::Bool
    opnorm::O
    function AnalyticalJacVecOperator{T}(f,u=nothing,p=nothing,t::Union{Nothing,Number}=nothing;ishermitian=false,opnorm=true) where T
        u===nothing ? uType = Any : uType = typeof(u)
        p===nothing ? P = Any : P = typeof(p)
        t===nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),uType,P,tType,typeof(opnorm)}(f,u,p,t,ishermitian,opnorm)
    end
    function AnalyticalJacVecOperator(f,u,args...;kwargs...)
        AnalyticalJacVecOperator{eltype(u)}(f,u,args...;kwargs...)
    end
end

LinearAlgebra.opnorm(L::AnalyticalJacVecOperator, p::Real=2) = L.opnorm
LinearAlgebra.ishermitian(L::AnalyticalJacVecOperator) = L.ishermitian

Base.size(L::AnalyticalJacVecOperator) = (length(L.u),length(L.u))
Base.size(L::AnalyticalJacVecOperator,i::Int) = length(L.u)
function update_coefficients!(L::AnalyticalJacVecOperator,u,p,t)
    L.u = u
    L.p = p
    L.t = t
end

# Interpret the call as df/du * u

function (L::AnalyticalJacVecOperator)(u,p,t::Number)
    update_coefficients!(L,u,p,t)
    L*u
end

function (L::AnalyticalJacVecOperator)(du,u,p,t::Number)
    update_coefficients!(L,u,p,t)
    mul!(du,L,u)
end

Base.:*(L::AnalyticalJacVecOperator,x::AbstractVector) = L.f(x,L.u,L.p,L.t)

function LinearAlgebra.mul!(du::AbstractVector,L::AnalyticalJacVecOperator,x::AbstractVector)
    L.f(du,x,L.u,L.p,L.t)
end


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq

using Distributions: Exponential, Geometric, cdf
using Random: randsubseq
using StatsBase: sample

using Plots
using ProgressBars

include("utils.jl")

# schema
@present ThSIR(FreeSchema) begin
    (S,I,R)::Ob
end

@acset_type SIR(ThSIR)

# infection rules
I = @acset SIR begin I=1 end
I2 = @acset SIR begin I=2 end
SI = @acset SIR begin S=1; I=1 end
L_infect =  ACSetTransformation(I, SI; I=[1]) # fn from dom I to codom SI, I is the injective function mapping stuff in dom to codom
R_infect = ACSetTransformation(I, I2; I=[1])

# recovery rules
R = @acset SIR begin R=1 end
L_recovery = ACSetTransformation(SIR(), I)
R_recovery = ACSetTransformation(SIR(), R)

# parameters
N = 1000
I0 = 5
S0 = N - I0

Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5
β = R0 * γ # R0 for corresponding ODEs

state = @acset SIR begin S=S0; I=I0; R=N-S0-I0 end

# simulation loop
out = fill(-1, (steps, 3))

for t = ProgressBar(1:steps)
    # matches
    infections_m = homomorphisms(SI, state)
    recoveries_m = homomorphisms(I, state)

    # queued updates objects
    infections = queued_updates(infections_m, L_infect, R_infect)
    recoveries = queued_updates(recoveries_m, L_recovery, R_recovery)

    # sample occurances
    sample_matches(infections, β/N, Δt)
    sample_matches(recoveries, γ, Δt)

    global state = fire_events(state, [infections, recoveries])

    # write output
    out[t, :] = [nparts(state, x) for x in [:S, :I, :R]]
end

plot(
    (1:steps) * Δt,
    out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)


text:
# definition of some distribution types

## Distribution & Make

abstract type Distribution{T} end

Base.eltype(::Type{<:Distribution{T}}) where {T} = T

struct Make{T,X<:Tuple, XX<:Tuple} <: Distribution{T}
    # XX is X with Type{A} replaced by Nothing (otherwise, X (e.g. Tuple{Type{Int}}) can't have instances)
    x::XX
end

# @inline necessary for one @inferred test on arrays
@inline function Base.getindex(m::Make{T,X}, i::Integer) where {T,X}
    i = Int(i)
    fieldtype(X, i) <: Type ?
        fieldtype(X, i).parameters[1] :
        m.x[i]
end

@inline Base.getindex(m::Make, idxs::AbstractVector{<:Integer}) =
    ntuple(i->m[idxs[i]], length(idxs))
# seems faster than `Tuple(m[i] for i in idxs)`

Base.lastindex(m::Make) = lastindex(m.x)

@generated function Make{T}(X...) where T
    XX = Tuple{(x <: Type ? Nothing : x for x in X)...}
    y = [x <: Type ? nothing : :(X[$i]) for (i, x) in enumerate(X)]
    :(Make{T,Tuple{$X...},$XX}(tuple($(y...))))
end

Make0{T} = Make{T,Tuple{}}
Make1{T} = Make{T,Tuple{X}}         where X
Make2{T} = Make{T,Tuple{X,Y}}       where {X, Y}
Make3{T} = Make{T,Tuple{X,Y,Z}}     where {X, Y, Z}
Make4{T} = Make{T,Tuple{X,Y,Z,U}}   where {X, Y, Z, U}
Make5{T} = Make{T,Tuple{X,Y,Z,U,V}} where {X, Y, Z, U, V}

Make0{T}()              where {T} = Make{T}()
Make1{T}(x)             where {T} = Make{T}(x)
Make2{T}(x, y)          where {T} = Make{T}(x, y)
Make3{T}(x, y, z)       where {T} = Make{T}(x, y, z)
Make4{T}(x, y, z, u)    where {T} = Make{T}(x, y, z, u)
Make5{T}(x, y, z, u, v) where {T} = Make{T}(x, y, z, u, v)

# default maketype & make & Make(...)

# Make(...) is not meant to be specialized, i.e. Make(a, b, c) always create a Make3,
# and is equal to the *default* make(...)
# (it's a fall-back for client code which can help break recursivity)
# TODO: add tests for Make(...)

maketype(::Type{T}, x...) where {T} = T

Make(::Type{T}, x...) where {T} = Make{maketype(T, x...)}(x...)
make(::Type{T}, x...) where {T} = Make{maketype(T, x...)}(x...)

# make(x) is defined in sampling.jl, and is a special case wrapping already valid
# distributions (explicit or implicit)
Make(x1, x2, xs...) = Make{maketype(x1, x2, xs...)}(x1, x2, xs...)
make(x1, x2, xs...) = Make{maketype(x1, x2, xs...)}(x1, x2, xs...)

find_deduced_type(::Type{T}, ::X,     ) where {T,X} = deduce_type(T, gentype(X))
find_deduced_type(::Type{T}, ::Type{X}) where {T,X} = deduce_type(T, X)

find_deduced_type(::Type{T}, ::X,       ::Y)       where {T,X,Y} = deduce_type(T, gentype(X), gentype(Y))
find_deduced_type(::Type{T}, ::Type{X}, ::Y)       where {T,X,Y} = deduce_type(T, X,          gentype(Y))
find_deduced_type(::Type{T}, ::X,       ::Type{Y}) where {T,X,Y} = deduce_type(T, gentype(X), Y)
find_deduced_type(::Type{T}, ::Type{X}, ::Type{Y}) where {T,X,Y} = deduce_type(T, X,          Y)

deduce_type(::Type{T}, ::Type{X}, ::Type{Y}) where {T,X,Y} = _deduce_type(T, Val(isconcretetype(T)), X, Y)
deduce_type(::Type{T}, ::Type{X}) where {T,X} = _deduce_type(T, Val(isconcretetype(T)), X)

_deduce_type(::Type{T}, ::Val{true},  ::Type{X}, ::Type{Y}) where {T,X,Y} = T
_deduce_type(::Type{T}, ::Val{false}, ::Type{X}, ::Type{Y}) where {T,X,Y} = deduce_type(T{X}, Y)

_deduce_type(::Type{T}, ::Val{true},  ::Type{X}) where {T,X} = T
_deduce_type(::Type{T}, ::Val{false}, ::Type{X}) where {T,X} = T{X}

# show(::Make)

function Base.show(io::IO, m::Make{T}) where {T}
    M = typeof(m)
    P = M.parameters[2].parameters
    t = ntuple(length(m.x)) do i
        P[i] isa Type{<:Type} ? P[i].parameters[1] : m.x[i]
    end
    Base.show_type_name(io, M.name)
    print(io, "{", T, "}", t)
end

## Uniform

abstract type Uniform{T} <: Distribution{T} end

struct UniformType{T} <: Uniform{T} end

Uniform(::Type{T}) where {T} = UniformType{T}()

Base.getindex(::UniformType{T}) where {T} = T

struct UniformWrap{T,E} <: Uniform{E}
    val::T
end

Uniform(x::T) where {T} = UniformWrap{T,gentype(T)}(x)

Base.getindex(x::UniformWrap) = x.val

## Normal & Exponential

abstract type Normal{T} <: Distribution{T} end

struct Normal01{T} <: Normal{T} end

struct Normalμσ{T} <: Normal{T}
    μ::T
    σ::T
end

const NormalTypes = Union{AbstractFloat,Complex{<:AbstractFloat}}

Normal(::Type{T}=Float64) where {T<:NormalTypes} = Normal01{T}()
Normal(μ::T, σ::T) where {T<:NormalTypes} = Normalμσ(μ, σ)
Normal(μ::T, σ::T) where {T<:Real} = Normalμσ(AbstractFloat(μ), AbstractFloat(σ))
Normal(μ, σ) = Normal(promote(μ, σ)...)

abstract type Exponential{T} <: Distribution{T} end

struct Exponential1{T} <: Exponential{T} end

struct Exponentialθ{T} <: Exponential{T}
    θ::T
end

Exponential(::Type{T}=Float64) where {T<:AbstractFloat} = Exponential1{T}()
Exponential(θ::T) where {T<:AbstractFloat} = Exponentialθ(θ)
Exponential(θ::Real) = Exponentialθ(AbstractFloat(θ))

## floats

abstract type FloatInterval{T<:AbstractFloat} <: Uniform{T} end

abstract type CloseOpen{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenClose{ T<:AbstractFloat} <: FloatInterval{T} end
abstract type CloseClose{T<:AbstractFloat} <: FloatInterval{T} end
abstract type OpenOpen{  T<:AbstractFloat} <: FloatInterval{T} end

struct CloseOpen12{T<:AbstractFloat} <: CloseOpen{T} end # interval [1,2)

struct CloseOpen01{ T<:AbstractFloat} <: CloseOpen{T}  end # interval [0,1)
struct OpenClose01{ T<:AbstractFloat} <: OpenClose{T}  end # interval (0,1]
struct CloseClose01{T<:AbstractFloat} <: CloseClose{T} end # interval [0,1]
struct OpenOpen01{  T<:AbstractFloat} <: OpenOpen{T}   end # interval (0,1)

struct CloseOpenAB{T<:AbstractFloat} <: CloseOpen{T} # interval [a,b)
    a::T
    b::T

    CloseOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenCloseAB{T<:AbstractFloat} <: OpenClose{T} # interval (a,b]
    a::T
    b::T

    OpenCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct CloseCloseAB{T<:AbstractFloat} <: CloseClose{T} # interval [a,b]
    a::T
    b::T

    CloseCloseAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

struct OpenOpenAB{T<:AbstractFloat} <: OpenOpen{T} # interval (a,b)
    a::T
    b::T

    OpenOpenAB{T}(a::T, b::T) where {T} = (check_interval(a, b); new{T}(a, b))
end

check_interval(a, b) = a >= b && throw(ArgumentError("invalid interval specification"))

const FloatInterval_64 = FloatInterval{Float64}
const CloseOpen01_64   = CloseOpen01{Float64}
const CloseOpen12_64   = CloseOpen12{Float64}

CloseOpen01(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen12(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen12{T}()

CloseOpen(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen(a::T, b::T) where {T<:AbstractFloat} = CloseOpenAB{T}(a, b)

OpenClose(::Type{T}=Float64) where {T<:AbstractFloat} = OpenClose01{T}()
OpenClose(a::T, b::T) where {T<:AbstractFloat} = OpenCloseAB{T}(a, b)

CloseClose(::Type{T}=Float64) where {T<:AbstractFloat} = CloseClose01{T}()
CloseClose(a::T, b::T) where {T<:AbstractFloat} = CloseCloseAB{T}(a, b)

OpenOpen(::Type{T}=Float64) where {T<:AbstractFloat} = OpenOpen01{T}()
OpenOpen(a::T, b::T) where {T<:AbstractFloat} = OpenOpenAB{T}(a, b)

# convenience functions

CloseOpen(a, b) = CloseOpen(promote(a, b)...)
CloseOpen(a::T, b::T) where {T} = CloseOpen(AbstractFloat(a), AbstractFloat(b))

OpenClose(a, b) = OpenClose(promote(a, b)...)
OpenClose(a::T, b::T) where {T} = OpenClose(AbstractFloat(a), AbstractFloat(b))

CloseClose(a, b) = CloseClose(promote(a, b)...)
CloseClose(a::T, b::T) where {T} = CloseClose(AbstractFloat(a), AbstractFloat(b))

OpenOpen(a, b) = OpenOpen(promote(a, b)...)
OpenOpen(a::T, b::T) where {T} = OpenOpen(AbstractFloat(a), AbstractFloat(b))

## Bernoulli

struct Bernoulli{T<:Number} <: Distribution{T}
    p::Float64

    Bernoulli{T}(p::Real) where {T} = let pf = Float64(p)
        0.0 <= pf <= 1.0 ? new(pf) :
            throw(DomainError(p, "Bernoulli: parameter p must satisfy 0.0 <= p <= 1.0"))
    end
end

Bernoulli(p::Real=0.5) = Bernoulli(Int, p)
Bernoulli(::Type{T}, p::Real=0.5) where {T<:Number} = Bernoulli{T}(p)

## Categorical

struct Categorical{T<:Number} <: Distribution{T}
    cdf::Vector{Float64}

    function Categorical{T}(weigths) where T
        if !isa(weigths, AbstractArray)
            # necessary for accumulate
            # TODO: will not be necessary anymore in Julia 1.5
            weigths = collect(weigths)
        end
        weigths = vec(weigths)

        isempty(weigths) &&
            throw(ArgumentError("Categorical requires at least one category"))

        s = Float64(sum(weigths))
        cdf = accumulate(weigths; init=0.0) do x, y
            x + Float64(y) / s
        end
        @assert isapprox(cdf[end], 1.0) # really?
        cdf[end] = 1.0 # to be sure the algo terminates
        new{T}(cdf)
    end
end

Categorical(weigths) = Categorical{Int}(weigths)

Categorical(n::Number) =
    Categorical{typeof(n)}(Iterators.repeated(1.0 / Float64(n), Int(n)))


text:
A = wild_symbol("WildA", Sets)
B = wild_symbol("WildB", Sets)
f = wild_symbol("Wildf", Functions(Statements))

P = wild_symbol("WildP", Statements)
Q = wild_symbol("WildQ", Statements)
R = wild_symbol("WildR", Statements)

x = wild_symbol("Wildx", AΩ)
y = wild_symbol("Wildy", AΩ)
z = wild_symbol("Wildz", AΩ)

_x = symbol("_x", AΩ)
_y = symbol("_y", AΩ)

truths = ATruth[
    Axiom("Given", P, P),
]

include("Properties.jl")

include("Logic.jl")
append!(truths, logic)

include("SetTheory.jl")
append!(truths, set_theory)


text:
# Please run from the `basic_gppp` directory.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Stheno, Plots, Random, ColorTypes, FixedPointNumbers
using Stheno: @model

###########################  Define and inspect our model  ###########################

#=
g1 and g2 are time-varying basis functions. In a real application, these might be specified
by data.

w1 and w2 are the time-varying basis functions for these bases, specified by slowly
varying GPs.

f is the prediction of the regressor.

y is the addition of f and rough temporally-correlated "noise".
=#
@model function model()
    g1, g2 = x->x / 4, cos
    w1, w2 = stretch(GP(EQ()), 0.2), stretch(GP(EQ()), 1)
    f = g1 * w1 + g2 * w2
    y = f + 0.3 * GP(Matern12())
    return w1, w2, f, y
end

# Sample from the prior from plotting and for conditioning.
rng, N, Nplot, S = MersenneTwister(123456), 250, 500, 100;
X, Xp = sort(rand(rng, N) * 10), range(-2.5, stop=12.5, length=Nplot);
w1, w2, f, y = model();
w1s, w2s, fs, ŷ = rand(rng, [w1(Xp, 1e-9), w2(Xp, 1e-9), f(Xp, 1e-9), y(X, 1e-9)]);

# Compute posterior distribution over f′.
w1′, w2′, f′, y′ = (w1, w2, f, y) | (y(X) ← ŷ);

# Sample from the posterior and write to file.
w1′s, w2′s, y′s = rand(rng, [w1′(Xp, 1e-9), w2′(Xp, 1e-9), y′(Xp, 1e-9)], S);

# Get posterior mean and marginals f′ and y′ and write them for plotting.

ms_w1′ = marginals(w1′(Xp));
ms_w2′ = marginals(w2′(Xp));
ms_y′ = marginals(y′(Xp));

μw1′, σw1′ = mean.(ms_w1′), std.(ms_w1′);
μw2′, σw2′ = mean.(ms_w2′), std.(ms_w2′);
μy′, σy′ = mean.(ms_y′), std.(ms_y′);

###########################  Plot results ###########################

gr();
posterior_plot = plot(
    legend=:topright,
    legendfont=Plots.Font(
        "sans-serif",
        10,
        :hcenter,
        :vcenter,
        0.0,
        RGB{Normed{UInt8, 8}}(0.0,0.0,0.0)
    ),
    background_color_legend=RGBA(1, 1, 1, 0),
    foreground_color_legend=RGBA(1, 1, 1, 0),
);

# Plot posterior over w1.
plot!(posterior_plot, Xp, μw1′;
    linecolor=:green,
    linewidth=2.0,
    label="w1",
);
plot!(posterior_plot, Xp, [μw1′ μw1′];
    linewidth=0.0,
    fillrange=[μw1′ .- 3 .* σw1′, μw1′ .+ 3 * σw1′],
    fillalpha=0.2,
    fillcolor=:green,
    label="",
);
plot!(posterior_plot, Xp, w1′s;
    linecolor=:green,
    linealpha=0.1,
    label="",
);

# Plot posterior over w2.
plot!(posterior_plot, Xp, μw2′;
    linecolor=:magenta,
    linewidth=2.0,
    label="w2",
);
plot!(posterior_plot, Xp, [μw2′ μw2′];
    linewidth=0.0,
    fillrange=[μw2′ .- 3 .* σw2′, μw2′ .+ 3 * σw2′],
    fillalpha=0.2,
    fillcolor=:magenta,
    label="",
);
plot!(posterior_plot, Xp, w2′s;
    linecolor=:magenta,
    linealpha=0.1,
    label="",
);

# Plot x1 and x2
plot!(posterior_plot, Xp, (x->x / 4).(Xp);
    linecolor=:black,
    linewidth=1.0,
    label="sin",
);
plot!(posterior_plot, Xp, cos.(Xp);
    linecolor=:black,
    linewidth=1.0,
    linestyle=:dash,
    label="cos",
);

# Plot samples against which we're regressing.
scatter!(posterior_plot, X, ŷ;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

# Plot posterior over `y`.
plot!(posterior_plot, Xp, μy′;
    linecolor=:blue,
    linewidth=2.0,
    label="y",
);
plot!(posterior_plot, Xp, [μy′ μy′];
    linewidth=0.0,
    fillrange=[μy′ .- 3 .* σy′, μy′ .+ 3 * σy′],
    fillalpha=0.3,
    fillcolor=:blue,
    label="",
);

# savefig(posterior_plot, "plotting_research_talk/tv_blr.pdf");

display(posterior_plot);


text:
export first_permuted_prime
export isprime
export sieve_eratosthenes

"""
Creates permutations of numbers and returns when the first prime value is reached.
"""
function first_permuted_prime(elementlist::Array, multiplier::Array, k::Integer=1)
    if k == length(elementlist)
        if elementlist[end] % 2 == 1
            val = sum(elementlist .* multiplier)
            if isprime(val)
                return val
            else
                return -1
            end
        else
            return -1
        end
    else
        for i in k:length(elementlist)
            elementlist[k], elementlist[i] = elementlist[i], elementlist[k]
            v = first_permuted_prime(elementlist, multiplier, k + 1)
            if v > 0
                return v
            end
            elementlist[k], elementlist[i] = elementlist[i], elementlist[k]
        end
    end
    return -1
end

"""
Reasonably fast test for primality.
"""
function isprime(n)
    if n == 2 || n == 3
        return true
    end
    if n < 2 || n % 2 == 0
        return false
    end
    if n < 9
        return true
    end
    if n % 3 == 0
        return false
    end
    r = typeof(n)(floor(sqrt(n)))
    f = 5
    while f <= r
        if n % f == 0
            return false
        end
        if n % (f + 2) == 0
            return false
        end
        f += 6
    end
    return true
end

"""
Implements the Sieve of Eratosthenes to return all prime numbers up "n".
https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
"""
function sieve_eratosthenes(n::Integer=500, nstart::Integer=2)::Array{Integer}
    primes = trues(n)
    primes[1] = false

    for i in 2:floor(Integer, sqrt(n))
        if primes[i]
            @inbounds primes[(i+i):i:n] .= false
        end
    end
    primes[1:nstart-1] .= false
    return findall(primes)
end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
export TCPSCPIInstrument

mutable struct TCPSCPIInstrument <: SCPIInstrument
    socket::TCPSocket
    is_connected::Bool
    delimiter::String

    function TCPSCPIInstrument(ip::IPAddr, port::Integer)
        return new(connect(ip, port), true, "\n")
    end
end

"""
Constructor interpreting the given `ip` string as a IPv4 address
"""
TCPSCPIInstrument(ip::String, port::Integer) = TCPSCPIInstrument(IPv4(ip), port)

function Base.close(inst::TCPSCPIInstrument)
    close(inst.socket)
    inst.is_connected = false
end

"""
Send a command to the instrument
"""
function command(inst::TCPSCPIInstrument, cmd::String)
  write(inst.socket, cmd*inst.delimiter)
end

"""
Perform a query to the instrument. Return String.
"""
function query(inst::TCPSCPIInstrument, cmd::String)
  command(inst, cmd)
  return readline(inst.socket)[1:end]
end

"""
Perform a query to the iunstrument. Parse result as type T.
"""
function query(inst::TCPSCPIInstrument, cmd::String, T::Type)
  a = query(inst, cmd)
  return parse(T, a)
end


text:
module FITSIO

export FITS,
       HDU,
       ImageHDU,
       TableHDU,
       ASCIITableHDU,
       FITSHeader,
       read_key,
       write_key,
       read_header,
       get_comment,
       set_comment!,
       copy_section,
       default_header,
       deletekey!

import Base: getindex,
             setindex!,
             length,
             show,
             read,
             read!,
             write,
             close,
             ndims,
             size,
             getkey,
             haskey,
             keys,
             values,
             eltype,
             deleteat!

# Deal with compatibility issues.
using Printf
import Base: iterate, lastindex
# Libcfitsio submodule

## Have to manually import while deprecating `libcfitsio_version`.
## After removing that, can return to using CFITSIO
import CFITSIO: FITSFile,
                FITSMemoryHandle,
                fits_open_file,
                fits_open_diskfile,
                fits_create_file,
                fits_create_diskfile,
                fits_assert_open,
                fits_file_mode,
                fits_create_img,
                fits_close_file,
                fits_write_pix,
                fits_get_num_hdus,
                fits_movabs_hdu,
                fits_delete_hdu,
                fits_get_img_size,
                fits_get_img_type,
                fits_get_img_equivtype,
                type_from_bitpix,
                fits_read_pix,
                fits_read_subset,
                fits_copy_image_section,
                fits_file_name,
                fits_get_img_dim,
                fits_read_tdim,
                fits_write_tdim,
                fits_read_col,
                fits_write_col,
                fits_get_num_cols,
                fits_get_num_rows,
                fits_get_colnum,
                fits_get_eqcoltype,
                fits_read_descript,
                fits_update_key,
                fits_write_comment,
                fits_write_history,
                fits_get_hdrspace,
                fits_hdr2str,
                fits_read_keyword,
                fits_read_keyn,
                fits_open_memfile,
                fits_movnam_hdu,
                fits_get_hdu_num,
                fits_get_hdu_type

# There are a few direct `ccall`s to libcfitsio in this module. For this, we
# need a few non-exported things from Libcfitsio: the shared library handle,
# and a helper function for raising errors.
import CFITSIO: libcfitsio,
                fits_assert_ok,
                fits_assert_isascii

import CFITSIO
@deprecate libcfitsio_version CFITSIO.libcfitsio_version

import Tables

## DEPRECATED
module Libcfitsio
    using Reexport
    @reexport using CFITSIO
end

# HDU Types
abstract type HDU end

mutable struct ImageHDU{T<:Real,N} <: HDU
    fitsfile::FITSFile
    ext::Int

    function ImageHDU(fitsfile::FITSFile, ext::Integer)
        fits_assert_open(fitsfile)
        fits_movabs_hdu(fitsfile, ext)
        N = Int(fits_get_img_dim(fitsfile))
        bitpix = fits_get_img_equivtype(fitsfile)
        T = type_from_bitpix(bitpix)
        new{T,N}(fitsfile, ext)
    end
end

mutable struct TableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

mutable struct ASCIITableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

# FITS
#
# The FITS type represents a FITS file. It holds a reference to a
# FITSFile object (basically the low-level CFITSIO pointer). It also
# holds a reference to all of the previously accessed HDU
# objects. This is so that only a single HDU object is created for
# each extension in the file. It also allows a FITS object to tell
# previously created HDUs about events that happen to the file, such
# as deleting extensions. This could be done by, e.g., setting ext=-1
# in the HDU object.
"""
    FITS(filename::String[, mode::String = "r"]; extendedparser = true)

Open or create a FITS file. `mode` can be one of `"r"` (read-only),
`"r+"` (read-write) or `"w"` (write). In "write" mode, any existing
file of the same name is overwritten.

A `FITS` object is a collection of "Header-Data Units" (HDUs) and
supports the following operations:

* `f[i]`: Return the `i`-th HDU.

* `f[name]` or `f[name, ver]`: Return the HDU containing the given the
  given EXTNAME (or HDUNAME) keyword (a String), and optionally the
  given EXTVER (or HDUVER) number (an Integer).

* Iteration:
  ```
  for hdu in f
      ...
  end
  ```

The keyword argument `extendedparser` may be used to enable or disable the
[extended filename parser](https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node83.html).
If disabled, `filename` is treated exactly as the name of the file and is not tokenized into
parameters.
"""
mutable struct FITS
    fitsfile::FITSFile
    filename::AbstractString
    mode::AbstractString
    hdus::Dict{Int, HDU}

    # Hold on to memory if backed by a julia buffer
    memhandle::FITSMemoryHandle
    hidden::Any

    function FITS(filename::AbstractString, mode::AbstractString="r"; extendedparser = true)
        openfn = extendedparser ? fits_open_file : fits_open_diskfile
        createfn = extendedparser ? fits_create_file : fits_create_diskfile
        f = (mode == "r"                      ? openfn(filename, 0)    :
             mode == "r+" && isfile(filename) ? openfn(filename, 1)    :
             mode == "r+"                     ? createfn(filename)     :
             mode == "w"                      ?
                (rm(filename, force = true); createfn(filename)) :
             error("invalid open mode: $mode"))

        new(f, filename, mode, Dict{Int, HDU}(), FITSMemoryHandle(), nothing)
    end

    function FITS(data::Vector{UInt8}, mode::AbstractString="r", filename = "")
        @assert mode == "r"
        f, handle = fits_open_memfile(data, 0)
        new(f, filename, mode, Dict{Int, HDU}(), handle, data)
    end
end

"""
    FITS(f::Function, args...; kwargs...)

Apply the function `f` to the result of `FITS(args...; kwargs...)` and close the
resulting file descriptor upon completion.
"""
function FITS(f::Function, args...; kwargs...)
    io = FITS(args...; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

"""
    FITSHeader(keys::Vector{String}, values::Vector, comments::Vector{String})

An in-memory representation of the header of an HDU. It stores the
(key, value, comment) information for each 80-character "card" in a header.

Note that this structure is not linked to a FITS file in any way; it
is just a convenient structure for storing the header contents after
reading from a file. (This is similar to how an `Array` returned by
`read(f[1])` is not linked to the FITS file `f`.)  Manipulating a
`FITSHeader` will therefore have no immediate impact on any file, even
if it was created by `read_header(::HDU)`.  You can, however, write a
`FITSHeader` to a file using the `write(::FITS, ...)` methods that
append a new HDU to a file.
"""
mutable struct FITSHeader
    keys::Vector{String}
    values::Vector{Any}
    comments::Vector{String}
    map::Dict{String, Int}

    function FITSHeader(keys::Vector{String}, values::Vector,
                        comments::Vector{String})
        if ((length(keys) != length(values)) ||
            (length(keys) != length(comments)))
            error("keys, values, comments must be same length")
        end
        map = Dict{String, Int}()
        for i in 1:length(keys)
          map[keys[i]] = i
        end
        new(keys, Vector{Any}(values), comments, map)
    end
end

include("fits.jl")  # FITS methods
include("header.jl")  # FITSHeader methods
include("image.jl")  # ImageHDU methods
include("table.jl")  # TableHDU & ASCIITableHDU methods

end # module


text:
module TestMCTools

using QuantEcon
using Base.Test
using FactCheck
srand(42)

# these matrices come from RMT4 section 2.2.1
mc1 = [1 0 0; .2 .5 .3; 0 0 1]
mc2 = [.7 .3 0; 0 .5 .5; 0 .9 .1]
mc3 = [0.4 0.6; 0.2 0.8]
mc4 = eye(2)
mc5 = [
     0 1 0 0 0 0
     1 0 0 0 0 0
     1//2 0 0 1//2 0 0
     0 0 0 0 1 0
     0 0 0 0 0 1
     0 0 0 1 0 0
     ]
mc5_stationary = zeros(Rational,6,2)
mc5_stationary[[1,2]] = 1//2; mc5_stationary[[10,11,12]] = 1//3

mc6 = [2//3 1//3; 1//4 3//4]  # Rational elements
mc6_stationary = [3//7, 4//7]

mc7 = [1 0; 0 1]
mc7_stationary = [1 0;0 1]

# Reducible mc with a unique recurrent class,
# where n=2 is a transient state
mc10 = [1. 0; 1. 0]
mc10_stationary = [1., 0]

mc1 = MarkovChain(mc1)
mc2 = MarkovChain(mc2)
mc3 = MarkovChain(mc3)
mc4 = MarkovChain(mc4)
mc5 = MarkovChain(mc5)
mc6 = MarkovChain(mc6)
mc7 = MarkovChain(mc7)
mc10 = MarkovChain(mc10)

# examples from
# Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

fig1_p = zeros(Rational{Int64}, 5, 5)
fig1_p[[3, 4, 9, 10, 11, 13, 18, 19, 22, 24]] =
    [1//2, 2//5, 1//10, 1, 1, 1//5, 3//10, 1//5, 1, 3//10]
fig2_p = zeros(Rational{Int64}, 5, 5)
fig2_p[[3, 10, 11, 13, 14, 17, 18, 19, 22]] =
    [1//3, 1, 1, 1//3, 1//2, 1//2, 1//3, 1//2, 1//2]

fig1 = MarkovChain(convert(Matrix{Float64}, fig1_p))
fig1_rat = MarkovChain(fig1_p)

fig2 = MarkovChain(convert(Matrix{Float64}, fig2_p))
fig2_rat = MarkovChain(fig2_p)

function kmr_markov_matrix_sequential{T<:Real}(n::Integer, p::T, ε::T)
    """
    Generate the MarkovChain with the associated transition matrix from the KMR model with *sequential* move

    n: number of players
    p: level of p-dominance for action 1
       = the value of p such that action 1 is the BR for (1-q, q) for any q > p,
         where q (1-q, resp.) is the prob that the opponent plays action 1 (0, resp.)
    ε: mutation probability

    References:
        KMRMarkovMatrixSequential is contributed from https://github.com/oyamad
    """
    x = zeros(T, n+1, n+1)

    x[1, 1], x[1, 2] = 1 - ε/2, ε/2
    @inbounds for i = 1:n-1
        x[i+1, i] = (i/n) * (ε/2 + (1 - ε) *
                             (((i-1)/(n-1) < p) + ((i-1)/(n-1) == p)/2))
        x[i+1, i+2] = ((n-i)/n) * (ε/2 + (1 - ε) *
                                 ((i/(n-1) > p) + (i/(n-1) == p)/2))
        x[i+1, i+1] = 1 - x[i+1, i] - x[i+1, i+2]
    end
    x[end, end-1], x[end, end] = ε/2, 1 - ε/2
    return MarkovChain(x)
end

mc8 = kmr_markov_matrix_sequential(27, 1/3, 1e-2)
mc9 = kmr_markov_matrix_sequential(3, 1/3, 1e-14)

tol = 1e-15

facts("Testing mc_tools.jl") do
    context("test mc_compute_stationary using exact solutions") do
        @fact mc_compute_stationary(mc1) --> eye(3)[:, [1, 3]]
        @fact mc_compute_stationary(mc2) --> roughly([0, 9/14, 5/14])
        @fact mc_compute_stationary(mc3) --> roughly([1/4, 3/4])
        @fact mc_compute_stationary(mc4) --> eye(2)
        @fact mc_compute_stationary(mc5) --> mc5_stationary
        @fact mc_compute_stationary(mc6) --> mc6_stationary
        @fact mc_compute_stationary(mc7) --> mc7_stationary
        @fact mc_compute_stationary(mc10) --> mc10_stationary
    end

    context("test gth_solve with KMR matrices") do
        for d in [mc8,mc9]
            x = mc_compute_stationary(d)

            # Check elements sum to one
            @fact sum(x) --> roughly(1; atol=tol)

            # Check elements are nonnegative
            for i in 1:length(x)
                @fact x[i] --> greater_than_or_equal(-tol)
            end

            # Check x is a left eigenvector of P
            @fact vec(x'*d.p) --> roughly(x; atol=tol)
        end
    end

    context("test MarkovChain throws errors") do
        @fact_throws MarkovChain(rand(4, 5))  # not square
        @fact_throws MarkovChain([0.0 0.5; 0.2 0.8])  # first row doesn't sum to 1
        @fact_throws MarkovChain([-1 1; 0.2 0.8])  # negative element, but sums to 1
    end

    context("test graph theoretic algorithms") do
        for fig in [fig1, fig1_rat]
            @fact recurrent_classes(fig) --> Vector[[2, 5]]
            @fact communication_classes(fig) --> Vector[[2, 5], [1, 3, 4]]
            @fact is_irreducible(fig) --> false
            # @fact period(fig) --> 2
            # @fact is_aperiodic(fig) --> false
        end

        for fig in [fig2, fig2_rat]
            @fact recurrent_classes(fig) --> Vector[[1, 3, 4]]
            @fact communication_classes(fig) --> Vector[[1, 3, 4], [2, 5]]
            @fact is_irreducible(fig) --> false
            # @fact period(fig) --> 1
            # @fact is_aperiodic(fig) --> true
        end
    end
end  # facts

end  # module


text:
## Triangular
immutable Triangular{T<:Number} <: AbstractMatrix{T}
    UL::Matrix{T}
    uplo::Char
    unitdiag::Char
end
function Triangular{T<:Number}(A::Matrix{T}, uplo::Symbol, unitdiag::Bool)
    if size(A, 1) != size(A, 2) throw(DimensionMismatch("matrix must be square")) end
    return Triangular(A, string(uplo)[1], unitdiag ? 'U' : 'N')
end
Triangular(A::Matrix, uplo::Symbol) = Triangular(A, uplo, all(diag(A) .== 1) ? true : false)
function Triangular(A::Matrix)
    if istriu(A) return Triangular(A, :U) end
    if istril(A) return Triangular(A, :L) end
    throw(ArgumentError("matrix is not triangular"))
end

######################
# BlasFloat routines #
######################

# Vector multiplication
*{T<:BlasFloat}(A::Triangular{T}, b::Vector{T}) = BLAS.trmv(A.uplo, 'N', A.unitdiag, A.UL, b)
Ac_mul_B{T<:BlasComplex}(A::Triangular{T}, b::Vector{T}) = BLAS.trmv(A.uplo, 'C', A.unitdiag, A.UL, b)
At_mul_B{T<:BlasReal}(A::Triangular{T}, b::Vector{T}) = BLAS.trmv(A.uplo, 'T', A.unitdiag, A.UL, b)

# Matrix multiplication
*{T<:BlasFloat}(A::Triangular{T}, B::StridedMatrix{T}) = BLAS.trmm('L', A.uplo, 'N', A.unitdiag, one(T), A.UL, B)
*{T<:BlasFloat}(A::StridedMatrix{T}, B::Triangular{T}) = BLAS.trmm('R', B.uplo, 'N', B.unitdiag, one(T), A, B.UL)
A_mul_B!{T<:BlasFloat}(A::Triangular{T},B::Matrix{T}) = BLAS.trmm!('L',A.uplo,'N',A.unitdiag,one(eltype(A)),A.UL,B)
Ac_mul_B!{T<:BlasComplex}(A::Triangular{T}, B::StridedMatrix{T}) = BLAS.trmm('L', A.uplo, 'C', A.unitdiag, one(T), A.UL, B)
Ac_mul_B!{T<:BlasReal}(A::Triangular{T}, B::StridedMatrix{T}) = BLAS.trmm('L', A.uplo, 'T', A.unitdiag, one(T), A.UL, B)
A_mul_Bc!{T<:BlasComplex}(A::StridedMatrix{T}, B::Triangular{T}) = BLAS.trmm('R', B.uplo, 'C', B.unitdiag, one(T), B.UL, A)
A_mul_Bc!{T<:BlasReal}(A::StridedMatrix{T}, B::Triangular{T}) = BLAS.trmm('R', B.uplo, 'T', B.unitdiag, one(T), B.UL, A)

function \{T<:BlasFloat}(A::Triangular{T}, B::StridedVecOrMat{T})
    x = LAPACK.trtrs!(A.uplo, 'N', A.unitdiag, A.UL, copy(B))
    errors=LAPACK.trrfs!(A.uplo, 'N', A.unitdiag, A.UL, B, x)
    all(isfinite, [errors...]) || all([errors...] .< one(T)/eps(T)) || warn("""Unreasonably large error in computed solution:
forward errors:
$(errors[1])
backward errors:
$(errors[2])""")
    x
end
Ac_ldiv_B{T<:BlasReal}(A::Triangular{T}, B::StridedVecOrMat{T}) = LAPACK.trtrs!(A.uplo, 'T', A.unitdiag, A.UL, copy(B))
Ac_ldiv_B{T<:BlasComplex}(A::Triangular{T}, B::StridedVecOrMat{T}) = LAPACK.trtrs!(A.uplo, 'C', A.unitdiag, A.UL, copy(B))

/{T<:BlasFloat}(A::StridedVecOrMat{T}, B::Triangular{T}) = BLAS.trsm!('R', B.uplo, 'N', B.unitdiag, one(T), B.UL, copy(A))
A_rdiv_Bc{T<:BlasReal}(A::StridedVecOrMat{T}, B::Triangular{T}) = BLAS.trsm!('R', B.uplo, 'T', B.unitdiag, one(T), B.UL, copy(A))
A_rdiv_Bc{T<:BlasComplex}(A::StridedVecOrMat{T}, B::Triangular{T}) = BLAS.trsm!('R', B.uplo, 'C', B.unitdiag, one(T), B.UL, copy(A))

inv{T<:BlasFloat}(A::Triangular{T}) = LAPACK.trtri!(A.uplo, A.unitdiag, copy(A.UL))

#Eigensystems
function eigvecs{T<:BlasFloat}(A::Triangular{T})
    if A.uplo=='U'
        V = LAPACK.trevc!('R', 'A', Array(Bool,1), A.UL)
    else #A.uplo=='L'
        V = LAPACK.trevc!('L', 'A', Array(Bool,1), A.UL')
    end
    for i=1:size(V,2) #Normalize
        V[:,i] /= norm(V[:,i])
    end
    V
end

function cond{T<:BlasFloat}(A::Triangular{T}, p::Real=2)
    chksquare(A)
    if p==1
        return inv(LAPACK.trcon!('O', A.uplo, A.unitdiag, A.UL))
    elseif p==Inf
        return inv(LAPACK.trcon!('I', A.uplo, A.unitdiag, A.UL))
    else #use fallback
        return cond(full(A), p)
    end
end

####################
# Generic routines #
####################

size(A::Triangular, args...) = size(A.UL, args...)

convert(::Type{Matrix}, A::Triangular) = full(A)
convert{T}(::Type{Triangular{T}}, A::Triangular{T}) = A
convert{T}(::Type{Triangular{T}}, A::Triangular) = Triangular(convert(AbstractMatrix{T}, A.UL), A.uplo, A.unitdiag)
convert{T}(::Type{AbstractMatrix{T}}, A::Triangular) = convert(Triangular{T}, A)

full(A::Triangular) = A.uplo == 'U' ? triu!(A.UL) : tril!(A.UL)

getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? (A.unitdiag == 'U' ? one(T) : A.UL[i,j]) : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T))

istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL)
istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL)

transpose(A::Triangular) = Triangular(transpose(A.UL), A.uplo=='U'?'L':'U', A.unitdiag)
ctranspose(A::Triangular) = Triangular(ctranspose(A.UL), A.uplo=='U'?'L':'U', A.unitdiag)
transpose!(A::Triangular) = Triangular(copytri!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag)
ctranspose!(A::Triangular) = Triangular(copytri!(A.UL, A.uplo, true), A.uplo=='U'?'L':'U', A.unitdiag)
diag(A::Triangular) = diag(A.UL)
big(A::Triangular) = Triangular(big(A.UL), A.uplo, A.unitdiag)

#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
    @eval begin
        ($func){T}(A::Triangular{T}, B::AbstractVector{T}) = ($func)(full(A), B)
        #($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B)
    end
end

function sqrtm{T}(A::Triangular{T})
    n = size(A, 1)
    R = zeros(T, n, n)
    if A.uplo == 'U'
        for j = 1:n
            (T<:Complex || A[j,j]>=0) ? (R[j,j]=sqrt(A[j,j])) : throw(SingularException(j))
            for i = j-1:-1:1
                r = A[i,j]
                for k = i+1:j-1
                    r -= R[i,k]*R[k,j]
                end
                r==0 || (R[i,j] = r / (R[i,i] + R[j,j]))
            end
        end
        return Triangular(R)
    else #A.uplo == 'L' #Not the usual case
        return sqrtm(A.').'
    end
end

#Generic solver using naive substitution
function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b)
    N = size(A, 2)
    N==length(b)==length(x) || throw(DimensionMismatch(""))

    if A.uplo == 'L' #do forward substitution
        for j = 1:N
            x[j] = b[j]
            for k = 1:j-1
                x[j] -= A[j,k] * x[k]
            end
            if A.unitdiag=='N'
                x[j]/= A[j,j]==0 ? throw(SingularException(j)) : A[j,j]
            end
        end
    elseif A.uplo == 'U' #do backward substitution
        for j = N:-1:1
            x[j] = b[j]
            for k = j+1:1:N
                x[j] -= A[j,k] * x[k]
            end
            if A.unitdiag=='N'
                x[j]/= A[j,j]==0 ? throw(SingularException(j)) : A[j,j]
            end
        end
    else
        throw(ArgumentError("Unknown uplo=$(A.uplo)"))
    end
    x
end

#Generic eigensystems
eigvals(A::Triangular) = diag(A.UL)
det(A::Triangular) = prod(eigvals(A))

function eigvecs{T}(A::Triangular{T})
    evecs = zeros(A)
    N = size(A,1)
    if A.unitdiag == 'U' #Trivial
        return eye(A)
    elseif A.uplo == 'L' #do forward substitution
        for i=1:N
            evecs[i,i] = one(T)
            for j = i+1:N
                for k = i:j-1
                    evecs[j,i] -= A[j,k] * evecs[k,i]
                end
                evecs[j,i] /= A[j,j]-A[i,i]
            end
            evecs[i:N, i] /= norm(evecs[i:N, i])
        end
        evecs
    elseif A.uplo == 'U' #do backward substitution
        for i=1:N
            evecs[i,i] = one(T)
            for j = i-1:-1:1
                for k = j+1:i
                    evecs[j,i] -= A[j,k] * evecs[k,i]
                end
                evecs[j,i] /= A[j,j]-A[i,i]
            end
            evecs[1:i, i] /= norm(evecs[1:i, i])
        end
    else
        throw(ArgumentError("Unknown uplo=$(A.uplo)"))
    end
    evecs
end

eigfact(A::Triangular) = Eigen(eigvals(A), eigvecs(A))

#Generic singular systems
for func in (:svd, :svdfact, :svdfact!, :svdvals, :svdvecs)
    @eval begin
        ($func)(A::Triangular) = ($func)(full(A))
    end
end

max_stars_repo_name: Maple-Wang/julia

max_issues_repo_name: Maple-Wang/julia

max_forks_repo_name: Maple-Wang/julia