text:
export SSPRK22Heuns, SSPRK22Ralstons, SSPRK33ShuOsher, SSPRK34SpiteriRuuth

"""
    StrongStabilityPreservingRungeKutta <: DistributedODEAlgorithm

A class of Strong Stability Preserving Runge--Kutta methods.
These require two additional copies of the state vector.

The available implementations are:

  - [`SSPRK22Heuns`](@ref)
  - [`SSPRK22Ralstons`](@ref)
  - [`SSPRK33ShuOsher`](@ref)
  - [`SSPRK34SpiteriRuuth`](@ref)
"""
abstract type StrongStabilityPreservingRungeKutta <: DistributedODEAlgorithm end

struct StrongStabilityPreservingRungeKuttaTableau{Nstages, RT}
    "Storage RK coefficient vector A1 (rhs scaling of u)"
    A1::NTuple{Nstages, RT}
    "Storage RK coefficient vector A2 (rhs scaling of U)"
    A2::NTuple{Nstages, RT}
    "Storage RK coefficient vector B (rhs add in scaling)"
    B::NTuple{Nstages, RT}
    "Storage RK coefficient vector C (time scaling)"
    C::NTuple{Nstages, RT}
end

struct StrongStabilityPreservingRungeKuttaCache{Nstages, RT, A}
    tableau::StrongStabilityPreservingRungeKuttaTableau{Nstages, RT}
    "Storage for RHS during the `StrongStabilityPreservingRungeKutta` update"
    fU::A
    "Storage for the stage state during the `StrongStabilityPreservingRungeKutta` update"
    U::A
end

function cache(
    prob::DiffEqBase.AbstractODEProblem{uType, tType, true},
    alg::StrongStabilityPreservingRungeKutta; kwargs...) where {uType,tType}

    tab = tableau(alg, eltype(prob.u0))
    # can't use Vector{T}(undef) as need to ensure no NaNs
    fU = zero(prob.u0)
    U = zero(prob.u0)
    return StrongStabilityPreservingRungeKuttaCache(tab, fU, U)
end

function step_u!(int, cache::StrongStabilityPreservingRungeKuttaCache{Nstages, RT, A}) where {Nstages, RT, A}
    tab = cache.tableau

    f! = int.prob.f
    u = int.u
    p = int.prob.p
    t = int.t
    dt = int.dt

    for s in 1:Nstages
        if s == 1
            f!(cache.fU, u, p, t + tab.C[s]*dt)
        else
            f!(cache.fU, cache.U, p, t + tab.C[s]*dt)
        end
        if s < Nstages
            cache.U .= tab.A1[s] .* u  .+  tab.A2[s] .* cache.U .+ (dt * tab.B[s]) .* cache.fU
        else
            u .= tab.A1[s] .* u  .+  tab.A2[s] .* cache.U .+ (dt * tab.B[s]) .* cache.fU
        end
    end
end

"""
    SSPRK22Heuns()

The second-order, 2-stage, strong-stability-preserving, Runge--Kutta scheme of
[SO1988](@cite), also known as Heun's method ([Heun1900](@ref).

Exact choice of coefficients from wikipedia page for Heun's method :)
"""
struct SSPRK22Heuns <: StrongStabilityPreservingRungeKutta end

function tableau(::SSPRK22Heuns, RT)
    RKA1 = (RT(1), RT(1 // 2))
    RKA2 = (RT(0), RT(1 // 2))
    RKB  = (RT(1), RT(1 // 2))
    RKC  = (RT(0), RT(1))
    StrongStabilityPreservingRungeKuttaTableau(RKA1, RKA2, RKB, RKC)
end

"""
    SSPRK22Ralstons()

The second-order, 2-stage, strong-stability-preserving, Runge--Kutta scheme
of [SO1988](@cite), also known as Ralstons's method ([Rals1962](@cite).

Exact choice of coefficients from wikipedia page for Heun's method :)
"""
struct SSPRK22Ralstons <: StrongStabilityPreservingRungeKutta end

function tableau(::SSPRK22Ralstons, RT)
    RKA1 = (RT(1), RT(5 // 8))
    RKA2 = (RT(0), RT(3 // 8))
    RKB  = (RT(2 // 3), RT(3 // 4))
    RKC  = (RT(0), RT(2 // 3))
    StrongStabilityPreservingRungeKuttaTableau(RKA1, RKA2, RKB, RKC)
end

"""
    SSPRK33ShuOsher()

The third-order, 3-stage, strong-stability-preserving, Runge--Kutta scheme
of [SO1988](@cite).
"""
struct SSPRK33ShuOsher <: StrongStabilityPreservingRungeKutta end

function tableau(::SSPRK33ShuOsher, RT)
    RKA1 = (RT(1), RT(3 // 4), RT(1 // 3))
    RKA2 = (RT(0), RT(1 // 4), RT(2 // 3))
    RKB = (RT(1), RT(1 // 4), RT(2 // 3))
    RKC = (RT(0), RT(1), RT(1 // 2))
    StrongStabilityPreservingRungeKuttaTableau(RKA1, RKA2, RKB, RKC)
end

"""
    SSPRK34SpiteriRuuth()

The third-order, 4-stage, strong-stability-preserving, Runge--Kutta scheme
of [SR2002](@cite).
"""
struct SSPRK34SpiteriRuuth <: StrongStabilityPreservingRungeKutta end

function tableau(::SSPRK34SpiteriRuuth, RT)
    RKA1 = (RT(1), RT(0), RT(2 // 3), RT(0))
    RKA2 = (RT(0), RT(1), RT(1 // 3), RT(1))
    RKB  = (RT(1 // 2), RT(1 // 2), RT(1 // 6), RT(1 // 2))
    RKC  = (RT(0), RT(1 // 2), RT(1), RT(1 // 2))
    StrongStabilityPreservingRungeKuttaTableau(RKA1, RKA2, RKB, RKC)
end


text:
using Multiplexers
using Sockets
using Serialization

N = 1
@sync begin
    @async let
        server = listen(2014)
        try
            socket = accept(server)
            mux = Multiplexer(socket,N)

            # Tesing Serializer
            serialize(mux.lines[1],"from serializer")
            @show deserialize(mux.lines[1])

            # Testing asynchronous conection
            @async serialize(mux.lines[1],"Hello from here")
            @show deserialize(mux.lines[1])

            close(mux)
        finally
            close(server)
        end
    end

    @async let
        socket = connect(2014)
        mux = Multiplexer(socket,N)

        # Testing serializer
        @show deserialize(mux.lines[1])
        serialize(mux.lines[1],"HEllo")

        # Testing asynchronous conection
        @async serialize(mux.lines[1],"Hello from there")
        @show deserialize(mux.lines[1])

        wait(mux)
    end
end

ext: jl


text:
# Using matplotlib, make a stick figure showing the disk inclination
using consonance

push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

using constants
using model

import PyPlot.plt
using LaTeXStrings

# let's say the x axis stretches from -70 to + 70 "AU"

# Based upon the mass and inclination of the disk, caluclate length of the
# velocity vector.
# Say that max vel = 0.5 r_c, for .
# Rather, can we color in dots showing projected velocity? Shades of blue/red?
# Alternative would be to shade the halves of the disk based upon the velocity.

# We could also show the central star in the disk, as an orangy thing getting
# bigger and smaller. That's probably fine.

# If we're nifty, maybe we can make a rotation arrow along the side of the
# disk that can grow larger when the mass increases

function plot_stick(pars::Parameters, id::Int)

    # Projected arrow length
    al = pars.r_c * sind(pars.incl)

    major_axis = 2 * pars.r_c
    minor_axis = 2 * pars.r_c * cosd(pars.incl)

    fig, ax = plt.subplots(nrows=1, figsize=(3,3))

    ellipse = PyPlot.matplotlib[:patches][:Ellipse]((0.,0.), major_axis, minor_axis)
    ax[:add_artist](ellipse)
    ellipse[:set_facecolor]("Chocolate")

    ellipse = PyPlot.matplotlib[:patches][:Ellipse]((0.,0.), 12/cosd(pars.incl), 12)
    ax[:add_artist](ellipse)
    ellipse[:set_facecolor]("White")

    smidge = 0.3
    space = 3.

    # This should be an ellipse
    # ax[:plot](0., 0., "o", mfc="White", mec="Black", markersize=12, zorder=100)
    ax[:plot](0., 0., "o", mfc="Gold", mec="Black", markersize=10)

    ax[:annotate]("", xy=(smidge, al), xytext=(smidge, 0.), arrowprops={"width"=>0.08, "frac"=>0.2, "headwidth"=>4., "shrink"=>0.0, "facecolor"=>"black"})

    # ax[:annotate]("", xy=(-pars.r_c - space, -al), xytext=(-pars.r_c - space, 0.), arrowprops={"width"=>0.8, "frac"=>0.2, "headwidth"=>4., "shrink"=>0.0, "facecolor"=>"blue", "edgecolor"=>"blue"})

    # ax[:annotate]("", xy=(pars.r_c + space, al), xytext=(pars.r_c + space, 0.), arrowprops={"width"=>0.8, "frac"=>0.2, "headwidth"=>4., "shrink"=>0.0, "facecolor"=>"red", "edgecolor"=>"red"})

    width = 70.
    ax[:set_xlim](-width, width)
    ax[:set_ylim](-width, width)

    plt.savefig("stick.png")
    # plt.savefig(outdir * @sprintf("stick%04d.png", id))
end

p = pars[1]
p.incl = 0.
plot_stick(p, 1)

# How many frames per process?
# start = run_index * nframes_per_proc + 1
#
# ids = Int[i for i=start:(start + nframes_per_proc)]
#
# for id in ids
#     if id <= nframes
#         plot_stick(pars[id] , id)
#     else
#         break
#     end
# end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# port form musl library, commit e216951f509b71da193da2fc63e25b998740d58b
# see http://git.musl-libc.org/cgit/musl/tree/src/math/erf.c for implementation details
# origin: FreeBSD /usr/src/lib/msun/src/s_erf.c Copyright (c) 1993: Sun Microsystems, Inc.

let
global erf
global erfc

const erx  = 8.45062911510467529297e-01 # 0x3FEB0AC1, 0x60000000

# Coefficients for approximation to  erf on [0,0.84375]
const efx8 =  1.02703333676410069053e+00 # 0x3FF06EBA, 0x8214DB69
const pp0  =  1.28379167095512558561e-01 # 0x3FC06EBA, 0x8214DB68
const pp1  = -3.25042107247001499370e-01 # 0xBFD4CD7D, 0x691CB913
const pp2  = -2.84817495755985104766e-02 # 0xBF9D2A51, 0xDBD7194F
const pp3  = -5.77027029648944159157e-03 # 0xBF77A291, 0x236668E4
const pp4  = -2.37630166566501626084e-05 # 0xBEF8EAD6, 0x120016AC
const qq1  =  3.97917223959155352819e-01 # 0x3FD97779, 0xCDDADC09
const qq2  =  6.50222499887672944485e-02 # 0x3FB0A54C, 0x5536CEBA
const qq3  =  5.08130628187576562776e-03 # 0x3F74D022, 0xC4D36B0F
const qq4  =  1.32494738004321644526e-04 # 0x3F215DC9, 0x221C1A10
const qq5  = -3.96022827877536812320e-06 # 0xBED09C43, 0x42A26120

# Coefficients for approximation to  erf  in [0.84375,1.25]
const pa0  = -2.36211856075265944077e-03 # 0xBF6359B8, 0xBEF77538
const pa1  =  4.14856118683748331666e-01 # 0x3FDA8D00, 0xAD92B34D
const pa2  = -3.72207876035701323847e-01 # 0xBFD7D240, 0xFBB8C3F1
const pa3  =  3.18346619901161753674e-01 # 0x3FD45FCA, 0x805120E4
const pa4  = -1.10894694282396677476e-01 # 0xBFBC6398, 0x3D3E28EC
const pa5  =  3.54783043256182359371e-02 # 0x3FA22A36, 0x599795EB
const pa6  = -2.16637559486879084300e-03 # 0xBF61BF38, 0x0A96073F
const qa1  =  1.06420880400844228286e-01 # 0x3FBB3E66, 0x18EEE323
const qa2  =  5.40397917702171048937e-01 # 0x3FE14AF0, 0x92EB6F33
const qa3  =  7.18286544141962662868e-02 # 0x3FB2635C, 0xD99FE9A7
const qa4  =  1.26171219808761642112e-01 # 0x3FC02660, 0xE763351F
const qa5  =  1.36370839120290507362e-02 # 0x3F8BEDC2, 0x6B51DD1C
const qa6  =  1.19844998467991074170e-02 # 0x3F888B54, 0x5735151D

# Coefficients for approximation to  erfc in [1.25,1/0.35]
const ra0  = -9.86494403484714822705e-03 # 0xBF843412, 0x600D6435
const ra1  = -6.93858572707181764372e-01 # 0xBFE63416, 0xE4BA7360
const ra2  = -1.05586262253232909814e+01 # 0xC0251E04, 0x41B0E726
const ra3  = -6.23753324503260060396e+01 # 0xC04F300A, 0xE4CBA38D
const ra4  = -1.62396669462573470355e+02 # 0xC0644CB1, 0x84282266
const ra5  = -1.84605092906711035994e+02 # 0xC067135C, 0xEBCCABB2
const ra6  = -8.12874355063065934246e+01 # 0xC0545265, 0x57E4D2F2
const ra7  = -9.81432934416914548592e+00 # 0xC023A0EF, 0xC69AC25C
const sa1  =  1.96512716674392571292e+01 # 0x4033A6B9, 0xBD707687
const sa2  =  1.37657754143519042600e+02 # 0x4061350C, 0x526AE721
const sa3  =  4.34565877475229228821e+02 # 0x407B290D, 0xD58A1A71
const sa4  =  6.45387271733267880336e+02 # 0x40842B19, 0x21EC2868
const sa5  =  4.29008140027567833386e+02 # 0x407AD021, 0x57700314
const sa6  =  1.08635005541779435134e+02 # 0x405B28A3, 0xEE48AE2C
const sa7  =  6.57024977031928170135e+00 # 0x401A47EF, 0x8E484A93
const sa8  = -6.04244152148580987438e-02 # 0xBFAEEFF2, 0xEE749A62

# Coefficients for approximation to  erfc in [1/.35,28]
const rb0  = -9.86494292470009928597e-03 # 0xBF843412, 0x39E86F4A
const rb1  = -7.99283237680523006574e-01 # 0xBFE993BA, 0x70C285DE
const rb2  = -1.77579549177547519889e+01 # 0xC031C209, 0x555F995A
const rb3  = -1.60636384855821916062e+02 # 0xC064145D, 0x43C5ED98
const rb4  = -6.37566443368389627722e+02 # 0xC083EC88, 0x1375F228
const rb5  = -1.02509513161107724954e+03 # 0xC0900461, 0x6A2E5992
const rb6  = -4.83519191608651397019e+02 # 0xC07E384E, 0x9BDC383F
const sb1  =  3.03380607434824582924e+01 # 0x403E568B, 0x261D5190
const sb2  =  3.25792512996573918826e+02 # 0x40745CAE, 0x221B9F0A
const sb3  =  1.53672958608443695994e+03 # 0x409802EB, 0x189D5118
const sb4  =  3.19985821950859553908e+03 # 0x40A8FFB7, 0x688C246A
const sb5  =  2.55305040643316442583e+03 # 0x40A3F219, 0xCEDF3BE6
const sb6  =  4.74528541206955367215e+02 # 0x407DA874, 0xE79FE763
const sb7  = -2.24409524465858183362e+01 # 0xC03670E2, 0x42712D62

function erfc1{T<:FloatTypes}(x::T)
    s = abs(x) - 1
    P = @horner_oftype s pa0 pa1 pa2 pa3 pa4 pa5 pa6
    Q = @horner_oftype s 1.0 qa1 qa2 qa3 qa4 qa5 qa6
    return 1 - T(erx) - P/Q
end

function erfc2{T<:FloatTypes}(ix::UInt32, x::T)
    if ix < highword(T(1.25))
        # 0.84375 <= |x| < 1.25
        return erfc1(x)
    end
    # 1.25 <= |x| < 28
    x = abs(x)
    s = 1/(x*x)
    if ix < highword(T(1/0.35000001))
        # 1.25 <= |x| < 1/.35 ~ 2.85714
        R = @horner_oftype s ra0 ra1 ra2 ra3 ra4 ra5 ra6 ra7
        S = @horner_oftype s 1.0 sa1 sa2 sa3 sa4 sa5 sa6 sa7 sa8
    else
        # 1/.35 <= |x| < 28
        R = @horner_oftype s rb0 rb1 rb2 rb3 rb4 rb5 rb6
        S = @horner_oftype s 1.0 sb1 sb2 sb3 sb4 sb5 sb6 sb7
    end
    z = trunclo(x)
    return exp(-z*z-T(0.5625))*exp((z-x)*(z+x)+R/S)/x
end

function erf{T<:FloatTypes}(x::T)
    ix = highword(x)
    sign = (ix>>31) % Bool
    ix &= 0x7fffffff
    if ix >= highword(T(Inf)) # erf(nan)=nan, erf(+-inf)=+-1
        return 1-2*sign + 1/x
    end
    if ix < highword(T(0.84375)) # |x| < 0.84375
        if ix < highword(T(2)^-28) #|x| < 2**-28  avoid underflow
            return (8*x +T(efx8)*x)/8
        end
        z = x*x
        r = @horner_oftype z pp0 pp1 pp2 pp3 pp4
        s = @horner_oftype z 1 qq1 qq2 qq3 qq4 qq5
        y = r/s
        return x + x*y
    end
    if ix < highword(T(6)) # 0.84375 <= |x| < 6
        y = 1 - erfc2(ix,x)
    else
        y = 1 - realmin(T)
    end
    return sign ? -y : y
end

function erfc{T<:FloatTypes}(x::T)
    ix = highword(x)
    sign = (ix>>31) % Bool
    ix &= 0x7fffffff
    if ix >= highword(T(Inf)) # erfc(nan)=nan, erfc(+-inf)=0,2
        return 2*sign + 1/x
    end
    if ix < highword(T(0.84375)) # |x| < 0.84375
        if ix < highword(T(2)^-56)  # |x| < 2**-56
            return 1 - x
        end
        z = x*x
        r = @horner_oftype z pp0 pp1 pp2 pp3 pp4
        s = @horner_oftype z 1 qq1 qq2 qq3 qq4 qq5
        y = r/s
        if sign || ix < highword(T(0.25)) # x < 1/4
            return 1 - (x+x*y)
        end
        return T(0.5) - (x - T(0.5) + x*y)
    end
    if ix < highword(T(28)) # 0.84375 <= |x| < 28
        return sign ? 2 - erfc2(ix,x) : erfc2(ix,x)
    end
    return sign ? 2 - realmin(T) : realmin(T)*realmin(T)
end

end


text:
module GaussLegendre
    using FastGaussQuadrature

    function constract(integtablenum)
        return gausslegendre(integtablenum)
    end

    function gl_integ(f, x1, x2, solve_tf_val)
        xm = 0.5 * (x1 + x2)
        xr = 0.5 * (x2 - x1)

        s = sum(i -> solve_tf_val.w[i] * f(xm + xr * solve_tf_val.x[i]), eachindex(solve_tf_val.x))

        return s * xr
    end
end


text:
# ---
# title: 368. Largest Divisible Subset
# id: problem368
# author: Indigo
# date: 2021-06-14
# difficulty: Medium
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/largest-divisible-subset/description/>
# hidden: true
# ---
#
# Given a set of **distinct** positive integers, find the largest subset such
# that every pair (Si, Sj) of elements in this subset satisfies:
#
# Si % Sj = 0 or Sj % Si = 0.
#
# If there are multiple solutions, return any subset is fine.
#
# **Example 1:**
#
#
#
#     Input: [1,2,3]
#     Output: [1,2] (of course, [1,3] will also be ok)
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,4,8]
#     Output: [1,2,4,8]
#
#
#
## @lc code=start
using LeetCode

function largest_divisible_subset(nums::Vector{Int})
    sort!(nums)
    len = length(nums)
    dp = fill(1, len)
    prev = fill(0, len)
    for i in 2:len, j in 1:(i - 1)
        if nums[i] % nums[j] == 0 && dp[i] < dp[j] + 1
            dp[i] = dp[j] + 1
            prev[i] = j
        end
    end
    idx = argmax(dp)
    res =  Int[]
    while idx != 0
        push!(res, nums[idx])
        idx = prev[idx]
    end
    return res
end
## @lc code=end


text:
using NIPALS_PCA
using Test
using CSV
using DataFrames
using Pipe
# include("../src/structs.jl")

@testset "PCA" begin

    isnumcol(type) = type <: Union{Missing,Number}

    pwd()|> println
    pathof(NIPALS_PCA) |> println

    irisdir = @pipe pathof(NIPALS_PCA) |> splitpath |> _[1:end - 2] |> joinpath(_..., "test", "data", "iris")

    erange = range(1,length=4,step=2)

    x_df = CSV.File(joinpath(irisdir, "Dataset_01.txt"), header=false) |> DataFrame
    #x_df = CSV.File(joinpath(irisdir, "iris.csv"), header=false) |> DataFrame

    #x_df = x_df[:,erange]
    select!(x_df, eltype.(eachcol(x_df)) |> coltypes -> isnumcol.(coltypes) |> findall)
    #xdataset = x_df |> parseDataFrame |> dataset -> filterDataset(dataset, filters=[dataset -> dataset.mvs .< 0.25]) |> normalize
    xdataset = x_df |> parseDataFrame |> ds -> normalize(ds,doscale=true)

    pca = calcPCA(xdataset, 4)

    aeb_r2x = [0.7277,0.2303,0.0368,0.0052]

    @test isapprox(pca.r2x,aeb_r2x,atol=0.0001)
end

@testset "PLS" begin
    @test true
end


text:
"""
# Summary statistics

```julia
describe(ts::TS)
```

Compute summary statistics of `ts`. The output is a `DataFrame`
containing standard statistics along with number of missing values and
data types of columns.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;
julia> random(x) = rand(MersenneTwister(123), x...);
julia> ts = TS(random(([1, 2, 3, 4, missing], 10)))
julia> describe(ts)
2×7 DataFrame
 Row │ variable  mean     min    median   max    nmissing  eltype
     │ Symbol    Float64  Int64  Float64  Int64  Int64     Type
─────┼───────────────────────────────────────────────────────────────────────────
   1 │ Index        5.5       1      5.5     10         0  Int64
   2 │ x1           2.75      2      3.0      4         2  Union{Missing, Int64}

```
"""
function describe(io::IO, ts::TS)
    DataFrames.describe(ts.coredata)
end
TSx.describe(ts::TS) = TSx.describe(stdout, ts)

function Base.show(io::IO, ts::TS)
    println("(", TSx.nrow(ts), " x ", TSx.ncol(ts), ") TS with ", eltype(index(ts)), " Index")
    println("")
    DataFrames.show(ts.coredata, show_row_number=false, summary=false)
    return nothing
end
Base.show(ts::TS) = show(stdout, ts)

function Base.summary(io::IO, ts::TS)
    println("(", nr(ts), " x ", nc(ts), ") TS")
end

"""
# Size methods

```julia
nrow(ts::TS)
nr(ts::TS)
```

Return the number of rows of `ts`. `nr` is an alias for `nrow`.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> ts = TS(collect(1:10))
julia> TSx.nrow(ts)
10
```
"""
function nrow(ts::TS)
    DataFrames.size(ts.coredata)[1]
end
# alias
nr = TSx.nrow

function Base.lastindex(ts::TS)
    lastindex(index(ts))
end

function Base.length(ts::TS)
    TSx.nrow(ts)
end

# Number of columns
"""
# Size methods

```julia
ncol(ts::TS)
```

Return the number of columns of `ts`. `nc` is an alias for `ncol`.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;

julia> random(x) = rand(MersenneTwister(123), x);

julia> TSx.ncol(TS([random(100) random(100) random(100)]))
3

julia> nc(TS([random(100) random(100) random(100)]))
3
```
"""
function ncol(ts::TS)
    DataFrames.size(ts.coredata)[2] - 1
end
# alias
nc = TSx.ncol

# Size of
"""
# Size methods
```julia
size(ts::TS)
```

Return the number of rows and columns of `ts` as a tuple.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> TSx.size(TS([collect(1:100) collect(1:100) collect(1:100)]))
(100, 3)
```
"""
function size(ts::TS)
    nr = TSx.nrow(ts)
    nc = TSx.ncol(ts)
    (nr, nc)
end

# Return index column
"""
# Index column

```julia
index(ts::TS)
```

Return the index vector from the `coredata` DataFrame.

# Examples

```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;

julia> random(x) = rand(MersenneTwister(123), x);

julia> ts = TS(random(10), Date("2022-02-01"):Month(1):Date("2022-02-01")+Month(9));

julia> show(ts)
(10 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
───────────────────────
 2022-02-01  0.768448
 2022-03-01  0.940515
 2022-04-01  0.673959
 2022-05-01  0.395453
 2022-06-01  0.313244
 2022-07-01  0.662555
 2022-08-01  0.586022
 2022-09-01  0.0521332
 2022-10-01  0.26864
 2022-11-01  0.108871

julia> index(ts)
10-element Vector{Date}:
 2022-02-01
 2022-03-01
 2022-04-01
 2022-05-01
 2022-06-01
 2022-07-01
 2022-08-01
 2022-09-01
 2022-10-01
 2022-11-01

julia>  eltype(index(ts))
Date
```
"""
function index(ts::TS)
    ts.coredata[!, :Index]
end

"""
# Column names
```julia
names(ts::TS)
```

Return a `Vector{String}` containing column names of the TS object,
excludes index.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> names(TS([1:10 11:20]))
2-element Vector{String}:
 "x1"
 "x2"
```
"""

function names(ts::TS)
    names(ts.coredata[!, Not(:Index)])
end

"""
# First Row
```julia
first(ts::TS)
```

Returns the first row of `ts` as a TS object.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random)
julia> first(TS(1:10))
(10 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
───────────────────────
 2022-02-01  0.768448

```
"""
function Base.first(ts::TS)
    TS(Base.first(ts.coredata,1))
end

"""
# Head
```julia
head(ts::TS, n::Int = 10)
```
Returns the first `n` rows of `ts`.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random)
julia> head(TS(1:100))
(10 x 1) TS with Int64 Index

 Index  x1
 Int64  Int64
──────────────
     1      1
     2      2
     3      3
     4      4
     5      5
     6      6
     7      7
     8      8
     9      9
    10     10
```
"""
function head(ts::TS, n::Int = 10)
    TS(Base.first(ts.coredata, n))
end

"""
# Tail
```julia
tail(ts::TS, n::Int = 10)
```

Returns the last `n` rows of `ts`.

```jldoctest; setup = :(using TSx, DataFrames, Dates, Random)
julia> tail(TS(1:100))
(10 x 1) TS with Int64 Index

 Index  x1
 Int64  Int64
──────────────
    91     91
    92     92
    93     93
    94     94
    95     95
    96     96
    97     97
    98     98
    99     99
   100    100
```
"""
function tail(ts::TS, n::Int = 10)
    TS(DataFrames.last(ts.coredata, n))
end


text:
@testset "StudentT" begin
    test_interface(StudentTLikelihood(3.0, 1.5), Distributions.AffineDistribution{Float64,Continuous,<:TDist})
    test_auglik(StudentTLikelihood(3.0, 1.5))
end


text:
using Documenter, ModelingToolkit

include("pages.jl")

makedocs(sitename = "ModelingToolkit.jl",
         authors = "Chris Rackauckas",
         modules = [ModelingToolkit],
         clean = true, doctest = false,
         strict = [
             :doctest,
             :linkcheck,
             :parse_error,
             :example_block,
             # Other available options are
             # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
         ],
         format = Documenter.HTML(analytics = "UA-90474609-3",
                                  assets = ["assets/favicon.ico"],
                                  canonical = "https://mtk.sciml.ai/stable/"),
         pages = pages)

deploydocs(repo = "github.com/SciML/ModelingToolkit.jl.git";
           push_preview = true)


text:
@testset "modules" begin
    @testset "find module definition" begin
        using Atom: moduledefinition

        let (path, line) = moduledefinition(Atom)
            @test path == atommodfile
            @test line == 4
        end
        let (path, line) = moduledefinition(Junk)
            @test path == junkpath
            @test line == 1
        end
        let (path, line) = moduledefinition(Junk.SubJunk)
            @test path == subjunkspath
            @test line == 4
        end
    end

    @testset "find module files" begin
        using Atom: modulefiles

        ## Revise-like module file detection
        # works for precompiled packages
        let (parentfile, included_files) = modulefiles(Atom)
            expected = Set(atommodfiles)
            actual = Set((parentfile, included_files...))
            @test actual == expected

            # can't detect display/webio.jl
            @test_broken webiofile in modfiles
        end

        # fails for non-precompiled packages
        @test_broken junkpath == modulefiles(Junk)[1]

        ## CSTPraser-based module file detection
        let included_files = normpath.(modulefiles("Atom", atommodfile))
            # finds all the files in Atom module except display/webio.jl
            for f in atommodfiles
                f == webiofile && continue
                @test f in included_files
            end

            # only finds files in a module -- exclude files in the submodules
            @test length(atommodfiles) == length(included_files)

            # can't look for non-toplevel `include` calls
            @test_broken webiofile in included_files
        end
    end
end


text:
module PartitionedTableaus

    using ..RungeKutta

    using .Tableaus

    include("tableaus/prk.jl")

    export TableauLobattoIIIAIIIB,
           TableauLobattoIIIBIIIA,
           TableauLobattoIIICIIIC̄,
           TableauLobattoIIIC̄IIIC,
           PartitionedTableauGauss

    PartitionedTableauList = (
    )

    export PartitionedTableauList

end


text:
using TimeZones: Class, Transition

@testset "Class" begin
    @testset "construct" begin
        @test Class(:NONE) == Class(0x00)
        @test Class(:FIXED) == Class(0x01)
        @test Class(:STANDARD) == Class(0x02)
        @test Class(:LEGACY) == Class(0x04)

        @test Class(:DEFAULT) == Class(0x01 | 0x02)
        @test Class(:ALL) == Class(0x07)
    end

    @testset "getproperty field fallback" begin
        @test Class isa DataType
        @test sprint(show, Class) == "Class"
    end

    @testset "classify name/regions" begin
        @test Class("Foobar", []) == Class(:NONE)
        @test Class("UTC+1", []) == Class(:FIXED)
        @test Class("Europe/Warsaw", ["europe"]) == Class(:STANDARD)
        @test Class("US/Pacific", ["backward"]) == Class(:LEGACY)
        @test Class("Etc/GMT-14", ["etcetera"]) == Class(:LEGACY)
        @test Class("UTC", ["utc", "backward"]) == Class(:FIXED) | Class(:STANDARD) | Class(:LEGACY)
    end

    @testset "bitwise-or" begin
        @test Class(0x00) | Class(0x00) == Class(0x00)
        @test Class(0x00) | Class(0x01) == Class(0x01)
        @test Class(0x01) | Class(0x00) == Class(0x01)
        @test Class(0x01) | Class(0x01) == Class(0x01)
    end

    @testset "bitwise-and" begin
        @test Class(0x00) & Class(0x00) == Class(0x00)
        @test Class(0x00) & Class(0x01) == Class(0x00)
        @test Class(0x01) & Class(0x00) == Class(0x00)
        @test Class(0x01) & Class(0x01) == Class(0x01)
    end

    @testset "repr" begin
        @test repr(Class(:NONE)) == "Class(:NONE)"
        @test repr(Class(:FIXED)) == "Class(:FIXED)"
        @test repr(Class(:STANDARD)) == "Class(:STANDARD)"
        @test repr(Class(:LEGACY)) == "Class(:LEGACY)"

        @test repr(Class(:DEFAULT)) == "Class(:FIXED) | Class(:STANDARD)"
        @test repr(Class(:ALL)) == "Class(:FIXED) | Class(:STANDARD) | Class(:LEGACY)"

        @test repr(Class(0x08)) == "Class(0x08)"
        @test repr(Class(0x09)) == "Class(:FIXED) | Class(0x08)"
    end
end


text:
using YaoEGraph
using Test

using Metatheory, Metatheory.EGraphs
using YaoHIR
using YaoLocations
using YaoHIR: X

@metatheory_init

naive = @theory begin
    Chain(Gate(X, locs), Gate(X, locs)) |> Chain([])
end

circ = Chain(Gate(X, Locations(1)), Gate(X, Locations(1)))
egraph = EGraph(circ)
saturate!(egraph, naive)
extract!(egraph, astsize)


text:
@testset "Inner products" begin
    rₘₐₓ = 300
    ρ = rand()
    N = ceil(Int, rₘₐₓ/ρ + 1/2)

    for B in [StaggeredFiniteDifferences(N, ρ), FiniteDifferences(N, 1.0)]
        S = B'B
        for T in [Float64,ComplexF64]
            vv = rand(T, size(B,2))
            v = B*vv

            normalize!(v)

            @test norm(v) ≈ 1.0
            @test *(v'.args..., v.args...) ≈ 1.0
            @test v'v ≈ 1.0
            @test dot(vv,S,vv) ≈ 1.0

            # lazyip = Dot(v, v)
            # @test materialize(lazyip) ≈ 1.0
        end
    end

    @testset "Restricted bases" begin
        N = 100
        n = 10
        dr = 0.1

        R = StaggeredFiniteDifferences(N, dr)
        R̃ = R[:,1:n]

        f = 3

        ϕ = R*(f*ones(N))
        ϕ̃ = R̃*(f*ones(n))

        @test norm(ϕ) ≈ f*√(N*dr)
        @test norm(ϕ̃) ≈ f*√(n*dr)

        @test ϕ'ϕ ≈ (f^2*N*dr)
        @test ϕ̃'ϕ̃ ≈ (f^2*n*dr)

        @test apply(*, ϕ', ϕ) ≈ (f^2*N*dr)
        @test apply(*, ϕ̃', ϕ̃) ≈ (f^2*n*dr)

        normalize!(ϕ)
        @test norm(ϕ) ≈ 1.0
        normalize!(ϕ̃)
        @test norm(ϕ̃) ≈ 1.0
    end
end


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

abstract type AbstractDatabase{D<:AbstractDict} end
abstract type AbstractSubDatabase{D<:AbstractDict} <: AbstractDatabase{D} end

sources(database::AbstractDatabase) = database.sources
sources!(database::AbstractDatabase, value) = database.sources = value

default_kwargs(database::AbstractDatabase) = database.kwargs
default_kwargs!(database::AbstractDatabase, value) = database.kwargs = value

computing_functions(database::AbstractDatabase) = database.computing_functions
computing_functions!(database::AbstractDatabase, value) =
    database.computing_functions = value

datadict(database::AbstractDatabase) = database.datadict
datadict!(database::AbstractDatabase, value) = database.datadict = value

modeldict(database::AbstractDatabase) = database.modeldict
modeldict!(database::AbstractDatabase, value) = database.modeldict = value

name(database::AbstractDatabase) = database.name
name!(database::AbstractDatabase, value) = database.name = value

Base.iterate(database::AbstractDatabase) = iterate(datadict(database))
Base.keys(database::AbstractDatabase) = keys(datadict(database))
Base.values(database::AbstractDatabase) = values(datadict(database))
Base.pairs(database::AbstractDatabase) = pairs(datadict(database))

function Base.summary(database::AbstractDatabase{D}) where {D}
    data::D = datadict(database)
    len = length(data)
    nsubdata = count(v -> v isa D, values(data))
    ntables = len - nsubdata
    type = typeof(database).name
    "$(len)-element $(type) [$(ntables) table(s), $(nsubdata) subdatabase(s)]"
end

mutable struct Database{D<:AbstractDict} <: AbstractDatabase{D}
    sources::Vector{Symbol}
    kwargs::Dict{Symbol,Any}
    computing_functions::Vector{Function}
    datadict::D
    modeldict::D
    paramdict::D
    name::String
    function Database{D}(
        sources::Vector{Symbol},
        kwargs::Dict{Symbol,Any},
        computing_functions::Vector{Function},
        datadict::Union{D,Nothing},
        name::String = "Database",
    ) where {D<:AbstractDict}
        data = something(datadict, D())
        modeldict, paramdict = D(), D()
        db = new(sources, kwargs, computing_functions, data, modeldict, paramdict, name)
        if isnothing(datadict)
            data = import_data!(db; kwargs...)
            eltdata = length(data) == 1 ? first(data)[2] : data
            eltdata isa D && (db.datadict = data)
        end
        db
    end
end

root(database::Database) = database
Base.parent(database::Database) = database

function Database{D}(
    sources,
    kwargs,
    computing_functions;
    name::String = "Database",
) where {D<:AbstractDict}
    sources = collect(Symbol.(sources))
    kwargs, cmp = Dict{Symbol,Any}(kwargs), Vector{Function}(computing_functions)
    Database{D}(sources, kwargs, cmp, nothing, name)
end

Database(args...) = Database{Dict{Symbol,Any}}(args...)

function Database{D}(
    sources::Union{AbstractVector,Tuple},
    computing_functions::Vector{Function} = Function[];
    name::String = "Database",
    kwargs...,
) where {D<:AbstractDict}
    sources = collect(Symbol.(sources))
    Database{D}(sources, Dict{Symbol,Any}(kwargs), computing_functions; name = name)
end

function Database{D}(
    source::Union{Symbol,AbstractString},
    args...;
    name::String = "Database",
    kwargs...,
) where {D}
    Database{D}(Vector{Symbol}([Symbol(source)]), args...; name = name, kwargs...)
end

function Base.show(io::IO, database::AbstractDatabase{D}) where {D}
    print(io, summary(database), "\n")
end

function Base.show(io::IO, ::MIME"text/plain", database::AbstractDatabase{D}) where {D}
    show(io, database)
end

function Base.show(database::AbstractDatabase{D}) where {D}
    show(stdout, database)
end

function create_column!(
    df::AbstractDataFrame,
    colname::Symbol,
    columnsyms,
    func::Function;
    kwargs...,
)
    dfcols = names(df)
    if hasproperty(df, colname) || any(col -> !hasproperty(df, col), columnsyms)
        (df, false)
    else
        ind = ncol(df) + 1
        columns = [df[!, col] for col ∈ columnsyms]
        newcol = func(columns...)
        df = insertcols!(df, ind, colname => newcol)
        (df, true)
    end
end

function create_column!(
    data::AbstractDict,
    colname::Symbol,
    columnsyms,
    func::Function;
    kwargs...,
)
    updated = false
    for (_, elt) ∈ data
        (_, elt_updt) = create_column!(elt, colname, columnsyms, func; kwargs...)
        updated = updated || elt_updt
    end
    (data, updated)
end

function create_column!(
    database::AbstractDatabase,
    colname::Symbol,
    columnsyms,
    func::Function;
    kwargs...,
)
    data = datadict(database)
    (_, updated) = create_column!(data, colname, columnsyms, func; kwargs...)
    (database, updated)
end

function create_table!(
    data::AbstractDict,
    tblname::Symbol,
    tablesyms,
    func::Function;
    kwargs...,
)
    updated = false
    for (key, elt) ∈ data
        if elt isa AbstractDict
            (df, elt_updt) = create_table!(elt, tblname, tablesyms, func; kwargs...)
            updated = updated || elt_updt
        end
    end
    if !haskey(data, tblname) && all(tbl -> haskey(data, tbl), tablesyms)
        updated = true
        tables = [data[tbl] for tbl ∈ tablesyms]
        data[tblname] = func(tables...)
    else
        updated = false
    end
    (data, updated)
end

function create_table!(
    database::AbstractDatabase,
    tblname::Symbol,
    tablesyms,
    func::Function;
    kwargs...,
)
    (_, updated) = create_table!(datadict(database), tblname, tablesyms, func; kwargs...)
    (database, updated)
end

function group_table!(
    data::AbstractDict,
    tblname::Symbol,
    sourcetbl::Symbol,
    keycols,
    cols_and_funcs::Pair...;
    kwargs...,
)
    if !haskey(data, tblname) && haskey(data, sourcetbl)
        df = data[sourcetbl]
        keycols = Vector{Symbol}(keycols)
        combargs = []
        for (vcol, gfunc) ∈ cols_and_funcs
            append!(combargs, vcol .=> (col -> (; vcol => gfunc(col))))
        end
        groupdf = groupby(df, keycols; sort = true)
        groupdict = Dict{Symbol,Any}()
        for (gkey, df) ∈ pairs(groupdf)
            strs = map(string ∘ simplify ∘ Symbol, gkey) |> collect
            key = Symbol(join(strs, "__"))
            key = ensure_unique(groupdict, key)
            groupdict[key] = DataFrame(df)
        end
        combdf = combine(groupdf) do subdf
            colpairs = []
            for (vcols, gfunc) ∈ cols_and_funcs, vcol ∈ vcols
                push!(colpairs, vcol => gfunc(subdf[!, vcol]))
            end
            DataFrame(; colpairs...)
        end
        tblname_groups = ensure_unique(data, Symbol(tblname, :_group))
        data[tblname_groups] = groupdict
        data[tblname] = combdf
        updated = true
    else
        updated = false
    end
    for (_, elt) ∈ data
        if elt isa AbstractDict
            (df, elt_updt) =
                group_table!(elt, tblname, sourcetbl, keycols, cols_and_funcs...; kwargs...)
            updated = updated || elt_updt
        end
    end
    (data, updated)
end

function group_table!(
    database::AbstractDatabase,
    tblname::Symbol,
    sourcetbl::Symbol,
    keycols,
    cols_and_funcs::Pair...;
    kwargs...,
)
    (_, updated) = group_table!(
        datadict(database),
        tblname,
        sourcetbl,
        keycols,
        cols_and_funcs...;
        kwargs...,
    )
    (database, updated)
end

function group_table!(data, tblname, sourcetbls, keycols, cols_and_funcs...; kwargs...)
    updated = false
    for sourcetbl ∈ sourcetbls
        (_, updt) =
            group_table!(data, tblname, sourcetbl, keycols, cols_and_funcs...; kwargs...)
        updated = updated || updt
    end
    (data, updated)
end

macro defcolumn(fcall::Expr, body)
    @argcheck fcall.head == :call
    colname = fcall.args[1]
    funname = Symbol(:column_, colname)
    columnnames = fcall.args[2:end]
    args = Tuple(columnnames)
    colsym = QuoteNode(colname)
    quote
        function $(esc(funname))(data)
            func = ($(esc.(args)...),) -> $(esc(body))
            create_column!(data, $(colsym), $(args), func)[2]
        end
    end
end

macro deftable(fcall::Expr, body)
    @argcheck fcall.head == :call
    tblname = fcall.args[1]
    funname = Symbol(:table_, tblname)
    tablenames = fcall.args[2:end]
    args = Tuple(tablenames)
    tblsym = QuoteNode(tblname)
    quote
        function $(esc(funname))(data)
            func = ($(esc.(args)...),) -> $(esc(body))
            create_table!(data, $(tblsym), $(args), func)[2]
        end
    end
end

macro defgroup(fcall::Expr, keynames::Expr, cols_and_funcs::Expr)
    @argcheck fcall.head == :call
    @argcheck keynames.head ∈ (:tuple, :vect)
    @argcheck cols_and_funcs.head ∈ (:tuple, :vect, :block)

    tblname = fcall.args[1]
    funname = Symbol(:group_, tblname)
    tablesyms = fcall.args[2:end]
    args = Tuple(tablesyms)
    tblsym = QuoteNode(tblname)
    keysyms = [k isa QuoteNode ? k.value : k for k ∈ keynames.args]

    cols_funcs = []
    for pair ∈ cols_and_funcs.args
        if !(pair isa LineNumberNode)
            @argcheck pair isa Expr && pair.head == :call && pair.args[1] == :(=>)
            (cols, func) = filter(elt -> !(elt isa LineNumberNode), pair.args[2:end])
            @argcheck cols isa Expr && cols.head ∈ (:tuple, :vect)

            push!(cols_funcs, :($(Tuple(cols.args)) => $(esc(func))))
        end
    end

    quote
        function $(esc(funname))(data)
            group_table!(data, $(tblsym), $(args), $keysyms, $(cols_funcs...))[2]
        end
    end
end


text:
# functions related to F distribution

import .RFunctions:
    fdistpdf,
    fdistlogpdf,
    fdistcdf,
    fdistccdf,
    fdistlogcdf,
    fdistlogccdf,
    fdistinvcdf,
    fdistinvccdf,
    fdistinvlogcdf,
    fdistinvlogccdf


text:
push!(LOAD_PATH, "../")

using Documenter, JavisNB

makedocs(sitename = "JavisNB.jl")

makedocs(;
    modules = [JavisNB],
    authors = "Giovanni Puccetti <g.puccetti92@gmail.com>, Jacob Zelko <jacobszelko@gmail.com> and contributors",
    repo = "https://github.com/JuliaAnimators/JavisNB.jl/blob/{commit}{path}#L{line}",
    sitename = "JavisNB.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://JuliaAnimators.github.io/JavisNB.jl",
        assets = String[],
    ),
    pages = ["Home" => "index.md", "References" => "references.md"],
)

deploydocs(;
    repo = "github.com/JuliaAnimators/JavisNB.jl",
    push_preview = true,
    devbranch = "main",
)


text:
function remove_zero_cols(M::Array{Float64,2})
    ## find out the indices of the zero columns of Mi
    ind_zero_column = []
    for j = 1:size(M,2)
        if M[:,j]==zeros(size(M,1))
            push!(ind_zero_column, j)
        end
    end
    ind_zero_column
    println("Number of zero columns: $(length(ind_zero_column))")
    return M[:, setdiff(1:size(M,2), ind_zero_column)]
end

function remove_zero_rows(M::Array{Float64,2})
    ## find out the indices of the zero columns of Mi
    ind_zero_rows = []
    for i = 1:size(M,1)
        if M[i,:]==zeros(size(M,2))
            push!(ind_zero_rows, i)
        end
    end
    ind_zero_rows
    println("Number of zero rows: $(length(ind_zero_rows))")
    return M[setdiff(1:size(M,1), ind_zero_rows), :]
end

using Simplicial
#############################################################
## Idea: use normalized speed to construct a filtration

function FiltSparseMatrix(M::Array{Float64,2}, cutoff_size::Int)
    Ord_dic = Dict()
    for i = 1:size(M,1)
        ind_nz_i = findall(M[i,:].!=0)
        Ord_i = Array{Int,1}(sortslices([M[i,ind_nz_i] ind_nz_i], dims = 1)[:,2])
        Ord_dic[i] = Ord_i
    end
    Ord_dic

    Rho = zeros(Int,size(M,1))
    for i = 1:size(M,1)
        Rho[i] = length(Ord_dic[i])
    end

    ind_rem = [i for i = 1:size(M,1) if Rho[i]>=cutoff_size] ## rem = remain
    println("number of rows remained: ", length(ind_rem))
    M_rem = M[ind_rem, :];

    ## Construct the corresponding order matrices.
    ## Note: the indices are reset to 1:length(ind_rem)
    Ord_dic_rem = Dict()
    for i = 1:length(ind_rem)
        Ord_dic_rem[i] = Ord_dic[ind_rem[i]]
    end
    Ord_dic_rem

    Rho_rem = BigInt[length(Ord_dic_rem[i]) for i = 1:length(ind_rem)];

    ## find out the filtration indices
    ind_filt = BigInt[]
    ind_top = lcm(Rho_rem)

    for i = 1:length(ind_rem)
        temp = Array{BigInt, 1}(Array{BigInt, 1}(collect(1:Rho_rem[i])).*ind_top./Rho_rem[i])
        ind_filt = BigInt[ind_filt; temp]
    end
    ind_filt = sort(collect(Set(ind_filt)));

    ## To efficiently compute the filtration,
    ## compute the dictionary
    ## that maps each filtration index to the matrix position
    ind_filt_2_pos = Dict{BigInt,Array{Tuple{Int64,Int64}, 1}}(ind=>[] for ind in ind_filt)
    for i = 1:length(ind_rem)
        diff_i = BigInt(ind_top/Rho_rem[i])
        ##println(diff_i)
        ##println(Rho_rem[i])
        for l = 1:Rho_rem[i]
            ##println("l:", l)
            ##println("diff_i*l: ", diff_i*l)
            ##println("diff_i*l: ", BigInt(diff_i)*l)
            push!(ind_filt_2_pos[diff_i*l], (i,Ord_dic_rem[i][l]))
        end
    end
    ind_filt_2_pos;

    ## Construct the filtration
    Birth = BigInt[]
    Faces = CodeWord[]
    column_2_face = Dict(j=>CodeWord() for j = 1:size(M_rem,2))
    for t = 1:length(ind_filt)
        ind_t = ind_filt[t]
        for pair in ind_filt_2_pos[ind_t]
            push!(Birth, ind_filt[t])
            push!(column_2_face[pair[2]], pair[1])
            push!(Faces, copy(column_2_face[pair[2]]))
        end
    end

    ## This block is run in case
    ## the indices in Birth is too big to be handled by FiltrationOfSimplicialComplexes
    true_birth_dic = Dict{Int, BigInt}(1=>Birth[1])
    new_birth = [1]
    i = 1
    for j = 2:length(Birth)
        if Birth[j]!=Birth[j-1]
            i+=1
            true_birth_dic[i] = Birth[j]
        end
        push!(new_birth, i)
    end

    FS = FiltrationOfSimplicialComplexes(Faces, new_birth);
    return (FS, ind_rem, true_birth_dic)
end

function FiltSparseMatrix2(M::Array{Float64,2}, cutoff_size::Int)
    Ord_dic = Dict()
    for i = 1:size(M,1)
        ind_nz_i = findall(M[i,:].!=0)
        Ord_i = Array{Int,1}(sortslices([M[i,ind_nz_i] ind_nz_i], dims = 1)[:,2])
        Ord_dic[i] = Ord_i
    end
    Ord_dic

    Rho = zeros(Int,size(M,1))
    for i = 1:size(M,1)
        Rho[i] = length(Ord_dic[i])
    end

    ind_rem = [i for i = 1:size(M,1) if Rho[i]>=cutoff_size] ## rem = remain
    println("number of rows remained: ", length(ind_rem))
    M_rem = M[ind_rem, :];

    ## Construct the corresponding order matrices.
    ## Note: the indices are reset to 1:length(ind_rem)
    Ord_dic_rem = Dict()
    for i = 1:length(ind_rem)
        Ord_dic_rem[i] = Ord_dic[ind_rem[i]]
    end
    Ord_dic_rem

    Rho_rem = BigInt[length(Ord_dic_rem[i]) for i = 1:length(ind_rem)];

    ## find out the filtration indices
    ind_filt = BigInt[]
    ind_top = lcm(Rho_rem)

    for i = 1:length(ind_rem)
        temp = Array{BigInt, 1}(Array{BigInt, 1}(collect(1:Rho_rem[i])).*ind_top./Rho_rem[i])
        ind_filt = BigInt[ind_filt; temp]
    end
    ind_filt = sort(collect(Set(ind_filt)));

    ## To efficiently compute the filtration,
    ## compute the dictionary
    ## that maps each filtration index to the matrix position
    ind_filt_2_pos = Dict{BigInt,Array{Tuple{Int64,Int64}, 1}}(ind=>[] for ind in ind_filt)
    for i = 1:length(ind_rem)
        diff_i = BigInt(ind_top/Rho_rem[i])
        ##println(diff_i)
        ##println(Rho_rem[i])
        for l = 1:Rho_rem[i]
            ##println("l:", l)
            ##println("diff_i*l: ", diff_i*l)
            ##println("diff_i*l: ", BigInt(diff_i)*l)
            push!(ind_filt_2_pos[diff_i*l], (i,Ord_dic_rem[i][l]))
        end
    end
    ind_filt_2_pos;

    ## Construct the filtration
    Birth = BigInt[]
    Faces = CodeWord[]
    column_2_face = Dict(j=>CodeWord() for j = 1:size(M_rem,2))
    for t = 1:length(ind_filt)
        ind_t = ind_filt[t]
        for pair in ind_filt_2_pos[ind_t]
            push!(Birth, ind_filt[t])
            push!(column_2_face[pair[2]], pair[1])
            push!(Faces, copy(column_2_face[pair[2]]))
        end
    end

    ## This block is run in case
    ## the indices in Birth is too big to be handled by FiltrationOfSimplicialComplexes
    true_birth_dic = Dict{Int, BigInt}(1=>Birth[1])
    new_birth = [1]
    i = 1
    for j = 2:length(Birth)
        if Birth[j]!=Birth[j-1]
            i+=1
            true_birth_dic[i] = Birth[j]
        end
        push!(new_birth, i)
    end

    FS = FiltrationOfSimplicialComplexes(Faces, new_birth);
    return (FS, ind_rem, true_birth_dic, ind_top)
end


text:
export RedisDict

struct RedisDict{K,V} <: AbstractRedisCollection
    conn::AbstractRedisConnection
    key::String

    RedisDict{K,V}(conn, key) where {K,V} = if serializeable(K) && serializeable(V)
        new(conn, key)
    elseif V == Nothing
        throw(ArgumentError("Sorry RedisDict don't support Nothing as value"))
    else
        throw(ArgumentError("RedisDict currently not supports arbitrary element type"))
    end

    RedisDict{K,V}(key) where {K,V} = RedisDict{K,V}(default_connection, key)
end

struct RedisDictKeyIterator{K,V}
    rd::RedisDict{K,V}
end

struct RedisDictValueIterator{K,V}
    rd::RedisDict{K,V}
end

function getindex(rd::RedisDict{K,V}, key) where {K,V}
    res = exec(rd.conn, "hget", rd.key, serialize(K, key))
    deserialize(V, @some(res))
end

function get(rd::RedisDict{K,V}, key, default) where {K,V}
    res = exec(rd.conn, "hget", rd.key, serialize(K, key))
    res == nothing ? default : deserialize(V, res)
end

function setindex!(rd::RedisDict{K,V}, value, key) where {K,V}
    exec(rd.conn, "hset", rd.key, serialize(K, key), serialize(V, value))
end

function keys(rd::RedisDict)
    RedisDictKeyIterator(rd)
end

function values(rd::RedisDict)
    RedisDictValueIterator(rd)
end

function keytype(::RedisDict{K,V}) where {K,V}
    K
end

function valtype(::RedisDict{K,V}) where {K,V}
    V
end

function eltype(::RedisDict{K,V}) where {K,V}
    Pair{K,V}
end

function eltype(::RedisDictKeyIterator{K,V}) where {K,V}
    K
end

function eltype(::RedisDictValueIterator{K,V}) where {K,V}
    V
end

function length(rd::RedisDict)
    exec(rd.conn, "hlen", rd.key)::Int64
end

function length(ki::RedisDictKeyIterator)
    length(ki.rd)
end

function length(vi::RedisDictValueIterator)
    length(vi.rd)
end

function iterate(rd::RedisDict)
    handle, cache = exec(rd.conn, "hscan", rd.key, 0)
    iterate(rd, (cache, 1, parse(Int, String(handle))))
end

function iterate(rd::RedisDict{K,V}, (cache, ci, handle)) where {K,V}
    ci > length(cache) && return nothing

    cache[ci], cache[ci+1], if ci+2 > length(cache) && handle != 0
        handle, cache = exec(rd.conn, "hscan", rd.key, handle)
        cache, 1, parse(Int, String(handle))
    else
        cache, ci+2, handle
    end

    Pair(deserialize(K, key), deserialize(V, val)), iter
end

"simply delegate to the Dict, will do extra deserializations for value"
function iterate(ki::RedisDictKeyIterator)
    res, state = @some iterate(ki.rd)
    car(res), state
end

function iterate(ki::RedisDictKeyIterator, state)
    res, state = @some iterate(ki.rd, state)
    car(res), state
end

"simply delegate to the Dict, will do extra deserializations for key"
function iterate(vi::RedisDictValueIterator)
    res, state = @some iterate(vi.rd)
    cadr(res), state
end

function iterate(vi::RedisDictValueIterator, iter)
    res, state = @some iterate(vi.rd, state)
    cadr(res), state
end

function collect(rd::RedisDict{K,V}) where {K,V}
    raw = exec(rd.conn, "hgetall", rd.key)
    res = Vector{Pair{K,V}}(length(raw)÷2)

    for i in 1:length(raw)÷2
        res[i] = Pair(deserialize(K, raw[2i-1]), deserialize(V, raw[2i]))
    end

    res
end

function collect(ki::RedisDictKeyIterator{K,V}) where {K,V}
    map(deserialize(K), exec(ki.rd.conn, "hkeys", ki.rd.key))
end

function in(x::Pair, rd::RedisDict)
    cadr(x) != nothing && rd[car(x)] == cadr(x)
end

function in(x, rd::RedisDict)
    throw(ArgumentError("Dicts only contain Pairs"))
end

function in(x, ki::RedisDictKeyIterator{K,V}) where {K,V}
    exec(ki.rd.conn, "hexists", ki.rd.key, serialize(K, x)) |> Bool
end


text:
using Primes
using OffsetArrays

struct primePower
    prime::Int
    power::Int
end

function memoizeany(func)
    thedict = Dict()
    return (a)-> memoizeanyaux(a, func, thedict)
end

function memoizeanyaux(a, func, thedict)
    if haskey(thedict, a)
        return thedict[a]
    end
    res = func(a)
    thedict[a] =res
    return res
end

function maxpow(n::Int, p::Int)
    pow = 1
    n /= p
    while true
        rem = n%p
        if rem != 0
            return pow
        end
        pow+=1
        n /=p
    end
end

function factorall(n)
    dictlist = [Dict{Int, Int}() for i in 1:n]
    theprimes = primes(Int(ceil(sqrt(n))))
    for i in theprimes
        for j in i:i:n
            @inbounds dictlist[j][i]=maxpow(j, i)
        end
    end
    for i in 2:n
        dl = dictlist[i]
        if isempty(dl)
            dl[i] = 1
        end
        pfacs = prod(i^dl[i] for i in keys(dl))
        therat = i/pfacs
        if therat > 1
            dl[therat]=1
        end
    end
    return dictlist
end

function factorall(m,n)
    dictlist = [Dict{Int, Int}() for i in m:n]
    dictlist = OffsetArray(dictlist, m:n)
    theprimes = primes(Int(ceil(sqrt(n))))
    for i in theprimes
        begprime = Int(ceil(m/i))*i
        for j in begprime:i:n
            @inbounds dictlist[j][i]=maxpow(j, i)
        end
    end
    thetwo = max(2, m)
    for i in thetwo:n
        dl = dictlist[i]
        if isempty(dl)
            dl[i] = 1
        end
        pfacs = prod(i^dl[i] for i in keys(dl))
        therat = i/pfacs
        if therat > 1
            dl[therat]=1
        end
    end
    return dictlist
end

function domult(ndict, func)
    if isempty(ndict)
        return 1
    end
    return prod(func(i, ndict[i]) for i in keys(ndict))
end

newfact = memoizeany(factorall)

function doallmult(n, func)
    return [domult(i, func) for i in newfact(n)]
end

function doallmult(m, n, func)
    tmpar = collect(factorall(m, n)) #compiler bug
    res = [domult(i, func) for i in tmpar]
    return OffsetArray(res, m:n)
end

function makesigmab(k)
    thedict = Dict((0, 0)=>BigInt(0))
    return (p, n) -> sigkb(p, n, k, thedict)
end

function makesigma(k)
    thedict = Dict((0, 0)=>0)
    if k==0
        return (p, n)-> n+1
    end
    return (p, n) -> sigk(p, n, k, thedict)
end

function sigkb(p, n, k,thedict)
    if haskey(thedict, (p, n))
        return thedict[(p, n)]
    end
    bp = BigInt(p)
    tmp = sum(bp^(i*k) for i in 0:n)
    thedict[(p, n)] = tmp
    return tmp
end

function sigk(p, n, k,thedict)
    if haskey(thedict, (p, n))
        return thedict[(p, n)]
    end
    @inbounds tmp = sum(p^(i*k) for i in 0:n)
    thedict[(p, n)] = tmp
    return tmp
end

function memoize( func)
    thedict = Dict((0, 0)=>0)
    return (p, n)-> memoizeaux(p, n, func, thedict)
end

function memoizeaux(p, n, func, thedict)
    if haskey(thedict, (p, n))
        return thedict[(p, n)]
    end
    res = func(p, n)
    thedict[(p, n)] = res
    return res
end

mobius(p, n) = n==1?-1:0

phifunc(p, n) = p^(n-1)*(p-1)

function tally(xs)
  freqs = Dict{eltype(xs),Int}()
  for x in xs
    freqs[x] = get(freqs, x, 0) + 1
  end
  return freqs
end

max_stars_repo_name: igorrivin/juliant

max_issues_repo_name: igorrivin/juliant

max_forks_repo_name: igorrivin/juliant


text:
"""
    Wrap(value)
    Wrap{T}()

A wrapper which makes it convenient to instantiate without unitializing memory.
Used internally by `Switch`.
"""
immutable Wrap{T}
    value::T

    Wrap(v::T) = new(v)
    Wrap(v) = new(convert(T,v))
    Wrap() = new()
end

(::Type{Wrap}){T}(v::T) = Wrap{T}(v)

Base.show{T}(io::IO, x::Wrap{T}) = print(io, "Wrap(", x.value, ")")

"""
    Switch{A,B}(value::Union{A,B})

Create an single-element storage container which can contain one of two possible
input types. Semantically this type has similarities with `Union{A,B}`, but is
optimized for faster dispatch.

Using the `Switch` instance under `broadcast` is efficient and type-stable, and
is convenient to use with "dot-call" syntax. For example

    x = Switch{Int, Float64}(1)
    x .+ 1    # = Switch{Int, Float64}(2)
    x .+ 1.0  # = 2.0 (automatically fall back to simple type where possible)
"""
immutable Switch{A,B}
    which::Bool
    a::Wrap{A}
    b::Wrap{B}

    Switch(a::A) = new(true,  Wrap{A}(a), Wrap{B}())
    Switch(b::B) = new(false, Wrap{A}(),  Wrap{B}(b))
end

Base.show{A,B}(io::IO, x::Switch{A,B}) = x.which ? print(io, "Switch{A,B}(", x.a.value, ")") : print(io, "Switch{A,B}(", x.b.value, ")")

# broadcast methods
size(::Switch) = ()

getindex(x::Switch) = x.which ? x.a.value : x.b.value # not type-stable

@pure switchtype{A}(::Type{A}, ::Type{A}) = A
@pure switchtype{A,B}(::Type{A}, ::Type{B}) = Switch{A,B}
@pure switchtype{A,B}(::Type{Switch{A,B}}, ::Type{A}) = Switch{A,B}
@pure switchtype{A,B}(::Type{Switch{A,B}}, ::Type{B}) = Switch{A,B}
@pure switchtype{A,B}(::Type{A}, ::Type{Switch{A,B}}) = Switch{A,B}
@pure switchtype{A,B}(::Type{B}, ::Type{Switch{A,B}}) = Switch{A,B}
@pure switchtype{A,B}(::Type{Switch{A,B}}, ::Type{Switch{A,B}}) = Switch{A,B}
@pure switchtype{A,B}(::Type{Switch{A,B}}, ::Type{Switch{B,A}}) = Switch{A,B}

switchtype{A,B,C}(::Type{Switch{A,B}}, ::Type{C}) = error("Switches with more than two types is not supported")
switchtype{A,B,C}(::Type{C}, ::Type{Switch{A,B}}) = error("Switches with more than two types is not supported")
@pure switchtype{A,B,C,D}(::Type{Switch{A,B}}, ::Type{Switch{C,D}}) = error("Switches with more than two types is not supported")

@inline function broadcast{A,B}(f, x::Switch{A,B})
    A2 = promote_op(f, A)
    B2 = promote_op(f, B)
    out_type = switchtype(A2,B2)

    if out_type <: Switch
        return x.which ? out_type(f(x.a.value)) : out_type(f(x.b.value))
    else
        return x.which ? f(x.a.value) : f(x.b.value)
    end
end

@inline function broadcast{A,B,C,D}(f, x::Switch{A,B}, y::Switch{C,D})
    T1 = promote_op(f, A, C)
    T2 = promote_op(f, A, D)
    T3 = promote_op(f, B, C)
    T4 = promote_op(f, B, D)
    out_type = switchtype(switchtype(switchtype(T1,T2),T3),T4)

    if out_type <: Switch
        return x.which ? (y.which ? out_type(f(x.a.value, y.a.value)) :
                                    out_type(f(x.a.value, y.b.value))) :
                         (y.which ? out_type(f(x.b.value, y.a.value)) :
                                    out_type(f(x.b.value, y.b.value)))
    else
        return x.which ? (y.which ? f(x.a.value, y.a.value) :
                                    f(x.a.value, y.b.value)) :
                         (y.which ? f(x.b.value, y.a.value) :
                                    f(x.b.value, y.b.value))
    end
end

@inline function broadcast{A,B}(f, x::Switch{A,B}, y)
    T1 = promote_op(f, A, typeof(y))
    T2 = promote_op(f, B, typeof(y))
    out_type = switchtype(T1,T2)

    if out_type <: Switch
        return x.which ? out_type(broadcast(yi -> f(x.a.value, yi), y)) : out_type(broadcast(f(x.b.value, yi), y)))
    else
        return x.which ? broadcast(yi -> f(x.a.value, yi), y) : yi -> broadcast(f(x.b.value, yi), y))
    end
end

@inline function broadcast{A,B}(f, x, y::Switch{A,B})
    T1 = promote_op(f, typeof(x), A)
    T2 = promote_op(f, typeof(x), B)
    out_type = switchtype(T1,T2)

    if out_type <: Switch
        return y.which ? out_type(broadcast(xi -> f(xi, y.a.value), x)) : out_type(broadcast(xi -> f(xi, y.b.value), x))
    else
        return y.which ? broadcast(xi -> f(xi, y.a.value), x) : broadcast(xi -> f(xi, y.b.value), x)
    end
end

@inline Base.:(.+){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(+, x, y)
@inline Base.:(.+){A,B}(x::Switch{A,B}, y) = broadcast(+, x, y)
@inline Base.:(.+){C,D}(x, y::Switch{C,D}) = broadcast(+, x, y)

@inline Base.:(.*){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(*, x, y)
@inline Base.:(.*){A,B}(x::Switch{A,B}, y) = broadcast(*, x, y)
@inline Base.:(.*){C,D}(x, y::Switch{C,D}) = broadcast(*, x, y)

@inline Base.:(.-){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(-, x, y)
@inline Base.:(.-){A,B}(x::Switch{A,B}, y) = broadcast(-, x, y)
@inline Base.:(.-){C,D}(x, y::Switch{C,D}) = broadcast(-, x, y)

@inline Base.:(./){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(/, x, y)
@inline Base.:(./){A,B}(x::Switch{A,B}, y) = broadcast(/, x, y)
@inline Base.:(./){C,D}(x, y::Switch{C,D}) = broadcast(/, x, y)

@inline Base.:(.//){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(//, x, y)
@inline Base.:(.//){A,B}(x::Switch{A,B}, y) = broadcast(//, x, y)
@inline Base.:(.//){C,D}(x, y::Switch{C,D}) = broadcast(//, x, y)

@inline Base.:(.%){A,B,C,D}(x::Switch{A,B}, y::Switch{C,D}) = broadcast(%, x, y)
@inline Base.:(.%){A,B}(x::Switch{A,B}, y) = broadcast(%, x, y)
@inline Base.:(.%){C,D}(x, y::Switch{C,D}) = broadcast(%, x, y)


text:
const PISpins = Matrix{Int16}

"""
generate a random spin configuration.
"""
rand_pispins(nsite::Int, ntau::Int) = rand([-Int16(1), Int16(1)], nsite, ntau)

"""
number of imaginary time slice.
"""
ntau(spins::PISpins) = size(spins, 2)

nsite(spins::PISpins) = size(spins, 1)

"""
    flip!(spins::PISpins, i::Int, j::Int) -> PISpins

flip a spin at specific position inplace.
"""
flip!(spins::PISpins, i::Int, j::Int) = (spins[i, j] = -spins[i, j]; spins)

"""
    PathIntegralIter{MT<:AbstractModel}
    PathIntegralIter(model::MT, spins::PISpins) where {MT<:AbstractModel} -> PathIntegralIter{MT}

An iterator for sweeps in Path Integral MonteCarlo.
"""
struct PathIntegralIter{MT<:AbstractModel}
    model::MT
    spins::PISpins
    PathIntegralIter(model::MT, spins::PISpins) where {MT<:AbstractModel} = new{MT}(model, spins)
end

function Base.iterate(si::PathIntegralIter, state::Int=1)
    model = si.model
    spins = si.spins
    lt = model.lattice
    NN = nsite(lt)
    NTAU = ntau(spins)

    beta = 1/model.temp
    dtau = beta/NTAU
    J_spatial = model.J0/NTAU
    J_Trotter = - log(tanh(dtau*model.Γ)) * model.temp / 2

    # visit each spin on the space-time lattice in sequential order
    for itau = 1:NTAU
        for ir = 1:NN
            # interaction energy with the neighbouring spins in real space and
            # the neighbouring spins in imaginary time
            exchange_field = J_spatial * sum(view(spins, neighbors(lt, ir), itau))
            # periodic boundary conditions in imaginary time
            itau_up = mod(itau, NTAU) + 1
            itau_down = mod(itau-2, NTAU) + 1
            exchange_field += J_Trotter * (spins[ir,itau_up] + spins[ir,itau_down])

            energy_diff =  2 * spins[ir,itau] * exchange_field   # Note: J_spatial > 0 corresponds fo FM interactions
            # Metropolis-Hastings update: flip the spin with probability min(1, exp(-beta*energy_diff))
            (energy_diff < 0 || rand() < exp(-beta*energy_diff)) && flip!(spins, ir, itau)
        end
    end
    ((state, model, spins), state+1)
end

"""
    measure(model::AbstractModel, spins::PISpins) -> Tuple

Measure observables, returns a tuple of (E, Mz, Mz², Mz⁴).
"""
function measure(model::AbstractModel, spins::PISpins)
    lt = model.lattice
    NN = nsite(lt)
    NTAU = ntau(spins)
    beta = 1/model.temp
    dtau = beta/NTAU

    J_spatial = model.J0/NTAU
    J_Trotter = - log(tanh(dtau*model.Γ)) * model.temp / 2
    # Measure the total energy
    energy_tot = 0.0
    magnz = 0.0
    for itau = 1:NTAU
        for ir = 1:NN
            magnz += spins[ir,itau]
            energy_tot -= J_spatial * spins[ir,itau] * sum(view(spins, neighbors(lt, ir), itau)) / 2  # compenstate for double counting of bonds
            # periodic boundary conditions in imaginary time
            itau_up = mod(itau, NTAU) + 1
            energy_tot -= J_Trotter * spins[ir,itau] * spins[ir,itau_up]
        end
    end
    magnz /= (NTAU * NN)
    energy_tot / (NTAU * NN), magnz, magnz^2, magnz^4
end

"""
    runsse(sseiter, ntherm::Int, nmeas::Int; binsize::Int=200, seed::Int=2) -> DynamicBin

Run an PathIntegralIter application for calculating mean energy and magnetization, where
    * sseiter: PathIntegralIter instance.
    * ntherm: number of steps for heat bath.
    * nmeas: number of measures.
    * binsize: size of bin.
"""
function runsse(sseiter, ntherm::Int, nmeas::Int; binsize::Int=200)
    println("thermalizing ...")
    for (k, model, spins) in sseiter
        k == ntherm && break
    end

    println("measuring ...")
    println("      E          Mz         Mz²        Mz⁴")
    bin = DynamicBin(binsize, Float64, Float64, Float64, Float64)

    for (k, model, spins) in sseiter
        res = measure(model, spins)
        push!(bin, res)

        if fitted(bin)
            print_lastbin(bin); println()
        end
        k == nmeas && break
    end
    bin
end


text:
__precompile__(true)
module ADCME

    export tf,
    AUTO_REUSE,
    GLOBAL_VARIABLES,
    TRAINABLE_VARIABLES,
    UPDATE_OPS

    using PyCall
    using Random
    using LinearAlgebra
    import Optim
    import Conda
    using SparseArrays
    using LibGit2

    ENV["KMP_DUPLICATE_LIB_OK"]="TRUE"

    tf = PyNULL()
    DTYPE = Dict{Type, PyObject}()
    # a list of custom operators
    COLIB = Dict{String, Tuple{String, String, String, Bool}}()

    libSuffix = Sys.isapple() ? "dylib" : (Sys.islinux() ? "so" : "dll")

    if isfile("$(@__DIR__)/../deps/deps.jl")
        include("$(@__DIR__)/../deps/deps.jl")
    else
        error("ADCME is not properly built; run `Pkg.build(\"ADCME\")` to fix the problem.")
    end
    run_metadata = nothing

    # global options to control ADCME behavior
    mutable struct OptionsSparse
        auto_reorder::Bool
        solver::String
        OptionsSparse() = new(true, "SparseLU")
    end
    mutable struct Options
        sparse::OptionsSparse
        Options() = new(OptionsSparse())
    end
    options = Options()

    function __init__()
        # install_custom_op_dependency() # always install dependencies
        global AUTO_REUSE, GLOBAL_VARIABLES, TRAINABLE_VARIABLES, UPDATE_OPS, DTYPE, COLIB
        copy!(tf, pyimport("tensorflow"))
        DTYPE = Dict(Float64=>tf.float64,
            Float32=>tf.float32,
            Int64=>tf.int64,
            Int32=>tf.int32,
            Bool=>tf.bool,
            ComplexF64=>tf.complex128,
            ComplexF32=>tf.complex64,
            String=>tf.string)
        AUTO_REUSE = tf.compat.v1.AUTO_REUSE
        GLOBAL_VARIABLES = tf.compat.v1.GraphKeys.GLOBAL_VARIABLES
        TRAINABLE_VARIABLES = tf.compat.v1.GraphKeys.TRAINABLE_VARIABLES
        UPDATE_OPS = tf.compat.v1.GraphKeys.UPDATE_OPS

        colibs = include("$(@__DIR__)/../deps/CustomOps/default_formulas.jl")
        for c in colibs
            push!(COLIB, c)
        end

        if isfile("$(@__DIR__)/../deps/CustomOps/formulas.jl")
            colibs = include("$(@__DIR__)/../deps/CustomOps/formulas.jl")
            for c in colibs
                push!(COLIB, c)
            end
        end
    end

    include("core.jl")
    include("io.jl")
    include("optim.jl")
    include("run.jl")
    include("variable.jl")
    include("ops.jl")
    include("layers.jl")
    include("extra.jl")
    include("sparse.jl")
    include("random.jl")
    include("gan.jl")
    include("ot.jl")
    include("ode.jl")
end


text:
using Pkg

deps = [
    Pkg.PackageSpec(url = "https://github.com/JuliaGPU/GemmKernels.jl", rev = "master"),
]
Pkg.add(deps)
Pkg.build()
Pkg.precompile()


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

using Distributions: Exponential, Geometric, cdf

using Plots
using LinearAlgebra
using ProgressBars

include("utils.jl")

# so when someone goes from S->I, we want to schedule their I->R rewrite. But how do we store that?
# right now, we just update them to R in X time steps, regardless of their current state at the time.


text:
#=
            # P is the precision used with the typed occurance
            #
type ArbFloat{P}  <: Real     # field and struct names from arb.h
  exponentOf2        ::Int    #           fmpz
  nwords_sign        ::UInt   #           mp_size_t
  significand1       ::UInt   #           significand_struct
  significand2       ::UInt   #
  radius_exponentOf2 ::Int    #           fmpz
  radius_significand ::UInt   #
end
=#

# get and set working precision for ArbFloat

const ArbFloatPrecision = [116,]

precision{P}(x::ArbFloat{P}) = P
precision{P}(::Type{ArbFloat{P}}) = P
precision(::Type{ArbFloat}) = ArbFloatPrecision[1]
# allow inquiring the precision of the module: precision(ArbFloats)
precision(::Type{Type{Val{:ArbFloats}}}) = precision(ArbFloat)
precision(m::Module) = precision(Type{Val{Symbol(m)}})

function setprecision(::Type{ArbFloat}, x::Int)
    x = max(11, abs(x))
    x > 4095 && warn("ArbFloats are designed to work best at precisions < 4096 bits")
    ArbFloatPrecision[1] = x
    return x
end

function setprecisionAugmented(::Type{ArbFloat}, x::Int, offset::Int=10)
    return setprecision(ArbFloat, x+offset)
end

# a type specific hash function helps the type to 'just work'
const hash_arbfloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4
const hash_0_arbfloat_lo = hash(zero(UInt), hash_arbfloat_lo)
# two values of the same precision
#    with identical midpoint significands and identical radus_exponentOf2s hash equal
# they are the same value, one is less accurate yet centered about the other
hash{P}(z::ArbFloat{P}, h::UInt) =
    hash(z.significand1$z.exponentOf2,
         (h $ hash(z.significand2$(~reinterpret(UInt,P)), hash_arbfloat_lo)
            $ hash_0_arbfloat_lo))

# empty constructor
ArbFloat() = ArbFloat{precision(ArbFloat)}()

# typemax,realmax realmax,realmin

typemax{P}(::Type{ArbFloat{P}}) = ArbFloat{P}("Inf")
typemin{P}(::Type{ArbFloat{P}}) = ArbFloat{P}("-Inf")
realmax{P}(::Type{ArbFloat{P}}) = ArbFloat{P}(2)^(P+29)
realmin{P}(::Type{ArbFloat{P}}) = ArbFloat{P}(2)^(-P-29)

function zero{T<:ArbFloat}(::Type{T})
   z = T()
   return z
end

#zero{P}(x::ArbFloat{P}) = zero(ArbFloat{P})
function zero{T<:ArbFloat}(x::T)
    P = precision(T)
    return zero( ArbFloat{P} )
end

function one{T<:ArbFloat}(::Type{T})
    z = T()
    ccall(@libarb(arb_set_si), Void, (Ptr{ArfFloat}, Int), &z, one(Int))
    return z
end

#one{P}(x::ArbFloat{P}) = one(ArbFloat{P})
function one{T<:ArbFloat}(x::T)
    P = precision(T)
    return one( ArbFloat{P} )
end

# parts and aspects
# midpoint, radius, lowerbound, upperbound

@inline function ptr_to_midpoint{T<:ArbFloat}(x::T) # Ptr{ArfFloat}
    return ccall(@libarb(arb_mid_ptr), Ptr{ArfFloat}, (Ptr{T}, ), &x)
end
@inline function ptr_to_radius{T<:ArbFloat}(x::T) # Ptr{ArfFloat}
    return ccall(@libarb(arb_rad_ptr), Ptr{ArfFloat}, (Ptr{T}, ), &x)
end

function midpoint{P}(x::ArbFloat{P})
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_mid_arb), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    return z
end

function radius{P}(x::ArbFloat{P})
    z = zero(ArbFloat{P})
    if !isexact(x)
        ccall(@libarb(arb_get_rad_arb), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    end
    return z
end

function diameter{T<:ArbFloat}(x::T)
    return 2.0*radius(x)
end

function upperbound{T<:ArbFloat}(x::T)
    P = precision(T)
    a = ArfFloat{P}()
    z = T()
    ccall(@libarb(arb_get_ubound_arf), Void, (Ptr{ArfFloat}, Ptr{ArbFloat}, Int), &a, &x, P)
    ccall(@libarb(arb_set_arf), Void, (Ptr{ArbFloat}, Ptr{ArfFloat}), &z, &a)
    return z
end

function lowerbound{T<:ArbFloat}(x::T)
    P = precision(T)
    a = ArfFloat{P}()
    z = T()
    ccall(@libarb(arb_get_lbound_arf), Void, (Ptr{ArfFloat}, Ptr{ArbFloat}, Int), &a, &x, P)
    ccall(@libarb(arb_set_arf), Void, (Ptr{ArbFloat}, Ptr{ArfFloat}), &z, &a)
    return z
end

bounds{T<:ArbFloat}(x::T) = ( lowerbound(x), upperbound(x) )

"""
isolate_nonnegative_content(x::ArbFloat)
returns x without any content < 0
if x is strictly < 0, returns ArbFloat's NaN
"""
function isolate_nonnegative_content{T<:ArbFloat}(x::T)
    lo, hi = bounds(x)
    z = if lo > 0
              x
          elseif hi < 0
              T(NaN)
          else
              mid = hi * 0.5
              r = Float32(mid)
              dr = eps(r) * 0.125
              m = Float64(mid)
              while Float64(r) > m
                  r = r - dr
              end
              midpoint_radius(mid, r)
          end
    return z
end

"""
isolate_positive_content(x::ArbFloat)
returns x without any content <= 0
if x is strictly <= 0, returns ArbFloats' NaN
"""
function isolate_positive_content{T<:ArbFloat}(x::T)
    lo, hi = bounds(x)
    z = if lo > 0
              x
          elseif hi <= 0
              T(NaN)
          else
              mid = hi * 0.5
              r = Float32(mid)
              dr = 0.125 * eps(r)
              m = Float64(mid)
              while Float64(r) >= m
                  r = r - dr
              end
              midpoint_radius(mid, r)
          end
    return z
end

"""
force_nonnegative_content(x::ArbFloat)
returns x without any content < 0
if x is strictly < 0, returns 0
"""
function force_nonnegative_content{T<:ArbFloat}(x::T)
    lo, hi = bounds(x)
    z = if lo >= 0
              x
          elseif hi < 0
              zero(T)
          else
              isolate_nonnegative_content(x)
          end
    return z
end

"""
force_positive_content(x::ArbFloat)
returns x without any content <= 0
if x is strictly <= 0, returns eps(lowerbound(x))
"""
function force_positive_content{T<:ArbFloat}(x::T)
    lo, hi = bounds(x)
    z = if lo > 0
              x
          elseif hi < 0
              eps(lo)
          else
              isolate_positive_content(x)
          end
    return z
end

function max{T<:ArbFloat}(x::T, y::T)
    return (x + y + abs(x - y))/2
end

function min{T<:ArbFloat}(x::T, y::T)
    return (x + y - abs(x - y))/2
end
#=
function min2{T<:ArbFloat}(x::T, y::T)
    return
        if donotoverlap(x,y)
            x < y ? x : y
        else
            xlo, xhi = bounds(x)
            ylo, yhi = bounds(y)
            lo,hi = min(xlo, ylo), min(xhi, yhi)
            md = (hi+lo)/2
            rd = (hi-lo)/2
            midpoint_radius(md, rd)
        end
end

function max2{T<:ArbFloat}(x::T, y::T)
    return
        if donotoverlap(x,y)
            return x > y ? x : y
        else
            xlo, xhi = bounds(x)
            ylo, yhi = bounds(y)
            lo,hi = max(xlo, ylo), max(xhi, yhi)
            md = (hi+lo)/2
            rd = (hi-lo)/2
            return midpoint_radius(md, rd)
        end
end
=#
#=
function max{T<:ArbFloat}(x::T, y::T)
    return ((x>=y) | !(y<x)) ? x : y
end
=#

function minmax{T<:ArbFloat}(x::T, y::T)
   return min(x,y), max(x,y) # ((x<=y) | !(y>x)) ? (x,y) : (y,x)
end

"""
Returns the effective relative error of x measured in bits,
  defined as the difference between the position of the
  top bit in the radius and the top bit in the midpoint, plus one.
  The result is clamped between plus/minus ARF_PREC_EXACT.
"""
function relativeError{T<:ArbFloat}(x::T)
    re_bits = ccall(@libarb(arb_rel_error_bits), Int, (Ptr{T},), &x)
    return re_bits
end

"""
Returns the effective relative accuracy of x measured in bits,
  equal to the negative of the return value from relativeError().
"""
function relativeAccuracy{T<:ArbFloat}(x::T)
    ra_bits = ccall(@libarb(arb_rel_accuracy_bits), Int, (Ptr{T},), &x)
    return ra_bits
end

"""
Returns the number of bits needed to represent the absolute value
  of the mantissa of the midpoint of x, i.e. the minimum precision
  sufficient to represent x exactly.
  Returns 0 if the midpoint of x is a special value.
"""
function midpointPrecision{T<:ArbFloat}(x::T)
    mp_bits = ccall(@libarb(arb_bits), Int, (Ptr{T},), &x)
    return mp_bits
end

"""
Sets y to a trimmed copy of x: rounds x to a number of bits equal
  to the accuracy of x (as indicated by its radius), plus a few
  guard bits. The resulting ball is guaranteed to contain x,
  but is more economical if x has less than full accuracy.
"""
function trimmed{T<:ArbFloat}(x::T)
    z = T()
    ccall(@libarb(arb_trim), Void, (Ptr{T}, Ptr{T}), &z, &x)
    return z
end


text:
@coeff function D(reactor::Reactor)
  ( 1 + zeta ) / ( zeta^0.25 ) / sqrt( 1 + zeta + zeta^2 / 3 )
end


text:
# Test the inf_alleles() function from src/Evolve.jl
# to run from the test folder:  julia -p 4 -L CGP.jl test_inf_alleles.jl

#include("../src/CGP.jl")
include("../src/Inf_alleles.jl")
include("../src/Propsel.jl")
using Main.CGP
MyInt = Main.CGP.MyInt

nreps = 4
numinputs = 3
numoutputs = 2
numints = 10
levsback = 10
popsize = 5:5:20
max_pop_gens = 50000
ngoals = 8   # number goals in goallist
tourn_size = 0
#println("pwd: ",pwd())
# The following lines allow this file to included from CGP.jl, CGP.jl/test, and CGP.jl/src.
cwd = pwd()
if cwd[end-5:end] == "CGP.jl"
  csvfile = "$cwd/test/data/test_inf_alleles.csv"
else
  csvfile = "../test/data/test_inf_alleles.csv"
end
#println("csvfile: ",csvfile)
#inf_alleles( p, popsize, max_pop_gens, gl, tourn_size=tourn_size )
df = run_inf_alleles(nreps,numinputs,numoutputs,numints,ngoals,levsback,popsize,max_pop_gens,tourn_size,csvfile)


text:
using Statistics
using FFTW
#FFTW.set_provider!("mkl")
#import Pkg; Pkg.add("FileIO")
using FileIO
using Base.Threads

include("aux_functions.jl")

run_name = "try_"
order = 5
const L = 5
#N = 80000
const N = 20000
const J = 50
exp_Theta = 2
exp_t = 2
θ = 10.0^(-exp_Theta)
t_i = 0.0
t_f = 10.0^(exp_t)
M = 1001
M_g = 100 + 1 #number of outputs, starting from the initial data
dt = t_f / (M-1)
t = t_i
#M = convert(Int64,t_f/dt)
const κ = 2π/L # for Fourier Transform
dx = L/J
x = [dx*(i-1) for i in 1:J] ;
p = (L, N, J, κ, dx, order)

animation = false
display = true

println("t_f = $(t_f), M = $(M), dt = $(dt), exp_Theta = $(exp_Theta)")

run_name = run_name * "_L$(L)_N2_5_J$(J)_M$(M)_o$(order)_T$(exp_Theta)"

E = zeros(J)
ϕ = zeros(J)
n = zeros(J) #charge density
S = zeros(J) #carge current
du = zeros(2*N+J); # contains r, v and E

par_dis = load("Initial_Distributions/par_dis_L5_N2_5_theta01.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N2_5_theta001.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N2_5_theta0001.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N2_5_theta00001.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta01.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta001.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta0001.jld2", "par_dis");
#par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta00001.jld2", "par_dis");

@assert length(par_dis) ÷ 2 == N

@time get_density!(par_dis, n, p)
n0 = N/L
#get_ϕ!(ϕ, -n/n0 .+ 1., κ) # chenge the sign here to make it consistent with charge conservation and the time derivative of E
@time get_ϕ!(ϕ, n/n0 .- 1., κ)
@time get_E_from_ϕ!(ϕ,E,dx)
global u = [par_dis;E];
length(u)

@time get_current!(u, S, p)

@time Coordinate_test(u[1:N],L)

println("n_total = $(sum(n .- n0))")
println("v_total = $(sum(u[N+1:2N]))")
println("E_total = $(sum(E))")

println("S_total = $(sum(S))")

if nthreads() > 1
    TS = zeros(J, nthreads())
    p_RHSC = (N, J, L, dx, order, n, S, du, get_density!, get_current_threads!, Interpolate_2, TS) ;
else
    p_RHSC = (N, J, L, dx, order, n, S, du, get_density!, get_current!, Interpolate_2) ;
end

Energy_K = zeros(M_g)
Energy_E = zeros(M_g)
E_T = zeros(M_g)
v_T = zeros(M_g)
D_T = zeros(M_g)
S_T = zeros(M_g)
T = zeros(M_g)
if animation
    par = zeros(M_g,2N)
end
Energy_K[1], Energy_E[1]  = get_energy(u,(L, N, J))
E_T[1] = sum(u[2N+1:end])
v_T[1] = sum(u[N+1:2N])
get_density!(u, n, p)
get_current!(u, S, p)
D_T[1] = sum(n)/n0/J - 1
S_T[1] = sum(S)/n0/J
T[1] = var(u[N+1:2N])
if animation
    par[1,:] = u[1:2N]
end

global t = 0.0
#global j = 1

    for k in 2:M
        RK4_Step!(RHSC,u,t,dt,p_RHSC)
        global u = [mod1.(u[1:N],L); u[N+1:end]]
        #global u = [make_periodic!(u[1:N],L); u[N+1:end]]
        #filter_constant!(u[2N+1:end])
        global t = t + dt
        if display
            if (k-1) % (M÷(M_g-1)) == 0
                j = (k-1)÷(M÷(M_g-1))+1
                #scatter(plt, u[1:N], u[N+1:2*N])
                Energy_K[j], Energy_E[j] = get_energy(u,(L, N, J))
                E_T[j] = sum(u[2N+1:end])
                v_T[j] = sum(u[N+1:2N])
                get_density!(u, n, p)
                get_current!(u, S, p)
                D_T[j] = sum(n)/n0/J - 1
                S_T[j] = sum(S)/n0/J
                T[j] = var(u[N+1:2N])
                println("j = $j , t = $t, k = $k, threads = $(nthreads()), E_K = $(Energy_K[j])")
                if animation
                    par[j,:] = u[1:2N]
                end
            end
        end
    end

n_F = zeros(J)
S_F = zeros(J)
get_density!(u, n_F, p)
get_current!(u, S_F, p)
save(run_name * "th$(nthreads())_results.jld2", Dict("p" => p, "Energy_E" => Energy_E, "Energy_K" => Energy_K, "E_f" => u[2N+1:end], "n_F" => n_F, "S_F" => S_F, "E_T"=> E_T, "v_T" => v_T, "S_T" => S_T, "D_T" => D_T, "T" => T, "u" => u))

#plot(layout=(2,2))
#plot!(subplot=1,E_T,title="Total Electric Field")
#plot!(subplot=2,v_T./N, title="Total velocity")
#plot!(subplot=3,D_T,title = "Total density")
#plot!(subplot=4,S_T,title = "Total Current")
#png("$run_name * totals")

E_F = zeros(J)
ϕ_F = zeros(J)
n_F = zeros(J)
get_density!(u, n_F, p)
n0 = N/L
get_ϕ!(ϕ_F, n_F/n0 .+ 1, κ)
#get_ϕ!(ϕ_F, n_F .- n0, κ)
get_E_from_ϕ!(ϕ_F,E_F,dx)
#plot(x,u[2N+1:end], label = "dynamical", title="Electric Field with order = $(order)")
#plot!(x,E_F,label="from constraint", ls=:dash)
#t_f = 40
#png("$run_name * E")

println("averaged total E field = $(sum(E_F))")
println("E_total = $(sum(u[2N+1:end])/J)")
println("Total velocity = $(sum(u[N+1:2N])/N)")
println("Total Charge = $(sum(n_F .- n0))")
println("Final Energy = $(get_energy(u,(L, N, J)))")
get_current!(u, S, p)
println("Total_current = $(sum(S)/J)")

#plot(x,u[2N+1:end]-E_F, label = "difference", title="Electric Field")
#png("$run_name * electric_diff")

#scatter(u[1:N],u[N+1:2N]
#, markersize = 0.3
#, title = "phase-space", legend =:false)
#png("run_name * ps")

#load(run_name * "results.jld2", "S_F")


text:
using Revise
using ZigZagBoomerang
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using Profile
const ZZB = ZigZagBoomerang
PLOT = true
PROFILE = false
n = 10 # 10-dimensional problem

# Create sparse 10-dimensional precision matrix
GS = [1=>2:10, 2=>3:5, 3=>5:10, 4=>[5], 5=>[], 6=>7:8, 7=>[], 8=>[], 9=>[10], 10=>[]] # start from conditional dependence tree
S = Matrix(1.0I, n, n)
for (i, nb) in GS, j in nb
    S[i, j] = 0.3randn()
end

Γ = sparse(S * S')

ϕ(x) = 0.5*x'*Γ*x # potential or -log(target density)

#using ForwardDiff
#∇ϕ(x) = ForwardDiff.gradient(ϕ, x)

∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # partial derivative of ϕ(x) with respect to x[i]

t0 = 0.0
x0 = rand(n)
θ0 = rand([-1.0,1.0], n)

c = [norm(Γ[:, i], 2) for i in 1:n]

Z = ZigZag(Γ, x0*0)
T = 200.0

if PROFILE
    Profile.init()
    @time trace, (tT, xT, θT), (acc, num) = @profile pdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)
    Profile.clear()
end
@time trace, (tT, xT, θT), (acc, num) = @profile pdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)

xs = last.(collect(trace))
@show acc, num, acc/num

if PLOT
    using Makie
    p1 = Makie.lines(first.(xs), getindex.(xs, 2))
    save("figures/localzigzag.png", title(p1, "ZigZag"))
    p2 = Makie.lines(getindex.(xs, 1), getindex.(xs, 2), getindex.(xs, 3))
    save("figures/localzigzag3d.png", title(p2, "ZigZag 3d"))

    xs2 = last.(collect(discretize(trace, 0.1)))
    p3 = Makie.lines(first.(xs), getindex.(xs, 2), linewidth=0.5)
    Makie.scatter!(p3, first.(xs2), getindex.(xs2, 2), markersize=0.05)
    save("figures/localzigzagdis.png", title(p3, "ZigZag discretized"))
end

λref = 1.0
x0 = randn(n)
θ0 = randn(n)
B = FactBoomerang(Γ, x0*0, λref)
@time trace, (tT, xT, θT), (acc, num) = pdmp(∇ϕ, t0, x0, θ0, T, c, B, Γ)

xs = last.(collect(trace))
@show acc, num, acc/num

if PLOT
    using Makie
    p1 = Makie.lines(first.(xs), getindex.(xs, 2))
    save("figures/factboom.png", title(p1, "Factorised Boomerang"))
    p2 = Makie.lines(getindex.(xs, 1), getindex.(xs, 2), getindex.(xs, 3))
    save("figures/factboom3d.png", title(p2, "Factorised Boomerang 3d"))

    xs2 = last.(collect(discretize(trace, 0.1)))
    p3 = Makie.lines(first.(xs), getindex.(xs, 2), linewidth=0.5)
    Makie.scatter!(p3, first.(xs2), getindex.(xs2, 2), markersize=0.05)
    save("figures/factboomdis.png", title(p3, "Factorised Boomerang discretized"))
end


text:
# T = Eltype, K = Labelcount, M = Arraydimensions
abstract type LabelEncoding{T,K,M} end
# act as scalar in broadcast, see julia #18618
Base.Broadcast.broadcastable(x::LabelEncoding) = Ref(x)

const BinaryLabelEncoding{T,M} = LabelEncoding{T,2,M}
const VectorLabelEncoding{T,K} = LabelEncoding{T,K,1}
const MatrixLabelEncoding{T,K} = LabelEncoding{T,K,2}

"""
    nlabel(obj)::Int

Returns the number of labels represented in the given object `obj`.
"""
function nlabel end

"""
    label(obj)::Vector

Returns the labels represented in the given object `obj`.
Note that the order of the labels matters.
In the case of two labels, the first element represents the positive
label and the second element the negative label.
"""
function label end

"""
    labeltype(::Type{<:LabelEncoding})

Determine the type of the labels represented by the label encoding.
"""
function labeltype end

"""
    ind2label(index, encoding)

Converts the given `index` into the corresponding label defined by the `encoding`.
Note that in the binary case, `index == 1` represents the positive label
and `index == 2` the negative label.
"""
function ind2label end

"""
    label2ind(label, encoding)::Int

Converts the given `label` into the corresponding index defined by the encoding.
Note that in the binary case, the positive label
will result in the index `1` and the negative label in the index `2`, respectively.
"""
function label2ind end

"""
    poslabel(encoding)

If the encoding is binary it will return the positive label of it.
The function will throw an error otherwise.
"""
function poslabel end

"""
    neglabel(encoding)

If the encoding is binary it will return the negative label of it.
The function will throw an error otherwise.
"""
function neglabel end

"""
    labelenc(obj) -> LabelEncoding

    Tries to determine the most approriate label-encoding to describe the
given object `obj` based on the result of `label(obj)`. Note that in
most cases this function is not typestable.

```juliarepl
julia> labelenc([:yes,:no,:no,:yes,:maybe])
MLLabelUtils.LabelEnc.NativeLabels{Symbol,3}(Symbol[:yes,:no,:maybe],Dict(:yes=>1,:maybe=>3,:no=>2))

julia> labelenc([1,0,0,1,0,1])
MLLabelUtils.LabelEnc.ZeroOne{Int64,Float64}(0.5)

julia> labelenc(Int8[1,-1,-1,1,-1,1])
MLLabelUtils.LabelEnc.MarginBased{Int8}()
```
"""
function labelenc end

"""
    islabelenc(obj, encoding)

Checks is the given object `obj` can be described as being produced
by the given `encoding` in which case the function returns true, or false otherwise.
"""
function islabelenc end

"""
    isposlabel(x, encoding)

Checks if the given value `x` can be interpreted as the positive label
given the `encoding`. This function takes potential classification rules into account.
"""
function isposlabel end

"""
    isneglabel(x, encoding)

Checks if the given value `x` can be interpreted as the negative label
given the `encoding`. This function takes potential classification rules into account.
"""
function isneglabel end

"""
    classify(x, encoding)

Returns the classified version of `x` given the `encoding`.
Which means that if `x` can be interpreted as a positive label,
the positive label of `encoding` is returned; the negative otherwise.
"""
function classify end

"""
    classify!(out, x, encoding)

Same as `classify`, but uses `out` to store the result.
"""
function classify! end

"""
    convertlabel(new_encoding, x, [old_encoding])
Converts the given value/array `x` from the `old_encoding` into the
`new_encoding`. Note that if `old_encoding` is not specified it will
be derived automaticaly using `labelenc`.
    julia> convertlabel(LabelEnc.MarginBased, [0, 1, 1, 0, 0])
    5-element Array{Int64,1}:
     -1
      1
      1
     -1
     -1
    julia> convertlabel([:yes,:no], [0, 1, 1, 0, 0])
    5-element Array{Symbol,1}:
     :no
     :yes
     :yes
     :no
     :no
For more information on the available encodings, see `?LabelEnc`.
    convertlabel(new_encoding, x, [old_encoding], [obsdim])
When working with `OneOfK` one can additionally specifify which
dimension of the array denotes the observations using `obsdim`
    julia> convertlabel(LabelEnc.OneOfK, [0, 1, 1, 0, 0], obsdim = 2)
    2×5 Array{Int64,2}:
     0  1  1  0  0
     1  0  0  1  1
"""
function convertlabel end
function convertlabel! end

"""
    convertlabel(new_encoding, vec::AbstractVector, [old_encoding]) -> (Readonly)MappedArray

Creates a lazy view into `vec` that makes it look like it is
in the encoding specified by `new_encoding`, while it is actually
preserved as being of `old_encoding`.
This method only works for label-encodings that are vector-based
(i.e. pretty much all but `OneOfK`). The resulting MappedArray
will be writeable unless `old_encoding` is of type `OneVsRest`,
in which case the result will be a `ReadonlyMappedArray`.
"""
function convertlabelview end

"""
    labelmap(obj) -> Dict

Computes a mapping from the labels in `obj` to all the individual
element-indices in `obj` that correspond to that label
    julia> labelmap([0, 1, 1, 0, 0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5]
      1 => [2,3]
"""
function labelmap end

"""
    labelmap!(dict, idx, elem) -> Dict

Updates the given label-map `dict` with the new element `elem`,
which is assumed to be associated with the index `idx`.
    julia> lm = labelmap([0, 1, 1, 0, 0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5]
      1 => [2,3]
    julia> labelmap!(lm, 6, 0)
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5,6]
      1 => [2,3]
    julia> labelmap!(lm, 7:8, [1,0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5,6,8]
      1 => [2,3,7]
"""
function labelmap! end

"""
    labelfreq(obj) -> Dict

Computes the absolute frequencies for each label in `obj`.
    julia> labelfreq([0, 1, 1, 0, 0])
    Dict{Int64,Int64} with 2 entries:
      0 => 3
      1 => 2
"""
function labelfreq end

"""
    labelfreq!(dict, obj) -> Dict

Updates the given label-frequency-map `dict` with the absolute
frequencies for each label in `obj`
    julia> ld = labelfreq([0, 1, 1, 0, 0])
    Dict{Int64,Int64} with 2 entries:
      0 => 3
      1 => 2
    julia> labelfreq!(ld, [1,0,0])
    Dict{Int64,Int64} with 2 entries:
      0 => 5
      1 => 3
"""
function labelfreq! end

"""
    labelmap2vec(dict) -> Vector

Inverse function of labelmap.
Computes an `array` of labels by element-wise
traversal of the entries in `dict`.
    julia> labelvec = [:yes,:no,:no,:yes,:yes]
    julia> lm = labelmap(labelvec)
    Dict{Symbol,Array{Int64,1}} with 2 entries:
        :yes => [1, 4, 5]
        :no  => [2, 3]
    julia> labelmap2vec(lm)
    5-element Array{Symbol,1}:
        :yes
        :no
        :no
        :yes
        :yes
"""
function labelmap2vec end


text:
using LinearAlgebra, Jets, JetPack, Test

n1,n2 = 33,44

@testset "JopShift, correctness T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    b = rand(T,1)
    F = JopShift(JetSpace(T,n1,n2),b)
    x1 = rand(domain(F))
    @test F*x1 == x1 .+ b

    b1 = rand(T,n1,n2)
    F = JopShift(b1)
    @test size(domain(F)) == (n1,n2)
    @test size(range(F)) == (n1,n2)

    @test F*x1 == x1 .+ b1
end

@testset "JopShift, linearity test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    @info "Testing scalar b"
    b = rand(T,1)
    F = JopShift(JetSpace(T,n1,n2),b)
    J = jacobian(F, rand(domain(F)))
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs

    @info "Testing vector b"
    b = rand(T,n1,n2)
    F = JopShift(b)
    J = jacobian(F, rand(domain(F)))
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs
end

@testset "JopShift, dot product test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    @info "Testing scalar b"
    b = rand(T,1)
    F = JopShift(JetSpace(T,n1,n2),b)
    J  = jacobian!(F, rand(domain(F)))
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(J)), -1 .+ 2 .* rand(range(J)))
    @test lhs ≈ rhs

    @info "Testing vector b"
    b = rand(T,n1,n2)
    F = JopShift(b)
     J  = jacobian!(F, rand(domain(F)))
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(J)), -1 .+ 2 .* rand(range(J)))
    @test lhs ≈ rhs
end

#Linearization test modified since the linearization of an affine operator is exactly linear.
#Usual linearization test will fail since ||F(m + dm) - F(m) - Jdm|| will be close to zero and observed quantity will have 0/0
@testset "JopShift, linearization test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    @info "Testing scalar b"
    b = rand(T,1)
    F = JopShift(JetSpace(T,n1,n2),b)
    m0 = rand(domain(F))
    μ  = sqrt.([1/1,1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256,1/512,1/1024,1/2048,1/4096,1/8192])
    δm = rand(domain(F))
    J = jacobian!(F,m0)
    phi = ones(size(μ))
    for i=1:length(μ)
        dnon = F * (m0 .+ μ[i].*δm)
        dlin = F*m0 .+ (μ[i] .* (J*δm))
        phi[i] = norm(dnon .- dlin)
    end
    @show phi
    @test minimum(phi) <= 100*eps(eltype(abs.(b[1])))

    @info "Testing vector b"
    b = rand(T,n1,n2)
    F = JopShift(b)
    m0 = rand(domain(F))
    μ  = sqrt.([1/1,1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256,1/512,1/1024,1/2048,1/4096,1/8192])
    δm = rand(domain(F))
    J = jacobian!(F,m0)
    phi = ones(size(μ))
    for i=1:length(μ)
        dnon = F * (m0 .+ μ[i].*δm)
        dlin = F*m0 .+ (μ[i] .* (J*δm))
        phi[i] = norm(dnon .- dlin)
    end
    @show phi
    @test minimum(phi) <= 100*eps(eltype(abs.(b[1])))
end


text:
using Printf

function exception_examples()
    println("--------------------------------------")
    println("Exception examples")

    println()

    print("Enter a number ")
    num1 = chomp(readline())
    print("Enter a number ")
    num2 = chomp(readline())

    try
        val = (parse(Int32, num1)) / (parse(Int32, num2))
        if (val == Inf)
            error("Can't divide by zero")
        else
            println(val)
        end
    catch e
        println(e)
    end
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using DifferentialEquations
using Distributions

function getsolution(sol, times)
  #function to get species abundance from ODE solution
  if VERSION < v"0.6-"
    x1 = map(x -> x[1], sol(times))
    y1 = map(x -> x[2], sol(times))
  else
    x1 = map(x -> x[1], sol(times).u)
    y1 = map(x -> x[2], sol(times).u) #Need .u syntax as of 0.6
  end
  return x1, y1
end

srand(1)
println("Test ABC SMC with ODE model")

f = @ode_def LV begin
  dx = a*x - b*x*y
  dy = b*x*y - y
end a=>1.0 b=>1.0

inputa = 0.8
inputb = 2.2
g = LV(a = inputa, b = inputb)
x0 = [1.0; 0.5]
tspan = (0.0, 15.0)
prob = ODEProblem(g, x0, tspan)
sol = solve(prob)

#generate target data
times = 1.0:1.0:15.0
x, y = getsolution(sol, times)
x += rand(Normal(0.0, 0.25^2), length(times))
y += rand(Normal(0.0, 0.25^2), length(times))
targetdata = [x, y]

function simLV(params, constants, targetdata)
  a = params[1]
  b = params[2]
  x0 = [1.0; 0.5]
  tspan = (0.0, 15.0)
  times = 1.0:1.0:15.0
  h = LV(a = a, b = b)
  prob = ODEProblem(h, x0, tspan)
  sol = solve(prob)
  x1, y1 = getsolution(sol, times)
  d = sum((x1 .- targetdata[1]).^2 + (y1 .- targetdata[2]).^2)
  return d, sol
end

setup = ABCSMC(simLV,
  2,
  0.1,
  Prior([Uniform(0.0, 5.0), Uniform(0.0, 5.0)]),
  maxiterations = 10^6,
  convergence = 0.05,
  nparticles = 500
  )
@time ressmc = runabc(setup, targetdata, verbose = false, progress = false);

println("\t Check parameters are inferred correctly for Lotka Volterra model")
@test isapprox(mean(ressmc.parameters, weights(ressmc.weights), 1)[1], inputa, rtol = 0.05)
@test isapprox(mean(ressmc.parameters, weights(ressmc.weights), 1)[2], inputb, rtol = 0.05)


text:
using Memento
export m_debug, m_info, m_notice, m_warn, m_error

# NOTE: It isn't generally recommended to configure your logger at the package/library level.
push!(getlogger(Mocha),
      DefaultHandler(STDOUT, DefaultFormatter("[{date} | {level} | {name}]: {msg}")))
setlevel!(getlogger(Mocha), "info")
setpropagating!(getlogger(Mocha), false)

function m_debug(msg :: AbstractString...)
  debug(getlogger(Mocha), prod(msg))
end

function m_info(msg :: AbstractString...)
  info(getlogger(Mocha), prod(msg))
end

function m_notice(msg :: AbstractString...)
  notice(getlogger(Mocha), prod(msg))
end

function m_warn(msg :: AbstractString...)
  warn(getlogger(Mocha), prod(msg))
end

function m_error(msg :: AbstractString...)
  error(getlogger(Mocha), prod(msg))
end