text:
"""
Describes a single water reservoir, using volumes (m^3) and discharges (m^3/s).
It is linked to a series of input rivers that provied water, and has to
fulfil some requirements with the output (tap or industry water, environmental
flow, with or without uncertainty).

TODO:

  * what about topology? (where inputs and diversions can be other reservoirs!).
    Store the rivers inside the basin, and only use names here? Could also get
    rid of different arrays for rivers. Could also store information about the links
    to improve the modelling (e.g., add the wave propagation along this edge).

  * consider uncertainty about purposes (like drinking water consumption)?
"""
immutable Reservoir
  name::AbstractString
  variant::AbstractString
  purposes::Array{Purpose, 1}
  hydropower::Array{HydropowerUnit, 1}
  outputs::Array{DamOutput, 1}

  # Reservoir bathymetry.
  capacity::Tuple{Float64, Float64} # 10^6 m^3; minimum is enough to ensure drinking water intake, by hypothesis.
  bathymetry::ReservoirBathymetry

  # Topography.
  rivers_in::Array{NaturalRiver, 1}
  rivers_diverted::Array{DivertedRiver, 1}

  # Private constructor, without units nor defaut values.
  function Reservoir(name::AbstractString, variant::AbstractString, purposes::Array{Purpose, 1},
                     hydropower::Array{HydropowerUnit, 1}, outputs::Array{DamOutput, 1},
                     capacity::Tuple{Float64, Float64}, bathymetry::ReservoirBathymetry,
                     rivers_in::Array{NaturalRiver, 1}, rivers_diverted::Array{DivertedRiver, 1}
                    )
    ## Some consistency check for the inputs.
    if length(rivers_in) + length(rivers_diverted) == 0
      error("Must have at least one input river!")
    end

    ## If there is hydropower, complement the purposes and outputs.
    # First check if these hydropower units are notalready in the given arrays!
    purposeHPs = HydropowerUnit[getHydropowerUnit(p) for p in filter(p -> isHydropower(p), purposes)]
    outputHPs = HydropowerUnit[getHydropowerUnit(o) for o in filter(o -> isHydropower(o), outputs)]
    for hp in hydropower
      if ! in(hp, purposeHPs)
        push!(purposes, DeterministicHydropower(hp))
      end
      if ! in(hp, outputHPs)
        push!(outputs, HydropowerDamOutput(hp))
      end
    end

    ## Build the object.
    obj = new(name, variant, purposes, hydropower, outputs,
              capacity, bathymetry,
              rivers_in, rivers_diverted)

    ## Check consistency.
    # Check consistency of the bathymetry.
    if capacity[1] > capacity[2]
      error("Inconsistent capacities: the minimum is higher than the maximum!")
    end

    # Check consistency of scenarios.
    countScenariosWithConsistency(obj)
    countTimeStepsWithConsistency(obj)
    getPeriodWithConsistency(obj)
    getTimePointsWithConsistency(obj)

    # Check consistency of the ids.
    ids = vec(vcat([getId(p) for p in purposes], [getId(hp) for hp in hydropower], [getId(o) for o in outputs]))
    ids_cnt = filter((k,v) -> v > 1, countmap(ids))
    if length(ids_cnt) > 0
      error("IDs are not consistent! Already met: " * string(ids_cnt) * ".")
    end

    ## Done!
    return obj
  end
end

# Public constructor, with units and keyword arguments.
Reservoir(;name::AbstractString=error("Missing reservoir name"),
           variant::AbstractString="vanilla",
           purposes::Array{Purpose, 1}=Purpose[],
           hydropower::Array{HydropowerUnit, 1}=HydropowerUnit[],
           outputs::Array{DamOutput, 1}=DamOutput[],

           capacity::Tuple{Volume, Volume}=error("Missing reservoir capacities"),
           bathymetry::ReservoirBathymetry=NullReservoirBathymetry(),

           rivers_in::Array{NaturalRiver, 1}=NaturalRiver[],
           rivers_diverted::Array{DivertedRiver, 1}=DivertedRiver[]
          ) = Reservoir(name, variant, purposes, hydropower, outputs,
                        _from_unitful(capacity), bathymetry,
                        rivers_in, rivers_diverted)

copy(r::Reservoir; name::AbstractString=r.name,
                   variant::AbstractString=r.variant,
                   purposes::Array{Purpose, 1}=r.purposes,
                   hydropower::Array{HydropowerUnit, 1}=r.hydropower,
                   outputs::Array{DamOutput, 1}=r.outputs,
                   capacity::Tuple{Union{Float64, Volume}, Union{Float64, Volume}}=r.capacity,
                   bathymetry::ReservoirBathymetry=r.bathymetry,
                   rivers_in::Array{NaturalRiver, 1}=r.rivers_in,
                   rivers_diverted::Array{DivertedRiver, 1}=r.rivers_diverted
    ) =
  Reservoir(name, variant, purposes, hydropower, outputs,
            _from_unitful(capacity), bathymetry,
            rivers_in, rivers_diverted)

Base.isempty(r::Reservoir) = false

### Generic getters.

getName(r::Reservoir) = r.name
getVariant(r::Reservoir) = r.variant
getMinimumCapacity(r::Reservoir) = r.capacity[1]
getMaximumCapacity(r::Reservoir) = r.capacity[2]

getHydropowerUnits(r::Reservoir) = r.hydropower
hasHydropower(r::Reservoir) = length(getHydropowerUnits(r)) > 0
getHydropowerIds(r::Reservoir) = [getId(hp) for hp in getHydropowerUnits(r)]
getPenstockHydropower(r::Reservoir) = sum([hp.max_discharge for hp in getHydropowerUnits(r)])

getPurposes(r::Reservoir) = r.purposes
getOutputs(r::Reservoir) = r.outputs
getBathymetry(r::Reservoir) = r.bathymetry

function _check_purpose_subtype(purposeType::DataType)
  if ! (purposeType <: Purpose)
    error("Wrong type argument: must be a subtype of Purpose; got " * string(purposeType) * ".")
  end
end
function hasPurpose(r::Reservoir, purposeType::DataType)
  _check_purpose_subtype(purposeType)
  return length(filter(t -> isa(t, purposeType), getPurposes(r))) > 0
end
function getPurpose(r::Reservoir, purposeType::DataType)
  _check_purpose_subtype(purposeType)
  # Return an array of purposes, as the type has no unicity guarantee (as opposed to IDs).
  return filter(t -> isa(t, purposeType), getPurposes(r))
end
getPurposeIds(r::Reservoir) = Symbol[getId(o) for o in getPurposes(r)]
hasPurpose(r::Reservoir, s::Symbol) = length(filter(o -> getId(o) == s, getPurposes(r))) > 0
function getPurpose(r::Reservoir, s::Symbol; nullable::Bool=false)
  elt = filter(o -> getId(o) == s, getPurposes(r))
  if nullable
    if length(elt) == 0
      return Nullable{Symbol}()
    else
      return Nullable(elt[1])
    end
  else
    if length(elt) == 0
      error("Purpose " * string(s) * " not found.")
    else
      return elt[1]
    end
  end
end

getOutputIds(r::Reservoir) = Symbol[getId(o) for o in getOutputs(r)]
hasOutput(r::Reservoir, s::Symbol) = length(filter(o -> getId(o) == s, getOutputs(r))) > 0
function getOutput(r::Reservoir, s::Symbol; nullable::Bool=false)
  elt = filter(o -> getId(o) == s, getOutputs(r))
  if nullable
    if length(elt) == 0
      return Nullable{Symbol}()
    else
      return Nullable(elt[1])
    end
  else
    if length(elt) == 0
      error("Purpose " * string(s) * " not found.")
    else
      return elt[1]
    end
  end
end

name(r::Reservoir) = getName(r)
variant(r::Reservoir) = getVariant(r)
minCapacity(r::Reservoir) = getMinimumCapacity(r)
maxCapacity(r::Reservoir) = getMaximumCapacity(r)
penstockHydropower(r::Reservoir) = getPenstockHydropower(r)
hydropowerUnits(r::Reservoir) = getHydropowerUnits(r)
outputs(r::Reservoir) = getOutputs(r)
outputsId(r::Reservoir) = getOutputIds(r)
output(r::Reservoir, s::Symbol; kwargs...) = getOutput(r, s; kwargs...)
bathymetry(r::Reservoir) = getBathymetry(r)
purpose(r::Reservoir, s::Symbol; kwargs...) = getPurpose(r, s; kwargs...)
purpose(r::Reservoir, purposeType::DataType) = getPurpose(r, purposeType)

# Counting things.

"""
Counts the number of scenarios within this reservoir.
"""
countScenarios(r::Reservoir) = countScenarios(first(rivers(r)))

"""
Counts the number of scenarios within this reservoir and ensures the whole data structure is consistent.
"""
function countScenariosWithConsistency(r::Reservoir)
  cnt = -1

  for river in rivers(r)
    if cnt < 0
      cnt = countScenariosWithConsistency(river)
    elseif cnt != countScenariosWithConsistency(river)
      error("Inconsistent data in reservoir " * getName(r) * ": different number of scenarios within the reservoir (while checking river " * getName(river) * ").")
    end
  end

  if cnt < 0
    error("Inconsistent data in reservoir " * getName(r) * ": no time series found in the rivers.")
  end

  return cnt
end

"""
Counts the number of time steps within this reservoir.
"""
countTimeSteps(r::Reservoir) = countTimeSteps(first(rivers(r)))

"""
Counts the number of time steps within this reservoir and ensures the whole data structure is consistent.
"""
function countTimeStepsWithConsistency(r::Reservoir)
  cnt = -1

  for river in rivers(r)
    if cnt < 0
      cnt = countTimeStepsWithConsistency(river)
    elseif cnt != countTimeStepsWithConsistency(river)
      error("Inconsistent data in reservoir " * getName(reservoir) * ": different number of time steps within the reservoir " *
            "(while checking river " * getName(river) * ").")
    end
  end

  if cnt < 0
    error("Inconsistent data in reservoir " * getName(reservoir) * ": no time series found in the rivers.")
  end

  return cnt
end

"""
Retrieves the periodicity of this reservoir data.
"""
getPeriod(r::Reservoir) = getPeriod(getRivers(r)[1])

"""
Retrieves the periodicity of this reservoir data and ensures the whole data structure is consistent.
"""
function getPeriodWithConsistency(r::Reservoir)
  period = :Unknown

  for river in rivers(r)
    if period == :Unknown
      period = getPeriodWithConsistency(river)
    elseif period != getPeriodWithConsistency(river)
      error("Inconsistent data in reservoir " * getName(reservoir) * ": different period within the reservoir " *
            "(while checking river " * getName(river) * ").")
    end
  end

  if period == :Unknown
    error("Inconsistent data in reservoir " * getName(reservoir) * ": no time series found in the rivers.")
  end

  return period::Period
end

"""
Retrieves the time series of this reservoir data, i.e. the time points where historical information is known.
"""
getTimePoints(r::Reservoir) = getTimePoints(getRivers(r)[1])

"""
Retrieves the time instants of this reservoir data, i.e. the time points where historical information is known.
This function also ensures the whole data structure is consistent.
"""
function getTimePointsWithConsistency(r::Reservoir)
  ts = nothing

  for river in rivers(r)
    if ts == nothing
      ts = getTimePointsWithConsistency(river)
    elseif ts != getTimePointsWithConsistency(river)
      error("Inconsistent data in reservoir " * getName(reservoir) * ": different time series within the reservoir " *
            "(while checking river " * getName(river) * ").")
    end
  end

  if ts == nothing
    error("Inconsistent data in reservoir " * getName(r) * ": no time series found in the rivers.")
  end

  return ts::Array{DateTime, 1}
end

### Getters for rivers.

# Iterate through all rivers.
rivers(r::Reservoir) = RiverIterator(r)
naturalRivers(r::Reservoir) = NaturalRiverIterator(r)
divertedRivers(r::Reservoir) = DivertedRiverIterator(r)

getNaturalRivers(r::Reservoir) = r.rivers_in
getDivertedRivers(r::Reservoir) = r.rivers_diverted
getRivers(r::Reservoir) = vcat(getNaturalRivers(r), getDivertedRivers(r))

getNaturalRiver(r::Reservoir, naturalRiver::Int) = r.rivers_in[naturalRiver]
getDivertedRiver(r::Reservoir, divertedRiver::Int) = r.rivers_diverted[divertedRiver]

hasRiver(r::Reservoir, river::AbstractString) = length(filter(nr -> getName(nr) == river, getRivers(r))) >= 1
function getRiver(r::Reservoir, river::AbstractString)
  elt = filter(nr -> getName(nr) == river, getRivers(r))
  if length(elt) < 1
    return Nullable{River}()
  else
    return Nullable(elt[1])
  end
end

# Get river name.
getDivertedRiverNames(reservoir::Reservoir) = UTF8String[r.name for r in divertedRivers(reservoir)]
getNaturalRiverNames(reservoir::Reservoir) = UTF8String[r.name for r in naturalRivers(reservoir)]

function getRiverNames(reservoir::Reservoir)
  names = getNaturalRiverNames(reservoir)
  append!(names, getDivertedRiverNames(reservoir))
  return names
end

### Agregators.

# Total inflow.

"""
Compute the total inflow for this reservoir for the given scenario,
i.e. this function agregates all tributaries to the reservoir.
The output is given along time, either as an array (`asTimeSeries=false`, default),
or as a time series object  (`asTimeSeries=true`).
"""
function totalInflow(r::Reservoir, scenario::Int; asTimeSeries::Bool=false)
  if asTimeSeries
    anyRiver = getRivers(r)[1].scenarios[1]
    return TimeArray(timestamp(anyRiver), totalInflow(r, scenario, asTimeSeries=false), colnames(anyRiver), meta(anyRiver))
  else
    return sum([values(riv.scenarios[scenario]) for riv in eachriver(r)])
  end
end

"""
Compute the total inflow for this reservoir,
i.e. this function agregates all tributaries to the reservoir.
The output is given along time, either as an array (`asTimeSeries=false`, default),
or as a time series object  (`asTimeSeries=true`).
"""
function totalInflow(r::Reservoir; kwargs...)
  if countScenarios(r) > 1
    error("This reservoir has multiple scenarios, please indicate which one to use.")
  end
  return totalInflow(r, 1; kwargs...)
end

"""
Compute the total inflow for this reservoir for the given range of scenarios,
i.e. this function agregates all tributaries to the reservoir.
The output is given along time, either as an array (`asTimeSeries=false`, default),
or as a time series object  (`asTimeSeries=true`).
"""
function totalInflow{T}(r::Reservoir, scenarios::Range{T}; kwargs...)
  kwargs = Dict(kwargs)
  if haskey(kwargs, :asTimeSeries) && kwargs[:asTimeSeries]
    return TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}[totalInflow(r, scenario; kwargs...) for scenario in scenarios]
  else
    return Float64[totalInflow(r, scenario; kwargs...) for scenario in scenarios]
  end
end

# No agregator for outputs or purposes: they have a high variability in the way to take them into account.
# Such a method would only make sense for a subset of all purposes.

### Altering the set of scenarios.

"""
Removes the given set of scenarios from the reservoir, of which a modified copy is given.
"""
function removeScenarios(reservoir::Reservoir, toRemove::IntSet)
  nScenarios = countScenarios(reservoir)
  indices = trues(nScenarios)
  for i in toRemove
    indices[i] = false
  end

  newRivers = NaturalRiver[copy(r, scenarios=r.scenarios[indices]) for r in eachnaturalriver(reservoir)]
  newDivertedRivers = DivertedRiver[copy(r, scenarios=r.scenarios[indices]) for r in eachdivertedriver(reservoir)]

  return copy(reservoir, rivers_in=newRivers, rivers_diverted=newDivertedRivers)
end

removeScenarios(r::Reservoir, toRemove::Int) = removeScenarios(r, IntSet(toRemove))
removeScenarios(r::Reservoir, toRemove::Array{Int, 1}) = removeScenarios(r, IntSet(toRemove))

"""
Limits the scenarios in the reservoir to the given time steps.
"""
function keepTimeSteps(reservoir::Reservoir, timeSteps)
  newRivers = NaturalRiver[keepTimeSteps(r, timeSteps) for r in eachnaturalriver(reservoir)]
  newDivertedRivers = DivertedRiver[keepTimeSteps(r, timeSteps) for r in eachdivertedriver(reservoir)]

  return copy(reservoir, rivers_in=newRivers, rivers_diverted=newDivertedRivers)
end

"""
Limits the scenarios in the river to the given time steps.
"""
keepTimeSteps(r::River, timeSteps) = copy(r, scenarios=[s[timeSteps] for s in r.scenarios])

### River iterators.

# Natural rivers.
immutable NaturalRiverIterator
  r::Reservoir
end

eachnaturalriver(r::Reservoir) = NaturalRiverIterator(r)
Base.start(::NaturalRiverIterator) = 1
Base.next(i::NaturalRiverIterator, s) = (i.r.rivers_in[s], s + 1)
Base.done(i::NaturalRiverIterator, s) = s > length(i.r.rivers_in)
Base.eltype(i::NaturalRiverIterator) = NaturalRiver
Base.length(i::NaturalRiverIterator) = length(i.r.rivers_in)

# Diverted rivers.
immutable DivertedRiverIterator
  r::Reservoir
end

eachdivertedriver(r::Reservoir) = DivertedRiverIterator(r)
Base.start(::DivertedRiverIterator) = 1
Base.next(i::DivertedRiverIterator, s) = (i.r.rivers_diverted[s], s + 1)
Base.done(i::DivertedRiverIterator, s) = s > length(i.r.rivers_diverted)
Base.eltype(i::DivertedRiverIterator) = DivertedRiver
Base.length(i::DivertedRiverIterator) = length(i.r.rivers_diverted)

# Finally for both at once.
# Implementation detail: first natural rivers, then diverted.
# TODO: redo with Iterators.chain()?
immutable RiverIterator
  r::Reservoir
  nIn::Int

  function RiverIterator(reservoir::Reservoir)
    return new(reservoir, length(reservoir.rivers_in))
  end
end

eachriver(r::Reservoir) = RiverIterator(r)
Base.start(::RiverIterator) = 1
Base.next(i::RiverIterator, s) = ((s <= i.nIn) ? i.r.rivers_in[s] : i.r.rivers_diverted[s - i.nIn], s + 1)
Base.done(i::RiverIterator, s) = s > i.nIn + length(i.r.rivers_diverted)
Base.eltype(i::RiverIterator) = River
Base.length(i::RiverIterator) = i.nIn + length(i.r.rivers_diverted)

### Iterate over scenarios in a reservoir.
# Provides a copy of the reservoir with a single scenario.
immutable ScenarioReservoirIterator
  r::Reservoir
  max::Int

  function ScenarioReservoirIterator(reservoir::Reservoir)
    return new(reservoir, countScenarios(reservoir))
  end
end

eachscenario(r::Reservoir) = ScenarioReservoirIterator(r)
Base.start(::ScenarioReservoirIterator) = 1
Base.next(i::ScenarioReservoirIterator, s) = (
  copy(i.r,
       rivers_in=NaturalRiver[copy(riv, scenarios=[riv.scenarios[s]]) for riv in eachnaturalriver(i.r)],
       rivers_diverted=DivertedRiver[copy(riv, scenarios=[riv.scenarios[s]]) for riv in eachdivertedriver(i.r)]),
  s + 1
)
Base.done(i::ScenarioReservoirIterator, s) = s > i.max
Base.eltype(i::ScenarioReservoirIterator) = Reservoir
Base.length(i::ScenarioReservoirIterator) = i.max


text:
module Schemata

export Schema, TableSchema, ColumnSchema, # types
       compare,                           # core function
       readschema, writeschema            # read/write schema from/to config file

include("handle_validvalues.jl")
include("types.jl")
include("compare_table_to_schema.jl")
include("readwrite.jl")

using .handle_validvalues
using .types
using .compare_table_to_schema
using .readwrite

end

ext: jl


text:
function randsdp_data(seed, m, n)
    Random.seed!(seed)
    # n = 15  # Instance size
    # m = 10  # Number of constraints
    # Objective function
    c_sqrt = rand(Float64, (n, n))
    C = c_sqrt * c_sqrt'
    # C[1, 2] *= 0.5
    # C[2, 1] = C[1, 2]
    # Generate m-dimensional feasible system
    A, b = Dict(), Dict()
    X_ = randn((n, n))
    X_ = X_ * X_'
    for i in 1:m
        A[i] = rand(Float64, (n, n))
        A[i] = A[i] * A[i]'
        b[i] = tr(A[i] * X_)
    end
    return A, b, C
end
function randsdp_eval(A,b,C,n,m,XX)
    @show minus_rank = length([eig for eig in eigen(XX).values if eig < -1e-10])

    @show rank = length([eig for eig in eigen(XX).values if eig > 1e-10])

    @show tr(C * XX)
    for i in 1:m
        @show tr(A[i] * XX)-b[i]
    end

    nothing
end


text:
@testset "Q" begin
    @testset "q2m" begin
        exp = [0.607843137254902 0.27450980392156854 0.7450980392156862;
               0.6666666666666666 0.33333333333333326 -0.6666666666666666;
               -0.43137254901960775 0.9019607843137255 0.019607843137254832]
        q = [0.5, 0.4, 0.3, 0.1]
        act = q2m(q)
        @testset for i in eachindex(act, exp)
            @test act[i] ≈ exp[i]
        end
    end
    @testset "qdq2av" begin
        angle = deg2rad.([-20.0, 50.0, -60.0])
        m = eul2m(angle[3], angle[2], angle[1], 3, 1, 3)
        q = m2q(m)
        exp = [1.0, 2.0, 3.0]
        qav = [0.0, 1.0, 2.0, 3.0]
        dq = -0.5 * qxq(q, qav)
        act = qdq2av(q, dq)
        @testset for i in eachindex(act, exp)
            @test act[i] ≈ exp[i]
        end
    end
    @testset "qxq" begin
        q_id = [1.0, 0.0, 0.0, 0.0]
        nq_id = [-1.0, 0.0, 0.0, 0.0]
        q_i = [0.0, 1.0, 0.0, 0.0]
        q_j = [0.0, 0.0, 1.0, 0.0]
        q_k = [0.0, 0.0, 0.0, 1.0]
        q_ij = qxq(q_i, q_j)
        @testset for i in eachindex(q_ij, q_k)
            @test q_ij[i] ≈ q_k[i]
        end
        q_jk = qxq(q_j, q_k)
        @testset for i in eachindex(q_jk, q_i)
            @test q_jk[i] ≈ q_i[i]
        end
        q_ki = qxq(q_k, q_i)
        @testset for i in eachindex(q_ki, q_j)
            @test q_ki[i] ≈ q_j[i]
        end
        q_ii = qxq(q_i, q_i)
        @testset for i in eachindex(q_ii, nq_id)
            @test q_ii[i] ≈ nq_id[i]
        end
        q_jj = qxq(q_j, q_j)
        @testset for i in eachindex(q_jj, nq_id)
            @test q_jj[i] ≈ nq_id[i]
        end
        q_kk = qxq(q_k, q_k)
        @testset for i in eachindex(q_kk, nq_id)
            @test q_kk[i] ≈ nq_id[i]
        end
        q_idi = qxq(q_id, q_i)
        @testset for i in eachindex(q_idi, q_i)
            @test q_idi[i] ≈ q_i[i]
        end
        q_iid = qxq(q_i, q_id)
        @testset for i in eachindex(q_iid, q_i)
            @test q_iid[i] ≈ q_i[i]
        end
    end
end


text:
function OffsetModel(arg0::DateComponents, arg1::jint)
    return OffsetModel((DateComponents, jint), arg0, arg1)
end

function OffsetModel(arg0::DateComponents, arg1::jint, arg2::jdouble, arg3::jdouble)
    return OffsetModel((DateComponents, jint, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function get_mjd_ref(obj::OffsetModel)
    return jcall(obj, "getMJDRef", jint, ())
end

function get_offset(obj::OffsetModel)
    return jcall(obj, "getOffset", jdouble, ())
end

function get_slope(obj::OffsetModel)
    return jcall(obj, "getSlope", jdouble, ())
end

function get_start(obj::OffsetModel)
    return jcall(obj, "getStart", DateComponents, ())
end


text:
function rank(obj::RankingAlgorithm, arg0::Vector{jdouble})
    return jcall(obj, "rank", Vector{jdouble}, (Vector{jdouble},), arg0)
end


text:
module B
function myfunc()
    println("hello")
end
end


text:
using GNSSTools
using PyPlot
pygui(true)

fig = figure(figsize=(6.5,2))
ax1 = fig.add_subplot(1,2,1)
ax1.step(Array(1:10) .- 0.5, nh10, "k-")
ax1.set_yticks([0, 1])
ax1.set_ylim([-0.1, 1.1])
ax1.set_xticks([0, 5, 10])
title("(a)")
xlabel("Time (ms)")
ax2 = fig.add_subplot(1,2,2)
ax2.step(Array(1:20) .- 0.5, nh20, "k-")
ax2.set_yticks([0, 1])
ax2.set_ylim([-0.1, 1.1])
ax2.set_xticks([0, 10, 20])
title("(b)")
xlabel("Time (ms)")
fig.subplots_adjust(hspace=0.2, wspace=0.3, top=0.88, left=0.08, right=0.93, bottom=0.22)
savefig("figures/ch2_neuman_sequence.pdf", dpi=300)


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

import .Base: unsafe_convert, lock, trylock, unlock, islocked, wait, notify, AbstractLock

# Important Note: these low-level primitives defined here
#   are typically not for general usage

##########################################
# Atomic Locks
##########################################

# Test-and-test-and-set spin locks are quickest up to about 30ish
# contending threads. If you have more contention than that, perhaps
# a lock is the wrong way to synchronize.
"""
    SpinLock()

Create a non-reentrant lock.
Recursive use will result in a deadlock.
Each [`lock`](@ref) must be matched with an [`unlock`](@ref).

Test-and-test-and-set spin locks are quickest up to about 30ish
contending threads. If you have more contention than that, perhaps
a lock is the wrong way to synchronize.

See also [`Mutex`](@ref) for a more efficient version on one core or if the
lock may be held for a considerable length of time.
"""
struct SpinLock <: AbstractLock
    handle::Atomic{Int}
    SpinLock() = new(Atomic{Int}(0))
end

function lock(l::SpinLock)
    while true
        if l.handle[] == 0
            p = atomic_xchg!(l.handle, 1)
            if p == 0
                return
            end
        end
        ccall(:jl_cpu_pause, Cvoid, ())
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Cvoid, ())
    end
end

function trylock(l::SpinLock)
    if l.handle[] == 0
        return atomic_xchg!(l.handle, 1) == 0
    end
    return false
end

function unlock(l::SpinLock)
    l.handle[] = 0
    ccall(:jl_cpu_wake, Cvoid, ())
    return
end

function islocked(l::SpinLock)
    return l.handle[] != 0
end

##########################################
# System Mutexes
##########################################

# These are mutexes from libuv.
const UV_MUTEX_SIZE = ccall(:jl_sizeof_uv_mutex, Cint, ())

"""
    Mutex()

These are standard system mutexes for locking critical sections of logic.

On Windows, this is a critical section object,
on pthreads, this is a `pthread_mutex_t`.

See also [`SpinLock`](@ref) for a lighter-weight lock.
"""
mutable struct Mutex <: AbstractLock
    ownertid::Int16
    handle::Ptr{Cvoid}
    function Mutex()
        m = new(zero(Int16), Libc.malloc(UV_MUTEX_SIZE))
        ccall(:uv_mutex_init, Cvoid, (Ptr{Cvoid},), m.handle)
        finalizer(mutex_destroy, m)
        return m
    end
end

unsafe_convert(::Type{Ptr{Cvoid}}, m::Mutex) = m.handle

function mutex_destroy(x::Mutex)
    h = x.handle
    if h != C_NULL
        x.handle = C_NULL
        ccall(:uv_mutex_destroy, Cvoid, (Ptr{Cvoid},), h)
        Libc.free(h)
        nothing
    end
end

function lock(m::Mutex)
    m.ownertid == threadid() && concurrency_violation() # deadlock
    # Temporary solution before we have gc transition support in codegen.
    # This could mess up gc state when we add codegen support.
    gc_state = ccall(:jl_gc_safe_enter, Int8, ())
    ccall(:uv_mutex_lock, Cvoid, (Ptr{Cvoid},), m)
    ccall(:jl_gc_safe_leave, Cvoid, (Int8,), gc_state)
    m.ownertid = threadid()
    return
end

function trylock(m::Mutex)
    m.ownertid == threadid() && concurrency_violation() # deadlock
    r = ccall(:uv_mutex_trylock, Cint, (Ptr{Cvoid},), m)
    if r == 0
        m.ownertid = threadid()
    end
    return r == 0
end

function unlock(m::Mutex)
    m.ownertid == threadid() || concurrency_violation()
    m.ownertid = 0
    ccall(:uv_mutex_unlock, Cvoid, (Ptr{Cvoid},), m)
    return
end

function islocked(m::Mutex)
    return m.ownertid != 0
end

max_stars_repo_name: terasakisatoshi/julia

max_issues_repo_name: terasakisatoshi/julia

max_forks_repo_name: terasakisatoshi/julia


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

mutable struct BgpPeerStatus <: SwaggerModel
    localAddress::Any # spec type: Union{ Nothing, String } # spec name: localAddress
    neighbor::Any # spec type: Union{ Nothing, String } # spec name: neighbor
    asn::Any # spec type: Union{ Nothing, Int32 } # spec name: asn
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    connectedDuration::Any # spec type: Union{ Nothing, String } # spec name: connectedDuration
    routesReceived::Any # spec type: Union{ Nothing, Int64 } # spec name: routesReceived
    messagesSent::Any # spec type: Union{ Nothing, Int64 } # spec name: messagesSent
    messagesReceived::Any # spec type: Union{ Nothing, Int64 } # spec name: messagesReceived

    function BgpPeerStatus(;localAddress=nothing, neighbor=nothing, asn=nothing, state=nothing, connectedDuration=nothing, routesReceived=nothing, messagesSent=nothing, messagesReceived=nothing)
        o = new()
        validate_property(BgpPeerStatus, Symbol("localAddress"), localAddress)
        setfield!(o, Symbol("localAddress"), localAddress)
        validate_property(BgpPeerStatus, Symbol("neighbor"), neighbor)
        setfield!(o, Symbol("neighbor"), neighbor)
        validate_property(BgpPeerStatus, Symbol("asn"), asn)
        setfield!(o, Symbol("asn"), asn)
        validate_property(BgpPeerStatus, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(BgpPeerStatus, Symbol("connectedDuration"), connectedDuration)
        setfield!(o, Symbol("connectedDuration"), connectedDuration)
        validate_property(BgpPeerStatus, Symbol("routesReceived"), routesReceived)
        setfield!(o, Symbol("routesReceived"), routesReceived)
        validate_property(BgpPeerStatus, Symbol("messagesSent"), messagesSent)
        setfield!(o, Symbol("messagesSent"), messagesSent)
        validate_property(BgpPeerStatus, Symbol("messagesReceived"), messagesReceived)
        setfield!(o, Symbol("messagesReceived"), messagesReceived)
        o
    end
end # type BgpPeerStatus

const _property_map_BgpPeerStatus = Dict{Symbol,Symbol}(Symbol("localAddress")=>Symbol("localAddress"), Symbol("neighbor")=>Symbol("neighbor"), Symbol("asn")=>Symbol("asn"), Symbol("state")=>Symbol("state"), Symbol("connectedDuration")=>Symbol("connectedDuration"), Symbol("routesReceived")=>Symbol("routesReceived"), Symbol("messagesSent")=>Symbol("messagesSent"), Symbol("messagesReceived")=>Symbol("messagesReceived"))
const _property_types_BgpPeerStatus = Dict{Symbol,String}(Symbol("localAddress")=>"String", Symbol("neighbor")=>"String", Symbol("asn")=>"Int32", Symbol("state")=>"String", Symbol("connectedDuration")=>"String", Symbol("routesReceived")=>"Int64", Symbol("messagesSent")=>"Int64", Symbol("messagesReceived")=>"Int64")
Base.propertynames(::Type{ BgpPeerStatus }) = collect(keys(_property_map_BgpPeerStatus))
Swagger.property_type(::Type{ BgpPeerStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_BgpPeerStatus[name]))}
Swagger.field_name(::Type{ BgpPeerStatus }, property_name::Symbol) =  _property_map_BgpPeerStatus[property_name]

const _allowed_BgpPeerStatus_state = ["Unknown", "Stopped", "Idle", "Connecting", "Connected"]

function check_required(o::BgpPeerStatus)
    true
end

function validate_property(::Type{ BgpPeerStatus }, name::Symbol, val)
    if name === Symbol("state")
        Swagger.validate_param(name, "BgpPeerStatus", :enum, val, _allowed_BgpPeerStatus_state)
    end
end


text:
# SAM Reader
# =========

mutable struct Reader <: BioCore.IO.AbstractReader
    state::BioCore.Ragel.State
    header::Header

    function Reader(input::BufferedStreams.BufferedInputStream)
        reader = new(BioCore.Ragel.State(sam_header_machine.start_state, input), Header())
        readheader!(reader)
        reader.state.cs = sam_body_machine.start_state
        return reader
    end
end

"""
    SAM.Reader(input::IO)

Create a data reader of the SAM file format.

# Arguments
* `input`: data source
"""
function Reader(input::IO)
    return Reader(BufferedStreams.BufferedInputStream(input))
end

function BioCore.IO.stream(reader::Reader)
    return reader.state.stream
end

"""
    header(reader::Reader)::Header

Get the header of `reader`.
"""
function header(reader::Reader)::Header
    return reader.header
end

function BioCore.header(reader::Reader)
    return header(reader)
end

function Base.eltype(::Type{Reader})
    return Record
end

# file   = header . body
# header = metainfo*
# body   = record*
isinteractive() && info("compiling SAM")
const sam_metainfo_machine, sam_record_machine, sam_header_machine, sam_body_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt
    any = Automa.RegExp.any

    metainfo = let
        tag = re"[A-Z][A-Z]" \ cat("CO")
        tag.actions[:enter] = [:mark1]
        tag.actions[:exit]  = [:metainfo_tag]

        dict = let
            key = re"[A-Za-z][A-Za-z0-9]"
            key.actions[:enter] = [:mark2]
            key.actions[:exit]  = [:metainfo_dict_key]
            val = re"[ -~]+"
            val.actions[:enter] = [:mark2]
            val.actions[:exit]  = [:metainfo_dict_val]
            keyval = cat(key, ':', val)

            cat(keyval, rep(cat('\t', keyval)))
        end
        dict.actions[:enter] = [:mark1]
        dict.actions[:exit]  = [:metainfo_val]

        co = cat("CO")
        co.actions[:enter] = [:mark1]
        co.actions[:exit]  = [:metainfo_tag]

        comment = re"[^\r\n]*"
        comment.actions[:enter] = [:mark1]
        comment.actions[:exit]  = [:metainfo_val]

        cat('@', alt(cat(tag, '\t', dict), cat(co, '\t', comment)))
    end
    metainfo.actions[:enter] = [:anchor]
    metainfo.actions[:exit]  = [:metainfo]

    record = let
        qname = re"[!-?A-~]+"
        qname.actions[:enter] = [:mark]
        qname.actions[:exit]  = [:record_qname]

        flag = re"[0-9]+"
        flag.actions[:enter] = [:mark]
        flag.actions[:exit]  = [:record_flag]

        rname = re"\*|[!-()+-<>-~][!-~]*"
        rname.actions[:enter] = [:mark]
        rname.actions[:exit]  = [:record_rname]

        pos = re"[0-9]+"
        pos.actions[:enter] = [:mark]
        pos.actions[:exit]  = [:record_pos]

        mapq = re"[0-9]+"
        mapq.actions[:enter] = [:mark]
        mapq.actions[:exit]  = [:record_mapq]

        cigar = re"\*|([0-9]+[MIDNSHPX=])+"
        cigar.actions[:enter] = [:mark]
        cigar.actions[:exit]  = [:record_cigar]

        rnext = re"\*|=|[!-()+-<>-~][!-~]*"
        rnext.actions[:enter] = [:mark]
        rnext.actions[:exit]  = [:record_rnext]

        pnext = re"[0-9]+"
        pnext.actions[:enter] = [:mark]
        pnext.actions[:exit]  = [:record_pnext]

        tlen = re"[-+]?[0-9]+"
        tlen.actions[:enter] = [:mark]
        tlen.actions[:exit]  = [:record_tlen]

        seq = re"\*|[A-Za-z=.]+"
        seq.actions[:enter] = [:mark]
        seq.actions[:exit]  = [:record_seq]

        qual = re"[!-~]+"
        qual.actions[:enter] = [:mark]
        qual.actions[:exit]  = [:record_qual]

        field = let
            tag = re"[A-Za-z][A-Za-z0-9]"
            val = alt(
                re"A:[!-~]",
                re"i:[-+]?[0-9]+",
                re"f:[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?",
                re"Z:[ !-~]*",
                re"H:([0-9A-F][0-9A-F])*",
                re"B:[cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+")

            cat(tag, ':', val)
        end
        field.actions[:enter] = [:mark]
        field.actions[:exit]  = [:record_field]

        cat(
            qname, '\t',
            flag,  '\t',
            rname, '\t',
            pos,   '\t',
            mapq,  '\t',
            cigar, '\t',
            rnext, '\t',
            pnext, '\t',
            tlen,  '\t',
            seq,   '\t',
            qual,
            rep(cat('\t', field)))
    end
    record.actions[:enter] = [:anchor]
    record.actions[:exit]  = [:record]

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(re"\r?", lf)
    end

    header′ = rep(cat(metainfo, newline))
    header′.actions[:exit] = [:header]
    header = cat(header′, opt(any() \ cat('@')))  # look ahead

    body = rep(cat(record, newline))

    return map(Automa.compile, (metainfo, record, header, body))
end)()

const sam_metainfo_actions = Dict(
    :metainfo_tag => :(record.tag = (mark1:p-1) .- offset),
    :metainfo_val => :(record.val = (mark1:p-1) .- offset),
    :metainfo_dict_key => :(push!(record.dictkey, (mark2:p-1) .- offset)),
    :metainfo_dict_val => :(push!(record.dictval, (mark2:p-1) .- offset)),
    :metainfo => quote
        BioCore.ReaderHelper.resize_and_copy!(record.data, data, offset+1:p-1)
        record.filled = (offset+1:p-1) .- offset
    end,
    :anchor => :(),
    :mark1  => :(mark1 = p),
    :mark2  => :(mark2 = p))
eval(
    BioCore.ReaderHelper.generate_index_function(
        MetaInfo,
        sam_metainfo_machine,
        :(mark1 = mark2 = offset = 0),
        sam_metainfo_actions))
eval(
    BioCore.ReaderHelper.generate_readheader_function(
        Reader,
        MetaInfo,
        sam_header_machine,
        :(mark1 = mark2 = offset = 0),
        merge(sam_metainfo_actions, Dict(
            :metainfo => quote
                BioCore.ReaderHelper.resize_and_copy!(record.data, data, BioCore.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) .- offset
                @assert isfilled(record)
                push!(reader.header.metainfo, record)
                BioCore.ReaderHelper.ensure_margin!(stream)
                record = MetaInfo()
            end,
            :header => :(finish_header = true; @escape),
            :countline => :(linenum += 1),
            :anchor => :(BioCore.ReaderHelper.anchor!(stream, p); offset = p - 1))),
        quote
            if !eof(stream)
                stream.position -= 1  # cancel look-ahead
            end
        end))

const sam_record_actions = Dict(
    :record_qname => :(record.qname = (mark:p-1) .- offset),
    :record_flag  => :(record.flag  = (mark:p-1) .- offset),
    :record_rname => :(record.rname = (mark:p-1) .- offset),
    :record_pos   => :(record.pos   = (mark:p-1) .- offset),
    :record_mapq  => :(record.mapq  = (mark:p-1) .- offset),
    :record_cigar => :(record.cigar = (mark:p-1) .- offset),
    :record_rnext => :(record.rnext = (mark:p-1) .- offset),
    :record_pnext => :(record.pnext = (mark:p-1) .- offset),
    :record_tlen  => :(record.tlen  = (mark:p-1) .- offset),
    :record_seq   => :(record.seq   = (mark:p-1) .- offset),
    :record_qual  => :(record.qual  = (mark:p-1) .- offset),
    :record_field => :(push!(record.fields, (mark:p-1) .- offset)),
    :record       => quote
        BioCore.ReaderHelper.resize_and_copy!(record.data, data, 1:p-1)
        record.filled = (offset+1:p-1) .- offset
    end,
    :anchor       => :(),
    :mark         => :(mark = p))
eval(
    BioCore.ReaderHelper.generate_index_function(
        Record,
        sam_record_machine,
        :(mark = offset = 0),
        sam_record_actions))
eval(
    BioCore.ReaderHelper.generate_read_function(
        Reader,
        sam_body_machine,
        :(mark = offset = 0),
        merge(sam_record_actions, Dict(
            :record    => quote
                BioCore.ReaderHelper.resize_and_copy!(record.data, data, BioCore.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) .- offset
                found_record = true
                @escape
            end,
            :countline => :(linenum += 1),
            :anchor    => :(BioCore.ReaderHelper.anchor!(stream, p); offset = p - 1)))))


text:
"""
    poisson_eigenvalues(N, L, dim, ::Periodic)

Return the eigenvalues satisfying the discrete form of Poisson's equation
with periodic boundary conditions along the dimension `dim` with `N` grid
points and domain extent `L`.
"""
function poisson_eigenvalues(N, L, dim, ::Periodic)
    inds = reshape(1:N, reshaped_size(N, dim)...)
    return @. (2sin((inds - 1) * π / N) / (L / N))^2
end

"""
    poisson_eigenvalues(N, L, dim, ::Bounded)

Return the eigenvalues satisfying the discrete form of Poisson's equation
with staggered Neumann boundary conditions along the dimension `dim` with
`N` grid points and domain extent `L`.
"""
function poisson_eigenvalues(N, L, dim, ::Bounded)
    inds = reshape(1:N, reshaped_size(N, dim)...)
    return @. (2sin((inds - 1) * π / 2N) / (L / N))^2
end

"""
    poisson_eigenvalues(N, L, dim, ::Flat)

No need to do any transforms along `Flat` dimensions so just return zeros
since they will be added to the eigenvalues for other dimensions.
"""
poisson_eigenvalues(N, L, dim, ::Flat) = zeros(1, 1, 1)


text: AbstractLDA = Union{LDA, gpuLDA}


text:
using Flux.Tracker, Statistics, DelimitedFiles
using Flux.Tracker: Params, gradient, update!
using DelimitedFiles, Statistics
using Flux: gpu

# This replicates the housing data example from the Knet.jl readme. Although we
# could have reused more of Flux (see the mnist example), the library's
# abstractions are very lightweight and don't force you into any particular
# strategy.

cd(@__DIR__)

isfile("housing.data") ||
  download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
           "housing.data")

rawdata = readdlm("housing.data")'

# The last feature is our target -- the price of the house.
split_ratio = 0.1 # For the train test split

x = rawdata[1:13,:] |> gpu
y = rawdata[14:14,:] |> gpu

# Normalise the data
x = (x .- mean(x, dims = 2)) ./ std(x, dims = 2)

# Split into train and test sets
split_index = floor(Int,size(x,2)*split_ratio)
x_train = x[:,1:split_index]
y_train = y[:,1:split_index]
x_test = x[:,split_index+1:size(x,2)]
y_test = y[:,split_index+1:size(x,2)]

# The model
W = param(randn(1,13)/10) |> gpu
b = param([0.]) |> gpu

predict(x) = W*x .+ b
meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)
loss(x, y) = meansquarederror(predict(x), y)

η = 0.1
θ = Params([W, b])

for i = 1:10
  g = gradient(() -> loss(x_train, y_train), θ)
  for x in θ
    update!(x, -g[x]*η)
  end
  @show loss(x_train, y_train)
end

# Predict the RMSE on the test set
err = meansquarederror(predict(x_test),y_test)
println(err)

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using Test

using MagicUnderscores
using MagicUnderscores: lower_underscores

@testset "MagicUnderscores" begin
    @test @_(map(_^2, [1,2])) == [1,4]
    @test @_(map(sqrt(abs(_)) + 1, [-1,2,3])) == sqrt.(abs.([-1,2,3])) .+ 1
    @test @_(map(map(x->x^2, _), [[1,2]])) == [[1,4]]

    @test @_(filter(_>2, [1,2])) == []

    @test @_([1,2] |> length) == 2
    @test @_([1,2,3,4] |> filter(_>2, _)) == [3,4]
    @test @_([1,2,3,4] |> filter(_>2, _) |> length) == 2

    @test @_(map(_.a == 0, [(a=1, b=2), (a=0, b=2)])) == [false, true]
    @test @_(findall(2*_ in (4,6), [1,2,3])) == [2,3]
end


text:
function update_edge_trees!(edge_trees, tree, σ, α, ε, fs, rates_all, branch_lengths_all ; do_tips = true, max_try = 1_000,
    max_node_number = 50, keep_if_any = true, rep = 1)
    n_edges = length(edge_trees)
    branch_lengths = branch_lengths_all[2:end]
    left = n_left(tree)

    for edge_id in 1:n_edges#:-1:1
        f = fs[edge_id]
        change_edge_trees_rrMH!(edge_trees, tree, edge_id, σ, α, ε, f, rates_all, branch_lengths,left,
            max_try = max_try, max_node_number = max_node_number, rep = rep);
    end

    update_rates!(tree,rates_all)
end

function update_edge_trees!(edge_trees, tree, σ, α, ε, fs; max_try = 1_000, max_node_number = 50, keep_if_any = true, do_tips = true)
    rates_all = extract_rates(tree)
    branch_lengths_all = extract_branch_lengths(tree)

    update_edge_trees!(edge_trees, tree, σ, α, ε, fs, rates_all, branch_lengths_all, max_try = max_try,
        max_node_number = max_node_number, keep_if_any = keep_if_any, do_tips = do_tips)
end

function change_edge_trees_rrMH!(edge_trees, tree, edge_id, σ, α, ε, f, rates, branch_lengths, lefts;
    max_try = 1_000, max_node_number = 50, rep = 1)

    log_α = log(α)
    lambda_law = LogNormal(log_α,σ)                                                        # relative new rate density

    trial = 0
    parent_rate = get_parent_rate(tree, edge_id, edge_trees, rates)
    #pr2 = rates[1]
    #if edge_trees[edge_id].parent_edge > 0
    #    pr2 = edge_trees[edge_trees[edge_id].parent_edge].stem_rate[1] *edge_trees[edge_trees[edge_id].parent_edge].tip_rate             # the parent rate in the augmented tree
    #end
    #println("$pr2, $parent_rate")
    depth = edge_trees[edge_id].depth
    sd = edge_trees[edge_id].stem_depth

    parent_edge = edge_trees[edge_id].parent_edge

    if lefts[edge_id+1] == 0                                                                   # meaning it's a tip branch

        function new_rates_tip(l)                                                                   # new rate density
            l * rand(lambda_law, 1)[1]
        end

        int_tree = Tree()                                                                       # initialize the tree
        n_alive = edge_trees[edge_id].tip_number                                                    # tip number in the current MCMC state

        if n_alive == 0                                                                         # compute the former likelihood of having exactly one tip in the reconstructed tree
            former_lik = -Inf
        elseif n_alive == 1
            former_lik = log(f)
        else
            former_lik = log(n_alive) + log(f) + (n_alive - 1) * log(1-f)
        end

        while trial < rep

            trial += 1
            ext = 0

            while ext < max_try                                                                 # conditionned to the non-extinction (if you don't want it put max_try to 1)
                ext += 1

                u = log(rand())
                λ = new_rates_tip(parent_rate)                                                      # draw the new stem_rate
                int_tree, unsampled, nalive, tip_rates =
                    sim_ClaDS2_time_unsampled_rates(branch_lengths[edge_id]*λ, σ, α, ε,
                        1.,u,-1*former_lik,
                        max_node_number = max_node_number,
                        sampling_proba = f, return_if_extinct = true)
                #println("$edge_id ; $trial , $ext , $nalive, $unsampled, $(unsampled && nalive > 0 )")

                if (unsampled && nalive > 0 )                                                       # one tip sampled exactly
                    rates[edge_id+1] = λ
                    edge_trees[edge_id] =
                        EdgeTreeRates2(int_tree, λ, 1., 1,
                            tip_rates, nalive, depth,sd, parent_edge)

                    #former_lik = log(nalive) + log(f) + (nalive - 1) * log(1-f)               # new likelihood
                    break
                elseif nalive > 0
                    break
                end
                #println(ext)
            end
        end

        return edge_trees[edge_id]

    else                                                                                # meaning it's an internal branch
        daughter_rates = [rates[edge_id+2]; rates[edge_id+2+lefts[edge_id+1]]]                  # the daugther rates in that case
        stem_rate = rates[edge_id+1]
        function new_rates_int(l)                                                                   # new rate density
            l * rand(lambda_law, 1)[1]
        end

        function nu(former_rate, former_ratio)                                                  # speciating at t_i and giving birth to the new  daugther rates
            if isnan(former_rate)
                return 0.
            end
            fr = former_rate*former_ratio
            log_fr = log(fr)
            return fr * exp((-(log_fr + log_α - log(daughter_rates[1]))^2 -
                (log_fr + log_α - log(daughter_rates[2]))^2)/(2 * σ^2))
        end

        sum_rates = 0.
        for r in edge_trees[edge_id].tip_rates
            sum_rates += nu(r, stem_rate)
        end
        former_lik_rate = log(sum_rates)                                                        # Σ λνν
        former_tip_number = edge_trees[edge_id].tip_number
        former_lik = former_lik_rate

        if former_tip_number == 0
            former_lik += 0
        else
            former_lik += former_tip_number * log(1-f)                                          # no sampled tip
        end

        it = 0
        λ = 1.

        while it < rep
            it += 1
            trial = 0

            while trial < max_try
                nalive = 0
                int_tree = Tree()
                int_tip_rates = Array{Float64,1}(undef,0)
                which_alive = Array{Bool,1}(undef,0)

                while (nalive == 0) && (trial < max_try)

                    trial += 1
                    λ = new_rates_int(parent_rate)

                    int_tree, int_tip_rates, which_alive =
                        sim_ClaDS2_time_rates(branch_lengths[edge_id]*λ, σ, α, ε, 1.,
                            prune_extinct = false, max_node_number = max_node_number,
                            return_if_max = true, return_na = true)

                    if int_tree.branch_length<=0
                        nalive = Inf
                    else
                        nalive = length(int_tip_rates)
                    end
                    #=if nalive>1
                        println("$edge_id ; $trial $nalive $ε $(branch_lengths[edge_id]*λ) $σ $α")
                        println(int_tip_rates)
                        plot_ClaDS(int_tree)
                    end=#
                end

                if nalive == Inf
                    break
                end

               if (0 < nalive)

                   w = Array{Float64,1}(undef,0)
                   Ws = 0.
                   for r in 1:length(int_tip_rates)
                       push!(w,nu(int_tip_rates[r],λ))
                       Ws += w[end]
                   end
                   new_lik = log(Ws)                                                        # Σ λνν

                   new_tip = sample(1:nalive, Weights(w))

                   new_tip_id = NaN
                   new_parent_rate = 0.
                   former_tip_number = 0

                   k = 0
                   for i in 1:length(which_alive)
                       if which_alive[i]
                           k += 1
                           if k == new_tip
                               new_parent_rate = int_tip_rates[k]
                               new_tip_id = i
                               #which_alive[i] = false
                           end
                       end
                   end

                   u = log(rand())

                   keep = u < (new_lik-former_lik)

                   ntips = 0
                   n_alive = 0
                   current_lik = new_lik - former_lik

                   k = 0
                   i = 0
                   tip_trees = []
                   graft_tips = []
                   add_to_tip_id = 0
                   new_tip_number = 0

                   while keep & (i < length(which_alive))
                       i += 1

                       if which_alive[i]
                           k += 1
                           if k != new_tip
                               tip_tree, current_lik =
                                    sim_ClaDS2_time(depth*λ, σ, α, ε,
                                        int_tip_rates[k],u,
                                        current_lik,log(1-f),
                                        max_node_number = max_node_number)

                                   if tip_tree.branch_length < 0                            # more than max_node_number nodes, or sampled
                                       keep = false
                                   else
                                       push!(tip_trees, tip_tree)
                                       push!(graft_tips, i)
                                       nt = n_extant_tips(tip_tree)
                                       new_tip_number += nt

                                       if i < new_tip_id
                                           add_to_tip_id += n_tip(tip_tree) - 1
                                       end
                                   end
                              end
                         end

                   end

                   if keep
                       rates[edge_id+1] = λ
                       to_graft = graft_to_tips(int_tree,
                            [Tree(Array{Tree,1}(undef,0),0., [0.], 0, true); tip_trees],    # the tree
                            [new_tip_id; graft_tips])
                       edge_trees[edge_id] =
                            EdgeTreeRates2(to_graft,
                                λ,                                                              # stem_rate
                                new_parent_rate,                                                # tip_rate
                                new_tip_id + add_to_tip_id,                                     # tip_id
                                int_tip_rates,                                                  # tip_rates
                                new_tip_number,
                                depth, sd, parent_edge)

                        if it < rep
                            former_lik = log(Ws)

                            if new_tip_number == 0
                                former_lik += 0
                            else
                                former_lik += new_tip_number * log(1-f)                                          # no sampled tip
                            end
                        end
                   end

                   trial = max_try
               end
           end
        end
        return edge_trees[edge_id]
    end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


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

#Macros
macro OK(func)
    :($(esc(func)) == SQLITE_OK)
end

macro CHECK(db,ex)
    esc(quote
        if !(@OK $ex)
            sqliteerror($db)
        end
        SQLITE_OK
    end)
end

const SQLNullPtrError = SQLiteException("Cannot operate on null pointer")
macro NULLCHECK(ptr)
    esc(quote
        if $ptr == C_NULL
            throw(SQLNullPtrError)
        end
    end)
end

#Return codes
const SQLITE_OK =           0   # /* Successful result */
const SQLITE_ERROR =        1   # /* SQL error or missing database */
const SQLITE_INTERNAL =     2   # /* Internal logic error in SQLite */
const SQLITE_PERM =         3   # /* Access permission denied */
const SQLITE_ABORT =        4   # /* Callback routine requested an abort */
const SQLITE_BUSY =         5   # /* The database file is locked */
const SQLITE_LOCKED =       6   # /* A table in the database is locked */
const SQLITE_NOMEM =        7   # /* A malloc() failed */
const SQLITE_READONLY =     8   # /* Attempt to write a readonly database */
const SQLITE_INTERRUPT =    9   # /* Operation terminated by sqlite3_interrupt()*/
const SQLITE_IOERR =       10   # /* Some kind of disk I/O error occurred */
const SQLITE_CORRUPT =     11   # /* The database disk image is malformed */
const SQLITE_NOTFOUND =    12   # /* Unknown opcode in sqlite3_file_control() */
const SQLITE_FULL =        13   # /* Insertion failed because database is full */
const SQLITE_CANTOPEN =    14   # /* Unable to open the database file */
const SQLITE_PROTOCOL =    15   # /* Database lock protocol error */
const SQLITE_EMPTY =       16   # /* Database is empty */
const SQLITE_SCHEMA =      17   # /* The database schema changed */
const SQLITE_TOOBIG =      18   # /* String or BLOB exceeds size limit */
const SQLITE_CONSTRAINT =  19   # /* Abort due to constraint violation */
const SQLITE_MISMATCH =    20   # /* Data type mismatch */
const SQLITE_MISUSE =      21   # /* Library used incorrectly */
const SQLITE_NOLFS =       22   # /* Uses OS features not supported on host */
const SQLITE_AUTH =        23   # /* Authorization denied */
const SQLITE_FORMAT =      24   # /* Auxiliary database format error */
const SQLITE_RANGE =       25   # /* 2nd parameter to sqlite3_bind out of range */
const SQLITE_NOTADB =      26   # /* File opened that is not a database file */
const SQLITE_NOTICE =      27   # /* Notifications from sqlite3_log() */
const SQLITE_WARNING =     28   # /* Warnings from sqlite3_log() */
const SQLITE_ROW =         100  # /* sqlite3_step() has another row ready */
const SQLITE_DONE =        101  # /* sqlite3_step() has finished executing */
#Extended Return codes
const SQLITE_IOERR_READ =              (SQLITE_IOERR | (1<<8))
const SQLITE_IOERR_SHORT_READ =        (SQLITE_IOERR | (2<<8))
const SQLITE_IOERR_WRITE =             (SQLITE_IOERR | (3<<8))
const SQLITE_IOERR_FSYNC =             (SQLITE_IOERR | (4<<8))
const SQLITE_IOERR_DIR_FSYNC =         (SQLITE_IOERR | (5<<8))
const SQLITE_IOERR_TRUNCATE =          (SQLITE_IOERR | (6<<8))
const SQLITE_IOERR_FSTAT =             (SQLITE_IOERR | (7<<8))
const SQLITE_IOERR_UNLOCK =            (SQLITE_IOERR | (8<<8))
const SQLITE_IOERR_RDLOCK =            (SQLITE_IOERR | (9<<8))
const SQLITE_IOERR_DELETE =            (SQLITE_IOERR | (10<<8))
const SQLITE_IOERR_BLOCKED =           (SQLITE_IOERR | (11<<8))
const SQLITE_IOERR_NOMEM =             (SQLITE_IOERR | (12<<8))
const SQLITE_IOERR_ACCESS =            (SQLITE_IOERR | (13<<8))
const SQLITE_IOERR_CHECKRESERVEDLOCK = (SQLITE_IOERR | (14<<8))
const SQLITE_IOERR_LOCK =              (SQLITE_IOERR | (15<<8))
const SQLITE_IOERR_CLOSE =             (SQLITE_IOERR | (16<<8))
const SQLITE_IOERR_DIR_CLOSE =         (SQLITE_IOERR | (17<<8))
const SQLITE_IOERR_SHMOPEN =           (SQLITE_IOERR | (18<<8))
const SQLITE_IOERR_SHMSIZE =           (SQLITE_IOERR | (19<<8))
const SQLITE_IOERR_SHMLOCK =           (SQLITE_IOERR | (20<<8))
const SQLITE_IOERR_SHMMAP =            (SQLITE_IOERR | (21<<8))
const SQLITE_IOERR_SEEK =              (SQLITE_IOERR | (22<<8))
const SQLITE_IOERR_DELETE_NOENT =      (SQLITE_IOERR | (23<<8))
const SQLITE_IOERR_MMAP =              (SQLITE_IOERR | (24<<8))
const SQLITE_LOCKED_SHAREDCACHE =      (SQLITE_LOCKED |  (1<<8))
const SQLITE_BUSY_RECOVERY =           (SQLITE_BUSY   |  (1<<8))
const SQLITE_CANTOPEN_NOTEMPDIR =      (SQLITE_CANTOPEN | (1<<8))
const SQLITE_CANTOPEN_ISDIR =          (SQLITE_CANTOPEN | (2<<8))
const SQLITE_CANTOPEN_FULLPATH =       (SQLITE_CANTOPEN | (3<<8))
const SQLITE_CORRUPT_VTAB =            (SQLITE_CORRUPT | (1<<8))
const SQLITE_READONLY_RECOVERY =       (SQLITE_READONLY | (1<<8))
const SQLITE_READONLY_CANTLOCK =       (SQLITE_READONLY | (2<<8))
const SQLITE_READONLY_ROLLBACK =       (SQLITE_READONLY | (3<<8))
const SQLITE_ABORT_ROLLBACK =          (SQLITE_ABORT | (2<<8))
const SQLITE_CONSTRAINT_CHECK =        (SQLITE_CONSTRAINT | (1<<8))
const SQLITE_CONSTRAINT_COMMITHOOK =   (SQLITE_CONSTRAINT | (2<<8))
const SQLITE_CONSTRAINT_FOREIGNKEY =   (SQLITE_CONSTRAINT | (3<<8))
const SQLITE_CONSTRAINT_FUNCTION =     (SQLITE_CONSTRAINT | (4<<8))
const SQLITE_CONSTRAINT_NOTNULL =      (SQLITE_CONSTRAINT | (5<<8))
const SQLITE_CONSTRAINT_PRIMARYKEY =   (SQLITE_CONSTRAINT | (6<<8))
const SQLITE_CONSTRAINT_TRIGGER =      (SQLITE_CONSTRAINT | (7<<8))
const SQLITE_CONSTRAINT_UNIQUE =       (SQLITE_CONSTRAINT | (8<<8))
const SQLITE_CONSTRAINT_VTAB =         (SQLITE_CONSTRAINT | (9<<8))
const SQLITE_NOTICE_RECOVER_WAL =      (SQLITE_NOTICE | (1<<8))
const SQLITE_NOTICE_RECOVER_ROLLBACK = (SQLITE_NOTICE | (2<<8))
#Text Encodings
const SQLITE_UTF8 =            1 #
const SQLITE_UTF16LE =         2 #
const SQLITE_UTF16BE =         3 #
const SQLITE_UTF16 =           4 #    /* Use native byte order */
const SQLITE_ANY =             5 #    /* DEPRECATED */
const SQLITE_UTF16_ALIGNED =   8 #    /* sqlite3_create_collation only */

#Fundamental Data Types
const SQLITE_INTEGER = 1
const SQLITE_FLOAT   = 2
const SQLITE_TEXT    = 3
const SQLITE_BLOB    = 4
const SQLITE_NULL    = 5

const SQLITE3_TEXT   = 3

#Checkpoint operation parameters
const SQLITE_CHECKPOINT_PASSIVE =  0 #
const SQLITE_CHECKPOINT_FULL =     1 #
const SQLITE_CHECKPOINT_RESTART =  2 #

#Configuration Options
const SQLITE_CONFIG_SINGLETHREAD =   1 #  /* nil */
const SQLITE_CONFIG_MULTITHREAD =    2 #  /* nil */
const SQLITE_CONFIG_SERIALIZED =     3 #  /* nil */
const SQLITE_CONFIG_MALLOC =         4 #  /* sqlite3_mem_methods* */
const SQLITE_CONFIG_GETMALLOC =      5 #  /* sqlite3_mem_methods* */
const SQLITE_CONFIG_SCRATCH =        6 #  /* void*, int sz, int N */
const SQLITE_CONFIG_PAGECACHE =      7 #  /* void*, int sz, int N */
const SQLITE_CONFIG_HEAP =           8 #  /* void*, int nByte, int min */
const SQLITE_CONFIG_MEMSTATUS =      9 #  /* boolean */
const SQLITE_CONFIG_MUTEX =         10 #  /* sqlite3_mutex_methods* */
const SQLITE_CONFIG_GETMUTEX =      11 #  /* sqlite3_mutex_methods* */
#/* previously SQLITE_CONFIG_CHUNKALLOC 12 which is now unused. */
const SQLITE_CONFIG_LOOKASIDE =     13 #  /* int int */
const SQLITE_CONFIG_PCACHE =        14 #  /* no-op */
const SQLITE_CONFIG_GETPCACHE =     15 #  /* no-op */
const SQLITE_CONFIG_LOG =           16 #  /* xFunc, void* */
const SQLITE_CONFIG_URI =           17 #  /* int */
const SQLITE_CONFIG_PCACHE2 =       18 #  /* sqlite3_pcache_methods2* */
const SQLITE_CONFIG_GETPCACHE2 =    19 #  /* sqlite3_pcache_methods2* */
const SQLITE_CONFIG_COVERING_INDEX_SCAN =  20 #  /* int */
const SQLITE_CONFIG_SQLLOG =        21 #  /* xSqllog, void* */
const SQLITE_CONFIG_MMAP_SIZE =     22 #  /* sqlite3_int64, sqlite3_int64 */

#Database Connection Configuration Options
const SQLITE_DBCONFIG_LOOKASIDE =        1001 #  /* void* int int */
const SQLITE_DBCONFIG_ENABLE_FKEY =      1002 #  /* int int* */
const SQLITE_DBCONFIG_ENABLE_TRIGGER =   1003 #  /* int int* */

#Status Parameters for database connections
const SQLITE_DBSTATUS_LOOKASIDE_USED =        0 #
const SQLITE_DBSTATUS_CACHE_USED =            1 #
const SQLITE_DBSTATUS_SCHEMA_USED =           2 #
const SQLITE_DBSTATUS_STMT_USED =             3 #
const SQLITE_DBSTATUS_LOOKASIDE_HIT =         4 #
const SQLITE_DBSTATUS_LOOKASIDE_MISS_SIZE =   5 #
const SQLITE_DBSTATUS_LOOKASIDE_MISS_FULL =   6 #
const SQLITE_DBSTATUS_CACHE_HIT =             7 #
const SQLITE_DBSTATUS_CACHE_MISS =            8 #
const SQLITE_DBSTATUS_CACHE_WRITE =           9 #
const SQLITE_DBSTATUS_MAX =                   9 #   /* Largest defined DBSTATUS */

#Authorizer Return Codes
const SQLITE_DENY =    1 #   /* Abort the SQL statement with an error */
const SQLITE_IGNORE =  2 #   /* Don't allow access, but don't generate an error */

#Conflict resolution modes
const SQLITE_ROLLBACK =  1 #
#/* const SQLITE_IGNORE =  2 # // Also used by sqlite3_authorizer() callback */
const SQLITE_FAIL =      3 #
#/* const SQLITE_ABORT =  4 #  // Also an error code */
const SQLITE_REPLACE =   5 #

#Standard File Control Opcodes
const SQLITE_FCNTL_LOCKSTATE =                1 #
const SQLITE_GET_LOCKPROXYFILE =              2 #
const SQLITE_SET_LOCKPROXYFILE =              3 #
const SQLITE_LAST_ERRNO =                     4 #
const SQLITE_FCNTL_SIZE_HINT =                5 #
const SQLITE_FCNTL_CHUNK_SIZE =               6 #
const SQLITE_FCNTL_FILE_POINTER =             7 #
const SQLITE_FCNTL_SYNC_OMITTED =             8 #
const SQLITE_FCNTL_WIN32_AV_RETRY =           9 #
const SQLITE_FCNTL_PERSIST_WAL =             10 #
const SQLITE_FCNTL_OVERWRITE =               11 #
const SQLITE_FCNTL_VFSNAME =                 12 #
const SQLITE_FCNTL_POWERSAFE_OVERWRITE =     13 #
const SQLITE_FCNTL_PRAGMA =                  14 #
const SQLITE_FCNTL_BUSYHANDLER =             15 #
const SQLITE_FCNTL_TEMPFILENAME =            16 #
const SQLITE_FCNTL_MMAP_SIZE =               18 #

#Device Characteristics
const SQLITE_IOCAP_ATOMIC =                  0x00000001 #
const SQLITE_IOCAP_ATOMIC512 =               0x00000002 #
const SQLITE_IOCAP_ATOMIC1K =                0x00000004 #
const SQLITE_IOCAP_ATOMIC2K =                0x00000008 #
const SQLITE_IOCAP_ATOMIC4K =                0x00000010 #
const SQLITE_IOCAP_ATOMIC8K =                0x00000020 #
const SQLITE_IOCAP_ATOMIC16K =               0x00000040 #
const SQLITE_IOCAP_ATOMIC32K =               0x00000080 #
const SQLITE_IOCAP_ATOMIC64K =               0x00000100 #
const SQLITE_IOCAP_SAFE_APPEND =             0x00000200 #
const SQLITE_IOCAP_SEQUENTIAL =              0x00000400 #
const SQLITE_IOCAP_UNDELETABLE_WHEN_OPEN =   0x00000800 #
const SQLITE_IOCAP_POWERSAFE_OVERWRITE =     0x00001000 #

#Run-Time Limit Categories
const SQLITE_LIMIT_LENGTH =                     0 #
const SQLITE_LIMIT_SQL_LENGTH =                 1 #
const SQLITE_LIMIT_COLUMN =                     2 #
const SQLITE_LIMIT_EXPR_DEPTH =                 3 #
const SQLITE_LIMIT_COMPOUND_SELECT =            4 #
const SQLITE_LIMIT_VDBE_OP =                    5 #
const SQLITE_LIMIT_FUNCTION_ARG =               6 #
const SQLITE_LIMIT_ATTACHED =                   7 #
const SQLITE_LIMIT_LIKE_PATTERN_LENGTH =        8 #
const SQLITE_LIMIT_VARIABLE_NUMBER =            9 #
const SQLITE_LIMIT_TRIGGER_DEPTH =             10 #

#File Locking Levels
const SQLITE_LOCK_NONE =           0 #
const SQLITE_LOCK_SHARED =         1 #
const SQLITE_LOCK_RESERVED =       2 #
const SQLITE_LOCK_PENDING =        3 #
const SQLITE_LOCK_EXCLUSIVE =      4 #

#Mutex Types
const SQLITE_MUTEX_FAST =              0 #
const SQLITE_MUTEX_RECURSIVE =         1 #
const SQLITE_MUTEX_STATIC_MASTER =     2 #
const SQLITE_MUTEX_STATIC_MEM =        3 #  /* sqlite3_malloc() */
const SQLITE_MUTEX_STATIC_MEM2 =       4 #  /* NOT USED */
const SQLITE_MUTEX_STATIC_OPEN =       4 #  /* sqlite3BtreeOpen() */
const SQLITE_MUTEX_STATIC_PRNG =       5 #  /* sqlite3_random() */
const SQLITE_MUTEX_STATIC_LRU =        6 #  /* lru page list */
const SQLITE_MUTEX_STATIC_LRU2 =       7 #  /* NOT USED */
const SQLITE_MUTEX_STATIC_PMEM =       7 #  /* sqlite3PageMalloc() */

#Flags for the xShmLock VFS method
const SQLITE_SHM_UNLOCK =        1 #
const SQLITE_SHM_LOCK =          2 #
const SQLITE_SHM_SHARED =        4 #
const SQLITE_SHM_EXCLUSIVE =     8 #

#Constants Defining Special Destructor Behavior
# typedef void (*sqlite3_destructor_type)(void*);
const SQLITE_STATIC = reinterpret(Ptr{Void},0)
const SQLITE_TRANSIENT = reinterpret(Ptr{Void},-1)

#Function Flags
const SQLITE_DETERMINISTIC = 0x800

#Maximum xShmLock index
const SQLITE_SHM_NLOCK =         8 #

#Status Parameters
const SQLITE_STATUS_MEMORY_USED =           0 #
const SQLITE_STATUS_PAGECACHE_USED =        1 #
const SQLITE_STATUS_PAGECACHE_OVERFLOW =    2 #
const SQLITE_STATUS_SCRATCH_USED =          3 #
const SQLITE_STATUS_SCRATCH_OVERFLOW =      4 #
const SQLITE_STATUS_MALLOC_SIZE =           5 #
const SQLITE_STATUS_PARSER_STACK =          6 #
const SQLITE_STATUS_PAGECACHE_SIZE =        7 #
const SQLITE_STATUS_SCRATCH_SIZE =          8 #
const SQLITE_STATUS_MALLOC_COUNT =          9 #

#Status Parameters for prepared statements
const SQLITE_STMTSTATUS_FULLSCAN_STEP =      1 #
const SQLITE_STMTSTATUS_SORT =               2 #
const SQLITE_STMTSTATUS_AUTOINDEX =          3 #

#Synchronization Type Flags
const SQLITE_SYNC_NORMAL =         0x00002 #
const SQLITE_SYNC_FULL =           0x00003 #
const SQLITE_SYNC_DATAONLY =       0x00010 #

#Testing Interface Operation Codes
const SQLITE_TESTCTRL_FIRST =                     5 #
const SQLITE_TESTCTRL_PRNG_SAVE =                 5 #
const SQLITE_TESTCTRL_PRNG_RESTORE =              6 #
const SQLITE_TESTCTRL_PRNG_RESET =                7 #
const SQLITE_TESTCTRL_BITVEC_TEST =               8 #
const SQLITE_TESTCTRL_FAULT_INSTALL =             9 #
const SQLITE_TESTCTRL_BENIGN_MALLOC_HOOKS =      10 #
const SQLITE_TESTCTRL_PENDING_BYTE =             11 #
const SQLITE_TESTCTRL_ASSERT =                   12 #
const SQLITE_TESTCTRL_ALWAYS =                   13 #
const SQLITE_TESTCTRL_RESERVE =                  14 #
const SQLITE_TESTCTRL_OPTIMIZATIONS =            15 #
const SQLITE_TESTCTRL_ISKEYWORD =                16 #
const SQLITE_TESTCTRL_SCRATCHMALLOC =            17 #
const SQLITE_TESTCTRL_LOCALTIME_FAULT =          18 #
const SQLITE_TESTCTRL_EXPLAIN_STMT =             19 #
const SQLITE_TESTCTRL_LAST =                     19 #

#Virtual Table Configuration Options
const SQLITE_VTAB_CONSTRAINT_SUPPORT =  1 #

max_stars_repo_name: beffaCo/updateSqlite

max_issues_repo_name: beffaCo/updateSqlite

max_forks_repo_name: beffaCo/updateSqlite


text:
function optimize(obj::LeastSquaresOptimizer, arg0::LeastSquaresProblem)
    return jcall(obj, "optimize", LeastSquaresOptimizer_Optimum, (LeastSquaresProblem,), arg0)
end


text:
module BasisTest
using Base.Test
using Quantum

xb = basis(:X,[1:5])
xbplus = basis(:X,[2:6])
xb2 = xb*xb
@test separate(xb2)==[xb,xb]
@test separate(xb)==[xb]
@test basis([Ket{:X,Int}(i) for i=1:5])==xb
@test basis([xb[j]*xb[i] for i=1:5, j=1:5])==xb2
@test isdual(xb, xb')
@test basis(:X, [1,3,5]) == filter(s->label(s)%2!=0, xb)
@test basis(xb[1]*xb[1]) == filter(s->s==xb[1]*xb[1], xb2)
@test map(s->ket(:X,label(s)+1), xb) == basis(:X,[2:6])
@test map(s->ket(:X,label(s[1])+1)*s[2], xb2) == xbplus*xb
@test setdiff(xbplus, xb) == [xbplus[5]]
@test xb+xb == xb
@test xbplus[5]+xb == basis(:X,[6, 1:5])
@test xb+xbplus[5] == basis(:X,[1:6])

end


text:
abstract type AbstractState end

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

"A StateVector holds both the functions which will populate the state, and the most recent state."
mutable struct StateVector{S} <: AbstractState
	queries::Vector{Function}
	state::Vector{S}
	names::Vector{String}
end

function StateVector(queries::AbstractVector{Function}; names=fill("",length(queries)))
	StateVector(queries, zeros(length(queries)), names)
end
function StateVector(queries::Function...; names=fill("",length(queries)))
	StateVector(Function[f for f in queries], names=names)
end

Base.length(sv::StateVector) = length(sv.queries)

"retreive the most recently calculated state"
state(sv::StateVector) = sv.state

"update the state, then return it"
function state!(sv::StateVector)
	for (i,f) in enumerate(sv.queries)
		sv.state[i] = f()
	end
	sv.state
end

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

mutable struct History{T}
    sv::StateVector{T}
    states::Matrix{T}
end
History(sv::StateVector{T}) where {T} = History(sv, Matrix{T}(length(sv),0))

function state!(hist::History)
	s = state!(hist.sv)
	hist.states = hcat(hist.states, s)
	s
end

StatsBase.nobs(hist::History) = size(hist.states, 2)

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


text:
#!/usr/bin/env julia

#Benchmark iterative methods against matrices from the University of Florida
#sparse collection

#Root URL to the matrix collection
UFL_URL_ROOT = "http://www.cise.ufl.edu/research/sparse/matrices"
#Plain text file containing list of matrices
MASTERLIST = "matrices.txt"
#Download UFL collection to this directory
BASEDIR = "florida"

# 1. Read in master list of matrices from BASEDIR/MASTERLIST
#    If absent, generate this file. Requires Gumbo.jl
Pkg.installed("Gumbo")==nothing || using Gumbo
isdir(BASEDIR) || mkdir(BASEDIR)
if !isfile(joinpath(BASEDIR, MASTERLIST))
    Pkg.installed("Gumbo")===nothing && error("Parsing list from UFL website requires Gumbo.jl")
    #Parse the list on the UFL website. Requires the Gumbo HTML parser

    dlfilename = joinpath(BASEDIR, "florida.html")
    isfile(dlfilename) || download(UFL_URL_ROOT*"/list_by_id.html", dlfilename)

    doc = parsehtml(open(readall, dlfilename))

    matrices = Tuple{AbstractString, AbstractString}[]
    for elem in preorder(doc.root)
        if isa(elem, HTMLElement) && tag(elem)==:a
            relurl = getattr(elem, "href")
            if endswith(relurl, ".mat")
                tokens = split(relurl, '/')
                group = tokens[end-1]
                name = split(tokens[end], '.')[1]

                push!(matrices, (group, name))
            end
        end
    end

    #Write master list
    open(joinpath(BASEDIR, MASTERLIST), "w") do io
        for (group, name) in matrices
            println(io, group, '\t', name)
        end
    end
else #Read list of matrices from saved file
    matrices = Tuple{AbstractString, AbstractString}[]
    open(joinpath(BASEDIR, MASTERLIST)) do io
        for l in eachline(io)
            tokens = split(l)
            push!(matrices, (tokens[1], tokens[2]))
        end
    end
end

info("Listing contains $(length(matrices)) matrices")

# 2. Download all the matrices
for (group, matrix) in matrices
    groupdir = joinpath(BASEDIR, group)
    isdir(groupdir) || mkdir(groupdir)

    if !isfile(joinpath(groupdir, matrix*".mat"))
        info("Downloading $group/$matrix.mat")
        download(joinpath(UFL_URL_ROOT, "..", "mat", group, matrix*".mat"),
            joinpath(BASEDIR, group, matrix*".mat"))
    end
end
info("Downloaded $(length(matrices)) matrices")


text:
using Raylib
using Test

@testset "Raylib.jl" begin
    include("core.jl")
    include("easing.jl")
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export WFST

struct WFST{W, A<:Arc{W}, D, I, O,
            IS <: AbstractDict{I,D}, OS <: AbstractDict{O,D},
            IF <: AbstractDict{D,W},
            S  <: AbstractVector{<:AbstractVector{A}}
           }
  states::S
  initials::IF # index of initial state and corresponding weights
  finals::IF   # index of final state and corresponding weights
  isym::IS     # input  symbol table idx::Int => sym::I
  osym::OS     # output symbol table idx::Int => sym::O
end

"""
  `WFST(isym, [osym]; W = TropicalWeight{Float32})`

Construct an empty Weighted Finite State Transducer (WFST)

* `isym` input symbol table (can be a dictionary or a vector)
* `osym` output symbol dictionary (optional)
* `W` weight type

If the symbol table are `AbstractDict{I,D}`, `D` must be an integer and `I` the type of the symbol.
"""
function WFST(isym::Dict{I,D}, osym::Dict{O,D}; W=TropicalWeight{Float32}) where {D<:Integer, I, O}
  A = Arc{W,D}
  states=Vector{Vector{A}}()
  initials=Dict{D,W}()
  finals=Dict{D,W}()
  WFST(states,initials,finals,isym,osym)
end

parse_sym(sym::AbstractArray) = Dict(s => i for (i,s) in enumerate(sym))
parse_sym(sym::AbstractDict)  = sym
parse_sym(sym) = throw(ErrorException("Symbol table must be either a Vector or a Dict"))

WFST(isym;kwargs...) = WFST(isym,isym; kwargs...)
WFST(isym,osym;kwargs...) = WFST(parse_sym(isym),parse_sym(osym);kwargs...)

# methods
export add_states!

"""
`add_states!(fst,N)`

Add `N` empty states to the FST.
"""
function add_states!(fst::WFST{W,A}, N::Int) where {W,A}
  s = get_states(fst)
  for i = 1:N
    push!(s, Vector{A}() )
  end
  return fst
end

export get_initial
"""
`get_initial(fst::WFST, [single=false])`

Returns the initial state id of the `fst`. If `single` is `true` the first initial state is returned.
"""
function get_initial(fst::WFST; single=true)
  if isempty(fst.initials)
    throw(ErrorException("No initial state is present"))
  end
  if single
    return first(keys(fst.initials))
  else
    return fst.initials
  end
end

export get_initial_weight
"""
`get_initial_weight(fst::WFST,i)`

Returns the weight of initial `i`-th state of `fst`.
"""
get_initial_weight(fst::WFST,i) = get_initial(fst;single=false)[i]

export isinitial
"""
`isinitial(fst::WFST, i)`

Check if the `i`-th state is an initial state.
"""
isinitial(fst::WFST,i) = i in keys(fst.initials)

export initial!
"""
`initial!(fst::WFST{W},i, [w=one(W)])`

Set the `i`-th state as a initial state. A weight `w` can also be provided.
"""
function initial!(fst::WFST{W}, i, w=one(W)) where {W}
  fst.initials[i] = W(w)
  return fst
end

export rm_initial!
"""
`rm_initial!(fst,i)`

Remove initial state labeling of the state `i`.
"""
rm_initial!(fst,i) = delete!(fst.initials,i)

export get_final
"""
`get_final(fst::WFST, [single=false])`

Returns the final state id of the `fst`. If `single` is `true` the first final state is returned.
"""
function get_final(fst::WFST; single=false)
  if isempty(fst.finals)
    throw(ErrorException("No final state is present"))
  end
  if single
    return first(keys(fst.finals))
  else
    return fst.finals
  end
end

export get_final_weight
"""
`get_final_weight(fst::WFST,i)`

Returns the weight of `i`-th final state of `fst`.
"""
get_final_weight(fst::WFST,i) = get_final(fst;single=false)[i]

export isfinal
"""
`isfinal(fst::WFST, i)`

Check if the `i`-th state is an final state.
"""
isfinal(fst::WFST,i) = i in keys(fst.finals)

export final!
"""
`final!(fst::WFST{W},i, [w=one(W)])`

Set the `i`-th state as a final state. A weight `w` can also be provided.
"""
function final!(fst::WFST{W}, i, w) where {W}
  fst.finals[i] = W(w)
  return fst
end
final!(fst::WFST{W},i) where {W} = final!(fst,i,one(W))

export rm_final!
"""
`rm_final!(fst,i)`

Remove final state labeling of the `i`-th state.
"""
rm_final!(fst::WFST,i) = delete!(fst.finals,i)

export get_states
"""
`get_states(fst::WFST)`

Returns the states the `fst`.
"""
get_states(fst::WFST) = fst.states

export add_arc!
"""
`add_arc!(fst, src, dest, ilabel, olabel[, w=one(W)])`

`add_arc!(fst, srcdest::Pair, ilabelolabel::Pair[, w=one(W)])`

Adds an arc from state `src` to state `dest` with input label `ilabel`, output label `olabel` and weight `w`.
Alternative notation utilizes `Pair`s.

If `w` is not provided this defaults to `one(W)` where `W` is the weight type of `fst`.

"""
function add_arc!(fst::WFST{W},
                  srcdest::Pair,
                  iolabel::Pair,
                  w=one(W)) where {W}
  src, dest = srcdest.first, srcdest.second
  ilabel, olabel = iolabel.first, iolabel.second
  add_arc!(fst,src,dest,ilabel,olabel,w)
end

function parse_label(label,sym)
  if label in keys(sym)
    idx = sym[label]
  elseif iseps(label)
    idx = 0
  else
    throw(ErrorException("Symbol $label not found in symbol table $sym"))
  end
  return idx
end

function add_arc!(fst::WFST{W,A},
                  src::D,dest::D,
                  ilabel, olabel, w=one(W)) where {W,D<:Integer,A<:Arc{W,D}}
  if max(src,dest) - length(fst) > 0
    add_states!(fst, max(src,dest) - length(fst) )  # eventually expand states
  end
  i = parse_label(ilabel, fst.isym)
  o = parse_label(olabel, fst.osym)
  arc = Arc{W,D}(i,o,W(w),dest)
  push_arc!(fst, src, arc)
  return fst
end

# Note that this avoids checking if symbol is already present in the fst
function push_arc!(fst::WFST{W,A,D}, src::D, arc::A) where {W,A,D}
  s = get_states(fst)
  push!(s[src], arc)
  return fst
end

count_arcs(fst::WFST) = isempty(get_states(fst)) ? 0 : sum(length(s) for s in get_states(fst))

Base.getindex(fst::WFST, i) = getindex(get_states(fst),i)

Base.iterate(fst::WFST, args...) = iterate(get_states(fst), args...)

Base.isempty(fst::WFST) = isempty(get_states(fst))

struct ArcIterator
  fst::WFST
end

function Base.iterate(iter::ArcIterator, it_st=(1,1))
  i, j = it_st     # state id, arc id
  s = get_states(iter.fst)[i]
  if j > length(s)
    while true
      i += 1 # next state
      if i > length(iter.fst) return nothing; end
      s = get_states(iter.fst)[i]
      if !isempty(s) break; end
    end
    j = 1  # reinit arc
  end
  return s[j], (i,j+1)
end

Base.length(iter::ArcIterator) = count_arcs(iter.fst)

export get_arcs
"""
`get_arcs(fst::WFST) = ArcIterator(fst)`

Returns an iterator that can be used to loop through all of the arcs of `fst`.
"""
get_arcs(fst::WFST) = ArcIterator(fst)

export get_weight
"""
`get_weight(fst::WFST,i)`

Returns the weight of the `i`-th state of `fst`. If there is no weight `nothing` is returned.
"""
function get_weight(fst::WFST,i)
  if isinitial(fst,i)
    return get_initial_weight(fst,i)
  elseif isfinal(fst,i)
    return get_final_weight(fst,i)
  else
    return nothing
  end
end

# arc inspectors
function get_ilabel(fst::WFST{W,A,D},s::D,a::D) where {W,A,D}
  return get_iisym(fst)[get_ilabel(fst[s][a])]
end
function get_olabel(fst::WFST{W,A,D},s::D,a::D) where {W,A,D}
  return get_iosym(fst)[get_olabel(fst[s][a])]
end

function (fst::WFST{W,A,D,I,O})(labels::Vector{I}) where {W,A,D,I,O}
  i = get_initial(fst)
  w = get_initial_weight(fst,i)
  olabels = Vector{O}()
  isym = get_isym(fst)
  osym = get_iosym(fst)
  cnt = 1
  while cnt <= length(labels)
    n = 0
    lab = labels[cnt]
    if iseps(lab)
      cnt += 1
    else
      for arc in fst[i]
        ilab = get_ilabel(arc)
        if ilab == isym[lab] || iseps(ilab)
          if !iseps(ilab)
            cnt += 1
          end
          olab = get_olabel(arc)
          if !iseps(olab)
            push!(olabels,osym[olab])
          end
          n = get_nextstate(arc)
          w *= get_weight(arc)
          break
        end
      end
      if n == 0 # there is no state with lab
        w = zero(W)
        break
      else
        i = n
      end
    end
  end

  # finished labels but there could still be ilab == ϵ
  while !isfinal(fst,i)
    n = 0
    for arc in fst[i]
      ilab = get_ilabel(arc)
      olab = get_olabel(arc)
      if iseps(ilab)
        if !iseps(olab)
          push!(olabels,osym[olab])
        end
        n = get_nextstate(arc)
        w *= get_weight(arc)
        break
      end
    end
    if n == 0 # there is no state with lab
      w = zero(W)
      break
    else
      i = n
    end
  end

  if isfinal(fst,i)
    w *= get_final_weight(fst,i)
  else
    w = zero(W)
  end
  return olabels, w
end


text:
# This file was copied from Transducers.jl
# which is available under an MIT license (see LICENSE).
using PkgBenchmark

function mkconfig(; kwargs...)
    return BenchmarkConfig(env = Dict("JULIA_NUM_THREADS" => "1"); kwargs...)
end

group_target = benchmarkpkg(
    dirname(@__DIR__),
    mkconfig(),
    resultfile = joinpath(@__DIR__, "result-target.json"),
)

group_baseline = benchmarkpkg(
    dirname(@__DIR__),
    mkconfig(id = "baseline"),
    resultfile = joinpath(@__DIR__, "result-baseline.json"),
)

judgement = judge(group_target, group_baseline)


text:
module Api

using HTTP, JSON3
using ..MangroveBase: ReviewInfo, ReviewContent, Reviews

"Fetch all possible data from the Mangrove Api."
all_data() = HTTP.get("https://api.mangrove.reviews/reviews?issuers=true&maresi_subjects=true").body |> JSON3.read

"Extract internal datastructures from the Api output."
function extract(all::JSON3.Object)::Reviews
  reviews = all[:reviews]
  Dict(
    ReviewInfo(r[:payload][:sub], r[:kid]) =>
      ReviewContent(r[:payload][:rating], get(r[:payload], :opinion, missing))
    for r in reviews if haskey(r[:payload], :rating)
  )
end

"Fetch and extract data."
function mangrove_reviews()::Reviews
  extract(all_data())
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CookieCutter(master, others)

A cookie-cutter simulation solver.

## Parameters

* `master` - Master simulation solver (a.k.a. facies solver)
* `others` - Solvers for each value of the master realization

## Examples

Simulate lithofacies with image quilting and fill porosity
values with LU Gaussian simulation:

```julia
julia> f  = IQ(:facies => (TI=IMG, template=(30,30,1)))
julia> p₀ = LUGS(:poro => (variogram=SphericalVariogram(range=10.),))
julia> p₁ = LUGS(:poro => (variogram=SphericalVariogram(range=20.),))
julia> CookieCutter(f, Dict(0=>p₀, 1=>p₁))
```
"""
struct CookieCutter <: SimulationSolver
  master::SimulationSolver
  others::Dict{Number,SimulationSolver}
end

function solve(problem::SimulationProblem, solver::CookieCutter)
  # retrieve problem info
  pdata   = data(problem)
  pdomain = domain(problem)
  preals  = nreals(problem)

  # master variable
  pvars = Dict(name(v) => mactype(v) for v in variables(problem))
  mvars = targets(solver.master)
  @assert length(mvars) == 1 "one single variable must be specified in master solver"
  mvar = mvars[1]
  @assert mvar ∈ keys(pvars) "invalid variable in master solver"
  mtype = pvars[mvar]

  # other variables
  ovars = Tuple(var => V for (var, V) in pvars if var ≠ mvar)
  @assert length(ovars) > 0 "cookie-cutter requires problem with more than one target variable"

  # define master and others problem
  if hasdata(problem)
    mproblem = SimulationProblem(pdata, pdomain, mvar, preals)
    oproblem = SimulationProblem(pdata, pdomain, first.(ovars), preals)
  else
    mproblem = SimulationProblem(pdomain, mvar => mtype, preals)
    oproblem = SimulationProblem(pdomain, ovars, preals)
  end

  # realizations of master variable
  msol = solve(mproblem, solver.master)
  mreals = msol[mvar]

  # pre-allocate memory for realizations
  reals = Dict(var => [Vector{V}(undef, nelements(pdomain)) for i in 1:preals] for (var,V) in ovars)

  # solve other problems
  for (mval, osolver) in solver.others
    osol = solve(oproblem, osolver)

    # cookie-cutter step
    for (var, V) in ovars
      vreals = osol[var]
      for i in 1:preals
        mask = mreals[i] .== mval
        reals[var][i][mask] .= vreals[i][mask]
      end
    end
  end

  realizations = merge(reals, Dict(mvar => mreals))

  Ensemble(pdomain, realizations)
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, solver::CookieCutter)
  print(io, "CookieCutter")
end

function Base.show(io::IO, ::MIME"text/plain", solver::CookieCutter)
  mvar = targets(solver.master)[1]
  println(io, solver)
  println(io, "  └─", mvar, " ⇨ ", solver.master)
  lines = ["    └─$val ⇨ $osolver" for (val,osolver) in solver.others]
  print(io, join(lines, "\n"))
end


text:
# abstract type AbstractAllocator end
# abstract type StaticAllocator <: AbstractAllocator end
# include("PseudoInverseAllocator.jl")
include("ConstrainedAllocator.jl")
include("AdaptiveAllocator.jl")


text:
module SeqLoggers

using Logging
using LoggingExtras
using Dates
using HTTP
using JSON3

include("utils.jl")

include("loggers.jl")
export SeqLogger
export run_with_logger

include("advanced_file_logger.jl")
export AdvancedFileLogger

include("logging_extras.jl")

include("load_from_config.jl")
export load_logger_from_config

end


text:
function test_reification_singlevariable(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T}) # c1
    @MOIT.requires MOI.supports_constraint(model, MOI.VectorOfVariables, CP.Reification{MOI.EqualTo{T}}) # c2

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x2, MOI.EqualTo(1))
    c2 = MOI.add_constraint(model, MOI.VectorOfVariables([x1, x2]), CP.Reification(MOI.EqualTo(2)))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if MOIT._supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

function MOIT.setup_test(
    ::typeof(test_reification_singlevariable),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[0, 1]),
    )
    return
end

function test_reification_scalaraffinefunction(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2
    @MOIT.requires MOI.supports_constraint(model, MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}) # c1
    @MOIT.requires MOI.supports_constraint(model, MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}) # c2

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, one(T) * x2, MOI.EqualTo(1))
    c2 = MOI.add_constraint(model, MOIU.vectorize(one(T) .* [x1, x2]), CP.Reification(MOI.EqualTo(2)))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if MOIT._supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

function MOIT.setup_test(
    ::typeof(test_reification_scalaraffinefunction),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[0, 1]),
    )
    return
end