text:
# Exact Search
# ============
#
# Exact sequence search tools.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Reference Implementation
# ------------------------
#
# `searchindex` and `rsearchindex` should work exactly the same way as the
# following reference implementation:
#=
# Return the index of the first occurrence of `pat` in `seq[start:stop]`.
function searchindex(seq, pat, start=1, stop=endof(seq))
    m = length(pat)
    n = length(seq)
    for s in max(start-1, 0):min(stop, n)-m
        if occurs_with_shift(pat, seq, s)
            return s+1  # found
        end
    end
    return 0  # not found
end

# Return the index of the last occurrence of `pat` in `seq[stop:start]`.
function rsearchindex(seq, pat, start=endof(seq), stop=1)
    n = length(seq)
    m = length(pat)
    for s in min(start, n)-m:-1:max(stop-1, 0)
        if occurs_with_shift(pat, seq, s)
            return s+1  # found
        end
    end
    return 0  # not found
end

function occurs_with_shift(pat, seq, shift)
    for i in 1:endof(pat)
        if !iscompatible(pat[i], seq[shift+i])
            return false
        end
    end
    return true
end
=#

# Actual Implementation
# ---------------------

"""
Query type for exact sequence search.
"""
immutable ExactSearchQuery{S<:Sequence}
    seq::S         # query sequence
    cbits::UInt32  # compatibility bits
    fshift::Int    # shift length for forward search
    bshift::Int    # shift length for backward search
end

function ExactSearchQuery(query::Sequence)
    cbits, fshift, bshift = preprocess(query)
    return ExactSearchQuery(query, cbits, fshift, bshift)
end

function preprocess(query)
    if length(query) == 0
        return UInt32(0), 0, 0
    end
    m = length(query)
    first = query[1]
    last = query[end]
    cbits::UInt32 = 0
    fshift = bshift = m
    for i in 1:endof(query)
        x = query[i]
        cbits |= compatbits(x)
        if iscompatible(x, last) && i < m
            fshift = m - i
        end
    end
    for i in endof(query):-1:1
        x = query[i]
        if iscompatible(x, first) && i > 1
            bshift = i - 1
        end
    end
    return cbits, fshift, bshift
end

function checkeltype(seq1, seq2)
    if eltype(seq1) != eltype(seq2)
        throw(ArgumentError("the element type of two sequences must match"))
    end
end

# Forward
# -------

# NOTE: the return value is a range while the counterpart for strings is an index
"""
    search(seq::Sequence, pat[, start=1[, stop=endof(seq)]])

Return the range of the first occurrence of `pat` in `seq[start:stop]`; symbol
comparison is done using `Bio.Seq.iscompatible`.
"""
function Base.search(seq::Sequence, val,
                     start::Integer=1, stop::Integer=endof(seq))
    s = searchindex(seq, val, start, stop)
    if s == 0
        return 0:-1
    else
        return s:s
    end
end

function Base.search(seq::Sequence, pat::Sequence,
                     start::Integer=1, stop::Integer=endof(seq))
    return search(seq, ExactSearchQuery(pat), start, stop)
end

function Base.search(seq::Sequence, query::ExactSearchQuery,
                     start::Integer=1, stop::Integer=endof(seq))
    s = searchindex(seq, query, start, stop)
    if s == 0
        return 0:-1
    end
    return s:s+length(query.seq)-1
end

"""
    searchindex(seq::Sequence, pat[, start=1[, stop=endof(seq)]])

Return the index of the first occurrence of `pat` in `seq[start:stop]`; symbol
comparison is done using `Bio.Seq.iscompatible`.
"""
function Base.searchindex(seq::Sequence, val,
                          start::Integer=1, stop::Integer=endof(seq))
    x = convert(eltype(seq), val)
    for i in max(start, 1):min(stop, endof(seq))
        if iscompatible(seq[i], x)
            return i  # found
        end
    end
    return 0  # not found
end

function Base.searchindex(seq::Sequence, pat::Sequence,
                          start::Integer=1, stop::Integer=endof(seq))
    return searchindex(seq, ExactSearchQuery(pat), start, stop)
end

function Base.searchindex(seq::Sequence, query::ExactSearchQuery,
                          start::Integer=1, stop::Integer=endof(seq))
    checkeltype(seq, query.seq)
    return quicksearch(query, seq, start, stop)
end

# This algorithm is borrowed from base/strings/search.jl, which looks similar to
# Sunday's Quick Search algorithm.
function quicksearch(query, seq, start, stop)
    pat = query.seq
    m = length(pat)
    n = length(seq)
    stop′ = min(stop, n) - m
    s::Int = max(start - 1, 0)

    if m == 0  # empty query
        if s ≤ stop′
            return s + 1  # found
        else
            return 0  # not found
        end
    end

    while s ≤ stop′
        if iscompatible(pat[m], seq[s+m])
            i = m - 1
            while i > 0
                if !iscompatible(pat[i], seq[s+i])
                    break
                end
                i -= 1
            end
            if i == 0
                return s + 1  # found
            elseif s < stop′ && (query.cbits & compatbits(seq[s+m+1]) == 0)
                s += m + 1
            elseif isambiguous(seq[s+m])
                s += 1
            else
                s += query.fshift
            end
        elseif s < stop′ && (query.cbits & compatbits(seq[s+m+1]) == 0)
            s += m + 1
        else
            s += 1
        end
    end

    return 0  # not found
end

# Backward
# --------

# NOTE: the return value is a range while the counterpart for strings is an index
"""
    rsearch(seq::Sequence, pat[, start=endof(seq)[, stop=1]])

Return the range of the last occurrence of `pat` in `seq[stop:start]`; symbol
comparison is done using `Bio.Seq.iscompatible`.
"""
function Base.rsearch(seq::Sequence, val,
                      start::Integer=endof(seq), stop::Integer=1)
    s = rsearchindex(seq, val, start, stop)
    if s == 0
        return 0:-1
    else
        return s:s
    end
end

function Base.rsearch(seq::Sequence, pat::Sequence,
                      start::Integer=endof(seq), stop::Integer=1)
    return rsearch(seq, ExactSearchQuery(pat), start, stop)
end

function Base.rsearch(seq::Sequence, query::ExactSearchQuery,
                      start::Integer=endof(seq), stop::Integer=1)
    s = rsearchindex(seq, query, start, stop)
    if s == 0
        return 0:-1
    end
    return s:s+length(query.seq)-1
end

"""
    rsearchindex(seq::Sequence, pat[, start=1[, stop=endof(seq)]])

Return the index of the last occurrence of `pat` in `seq[start:stop]`; symbol
comparison is done using `Bio.Seq.iscompatible`.
"""
function Base.rsearchindex(seq::Sequence, val,
                           start::Integer=endof(seq), stop::Integer=1)
    x = convert(eltype(seq), val)
    for i in min(start, endof(seq)):-1:max(stop, 1)
        if iscompatible(seq[i], x)
            return i  # found
        end
    end
    return 0  # not found
end

function Base.rsearchindex(seq::Sequence, pat::Sequence,
                           start::Integer=endof(seq), stop::Integer=1)
    return rsearchindex(seq, ExactSearchQuery(pat), start, stop)
end

function Base.rsearchindex(seq::Sequence, query::ExactSearchQuery,
                           start::Integer=endof(seq), stop::Integer=1)
    checkeltype(seq, query.seq)
    return quickrsearch(seq, query, start, stop)
end

function quickrsearch(seq, query, start, stop)
    pat = query.seq
    m = length(pat)
    n = length(seq)
    stop′ = max(stop - 1, 0)
    s::Int = min(start, n) - m

    if m == 0  # empty query
        if s ≥ stop′
            return s + 1  # found
        else
            return 0  # not found
        end
    end

    while s ≥ stop′
        if iscompatible(pat[1], seq[s+1])
            i = 2
            while i < m + 1
                if !iscompatible(pat[i], seq[s+i])
                    break
                end
                i += 1
            end
            if i == m + 1
                return s + 1  # found
            elseif s > stop′ && (query.cbits & compatbits(seq[s]) == 0)
                s -= m + 1
            elseif isambiguous(seq[s+1])
                s -= 1
            else
                s -= query.bshift
            end
        elseif s > stop′ && (query.cbits & compatbits(seq[s]) == 0)
            s -= m + 1
        else
            s -= 1
        end
    end

    return 0  # not found
end


text:
module SpatialMoran1D

    export  ETKomaP,
            ETKomaErfP,
            ETotalKoma,
            ETMRCAKoma,
            ESegSitesKoma,
            EPiKoma,
            PSegSitesMoran,
            ESegSitesMoran,
            TMoran,
            ETMoran,
            PiMoran

    using Distributions

    """
        ETKomaP(k, n)

    Expected value for the kth coalescent time under the spatial model with PBC.
    # Arguments
    * `k`: Lineages present
    * `n`: Population size
    """
    ETKomaP(k, n)   = n^3/pi * (1/(k-1)^2-1/k^2)*(1-1/k) + (1-2/pi)
    """
        ETKomaP(k, n, m)

    Expected value for the kth coalescent time under the spatial model with PBC.
    # Arguments
    * `k`: Lineages present
    * `n`: Population size
    * `m`: Tweak parameter
    """
    ETKomaP(k, n, m) = n^3/pi *(1/(k-m)^2-1/(k+1-m)^2)*((k-m)/(k-m+1))+(1-2/pi)
    """
        ETKomaErfP(k, n, C0::Real=1.)

    Expected value for the kth coalescent time under the spatial model with PBC.
    Exact solution via inverse of error function.
    # Arguments
    * `k`: Lineages present
    * `n`: Population size
    * `C0=1.0`: Initial concentration; 1.0 is fully occupied.
    """
    ETKomaErfP(k, n, C0::Real=1.) = 1/4/C0^2*(1/erfinv((k-1)/n)^2-1/erfinv(k/n)^2)*(1-1/k)*n
    """
        ETotalKoma(n)

    Expected value for total length of a genealogy under the spatial model.
    # Arguments
    * `n`: Population size
    """
    ETotalKoma(n) = sum( k->ETKomaP(k,n)*k, 2:n )
    """
        ETMRCAKoma(n)

    Expected value for time to the MRCA under the spatial model.
    # Arguments
    * `n`: Population size
    """
    ETMRCAKoma(n) = sum( k->ETKomaP(k,n), 2:n )
    """
        ESegSitesKoma(n, μ)

    Expected value for the number of segregating sites under the spatial model.
    # Arguments
    * `n`: Population size
    * `μ`: Mutation rate
    """
    ESegSitesKoma(n, μ) = μ*ETotalKoma(n)
    """
        EPiKoma(n, μ)

    Expected value for the number of pairwise differences under the spatial model.
    # Arguments
    * `n`: Population size
    * `μ`: Mutation rate
    """
    EPiKoma(n, μ) = μ/12*n^2*(n+1)

    """
        PSegSitesMoran(k, n, θ)

    Probability Mass Function of the number of segregating sites under Moran.
    # Arguments
    * `n`: Population size
    * `θ`: Mutation rate
    """
    PSegSitesMoran(k, n, θ) = sum( i->(-1)^i*binomial(n-1,i-1)*(i-1)/(θ+i-1)*(θ/(θ+i-1))^k, 2:n)

    """
        ESegSitesMoran(n, μ)

    Expected number of segregrating sites under Moran.
    # Arguments
    * `n`: Population size
    * `θ`: Mutation rate
    """
    ESegSitesMoran(n, μ) = n*μ*sum( i->1/i, 1:n-1)

    """
        TMoran(t, k, n)

    Distribution of coalescence times according to Kingman. Moran-timescale `2/n^2`
    # Arguments
    * `k`: Lineages present
    * `n`: Population size
    """
    TMoran(t, k, n) = pdf( Exponential(binomial(k,2)), t/n^2*2)

    """
        ETMoran(k, n) = n^2/2 / binomial(k,2)

    Expected value for the kth coalescence time under Moran.
    # Arguments
    * `k`: Lineages present
    * `n`: Population size
    """
    ETMoran(k, n) = n^2/2 / binomial(k,2)

    """
        PiMoran(n, θ)

    Expected number of pairwise differences under Moran.
    # Arguments
    * `θ`: Mutation rate
    * `n`: Population size
    """
    PiMoran(n, θ) = n*θ

end # module

ext: jl


text:
#
# A simple steepest descent algorithm implementation
#
export initializeSolver!, doSolverStep!, getSolverResult
export steepestDescent
export DebugGradient, DebugGradientNorm, DebugStepsize
@doc doc"""
    steepestDescent(M, F, ∇F, x)
perform a steepestDescent $x_{k+1} = \mathrm{retr}_{x_k} s_k\nabla f(x_k)$ with
different choices of $s_k$ available (see `stepsize` option below).

# Input
* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `∇F`: the gradient $\nabla F\colon\mathcal M\to T\mathcal M$ of F
* `x` – an initial value $x\in\mathcal M$

# Optional
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`) specify a [`Stepsize`](@ref)
  functor.
* `retraction` – (`exp`) a `retraction(M,x,ξ)` to use.
* `stoppingCriterion` – (`[`stopWhenAny`](@ref)`(`[`stopAfterIteration`](@ref)`(200), `[`stopWhenGradientNormLess`](@ref)`(10.0^-8))`)
  a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop.

and the ones that are passed to [`decorateOptions`](@ref) for decorators.

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
* `record` - if activated (using the `record` key, see [`RecordOptions`](@ref)
  an array containing the recorded values.
"""
function steepestDescent(M::mT,
    F::Function, ∇F::Function, x::MP;
    stepsize::Stepsize = ConstantStepsize(1.0),
    retraction::Function = exp,
    stoppingCriterion::StoppingCriterion = stopWhenAny( stopAfterIteration(200), stopWhenGradientNormLess(10.0^-8)),
    kwargs... #collect rest
  ) where {mT <: Manifold, MP <: MPoint}
  p = GradientProblem(M,F,∇F)
  o = GradientDescentOptions(x,stoppingCriterion,stepsize,retraction)
  o = decorateOptions(o; kwargs...)
  resultO = solve(p,o)
  if hasRecord(resultO)
      return getSolverResult(p,getOptions(resultO)), getRecord(resultO)
  end
  return getSolverResult(p,resultO)
end
#
# Solver functions
#
function initializeSolver!(p::P,o::O) where {P <: GradientProblem, O <: GradientDescentOptions}
    o.∇ = getGradient(p,o.x)
end
function doSolverStep!(p::P,o::O,iter) where {P <: GradientProblem, O <: GradientDescentOptions}
    o.∇ = getGradient(p,o.x)
    o.x = o.retraction(p.M, o.x , -getStepsize!(p,o,iter) * o.∇)
end
getSolverResult(p::P,o::O) where {P <: GradientProblem, O <: GradientDescentOptions} = o.x


text:
function test_deprecated()
  @testset ExtendedTestSet "Deprecated methods" begin
    n = 10
    mem = 3
    T = Float16
    @test_deprecated LBFGSOperator(n, mem, scaling = false)
    @test_deprecated InverseLBFGSOperator(n, mem, scaling = false)
    @test_deprecated LBFGSOperator(T, n, mem, scaling = false)
    @test_deprecated InverseLBFGSOperator(T, n, mem, scaling = false)
    @test_deprecated LSR1Operator(n, mem, scaling = false)
    @test_deprecated LSR1Operator(T, n, mem, scaling = false)
  end
end

test_deprecated()


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

mutable struct ApplicationGatewayPrivateLinkResourceProperties <: SwaggerModel
    groupId::Any # spec type: Union{ Nothing, String } # spec name: groupId
    requiredMembers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: requiredMembers
    requiredZoneNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: requiredZoneNames

    function ApplicationGatewayPrivateLinkResourceProperties(;groupId=nothing, requiredMembers=nothing, requiredZoneNames=nothing)
        o = new()
        validate_property(ApplicationGatewayPrivateLinkResourceProperties, Symbol("groupId"), groupId)
        setfield!(o, Symbol("groupId"), groupId)
        validate_property(ApplicationGatewayPrivateLinkResourceProperties, Symbol("requiredMembers"), requiredMembers)
        setfield!(o, Symbol("requiredMembers"), requiredMembers)
        validate_property(ApplicationGatewayPrivateLinkResourceProperties, Symbol("requiredZoneNames"), requiredZoneNames)
        setfield!(o, Symbol("requiredZoneNames"), requiredZoneNames)
        o
    end
end # type ApplicationGatewayPrivateLinkResourceProperties

const _property_map_ApplicationGatewayPrivateLinkResourceProperties = Dict{Symbol,Symbol}(Symbol("groupId")=>Symbol("groupId"), Symbol("requiredMembers")=>Symbol("requiredMembers"), Symbol("requiredZoneNames")=>Symbol("requiredZoneNames"))
const _property_types_ApplicationGatewayPrivateLinkResourceProperties = Dict{Symbol,String}(Symbol("groupId")=>"String", Symbol("requiredMembers")=>"Vector{String}", Symbol("requiredZoneNames")=>"Vector{String}")
Base.propertynames(::Type{ ApplicationGatewayPrivateLinkResourceProperties }) = collect(keys(_property_map_ApplicationGatewayPrivateLinkResourceProperties))
Swagger.property_type(::Type{ ApplicationGatewayPrivateLinkResourceProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPrivateLinkResourceProperties[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPrivateLinkResourceProperties }, property_name::Symbol) =  _property_map_ApplicationGatewayPrivateLinkResourceProperties[property_name]

function check_required(o::ApplicationGatewayPrivateLinkResourceProperties)
    true
end

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


text:
# inspired by
# - https://github.com/dalum/InlineExports.jl/blob/master/src/InlineExports.jl
# - https://github.com/QuantumBFS/YaoBase.jl/blob/master/src/utils/interface.jl

# TODO: use MLStyle.jl instead!

const API = []
const ENV_API = []
const MULTI_AGENT_ENV_API = []

macro api(ex)
    interfacem(__module__, __source__, ex, API)
end

macro env_api(ex)
    interfacem(__module__, __source__, ex, ENV_API)
end

macro multi_agent_env_api(ex)
    interfacem(__module__, __source__, ex, MULTI_AGENT_ENV_API)
end

function interfacem(__module__::Module, __source__::LineNumberNode, ex::Expr, store)
    name, is_body_missing = handle(ex)
    if name === nothing
        :(error("unknown expression"))
    else
        push!(store, name)
        if is_body_missing
            quote
                export $(esc(name))
                Core.@__doc__ $(esc(ex)) = error("method not implemented")
            end
        else
            quote
                export $(esc(name))
                Core.@__doc__ $(esc(ex))
            end
        end
    end
end

handle(ex) = extract_name(ex), is_body_missing(ex)

extract_name(::Any) = nothing
extract_name(x::Symbol) = x
extract_name(x::QuoteNode) = x.value
extract_name(ex::Expr) = extract_name(Val(ex.head), ex)
extract_name(::Val{:abstract}, ex) = extract_name(ex.args[1])
extract_name(::Val{:(=)}, ex) = extract_name(ex.args[1])
extract_name(::Val{:call}, ex) = extract_name(ex.args[1])
extract_name(::Val{:(::)}, ex) = extract_name(ex.args[2])
extract_name(::Val{:(<:)}, ex) = extract_name(ex.args[1])
extract_name(::Val{:curly}, ex) = extract_name(ex.args[1])
extract_name(::Val{:.}, ex) = extract_name(ex.args[2])
extract_name(::Val{:where}, ex) = extract_name(ex.args[1])
extract_name(::Val{:function}, ex) = extract_name(ex.args[1])
extract_name(::Val{:struct}, ex) = extract_name(ex.args[2])
extract_name(::Val{:const}, ex) = extract_name(ex.args[1])

is_body_missing(::Any) = false

is_body_missing(ex::Expr) = is_body_missing(Val(ex.head), ex)
is_body_missing(::Any, ::Any) = false
is_body_missing(::Val{:call}, ::Any) = true


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using Base.Test

import FEMBase: solve!, can_solve

mutable struct LSSolver1 <: AbstractLinearSystemSolver
    a :: Int64
end

mutable struct LSSolver2 <: AbstractLinearSystemSolver
end

function LSSolver1()
    return LSSolver1(1)
end

function can_solve(ls::LinearSystem, solver::LSSolver1)
    return (false, "don't know how to do")
end

function solve!(ls::LinearSystem, solver::LSSolver2)
    fill!(ls.u, 1.0)
end

@testset "test creating new linear system solver" begin
    ls = LinearSystem(3)
    ls.u = spzeros(3)
    solver1 = LSSolver1()
    solver2 = LSSolver2()
    solve!(ls, solver1)
    @test can_solve(ls, solver1)[1] == false
    solve!(ls, [solver1, solver2])
    @test isapprox(ls.u, sparse(ones(3)))
    @test_throws ErrorException solve!(ls, [solver1])
end


text:
function time_integration_explicit(dt,Qbase_hat,RHS,cellxmax,cellymax)
    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax-1
            Qbase_hat[i,j] = Qbase_hat[i,j] + dt*RHS[i,j]
        end
    end
    return Qbase_hat
end


text:
#something
F(x,y) = 2x+3y;


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

name = "SQLiteBuilder"
version = v"0.1.0"

# Collection of sources required to build SQLiteBuilder
sources = [
    "https://www.sqlite.org/2018/sqlite-autoconf-3240000.tar.gz" =>
    "d9d14e88c6fb6d68de9ca0d1f9797477d82fc3aed613558f87ffbdbbc5ceb74a",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd sqlite-autoconf-3240000/
./configure --prefix=$prefix --host=$target && make && make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = BinaryBuilder.supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libsqlite3", :libsqlite)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using CoordinateTransformations

@testset "Reading rasters" begin
    ga = GeoArrays.read("data/utmsmall.tif")
    @test bbox(ga) == (min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6)
    @test bboxes(ga)[1] == (min_x = 440720.0, max_x = 440780.0, min_y = 3.75126e6, max_y = 3.75132e6)
    @test bboxes(ga)[end] == (min_x = 446660.0, max_x = 446720.0, min_y = 3.74532e6, max_y = 3.74538e6)
end

@testset "coords" begin
    straight = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 -1.0], [375000.03, 380000.03]), "")
    rot = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.5; 0.1 1.0], [0.0, 0.0]), "")
    @test coords(straight, :x) == collect(375000.03:1:375005.03)
    @test coords(straight, :y) == collect(380000.03:-1:379995.03)
    @test_throws ErrorException coords(straight, :z)
    @test_throws ErrorException coords(rot, :x)
    @test_throws ErrorException coords(rot, :y)
end

@testset "GeoArray constructors" begin
    x, y = range(4, stop = 8.0, length = 10), range(0, stop = 1, length = 9)
    ga2 = GeoArray(rand(10, 9), x, y)
    ga2 = GeoArray(rand(10, 9), x, y, "")
    ga3 = GeoArray(rand(10, 9, 8), x, y)
    ga3 = GeoArray(rand(10, 9, 8), x, y, "")
    for i in 1:length(x), j in 1:length(y)
        @test GeoArrays.centercoords(ga2, [i,j]) ≈ [x[i],y[j]]
    end
    for i in 1:length(x), j in 1:length(y)
        @test GeoArrays.centercoords(ga3, [i,j]) ≈ [x[i],y[j]]
    end
    x, y = range(4, stop = 8.0, length = 11), range(0, stop = 1, length = 9)
    @test_throws ErrorException GeoArray(rand(10, 9), x, y)
end

@testset "Conversions" begin
    ga = GeoArray(rand(1:32000, 5, 5))
    GeoArrays.write!("test_conversion.tif", ga)

    ga = GeoArray(rand(Bool, 5, 5))
    @test_throws ErrorException GeoArrays.write!("test_conversion.tif", ga)
end


text:
# This file is part of Kpax3. License is MIT.

# normalized distance between a generic partition 'a' and the partition
# {{1, ..., v}}
function distsjlower(n::Vector{T},
                     k::Int,
                     v::Real;
                     β::Real=3.0) where T <: Real
  w = v^β

  z = 0.0
  for g in 1:k
    z += n[g]^β
  end

  (w - z) / (w - v)
end

# normalized distance between a generic partition 'a' and the partition
# {{1}, ..., {v}}
function distsjupper(n::Vector{T},
                     k::Int,
                     v::Real;
                     β::Real=3.0) where T <: Real
  w = v^β

  z = 0.0
  for g in 1:k
    z += n[g]^β
  end

  (z - v) / (w - v)
end

# normalized distance between partition 'a' and partition 'b'
function distsj(a::Vector{Int},
                na::Vector{T},
                ka::Int,
                b::Vector{Int},
                nb::Vector{T},
                kb::Int,
                v::Int;
                β::Real=3.0) where T <: Real
  w = v^β

  count = zeros(Float64, ka, kb)
  for i in 1:v
    count[a[i], b[i]] += 1.0
  end

  za = 0.0
  zb = 0.0
  zab = 0.0

  for g in 1:ka
    za += na[g]^β
  end

  for g in 1:kb
    zb += nb[g]^β
  end

  for g in 1:kb, h in 1:ka
    zab += count[h, g]^β
  end

  (za + zb - 2 * zab) / (w - v)
end


text:
export color_material_physical_group_entities

"""
    color_material_physical_group_entities(materials::Vector{Material})

For a model with physical groups of the form "Material: X, Material: Y, ...", color
the entities in each material physical group according to corresponding Material in `materials`
"""
function color_material_physical_group_entities(materials::Vector{Material})
    material_names = ["Material: "*mat.name for mat in materials]
    for group in gmsh.model.get_physical_groups()
        gdim, gnum = group
        name = gmsh.model.get_physical_name(gdim, gnum)
        if startswith(name, "Material: ")
            tags = gmsh.model.get_entities_for_physical_group(gdim, gnum)
            mat_id = findfirst(mat_name->name == mat_name, material_names)
            color = materials[mat_id].color
            r = Int32(color.r.i)
            g = Int32(color.g.i)
            b = Int32(color.b.i)
            a = Int32(color.alpha.i)
            dtags = [(gdim, tag) for tag in tags]
            gmsh.model.set_color(dtags, r, g, b, a)
        end
    end
    return nothing
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
MPI.Barrier(comm)
comm2 = MPI.Comm_dup(comm)
MPI.Barrier(comm2)
comm3 = MPI.Comm_dup(comm2)
MPI.Barrier(comm3)
MPI.Comm_free(comm2)
MPI.Barrier(comm3)
# Don't free comm2

MPI.Finalize()


text:
#
# tools.jl --
#
# Tools for installing `AndorCameras` package: `parse_header` is a simple
# parser to fetch definitions from Andor SDK header file <atcore.h> and
# `make_deps` is called to generate "deps.jl".
#
module AndorInstallTools

using Libdl

# List of expected macros (except error code):
const MACROS = Dict(
    "INFINITE" => "Cuint",
    "TRUE" => "BOOL",
    "FALSE" => "BOOL",
    "HANDLE_UNINITIALISED" => "Handle",
    "HANDLE_SYSTEM" => "Handle",
    "SUCCESS" => "Cint",
    "CALLBACK_SUCCESS" => "Cint",
)

function parse_header(filename::String)
    code = ["# Constants.",]
    tail = String[]
    open(filename, "r") do io
        linenumber = 0
        while !eof(io)
            line = readline(io; keep=false)
            linenumber += 1
            m = match(r"^ *# *define +AT_([^ ]+) *([^ ]+) *$", line)
            m === nothing && continue
            name, value = m.captures
            if haskey(MACROS, name)
                type = MACROS[name]
            elseif startswith(name, "ERR_")
                type = "Cint"
            else
                continue
            end
            if tryparse(Int, value) === nothing
                println(stderr, "Warning: parsed value \"$value\" ",
                        "cannot be converted to an integer constant ",
                        "($filename, line $linenumber)")
                continue
            end
            push!((type == "Cint" ? tail : code),
                  "const $name = $(type)($value)")
        end
    end
    push!(code, "", "# Status codes.")
    append!(code, tail)
end

# Make sure to use a slash for directory separator.
const PATH_SEPARATOR_RE = Sys.iswindows() ? r"[/\\]+" : r"/+"
fix_path(str::AbstractString) = replace(str, PATH_SEPARATOR_RE => "/")

function make_deps(target::String;
                   header::String = "",
                   library::String = "",
                   compile::Bool =  Sys.islinux())
    AT_DIR = fix_path(get(ENV, "AT_DIR",
                          (Sys.iswindows() ? "C:/Program Files/AndorSDK3" :
                           "/usr/local/andor")))

    # Find Andor SDK header <atcore.h>
    if header != ""
        header = fix_path(header)
        incdir = fix_path(dirname(header))
    else
        incdir = ""
        hdrname = "atcore.h"
        if haskey(ENV, "AT_INCDIR")
            dir = fix_path(ENV["AT_INCDIR"])
            if !isfile(dir*"/"*hdrname)
                error("\nDirectory specified by environment variable ",
                      "\"AT_INCDIR\" (\"$dir\") does not contain file ",
                      "Andor SDK header file \"$hdrname\".\n",
                      "Fix the definition of \"AT_INCDIR\" and rebuild.")
            end
            incdir = dir
        else
            for dir in (AT_DIR*"/include", AT_DIR)
                if isdir(dir) && isfile(dir*"/"*hdrname)
                    incdir = dir
                    break
                end
            end
            if incdir == ""
                error("\nAndor SDK header file \"$hdrname\" not found.\n",
                      "Define environment variable \"AT_INCDIR\" with the ",
                      "directory containing this file and rebuild.")
            end
        end
        header = incdir*"/"*hdrname
    end

    # Find Andor SDK library.
    if library != ""
        library = fix_path(library)
        incdir = fix_path(dirname(library))
    else
        libdir = ""
        libname = (Sys.iswindows() ? "atcore."*Libdl.dlext :
                   "libatcore."*Libdl.dlext)
        libdir = ""
        if haskey(ENV, "AT_LIBDIR")
            dir = fix_path(ENV["AT_LIBDIR"])
            if !isfile(dir*"/"*libname)
                error("\nDirectory specified by environment variable ",
                      "\"AT_LIBDIR\" (\"$dir\") does not contain file ",
                      "Andor SDK library file \"$libname\".\n",
                      "Fix the definition of \"AT_LIBDIR\" and rebuild.")
            end
            libdir = dir
        else
            for dir in (AT_DIR*"/lib", AT_DIR)
                if isdir(dir) && isfile(dir*"/"*libname)
                    libdir = dir
                    break
                end
            end
            if libdir == ""
                error("\nAndor SDK library file \"$libname\" not found.\n",
                      "Define environment variable \"AT_LIBDIR\" with ",
                      "the directory containing this file and rebuild.")
            end
        end
        library = libdir*"/"*libname
    end

    if compile
        # FIXME: Building all require old usb.h to be available.
        #target = (Sys.islinux() ? "all" : "deps.jl")
        print(stderr, "\nI will run the following command:\n",
              "    make TARGET=$(repr(target)) AT_INCDIR=$(repr(incdir)) ",
              "AT_LIBDIR=$(repr(libdir)) AT_DLL=$(repr(library))\n\n")
        run(`make TARGET=$(repr(target)) AT_INCDIR=$(repr(incdir)) AT_LIBDIR=$(repr(libdir)) AT_DLL=$(repr(library))`)
    else
        code = parse_header(header)
        open(target, "w") do io
            println(io, "#")
            println(io, "# deps.jl --")
            println(io, "#")
            println(io, "# Definitions of types and constants for interfacing Andor cameras in Julia.")
            println(io, "#")
            println(io, "# *DO NOT EDIT* as this file has been automatically generated.")
            println(io, "#")
            println(io)
            println(io, "# Path to the dynamic library.")
            println(io, "const _DLL = $(repr(library))")
            println(io)
            foreach(line -> println(io, line), code)
        end
    end
end

end # module


text:
#=
  default run
=#
using HetaSimulator, Plots

model = load_jlmodel("./model.jl")

### default simulations
sim(model, tspan = (0,100)) |> plot

max_stars_repo_name: ivborissov/alco

max_issues_repo_name: ivborissov/alco

max_forks_repo_name: ivborissov/alco


text:
using SafeTestsets

@safetestset "FFT Tests" begin
    using FFTDists, StatsBase, Statistics, Distributions, Test, Random
    import StatsBase: normalize

    function trueXY(edges, data)
        h = fit(Histogram, data, edges)
        h = normalize(h)
        y = h.weights
        Δ = step(edges)*.5
        x = edges[1:end-1] .+ Δ
        return x,y
    end

    meanAbs(v1, v2) = mean(abs.(v1-v2))
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    Random.seed!(45484)
    model = LogNormal(-1,.7) + LogNormal(-1,.7)
    convolve!(model)
    edges = 0:.01:2
    trueVals = rand(model,10^6)
    x,y = trueXY(edges, trueVals)
    dens = pdf.(model, x)
    @test meanAbs(y, dens) ≈ 0 atol= .05
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    Random.seed!(4544)
    model = Gamma(5,.05) + LogNormal(-1,.7)
    convolve!(model)
    edges = 0:.01:2
    trueVals = rand(model, 10^6)
    x,y = trueXY(edges, trueVals)
    dens = pdf.(model, x)
    @test meanAbs(y, dens) ≈ 0 atol= .05
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    Random.seed!(4547)
    model = Uniform(0,1) + InverseGaussian(.5,.5)
    convolve!(model)
    edges = 0:.05:2
    trueVals = rand(model, 10^6)
    x,y = trueXY(edges, trueVals)
    dens = pdf.(model, x)
    @test meanAbs(y, dens) ≈ 0 atol= .05
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    Random.seed!(555)
    model = Uniform(0,.5) + Uniform(0,.5)
    convolve!(model)
    edges = 0:.01:2
    trueVals = rand(model, 10^6)
    x,y = trueXY(edges, trueVals)
    dens = pdf.(model, x)
    @test meanAbs(y, dens) ≈ 0 atol= .05
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    Random.seed!(985)
    model = Uniform(.5,1.5) + Uniform(0,.1^6)
    convolve!(model)
    edges = 0:.01:2
    trueVals = rand(model, 10^6)
    x,y = trueXY(edges, trueVals)
    dens = pdf.(model, x)
    @test meanAbs(y, dens) ≈ 0 atol= .05
    ############################################################
    #####################  Define Model  #######################
    ############################################################
    μ,σ = convolve_normal(motor=(μ=.21,N=1), cr=(μ=.05,N=2))
    x = rand(Uniform(.21 - .21/3, .21 + .21/3), 10^5) + rand(Uniform(.05 - .05/3, .05 + .05/3), 10^5) +
        rand(Uniform(.05 - .05/3, .05 + .05/3), 10^5)
    @test μ ≈ .21 + .05*2
    @test σ ≈ sqrt((.21/3*sqrt(1/3))^2 + 2*(.05/3*sqrt(1/3))^2)
    @test σ ≈ std(x) rtol = .01
end


text:
# Benchmark on pairwise distance evaluation

using Distances

macro bench_pairwise_dist(repeat, dist, x, y)
    quote
        println("bench ", typeof($dist))

        # warming up
        r1 = evaluate($dist, ($x)[:,1], ($y)[:,1])
        pairwise($dist, $x, $y)

        # timing

        t0 = @elapsed for k = 1 : $repeat
            nx = size($x, 2)
            ny = size($y, 2)
            r = Array(typeof(r1), (nx, ny))
            for j = 1 : ny
                for i = 1 : nx
                    r[i, j] = evaluate($dist, ($x)[:,i], ($y)[:,j])
                end
            end
        end
        @printf "    loop:      t = %9.6fs\n" (t0 / $repeat)

        t1 = @elapsed for k = 1 : $repeat
            r = pairwise($dist, $x, $y)
        end
        @printf "    pairwise:  t = %9.6fs  |  gain = %7.4fx\n" (t1 / $repeat) (t0 / t1)
        println()
    end
end

m = 100
nx = 200
ny = 250

x = rand(m, nx)
y = rand(m, ny)

w = rand(m)
Q = rand(m, m)
Q = Q' * Q

@bench_pairwise_dist 20 SqEuclidean() x y
@bench_pairwise_dist 20 Euclidean() x y
@bench_pairwise_dist 20 Cityblock() x y
@bench_pairwise_dist 20 Chebyshev() x y
@bench_pairwise_dist  5 Minkowski(3.0) x y
@bench_pairwise_dist 20 Hamming() x y

@bench_pairwise_dist 20 CosineDist() x y
@bench_pairwise_dist 10 CorrDist() x y
@bench_pairwise_dist 20 ChiSqDist() x y
@bench_pairwise_dist 10 KLDivergence() x y
@bench_pairwise_dist  5 JSDivergence() x y

@bench_pairwise_dist 10 BhattacharyyaDist() x y
@bench_pairwise_dist 10 HellingerDist() x y

@bench_pairwise_dist 20 WeightedSqEuclidean(w) x y
@bench_pairwise_dist 20 WeightedEuclidean(w) x y
@bench_pairwise_dist 20 WeightedCityblock(w) x y
@bench_pairwise_dist  5 WeightedMinkowski(w, 3.0) x y
@bench_pairwise_dist 20 WeightedHamming(w) x y

@bench_pairwise_dist 10 SqMahalanobis(Q) x y
@bench_pairwise_dist 10 Mahalanobis(Q) x y


text:
"""
    ChebyshevTransform(; modes)

Constructs a discrete Chebyshev transform transform with modes `modes` that
operates on the 2nd, 3rd, ...N-modes+1-st dimension of the data.
Input must be have even-numbered size.
# Example
```
julia> ChebyshevTransform(modes = (12, 3, 4, 12))
```
"""
Base.@kwdef struct ChebyshevTransform{N, T} <: AbstractTransform
    modes::NTuple{N, T}
end

function forward(::ChebyshevTransform{N}, x) where {N}
    return FFTW.r2r(x, FFTW.REDFT00, 2:(N + 1))
end

function inverse(::ChebyshevTransform{N}, x) where {N}
    return FFTW.r2r(
        x ./ (prod(2 .* (size(x)[2:(N + 1)] .- 1))),
        FFTW.REDFT00,
        2:(N + 1),
    )
end

function truncate_modes(tr::ChebyshevTransform{N}, coeff) where {N}
    # Returns a low-pass filtered version of coeff assuming
    # that coeff is a tensor of spectral weights.
    # Ex.: tr.modes = (3,)
    # [0, 1, 2, 3, 4, 5] -> [0, 1, 2]
    # [a, b, c, d, e, f] -> [a, b, c]

    # indices for the spectral coefficients that we need to retain
    inds = [collect(1:m) for (s, m) in zip(size(coeff)[2:(end - 1)], tr.modes)]

    coeff_truncated = OperatorFlux.mview(coeff, inds, Val(N))

    return coeff_truncated
end

function pad_modes(
    ::ChebyshevTransform,
    coeff,
    size_pad::NTuple{N, T},
) where {N, T}
    # return a padded-with-zeros version of coeff assuming
    # that coeff is a tensor of spectral weights, thereby inflating coeff.
    # generate a zero array and indices that point to the filled-in
    # locations
    size_padded = (size(coeff)[1], size_pad..., size(coeff)[end])
    coeff_padded = zeros(eltype(coeff), size_padded)
    inds = [collect(1:m) for (s, m) in zip(size_pad, size(coeff)[2:(end - 1)])]

    coeff_padded_view = OperatorFlux.mview(coeff_padded, inds, Val(N))
    coeff_padded_view .= coeff

    return coeff_padded
end

# Base extensions
Base.ndims(::ChebyshevTransform{N}) where {N} = N
Base.eltype(::ChebyshevTransform) = Float32
Base.size(tr::ChebyshevTransform) = [tr.modes...]

function Base.show(io::IO, tr::ChebyshevTransform)
    print(io, "ChebyshevTransform(modes = $(tr.modes), dims = $(tr.dims))")
end


text:
function starting_points!(pt0 :: Point{T}, fd :: Abstract_QM_FloatData{T}, id:: QM_IntData, itd :: IterData{T}) where {T<:Real}

    itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, :U), pt0.x)
    itd.ATy = mul!(itd.ATy, fd.AT, pt0.y)
    dual_val = itd.Qx .- itd.ATy .+ fd.c
    pt0.s_l = dual_val[id.ilow]
    pt0.s_u = -dual_val[id.iupp]

    # check distance to bounds δ for x, s_l and s_u
    itd.x_m_lvar .= @views pt0.x[id.ilow] .- fd.lvar[id.ilow]
    itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt0.x[id.iupp]
    if id.nlow == 0
        δx_l1, δs_l1 = zero(T), zero(T)
    else
        δx_l1 = max(-T(1.5)*minimum(itd.x_m_lvar), T(1.e-2))
        δs_l1 = max(-T(1.5)*minimum(pt0.s_l), T(1.e-4))
    end
    if id.nupp == 0
        δx_u1, δs_u1 = zero(T), zero(T)
    else
        δx_u1 = max(-T(1.5)*minimum(itd.uvar_m_x), T(1.e-2))
        δs_u1 = max(-T(1.5)*minimum(pt0.s_u), T(1.e-4))
    end

    # correct components that to not respect the bounds
    itd.x_m_lvar .+= δx_l1
    itd.uvar_m_x .+= δx_u1
    s0_l1 = pt0.s_l .+ δs_l1
    s0_u1 = pt0.s_u .+ δs_u1
    xs_l1, xs_u1 = dot(s0_l1, itd.x_m_lvar), dot(s0_u1, itd.uvar_m_x)
    if id.nlow == 0
        δx_l2, δs_l2 = zero(T), zero(T)
    else
        δx_l2 = δx_l1 + xs_l1 / sum(s0_l1) / 2
        δs_l2 = δs_l1 + xs_l1 / sum(itd.x_m_lvar) / 2
    end
    if id.nupp == 0
        δx_u2, δs_u2 = zero(T), zero(T)
    else
        δx_u2 = δx_u1 + xs_u1 / sum(s0_u1) / 2
        δs_u2 = δs_u1 + xs_u1 / sum(itd.uvar_m_x) / 2
    end
    δx = max(δx_l2, δx_u2)
    δs = max(δs_l2, δs_u2)
    pt0.x[id.ilow] .+= δx
    pt0.x[id.iupp] .-= δx
    pt0.s_l .= pt0.s_l .+ δs
    pt0.s_u .= pt0.s_u .+ δs

    # deal with the compensation phaenomenon in x if irng != []
    @inbounds @simd for i in id.irng
        if fd.lvar[i] >= pt0.x[i]
            pt0.x[i] = fd.lvar[i] + T(1e-4)
        end
        if pt0.x[i] >= fd.uvar[i]
            pt0.x[i] = fd.uvar[i] - T(1e-4)
        end
        if (fd.lvar[i] < pt0.x[i] < fd.uvar[i]) == false
            pt0.x[i] = (fd.lvar[i] + fd.uvar[i]) / 2
        end
    end

    # verify bounds
    @assert all(pt0.x .> fd.lvar) && all(pt0.x .< fd.uvar)
    @assert all(pt0.s_l .> zero(T)) && all(pt0.s_u .> zero(T))

    # update itd
    update_IterData!(itd, pt0, fd, id, false)

end


text:
"""
    Kmeans(K; max=1e+8, th=1e-4)
Kmeans method.

# Parameters:
- `K`: number of class
- `max`: maximum number of repitition
- `th`: converge threshold

# Example
```jldoctest
julia> model = Kmenas(3)
Kmeans{3}(Matrix{Float64}(undef, 0, 0), 100000000, 0.0001)

julia> using HorseML.Clustering: fit!

julia> fit!(model, x)

julia> model(x) |> size
(100, 3)
```
"""
mutable struct Kmeans{K}
    μ::Matrix{Float64}
    max::Int64
    th::Float64
    Kmeans(K; max = 1e+8, th = 1e-4) = new{K}(Array{Float64}(undef, 0, 0), max, th)
end

function clustering(μ, x, K)
    r = zeros(size(x, 1), K)
    for n in 1 : size(x, 1)
        diff = Matrix(undef, size(μ)...)
        for k in 1 : K
            diff[k, :] = x[n, :] - μ[k, :]
        end
        distance = @. abs(diff)
        r[n, argmin(vec(sum(distance, dims=2)))] = 1
    end
    return r
end

function update(μ, x, K)
    r = clustering(μ, x, K)
    μ = similar(μ)
    for k in 1 : K
        μ[k, :] = sum(x .* r[:, k], dims=1) / sum(r[:, k])
    end
    return μ
end

"""
    fit!(model, x)
fit Kmeans or GMM model with data.

size of `x` is (numer of data, number of features).
"""
function fit!(model::Kmeans{K}, x) where {K}
    coverge(x, th) = @. abs(x) < th
    ma = maximum(x, dims=1)
    mi = minimum(x, dims=1)
    interval = (ma-mi) / (K-1)
    μ = Array{Float64, 2}(undef, K, size(x, 2))
    for i in 1 : K-1
        μ[i, :] = mi
        mi += interval
    end
    for _ in 1 : model.max
        μ_new = update(μ, x, K)
        if coverge.(μ_new - μ, model.th) == trues(size(μ)...)
            model.μ = μ_new
            return
        end
        μ = μ_new
    end
    @warn "Not Converged!"
    model.μ = μ
end

(model::Kmeans{K})(x) where{K} = clustering(model.μ, x, K)


text:
using DifferentialEquations
using PyCall
using PyPlot

function aginglaw(v, theta, dc)
    return 1 - theta * v / dc
end

function sliplaw(v, theta, dc)
    return -v * theta / dc * log(v * theta / dc)
end

function calcdvthetaclassic!(du, u, p, t)
    dc, eta, sigman, a, b, mu, Vp, L, rho, statelaw = p
    theta = u[1]
    v = u[2]
    dtau1_dt = mu * (Vp - v) / L
    du[1] = statelaw(v, theta, dc)
    du[2] = 1 / (eta / sigman + a / v) * (dtau1_dt / sigman - b * du[1] / theta)
    return nothing
end

function calcdv2theta!(du, u, p, t)
    dc, eta, sigman, a, b1, b2, mu, Vp, L, rho, statelaw = p
    theta1 = u[1]
    theta2 = u[2]
    v = u[3]
    dtau1_dt = mu * (Vp - v) / L
    du[1] = statelaw(v, theta1, dc)
    du[2] = statelaw(v, theta2, dc)
    du[3] = 1 / (eta / sigman + a / v) * (dtau1_dt / sigman - b1 * du[1] / theta1 - b2 * du[2] / theta2)
    return nothing
end

function calcdvthetatwo!(du, u, p, t)
    dc, eta, sigman, a, b, mu, Vp, L, rho, statelaw = p

    theta1 = u[1]
    v1 = u[2]
    theta2 = u[3]
    v2 = u[4]
    theta3 = u[5]
    v3 = u[6]

    dtau1_dt = mu * (Vp - (v1 + 0.01*v2 + 0.04*v3)) / L
    dtau2_dt = mu * (Vp - (v2 + 0.05*v1 + 0.95*v3)) / L
    dtau3_dt = mu * (Vp - (v3 + 0.03*v1 - 0.95*v2)) / L

    du[1] = statelaw(v1, theta1, dc)
    du[2] = 1 / (eta / sigman + a / v1) * (dtau1_dt / sigman - b * du[1] / theta1)
    du[3] = statelaw(v2, theta2, dc)
    du[4] = 1 / (eta / sigman + a / v2) * (dtau2_dt / sigman - b * du[3] / theta2)
    du[5] = statelaw(v3, theta3, dc)
    du[6] = 1 / (eta / sigman + a / v3) * (dtau3_dt / sigman - b * du[5] / theta3)

    return nothing
end

function rs_two()
    close("all")

    # Model parameters
    siay = 365.25 * 24 * 60 * 60
    tspan = (0.0, siay * 20000.0) # This will run
    mu = 3e10
    nu = 0.25
    rho = 2700.0
    eta = mu / (2.0 * sqrt(mu / rho))
    L = 60 * 1e3
    a = 0.015
    b = 0.02
    vp = 1e-9
    sigman = 30e6
    dc = 0.1
    abstol = 1e-8
    reltol = 1e-8

    # Time integrate - classic
    icsclassic = [1e8; vp / 1000]
    pclassic = (dc, eta, sigman, a, b, mu, vp, L, rho, sliplaw)
    probclassic = ODEProblem(calcdvthetaclassic!, icsclassic, tspan, pclassic)
    solclassic = solve(probclassic, RK4(), abstol = abstol, reltol = reltol)

    # Plot velocity time series
    t1 = [x / siay for x in solclassic.t]
    theta1 = [x[1] for x in solclassic.u]
    v1 = [x[2] for x in solclassic.u]

    figure()
    plot(t1, v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    # Time integrate - two state parameters
    icsclassic = [1e8; 2e8; vp / 1000]
    pclassic = (dc, eta, sigman, a, b, 0.99*b, mu, vp, L, rho, sliplaw)
    probclassic = ODEProblem(calcdv2theta!, icsclassic, tspan, pclassic)
    solclassic = solve(probclassic, RK4(), abstol = abstol, reltol = reltol)

    # Plot velocity time series
    t1 = [x / siay for x in solclassic.t]
    theta1 = [x[1] for x in solclassic.u]
    v1 = [x[3] for x in solclassic.u]

    figure()
    plot(t1, v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")
    return

    # Time integrate - multi
    icstwo = [1e8; vp / 1000; 1e8; vp / 1000; 1e8 ; vp / 1000]

    ptwo = (dc, eta, sigman, a, b, mu, vp, L, rho, sliplaw)
    probtwo = ODEProblem(calcdvthetatwo!, icstwo, tspan, ptwo)
    soltwo = solve(probtwo, RK4(), abstol = abstol, reltol = reltol)

    # Plot velocity time series
    t1 = [x / siay for x in soltwo.t]
    theta1 = [x[1] for x in soltwo.u]
    v1 = [x[2] for x in soltwo.u]
    theta2 = [x[3] for x in soltwo.u]
    v2 = [x[4] for x in soltwo.u]
    v3 = [x[6] for x in soltwo.u]

    figure()
    plot(t1, v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    figure()
    plot(t1, v2, "--r", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    figure()
    plot(v3, "-g", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    figure()
    subplot(2, 2, 1)
    plot(log10.(v1), log10.(v1), "-k", linewidth = 0.1)
    subplot(2, 2, 2)
    plot(log10.(v1), log10.(v3), "-k", linewidth = 0.1)
    subplot(2, 2, 3)
    plot(log10.(v2), log10.(v3), "-k", linewidth = 0.1)

    d1 = zeros(length(v1))
    d2 = zeros(length(v1))
    d3 = zeros(length(v1))
    for i in 2:length(v1)
        d1[i] = d1[i-1] + (v1[i-1]) / (t1[i] - t1[i-1])
        d2[i] = d2[i-1] + (v2[i-1]) / (t1[i] - t1[i-1])
        d3[i] = d3[i-1] + (v3[i-1]) / (t1[i] - t1[i-1])
    end
    figure()
    plot(t1, d1, "-b", linewidth = 0.5)
    plot(t1, d2, "-r", linewidth = 0.5)
    plot(t1, d3, "-g", linewidth = 0.5)

    xlabel("time (years)")
    ylabel("d (m)")

    return nothing
end
rs_two()


text:
import Dates: Date, DateTime

"""NumPy array protocol type string (typestr) format

A string providing the basic type of the homogenous array. The basic string format
consists of 3 parts: a character describing the byteorder of the data
(<: little-endian, >: big-endian, |: not-relevant), a character code giving the basic
type of the array, and an integer providing the number of bytes the type uses.

https://zarr.readthedocs.io/en/stable/spec/v2.html#data-type-encoding
"""

include("MaxLengthStrings.jl")
using .MaxLengthStrings: MaxLengthString

primitive type ASCIIChar <: AbstractChar 8 end
ASCIIChar(x::UInt8) = reinterpret(ASCIIChar, x)
ASCIIChar(x::Integer) = ASCIIChar(UInt8(x))
UInt8(x::ASCIIChar) = reinterpret(UInt8, x)
Base.codepoint(x::ASCIIChar) = UInt8(x)
Base.show(io::IO, x::ASCIIChar) = print(io, Char(x))
Base.zero(::Union{ASCIIChar,Type{ASCIIChar}}) = ASCIIChar(Base.zero(UInt8))

using Dates: Period, TimeType, Date, DateTime, Dates
import Base.==
struct DateTime64{P} <: TimeType
    i::Int64
end
Base.convert(::Type{Date},t::DateTime64{P}) where P = Date(1970)+P(t.i)
Base.convert(::Type{DateTime},t::DateTime64{P}) where P = DateTime(1970)+P(t.i)
Base.show(io::IO,t::DateTime64{P}) where P = print(io,"DateTime64[",P,"]: ",string(DateTime(t)))
Base.isless(x::DateTime64{P}, y::DateTime64{P}) where P = isless(x.i, y.i)
==(x::DateTime64{P}, y::DateTime64{P}) where P = x.i == y.i
strpairs = [Dates.Year => "Y", Dates.Month => "M", Dates.Week => "W", Dates.Day=>"D",
    Dates.Hour => "h", Dates.Minute => "m", Dates.Second=>"s", Dates.Millisecond =>"ms",
    Dates.Microsecond => "us", Dates.Nanosecond => "ns"]
const jlperiod = Dict{String,Any}()
const pdt64string = Dict{Any, String}()
for p in strpairs
    jlperiod[p[2]] = p[1]
    pdt64string[p[1]] = p[2]
end
Base.convert(::Type{DateTime64{P}}, t::Date) where P = DateTime64{P}(Dates.value(P(t-Date(1970))))
Base.convert(::Type{DateTime64{P}}, t::DateTime) where P = DateTime64{P}(Dates.value(P(t-DateTime(1970))))
Base.convert(::Type{DateTime64{P}}, t::DateTime64{Q}) where {P,Q} = DateTime64{P}(Dates.value(P(Q(t.i))))
Base.zero(t::Union{DateTime64, Type{<:DateTime64}}) = t(0)
# Base.promote_rule(::Type{<:DateTime64{<:Dates.DatePeriod}}, ::Type{Date}) = Date
# Base.promote_rule(::Type{<:DateTime64{<:Dates.DatePeriod}}, ::Type{DateTime}) = DateTime
# Base.promote_rule(::Type{<:DateTime64{<:Dates.TimePeriod}}, ::Type{Date}) = DateTime
# Base.promote_rule(::Type{<:DateTime64{<:Dates.TimePeriod}}, ::Type{DateTime}) = DateTime

typestr(t::Type) = string('<', 'V', sizeof(t))
typestr(t::Type{>:Missing}) = typestr(Base.nonmissingtype(t))
typestr(t::Type{Bool}) = string('<', 'b', sizeof(t))
typestr(t::Type{<:Signed}) = string('<', 'i', sizeof(t))
typestr(t::Type{<:Unsigned}) = string('<', 'u', sizeof(t))
typestr(t::Type{Complex{T}} where T<:AbstractFloat) = string('<', 'c', sizeof(t))
typestr(t::Type{<:AbstractFloat}) = string('<', 'f', sizeof(t))
typestr(::Type{MaxLengthString{N,UInt32}}) where N = string('<', 'U', N)
typestr(::Type{MaxLengthString{N,UInt8}}) where N = string('<', 'S', N)
typestr(::Type{<:Array}) = "|O"
typestr(::Type{<:DateTime64{P}}) where P = "<M8[$(pdt64string[P])]"

const typestr_regex = r"^([<|>])([tbiufcmMOSUV])(\d*)(\[\w+\])?$"
const typemap = Dict{Tuple{Char, Int}, DataType}(
    ('b', 1) => Bool,
    ('S', 1) => ASCIIChar,
    ('U', 1) => Char,
)
sizemapf(x::Type{<:Number}) = sizeof(x)
sizemapf(x::Type{<:Complex{T}}) where T = sizeof(T)
typecharf(::Type{<:Signed}) = 'i'
typecharf(::Type{<:Unsigned}) = 'u'
typecharf(::Type{<:AbstractFloat}) = 'f'
typecharf(::Type{<:Complex}) = 'c'
foreach([Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128,
  UInt8,UInt16,UInt32,UInt64,UInt128,
  Complex{Float16},Complex{Float32},Complex{Float64}]) do t
    typemap[(typecharf(t),sizemapf(t))] = t
end

function typestr(s::AbstractString, filterlist=nothing)
    m = match(typestr_regex, s)
    if m === nothing
        throw(ArgumentError("$s is not a valid numpy typestr"))
    else

        byteorder, typecode, typesize, typespec = m.captures
        if byteorder == ">"
            throw(ArgumentError("Big-endian data not yet supported"))
        end
        if typecode == "O"
            if filterlist === nothing
                throw(ArgumentError("Object array can only be parsed when an appropriate filter is defined"))
            end
            return Vector{sourcetype(first(filterlist))}
        end
        isempty(typesize) && throw((ArgumentError("$s is not a valid numpy typestr")))
        tc, ts = first(typecode), parse(Int, typesize)
        if (tc in ('U','S')) && ts > 1
          return MaxLengthString{ts,tc=='U' ? UInt32 : UInt8}
        end
        if tc == 'M' && ts == 8
            #We have a datetime64 value
            return DateTime64{jlperiod[String(typespec)[2:end-1]]}
        end
        # convert typecode to Char and typesize to Int
        typemap[(tc,ts)]
    end
end

"""Metadata configuration of the stored array

Each array requires essential configuration metadata to be stored, enabling correct
interpretation of the stored data. This metadata is encoded using JSON and stored as the
value of the “.zarray” key within an array store.

https://zarr.readthedocs.io/en/stable/spec/v2.html#metadata
"""
struct Metadata{T, N, C, F}
    zarr_format::Int
    shape::Ref{NTuple{N, Int}}
    chunks::NTuple{N, Int}
    dtype::String  # structured data types not yet supported
    compressor::C
    fill_value::Union{T, Nothing}
    order::Char
    filters::F  # not yet supported
    function Metadata{T2, N, C, F}(zarr_format, shape, chunks, dtype, compressor,fill_value, order, filters) where {T2,N,C,F}
        #We currently only support version
        zarr_format == 2 || throw(ArgumentError("Zarr.jl currently only support v2 of the protocol"))
        #Do some sanity checks to make sure we have a sane array
        any(<(0), shape) && throw(ArgumentError("Size must be positive"))
        any(<(1), chunks) && throw(ArgumentError("Chunk size must be >= 1 along each dimension"))
        order === 'C' || throw(ArgumentError("Currently only 'C' storage order is supported"))
        new{T2, N, C, F}(zarr_format, shape, chunks, dtype, compressor,fill_value, order, filters)
    end
end

#To make unit tests pass with ref shape
import Base.==
function ==(m1::Metadata, m2::Metadata)
  m1.zarr_format == m2.zarr_format &&
  m1.shape[] == m2.shape[] &&
  m1.chunks == m2.chunks &&
  m1.dtype == m2.dtype &&
  m1.compressor == m2.compressor &&
  m1.fill_value == m2.fill_value &&
  m1.order == m2.order &&
  m1.filters == m2.filters
end

"Construct Metadata based on your data"
function Metadata(A::AbstractArray{T, N}, chunks::NTuple{N, Int};
        zarr_format::Integer=2,
        compressor::C=BloscCompressor(),
        fill_value::Union{T, Nothing}=nothing,
        order::Char='C',
        filters::Nothing=nothing,
        fill_as_missing = fill_value !== nothing && deprec_fillvalue(),
    ) where {T, N, C}
    T2 = (fill_value === nothing || !fill_as_missing) ? T : Union{T,Missing}
    Metadata{T2, N, C, typeof(filters)}(
        zarr_format,
        size(A),
        chunks,
        typestr(eltype(A)),
        compressor,
        fill_value,
        order,
        filters
    )
end

Metadata(s::Union{AbstractString, IO},fill_as_missing) = Metadata(JSON.parse(s),fill_as_missing)

function deprec_fillvalue()
    @warn """In the future, fill values will not be interpreted as missing values by default.
    Please set the keyword argument `fill_as_missing` to a boolean accordingly. Setting to `true`
    for now, but in the future `false` will be the default.
    """
    true
end

"Construct Metadata from Dict"
function Metadata(d::AbstractDict, fill_as_missing)
    # create a Metadata struct from it

    compdict = d["compressor"]
    compressor = getCompressor(compdict)

    filters = getfilters(d)

    T = typestr(d["dtype"], filters)
    N = length(d["shape"])
    C = typeof(compressor)
    F = typeof(filters)

    fv = fill_value_decoding(d["fill_value"], T)

    TU = (fv === nothing || !fill_as_missing) ? T : Union{T,Missing}

    Metadata{TU, N, C, F}(
        d["zarr_format"],
        NTuple{N, Int}(d["shape"]) |> reverse,
        NTuple{N, Int}(d["chunks"]) |> reverse,
        d["dtype"],
        compressor,
        fv,
        first(d["order"]),
        filters,
    )
end

"Describes how to lower Metadata to JSON, used in json(::Metadata)"
function JSON.lower(md::Metadata)
    Dict{String, Any}(
        "zarr_format" => md.zarr_format,
        "shape" => md.shape[] |> reverse,
        "chunks" => md.chunks |> reverse,
        "dtype" => md.dtype,
        "compressor" => md.compressor,
        "fill_value" => fill_value_encoding(md.fill_value),
        "order" => md.order,
        "filters" => md.filters
    )
end

# Fill value encoding and decoding as described in
# https://zarr.readthedocs.io/en/stable/spec/v2.html#fill-value-encoding

fill_value_encoding(v) = v
fill_value_encoding(::Nothing)=nothing
function fill_value_encoding(v::AbstractFloat)
    if isnan(v)
        "NaN"
    elseif isinf(v)
        v>0 ? "Infinity" : "-Infinity"
    else
        v
    end
end

Base.eltype(::Metadata{T}) where T = T

# this correctly parses "NaN" and "Infinity"
fill_value_decoding(v::AbstractString, T::Type{<:Number}) = parse(T, v)
fill_value_decoding(v::Nothing, ::Any) = v
fill_value_decoding(v, T) = T(v)
fill_value_decoding(v, ::Type{ASCIIChar}) = v == "" ? nothing : v


text:
using MicrostructureNoise, Distributions

using Random, LinearAlgebra

Random.seed!(10)
η = 0.01

n = 1000
tt = 2*[0.0; sort(rand(n-1)); 1.0]
T = tt[end]
σ(t) = 1.0 + sin(t)
# ad hoc euler
x = [0.0; cumsum(σ.(tt[2:end]).*randn(n).*sqrt.(diff(tt)))]
y = x + sqrt(η)*randn(n+1)
prior = MicrostructureNoise.Prior(
N = 20,

α1 = 0.0,
β1 = 0.0,

αη = 0.3,
βη = 0.3,

Πα = LogNormal(1., 0.5),
μ0 = 0.0,
C0 = 5.0
)

α = 0.3
σα = 0.1
td, θs, ηs, αs, p = MicrostructureNoise.MCMC(prior, tt, y, α, σα, 1500)

@test_throws DimensionMismatch MicrostructureNoise.MCMC(prior, tt, y[2:end], α, σα, 10)

@testset "skipfirst" begin
    @test_logs (:info, "skip observation y[1] at t[1] (skipfirst == true)") MicrostructureNoise.MCMC(prior, tt, y, α, σα, 10; skipfirst = true)
    @test_throws DimensionMismatch MicrostructureNoise.MCMC(prior, tt, y[3:end], α, σα, 10; skipfirst = true)

    Random.seed!(10)
    p1 = MicrostructureNoise.MCMC(prior, tt, [y[2];y[2:end]], α, σα, 10; skipfirst = true)
    Random.seed!(10)
    p2 = MicrostructureNoise.MCMC(prior, tt, y[2:end], α, σα, 10; skipfirst = true)
    @test p1[1] == p2[1]
    @test p1[2][end] == p2[2][end]
end

posterior = MicrostructureNoise.posterior_volatility(td, θs)

# plot(MicrostructureNoise.piecewise(posterior.post_t[1:end], posterior.post_mean_root[1:end])...)
# plot!(tt,1+sin.(tt), ylim=(0,2))

# Test posterior mean versus truth

@test abs(mean(ηs)-η) < 0.2*η
@test norm(diff(td).*(mean(θs, dims=2) - σ.(td[1:end-1]).^2)) < 1


text:
struct Concrete <: Name end

# mandatory methods
format(::Type{Concrete}) = TabularData
task(::Type{Concrete}) = Regression
source(::Type{Concrete}) = "https://archive.ics.uci.edu/ml/datasets/BlogFeedback"

function downloadlink(::Type{Concrete})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/concrete/compressive/Concrete_Data.xls"
end

function preprocess(N::Type{Concrete})
    return path -> preprocess_csv(N, path, :train; istarget = 9, header = true, excel = true)
end

nattributes(::Type{Concrete}) = (8, )
ninstances(::Type{Concrete}) = (1030, 0, 0)

# optional methods
function checksum(::Type{Concrete})
    return "710076c66b9ca3f8050e7942f3dcbdbe04013534daeb0077ffd3079a52d8e0c4"
end

name(::Type{Concrete}) = "Concrete Compressive Strength Data Set"
author(::Type{Concrete}) = ["I-Cheng Yeh"]

function citation(::Type{Concrete})
    return """
    @article{yeh1998modeling,
        title={Modeling of strength of high-performance concrete using artificial neural networks},
        author={Yeh, I-C},
        journal={Cement and Concrete research},
        volume={28},
        number={12},
        pages={1797--1808},
        year={1998},
        publisher={Elsevier}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
using Graphs

function MaxLambda2inEV(e::Int,v::Int)
    # finds the graph with maximum algebraic connectivity given v vertices and e edges
    # Input: Integer e number of edges, integer v number of vertices
    # Output: The value of maximum lambda and the graph object

    @assert v>=1
    if e > v^v
        # Check the number of edges
        error("Too many edges")
    end
    possible_edges = collect(combinations(collect([1:v]),2))
    possible_graphs = combinations(possible_edges,e)
    tot = length(possible_graphs)
    iter = 1
    λ_max = typemin(Float64)
    max_g = SimpleAdjacencyList
    for i in possible_graphs

        if mod(iter,25000) == 0
            # print progress
            @printf("processing graph %d/%d\n",iter,tot)
        end
        g = simple_adjlist(v,is_directed=false)
        for j in i
            add_edge!(g, j[1], j[2])   # add edge
        end

        laplacian = Symmetric(laplacian_matrix(g))
        eigs = eigvals(laplacian,v-1:v)

        e_vals = eigfact(L2,2:2)
        λ = e_vals.values
        ν = e_vals.vectors
        if λ > λ_max
            max_g = g
            λ_max = λ
        end
        iter += 1
    end
    return (λ_max, max_g)
end

#
# num_v = 3
# num_e = 8
#
# (λ_max, max_g) = MaxLambda2inEV(num_e,num_v)
#
# println(max_g)


text: include("Proving.jl")


text:
"""
 Tests for https://arxiv.org/abs/1303.5103 with test problems as in the paper.
 DRI1, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
using DiffEqDevTools
#using DiffEqGPU
seed = 103473
function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
"""

@info "Diagonal noise"

u₀ = [0.1,0.1]
function f3!(du,u,p,t)
  du[1] = 3//2*u[1]
  du[2] = 3//2*u[2]
end
function g3!(du,u,p,t)
  du[1] = 1//10*u[1]
  du[2] = 1//10*u[2]
end
dts = 1 .//2 .^(3:-1:0)
tspan = (0.0,1.0)

h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and  == 1//100**exp(301//100*t) if h3(z) = z^2 )

prob = SDEProblem(f3!,g3!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h3(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(5e6)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,DRI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
println("DRI1:", sim.𝒪est[:weak_final])

numtraj = Int(8e6)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,DRI1NM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.4 # order is 2.91
println("DRI1NM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 3.05
println("RI1:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI3(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.77
println("RI3:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI5(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.45
println("RI5:", sim.𝒪est[:weak_final])

numtraj = Int(6e6)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI6(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.5 # order is 1.55
println("R6:", sim.𝒪est[:weak_final])

numtraj = Int(4e6)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-1) < 0.5 # order is 1.44
println("RDI1WM:", sim.𝒪est[:weak_final])

numtraj = Int(6e6)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI2WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.5 # order is 1.55
println("RDI2WM:", sim.𝒪est[:weak_final])

numtraj = Int(5e4)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI3WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.84
println("RDI3WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI4WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
println("RDI4WM:", sim.𝒪est[:weak_final])

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
let
    mapf = MAPF(CBSEnv.LowLevelEnv(), [1,2,3,4], [5,6,7,8])
    @test num_agents(mapf) == 4
    @test num_goals(mapf) == 4

    action_type(mapf)
    state_type(mapf)
    cost_type(mapf)
    get_starts(mapf)
    get_goals(mapf)
    i = 1
    env = mapf.env
    get_start(mapf, i)
    get_goal(mapf, i)
    get_start(mapf, env, i)
end


text:
#=

# 160: Unipolar degenerate drift-diffusion
([source code](SOURCE_URL))

See: C. Cancès, C. Chainais-Hillairet, J. Fuhrmann, and B. Gaudeul, "A numerical-analysis-focused comparison of several finite volume schemes for a unipolar degenerate drift-diffusion model" IMA Journal of Numerical Analysis, vol. 41, no. 1, pp. 271–314, 2021.

Available from [https://doi.org/10.1093/imanum/draa002](https://doi.org/10.1093/imanum/draa002),
the preprint is on [arxiv1907.11126](https://arxiv.org/abs/1907.11126).

The problem consists of a Poisson equation for the electrostatic potential $\phi$:

```math
-\nabla \varepsilon \nabla \phi = z(2c-1)
```
and a degenerate drift-diffusion equation of the transport of a charged species $c$:

```math
\partial_t u  - \nabla\cdot \left(\nabla c  + c \nabla (\phi - \log (1-c) )\right)
```

In particular, the paper, among others, investigates the "sedan" flux discretization which is able to handle the degeneracy coming from the $\log (1-c)$ term. The earliest reference to this scheme we found in the source code of the [SEDAN III](http://www-tcad.stanford.edu/tcad/programs/sedan3.html) semiconductor device simulator.
=#

module Example160_UnipolarDriftDiffusion1D

using Printf

using VoronoiFVM
using ExtendableGrids
using GridVisualize

mutable struct Data
    eps::Float64
    z::Float64
    ic::Int32
    iphi::Int32
    Data()=new()
end

function classflux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    ic=data.ic
    iphi=data.iphi
    f[iphi]=data.eps*(u[iphi,1]-u[iphi,2])
    bp,bm=fbernoulli_pm(u[iphi,1]-u[iphi,2])
    f[ic]=bm*u[ic,1]-bp*u[ic,2]
end

function storage!(f,u,node,data)
    ic=data.ic
    iphi=data.iphi
    f[iphi]=0
    f[ic]=u[ic]
end

function reaction!(f,u,node,data)
    ic=data.ic
    iphi=data.iphi
    f[iphi]=data.z*(1-2*u[ic])
    f[ic]=0
end

function sedanflux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    ic=data.ic
    iphi=data.iphi
    f[iphi]=data.eps*(u[iphi,1]-u[iphi,2])
    mu1=-log(1-u[ic,1])
    mu2=-log(1-u[ic,2])
    bp,bm=fbernoulli_pm(data.z*2*(u[iphi,1]-u[iphi,2])+(mu1-mu2))
    f[ic]=bm*u[ic,1]-bp*u[ic,2]
end

function main(;n=20,Plotter=nothing,dlcap=false,verbose=false,unknown_storage=:sparse,DiffEq=nothing)

    h=1.0/convert(Float64,n)
    grid=VoronoiFVM.Grid(collect(0:h:1))

    data=Data()
    data.eps=1.0e-3
    data.z=-1
    data.iphi=1
    data.ic=2

    ic=data.ic
    iphi=data.iphi

    physics=VoronoiFVM.Physics(data=data,
                               num_species=2,
                               flux=sedanflux!,
                               reaction=reaction!,
                               storage=storage!
                               )
    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])

    boundary_dirichlet!(sys,iphi,1,5.0)
    boundary_dirichlet!(sys,iphi,2,0.0)
    boundary_dirichlet!(sys,ic,2,0.5)

    inival=unknowns(sys)
    @views inival[iphi,:].=2
    @views inival[ic,:].=0.5

    if !dlcap
        ## Create solver control info for constant time step size
        tstep=1.0e-3
        control=VoronoiFVM.NewtonControl()
        control.verbose=verbose
        control.Δt_min=tstep
        control.Δt=tstep
        control.Δt_grow=1.2
        control.Δt_max=0.1
        control.Δu_opt=100
        control.damp_initial=0.5
        if isnothing(DiffEq)
            tsol=solve(inival,sys,[0.0,10],control=control)
        else # does not work yet...
            tsol=solve(DiffEq,inival,sys,[0.0,10],
                       initializealg=DiffEq.NoInit(),
                       dt=tstep)
        end
        vis=GridVisualizer(Plotter=Plotter,layout=(1,1),fast=true)
        for log10t=-4:0.01:0
            time=10^(log10t)
            sol=tsol(time)
            scalarplot!(vis[1,1],grid,sol[iphi,:],label="ϕ",title=@sprintf("time=%.3g",time),flimits=(0,5))
            scalarplot!(vis[1,1],grid,sol[ic,:],label="c",flimits=(0,5))
            reveal(vis)
        end
        return sum(tsol[end])

    else  # Calculate double layer capacitance
        U=unknowns(sys)
        control=VoronoiFVM.NewtonControl()
        control.damp_initial=1.0e-5
        delta=1.0e-4
        @views inival[iphi,:].=0
        @views inival[ic,:].=0.5
        sys.boundary_values[iphi,1]=0

        dphi=1.0e-1
        phimax=5
        delta=1.0e-4
        vplus=zeros(0)
        cdlplus=zeros(0)
        vminus=zeros(0)
        cdlminus=zeros(0)
        cdl=0
        vis=GridVisualizer(Plotter=Plotter,layout=(1,1),fast=true)
        for dir in [1,-1]
            phi=0.0
            while phi<phimax
                sys.boundary_values[iphi,1]=dir*phi
                solve!(U,inival,sys,control=control)
                inival.=U
                Q=integrate(sys,physics.reaction,U)
                sys.boundary_values[iphi,1]=dir*phi+delta
                solve!(U,inival,sys,control=control)
                inival.=U

                scalarplot!(vis[1,1],grid,U[iphi,:],label="ϕ",title=@sprintf("Δϕ=%.3g",phi),flimits=(-5,5))
                scalarplot!(vis[1,1],grid,U[ic,:],label="c",flimits=(0,5))
                reveal(vis)

                Qdelta=integrate(sys,physics.reaction,U)
                cdl=(Qdelta[iphi]-Q[iphi])/delta
                if dir==1
                    push!(vplus,dir*phi)
                    push!(cdlplus,cdl)
                else
                    push!(vminus,dir*phi)
                    push!(cdlminus,cdl)
                end
                phi+=dphi
            end
        end

        scalarplot!(vis[1,1],vplus,cdlplus,color=:green,clear=true,flimits=(0,0.05))
        scalarplot!(vis[1,1],vminus,cdlminus,color=:green,clear=false,show=true)
        return cdl
    end
end

function test()

    isapprox(main(unknown_storage=:sparse,dlcap=false),18.721369939561963,rtol=1.0e-5)&&
    isapprox(main(unknown_storage=:sparse,dlcap=true),0.010759276468375045,rtol=1.0e-5)&&
    isapprox(main(unknown_storage=:dense,dlcap=false),18.721369939561963,rtol=1.0e-5)&&
    isapprox(main(unknown_storage=:dense,dlcap=true),0.010759276468375045,rtol=1.0e-5)

end
end


text:
"""
    defaults for Individuals constructor
"""

day=86400.0
mon=365/12*day
OneMonth=[-0.5*mon,0.5*mon]

solver_default(prob) = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
#solver_default(prob) = solve(prob,Euler(),dt=day)
param_default = ( 𝑇=OneMonth , 🔄=(x->x), u0=[], u1=[], v0=[], v1=[])
rec_default = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[],
                        lon=Float64[], lat=Float64[], fid=Int[])
postprocess_default = (x->x)

"""
    struct Individuals{T}

- Data:           📌 (position),   🔴(record),           🆔 (ID)
- Functions:      🚄 (velocity),   ∫ (integration), 🔧(postprocessing)
- NamedTuples:    𝑃  (parameters), 𝐷 (diagnostics),      𝑀 (metadata)

The velocity function 🚄 typically computes velocity at an arbitrary position within the
chosen space-time domain (📌 to start) by interpolating gridded variables obtained from 𝑃.

Default keyword constructor example:

```
df=DataFrame( ID=[], x=[], y=[], z=[], t = [])
𝐼=Individuals{Float64,2}(📌=zeros(3,10),🆔=1:10,🔴=deepcopy(df))
𝐼=Individuals(📌=zeros(3,2),🆔=collect(1:2),🔴=deepcopy(df))
```

Plain text (or no-unicode) constructor example:

```
df=DataFrame( ID=[], x=[], y=[], z=[], t = [])
I=(position=zeros(3,2),ID=1:2,record=deepcopy(df))
I=Individuals(I)
```

Keyword cheatsheet:

- 📌=`\\:pushpin:<tab>`,          🔴=`\\:red_circle:<tab>`, 🆔=`\\:id:<tab>`
- 🚄=`\\:bullettrain_side:<tab>`, ∫=`\\int<tab>`,          🔧=`\\wrench<tab>`
- 𝑃=`\\itP<tab>`,                 𝐷=`\\itD<tab>`,           𝑀=`\\itM<tab>`
"""
Base.@kwdef struct Individuals{T,N}
   📌  ::Array{T,N} = Array{T,N}(undef, Tuple(Int.(zeros(1,N)))) #\:pushpin:<tab>
   🔴  ::DataFrame = similar(rec_default) #\:red_circle:<tab>
   🆔   ::Array{Int,1} = Array{Int,1}(undef, 0) #\:id:<tab>
   🚄  ::Function = dxy_dt #\:bullettrain_side:<tab>
   ∫   ::Function = solver_default #\int<tab>
   🔧  ::Function = postprocess_default #\wrench<tab>
   𝑃   ::NamedTuple = param_default #\itP<tab>
   𝐷   ::NamedTuple = NamedTuple() #\itD<tab>
   𝑀   ::NamedTuple = NamedTuple() #\itM<tab>vec
end

"""
    Individuals(NT::NamedTuple)

Constructor that uses a NamedTuple with only plain text keywords (i.e. no-unicode needed).

```
df=DataFrame( ID=[], x=[], y=[], z=[], t = [])
I=(position=zeros(3,2),ID=1:2,record=deepcopy(df))
I=Individuals(I)
```
"""
function Individuals(NT::NamedTuple)

    haskey(NT,:position) ? 📌=NT.position : 📌=Array{Float64,2}(undef, Tuple(Int.(zeros(1,2))))
    haskey(NT,:record) ? 🔴=NT.record : 🔴=similar(rec_default)
    haskey(NT,:ID) ? 🆔=NT.ID : 🆔=collect(1:size(📌,2))
    haskey(NT,:velocity) ? 🚄=NT.velocity : 🚄=dxy_dt
    haskey(NT,:integration) ? ∫=NT.integration : ∫=solver_default
    haskey(NT,:postprocessing) ? 🔧=NT.postprocessing : 🔧=postprocess_default
    haskey(NT,:parameters) ? 𝑃=NT.parameters : 𝑃=param_default
    haskey(NT,:diagnostics) ? 𝐷=NT.diagnostics : 𝐷=NamedTuple()
    haskey(NT,:metadata) ? 𝑀=NT.metadata : 𝑀=NamedTuple()
    isa(📌,UnitRange) ? 📌=collect(📌) : nothing
    haskey(NT,:type) ? T=NT.type : T=eltype(📌)

    Individuals{T,ndims(📌)}(📌=📌,🔴=🔴,🆔=🆔,🚄=🚄,∫=∫,🔧=🔧,𝑃=𝑃,𝐷=𝐷,𝑀=𝑀)
end

"""
    ∫!(𝐼::Individuals,𝑇::Tuple)

Displace simulated individuals continuously through space over time period 𝑇 starting from position 📌.

- This is typically achieved by computing the cumulative integral of velocity experienced by each individual along its trajectory (∫ 🚄 dt).
- The current default is `solve(prob,Euler(),dt=day)` but all solver options from the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package are available.
- After this, `∫!` is also equiped to postprocess results recorded into 🔴 via the 🔧 workflow, and the last step in `∫!` consiste in updating 📌 to be ready for continuing in a subsequent call to `∫!`.
"""
function ∫!(𝐼::Individuals,𝑇::Tuple)
    @unpack 🚄,📌,𝑃, 🔧, 🆔, 🔴, ∫ = 𝐼

    prob = ODEProblem(🚄,📌, 𝑇 ,𝑃)
    sol = ∫(prob)

    tmp = 🔧(sol,𝑃, id=🆔, 𝑇=𝑇)

    isempty(🔴) ? np =0 : np=length(🆔)
    append!(🔴,tmp[np+1:end,:])

    nd=length(size(sol))
    nd==3 ? 📌[:,:] = deepcopy(sol[:,:,end]) : 📌[:] = deepcopy(sol[:,end])

end

"""
    ∫!(𝐼::Individuals)

Call ∫!(𝐼::Individuals,𝐼.𝑃.𝑇)
"""
∫!(𝐼::Individuals) = ∫!(𝐼::Individuals,(𝐼.𝑃.𝑇[1],𝐼.𝑃.𝑇[2]))

## Convenience Methods (size,show,similar)

Base.size(A::Individuals) = size(A.📌)

function Base.show(io::IO, 𝐼::Individuals) where {T}
    @unpack 🚄,📌,𝑃, 𝐷, 𝑀, 🔧, 🆔, 🔴, ∫ = 𝐼
    printstyled(io, "  📌 details     = ",color=:normal)
    printstyled(io, "$(size(📌)) $(typeof(𝐼).parameters[1])\n",color=:blue)
    printstyled(io, "  🔴 details     = ",color=:normal)
    printstyled(io, "$(size(🔴)) $(names(🔴))\n",color=:blue)
    printstyled(io, "  🆔 range       = ",color=:normal)
    printstyled(io, "$(extrema(🆔))\n",color=:blue)
    printstyled(io, "  🚄 function    = ",color=:normal)
    printstyled(io, "$(🚄)\n",color=:blue)
    printstyled(io, "  ∫  function    = ",color=:normal)
    printstyled(io, "$(∫)\n",color=:blue)
    printstyled(io, "  🔧 function    = ",color=:normal)
    printstyled(io, "$(🔧)\n",color=:blue)
    printstyled(io, "  𝑃  details     = ",color=:normal)
    printstyled(io, "$(keys(𝑃))\n",color=:blue)
  return
end

function Base.similar(𝐼::Individuals)
    @unpack 🚄,📌,𝑃, 𝐷, 𝑀, 🔧, 🆔, 🔴, ∫ = 𝐼
    T = typeof(𝐼).parameters[1]
    return Individuals{T}(📌=similar(📌),🔴=similar(🔴),🆔=similar(🆔),
                          🚄=🚄, ∫=∫, 🔧=🔧, 𝑃=𝑃, 𝐷=𝐷, 𝑀=𝑀)
end

function Base.diff(𝐼::Individuals)
    f(x)=last(x).-first(x)
    🔴_by_ID = groupby(𝐼.🔴, :ID)
    return combine(🔴_by_ID,nrow,:lat => f => :dlat,:lon => f => :dlon)
end


text:
export projector

"""
    projector(x)

Return projector on `0` or projector on `1`.
"""
projector(x) = code==0 ? mat(P0) : mat(P1)


text:
# Partial correlation

"""
    partialcor(x, y, Z)

Compute the partial correlation of the vectors `x` and `y` given `Z`, which can be
a vector or matrix.
"""
function partialcor(x::AbstractVector, y::AbstractVector, Z::AbstractVecOrMat)
    length(x) == length(y) == size(Z, 1) ||
        throw(DimensionMismatch("Inputs must have the same number of observations"))
    length(x) > 0 || throw(ArgumentError("Inputs must be non-empty"))
    return Base.clampcor(_partialcor(x, mean(x), y, mean(y), Z))
end

function _partialcor(x::AbstractVector, μx, y::AbstractVector, μy, Z::AbstractMatrix)
    p = size(Z, 2)
    p == 1 && return _partialcor(x, μx, y, μy, vec(Z))
    z₀   = view(Z, :, 1)
    Zmz₀ = view(Z, :, 2:p)
    μz₀ = mean(z₀)
    rxz = _partialcor(x,  μx,  z₀, μz₀, Zmz₀)
    rzy = _partialcor(z₀, μz₀, y,  μy,  Zmz₀)
    rxy = _partialcor(x,  μx,  y,  μy,  Zmz₀)::typeof(rxz)
    return (rxy - rxz * rzy) / (sqrt(1 - rxz^2) * sqrt(1 - rzy^2))
end

function _partialcor(x::AbstractVector, μx, y::AbstractVector, μy, z::AbstractVector)
    μz = mean(z)

    # Initialize all of the accumulators to 0 of the appropriate types
    Σxx = abs2(zero(eltype(x)) - zero(μx))
    Σyy = abs2(zero(eltype(y)) - zero(μy))
    Σzz = abs2(zero(eltype(z)) - zero(μz))
    Σxy = zero(Σxx * Σyy)
    Σxz = zero(Σxx * Σzz)
    Σzy = zero(Σzz * Σyy)

    # We only want to make one pass over all of the arrays
    @inbounds begin
        @simd for i in eachindex(x, y, z)
            xi = x[i] - μx
            yi = y[i] - μy
            zi = z[i] - μz

            Σxx += abs2(xi)
            Σyy += abs2(yi)
            Σzz += abs2(zi)

            Σxy += xi * yi
            Σxz += xi * zi
            Σzy += zi * yi
        end
    end

    # Individual pairwise correlations
    rxy = Σxy / sqrt(Σxx * Σyy)
    rxz = Σxz / sqrt(Σxx * Σzz)
    rzy = Σzy / sqrt(Σzz * Σyy)

    return (rxy - rxz * rzy) / (sqrt(1 - rxz^2) * sqrt(1 - rzy^2))
end

"""
    PartialCorTest(x, y, Z)

Perform a t-test for the hypothesis that ``\\text{Cor}(x,y|Z=z) = 0``, i.e. the partial
correlation of vectors `x` and `y` given the matrix `Z` is zero.

Implements `pvalue` and `confint`. See also [`partialcor`](@ref).
"""
struct PartialCorTest <: HypothesisTest
    r::Real
    n::Int
    k::Int
    t::Real

    # Error checking is done in `partialcor`
    function PartialCorTest(x::AbstractVector, y::AbstractVector, Z::AbstractMatrix)
        r = partialcor(x, y, Z)
        n, k = size(Z)
        t = r * sqrt((n - 2 - k) / (1 - r^2))
        return new(r, n, k, t)
    end

    function PartialCorTest(x::AbstractVector, y::AbstractVector, z::AbstractVector)
        r = partialcor(x, y, z)
        n = length(z)
        t = r * sqrt((n - 3) / (1 - r^2))
        return new(r, n, 1, t)
    end
end

HypothesisTests.testname(::PartialCorTest) = "Test for partial correlation"
HypothesisTests.population_param_of_interest(p::PartialCorTest) =
    ("Partial correlation", zero(p.r), p.r)

StatsBase.nobs(p::PartialCorTest) = p.n
StatsBase.dof(p::PartialCorTest) = p.n - 2 - p.k

function StatsBase.confint(test::PartialCorTest, alpha::Float64=0.05)
    dof(test) > 1 || return (-one(test.r), one(test.r)) # Otherwise we can get NaNs
    q = quantile(Normal(), 1 - alpha / 2)
    fisher = (log1p(test.r) - log1p(-test.r)) / 2
    bound = q / sqrt(test.n - 3 - test.k)
    lo = 2 * (fisher - bound)
    hi = 2 * (fisher + bound)
    elo = Base.clampcor(expm1(lo) / (exp(lo) + 1))
    ehi = Base.clampcor(expm1(hi) / (exp(hi) + 1))
    return (elo, ehi)
end

HypothesisTests.default_tail(::PartialCorTest) = :both
HypothesisTests.pvalue(test::PartialCorTest; tail=:both) =
    pvalue(TDist(dof(test)), test.t, tail=tail)

function HypothesisTests.show_params(io::IO, test::PartialCorTest, indent="")
    println(io, indent, "number of observations:          ", nobs(test))
    println(io, indent, "number of conditional variables: ", test.k)
    println(io, indent, "t-statistic:                     ", test.t)
    println(io, indent, "degrees of freedom:              ", dof(test))
end


text:
using PipeLayout.GndStr

using LinearAlgebra
using JuMP
using SCIP
using GLPK

include("gndstruct_discdiam_iter.jl")
include("gndstruct_discdiam_minlp.jl")
include("gndstruct_discdiam_cb.jl")


text:
lattice = load_lattice(joinpath(@__DIR__, "CrFe2Ga_qe.in"))
atoms   = load_atoms(  joinpath(@__DIR__, "CrFe2Ga_qe.in"))
atoms   = attach_pseudos(atoms, Cr="hgh/pbe/cr-q14.hgh",
                                Fe="hgh/pbe/fe-q16.hgh",
                                Ga="hgh/pbe/ga-q3.hgh")
magnetic_moments = build_magnetic_moments(atoms; Ga=0.0, Fe=5.0, Cr=5.0)

cases = Dict{String, Any}(
    "adaptive" => (damping=DFTK.AdaptiveDamping(), ),  # Parameter-free adaptive damping
)
for α in 5:5:20
    cases["standard_$(@sprintf "%03d" α)"] = (damping=DFTK.FixedDamping(α / 100), )
end
for i in [1, 3]
    cases["adaptive_0$i"] = (damping=DFTK.AdaptiveDamping(i / 10), )
end

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 45,
    kgrid            = [13, 13, 13],
    supersampling    = 2.0,
    potential        = nothing,
    kwargs_scf       = (acceleration=DFTK.AndersonAcceleration(;m=0), ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
using Requires, Reexport

@reexport using ReachabilityAnalysis

# unexported methods
using ReachabilityAnalysis: _check_dim, _get_tspan, _get_cpost, _default_cpost,
                            ReachSolution, InitialValueProblem, numtype, post,
                            AbstractContinuousPost, TimeInterval,
                            AbstractLazyReachSet, AbstractTaylorModelReachSet

using Parameters: @with_kw

# aliases
const RA = ReachabilityAnalysis
const IA = IntervalArithmetic

using LazySets: _leq, _geq, isapproxzero, remove_zero_generators

import CommonSolve: solve

# optional dependencies
function __init__()
    @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
        include("init_DifferentialEquations.jl")
        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            import .Plots
            include("init_DifferentialEquations_Plots.jl")
        end
    end

    @require MAT = "23992714-dd62-5051-b70f-ba57cb901cac" begin
        using .MAT: matread
    end
end


text:
function isgrad{T <: Number}(
    g::Function,
    f::Function,
    x::T;
    tol::Real = @centralrule(x),
)
    return maxdiff(g(x), grad(f, x)) <= tol
end

function isgrad!{T <: Number}(
    g!::Function,
    f::Function,
    x::Vector{T};
    tol::Real = @centralrule(maximum(abs(x))),
)
    n = length(x)
    gr = Array(T, n)
    g!(gr, x)
    return maxdiff(gr, grad(f, x)) <= tol
end

function isgrad{T <: Number}(
    g::Function,
    f::Function,
    x::Vector{T};
    tol::Real =  @centralrule(maximum(abs(x))),
)
    return maxdiff(g(x), grad(f, x)) <= tol
end

function ishess{T <: Number}(
    h::Function,
    f::Function,
    x::T;
    tol::Real = @hessianrule(x),
)
    return maxdiff(h(x), hess(f, x)) <= tol
end

function ishess!{T <: Number}(
    h!::Function,
    f::Function,
    x::Vector{T};
    tol::Real = @hessianrule(maximum(abs(x))),
)
    n = length(x)
    H = Array(T, n, n)
    h!(H, x)
    return maxdiff(H, hess(f, x)) <= tol
end

function ishess{T <: Number}(
    h::Function,
    f::Function,
    x::Vector{T};
    tol::Real = @hessianrule(maximum(abs(x))),
)
    return maxdiff(h(x), hess(f, x)) <= tol
end


text:
@testset "cl1" begin
    data = open(readdlm, joinpath(@__DIR__, "data", "Cl1.txt"))

    for r in 1:size(data, 1)
        row      = data[r, :]
        x        = row[1]
        expected = row[2]

        @test ClausenFunctions.cl1(x) ≈ expected atol=1e-14
        @test ClausenFunctions.cl1(-x) ≈ expected atol=1e-14
        @test ClausenFunctions.cl1(x - 2.0*pi) ≈ expected atol=1e-12
        @test ClausenFunctions.cl1(x + 2.0*pi) ≈ expected atol=1e-12
    end

    @test ClausenFunctions.cl1(0.0) == Inf
end


text:
module problem112

description = """
Bouncy numbers

Working from left-to-right if no digit is exceeded by the digit to its left it is called an increasing number; for example, 134468.
Similarly if no digit is exceeded by the digit to its right it is called a decreasing number; for example, 66420.
We shall call a positive integer that is neither increasing nor decreasing a "bouncy" number; for example, 155349.
Clearly there cannot be any bouncy numbers below one-hundred, but just over half of the numbers below one-thousand (525) are bouncy. In fact, the least number for which the proportion of bouncy numbers first reaches 50% is 538.
Surprisingly, bouncy numbers become more and more common and by the time we reach 21780 the proportion of bouncy numbers is equal to 90%.

Find the least number for which the proportion of bouncy numbers is exactly 99%.
"""

function increasing(n :: Integer)
  lastDigit = 9
  while n > 0
    d = n % 10
    if d > lastDigit
      return false
    end
    lastDigit = d;
    n = div(n, 10)
  end
  return true
end

function decreasing(n :: Integer)
  lastDigit = 0
  while n > 0
    d = n % 10
    if d < lastDigit
      return false
    end
    lastDigit = d;
    n = div(n, 10)
  end
  return true
end

function bouncy(n :: Integer)
  return !(increasing(n) || decreasing(n))
end

function bouncySearch(percentage :: Integer)
  n = 1
  b = 0
  while b*100 != n*percentage
    n += 1
    b += bouncy(n) ? 1 : 0
  end
  return n
end

using Base.Test

@test bouncySearch(90) == 21780

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using Flux.Data
using Test

@test cmudict()["CATASTROPHE"] == :[K,AH0,T,AE1,S,T,R,AH0,F,IY0].args

@test length(CMUDict.phones()) == 39

@test length(CMUDict.symbols()) == 84

@test MNIST.images()[1] isa Matrix
@test MNIST.labels() isa Vector{Int64}

@test FashionMNIST.images()[1] isa Matrix
@test FashionMNIST.labels() isa Vector{Int64}

@test Data.Sentiment.train() isa Vector{Data.Tree{Any}}

@test Iris.features() isa Matrix
@test size(Iris.features()) == (4,150)

@test Iris.labels() isa Vector{String}
@test size(Iris.labels()) == (150,)


text:
# prologue.jl
using Yosys_jll
using Clang_jll
using Libdl
libcxxrtl = missing;

inc = joinpath(Yosys_jll.find_artifact_dir(), "share/yosys/include")
capi = joinpath(inc, "backends/cxxrtl/cxxrtl_capi.cc")
# end

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


text:
q = linspace(0.01,5,100)
x = 4*q.^2

g = (3,2,1)
α = (1,1,1)
C = (1,1,1)

@everywhere N = (25,25,50)
h = 0.000001

@everywhere include("junction-hamiltonian.jl")

######################### Just a Test #####################

@everywhere function energies(hw,I,phi,N)
    println("hw = $hw")
    println("I = $I")

    BLAS.set_num_threads(1)
    p = Parameters(hw,I,phi)
    @time Hm = FillHamiltonian(N...,p)
    @time solution = eigs(Hm,nev=4,tol=1e-18,maxiter=1000,which=:SR)
    return real(solution[1])
end

### Making a list of parameters
s = map(x->((sqrt(3/x),sqrt(2/x),sqrt(1/x)),(1/sqrt(2)/(3*x)^(1/4),1/sqrt(2)/(2*x)^(1/4),1/sqrt(2)/(1*x)^(1/4))),x)

@sync begin
     @async global E0 = pmap(s->energies(s...,(0,0,0),N),s)
     @async global Ex = pmap(s->energies(s...,(h,0,0),N),s)
     @async global Ey = pmap(s->energies(s...,(0,h,0),N),s)
     @async global Ez = pmap(s->energies(s...,(0,0,h),N),s)
end

DEx = (Ex - E0)/h
DEy = (Ey - E0)/h
DEz = (Ez - E0)/h

using JLD
save("outdir/derivatives-$N-$q-$g.jld","q",q,"DEx",DEx,"DEy",DEy,"DEz",DEz)

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type

# using Random, Distributions

"""
    LinSysVecRowPropToNormSampler{T} <: LinSysVecRowSampler where T <: Categorical

A parametric mutable structure that specifies sampling from the rows of the equation where
the probability of selecting a given equation is proportional to the sum of squares of the
coefficients of the given equation. The solver will appropriately initialize the
distribution.

See Strohmer, T., Vershynin, R. A Randomized Kaczmarz Algorithm with Exponential
Convergence. J Fourier Anal Appl 15, 262 (2009). https://doi.org/10.1007/s00041-008-9030-4

# Aliases
- `LinSysVecRowSVSampler`

# Fields
- `dist::T`, a categorical probability distribution.

Calling `LinSysVecRowPropToNormSampler()` or `LinSysVecRowSVSampler()` defaults `dist` to
`Categorical(1.0)`.
"""
mutable struct LinSysVecRowPropToNormSampler{T} <: LinSysVecRowSampler where T <:
    Categorical

    dist::T
end
LinSysVecRowSVSampler = LinSysVecRowPropToNormSampler
LinSysVecRowPropToNormSampler() = LinSysVecRowPropToNormSampler(Categorical(1.0))

# Common sample interface for linear systems
function sample(
    type::LinSysVecRowSVSampler{T} where T <: Categorical,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)
    if iter == 1
        prob = sum(A.^2, dims=2)[:,1]
        type.dist = Categorical(prob ./ sum(prob))
    end

    eqn_ind = rand(type.dist)

    return A[eqn_ind, :], b[eqn_ind]
end

#export LinSysVecRowPropToNormSampler, LinSysVecRowSVSampler