text:
"""
`SolutionPDAE`: Solution of a partitioned differential algebraic equation

Contains all fields necessary to store the solution of an PDAE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nm`: dimension of the constraint submanifold
* `nt`: number of time steps to store
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ni]` with `q[:,0,:]` the initial conditions
* `p`:  solution `p[nd, nt+1, ni]` with `p[:,0,:]` the initial conditions
* `λ`:  Lagrange multiplier `λ[nd, nt+1, ni]`
* `ntime`: number of time steps to compute
* `nsave`: store every nsave'th time step (default: 1)
* `nwrite`: save data to disk after every nwrite'th time step (default: ntime)
* `counter`: counter for copied solution entries
* `woffset`: counter for file offset
* `h5`: HDF5 file for storage
"""
abstract type SolutionPDAE{dType, tType, N} <: DeterministicSolution{dType, tType, N} end

# Create SolutionPDAEs with serial and parallel data structures.
for (TSolution, TDataSeries, Tdocstring) in
    ((:SSolutionPDAE, :SDataSeries, "Serial Solution of a partitioned differential algebraic equation."),
     (:PSolutionPDAE, :PDataSeries, "Parallel Solution of a partitioned differential algebraic equation."))
    @eval begin
        $Tdocstring
        mutable struct $TSolution{dType, tType, N} <: SolutionPDAE{dType, tType, N}
            nd::Int
            nm::Int
            nt::Int
            ni::Int
            t::TimeSeries{tType}
            q::$TDataSeries{dType,N}
            p::$TDataSeries{dType,N}
            λ::$TDataSeries{dType,N}
            ntime::Int
            nsave::Int
            nwrite::Int
            counter::Vector{Int}
            woffset::Int
            h5::HDF5File

            function $TSolution{dType, tType, N}(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nwrite) where {dType <: Number, tType <: Real, N}
                new(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nwrite, zeros(Int, ni), 0)
            end
        end

        function $TSolution(equation::Union{IODE{DT,TT},VODE{DT,TT},PDAE{DT,TT},IDAE{DT,TT}}, Δt::TT, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE; filename=nothing) where {DT,TT}
            @assert nsave > 0
            @assert ntime == 0 || ntime ≥ nsave
            @assert nwrite == 0 || nwrite ≥ nsave
            @assert mod(ntime, nsave) == 0

            if nwrite > 0
                @assert mod(nwrite, nsave) == 0
                @assert mod(ntime, nwrite) == 0
            end

            N  = equation.n > 1 ? 3 : 2
            nd = equation.d
            nm = equation.m
            ni = equation.n
            nt = div(ntime, nsave)
            nt = (nwrite == 0 ? nt : div(nwrite, nsave))
            nw = (nwrite == 0 ? ntime : nwrite)

            @assert nd > 0
            @assert nm > 0
            @assert ni > 0
            @assert nw > 0

            t = TimeSeries{TT}(nt, Δt, nsave)
            q = $TDataSeries(DT, nd, nt, ni)
            p = $TDataSeries(DT, nd, nt, ni)
            λ = $TDataSeries(DT, nm, nt, ni)
            s = $TSolution{DT,TT,N}(nd, nm, nt, ni, t, q, p, λ, ntime, nsave, nw)
            set_initial_conditions!(s, equation)

            if !isnothing(filename)
                isfile(filename) ? @warn("Overwriting existing HDF5 file.") : nothing
                create_hdf5(s, filename)
            end

            return s
        end

        function $TSolution(t::TimeSeries{TT}, q::$TDataSeries{DT,N}, p::$TDataSeries{DT,N}, λ::$TDataSeries{DT,N}, ntime::Int) where {DT,TT,N}
            @assert q.nd == p.nd >= λ.nd
            @assert q.nt == p.nt == λ.nt
            @assert q.ni == p.ni == λ.ni

            # extract parameters
            nd = q.nd
            nm = λ.nd
            ni = q.ni
            nt = t.n
            ns = div(ntime, nt)

            @assert mod(ntime, nt) == 0

            # create solution
            $TSolution{DT,TT,N}(nd, nm, nt, ni, t, q, p, λ, ntime, ns, 0)
        end

        function $TSolution(file::String)
            # open HDF5 file
            get_config(:verbosity) > 1 ? @info("Reading HDF5 file ", file) : nothing
            h5 = h5open(file, "r")

            # read attributes
            ntime = read(attrs(h5)["ntime"])
            nsave = read(attrs(h5)["nsave"])

            # reading data arrays
            t = TimeSeries(read(h5["t"]), nsave)
            q = $TDataSeries(read(h5["q"]))
            p = $TDataSeries(read(h5["p"]))
            λ = $TDataSeries(read(h5["λ"]))

            # need to close the file
            close(h5)

            # create solution
            $TSolution(t, q, p, λ, ntime)
        end
    end
end

hdf5(sol::SolutionPDAE)  = sol.h5
time(sol::SolutionPDAE)  = sol.t.t
ntime(sol::SolutionPDAE) = sol.ntime
nsave(sol::SolutionPDAE) = sol.nsave
offset(sol::SolutionPDAE) = sol.woffset

function set_initial_conditions!(sol::SolutionPDAE{DT,TT}, equ::Union{IODE{DT,TT},VODE{DT,TT},PDAE{DT,TT},IDAE{DT,TT}}) where {DT,TT}
    set_initial_conditions!(sol, equ.t₀, equ.q₀, equ.p₀, equ.λ₀)
end

function set_initial_conditions!(sol::SolutionPDAE{DT,TT}, t₀::TT, q₀::Union{Array{DT}, Array{TwicePrecision{DT}}}, p₀::Union{Array{DT}, Array{TwicePrecision{DT}}}, λ₀::Union{Array{DT}, Array{TwicePrecision{DT}}}) where {DT,TT}
    set_data!(sol.q, q₀, 0)
    set_data!(sol.p, p₀, 0)
    set_data!(sol.λ, λ₀, 0)
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function get_initial_conditions!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, k, n=1) where {DT,TT}
    get_solution!(sol, q, p, λ, n-1, k)
end

function get_initial_conditions!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, k, n=1) where {DT,TT}
    get_solution!(sol, q, p, n-1, k)
end

function get_initial_conditions(sol::SolutionPDAE, k, n=1)
    get_solution(sol, n-1, k)
end

function CommonFunctions.get_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, n, k) where {DT,TT}
    q .= sol.q[:, n, k]
    p .= sol.p[:, n, k]
    λ .= sol.λ[:, n, k]
end

function CommonFunctions.get_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, n, k) where {DT,TT}
    q .= sol.q[:, n, k]
    p .= sol.p[:, n, k]
end

function CommonFunctions.get_solution(sol::SolutionPDAE, n, k)
    (sol.t[n], sol.q[:, n, k], sol.p[:, n, k], sol.λ[:, n, k])
end

function CommonFunctions.set_solution!(sol::SolutionPDAE, t, q, p, n, k)
    set_solution!(sol, q, p, n, k)
end

function CommonFunctions.set_solution!(sol::SolutionPDAE, t, q, p, λ, n, k)
    set_solution!(sol, q, p, λ, n, k)
end

function CommonFunctions.set_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, λ::SolutionVector{DT}, n, k) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.p, p, sol.counter[k], k)
        set_data!(sol.λ, λ, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function CommonFunctions.set_solution!(sol::SolutionPDAE{DT,TT}, q::SolutionVector{DT}, p::SolutionVector{DT}, n, k) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.p, p, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function CommonFunctions.reset!(sol::SolutionPDAE)
    reset!(sol.q)
    reset!(sol.p)
    reset!(sol.λ)
    compute_timeseries!(sol.t, sol.t[end])
    sol.counter .= 1
    sol.woffset += sol.nt
end

function Base.close(solution::SolutionPDAE)
    # close(solution.t)
    # close(solution.q)
    # close(solution.p)
    # close(solution.λ)
    close(solution.h5)
end

"Creates HDF5 file and initialises datasets for PDAE solution object with single initial condition."
function create_hdf5(solution::SolutionPDAE{DT,TT,2}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution)

    # create datasets
    nt = div(solution.ntime, solution.nsave)
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, nt+1), "chunk", (solution.nd,1))
    p = d_create(solution.h5, "p", datatype(DT), dataspace(solution.nd, nt+1), "chunk", (solution.nd,1))
    λ = d_create(solution.h5, "λ", datatype(DT), dataspace(solution.nm, nt+1), "chunk", (solution.nm,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1] = solution.q[:, 0]
    p[:, 1] = solution.p[:, 0]
    λ[:, 1] = solution.λ[:, 0]

    return solution.h5
end

"Creates HDF5 file and initialises datasets for PDAE solution object with multiple initial conditions."
function create_hdf5(solution::SolutionPDAE{DT,TT,3}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution, solution.h5)

    # create datasets
    nt = div(solution.ntime, solution.nsave)
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ni), "chunk", (solution.nd,1,1))
    p = d_create(solution.h5, "p", datatype(DT), dataspace(solution.nd, nt+1, solution.ni), "chunk", (solution.nd,1,1))
    λ = d_create(solution.h5, "λ", datatype(DT), dataspace(solution.nm, nt+1, solution.ni), "chunk", (solution.nm,1,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1, :] = solution.q[:, 0, :]
    p[:, 1, :] = solution.p[:, 0, :]
    λ[:, 1, :] = solution.λ[:, 0, :]

    return solution.h5
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionPDAE{DT,TT,2}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # copy data from solution to HDF5 dataset
    h5["q"][:, j1:j2] = solution.q[:, 1:n]
    h5["p"][:, j1:j2] = solution.p[:, 1:n]
    h5["λ"][:, j1:j2] = solution.λ[:, 1:n]

    return nothing
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionPDAE{DT,TT,3}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # copy data from solution to HDF5 dataset
    h5["q"][:, j1:j2, :] = solution.q[:, 1:n, :]
    h5["p"][:, j1:j2, :] = solution.p[:, 1:n, :]
    h5["λ"][:, j1:j2, :] = solution.λ[:, 1:n, :]

    return nothing
end


text:
const VarCostArgs = Union{Float64, NTuple{2,Float64}, Vector{NTuple{2,Float64}}}

abstract type OperationalCost <: TechnicalParams end

mutable struct VariableCost{T}
    cost::T
end

get_cost(vc::PowerSystems.VariableCost) = vc.cost
Base.length(vc::PowerSystems.VariableCost) = length(vc.cost)
Base.getindex(vc::PowerSystems.VariableCost, ix::Int64) = getindex(vc.cost, ix)

ext: jl


text:
using TextParse

import TextParse: tryparsenext, unwrap, failedat, AbstractToken, LocalOpts
import CodecZlib: GzipCompressorStream
using Test
using Dates, Random
using Nullables

# dumb way to compare two AbstractTokens
Base.:(==)(a::T, b::T) where {T<:AbstractToken} = string(a) == string(b)

@testset "TextParse" begin

import TextParse: eatnewlines
@testset "eatnewlines" begin
    @test eatnewlines("\n\r\nx") == (4, 2)
    @test eatnewlines("x\n\r\nx") == (1, 0)

    # Also test the AbstractString variant
    @test eatnewlines(SubString("\n\r\nx", 1)) == (4, 2)
    @test eatnewlines(SubString("x\n\r\nx", 1)) == (1, 0)
end

import TextParse: eatwhitespaces
@testset "eatwhitespaces" begin
    @test eatwhitespaces("  x") == 3
    @test eatwhitespaces("x  x") == 1

    # Also test the AbstractString variant
    @test eatwhitespaces(SubString("  x", 1)) == 3
    @test eatwhitespaces(SubString("x  x", 1)) == 1
end

import TextParse: getlineend
@testset "getlineend" begin
    @test getlineend("\nx") == 0
    @test getlineend("x\nx") == 1
    @test getlineend("x\ny", 2) == 1
    @test getlineend("x\nyz", 3) == 4
    @test getlineend("xβ\nyz") == 2
end

import TextParse: getrowend
@testset "getrowend" begin
    opts = LocalOpts(',', false, '"', '"', true, true)

    @test getrowend("\nx", 1, lastindex("\nx"), opts, ',') == 0
    @test getrowend("x\nx", 1, lastindex("x\nx"), opts, ',') == 1
    @test getrowend("x\ny", 2, lastindex("x\ny"), opts, ',') == 1
    @test getrowend("x\nyz", 3, lastindex("x\nyz"), opts, ',') == 4
    @test getrowend("\"x\"\nyz", 1, lastindex("x\nyz"), opts, ',') == 3
    @test_throws ErrorException getrowend("\"x\" er", 1, lastindex("\"x\" er"), opts, ',')
    @test getrowend("", 1, lastindex(""), opts, ',') == 0
    @test getrowend("  ", 1, lastindex("  "), opts, ',') == 2
    @test getrowend("a,b,c\nd,e,f", 1, lastindex("a,b,c\nd,e,f"), opts, ',') == 5
    @test getrowend("a,\"b\"\"ef\",c\nd,e,f", 1, lastindex("a,\"b\"\"ef\",c\nd,e,f"), opts, ',') == 11
    @test getrowend("a,\"b\"\"ef\"", 1, lastindex("a,\"b\"\"ef\""), opts, ',') == 9
    @test_throws ErrorException getrowend("\"xy", 1, lastindex("\"xy"), opts, ',')

    opts = LocalOpts(',', false, '"', '\\', true, true)
    @test getrowend("a,\"bef\",f\na,b", 1, lastindex("a,\"bef\",f\na,b"), opts, ',') == 9
    @test_throws ErrorException getrowend("\"xy\\", 1, lastindex("\"xy"), opts, ',')
end

import TextParse: fromtype, Percentage
@testset "Float parsing" begin

    @test tryparsenext(fromtype(Float64), "1", 1, 1) |> unwrap == (1.0, 2)
    @test tryparsenext(fromtype(Float64), "12", 1, 2) |> unwrap == (12.0, 3)
    @test tryparsenext(fromtype(Float64), ".1", 1, 2) |> unwrap == (0.1, 3)
    @test tryparsenext(fromtype(Float64), "1.1", 1, 3) |> unwrap == (1.1, 4)
    @test tryparsenext(fromtype(Float32), "1.", 1, 2) |> unwrap == (1f0,3)
    @test tryparsenext(fromtype(Float64), "-1.1", 1, 4) |> unwrap == (-1.1,5)
    @test tryparsenext(fromtype(Float64), "-1.0e-12", 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(fromtype(Float64), "-1e-12") |> unwrap == (-1.0e-12,7)
    @test tryparsenext(fromtype(Float64), "-1.0E-12", 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(fromtype(Float64), "5.e-3", 1, 5) |> unwrap == (5.0e-3,6) # 32
    @test tryparsenext(Percentage(), "33%") |> unwrap == (.33,4)
    @test tryparsenext(Percentage(), "3.3%") |> unwrap == (.033,5)

    # Also test AbstractString variant
    @test tryparsenext(fromtype(Float64), SubString("1", 1), 1, 1) |> unwrap == (1.0, 2)
    @test tryparsenext(fromtype(Float64), SubString("12", 1), 1, 2) |> unwrap == (12.0, 3)
    @test tryparsenext(fromtype(Float64), SubString(".1", 1), 1, 2) |> unwrap == (0.1, 3)
    @test tryparsenext(fromtype(Float64), SubString("1.1", 1), 1, 3) |> unwrap == (1.1, 4)
    @test tryparsenext(fromtype(Float32), SubString("1.", 1), 1, 2) |> unwrap == (1f0,3)
    @test tryparsenext(fromtype(Float64), SubString("-1.1", 1), 1, 4) |> unwrap == (-1.1,5)
    @test tryparsenext(fromtype(Float64), SubString("-1.0e-12", 1), 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(fromtype(Float64), SubString("-1e-12", 1)) |> unwrap == (-1.0e-12,7)
    @test tryparsenext(fromtype(Float64), SubString("-1.0E-12", 1), 1, 8) |> unwrap == (-1.0e-12,9)
    @test tryparsenext(fromtype(Float64), SubString("5.e-3", 1), 1, 5) |> unwrap == (5.0e-3,6) # 32
    @test tryparsenext(Percentage(), SubString("33%", 1)) |> unwrap == (.33,4)
    @test tryparsenext(Percentage(), SubString("3.3%", 1)) |> unwrap == (.033,5)

    rng = MersenneTwister(0)
    floats = rand(1_000)
    parsed_floats = map(i->get(tryparsenext(fromtype(Float64), i, 1, lastindex(i))[1]), string.(floats))
    @test parsed_floats == floats

    # Also test AbstractString variant
    parsed_floats = map(i->get(tryparsenext(fromtype(Float64), SubString(i,1), 1, lastindex(i))[1]), string.(floats))
    @test parsed_floats == floats
end

@testset "Int parsing" begin
    @test tryparsenext(fromtype(Int64), "1", 1, 1) |> unwrap == (1, 2)
    @test tryparsenext(fromtype(Int64), "01", 1, 2) |> unwrap == (1, 3)
    @test tryparsenext(fromtype(Int64), "0001", 1, 4) |> unwrap == (1, 5)
    @test tryparsenext(fromtype(Int64), "123", 1, 3) |> unwrap == (123, 4)
    @test tryparsenext(fromtype(Int64), "00123", 1, 5) |> unwrap == (123, 6)
    @test tryparsenext(fromtype(Int64), "9223372036854775807", 1, 19) |> unwrap == (9223372036854775807, 20)
    @test tryparsenext(fromtype(Int64), "9223372036854775808", 1, 19) |> failedat == 1
    @test tryparsenext(fromtype(Int64), "19223372036854775808", 1, 20) |> failedat == 1

end

import TextParse: StringToken
using WeakRefStrings
@testset "String parsing" begin

    # default options
    @test tryparsenext(StringToken(String), "") |> unwrap == ("", 1)
    x = "x"
    @test tryparsenext(StringToken(WeakRefString), WeakRefString(pointer(x), 1)) |> unwrap == ("x", 2)
    @test tryparsenext(StringToken(String), "x") |> unwrap == ("x", 2)
    @test tryparsenext(StringToken(String), "x ") |> unwrap == ("x ", 3)
    @test tryparsenext(StringToken(String), " x") |> unwrap == (" x", 3)
    @test tryparsenext(StringToken(String), "x\ny") |> unwrap == ("x", 2)
    @test tryparsenext(StringToken(String), "x,y") |> unwrap == ("x", 2) # test escape

    opts = LocalOpts(',', false, '"', '"', true, true)
    @test tryparsenext(StringToken(String), "", opts) |> unwrap == ("", 1)
    @test tryparsenext(StringToken(String), "\"\"", opts) |> unwrap == ("\"\"", 3)
    @test tryparsenext(StringToken(String), "x", opts) |> unwrap == ("x", 2)
    # test including new lines
    @test tryparsenext(StringToken(String), "x\ny", opts) |> unwrap == ("x\ny", 4)
    @test tryparsenext(StringToken(String), "\"x\ny\"", opts) |> unwrap == ("\"x\ny\"", 6)

    opts = LocalOpts(',', false, '"', '"', false, true)
    # test that includequotes option doesn't affect string
    @test tryparsenext(StringToken(String), "\"\"", opts) |> unwrap == ("\"\"", 3)

    opts = LocalOpts(',', false, '"', '\\', false, false)
    str =  "Owner 2 ”Vicepresident\"\""
    @test tryparsenext(Quoted(String, '"', '\\'), str) |> unwrap == (str, lastindex(str)+1)
    str1 =  "\"Owner 2 ”Vicepresident\"\"\""
    @test tryparsenext(Quoted(String,'"', '"'), str1) |> unwrap == (str, lastindex(str1)+1)
    @test tryparsenext(Quoted(String,'"', '\\'), "\"\tx\"") |> unwrap == ("\tx", 5)
    opts = LocalOpts(',', true, '"', '\\', false, false)
    @test tryparsenext(StringToken(String), "x y",1,3, opts) |> unwrap == ("x", 2)

    @test tryparsenext(StringToken(String), "abcβ") |> unwrap == ("abcβ", 6)
end

import TextParse: Quoted, NAToken, Unknown
@testset "Quoted string parsing" begin
    opts = LocalOpts(',', false, '"', '"', true, true)

    @test tryparsenext(Quoted(String, '"', '"'), "\"\"") |> unwrap == ("", 3)
    @test tryparsenext(Quoted(String, '"', '"'), "\"\" ", opts) |> unwrap == ("", 3)
    @test tryparsenext(Quoted(String, '"', '"'), "\"x\"") |> unwrap == ("x", 4)
    @test tryparsenext(Quoted(String, '"', '"', includequotes=true), "\"x\"") |> unwrap == ("\"x\"", 4)
    str2 =  "\"\"\"\""
    @test tryparsenext(Quoted(String, '"', '"'), str2, opts) |> unwrap == ("\"\"", lastindex(str2)+1)
    str1 =  "\"x”y\"\"\""
    @test tryparsenext(Quoted(StringToken(String), '"', '"', required=true), "x\"y\"") |> failedat == 1

    @test tryparsenext(Quoted(String, '"', '"'), str1) |> unwrap == ("x”y\"\"", lastindex(str1)+1)
    @test tryparsenext(Quoted(StringToken(String), '"', '\\'), "\"x\\\"yz\"") |> unwrap == ("x\\\"yz", 8)
    @test tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "1") |> unwrap == (1,2)

    t = tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 1

    t = tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "\"\"") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 3
    @test tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "\"1\"") |> unwrap == (1, 4)

    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "\"abc\"") |> unwrap == ("abc", 6)
    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "x\"abc\"") |> unwrap == ("x\"abc\"", 7)
    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "\"a\nbc\"") |> unwrap == ("a\nbc", 7)
    @test tryparsenext(Quoted(StringToken(String), '"', '"', required=true), "x\"abc\"") |> failedat == 1
    @test tryparsenext(Quoted(fromtype(Int), '"', '"'), "21") |> unwrap == (21,3)
    @test tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "21") |> unwrap == (21,3)

    t = tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 1

    t = tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "\"\"") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 3

    @test tryparsenext(Quoted(NAToken(fromtype(Int)), '"', '"'), "\"21\"") |> unwrap == (21, 5)
    @test ismissing(tryparsenext(Quoted(NAToken(Unknown()), '"', '"'), " ") |> unwrap |> first)
    opts = LocalOpts(',', false,'"', '"', false, false)
    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "x,", opts) |> unwrap == ("x", 2)

    # stripspaces
    @test tryparsenext(Quoted(Percentage(), '"', '"'), "\" 10%\",", opts) |> unwrap == (0.1, 7)
    @test tryparsenext(Quoted(String, '"', '"'), "\" 10%\",", opts) |> unwrap == (" 10%", 7)
    opts = LocalOpts(',', true,'"', '"', false, false)
    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "\"x y\" y", opts) |> unwrap == ("x y", 6)
    @test tryparsenext(Quoted(StringToken(String), '"', '"'), "x y", opts) |> unwrap == ("x", 2)
end

@testset "NA parsing" begin
    t = tryparsenext(NAToken(fromtype(Float64)), ",") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 1

    t = tryparsenext(NAToken(fromtype(Float64)), "NA,") |> unwrap
    @test ismissing(t[1])
    @test t[2] == 3

    @test tryparsenext(NAToken(fromtype(Float64)), "X,") |> failedat == 1
    @test tryparsenext(NAToken(fromtype(Float64)), "1.212,") |> unwrap == (1.212, 6)
end

import TextParse: Field
@testset "Field parsing" begin
    f = fromtype(Int)
    @test tryparsenext(Field(f), "12,3") |> unwrap == (12, 4)
    @test tryparsenext(Field(f), "12 ,3") |> unwrap == (12, 5)
    @test tryparsenext(Field(f), " 12 ,3") |> unwrap == (12, 6)
    opts = LocalOpts('\t', false, 'x','x',true,false)
    @test tryparsenext(Field(f), "12\t3", 1, 4, opts) |> unwrap == (12, 4)
    @test tryparsenext(Field(f), "12 \t3", 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f), " 12 \t 3", 1, 6, opts) |> unwrap == (12, 6)
    opts = LocalOpts('\t', true, 'x','x',true,false)
    @test tryparsenext(Field(f), " 12 3", 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f, ignore_end_whitespace=false), " 12 \t 3", 1,6, opts) |> unwrap == (12, 5)
    opts = LocalOpts(' ', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), "12 3", 1,4,opts) |> unwrap == (12, 4)
#    @test tryparsenext(Field(f,ignore_end_whitespace=false), "12 \t3", 1,5,opts) |> failedat == 3
    opts = LocalOpts('\t', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), " 12\t 3", 1, 6, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), " 12\n", 1, 4, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), " 12\n\r\n", 1, 5, opts) |> unwrap == (12,6)
    @test tryparsenext(Field(f,eoldelim=true), " 12") |> unwrap == (12,4)

    # Also test AbstractString variant
    @test tryparsenext(Field(f), SubString("12,3",1)) |> unwrap == (12, 4)
    @test tryparsenext(Field(f), SubString("12 ,3",1)) |> unwrap == (12, 5)
    @test tryparsenext(Field(f), SubString(" 12 ,3",1)) |> unwrap == (12, 6)
    opts = LocalOpts('\t', false, 'x','x',true,false)
    @test tryparsenext(Field(f), SubString("12\t3",1), 1, 4, opts) |> unwrap == (12, 4)
    @test tryparsenext(Field(f), SubString("12 \t3",1), 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f), SubString(" 12 \t 3",1), 1, 6, opts) |> unwrap == (12, 6)
    opts = LocalOpts('\t', true, 'x','x',true,false)
    @test tryparsenext(Field(f), SubString(" 12 3",1), 1, 5, opts) |> unwrap == (12, 5)
    @test tryparsenext(Field(f, ignore_end_whitespace=false), SubString(" 12 \t 3",1), 1,6, opts) |> unwrap == (12, 5)
    opts = LocalOpts(' ', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), SubString("12 3",1), 1,4,opts) |> unwrap == (12, 4)
#    @test tryparsenext(Field(f,ignore_end_whitespace=false), "12 \t3", 1,5,opts) |> failedat == 3
    opts = LocalOpts('\t', false, 'x','x',false, false)
    @test tryparsenext(Field(f,ignore_end_whitespace=false), SubString(" 12\t 3",1), 1, 6, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), SubString(" 12\n",1), 1, 4, opts) |> unwrap == (12,5)
    @test tryparsenext(Field(f,eoldelim=true), SubString(" 12\n\r\n",1), 1, 5, opts) |> unwrap == (12,6)
    @test tryparsenext(Field(f,eoldelim=true), SubString(" 12",1)) |> unwrap == (12,4)
end

import TextParse: Record
@testset "Record parsing" begin
    r=Record((Field(fromtype(Int)), Field(fromtype(UInt)), Field(fromtype(Float64))))
    @test tryparsenext(r, "12,21,21,", 1, 9) |> unwrap == ((12, UInt(21), 21.0), 10)
    @test tryparsenext(r, "12,21.0,21,", 1, 9) |> failedat == 6
    s = "12   ,  21,  21.23,"
    @test tryparsenext(r, s, 1, length(s)) |> unwrap == ((12, UInt(21), 21.23), length(s)+1)
end

import TextParse: UseOne
@testset "UseOne" begin
    f = UseOne((Field(fromtype(Int)), Field(fromtype(Float64)), Field(fromtype(Int), eoldelim=true)), 3)
    @test tryparsenext(f, "1, 33.21, 45", 1, 12) |> unwrap == (45, 13)
end

import TextParse: Repeated
@testset "Repeated" begin
    f = Repeated(Field(fromtype(Int)), 3)
    @test tryparsenext(f, "1, 33, 45,", 1, 12) |> unwrap == ((1,33,45), 11)

    inp = join(map(string, [1:45;]), ", ") * ", "
    out = ntuple(identity, 45)
    f2 = Repeated(Field(fromtype(Int)), 45)
    @test tryparsenext(f2, inp, 1, length(inp)) |> unwrap == (out, length(inp))
    #@benchmark tryparsenext($f2, $inp, 1, length($inp))
end

import TextParse: quotedsplit
@testset "quotedsplit" begin
    opts = LocalOpts(',', false, '"', '\\', false, false)
    @test quotedsplit("x", opts, false, 1, 1) == ["x"]
    @test quotedsplit("x, y", opts, false, 1, 4) == ["x", "y"]
    @test quotedsplit("\"x\", \"y\"", opts,false, 1, 8) == ["x", "y"]
    @test quotedsplit("\"x\", \"y\"", opts,true, 1, 8) == ["\"x\"", "\"y\""]
    str = """x\nx,"s,", "\\",x" """
    @test quotedsplit(str, opts, false, 3, length(str)) == ["x", "s,", "\\\",x"]
    @test quotedsplit(",", opts, true, 1, 1) == ["", ""]
    @test quotedsplit(", ", opts, false, 1, 2) == ["", ""]
    str = "1, \"x \"\"y\"\" z\", 1"
    qopts = LocalOpts(',', false,'"', '"', false, false)
    @test quotedsplit(str, qopts,true, 1, lastindex(str)) == ["1", "\"x \"\"y\"\" z\"", "1"]
end

import TextParse: LocalOpts, readcolnames
@testset "CSV column names" begin
    str1 = """
     a, b,c d, e
    x,1,1,1
    ,1,1,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """

    str2 = """
     a, " b", "c", "d\\" e "
    """
    opts = LocalOpts(',', false, '"', '\\', false, false)
    @test readcolnames(str1, opts, 1, String[]) == (["a", "b", "c d", "e"], 13)
    @test readcolnames("\n\r$str1", opts, 3, Dict(3=>"x")) == (["a", "b", "x", "e"], 15)
    #@test readcolnames("$str2", opts, 3, Dict(3=>"x")) == (["a", "b", "x", "d\" e"], 24)
end

import TextParse: guesstoken, Unknown, Numeric, DateTimeToken, StrRange
@testset "guesstoken" begin
    opts = LocalOpts(UInt8(','), false, UInt8('"'), UInt8('"'), false, false)
    # Test null values
    @test guesstoken("", opts, Unknown()) == NAToken(Unknown())
    @test guesstoken("null", opts, Unknown()) == NAToken(Unknown())
    @test guesstoken("", opts, NAToken(Unknown())) == NAToken(Unknown())
    @test guesstoken("null", opts, NAToken(Unknown())) == NAToken(Unknown())

    # Test NA
    @test guesstoken("1", opts, NAToken(Unknown())) == NAToken(Numeric(Int))
    @test guesstoken("1", opts, NAToken(Numeric(Int))) == NAToken(Numeric(Int))
    @test guesstoken("", opts, NAToken(Numeric(Int))) == NAToken(Numeric(Int))
    @test guesstoken("1%", opts, NAToken(Unknown())) == NAToken(Percentage())

    # Test non-null numeric
    @test guesstoken("1", opts, Unknown()) == Numeric(Int)
    @test guesstoken("1", opts, Numeric(Int)) == Numeric(Int)
    @test guesstoken("", opts, Numeric(Int)) == NAToken(Numeric(Int))
    @test guesstoken("1.0", opts, Numeric(Int)) == Numeric(Float64)

    # Test strings
    @test guesstoken("x", opts, Unknown()) == StringToken(StrRange)

    # Test nullable to string
    @test guesstoken("x", opts, NAToken(Unknown())) == StringToken(StrRange)

    # Test string to non-null (short circuit)
    @test guesstoken("1", opts, StringToken(StrRange)) == StringToken(StrRange)

    # Test quoting
    @test guesstoken("\"1\"", opts, Unknown()) == Quoted(Numeric(Int), opts.quotechar, opts.escapechar)
    @test guesstoken("\"1\"", opts, Quoted(Numeric(Int), opts.quotechar, opts.escapechar)) == Quoted(Numeric(Int), opts.quotechar, opts.escapechar)

    # Test quoting with Nullable tokens
    @test guesstoken("\"\"", opts, Quoted(Unknown(), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Unknown()), opts.quotechar, opts.escapechar)
    @test guesstoken("\"\"", opts, Quoted(NAToken(Unknown()), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Unknown()), opts.quotechar, opts.escapechar)
    @test guesstoken("\"\"", opts, Quoted(Numeric(Int), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)
    @test guesstoken("\"\"", opts, Unknown()) == Quoted(NAToken(Unknown()), opts.quotechar, opts.escapechar)
    @test guesstoken("\"\"", opts, Numeric(Int)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)
    @test guesstoken("", opts, Quoted(Numeric(Int), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)
    @test guesstoken("", opts, Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)
    @test guesstoken("1", opts, Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)
    @test guesstoken("\"1\"", opts, Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)) == Quoted(NAToken(Numeric(Int)), opts.quotechar, opts.escapechar)

    # Test DateTime detection:
    tok = guesstoken("2016-01-01 10:10:10.10", opts, Unknown())
    @test tok == DateTimeToken(DateTime, dateformat"yyyy-mm-dd HH:MM:SS.s")
    @test guesstoken("2016-01-01 10:10:10.10", opts, tok) == tok
    @test guesstoken("2016-01-01 10:10:10.10", opts, Quoted(NAToken(Unknown()), opts.quotechar, opts.escapechar)) == Quoted(NAToken(tok), opts.quotechar, opts.escapechar)
end

import TextParse: guesscolparsers
@testset "CSV type detect" begin
    str1 = """
     a, b,c d, e
    x,1,1,1
    x,1,1,1
    x,1,1.,1
    x y,1.0,1,
    ,1.0,,1
    """
    opts = LocalOpts(',', false, '"', '\\', false, false)
    _, pos = readcolnames(str1, opts, 1, String[])
    testtill(i, colparsers=[]) = guesscolparsers(str1, String[], opts, pos, i, colparsers)
    @test testtill(0) |> first == Any[]
    @test testtill(1) |> first == map(fromtype, [StrRange, Int, Int, Int])
    @test testtill(2) |> first == map(fromtype, [StrRange, Int, Int, Int])
    @test testtill(3) |> first == map(fromtype, [StrRange, Int, Float64, Int])
    @test testtill(4) |> first == vcat(map(fromtype, [StrRange, Float64, Float64]),
                                       NAToken(fromtype(Int)))
    @test testtill(5) |> first == vcat(map(fromtype, [StrRange, Float64]),
                                       NAToken(fromtype(Float64)),
                                       NAToken(fromtype(Int)))
end

import TextParse: getlineat
@testset "getlineat" begin
    str = "abc\ndefg"
    @test str[getlineat(str,1)] == "abc\n"
    @test str[getlineat(str,4)] == "abc\n"
    @test str[getlineat(str,5)] == "defg"
    @test str[getlineat(str,lastindex(str))] == "defg"
    @test getlineat("x", 5) == 1:1
end

import TextParse: guessdateformat
@testset "date detection" begin
    @test guessdateformat("2016") |> typeof == DateTimeToken(Date, dateformat"yyyy-mm-dd") |> typeof
    @test guessdateformat("09/09/2016") |> typeof == DateTimeToken(Date, dateformat"mm/dd/yyyy") |> typeof
    @test guessdateformat("24/09/2016") |> typeof == DateTimeToken(Date, dateformat"dd/mm/yyyy") |> typeof
end

@testset "date parsing" begin
    tok = DateTimeToken(DateTime, dateformat"yyyy-mm-dd HH:MM:SS")
    opts = LocalOpts('y', false, '"', '\\', false, false)
    str = "1970-02-02 02:20:20"
    @test tryparsenext(tok, str, 1, length(str), opts) |> unwrap == (DateTime("1970-02-02T02:20:20"), length(str)+1)
    @test tryparsenext(tok, str*"x", 1, length(str)+1, opts) |> unwrap == (DateTime("1970-02-02T02:20:20"), length(str)+1)
    @test tryparsenext(tok, str[1:end-3]*"x", 1, length(str)-2, opts) |> failedat == length(str)-2
    @test tryparsenext(tok, str[1:end-3]*"y", 1, length(str)-2, opts) |> unwrap == (DateTime("1970-02-02T02:20"), length(str)-2)
end

import TextParse: _csvread
@testset "csvread" begin

    str1 = """
     a, b,c d, e
    x,1,1,1
    ,1,,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """
    data = ((["x", "","x","x y","x"],
              ones(5),
              [1,missing,1,1,missing],
              [1,1,1,missing,1]),
              ["a", "b", "c d", "e"])
    @test isequal(_csvread(str1, ','), data)
    coltype_test1 = _csvread(str1,
                            colparsers=Dict("b"=>Union{Missing, Float64},
                                          "e"=>Union{Missing,Float64}))
    coltype_test2 = _csvread(str1,
                            colparsers=Dict(2=>Union{Missing, Float64},
                                          4=>Union{Missing,Float64}))

    str2 = """
    x,1,1,1
    ,1,,1
    x,1,1.,1
    x y,1.0,1,
    x,1.0,,1
    """
    coltype_test3 = _csvread(str2, header_exists=false,
                            colparsers=Dict(2=>Union{Missing,Float64},
                                          4=>Union{Missing,Float64}))
    @test eltype(coltype_test1[1][2]) == Union{Missing, Float64}
    @test eltype(coltype_test1[1][4]) == Union{Missing, Float64}
    @test eltype(coltype_test2[1][2]) == Union{Missing, Float64}
    @test eltype(coltype_test2[1][4]) == Union{Missing, Float64}
    @test eltype(coltype_test3[1][2]) == Union{Missing, Float64}
    @test eltype(coltype_test3[1][4]) == Union{Missing, Float64}

    @test isequal(data, _csvread(str1, type_detect_rows=1))
    @test isequal(data, _csvread(str1, type_detect_rows=2))
    @test isequal(data, _csvread(str1, type_detect_rows=3))
    @test isequal(data, _csvread(str1, type_detect_rows=4))
    # But we can't go from a non-string column to a string column
    str3 = """
     a, b,c d, e
    1,1,1,1
    x,1,1,1
    """
    @test_throws TextParse.CSVParseError _csvread(str3, type_detect_rows=1)

    # test growing of columns if prediction is too low
    @test _csvread("x,y\nabcd, defg\n,\n,\n", type_detect_rows=1) ==
        ((String["abcd", "", ""], String["defg", "", ""]), String["x", "y"])

    # #19
    s="""
    x,y,z
    1,1,x
    "2",2,x
    1,2,"x \"\"y\"\""
    """

    res = (([1, 2, 1], [1, 2, 2], String["x", "x", "x \"\"y\"\""]), String["x", "y", "z"])
    @test _csvread(s, type_detect_rows=1, escapechar='"') == res
    @test _csvread(s, type_detect_rows=2, escapechar='"') == res

    @test csvread(IOBuffer("x\n1")) == (([1],),["x"])

    @test _csvread("x\n1\n") == (([1],),["x"])

    # test detection of newlines in fields
    s = """x, y
        abc, def
        g
        hi,jkl
        mno,pqr
        """

    @test _csvread(s, type_detect_rows=1) == ((["abc", "g\nhi", "mno"], ["def", "jkl", "pqr"]), ["x", "y"])
    # test custom na strings
    s = """
    x,y
    1,2
    ?,3
    4,*
    """
    nullness = ([false, true, false], [false, false, true])
    @test map(x->map(ismissing, x), first(_csvread(s, nastrings=["?","*"]))) == nullness
    @test map(x->map(ismissing, x), first(_csvread(s, nastrings=["?","*"], type_detect_rows=1))) == nullness

    @test isequal(csvread(["data/a.csv", "data/b.csv"]),
                  (([1.0, 2.0, 1.0, 2.0, 3.0], [2, 2, missing, missing, missing],
                    [missing, missing, missing, 2, 1]), String["x", "y", "z"], [2, 3]))
    @test isequal(csvread(["data/a.csv", "data/b.csv"], samecols=[("y","z")]),
                  (([1.0, 2.0, 1.0, 2.0, 3.0], [2, 2, missing, 2, 1]), String["x", "y"], [2,3]))

    # shouldn't fail because y doesn't exist
    @test _csvread("x\n1", colparsers=Dict("y"=>String)) == (([1],), ["x"])

    # Don't try to guess type if it's provided by user. Issue JuliaDB.jl#109
    s="""
    time,value
    "2017-11-09T07:00:07.391101180",0
    """
    @test _csvread(s) == ((String["2017-11-09T07:00:07.391101180"], [0]), String["time", "value"])
    @test _csvread(s, colparsers=Dict(:time=>String)) == ((String["2017-11-09T07:00:07.391101180"], [0]), String["time", "value"])

    @test _csvread("") == ((), String[])

    @test _csvread("""x""y"", z
                   a""b"", 1""") == ((["a\"\"b\"\""], [1]), ["x\"\"y\"\"", "z"])
end

@testset "skiplines_begin" begin
    str1 = """
    hello

    world
    x,y,z
    1,1,1
    """
    @test _csvread(str1, skiplines_begin=3) == (([1], [1], [1]), String["x", "y","z"])

    s = """
    x,y z
    a,b 1
    e  	3
    """
    @test _csvread(s, spacedelim=true) == ((["a,b", "e"],[1,3]), ["x,y","z"])
end

import TextParse: eatwhitespaces
@testset "custom parser" begin
    floatparser = Numeric(Float64)
    percentparser = CustomParser(Float64) do str, i, len, opts
        num, ii = tryparsenext(floatparser, str, i, len, opts)
        if num === nothing
            return num, ii
        else
            # parse away the % char
            ii = eatwhitespaces(str, ii, len)
            c, k = iterate(str, ii)
            if c != '%'
                return Nullable{Float64}(), ii # failed to parse %
            else
                return num, k # the point after %
            end
        end
    end

    @test tryparsenext(percentparser, "10%")  |> unwrap == (10.0, 4)
    @test tryparsenext(percentparser, "10.32 %") |> unwrap == (10.32, 8)
    @test tryparsenext(percentparser, "2k%") |> failedat ==  2
end

@testset "read gzipped files" begin
    fn   = joinpath(@__DIR__, "data", "a.csv")
    fngz = fn*".gz"
    open(fn, "r") do ior
        open(GzipCompressorStream, fngz, "w") do iow
            write(iow, ior)
        end
    end
    @test csvread(fn)   == csvread(fngz)
    @test csvread([fn]) == csvread([fngz])
    if isfile(fngz)
        rm(fngz)
    end
end

include("test_vectorbackedstrings.jl")

end


text:
#=
Common parameters for estimating the probability of for a diffusion with X(0) =
a satisfying X(T) ∈ (b, ∞) for the double well potential V(x) = (x²-1)².
=#

using LinearAlgebra
using Random
# adjust this path as neccessary
using BasicMD

function V(X)
    return (X⋅X - 1)^2;
end
function gradV!(gradV, X)
    @. gradV = 4.0 * (X^2 - 1) * X;
    gradV
end

a = -1.0;   # starting point
x₀ = [a];
b = 0.5;    # target set, (b, ∞)
Δt = 1e-4;  # time step
β = 5.0;    # inverse temperature
T = 1.0;    # terminal time
nΔt = ceil(Int, T/Δt);

f = x-> Int(x[1]>b); # define observable

# define sampler
sampler = MALA(V, gradV!, β, Δt);


text:
module NewtonsMethod

using LinearAlgebra, Statistics, Compat, ForwardDiff

function newtonroot(f, f_prime; x_0, tolerance = 1E-7, maxiter = 1000)
    x_old = x_0
    error = Inf
    iter = 1
    while error > tolerance && iter <= maxiter
        if f_prime(x_old) ≈ 0.0
            return (root = nothing, error = nothing, iter = nothing)
        end
        x_new = x_old - f(x_old)/f_prime(x_old)
        error = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter == maxiter+1
        return (root = nothing, error = nothing, iter = nothing)
    else
        return (root = x_old, error = error, iter = iter)
    end
end

# Applying auto-differentiation
D(f) = x -> ForwardDiff.derivative(f, x)

newtonroot(f; x_0, tolerance = 1E-7, maxiter = 1000) = newtonroot(f, D(f), x_0 = x_0, tolerance = tolerance, maxiter = maxiter)

# Export function
export newtonroot

end


text:
# Autogenerated wrapper script for Xorg_xkeyboard_config_jll for aarch64-linux-musl
using Xorg_xkbcomp_jll
JLLWrappers.@generate_wrapper_header("Xorg_xkeyboard_config")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_xkbcomp_jll)
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DifferentialEquations, ModelingToolkit, DiffEqParamEstim
using Turing, Distributions, StatsPlots
using CSV, DataFrames

include("mcmcODE.jl")

#=
1. Set up the ODE
2. Import and process the data
3. Transform the ODE output to match the data
4. Generate a loss function
5. Run the MCMC algorithm + Postprocessing
=#

###############################################################
                    # 1. Set up the ODE
###############################################################

#Constants in the model
cellVol = 3e-12 #Cell Volume (liters)
Na = 6.02e23 #Avagadro's number

#Defining non zero species in the model
m2c(molecule) = @. 1e9*molecule/(cellVol*Na) #Converts molecules to nM

#Initial Conditions
intConditions = (
cGAS = m2c(1e3),
DNA = m2c(1e3),
Sting = m2c(1e3),
cCAMP = 0.0,
IRF3 = m2c(1e4),
IFNbm = 0.0,
IFNb = 0.0,
STAT = 0.0,
SOCSm = 0.0,
IRF7m = 0.0,
TREX1m = 0.0,
IRF7 = 0.0,
TREX1 = 0.0)

#Algebraic Equations (mass balances)
cGAStot, Stingtot, IRF3tot = (
intConditions[:cGAS],
intConditions[:Sting],
intConditions[:IRF3])

## Define the ODE model
@parameters t k1f k1r k3f k3r k4f kcat5 Km5 k5r kcat6 Km6 kcat7 Km7 kcat8 Km8 k8f k9f k10f1 k10f2 k11f k12f k13f k6f kcat2 Km2 τ4 τ6 τ7 τ8 τ9 τ10 τ11 τ12 τ13
@variables cGAS(t) DNA(t) Sting(t) cGAMP(t) IRF3(t) IFNβm(t) IFNβ(t) STAT(t) SOCSm(t) IRF7m(t) TREX1m(t) IRF7(t) TREX1(t)
@derivatives D'~t

#Write the equations for the states
eqs= [
  D(cGAS) ~ -k1f*cGAS*DNA + k1r*(cGAStot - cGAS)
  D(DNA) ~ -k1f*cGAS*DNA + k1r*(cGAStot - cGAS) - kcat2*TREX1*DNA / (Km2 + DNA)
  D(Sting) ~ -k3f*cGAMP*Sting + k3r*(Stingtot - Sting)
  D(cGAMP) ~ k4f*(cGAStot - cGAS) - k3f*cGAMP*Sting + k3f*(Stingtot - Sting) - τ4*cGAMP
  D(IRF3) ~ -kcat5*IRF3*(Stingtot - Sting) / (Km5 +IRF3) + k5r*(IRF3tot - IRF3)
  D(IFNβm) ~ kcat6*(IRF3tot - IRF3) / (Km6 + (IRF3tot - IRF3)) + k6f*IRF7 - τ6*IFNβm
  D(IFNβ) ~ kcat7*IFNβm / (Km7 + IFNβm) - τ7*IFNβ
  D(STAT) ~ kcat8*IFNβ / (Km8 + IFNβ) * 1.0/(1.0+k8f*SOCSm) - τ8*STAT
  D(SOCSm) ~ k9f*STAT - τ9*SOCSm
  D(IRF7m) ~ k10f1*STAT + k10f2*IRF7 - τ10*IRF7m
  D(TREX1m) ~ k11f*STAT - τ11*TREX1m
  D(IRF7) ~ k12f*IRF7m - τ12*IRF7
  D(TREX1) ~ k13f*TREX1m - τ13*TREX1]

#Generate a function that evaluates the ODE
de = ODESystem(eqs)
f = ODEFunction(de,[cGAS, DNA, Sting, cGAMP, IRF3, IFNβm, IFNβ, STAT, SOCSm, IRF7m, TREX1m, IRF7, TREX1],[ k1f, k1r, k3f, k3r, k4f, kcat5, Km5, k5r, kcat6, Km6, kcat7, Km7, kcat8, Km8, k8f, k9f, k10f1, k10f2, k11f, k12f, k13f, k6f, kcat2, Km2, τ4, τ6, τ7, τ8, τ9, τ10, τ11, τ12, τ13])

#Contruct the ODE Problem
#Time to run the simulation
tspan = (0.0, 48.0)
#Initial guesses for the parameters
pNames = (:k1f, :k1r, :k3f, :k3r, :k4f, :kcat5, :Km5, :k5r, :kcat6, :Km6, :kcat7, :Km7, :kcat8, :Km8, :k8f, :k9f, :k10f1, :k10f2, :k11f, :k12f, :k13f, :k6f, :kcat2, :Km2, :τ4, :τ6, :τ7, :τ8, :τ9, :τ10, :τ11, :τ12, :τ13)

p = [2.6899, 4.8505, 0.0356, 7.487, 517.4056, 22328.3852, 11226.3682,0.9341,
         206.9446, 10305.461, 47639702.95,3.8474, 13.006, 78.2048, 0.0209,
         0.0059, 0.001, 0.0112, 0.001, 99.9466, 15.1436,0.0276, 237539.3249,
         61688.259, 0.96, 1.347, 12242.8736,1.2399, 1.5101, 0.347, 0.165, 6.9295,
         0.0178]

#Get the values for the initial condtions
u0 = collect(intConditions)

#Contruct the problem
prob = ODEProblem(f,u0,tspan,p)

#Choose an ODE algorithm and solve for the solution
alg = Rodas5()
sol = solve(prob,alg)

###############################################################
                    # 2. Import the data
###############################################################

#Read CSV in
data = CSV.read("NormData.csv",missingstring = "-")

#Two seperate error functions are used, one comparing normalized simulations and data, and the other comparing desired peaks (max/min)

#Define desired peaks
desiredScale = m2c([0.8*1e3,eps(),0.8*1e3, 1e5, 3e3,100.0,1e5,1e4,
                    30.0,30.0,30.0,1e4,1e4])

#Transform dataframe into Array of Arrays
timePoints = [data[data.Experiment .==i,2] for i ∈ unique(data.Experiment)]
experiments = [data[data.Experiment .==i,3:end] for i ∈ unique(data.Experiment)]

#Convert into matrices and vectors
timePoints = convert.(Vector,timePoints)
experiments = convert.(Matrix,experiments)

#Scale data by the desired value
dataIndex = [6,7,10,9,11,5,3]

scaledExperiments = [ desiredScale[dataIndex]'.*exp for exp in experiments]

###############################################################
        # 3. Transform the ODE output to match the data
###############################################################

#Depending on the data you might want to convert your ODE solution to match the data (convert to log fold change, measuring different species in mass balance, etc)

#Example expIdx
#expIdx = @. !ismissing(experiments[1][1,:])

function Observables(sol,intConditions,dataIndex,t)
  obversedSpecies = sol(t)[dataIndex,:]

  if isa(t,Number)
    #Have data for activated protiens, ode returns inactive
    obversedSpecies[end] = intConditions[:Sting] - obversedSpecies[end]
    obversedSpecies[end-1] = intConditions[:IRF3] - obversedSpecies[end-1]
  else
    obversedSpecies[end,:] = intConditions[:Sting] .- obversedSpecies[end,:]
    obversedSpecies[end-1,:] = intConditions[:IRF3] .- obversedSpecies[end-1,:]
  end

  #What species are being measured?
  return obversedSpecies
end

function heuristic(sol,desiredScale)
  peaksModel = maximum(sol[:,:],dims=2)
  peaksModel[[1,3,5]] = minimum(sol[[1,3,5],:],dims=2)
  h = @. (1.0-peaksModel/desiredScale)^2
  h[2] = 0.0
  return sum(h)
end
#What parameters are we changing?
#Need to pass the parameter values and the indicies of the parameters you want to change
parVary = (:k1f, :k1r, :k3f, :k3r, :k4f, :kcat5, :Km5, :k5r, :kcat6, :Km6, :k11f, :k13f, :kcat2, :Km2, :τ7, :kcat7, :Km7, :kcat8, :Km8, :k8f, :k9f, :k10f1, :k10f2, :k12f, :k6f)

#Number of parameters in the model to vary
parNum = length(parVary)

parVaryIdx = findall(par -> par ∈ parVary, pNames)

parChange = NamedTuple{parVary}(parVaryIdx)

###############################################################
                    # 4. Run the MCMC algorithm
###############################################################

#Provide any prior knowledge for the parameters
priors = fill(Uniform(0.0,1.0),parNum)

#include("mcmcODE.jl")
result = mcmcODE(prob,alg,timePoints,experiments,priors,parVary,parVaryIdx;num_samples=10_000)

###############################################################
                # 5. Postprocessing & Plotting
###############################################################

#Summary of the MCMC results
#describe(result)

#Save the parameter chains and information about the MCMC run
chainParameters = DataFrame(result,:parameters)
CSV.write("Parameters.csv",chainParameters)

chainInternals = DataFrame(result,:internals)
CSV.write("Internals.csv",chainInternals)

#Make a plot of the chains
chainsPlot = plot(result)
savefig(chainsPlot,"ChainsODE.pdf")

#Plot the data against the simulation

dataPlots = Vector(undef,length(dataIndex))
#loop through the data
for (i,protein) in enumerate(names(data)[3:end])
  dataPlots[i] = @df data plot(:Time,data[!,protein].*desiredScale[dataIndex[i]],
                 group=:Experiment,title=protein,markershape=:auto)
end

plot(dataPlots...,legend=false)

#Get the best parameter through the lowest log Posterior
bestPar = copy(p)
bestParIdx = argmax(chainInternals.lp)

#Why are the parameters out of order, why do you hate my Turing.jl?
for par in names(chainParameters)
  if par ∈ pNames
  bestPar[findfirst(x-> x==par,pNames)] = chainParameters[bestParIdx,par]
  end
end

newProb = remake(prob, p=bestPar)
newSol = solve(newProb,alg)

t = range(tspan[1],tspan[2],length=100)
obs = Observables(sol,intConditions,dataIndex,t)

plot!(t,obs',layout=length(dataIndex),legend=false,framestyle=:box)

max_stars_repo_name: RobertGregg/TuringODE

max_issues_repo_name: RobertGregg/TuringODE

max_forks_repo_name: RobertGregg/TuringODE


text:
using DSGE, JLD
using Base.Test

path = dirname(@__FILE__)

# Set up arguments
m = AnSchorfheide(testing = true)

system = jldopen("$path/../reference/forecast_args.jld","r") do file
    read(file, "system")
end

# Run impulse responses
states, obs, pseudo = impulse_responses(m, system)

# Compare to expected output
exp_states, exp_obs, exp_pseudo =
    jldopen("$path/../reference/impulse_responses_out.jld", "r") do file
        read(file, "exp_states"), read(file, "exp_obs"), read(file, "exp_pseudo")
    end

@testset "Compare irfs to expected output" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
end

nothing


text:
using SFML

ball_radius = 15.0
ballmaxvelocity = 8.0

paddle_width = 80.0
paddle_height = 5.0
paddle_velocity = 6.0

window_width = 800
window_height = 600

type Score
    paddle1::Int
    paddle2::Int
end

type Ball
    shape::CircleShape
    velocity::Vector2f
    starting_pos::Vector2f
end

function Ball(x, y)
    ball = Ball(CircleShape(), Vector2f(-ballmaxvelocity, -ballmaxvelocity), Vector2f(x, y))

    set_position(ball.shape, ball.starting_pos)
    set_radius(ball.shape, ball_radius)
    set_fillcolor(ball.shape, SFML.red)
    set_origin(ball.shape, Vector2f(ball_radius, ball_radius))

    return ball
end

function reset(ball::Ball)
    set_position(ball.shape, ball.starting_pos)
    ball.velocity = Vector2f(-ballmaxvelocity, -ballmaxvelocity)
end

function update(ball::Ball, score::Score)
    move(ball.shape, ball.velocity)

    position = get_position(ball.shape)
    x = position.x; y = position.y

    radius = get_radius(ball.shape)

    if x - radius < 0
        ball.velocity = Vector2f(ballmaxvelocity, ball.velocity.y)
    elseif x + radius > window_width
        ball.velocity = Vector2f(-ballmaxvelocity, ball.velocity.y)
    end

    if y - radius < 0
        score.paddle2 += 1
        # reset(ball)
        # sleep(1)
        ball.velocity = Vector2f(ball.velocity.x, ballmaxvelocity)
    elseif y + radius > window_height
        score.paddle1 += 1
        # reset(ball)
        # sleep(1)
        ball.velocity = Vector2f(ball.velocity.x, -ballmaxvelocity)
    end
end

type Paddle
    shape::RectangleShape
    velocity::Vector2f
    left::Int
    right::Int
end

function Paddle(x, y, left, right)
    paddle = Paddle(RectangleShape(), Vector2f(0.0, 0.0), left, right)

    set_position(paddle.shape, Vector2f(x, y))
    set_size(paddle.shape, Vector2f(paddle_width, paddle_height))
    set_fillcolor(paddle.shape, SFML.red)
    set_origin(paddle.shape, Vector2f(paddle_width / 2, paddle_height / 2))

    return paddle
end

function update(paddle::Paddle)
    move(paddle.shape, paddle.velocity)

    position = get_position(paddle.shape)
    x = position.x; y = position.y
    vy = paddle.velocity.y
    if (is_key_pressed(paddle.left))
        paddle.velocity = Vector2f(-paddle_velocity, vy)
    elseif (is_key_pressed(paddle.right))
        paddle.velocity = Vector2f(paddle_velocity, vy)
    else
        paddle.velocity = Vector2f(0, vy)
    end
end

function collides(ball::Ball, paddle::Paddle)
    ballbounds = get_globalbounds(ball.shape)
    paddlebounds = get_globalbounds(paddle.shape)

    return intersects(ballbounds, paddlebounds)
end

function main()
    score = Score(0, 0)

    score_text = RenderText()
    set_string(score_text, "$(score.paddle1)\n\n\n$(score.paddle2)")
    set_color(score_text, Color(255, 0, 0))
    set_position(score_text, Vector2f(60, window_height / 4))
    set_charactersize(score_text, 50)

    ball = Ball(window_width / 2, window_height / 2)
    paddle1 = Paddle(window_width / 2, window_height - 50, KeyCode.LEFT, KeyCode.RIGHT)
    paddle2 = Paddle(window_width / 2, 50, KeyCode.A, KeyCode.D)

    paddles = Paddle[]
    push!(paddles, paddle1)
    push!(paddles, paddle2)

    window = RenderWindow("Pong", window_width, window_height)
    set_framerate_limit(window, 60)
    set_vsync_enabled(window, true)

    event = Event()

    while isopen(window)
        while pollevent(window, event)
            if get_type(event) == EventType.CLOSED
                close(window)
            end
        end

        set_string(score_text, "$(score.paddle1)\n\n\n$(score.paddle2)")

        clear(window, SFML.white)

        if collides(ball, paddle1) || collides(ball, paddle2)
            ball.velocity = Vector2f(ball.velocity.x, -ball.velocity.y)
        end

        update(ball, score)
        draw(window, ball.shape)
        for i in 1:length(paddles)
            update(paddles[i])
            draw(window, paddles[i].shape)
        end

        draw(window, score_text)
        display(window)
    end
end

main()


text:
@inline function ExaTronKernels.dgpnorm(n::Int, x::CuDeviceArray{Float64,1}, xl::CuDeviceArray{Float64,1},
                         xu::CuDeviceArray{Float64,1}, g::CuDeviceArray{Float64,1})
    tx = threadIdx().x

    v = 0.0
    if tx <= n
        @inbounds begin
            if xl[tx] != xu[tx]
                if x[tx] == xl[tx]
                    v = min(g[tx], 0.0)
                elseif x[tx] == xu[tx]
                    v = max(g[tx], 0.0)
                else
                    v = g[tx]
                end

                v = abs(v)
            end
        end
    end

    # shfl_down_sync() will automatically sync threads in a warp.

    offset = 16
    while offset > 0
        v = max(v, CUDA.shfl_down_sync(0xffffffff, v, offset))
        offset >>= 1
    end
    v = CUDA.shfl_sync(0xffffffff, v, 1)

    return v
end


text:
fake_metsid(M) = ["M$i" for i in 1:M]
fake_rxnsid(N) = ["r$i" for i in 1:N]

# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)
struct MetNet{T<:Real}

    # Fields important for modeling
    S::AbstractMatrix{T} # Stoichiometric matrix M x N sparse
    b::AbstractVector{T} # right hand side of equation  S ν = b
    lb::AbstractVector{T} # fluxes lower bound N elements vector
    ub::AbstractVector{T} # fluxes upper bound N elements vector
    c::AbstractVector{T} # reaction index of biomass
    rxns::Vector{String} # reactions short-name N elements
    mets::Vector{String} # metabolites short-name M elements
    intake_info::Dict # Information required for enforcing the Chemostat bounds

    # Mustly for compatibility with COBRA
    metNames # metabolites long-names M elements
    rxnNames # reactions long-names N elements
    metFormulas # metabolites formula M elements
    genes # gene names
    rxnGeneMat #
    grRules # gene-reaction rule N elements vector of strings (and / or allowed)
    rev # reversibility of reactions N elements
    subSystems # cellular component of fluxes N elements

    function MetNet{T}(;kwargs...) where {T<:Real}
        kwargs = Dict(kwargs)

        # Just for better errors printing
        function checkget(k, T2)
            !haskey(kwargs, k) && error("Mandatoty field '$k' not found!!!")
            dat = kwargs[k]
            !(dat isa T2) && error("'$k' is a '$(typeof(dat))', expected '$T2'")
            return dat
        end

        # LP fields
        S = checkget(:S, AbstractMatrix{T})
        M, N = size(S)
        b = checkget(:b, AbstractVector{T})
        lb = checkget(:lb, AbstractVector{T})
        ub = checkget(:ub, AbstractVector{T})
        c = checkget(:c, AbstractVector{T})
        rxns = get(kwargs, :rxns, fake_rxnsid(N))
        mets = get(kwargs, :mets, fake_metsid(N))
        intake_info = get(kwargs, :intake_info, Dict())

        # Others
        metNames = get(kwargs, :metNames, [])
        rxnNames = get(kwargs, :rxnNames, [])
        metFormulas = get(kwargs, :metFormulas, [])
        genes = get(kwargs, :genes, [])
        rxnGeneMat = get(kwargs, :rxnGeneMat, [])
        grRules = get(kwargs, :grRules, [])
        rev = get(kwargs, :rev, [])
        subSystems = get(kwargs, :subSystems, [])

        new{T}(S, b, lb, ub, c, rxns, mets, intake_info,
            metNames, rxnNames, metFormulas, genes,
            rxnGeneMat, grRules, rev, subSystems)

    end
end

# Helpers
MetNet(;kwargs...) = MetNet{Float64}(;kwargs...)

# For COBRA .mat files use reshape
function MetNet(model_dict::Dict; reshape = false)
    net = to_symbol_dict(model_dict)
    reshape && (net = reshape_mat_dict(net))
    return MetNet(;net...)
end

"""
    Create a new MetNet from a template but overwriting the fields
    of the template with the given as kwargs.
    The returned MetNet will share the non-overwritten fields.
"""
function MetNet(template::MetNet; to_overwrite...)
    new_metnet_dict = Dict{Symbol, Any}(to_overwrite)

    for field in fieldnames(typeof(template))
        haskey(new_metnet_dict, field) && continue # avoid use the template version
        new_metnet_dict[field] = getfield(template, field)
    end

    return MetNet(new_metnet_dict; reshape = false)
end

# For compatibility with cobra mat models
function reshape_mat_dict(mat_dict::Dict; S::DataType = String,
    F::DataType = Float64, N::DataType = Int, I::DataType = Int64)

    new_dict::Dict = Dict()

    function _reshape!(rshf::Function, k::Symbol, T::Type)
        !haskey(mat_dict, k) && return
        dat = mat_dict[k]
        new_dict[k] = dat isa T ? deepcopy(dat) : rshf(dat)
    end

    # ----------------- Typed -----------------
    # String vectors
    for k in [:comps, :metNames, :metFormulas,
            :rxnFrom, :rxnNames, :genes, :inchis,
            :grRules, :mets, :rxns]
        _reshape!(k, Vector{S}) do dat
            dat |> vec .|> S
        end
    end

    # Float vectors
    for k in [:b, :lb, :ub, :c]
        _reshape!(k, Vector{F}) do dat
            dat |> vec .|> F
        end
    end

    # Integer vectors
    for k in [:metComps]
        _reshape!(k, Vector{I}) do dat
            floor.(I, dat |> vec)
        end
    end

    # S
    _reshape!(:S, Matrix{F}) do dat
        Matrix{F}(dat)
    end

    # ----------------- Type free -----------------
    # Matrices
    for k in [:rxnGeneMat]
        _reshape!(k, Matrix) do dat
            Matrix(dat)
        end
    end

    # Vectors
    for k in [:subSystems, :rev]
        _reshape!(k, Matrix) do dat
            dat |> vec |> collect
        end
    end

    return new_dict
end


text:
module Radon

export radon
export backprojection, backprojection_anim
export iradon, fitered_backprojection, SART

"""
    radon(f::Matrix{T}; θ=Vector(0:(π/180):π-π/180)::Vector{T}, in_circle=true) where {T<:Real}

return Radon transform of `f`.

# Keyword arguments

- `f`: 2-dimension array.(phantom image)
- `θ`: projection angle.(radian)
- `in_circle`:
"""
function radon(f::Matrix{T}; θ=Vector(0:(π/180):π-π/180)::Vector{T}, in_circle=true) where {T<:Real}
    ℛf = zeros(maximum(size(f)), length(θ))

    slim, θlim = size(ℛf)
    ylim, xlim = size(f)

    if in_circle == true
        L = ylim
    else
        L = round(Int64, hypot(xlim, ylim))
    end

    cos_ = cos.(θ)
    sin_ = sin.(θ)

    s = LinRange(-L/2, L/2, slim)
    t = LinRange(-L/2, L/2, L)

    for sᵢ ∈ 1:slim, θᵢ ∈ 1:θlim
        for l ∈ 1:L
            x = round(Int64, xlim/2 + s[sᵢ]*cos_[θᵢ] + t[l]*sin_[θᵢ])
            y = round(Int64, ylim/2 - t[l]*cos_[θᵢ] + s[sᵢ]*sin_[θᵢ])
            if 1 ≤ x ≤ xlim && 1 ≤ y ≤ ylim
                ℛf[sᵢ, θᵢ] += f[end-y+1, x]
            end
        end
    end

    return ℛf
end

"""
    backprojection(ℛf)(f::Matrix{T}, θ::Vector{T}) where {T<:Real}

return [back projection](https://freshrimpsushi.github.io/posts/back-projection-dual-transform-of-radon-transform/) of `ℛf`

# Keyword arguments

- `ℛf`: 2d array w.r.t. Radon transform of 'f'
- `θ`: projection angle(radian)
"""
function backprojection(ℛf, θ = 0:π/180:π-π/180)
    slim, θlim = size(ℛf)
    reconstructed_image = zeros(slim, slim)

    ylim, xlim = size(reconstructed_image)
    L = round(Int64, hypot(xlim, ylim))

    x= range(-xlim/2+1, xlim/2, xlim)
    y= range(ylim/2, -ylim/2+1, ylim)

    X = x'.*ones(xlim)
    Y = ones(ylim)'.*y

    s = range(-L/2, L/2, length=slim)

    for (value, angle) ∈ zip(eachcol(ℛf), θ)
        S = X.*cos(angle) + Y.*sin(angle)
        interpolation = LinearInterpolation(s, value, extrapolation_bc=Line())
        reconstructed_image .+= interpolation.(S)
    end
    return reconstructed_image
end

"""
    backprojection_anim(f::Matrix{T}, θ::Vector{T}) where {T<:Real}

return animation for calculation backprojection of ℛf

# Keyword arguments

- `ℛf`: 2d array w.r.t. Radon transform of 'f'
- `θ`: projection angle(radian)
"""
function backprojection_anim(ℛf)
    slim, θlim = size(ℛf)
    reconstructed_image = zeros(slim, slim)

    slim, θlim= size(ℛf)
    ylim, xlim = size(reconstructed_image)
    L = round(Int64, hypot(xlim, ylim))

    x= range(-xlim/2+1, xlim/2, xlim)
    y= range(ylim/2, -ylim/2+1, ylim)

    X = x'.*ones(xlim)
    Y = ones(ylim)'.*y

    θ = 0 : π/θlim : π-π/θlim
    s = range(-L/2, L/2, length=slim)

    anim = Animation()
    for (value, angle) ∈ zip(eachcol(ℛf), θ)
        S = X.*cos(angle) + Y.*sin(angle)
        interpolation = LinearInterpolation(s, value, extrapolation_bc=Line())
        reconstructed_image .+= interpolation.(S)
        frame(anim, heatmap(reverse(reconstructed_image, dims=1)))
    end
    for _ ∈ 1:80
        frame(anim, heatmap(reverse(reconstructed_image, dims=1)))
    end
    return anim
end

"""
    iradon(ℛf, method="fbp")

return Radon transform of `f`

# Keyword arguments

- `f`: phantom
- `θ`: projection angle
"""
function iradon(ℛf, method="fbp")

    if method == "fbp"
        reconstructed_f = fitered_backprojection(ℛf)
    end

    return reconstructed_f
end

"""
    fitered_backprojection(ℛff::Matrix{T}, θ::Vector{T}) where {T<:Real}
return f

# Keyword arguments

- `ℛf`: 2d array w.r.t. Radon transform of `f`
- `θ`: projection angle
"""
function fitered_backprojection(ℛf)
    N = size(ℛf)[1]

    ℱℛf = fft(ℛf, [1])

    fourier_filter = abs.(fftfreq(N))
    filterd_ℱℛf = ℱℛf .* fourier_filter

    ℱ⁻¹_filterd_ℱℛf = real.(ifft(filterd_ℱℛf, [1]))

    reconstructed_f = (π/360)*backprojection(ℱ⁻¹_filterd_ℱℛf)
    return reconstructed_f
end

"""
    SART(ℛf)
return Radon transform of `f`

# Keyword arguments

- `ℛf`: 2d array w.r.t. Radon transform of 'f'
- `θ`: projection angle
"""
function SART(ℛf)
end

end # module


text:
using JLD2
m = OneAssetHANK()#=
m <= Setting(:reduce_state_vars, false)
m <= Setting(:reduce_v, false)
=#
T, R, ~, inverse_basis = solvect(m)

n_vars = get_setting(m, :n_vars)
@show n_vars

#TODO: n_expectational_errors does not exist.
#n_expectation_errors = get_setting(m, :n_expectation_errors)
#n_shocks = get_setting(m, :n_shocks)
x = zeros(Float64, n_vars)

nst = size(T,1) ###
#C = zeros(nst)

n = 40 #100
s = zeros(nst,n) #s = zeros(n,406)
s[:,1] = inverse_basis' * x
e = randn(n)

dt = 1/2

for i = 2:n
    s[:,i] = s[:,i-1] + T*s[:,i-1]*dt + R*e[i]*(sqrt(dt))
end
data = zeros(1,Int(n/2))
data[1,:] = [s[64,2*i] for i = 1:Int(n/2)]
save("simulated_data.jld2", "data", data)


text:
using Zotero
using Test
using JSON

@testset "Zotero.jl" begin
    # Write your tests here.
end

c = ZoteroClient("test")
dicts = Zotero.request_json(c, "GET", "collections/LKXGUR69/items")
dict = dicts[1]
Zotero.dict_to_doc(dict)

l = Zotero.get_library(c)
print_tree(l)

Zotero.find_col(x->Zotero.title(x)=="Optimization", c, l)


text:
using DualNumbers

struct DualQuaternion{T<:Real} <: Number
  q0::Quaternion{T}
  qe::Quaternion{T}
  norm::Bool
end

DualQuaternion(q0::Quaternion, qe::Quaternion, n::Bool = false) =
  DualQuaternion(promote(q0, qe)..., n)

DualQuaternion(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool = false) =
  DualQuaternion(Quaternion(d1.re, d2.re, d3.re, d4.re, n),
                  Quaternion(d1.du, d2.du, d3.du, d4.du), n)

DualQuaternion(x::Real) = DualQuaternion(Quaternion(x), Quaternion(zero(x)), abs(x) == one(x))

DualQuaternion(d::Dual) = DualQuaternion(d, zero(d), zero(d), zero(d), abs(d) == one(d.re))

DualQuaternion(q::Quaternion) = DualQuaternion(q, zero(q), q.norm)

DualQuaternion(a::Vector) = DualQuaternion(zero(Quaternion{typeof(a[1])}), Quaternion(a))

convert(::Type{DualQuaternion{T}}, x::Real) where {T} =
    DualQuaternion(convert(Quaternion{T}, x), convert(Quaternion{T}, 0))

convert(::Type{DualQuaternion{T}}, d::Dual) where {T} =
    DualQuaternion(convert(Dual{T}, d), convert(Dual{T}, 0), convert(Dual{T}, 0), convert(Dual{T}, 0))

convert(::Type{DualQuaternion{T}}, q::Quaternion) where {T} =
    DualQuaternion(convert(Quaternion{T}, q), convert(Quaternion{T}, 0), q.norm)

convert(::Type{DualQuaternion{T}}, q::DualQuaternion{T}) where {T <: Real} = q

convert(::Type{DualQuaternion{T}}, dq::DualQuaternion) where {T} =
    DualQuaternion(convert(Quaternion{T}, dq.q0), convert(Quaternion{T}, dq.qe), dq.norm)

promote_rule(::Type{DualQuaternion{T}}, ::Type{T}) where {T <: Real} = DualQuaternion{T}
promote_rule(::Type{DualQuaternion}, ::Type{T}) where {T <: Real} = DualQuaternion
promote_rule(::Type{DualQuaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
promote_rule(::Type{Quaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
promote_rule(::Type{DualQuaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}

dualquat(q1, q2, n=false) = DualQuaternion(q1, q2, n)
dualquat(d1, d2, d3, d4, n=false) = DualQuaternion(d1, d2, d3, d4, n)
dualquat(x) = DualQuaternion(x)

function show(io::IO, dq::DualQuaternion)
  show(io, dq.q0)
  print(io, " + ( ")
  show(io, dq.qe)
  print(io, " )du")
end

Q0(dq::DualQuaternion) = dq.q0
Qe(dq::DualQuaternion) = dq.qe

(/)(dq::DualQuaternion, x::Real) = DualQuaternion(dq.q0 / x, dq.qe / x)

(/)(dq::DualQuaternion, d::Dual) =
  DualQuaternion(dual(dq.q0.s, dq.qe.s) / d,
                  dual(dq.q0.v1, dq.qe.v1) / d,
                  dual(dq.q0.v2, dq.qe.v2) / d,
                  dual(dq.q0.v3, dq.qe.v3) / d)

abs2(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) :
  dual(abs2(dq.q0),
        2.0 * (dq.q0.s  * dq.qe.s  +
                dq.q0.v1 * dq.qe.v1 +
                dq.q0.v2 * dq.qe.v2 +
                dq.q0.v3 * dq.qe.v3))

abs(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) : sqrt(abs2(dq))
float(dq::DualQuaternion{T}) where T = convert(DualQuaternion{float(T)}, dq)

conj(dq::DualQuaternion) = DualQuaternion(conj(dq.q0), conj(dq.qe), dq.norm)
dconj(dq::DualQuaternion) = DualQuaternion(dq.q0, -dq.qe, dq.norm)

inv(dq::DualQuaternion) = dq.norm ? conj(dq) : conj(dq) / abs2(dq)

function normalize(dq::DualQuaternion)
  if (dq.norm)
    return dq
  end
  a = abs(dq)
  if abs(a) > 0
    qa = dq / a
    dualquat(qa.q0, qa.qe, true)
  else
    dq
  end
end

function normalizea(dq::DualQuaternion)
  if (dq.norm)
    return (dq, one(dual))
  end
  a = abs(dq)
  if abs(a) > 0
    qa = dq / a
    dualquat(qa.q0, qa.qe, true), a
  else
    dq, zero(dual)
  end
end

(-)(dq::DualQuaternion) = DualQuaternion(-dq.q0, -dq.qe, dq.norm)

(+)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 + dw.q0, dq.qe + dw.qe)
(-)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 - dw.q0, dq.qe - dw.qe)
(*)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 * dw.q0,
                                                              dq.q0 * dw.qe + dq.qe * dw.q0,
                                                              dq.norm && dw.norm)
(/)(dq::DualQuaternion, dw::DualQuaternion) = dq * inv(dw)

function angleaxis(dq::DualQuaternion)
  tq = dq.qe * conj(dq.q0)
  t = [2.0 * tq.v1, 2.0 * tq.v2, 2.0 * tq.v3]
  q0s = dq.q0.s
  th0, s0 = angleaxis(dq.q0)
  sq0 = quat(0.0, s0)
  if abs(abs(q0s) - one(q0s)) == 0
    th = dual(th0, 0.5 * abs(quat(0, t)))
    th, dualquat(sq0)
  else
    th = dual(th0, 0.5 * dot(t, s0))
    s0c1 = cross(s0, t)
    tanth = tan(th0)
    s0c2 = (s0c1 / tanth + t) * 0.5
    sqev = cross(s0c2, s0)
    th, dualquat(sq0, quat(0.0, sqev))
  end
end

function angle(dq::DualQuaternion)
  th, ax = angleaxis(dq)
  th
end

function axis(dq::DualQuaternion)
  th, ax = angleaxis(dq)
  ax
end

function exp(dq::DualQuaternion)
  se = dual(dq.q0.s, dq.qe.s)
  se = exp(se)
  dq = dualquat(quat(0.0, imag(dq.q0)), quat(0.0, imag(dq.qe)))
  dq, th = normalizea(dq)
  if dq.norm
    dualquat(se) * (dualquat(cos(th)) + dq * dualquat(sin(th)))
  else
    dualquat(se)
  end
end

function log(dq::DualQuaternion)
  dq, a = normalizea(dq)
  sl = log(a)
  th, s = angleaxis(dq)
  s * dualquat(th) + dualquat(sl)
end

(^)(dq::DualQuaternion, dw::DualQuaternion) = exp(dw * log(dq))

function sqrt(dq::DualQuaternion)
  exp(0.5 * log(dq))
end

dualquatrand() = dualquat(quatrand(), quatrand())
ndualquatrand() = normalize(dualquatrand())


text:
function fixxing(f::Function, model::MetNet, ider, val::Real;
        btol = 0.0
    )
    idx = rxnindex(model, ider)
    bk_lb, bk_ub = model.lb[idx], model.ub[idx]
    try
        dflx = abs(val * btol)
        model.lb[idx] = val - dflx
        model.ub[idx] = val + dflx
        return f()
    finally
        model.lb[idx], model.ub[idx] = bk_lb, bk_ub
    end
end


text:
function lazyCheck(mdp::MDP, vi::ValueIteration)
  checkArgumentOrder(mdp)
  checkTransition(mdp)
  checkDiscretize(mdp, vi)
  info(string("mdp and value iteration solver passed basic checks"))
end

# Check that all |mdp| RangeVar variables have valid discretization schemes in |vi|
function checkDiscretize(mdp::MDP, vi::ValueIteration)

  for lazyvar in values(mdp.statemap)
    if isa(lazyvar, RangeVar)
      if !haskey(vi.statemap, lazyvar.varname)
        error(string(
          "state variable ", lazyvar.varname,
          " does not have a discretization scheme"))
      elseif lazyvar.maxval - lazyvar.minval < vi.statemap[lazyvar.varname].step
        error(string(
          "state variable ", lazyvar.varname,
          " has a discretization step size larger than than its range"))
      end
    end
  end

  for lazyvar in values(mdp.actionmap)
    if isa(lazyvar, RangeVar)
      if !haskey(vi.actionmap, lazyvar.varname)
        error(string(
          "action variable ", lazyvar.varname,
          " does not have a discretization scheme"))
      elseif lazyvar.maxval - lazyvar.minval < vi.actionmap[lazyvar.varname].step
        error(string(
          "action variable ", lazyvar.varname,
          " has a discretization step size larger than than its range"))
      end
    end
  end

end


text:
struct Simulation{FT<:Real}
     δt :: FT
     tt :: Float64
     ft :: Float64
    fnc :: AbstractString
end

function CreateSimulation(
	FT = Float64;
	δt  :: Real,
	tt  :: Real,
	ft  :: Real,
	fnc :: AbstractString
)

	return Simulation{FT}(δt,tt,ft,fnc)

end

function show(io::IO, S::Simulation{FT}) where FT <: Real
    print(
		io,
		"Gill-Matsuno Simulation{$FT} Setup:\n",
		"           Time Step (δt) : ", S.δt, '\n',
		"           Stop Time (tt) : ", S.tt, '\n',
		"    Output Frequency (ft) : ", S.ft, '\n',
		"    NetCDF Filepath (fnc) : ", S.fnc, '\n'
	)
end


text:
"""
$(SIGNATURES)

Returns Nemo's `MatElem` paramatrising the Jacobian of the dynamical system associated to a biochemical reaction network. We use "convex parameters", see
```julia
rn = @reaction_network begin
    k1,  X2 → X1
    k2,  X1 → X2
end k1 k2

N = netstoichmat(rn)
Y = substoichmat(rn)

E = raysof(cone_positivenullspace(N))
W = conservationlaws(N)
```
"""
function Jacobian_cp(N::MatElem, Y::MatElem, E::AbstractMatrix, ls, hs)
    ## Diagonal with coordinates of E
    D = diagonal(base_ring(N), E * ls)
    ## Square matrix with h's diagonal
    H = diagonal(base_ring(N), hs)
    ## Return the jacobian under convex parameters
    return N * D * Y * H
end

function Jacobian_cp(S,
                     N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix;
                     ls=gens(S)[1:size(E, 2)],
                     hs=gens(S)[.+(size(E, 2), 1:size(N, 1))])
    return Jacobian_cp(matrix(S, N),
                       matrix(S, permutedims(Y)),
                       E,
                       ls,
hs)
end

Jacobian_cp(H, L,
            N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix) =
                Jacobian_cp(H,
                            N, Y, E;
                            ls=gens(L), hs=gens(H))

## Computing the coeffitient of the least degree term (ldc) of the characteristic polynomial of the Jacobian, parametrised using convex parameters.
##
## Computing the whole char polynomial
function ldc_charJacobian_cp(S,
                             N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix;
                             ls=gens(S)[1:size(E, 2)],
                             hs=gens(S)[.+(size(E, 2), 1:size(N, 1))],
                             mindeg=size(N, 1) - rank(matrix(FlintZZ, N)))
    return coeff_charpoly(Jacobian_cp(S, N, Y, E; ls=ls, hs=hs), mindeg)
end

function ldc_charJacobian_cp(H, L,
                             N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix;
                             mindeg=size(N, 1) - rank(matrix(FlintZZ, N)))
    return ldc_charJacobian_cp(H,
                               N, Y, E;
                               ls=gens(L), hs=gens(H),
                               mindeg=mindeg)
end

## Computing the coeffitient of the least degree term (ldc) of the characteristic polynomial of the Jacobian, parametrised using convex parameters.
##
## If we know W, we can use it, simply add it as a last argument!
function ldc_charJacobian_cp(S,
                             N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix, W::AbstractMatrix;
                             ls=gens(S)[1:size(E, 2)],
                             hs=gens(S)[.+(size(E, 2), 1:size(N, 1))],
                             mindeg=size(N, 1) - rank(matrix(FlintZZ, N)))
    J = Jacobian_cp(S, N, Y, E; ls=ls, hs=hs)
    for (i, j) in enumerate(findpivotsof(W)(1))
        # Change row j-th of J by i-th of W
        # Nemo do not allow broadcasting.
        for k in 1:size(J, 2)
            J[j,k] = W[i,k]
        end
    end
    return det(J)
end

function ldc_charJacobian_cp(H, L,
                             N::AbstractMatrix, Y::AbstractMatrix, E::AbstractMatrix, W::AbstractMatrix;
                             mindeg=size(N, 1) - rank(matrix(FlintZZ, N)))
    return ldc_charJacobian_cp(H,
                               N, Y, E, W;
                               ls=gens(L), hs=gens(H),
                               mindeg=mindeg)
end

function from_lh_to_kx(point, E, Y, nls=size(E, 2))
    xs = inv.(point[nls + 1:end])
    e = E * point[begin:nls]
    x_to_y = [prod(xs.^r) for r in eachcol(Y)]
    ks = e .* inv.(x_to_y)
    return ks, xs
end


text:
"""
    update_dataset(stat_list, save_folder; path_runoff)

Update an existing initial dataset.
"""
function update_dataset(stat_list, save_folder;
    path_runoff = "/home/jmg/flood_forecasting/runoff_observed")

    # Read metadata from excel sheet

    @info "Read metadata"

    df_meta = read_metadata()

    df_meta = subset_metadata(df_meta, stat_list)

    # Dictonary that links drainage basin keys to senorge indices

    @info "Link drainage basin keys to senorge indices"

    dbk_ind = read_dbk_ind()

    dbk_ind = subset_dbk_ind(dbk_ind, df_meta)

    # Dictonary that links drainage basin keys to dataframes with
    # geographic information (landuse, elevation...)

    @info "Link drainage basin keys to dataframes with geoinfo"

    df_geo = watershed_info(dbk_ind)

    # Read existing data

    @info "Read existing data"

    tair_old, time_old = read_old_data(save_folder, "tair")

    prec_old, time_old = read_old_data(save_folder, "prec")

    # Find last available senorge data

    time_first = time_old[end] + Dates.Day(1)
    time_last = find_last_senorge("V2.0")

    time_vec = time_first:Day(1):time_last

    # Read new meteorological data

    met_var = "tm"

    tair_new, time_new = read_senorge_data(time_vec, met_var, df_geo; scale_func = x -> (x.-2732.0)./10.0)

    met_var = "rr"

    prec_new, time_new = read_senorge_data(time_vec, met_var, df_geo; scale_func = x -> x./10.0)

    # Write data to files

    write_update(save_folder, "tair", time_old, tair_old, time_new, tair_new, df_meta)

    write_update(save_folder, "prec", time_old, prec_old, time_new, prec_new, df_meta)

    # Process runoff data

    time_first = time_old[1]

    time_vec = time_first:time_last

    process_runoff_data(stat_list, time_vec, df_meta, path_runoff, save_folder)

end

max_forks_repo_name: jmgnve/NveData


text:
using Metalhead, Test
using Flux

@testset "ViT" begin
  for mode in [:tiny, :small, :base, :large] #,:huge, :giant, :gigantic]
    m = ViT(mode)
    @test size(m(rand(Float32, 256, 256, 3, 1))) == (1000, 1)
    @test_skip gradtest(m, rand(Float32, 256, 256, 3, 1))
  end
  GC.gc()
end


text:
## Get vertex and edge data

"""
    get_data(g)

Retrieve graph-level metadata.
"""
get_data(g::DataDiGraph) = g.graph_data

"""
    get_data(g, v)

Retrieve vertex-level metadata for vertex `v`.
"""
get_data(g::DataDiGraph, v::Integer) = g.vertex_data[v]

"""
    get_data(g, s, d)

Retrieve edge-level metadata for edge `(s, d)`.
"""
function get_data(g::DataDiGraph, s::Integer, d::Integer)
    d_index = searchsortedfirst(g.fadjlist[s], d)
    return g.edge_data[s][d_index]
end

## Set vertex and edge data

"""
    set_data!(g, v, data)

Set vertex-level metadata to value `data` for vertex `v`.
"""
function set_data!(g::DataDiGraph{T,VL,VD}, v::Integer, data::VD) where {T,VL,VD}
    if has_vertex(g, v)
        g.vertex_data[v] = data
        return true
    else
        return false
    end
end

"""
    set_data!(g, s, d, data)

Set edge-level metadata to value `data` for edge `(s, d)`.
"""
function set_data!(
    g::DataDiGraph{T,VL,VD,ED}, s::Integer, d::Integer, data::ED
) where {T,VL,VD,ED}
    if has_edge(g, s, d)
        d_index = searchsortedfirst(g.fadjlist[s], d)
        g.edge_data[s][d_index] = data
        return true
    else
        return false
    end
end


text:
using ApproxFun, BandedMatrices, BlockBandedMatrices, LinearAlgebra, BlockArrays, MultivariateOrthogonalPolynomials, Test
using SparseVolterraExamples

#############################
## Tests for Volterra integral equations of first and second kind
## with non-Clenshaw monomial expansion of the kernel.
####
## Includes tests for 1 to x and 1 to 1-x limits of integration.
#############################

@testset "Volterra Integral Equation 1st kind, mon. kernels, flip=false" begin
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> 1-sin(x)^2, S)
    # K given in monomial basis
    K = [2.3,0.33,2.1,1.2,7.6]  # = K(x,y)=2.3+0.33*y+2.1*x+1.2*x*y+7.6*y^2
        V=triVolterraOpP01(K,N,false)
        x = Fun(identity,S)
        J = Multiplication(x, S)[1:N,1:N]
    VintFun1=Fun(Jacobi(0,1,0..1),(V-J*V)*pad(g.coefficients,N))
    VintFun2=triVolterraFunP01(g,K,N,false)
    Vtest=triVolterraEQ1Solver(VintFun1,K,N,false)
    @test VintFun1 ≈ VintFun2
    @test Vtest ≈ g
    #########################################################################
    # sin(1-x) = integral_0^(1-x) u(t) dt
    # analytic solution is u(x)=cos(x)
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> sin(1-x), S)
    K = [1.,0.,0.,0.]  # = K(x,y)=1
    @test triVolterraEQ1Solver(g,K,N,false) ≈ Fun(x -> cos(x), S)
    #########################################################################
    # 2-x-exp(1-x) = integral_0^(1-x) (y-1-x) u(y) dy
    # analytic solution is u(x)=exp(x)
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> 2-x-exp(1-x), S)
    K = [-1.0,1.0,1.0,0.]  # = K(x,y)=(y-1+x)
    @test triVolterraEQ1Solver(g,K,N,false) ≈ Fun(x -> exp(x), S)
    #########################################################################
    # 3-x-2*exp(1-x)+(1-x)*exp(1-x) = integral_0^(1-x) (1-x-y) u(y) dy
    # analytic solution is u(x)=x*exp(x)
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> 3-x-2*exp(1-x)+(1-x)*exp(1-x), S)
    K = [1.,-1.,-1.0,0.]  # = K(x,y)=(1-x-y)
    @test triVolterraEQ1Solver(g,K,N,false) ≈ Fun(x -> x*exp(x), S)
end

@testset "Volterra Integral Equation 1st kind, mon. kernels, flip=true" begin
    g=Fun(x -> sin(x)-(x)*cos(x), Jacobi(0,1,0..1))
    K = [0.,1.,0.,0.]
    N = 20
    @test triVolterraEQ1Solver(g,K,N,true)  ≈ Fun(x->sin(x),Jacobi(0,1,0..1))
    #########################################################################
    g=Fun(x -> 2+x-2*exp(x)+x*exp(x), Jacobi(0,1,0..1))
    K = [0.,-1.,1.,0.]
    N = 15
    @test triVolterraEQ1Solver(g,K,N,true) ≈ Fun(x->x*exp(x),Jacobi(0,1,0..1))
end

@testset "Volterra Integral Equation 2nd kind, mon. kernels, flip=false" begin
    # u(x) = (x-1)*cos(2*pi*(1-x))/(2*pi)+((4*pi^2+1)/(4*pi^2))*sin(2*pi*(1-x))+\int_0^(1-x) t*u(t)dt
    # solution is u(x) = sin(2*pi*(1-x))
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> (x-1)*cos(2*pi*(1-x))/(2*pi)+((4*pi^2+1)/(4*pi^2))*sin(2*pi*(1-x)), S)
    K = zeros(5)
    K[2]=1.
    @test triVolterraEQ2Solver(g,K,N,false) ≈ Fun(x -> sin(2*pi*(1-x)), S)
    #########################################################################
    # u(x) = (-17/30)*x^6+2*x^5-(3/2)*x^4-(11/2)*x^3+(15/2)*x^2-(21/10)*x+1/6 + \int_0^(1-x) (7xy+5y^2) u(y)dy
    # solution is u(x) = 1-x^3
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> (-17/30)*x^6+2*x^5-(3/2)*x^4-(11/2)*x^3+(15/2)*x^2-(21/10)*x+1/6, S)
    K = [0.,0.,0.,7.,5.,0.,0.,0.,0.]
    @test triVolterraEQ2Solver(g,K,N,false) ≈ Fun(x -> 1-x^3, S)
    #########################################################################
    # u(x) = -x+x^2/2+1/2 + \int_0^(1-x) u(t)dt
    # solution is u(x) = 1-x
    N=15
    S = Jacobi(0,1,0..1)
    g = Fun(x -> -x+x^2/2+1/2, S)
    K = [1.,0.,0.]
    @test triVolterraEQ2Solver(g,K,N,false) ≈ Fun(x -> 1-x, S)
end

@testset "Volterra Integral Equation 2nd kind, mon. kernels, flip=true" begin
    g=Fun(x -> 1, Jacobi(0,1,0..1))
    K = [0.,1.,-1.,0.]
    N = 10
    @test triVolterraEQ2Solver(g,K,N,true)  ≈ Fun(x->cos(x),Jacobi(0,1,0..1))
    #########################################################################
    g=Fun(x -> 1-x-0.5*x^2, Jacobi(0,1,0..1))
    K = [0.,-1.,1.,0.]
    N = 10
    @test triVolterraEQ2Solver(g,K,N,true)  ≈ Fun(x->1-sinh(x),Jacobi(0,1,0..1))
    #########################################################################
    g=Fun(x -> x+x^4+0.5*x^2+0.2*x^5, Jacobi(0,1,0..1))
    K = [-1.,0.,0.,0.]
    N = 25
    @test triVolterraEQ2Solver(g,K,N,true)  ≈ Fun(x->x+x^4,Jacobi(0,1,0..1))
end


text:
using SimpleCache
using Test

@testset "SimpleCache.jl" begin
    cachename = tempname()
    @test cached(()->42, cachename) == 42
    @test cached(()->43, cachename) == 42
    @test cached(()->44, cachename; reset=true) == 44
    rm(cachename)
end


text:
function iszero{T}(x::T)
    return x == zero(T)
end

function Maybe(T::Type)
    return Union(T, Nothing)
end

# Cycle-zip. Zip two or more arrays, cycling the short ones.
function cyclezip(xs::AbstractArray...)
    if any(map(isempty, xs))
        return {}
    end
    n = maximum([length(x) for x in xs])
    return takestrict(zip([cycle(x) for x in xs]...), n)
end

# This generates optimized code for a reoccuring pattern in forms and patters
# that looks like:
#
#   return Circle([CirclePrimitive(Point(x, y), x_measure(r))
#                  for (x, y, r) in cyclezip(xs, ys, rs)])
#
# This macro does the equivalent with
#
#   return Circle(@makeprimitives CirclePrimitive,
#                        (x in xs, y in ys, r in rs),
#                        CirclePrimitive(Point(x, y), x_measure(r)))
#
# but much more efficiently.
macro makeprimitives(args)
    @assert args.head == :tuple
    @assert length(args.args) == 3
    T, iterators, constructor = args.args

    maxlen_ex = quote begin n = 0 end end
    iter_ex = quote begin end end

    for iterator in iterators.args
        @assert iterator.head == :in
        var = iterator.args[1]
        arr = iterator.args[2]

        push!(maxlen_ex.args, quote
            if isempty($(arr))
                primitives = Array($(T), 0)
                @goto done
            end end)
        push!(maxlen_ex.args, quote n = max(n, length($(arr))) end)
        push!(iter_ex.args, quote
            $(var) = $(arr)[((i - 1) % length($(arr))) + 1]
        end)
    end

    quote
        $(maxlen_ex)
        primitives = Array($(T), n)
        for i in 1:n
            $(iter_ex)
            primitives[i] = $(constructor)
        end
        @label done
        primitives
    end
end


text:
# Interactive clickable FitzHughNagumo state space

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using GLMakie, DynamicalSystems, InteractiveDynamics

a = 3.
b = -0.05
ε = 0.01
I = 0.0

fh = Systems.fitzhugh_nagumo(zeros(2); a, b, ε, I)

# Layout plot:
fig = Figure(resolution = (1000, 800))
ax = fig[1,1] = Axis(fig)
display(fig)
ax.xlabel = "u"
ax.ylabel = "w"
ax.title = "FitzHugh Nagumo, a=$a, b=$b"

# plot the two nullclines
us = -0.4:0.01:1.2
w1 = us
w2 = @. a*us*(us-b)*(1-us) + I
lines!(ax, us, w1; color = COLORS[1], linewidth = 2, linestyle = :dash)
lines!(ax, us, w2; color = COLORS[2], linewidth = 2, linestyle = :dot)

# Create trajectories on point selection
spoint = select_point(ax.scene)
on(spoint) do pos
    tr = trajectory(fh, 10000, SVector(pos...))
    lines!(ax, columns(tr)...;
        color = InteractiveDynamics.randomcolor(), linewidth = 4.0
    )
end

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 10)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module GML

using ...ParserCombinator
using Compat

export parse_raw, parse_dict, GMLError

# symbol required in 0.3
# both work in 0.4
# symbol gives deprecation warning in 0.5
Symbol_ = VERSION >= v"0.5-" ? Symbol : symbol

function mk_parser(string_input)

    # this is such a simple grammar that we don't need backtracking, so we can
    # use Seq! et al, and Error for useful diagnostics.

    # the only tricky things are getting the spaces right so that matching
    # spaces doesn't commit us to anything unexpected, and placing errors only
    # when we're sure we're wrong (you can't have one in the definition of
    # key, for example, because that can fail...).

    @with_names begin

        expect(x) = Error("Expected $x")

        pint(x) = parse(Int64, x)
        pflt(x) = parse(Float64, x)

        comment = P"(#.*)?"

        # we need string input as we match multiple lines
        if string_input && ParserCombinator.FAST_REGEX

            wspace   = "([\t ]+|[\r\n]+(#.*)?)"
            wstar(x) = string(x, wspace, "*")
            wplus(x) = string(x, wspace, "+")
            space    = ~Pattern(wplus(""))
            spc      = ~Pattern(wstar(""))

            open     = ~Pattern(wstar("\\["))
            close    = ~Pattern(wstar("]"))

            key      = Pattern(wplus("([a-zA-Z][a-zA-Z0-9]*)"), 1)    > Symbol_
            int      = Pattern(wstar("((\\+|-)?\\d+)"), 1)            > pint
            real     = Pattern(wstar("((\\+|-)?\\d+.\\d+((E|e)(\\+|-)?\\d+)?)"), 1) > pflt
            str      = Pattern(wstar("\"([^\"]*)\""), 1)

        else

            wspace   = Alt!(P"[\t ]+", Seq!(P"[\r\n]+", comment))
            space    = wspace[1:end,:!]
            spc      = wspace[0:end,:!]

            open     = Seq!(E"[", spc)
            close    = Seq!(E"]", spc)

            key      = Seq!(p"[a-zA-Z][a-zA-Z0-9]*", space)           > Symbol_
            int      = Seq!(p"(\+|-)?\d+", spc)                       > pint
            real     = Seq!(p"(\+|-)?\d+.\d+((E|e)(\+|-)?\d+)?", spc) > pflt
            str      = Seq!(Pattern("\"([^\"]*)\"", 1), spc)

        end

        list    = Delayed()
        sublist = Seq!(open, list, Alt!(close, expect("]")))
        value   = Alt!(real, int, str, sublist, expect("value"))
        element = Seq!(key, value)                            > tuple

        list.matcher = Nullable{Matcher}(element[0:end,:!]    > vcat)

        # first line comment must be explicit (no previous linefeed)
        Seq!(comment, spc, list, Alt!(Seq!(spc, Eos()), expect("key")))

    end
end

# this returns the "natural" representation as nested arrays and tuples
function parse_raw(s; debug=false)
    parser = mk_parser(isa(s, AbstractString))
    try
        if ParserCombinator.FAST_REGEX
            (debug ? parse_one_dbg : parse_one)(s, Trace(parser); debug=debug)
        else
            (debug ? parse_lines_dbg : parse_lines)(s, Trace(parser); debug=debug)
        end
    catch x
        if (debug)
            Base.show_backtrace(STDOUT, catch_backtrace())
        end
        rethrow()
    end
end

# (semi) structured model of GML graph files

# the GML specs that i have found are really rather frustrating, because they
# don't seme to acknowledge a fundamental problem with this format, which is
# that you cannot tell, from the file alone, whether a particuar field is a
# list or a single value.

# obviously, a name that occurs multiple times in a single scope is a list.
# but the opposite - that an isolated name is a single value - is not
# necessarily true, because it may be a singleton list.

# one solution is to make a model that very closely follows the "predefined
# keys" part of himsolt's 1996 document, available as part of the tarball from
# http://www.fim.uni-passau.de/fileadmin/files/lehrstuhl/brandenburg/projekte/gml/gml-documentation.tar.gz
# but that seems very specific, perhaps dated, and still doesn't help with
# additional fields.

# another solution is to treat everything as a list, and use dictionaries of
# lists.  but that means that the idea of "keys" is messed up with additional
# [1] indexes into singleton lists.

# after some reflection, i've decide to take a list of names of lists, and to
# validate against that.  this isn't perfect - it doesn't allow for the same
# name to have different meanings in different contexts, for example - but it
# seems to be a good middle ground for "doing the right thing" in general
# cases.

# users with different requirements are free to take the "raw" parse and build
# their own object models.

type GMLError<:Exception
    msg::AbstractString
end

LISTS = [:graph,:node,:edge]

typealias GMLDict Dict{Symbol, Any}

function build_dict(raw; lists=LISTS, unsafe=false)
    root = GMLDict()
    if length(raw) > 0
        build_dict(root, raw[1]; lists=lists, unsafe=unsafe)
    end
    root
end

function build_dict(dict::GMLDict, raw; lists=LISTS, unsafe=false)
    for (name, value) in raw
        if isa(value, Vector)
            entry = GMLDict()
            build_dict(entry, value; lists=lists, unsafe=unsafe)
        else
            entry = value
        end
        if name in lists
            if !haskey(dict, name)
                dict[name] = Any[]
            end
            push!(dict[name], entry)
        else
            if !haskey(dict, name)
                dict[name] = entry
            elseif !unsafe
                throw(GMLError("$name is a list"))
            end
        end
    end
end

# lists describes which symbols should be modelled as lists
# if unsafe is false, multiple values for non-list symbols throw an error;
# if true they are silently discarded
parse_dict(s; debug=false, lists=LISTS, unsafe=false) = build_dict(parse_raw(s; debug=debug); lists=lists, unsafe=unsafe)

end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.type;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.thing.Entity;
#
# import javax.annotation.CheckReturnValue;
# import java.util.stream.Stream;
#
# public interface EntityType extends ThingType {
#
#     @Override
#     @CheckReturnValue
#     default boolean isEntityType() {
#         return true;
#     }
#
#     @Override
#     @CheckReturnValue
#     EntityType.Remote asRemote(TypeDBTransaction transaction);
#
#     interface Remote extends ThingType.Remote, EntityType {
#
#         @CheckReturnValue
#         Entity create();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends Entity> getInstances();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends EntityType> getSubtypes();
#
#         void setSupertype(EntityType superEntityType);
#     }
# }


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

function relaxation(distance::AbstractArray, auxdistances::AbstractArray, cutoff::Real)
  # patterns enabled in the training image
  enabled = .!isinf.(distance)
  npatterns = sum(enabled)

  # candidates with good overlap
  dbsize = all(distance[enabled] .== 0) ? npatterns : ceil(Int, cutoff*npatterns)
  overlapdb = partialsortperm(vec(distance), 1:dbsize)

  # candidates in accordance with auxiliary data
  naux = length(auxdistances)
  softdb = [Vector{Int}() for i in 1:naux]

  patterndb = []
  softdistance = [copy(auxdistances[i]) for i in 1:naux]
  frac = .1 * (dbsize / npatterns)
  while true
    softdbsize = ceil(Int, frac*npatterns)

    patterndb = overlapdb
    for n in 1:naux
      softdistance[n][softdb[n]] .= Inf
      softdb[n] = [softdb[n]; partialsortperm(vec(softdistance[n]), 1:softdbsize-length(softdb[n]))]

      patterndb = fast_intersect(patterndb, softdb[n], length(distance))

      isempty(patterndb) && break
    end

    !isempty(patterndb) && break
    frac = min(frac + .1, 1)
  end

  patterndb
end

function fast_intersect(A::AbstractVector{Int}, B::AbstractVector{Int},
                         nbits::Integer)
  bitsA = falses(nbits)
  bitsB = falses(nbits)
  bitsA[A] .= true
  bitsB[B] .= true

  findall(bitsA .& bitsB)
end


text:
function get_account_settings(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/account/settings.json", options))
end

function get_account_verify_credentials(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/account/verify_credentials.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
#=
 This file declares as much of the Plugin API (as dictated by
 "vim_get_api_info") as possible.
=#

function declare_err(api::Dict)
    # function Base.showerr(io::IO, err::NeoVimError)
    #  if err.id == 0
    #   print(io, "NeoVim ", "Validation", " Error: ", err.msg)
    #  else
    #   print(io, "NeoVim Unknown Error: ", err.msg)
    #  end
    # end
    # pretty much.
    # TODO: This works but I'm fairly sure it can be done much more simply.
    #       Do I need to use Expr(...) in the push! or can I use :(...)?
    ex = Expr(
        :function,
        :(Base.showerror(io::IO, err::NeoVimError))
    )
    currex = ex
    for (id, name) in (Int64,ByteString)[(v["id"], k) for (k, v) in api]
        push!(currex.args, Expr(
            :if,
            :(err.id == $id),
            :(print(io, "NeoVim ", $name, " Error: ", err.msg))
        ))
        currex = currex.args[end]
    end
    push!(currex.args, :(print(io, "NeoVim Unknown Error: ", err.msg)))
    eval(ex)
end

function declare_type(api::Dict)
    # function sanitize(n::Nvim, vt::MsgPack.Ext)
    #  if vt.typecode == 0
    #   Buffer(n, vt)
    #  else
    #   error()
    #  end
    # end
    # or words to that effect.
    ex = Expr(
        :function,
        :(sanitize(n::Nvim, vt::MsgPack.Ext))
    )
    currex = ex
    for (id, typ) in (Int64,Symbol)[(v["id"], k) for (k, v) in api]
        push!(currex.args, Expr(
            :if,
            :(vt.typecode == $id),
            :($(typ)(n, vt))
        ))
        currex = currex.args[end]
    end
    push!(currex.args, :(error("OH MY GOD WHAT HAVE THEY DONE!")))
    eval(ex)
end

function sanitize(n::Nvim, d::Dict)
    @compat Dict{UTF8String,Any}([bytestring(k) => sanitize(n, v) for (k, v) in d])
end
sanitize(n::Nvim, a::Vector) = [sanitize(n, i) for i in a]
sanitize(::Nvim, a::Vector{UInt8}) = bytestring(a)
sanitize(::Nvim, x) = x

api_info(n::Nvim) = request(n, "vim_get_api_info")[2]

function vimtojulia(s)
    if s == "Nil"
        return :Void
    elseif s == "Boolean"
        return :Bool
    elseif s == "Integer"
        return :Int64
    elseif s == "Float"
        return :Float64
    elseif s == "String"
        return :ByteString
    elseif s == "Array"
        return :Vector
    elseif s == "Dictionary"
        return :Dict
    else
        return :Any
    end
end

function declare_func(api::Vector)
    for d in api
        vimnm = d["name"]
        # need to know the type of the first argument
        @compat typenm = split(vimnm, '_', keep=false)
        typenm, funcnm = typenm[1], join(typenm[1:end], "_")

        # TODO: would ignoring the types make this simpler?
        args = (Symbol, Symbol)[]
        for (typ, nam) in d["parameters"]
            nam == "end" && (nam = "finish")
            push!(args, (nam, vimtojulia(typ)))
        end

        if typenm == "vim"
            typenm = "nvim"
        else
            # the type is already given for non "vim_" commands so remove it
            shift!(args)
        end
        unshift!(args, (typenm, ucfirst(typenm)))

        eval(Expr(
            :function,
            Expr(:call, symbol(funcnm), [:($nm::$tnm) for (nm, tnm) in args]...),
            Expr(
                :call,
                :request,
                if args[1][2] == :Nvim
                    args[1][1]
                else
                    :($(args[1][1]).nvim)
                end,
                vimnm,
                (if args[1][2] == :Nvim
                    [nm for (nm, _) in args[2:end]]
                else
                    cont = shift!(args)[1]
                    vcat(:($(cont).data), [nm for (nm, _) in args])
                end)...
            )
        ))
    end
end