text:
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:AbstractGray,T2<:AbstractGray} = Gray{promote_type(eltype(T1), eltype(T2))}
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:AbstractRGB,T2<:AbstractRGB} = RGB{promote_type(eltype(T1), eltype(T2))}

# no-op and element-type conversions, plus conversion to and from transparency
# Colorimetry conversions are in Colors.jl
convert(::Type{C}, c::C) where {C<:Colorant} = c
convert(::Type{C}, c) where {C<:Colorant} = cconvert(ccolor(C, typeof(c)), c)
cconvert(::Type{C}, c::C) where {C} = c
cconvert(::Type{C}, c) where {C}    = _convert(C, base_color_type(C), base_color_type(c), c)
convert(::Type{C}, c::Color, alpha) where {C<:TransparentColor} = cconvert(ccolor(C, typeof(c)), c, alpha)
cconvert(::Type{AlphaColor{C,T,N}}, c::C, alpha) where {C<:Color,T,N} = alphacolor(C)(c, alpha)
cconvert(::Type{ColorAlpha{C,T,N}}, c::C, alpha) where {C<:Color,T,N} = coloralpha(C)(c, alpha)
cconvert(::Type{C}, c::Color, alpha) where {C<:TransparentColor} =_convert(C, base_color_type(C), base_color_type(c), c, alpha)

# Fallback definitions that print nice error messages
_convert(::Type{C}, ::Any, ::Any, c) where {C} = error("No conversion of ", c, " to ", C, " has been defined")
_convert(::Type{C}, C1::Any, C2::Any, c, alpha) where {C} = error("No conversion of (", c, ",alpha=$alpha) to ", C, " with consistency-types $C1 and $C2 has been defined")

# Any AbstractRGB types can be interconverted
# (the first 2 are just for ambiguity resolution)
# Note: on julia 0.3 these have to be before the block below, or you
# get a spurious ambiguity warning.
_convert(::Type{Cout}, ::Type{C1}, ::Type{C1}, c) where {Cout<:AbstractRGB,C1<:AbstractRGB} = Cout(red(c), green(c), blue(c))
_convert(::Type{A}, ::Type{C1}, ::Type{C1}, c) where {A<:TransparentRGB,C1<:AbstractRGB} = A(red(c), green(c), blue(c), alpha(c))
_convert(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) where {Cout<:AbstractRGB,C1<:AbstractRGB,C2<:AbstractRGB} = Cout(red(c), green(c), blue(c))
_convert(::Type{A}, ::Type{C1}, ::Type{C2}, c) where {A<:TransparentRGB,C1<:AbstractRGB,C2<:AbstractRGB} = A(red(c), green(c), blue(c), alpha(c))

# Implementations for when the base color type is not changing
# These might trip/add transparency, however
_convert(::Type{Cout}, ::Type{Ccmp}, ::Type{Ccmp}, c) where {Cout<:Color3,Ccmp<:Color3} = Cout(comp1(c), comp2(c), comp3(c))
_convert(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c) where {A<:Transparent3,Ccmp<:Color3} = A(comp1(c), comp2(c), comp3(c), alpha(c))
_convert(::Type{Cout}, ::Type{Ccmp}, ::Type{Ccmp}, c) where {Cout<:AbstractGray,Ccmp<:AbstractGray} = Cout(gray(c))
_convert(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c) where {A<:TransparentGray,Ccmp<:AbstractGray} = A(gray(c), alpha(c))

# With user-supplied alpha
_convert(::Type{A}, ::Type{Ccmp}, ::Type{Ccmp}, c, alpha) where {A<:Transparent3,Ccmp<:Color3} = A(comp1(c), comp2(c), comp3(c), alpha)

# Grayscale
_convert(::Type{Cout}, ::Type{C1}, ::Type{C2}, c) where {Cout<:AbstractGray,C1<:AbstractGray,C2<:AbstractGray} = Cout(gray(c))
_convert(::Type{A}, ::Type{C1}, ::Type{C2}, c) where {A<:TransparentGray,C1<:AbstractGray,C2<:AbstractGray} = A(gray(c), alpha(c))
_convert(::Type{A}, ::Type{C1}, ::Type{C2}, c, alpha) where {A<:TransparentGray,C1<:AbstractGray,C2<:AbstractGray} = A(gray(c), alpha)

convert(::Type{RGB24},   x::Real) = RGB24(x, x, x)
convert(::Type{ARGB32},  x::Real) = ARGB32(x, x, x, 1)
convert(::Type{Gray24},  x::Real) = Gray24(x)
convert(::Type{AGray32}, x::Real) = AGray32(x, 1)
convert(::Type{AGray32}, x::Real, alpha) = AGray32(x, alpha)

convert(::Type{Gray{T}},  x::Real) where {T}    = Gray{T}(x)
convert(::Type{T},  x::Gray) where {T<:Real}    = convert(T, x.val)
convert(::Type{T},  x::Gray24) where {T<:Real}  = convert(T, gray(x))
convert(::Type{AGray{T}}, x::Real) where {T}    = AGray{T}(x)
convert(::Type{GrayA{T}}, x::Real) where {T}    = GrayA{T}(x)

(::Type{T})(x::AbstractGray) where {T<:Real}    = T(gray(x))

# Define some constructors that just call convert since the fallback constructor in Base
# is removed in Julia 0.7
# The parametric types are handled in @make_constructors and @make_alpha
for t in (:ARGB32, :Gray24, :RGB24)
    @eval $t(x) = convert($t, x)
end

# Generate the transparent analog of a color
alphacolor(c::C) where {C<:Color} = alphacolor(C)(c)
alphacolor(c::C,a) where {C<:Color} = alphacolor(C)(c,a)
alphacolor(c::C) where {C<:TransparentColor} = alphacolor(base_color_type(C))(color(c), alpha(c))
alphacolor(c::C,a) where {C<:TransparentColor} = alphacolor(base_color_type(C))(color(c), a)
coloralpha(c::C) where {C<:Color} = coloralpha(C)(c)
coloralpha(c::C,a) where {C<:Color} = coloralpha(C)(c,a)
coloralpha(c::C) where {C<:TransparentColor} = coloralpha(base_color_type(C))(color(c), alpha(c))
coloralpha(c::C,a) where {C<:TransparentColor} = coloralpha(base_color_type(C))(color(c), a)


text:
@enum(Event,
    L1_DCM = 0x80000000,  # Level 1 data cache misses
    L1_ICM,  # Level 1 instruction cache misses
    L2_DCM,  # Level 2 data cache misses
    L2_ICM,  # Level 2 instruction cache misses
    L3_DCM,  # Level 3 data cache misses
    L3_ICM,  # Level 3 instruction cache misses
    L1_TCM,  # Level 1 total cache misses
    L2_TCM,  # Level 2 total cache misses
    L3_TCM,  # Level 3 total cache misses
    CA_SNP,  # Snoops
    CA_SHR,  # Request for shared cache line (SMP)
    CA_CLN,  # Request for clean cache line (SMP)
    CA_INV,  # Request for cache line Invalidation (SMP)
    CA_ITV,  # Request for cache line Intervention (SMP)
    L3_LDM,  # Level 3 load misses
    L3_STM,  # Level 3 store misses
    BRU_IDL, # Cycles branch units are idle
    FXU_IDL, # Cycles integer units are idle
    FPU_IDL, # Cycles floating point units are idle
    LSU_IDL, # Cycles load/store units are idle
    TLB_DM,  # Data translation lookaside buffer misses
    TLB_IM,  # Instr translation lookaside buffer misses
    TLB_TL,  # Total translation lookaside buffer misses
    L1_LDM,  # Level 1 load misses
    L1_STM,  # Level 1 store misses
    L2_LDM,  # Level 2 load misses
    L2_STM,  # Level 2 store misses
    BTAC_M,  # BTAC miss
    PRF_DM,  # Prefetch data instruction caused a miss
    L3_DCH,  # Level 3 Data Cache Hit
    TLB_SD,  # Xlation loadookaside buffer shootdowns (SMP)
    CSR_FAL, # Failed store conditional instructions
    CSR_SUC, # Successful store conditional instructions
    CSR_TOT, # Total store conditional instructions
    MEM_SCY, # Cycles Stalled Waiting for Memory Access
    MEM_RCY, # Cycles Stalled Waiting for Memory Read
    MEM_WCY, # Cycles Stalled Waiting for Memory Write
    STL_ICY, # Cycles with No Instruction Issue
    FUL_ICY, # Cycles with Maximum Instruction Issue
    STL_CCY, # Cycles with No Instruction Completion
    FUL_CCY, # Cycles with Maximum Instruction Completion
    HW_INT,  # Hardware interrupts
    BR_UCN,  # Unconditional branch instructions executed
    BR_CN,   # Conditional branch instructions executed
    BR_TKN,  # Conditional branch instructions taken
    BR_NTK,  # Conditional branch instructionsructions not taken
    BR_MSP,  # Conditional branch instructionsructions mispred
    BR_PRC,  # Conditional branch instructions corr. pred
    FMA_INS, # FMA instructions completed
    TOT_IIS, # Total instructions issued
    TOT_INS, # Total instructions executed
    INT_INS, # Integer instructions executed executed
    FP_INS,  # Floating point instructions executed
    LD_INS,  # Load instructions executed
    SR_INS,  # Store instructions executed
    BR_INS,  # Total branch instructions executed
    VEC_INS, # Vector/SIMD instructions executed (could include integer)
    RES_STL, # Cycles processor is stalled on resource
    FP_STAL, # Cycles any FP units are stalled
    TOT_CYC, # Total cycles executed
    LST_INS, # Total load/store inst. executed
    SYC_INS, # Sync. inst. executed
    L1_DCH,  # L1 D Cache Hit
    L2_DCH,  # L2 D Cache Hit
    L1_DCA,  # L1 D Cache Access
    L2_DCA,  # L2 D Cache Access
    L3_DCA,  # L3 D Cache Access
    L1_DCR,  # L1 D Cache Read
    L2_DCR,  # L2 D Cache Read
    L3_DCR,  # L3 D Cachee Read
    L1_DCW,  # L1 D Cache Write
    L2_DCW,  # L2 D Cache Write
    L3_DCW,  # L3 D Cache Write
    L1_ICH,  # L1 instruction cache hits
    L2_ICH,  # L2 instruction cache hitss
    L3_ICH,  # L3 instruction cache hits
    L1_ICA,  # L1 instruction cache accesses
    L2_ICA,  # L2 instruction cache accesses
    L3_ICA,  # L3 instruction cache accesses
    L1_ICR,  # L1 instruction cache reads
    L2_ICR,  # L2 instruction cache reads
    L3_ICR,  # L3 instruction cache reads
    L1_ICW,  # L1 instruction cache writes
    L2_ICW,  # L2 instruction cache writes
    L3_ICW,  # L3 instruction cache writes
    L1_TCH,  # L1 total cache hits
    L2_TCH,  # L2 total cache hits
    L3_TCH,  # L3 total cache hits
    L1_TCA,  # L1 total cache accesses
    L2_TCA,  # L2 total cache accesses
    L3_TCA,  # L3 total cache accesses
    L1_TCR,  # L1 total cache reads
    L2_TCR,  # L2 total cache reads
    L3_TCR,  # L3 total cache reads
    L1_TCW,  # L1 total cachee writes
    L2_TCW,  # L2 total cache writes
    L3_TCW,  # L3 total cache writes
    FML_INS, # FM ins
    FAD_INS, # FA ins
    FDV_INS, # FD ins
    FSQ_INS, # FSq ins
    FNV_INS, # Finv ins
    FP_OPS,  # Floating point operations executed
    SP_OPS,  # Floating point operations executed; optimized to count scaled single precision vector operations
    DP_OPS,  # Floating point operations executed; optimized to count scaled double precision vector operations
    VEC_SP,  # Single precision vector/SIMD instructions
    VEC_DP,  # Double precision vector/SIMD instructions
    REF_CYC, # Reference clock cycles
    END      # This should always be last!
)

ext: jl


text:
module Percolation

export clustgyration, hklabel, hkmaxclust!, genclust, clustfractal

using DataStructures

"""
Finds the radius of gyration for the largest closed cluster in the percolation grid

Percolation is assumed to occur from top to bottom. Uses the Hoshen-Kopelman algorithm.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

The gyration radius of the largest closed cluster (the largest cluster that doesn't connect
the top and bottom of the grid)
"""
function clustgyration(grid::BitMatrix)
    labels, list, clusters = hklabel(grid)
    if length(clusters) == 0
        return 0
    end

    maxlabel = hkmaxclust!(labels, list, clusters)
    if maxlabel == 0
        return 0
    end

    maxclust = Tuple.(findall(x -> x != 0 && hkfind(list, x) == maxlabel, labels))
    mass = length(maxclust)
    if mass == 0
        return 0
    end

    centerofmass = reduce(.+, maxclust) ./ mass
    return √(sum(pos -> reduce(.+, (pos .- centerofmass).^2),
        maxclust) / mass)
end

"""
Uses the Hoshen-Kopelman algorithm to label the grid

returns the labels, label mapping list, and cluster lengths
"""
@inline function hklabel(grid::BitMatrix)
    height, width = size(grid)
    labels = zeros(Integer, height, width)
    list = Vector{Integer}()
    clusters = Vector{Integer}()

    # first row, special case because of top boundary
    hkfirstrow!(grid, labels, list, clusters)

    for row in 2:height
        # first column, special case because of left boundary
        if grid[row, 1]
            hkfirstcol!(labels, list, clusters, row)
        end

        # "middle" points, with appropriate neighbors
        for col in 2:width
            if grid[row, col]
                hkiter!(labels, list, clusters, row, col)
            end
        end
    end

    return labels, list, clusters
end

"""
Helper function for `hkpercolate`; Adds new label.
"""
@inline function hknew!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    labels[row, col] = length(list) + 1
    push!(list, labels[row, col])
    push!(clusters, 0)
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first row.
"""
@inline function hkfirstrow!(grid::BitMatrix, labels::Matrix{<:Integer},
        list::Vector{<:Integer}, clusters::Vector{<:Integer})
    width = size(labels)[2]
    if grid[1, 1]
        hknew!(labels, list, clusters, 1, 1)
        clusters[labels[1, 1]] += 1
    end
    for col in 2:width
        if grid[1, col]
            if labels[1, col - 1] != 0
                labels[1, col] = hkfind!(list, labels[1, col - 1])
            else
                hknew!(labels, list, clusters, 1, col)
            end
            clusters[labels[1, col]] += 1
        end
    end
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first column.
"""
@inline function hkfirstcol!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer)
    if labels[row - 1, 1] != 0
        labels[row, 1] = hkfind!(list, labels[row - 1, 1])
    else
        hknew!(labels, list, clusters, row, 1)
    end
    clusters[labels[row, 1]] += 1
end

"""
Helper function for `hkpercolate`; Applies one iteration of the Hoshen-Kopelman algorithm.
"""
@inline function hkiter!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    above, left = labels[row - 1, col], labels[row, col - 1]
    if above == 0
        if left == 0
            hknew!(labels, list, clusters, row, col)
        else
            labels[row, col] = hkfind!(list, left)
        end
    elseif left == 0
        labels[row, col] = hkfind!(list, above)

    else
        labels[row, col] = hkunion!(list, clusters, above, left)
    end
    clusters[labels[row, col]] += 1
end

"""
Union function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkunion!(list::Vector{<:Integer}, clusters::Vector{<:Integer},
        x::Integer, y::Integer)
    root = hkfind!(list, x)
    dest = hkfind!(list, y)
    if dest != root
        list[dest] = root
        clusters[root] += clusters[dest]
        clusters[dest] = 0
    end
    return root
end

"""
Find function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkfind!(list::Vector{<:Integer}, x::Integer)
    y = x;
    while list[y] != y
        y = list[y];
    end
    while list[x] != x
        z = list[x];
        list[x] = y;
        x = z;
    end
    return y;
end

"""
Same as hkfind!, but doesn't update labels
"""
@inline function hkfind(list::Vector{<:Integer}, x::Integer)
    while x != list[x]
        x = list[x];
    end
    return x
end

"""
Finds the label of the largest closed cluster (a.k.a. finite cluster).

Also, sets the size of open clusters (a.k.a. infinite cluster) to zero, hence the !
"""
@inline function hkmaxclust!(labels::Matrix{Integer}, list::Vector{Integer},
        clusters::Vector{Integer})
    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0
                    toplabel = hkfind(list, root)
                    bottomlabel = hkfind(list, dest)
                    if toplabel == bottomlabel
                        clusters[toplabel] = 0
                    end
                end
            end
        end
    end

    if maximum(clusters) > 0
        return argmax(clusters)
    else
        return 0
    end
end

"""
Uses a depth-first search algorithm to generate a cluster with the given site probability

# Arguments

- `probability`: the probability a site is generated when growing the cluster
- `maxsize::Integer` (optional, defaults to 10000): the maximum size of the cluster
    (breaks the loop after the cluster is grown to this size)

# Returns

an integer matrix (8-bit) with 0 where no site exists, 1 where a cell is generated, and 2
where the growth is stopped
"""
function genclust(probability::Float64; maxsize::Integer=10000, crop::Bool=true)
    cluster = zeros(Int8, maxsize, maxsize)
    # create "fences" to stop the cluster form reaching the bounds of the matrix
    cluster[1, :] .= 2
    cluster[:, 1] .= 2
    cluster[end, :] .= 2
    cluster[:, end] .= 2

    stack = Stack{Tuple{Integer, Integer}}()
    seed = (maxsize÷2, maxsize÷2)
    push!(stack, seed)
    cluster[seed...] = 1

    clustdfs!(stack, cluster, probability, maxsize)

    if crop
        return cropclust!(cluster)
    else
        return cluster
    end
end

"""
The depth-first search algorithm used for `genclust`
"""
@inline function clustdfs!(stack::Stack, cluster::Matrix{<:Integer},
        probability::Float64, maxsize::Integer)
    size = 1
    while !isempty(stack) && size < maxsize
        site = pop!(stack)
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = site .+ step
            if cluster[neighbor...] == 0
                if rand() <= probability
                    cluster[neighbor...] = 1
                    push!(stack, neighbor)
                    size += 1
                else
                    cluster[neighbor...] = 2
                end
            end
        end
    end
end

"""
finds the mass and radius of gyration of the cluster from the output of `genclust`
"""
function clustfractal(cluster::Matrix{Int8})
    clust = Tuple.(findall(x -> x == 1, cluster))
    mass = length(clust)
    if mass == 0
        return 0
    end

    centerofmass = reduce(.+, clust) ./ mass
    return mass, √(sum(pos -> reduce(.+, (pos .- centerofmass).^2), clust) / mass)
end

"""
crop the `genclust` function output to remove excess zeros
"""
@inline function cropclust!(cluster::Matrix{<:Integer})
    height, width = size(cluster)
    top = bottom = height ÷ 2
    left = right = width ÷ 2

    while 1 in cluster[top + 1, 2:end-1] || 2 in cluster[top + 1, 2:end-1]
        top += 1
    end
    while 1 in cluster[bottom - 1, 2:end-1] || 2 in cluster[bottom - 1, 2:end-1]
        bottom -= 1
    end
    while 1 in cluster[2:end-1, left - 1] || 2 in cluster[2:end-1, left - 1]
        left -= 1
    end
    while 1 in cluster[2:end-1, right + 1] || 2 in cluster[2:end-1, right + 1]
        right += 1
    end

    return cluster[bottom:top, left:right]
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
using AQFED, Test
using StatsBase
using AQFED.Black
using AQFED.Collocation
import Polynomials: coeffs
using AQFED.Math
#using Plots

@testset "FilterConvex" begin
    strikes =
        Float64.([
            20,
            25,
            50,
            55,
            75,
            100,
            120,
            125,
            140,
            150,
            160,
            175,
            180,
            195,
            200,
            210,
            230,
            240,
            250,
            255,
            260,
            270,
            275,
            280,
            285,
            290,
            300,
            310,
            315,
            320,
            325,
            330,
            335,
            340,
            350,
            360,
            370,
            380,
            390,
            400,
            410,
            420,
            430,
            440,
            450,
            460,
            470,
            480,
            490,
            500,
            510,
            520,
            550,
            580,
            590,
            600,
            650,
            670,
            680,
            690,
            700,
        ])
    vols = [
        1.21744983334323,
        1.1529735541872308,
        1.0013512993166844,
        1.0087013871410198,
        0.9055919576135943,
        0.8196499269009432,
        0.779704840770866,
        0.753927847741657,
        0.7255349986293694,
        0.7036962946028743,
        0.6870907597202961,
        0.6631489459500445,
        0.6542809839143336,
        0.6310048431894977,
        0.6231979513191437,
        0.6154526014300009,
        0.5866214834144697,
        0.5783751483731193,
        0.5625036590124762,
        0.5625539176150428,
        0.5572331684618123,
        0.5485212417739607,
        0.5456131657256524,
        0.540060895711996,
        0.5384776792271245,
        0.5325298112839504,
        0.5222410647552144,
        0.5202396738775005,
        0.5168414254536685,
        0.5127405490209541,
        0.5100440087558921,
        0.50711984442627,
        0.5042896073005682,
        0.5013959697030379,
        0.4961897259221575,
        0.4914478237113829,
        0.48571052433313705,
        0.4820982302575811,
        0.4776551485043659,
        0.4682253137830999,
        0.46912624306506934,
        0.4652049749994563,
        0.4621036693145566,
        0.45969798571592985,
        0.4561356005182957,
        0.45418189139835186,
        0.4515451651258398,
        0.44541885580442636,
        0.4452833907060621,
        0.44303755690672525,
        0.43939212779385645,
        0.4413175310749832,
        0.4336322023390991,
        0.4297053821023934,
        0.4284357423754355,
        0.4241077476619805,
        0.4222672729031064,
        0.4203436892852212,
        0.4193419518701644,
        0.41934732346075626,
        0.41758929420417745,
    ]
    forward = 356.73063159822254
    tte = 1.5917808219178082
    refPricesf = [
        337.9478782712897,
        333.2264008151271,
        310.74145637852655,
        306.2719712239312,
        288.3940506055495,
        266.219801810884,
        249.06397096223955,
        244.7750182500784,
        232.02358238485036,
        223.57092985061377,
        215.48527697171835,
        203.47284496330764,
        199.46870596050405,
        187.64870164975872,
        183.87579062925917,
        176.5074646794032,
        161.77083277969126,
        154.64429335245714,
        147.51776392522302,
        144.31801826019864,
        141.1182775951742,
        134.7188062651254,
        131.6082360685426,
        128.49767087195968,
        125.40239375195168,
        122.30712163194369,
        116.11658739192775,
        110.79059339934288,
        108.12760140305048,
        105.46461440675806,
        102.80163241046563,
        100.15226884869381,
        97.60223566948574,
        95.07799161470712,
        90.23581650058495,
        85.6198931055724,
        81.10708395817021,
        76.78697167463605,
        72.46686939110188,
        68.14677710756771,
        64.77986823082749,
        61.41296935408727,
        58.04608047734704,
        54.95662415902851,
        51.90375987138524,
        49.029996842756624,
        46.15624381412802,
        43.2825007854994,
        41.01535154006894,
        38.748212294638485,
        36.481083049208024,
        34.612533924958775,
        29.006916552211,
        24.412399354175687,
        22.88090362149725,
        21.349417888818817,
        16.324515763803383,
        14.545066178175531,
        13.74560498753116,
        12.997715379078853,
        12.249835770626547,
    ]
    w1 = ones(length(strikes))
    prices, weights = Collocation.weightedPrices(true, strikes, vols, w1, forward, 1.0, tte)
    strikes, pricesf = Collocation.filterConvexPrices(strikes, prices, w1, forward, tol = 1e-6)
    for (k, p, pf, rpf) in zip(strikes, prices, pricesf, refPricesf)
        println(k, " ", p, " ", pf, " ", pf - p)
        @test isapprox(rpf, pf, atol = 1e-5)
    end

    isoc, m = Collocation.makeIsotonicCollocation(strikes, pricesf, weights, tte, forward, 1.0, deg = 5)
    sol = Collocation.Polynomial(isoc)
    println("Solution ", sol, " ", coeffs(sol), " ", Collocation.stats(sol), " measure ", m)
    refCoeffs = [
        348.59726724651006,
        221.54300193865407,
        -24.275731591637737,
        -12.821670197699545,
        10.80303198111674,
        2.940110896855002,
    ]
    for (ciRef, ci) in zip(refCoeffs, coeffs(sol))
        @test isapprox(ciRef, ci, atol = 1e-2)
    end
end

@testset "BadRoots" begin
    strikes =
        Float64.([
            20,
            25,
            50,
            55,
            75,
            100,
            120,
            125,
            140,
            150,
            160,
            175,
            180,
            195,
            200,
            210,
            230,
            240,
            250,
            255,
            260,
            270,
            275,
            280,
            285,
            290,
            300,
            310,
            315,
            320,
            325,
            330,
            335,
            340,
            350,
            360,
            370,
            380,
            390,
            400,
            410,
            420,
            430,
            440,
            450,
            460,
            470,
            480,
            490,
            500,
            510,
            520,
            550,
            580,
            590,
            600,
            650,
            670,
            680,
            690,
            700,
        ])
    vols = [
        1.21744983334323,
        1.1529735541872308,
        1.0013512993166844,
        1.0087013871410198,
        0.9055919576135943,
        0.8196499269009432,
        0.779704840770866,
        0.753927847741657,
        0.7255349986293694,
        0.7036962946028743,
        0.6870907597202961,
        0.6631489459500445,
        0.6542809839143336,
        0.6310048431894977,
        0.6231979513191437,
        0.6154526014300009,
        0.5866214834144697,
        0.5783751483731193,
        0.5625036590124762,
        0.5625539176150428,
        0.5572331684618123,
        0.5485212417739607,
        0.5456131657256524,
        0.540060895711996,
        0.5384776792271245,
        0.5325298112839504,
        0.5222410647552144,
        0.5202396738775005,
        0.5168414254536685,
        0.5127405490209541,
        0.5100440087558921,
        0.50711984442627,
        0.5042896073005682,
        0.5013959697030379,
        0.4961897259221575,
        0.4914478237113829,
        0.48571052433313705,
        0.4820982302575811,
        0.4776551485043659,
        0.4682253137830999,
        0.46912624306506934,
        0.4652049749994563,
        0.4621036693145566,
        0.45969798571592985,
        0.4561356005182957,
        0.45418189139835186,
        0.4515451651258398,
        0.44541885580442636,
        0.4452833907060621,
        0.44303755690672525,
        0.43939212779385645,
        0.4413175310749832,
        0.4336322023390991,
        0.4297053821023934,
        0.4284357423754355,
        0.4241077476619805,
        0.4222672729031064,
        0.4203436892852212,
        0.4193419518701644,
        0.41934732346075626,
        0.41758929420417745,
    ]
    forward = 356.73063159822254
    tte = 1.5917808219178082
    w1 = ones(length(strikes))
    prices, weights = Collocation.weightedPrices(true, strikes, vols, w1, forward, 1.0, tte)
    isoc, m = Collocation.makeIsotonicCollocation(strikes, prices, weights, tte, forward, 1.0, deg = 11, degGuess = 1)
    sol = Collocation.Polynomial(isoc)
    k = collect(800:0.1:1000)
    #Halley/SuperHalley does not behave well between 900 and 950
    #LithBoonkkampIJzerman does not bethave well between 825 and 830 (requires ~ 64*3 iterations)
    ivk = @. Black.impliedVolatility(true, Collocation.priceEuropean(sol, true, k, forward, 1.0), forward, k, tte, 1.0)
    for iv in ivk
        @test iv > 0 && iv < 10.0
    end
end

@testset "jaeckel" begin
    strikes = [
        0.035123777453185,
        0.049095433048156,
        0.068624781300891,
        0.095922580089594,
        0.134078990076508,
        0.18741338653678,
        0.261963320525776,
        0.366167980681693,
        0.511823524787378,
        0.715418426368358,
        1.0,
        1.39778339939642,
        1.95379843162821,
        2.73098701349666,
        3.81732831143284,
        5.33579814376678,
        7.45829006788743,
        10.4250740447762,
        14.5719954372667,
        20.3684933182917,
        28.4707418310251,
    ]
    vols = [
        0.642412798191439,
        0.621682849924325,
        0.590577891369241,
        0.553137221952525,
        0.511398042127817,
        0.466699250819768,
        0.420225808661573,
        0.373296313420122,
        0.327557513727855,
        0.285106482185545,
        0.249328882881654,
        0.228967051575314,
        0.220857187809035,
        0.218762825294675,
        0.218742183617652,
        0.218432406892364,
        0.217198426268117,
        0.21573928902421,
        0.214619929462215,
        0.2141074555437,
        0.21457985392644,
    ]
    forward = 1.0
    tte = 5.07222222222222
    w1 = ones(length(strikes))
    prices, weights = Collocation.weightedPrices(true, strikes, vols, w1, forward, 1.0, tte, vegaFloor = 1e-5)
    isoc, m = Collocation.makeIsotonicCollocation(strikes, prices, weights, tte, forward, 1.0, deg = 7, degGuess = 1)
    sol = Collocation.Polynomial(isoc)
    ivstrikes = @. Black.impliedVolatility(
        true,
        Collocation.priceEuropean(sol, true, strikes, forward, 1.0),
        forward,
        strikes,
        tte,
        1.0,
    )
    rmse = StatsBase.rmsd(ivstrikes, vols)
    println("poly ", rmse)
    bspl, m = Collocation.makeExpBSplineCollocation(
        strikes,
        prices,
        weights,
        tte,
        forward,
        1.0,
        penalty = 0e-2,
        size = 0,
        minSlope = 1e-8,
        rawFit = true,
    )
    ivstrikes = @. Black.impliedVolatility(
        true,
        Collocation.priceEuropean(bspl, true, strikes, forward, 1.0),
        forward,
        strikes,
        tte,
        1.0,
    )
    rmse = StatsBase.rmsd(ivstrikes, vols)
    println("bspl ", rmse)
    allStrikes = vcat(0.0, strikes, 50.0)
    allPrices = vcat(forward ,prices, 0.0)
    leftB = Math.FirstDerivativeBoundary(-1.0)
    rightB = Math.FirstDerivativeBoundary(0.0)
    cs = Math.makeConvexSchabackRationalSpline(allStrikes, allPrices, leftB, rightB, iterations=128)
    ivstrikes = @. Black.impliedVolatility(
        true,
        cs(strikes),
        forward,
        strikes,
        tte,
        1.0,
    )
    rmse = StatsBase.rmsd(ivstrikes, vols)
    println("scha ", rmse)
end


text:
abstract type Component end

mutable struct Bus
    contents::UInt8
end
value(c::Bus) = c.contents
set_value(c::Bus, v::UInt8) = c.contents = v

mutable struct Register <: Component
    value::UInt8
    bus::Bus
end
value(c::Register) = c.value
set_value(c::Register, v::UInt8) = c.value = v

mutable struct ProgramCounter <: Component
    value::UInt8
    bus::Bus
end
value(c::ProgramCounter) = c.value
set_value(c::ProgramCounter, v::UInt8) = c.value = v
enable(c::ProgramCounter) = c.value += 1

mutable struct ALU <: Component
    in1::Register
    in2::Register
    bus::Bus
end
value(c::ALU) = value(c.in1) + value(c.in2)

mutable struct RAM <: Component
    contents::Vector{UInt8}
    mar::Register
    bus::Bus

    RAM(mar::Register, bus::Bus) = new(zeros(UInt8, 16), mar, bus)
end
value(c::RAM) = c.contents[value(c.mar) << 4 >> 4 + 1]
set_value(c::RAM, v::UInt8) = c.contents[value(c.mar) << 4 >> 4 + 1] = v

output(c::Component) = c.bus.contents = value(c)
input(c::Component) = set_value(c, c.bus.contents)


text:
include("test_common.jl")
using Dates

macro check_type(js_script, jl_type, jl_script=nothing, js_type=nothing)
    quote
        @test $js_script isa $jl_type
        if !isnothing($js_type)
            @test instanceof($jl_script, $js_type)
        end
    end
end

@testset "types" begin
    # Primitive values:
    # `undefined` and `null`: both nothing
    @test node"undefined" ≡ nothing
    @test node"null" ≡ nothing
    # `boolean`:
    @test node"false" || node"true"
    @test node"x => typeof x === 'boolean'"(false)
    # `number`: always `Float64`
    @test node"0" ≡ 0.0 && node"Infinity" ≡ Inf && node"NaN" ≡ NaN
    @test node"x => typeof x === 'number'"(0)
    # `string`
    @test node"'a string'" ≡ "a string"
    @test node"x => typeof x === 'string'".(["str", 'c']) |> all
    # `bigint`
    @test node"282278710948156123453635394322245338923n" == 282278710948156123453635394322245338923
    @test node"x => typeof x === 'bigint'"(BigInt(0))
    # `symbol`
    @test node"Symbol.for('sym')" ≡ :sym
    @test node"x => typeof x === 'symbol'"(:sym)

    # Objects:
    # `object`
    @test node"[]"o isa NodeObject
    @check_type node"{}" JsObject (x=1,) node"Object"
    # `function`
    @check_type node"eval" JsFunction{Nothing} println node"Function"

    # `Date`:
    @check_type node"new Date()" DateTime now() node"Date"
    # `Array`: try to convert
    @check_type node"[]"        Vector{Any} [] node"Array"
    @check_type node"[1]"       Vector{Float64}
    @check_type node"[2, 't']"  Vector{Any}
    @check_type node"[3, 123n]" Vector{Real}
    # Typed arrays:
    @check_type node"new Int8Array(1)"         Vector{Int8}     Int8[]    node"Int8Array"
    @check_type node"new Uint8Array(1)"        Vector{UInt8}    UInt8[]   node"Uint8Array"
    @check_type node"new Uint8ClampedArray(1)" Vector{UInt8}
    @check_type node"new Int16Array(1)"        Vector{Int16}    Int16[]   node"Int16Array"
    @check_type node"new Uint16Array(1)"       Vector{UInt16}   UInt16[]  node"Uint16Array"
    @check_type node"new Int32Array(1)"        Vector{Int32}    Int32[]   node"Int32Array"
    @check_type node"new Uint32Array(1)"       Vector{UInt32}   UInt32[]  node"Uint32Array"
    @check_type node"new Float32Array(1)"      Vector{Float32}  Float32[] node"Float32Array"
    @check_type node"new Float64Array(1)"      Vector{Float64}  Float64[] node"Float64Array"
    @check_type node"new BigInt64Array(1)"     Vector{Int64}    Int64[]   node"BigInt64Array"
    @check_type node"new BigUint64Array(1)"    Vector{UInt64}   UInt64[]  node"BigUint64Array"
    # `Map`, `WeekMap`, `Set`, `WeakSet`
    @check_type node"new Map()"     Dict{Any, Any}
    @check_type node"new WeakMap()" JsObject
    @check_type node"new Set()"     Set{Any}
    @check_type node"new WeakSet()" JsObject
    # `ArrayBuffer`, `SharedArrayBuffer`, `DataView`
    @check_type node"new ArrayBuffer(1)" Vector{UInt8}
    @check_type node"new SharedArrayBuffer(1)" JsObject # Vector{UInt8}
    @check_type node"new DataView(new ArrayBuffer(1))" Vector{UInt8}
    # `Promise`
    @check_type node"new Promise(()=>{})" JsPromise

    o = node"{}"o
    @test value(Rational{Int}, node"0.75"o) == 3//4
    @test value(String, o) == "[object Object]"
    @test value(Int, node"123"o) == 123
    @test value(Int, node"123n"r, is_bigint=true) == 123
    @test value(UInt, node"123n"o) == UInt(123)
end


text:
using GraphIO, LightGraphs, JSON, Plots, GraphRecipes, ParserCombinator
using GraphIO: GML
gr(fmt=:svg)
theme(:juno)

function generate_syn(n₁, n₂, N, k)
    g = barbell_graph(n₁, n₂)
    barabasi_albert!(g, n₁ + n₂ + N, k)
    barabasi_albert!(g, n₁ + n₂ + 2*N, k)
    labels = Dict(i => i <= n₁ + n₂ ? 0 : 1 for i ∈ 1:(n₁ + n₂ + 2*N))
    return g, labels
end

function plot_syn_graph(g, labels)
    graphplot(g, nodecolor=map((x) -> labels[x] == 0 ? :green : :blue, 1:length(labels)))
end

function write_to_file(g, labels, n)
    open("Desktop/syn_graph_labels_$(n).json", "w") do io
        write(io, JSON.json(labels))
    end;

    open("Desktop/syn_graph_$(n).gml", "w") do io
        GML.savegml(io, g)
    end;
end

g, labels = generate_syn(8, 8, 4, 3)
write_to_file(g, labels, 3)


text:
@testset "Games" begin

@testset "registered game attributes" begin
    games = Dict(short_name(game_info) => game_info for game_info in registered_games())
    @test dynamics(games["kuhn_poker"]) == SEQUENTIAL
    @test chance_mode(games["kuhn_poker"]) == EXPLICIT_STOCHASTIC
    @test information(games["kuhn_poker"]) == IMPERFECT_INFORMATION
    @test utility(games["kuhn_poker"]) == ZERO_SUM
    @test min_num_players(games["kuhn_poker"]) == 2
end

@testset "create game" begin
    game = load_game("kuhn_poker")
    game_info = get_type(game)
    @test information(game_info) == IMPERFECT_INFORMATION
    @test num_players(game) == 2
end

@testset "play kuhn_poker" begin
    game = load_game("kuhn_poker")
    @test game == load_game("kuhn_poker")
    @test hash(game) == hash(load_game("kuhn_poker"))

    state = new_initial_state(game)
    @test game == load_game("kuhn_poker")
    @test hash(game) == hash(load_game("kuhn_poker"))

    @test is_chance_node(state) == true
    @test chance_outcomes(state) == [0 => 1/3, 1 => 1/3, 2 => 1/3]

    apply_action(state, 1)
    @test is_chance_node(state) == true
    @test chance_outcomes(state) == [0 => 1/2, 2 => 1/2]

    apply_action(state, 2)
    @test is_chance_node(state) == false
    @test legal_actions(state) == [0, 1]
end

@testset "tic_tac_toe" begin
    game = load_game("tic_tac_toe")
    state = new_initial_state(game)
    @test is_chance_node(state) == false
    @test is_terminal(state) == false
    @test legal_actions(state) == 0:8
end

@testset "GameParameter" begin
    io = IOBuffer()
    print(io, GameParameter(true))
    @test String(take!(io)) == "GameParameter(bool_value=True)"
    print(io, GameParameter(false))
    @test String(take!(io)) == "GameParameter(bool_value=False)"
    print(io, GameParameter("one"))
    @test String(take!(io)) == "GameParameter(string_value='one')"
    print(io, GameParameter(1))
    @test String(take!(io)) == "GameParameter(int_value=1)"
    print(io, GameParameter(1.0))
    @test String(take!(io)) == "GameParameter(double_value=1)"
    print(io, GameParameter(1.2))
    @test String(take!(io)) == "GameParameter(double_value=1.2)"
end

@testset "simultaneous game history" begin
    game = load_game("coop_box_pushing")
    state = new_initial_state(game)
    apply_action(state, 0)
    state2 = new_initial_state(game)
    apply_action(state2, fill(0, num_players(game)))
    @test history(state) == history(state2)
end

@testset "Matrixrame" begin
    matrix_blotto = load_matrix_game("blotto")
    @test num_rows(matrix_blotto) == 66
    @test num_cols(matrix_blotto) == 66

    kuhn_game = load_game("kuhn_poker")
    kuhn_matrix_game = extensive_to_matrix_game(kuhn_game)
    @test num_rows(kuhn_matrix_game) == 64
    @test num_cols(kuhn_matrix_game) == 64
end

end


text:
mutable struct CustomCallback{E<:FitEvent,P<:AbstractFittingPhase} <: AbstractCallback
    fn::Any
    everyn::Int
    _current::Int
    CustomCallback{E,P}(fn, everyn = 1) where {E,P} = new{E,P}(fn, everyn, 0)
end

function on(
    ::E,
    ::P,
    cb::CustomCallback{E,P},
    learner,
) where {E<:FitEvent,P<:AbstractFittingPhase}
    cb._current += 1
    if cb._current % cb.everyn == 0
        cb.fn(learner)
    end
end


text:
module JuliaCtags

export get_config

function write_config(filename = "./julia_ctags")
    sourcefile = joinpath(@__DIR__, "ctags")
    cp(sourcefile, filename)
end

end # module

max_stars_repo_name: m-wells/julia-ctags

max_issues_repo_name: m-wells/julia-ctags

max_forks_repo_name: m-wells/julia-ctags


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

"""
    BoundedSearcher(searcher, nmax)

A method for searching at most `nmax` neighbors using `searcher`.
"""
struct BoundedSearcher{S<:AbstractNeighborSearcher} <: AbstractBoundedNeighborSearcher
  searcher::S
  nmax::Int
end

object(searcher::BoundedSearcher) = object(searcher.searcher)

maxneighbors(searcher::BoundedSearcher) = searcher.nmax

function search!(neighbors::AbstractVector{Int}, xₒ::AbstractVector,
                 searcher::BoundedSearcher; mask=nothing)
  locs = search(xₒ, searcher.searcher)
  nmax = searcher.nmax

  nneigh = 0
  @inbounds for loc in locs
    if mask[loc]
      nneigh += 1
      neighbors[nneigh] = loc
    end
    nneigh == nmax && break
  end

  nneigh
end


text:
#=
Functions to compute quadratic sums

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-08-19
=#

"""
Computes the expected discounted quadratic sum

    q(x_0) = E sum_{t=0}^{infty} beta^t x_t' H x_t

Here {x_t} is the VAR process x_{t+1} = A x_t + C w_t with {w_t}
standard normal and x_0 the initial condition.

##### Arguments
- `A::Union{Float64, Matrix{Float64}}` The n x n matrix described above (scalar)
if n = 1
- `C::Union{Float64, Matrix{Float64}}` The n x n matrix described above (scalar)
if n = 1
- `H::Union{Float64, Matrix{Float64}}` The n x n matrix described above (scalar)
if n = 1
- `beta::Float64`: Discount factor in (0, 1)
- `x_0::Union{Float64, Vector{Float64}}` The initial condtion. A conformable
array (of length n) or a scalar if n=1

##### Returns

- `q0::Float64` : Represents the value q(x_0)

##### Notes

The formula for computing q(x_0) is q(x_0) = x_0' Q x_0 + v where

- Q is the solution to Q = H + beta A' Q A and
- v = \trace(C' Q C) \beta / (1 - \beta)

"""
function var_quadratic_sum(A::ScalarOrArray, C::ScalarOrArray, H::ScalarOrArray,
                           bet::Real, x0::ScalarOrArray)
    n = size(A, 1)

    # coerce shapes
    A = reshape([A;], n, n)
    C = reshape([C;], n, n)
    H = reshape([H;], n, n)
    x0 = reshape([x0;], n)

    # solve system
    Q = solve_discrete_lyapunov(sqrt(bet) .* A', H)
    cq = C'*Q*C
    v = trace(cq) * bet / (1 - bet)
    q0 = x0'*Q*x0 + v
    return q0[1]
end

"""
Computes the quadratic sum

    V = sum_{j=0}^{infty} A^j B A^{j'}

V is computed by solving the corresponding discrete lyapunov equation using the
doubling algorithm.  See the documentation of `solve_discrete_lyapunov` for
more information.

##### Arguments

- `A::Matrix{Float64}` : An n x n matrix as described above.  We assume in order
for convergence that the eigenvalues of A have moduli bounded by unity
- `B::Matrix{Float64}` : An n x n matrix as described above.  We assume in order
for convergence that the eigenvalues of B have moduli bounded by unity
- `max_it::Int(50)` : Maximum number of iterations

##### Returns

- `gamma1::Matrix{Float64}` : Represents the value V

"""
function m_quadratic_sum(A::Matrix, B::Matrix; max_it=50)
    solve_discrete_lyapunov(A, B, max_it)
end


text:
############################################################################################
"""
$(TYPEDEF)

Contains temporary buffer values to avoid allocations during particle propagation.

# Fields
$(TYPEDFIELDS)
"""
struct ParticleBuffer{A,I,P}
    "Contains buffer values for particles and ancestor for one iteration."
    parameter::BaytesCore.ParameterBuffer{A, I}
    "Proposal trajectory and predicted latent variable."
    proposal::Vector{A}
    "Predicted latent and oberved data"
    prediction::Vector{P}
    "Stores boolean if resampled at each iteration."
    resampled::Vector{Bool}
    function ParticleBuffer(
        prediction::Vector{P}, reference::AbstractVector{A}, Nparticles::Integer, Ndata::Integer, F::Type{I}
    ) where {P,A,I<:Integer}
        parameter = BaytesCore.ParameterBuffer(reference, Nparticles, F)
        proposal = typeof(reference)(undef, Ndata)
        resampled = zeros(Bool, Ndata)
        return new{A,I,P}(parameter, proposal, prediction, resampled)
    end
end
function update!(buffer::ParticleBuffer, Nparticles::Integer, Ndata::Integer)
    resize!(buffer.parameter, Nparticles)
    resize!(buffer.proposal, Ndata)
    resize!(buffer.resampled, Ndata)
    return nothing
end

############################################################################################
#export
export ParticleBuffer, update!


text:
function Base.show(io::IO, x::ParameterGroup)
    s1 = maximum(length.(string.(fieldnames(typeof(x)),"Variable")))
    vals = Dict{Symbol,Any}()
    desc = Dict{Symbol,Any}()
    for name in fieldnames(typeof(x))
        if typeof(getfield(x,name))<:AbstractParameter
            if typeof(getfield(getfield(x,name),:value)) <: Real
                vals[name] = round.(getfield(getfield(x,name),:value),digits = 3)
                desc[name] =  getfield(getfield(x,name),:description)
            else
                vals[name] = getfield(getfield(x,name),:value)
                desc[name] =  getfield(getfield(x,name),:description)
            end
        else
            vals[name] = getfield(x,name)
            desc[name] = "No Description"
        end
    end
    s2 = maximum(length.(string.(values(desc))))
    s3 = length.(string.("Values      "))
    s3 = s3+s2+s1
    print(io,"$(rpad("Variable",s1," ")) | ")
    print(io,"$(rpad("Description",s2," ")) | ")
    print(io,"Value\n")
    print(io,rpad("-",s3,"-"),"\n")
    for name in fieldnames(typeof(x))
        if typeof(getfield(x,name))<:AbstractParameter
            print(io,rpad(name,s1," ")," | ")
            print(io,"$(rpad(getfield(getfield(x,name),:description),s2," ")) | ")
            print(IOContext(io, :compact => true,:displaysize => (1,2),:limit=>true),vals[name],"\n")
        else
            print(io,rpad(name,s1," ")," | ")
            print(io,"$(rpad("No Description",s2," ")) | ")
            print(IOContext(io, :compact => true,:displaysize => (1,2),:limit=>true),vals[name],"\n")
        end
    end
end

function Base.show(io::IO,::MIME"text/plain", x::ParameterGroup)
    s1 = maximum(length.(string.(fieldnames(typeof(x)),"Variable")))
    vals = Dict{Symbol,Any}()
    desc = Dict{Symbol,Any}()
    for name in fieldnames(typeof(x))
        if typeof(getfield(x,name))<:AbstractParameter
            if typeof(getfield(getfield(x,name),:value)) <: Real
                vals[name] = round.(getfield(getfield(x,name),:value),digits = 3)
                desc[name] =  getfield(getfield(x,name),:description)
            else
                vals[name] = getfield(getfield(x,name),:value)
                desc[name] =  getfield(getfield(x,name),:description)
            end
        else
            vals[name] = getfield(x,name)
            desc[name] = "No Description"
        end
    end
    s2 = maximum(length.(string.(values(desc))))
    s3 = length.(string.("Values      "))
    s3 = s3+s2+s1
    print(io,"$(rpad("Variable",s1," ")) | ")
    print(io,"$(rpad("Description",s2," ")) | ")
    print(io,"Value\n")
    print(io,rpad("-",s3,"-"),"\n")
    for name in fieldnames(typeof(x))
        if typeof(getfield(x,name))<:AbstractParameter
            print(io,rpad(name,s1," ")," | ")
            print(io,"$(rpad(getfield(getfield(x,name),:description),s2," ")) | ")
            print(IOContext(io, :compact => true,:displaysize => (1,2),:limit=>true),vals[name],"\n")
        else
            print(io,rpad(name,s1," ")," | ")
            print(io,"$(rpad("No Description",s2," ")) | ")
            print(IOContext(io, :compact => true,:displaysize => (1,2),:limit=>true),vals[name],"\n")
        end
    end
end

function Base.show(io::IO,::MIME"text/plain", x::ParameterSpace)
    for str in fieldnames(typeof(x))
        println(io,"Parameters for ",typeof(getfield(x,str)),":")
        println(io,getfield(x,str))
    end
end

function Base.show(io::IO, x::ParameterSpace)
    for str in fieldnames(typeof(x))
        println(io,"Parameters for ",typeof(getfield(x,str)),":")
        println(io,getfield(x,str))
    end
end


text:
################################################################################
# Dataset related utilities
################################################################################
function get_data_dir()
  data_dir = joinpath(Pkg.dir("MXNet"), "data")
  mkpath(data_dir)
  data_dir
end

function get_mnist_ubyte()
  data_dir  = get_data_dir()
  mnist_dir = joinpath(data_dir, "mnist")
  mkpath(mnist_dir)
  filenames = Dict(:train_data  => "train-images-idx3-ubyte",
                   :train_label => "train-labels-idx1-ubyte",
                   :test_data   => "t10k-images-idx3-ubyte",
                   :test_label  => "t10k-labels-idx1-ubyte")
  filenames = Dict([k => joinpath(mnist_dir, v) for (k,v) in filenames])
  if !all(isfile, values(filenames))
    cd(mnist_dir) do
      mnist_dir = download("http://data.dmlc.ml/mxnet/data/mnist.zip", "mnist.zip")
        try
          run(`unzip -u $mnist_dir`)
        catch
          try
            run(pipe(`7z x $mnist_dir`,stdout=DevNull))
          catch
            error("Extraction Failed:No extraction program found in path")
          end
      end
    end
  end
  return filenames
end

function get_cifar10()
  data_dir    = get_data_dir()
  cifar10_dir = joinpath(data_dir, "cifar10")
  mkpath(cifar10_dir)
  filenames = Dict(:train => "cifar/train.rec", :test => "cifar/test.rec")
  filenames = Dict([k => joinpath(cifar10_dir, v) for (k,v) in filenames])
  if !all(isfile, values(filenames))
    cd(cifar10_dir) do
      download("http://data.dmlc.ml/mxnet/data/cifar10.zip", "cifar10.zip")
        try
          run(`unzip -u cifar10.zip`)
        catch
          try
            run(pipeline(`7z x cifar10.zip`, stdout=DevNull))
          catch
            error("Extraction Failed:No extraction program found in path")
          end
      end
    end
  end

  filenames[:mean] = joinpath(cifar10_dir, "cifar/cifar_mean.bin")
  return filenames
end

################################################################################
# Internal Utilities
################################################################################
function _format_typestring(typestr :: String)
  replace(typestr, r"\bSymbol\b", "SymbolicNode")
end
function _format_docstring(narg::Int, arg_names::Ref{char_pp}, arg_types::Ref{char_pp}, arg_descs::Ref{char_pp}, remove_dup::Bool=true)
  param_keys = Set{String}()

  arg_names  = unsafe_wrap(Array, arg_names[], narg)
  arg_types  = unsafe_wrap(Array, arg_types[], narg)
  arg_descs  = unsafe_wrap(Array, arg_descs[], narg)
  docstrings = String[]

  for i = 1:narg
    arg_name = unsafe_string(arg_names[i])
    if arg_name ∈ param_keys && remove_dup
      continue
    end
    push!(param_keys, arg_name)

    arg_type = _format_typestring(unsafe_string(arg_types[i]))
    arg_desc = unsafe_string(arg_descs[i])
    push!(docstrings, "* `$arg_name::$arg_type`: $arg_desc\n")
  end
  return join(docstrings, "\n")
end

function _format_signature(narg::Int, arg_names::Ref{char_pp})
  arg_names  = unsafe_wrap(Array, arg_names[], narg)

  return join([unsafe_string(name) for name in arg_names] , ", ")
end


text:
get_total_fitness_variance(history::History) = var(get_fitness.(get_detailed_history(history)))

function get_inherent_variance(history::History)
    organisms = get_detailed_history(history)
    return mean(var(mean(organisms, [EnvironmentType(), PopulationType(), NumOffspringsType()]), [EnvironmentType(), PopulationType()]))
end

function get_mutational_variance(history::History)
    organisms = get_detailed_history(history)
    return mean(var(organisms, [EnvironmentType(), PopulationType(), NumOffspringsType()]))
end

function get_abiotic_variance(history::History)
    organisms = get_detailed_history(history)
    return var(mean(organisms, [EnvironmentType()]))
end

function get_biotic_variance(history::History)
    organisms = get_detailed_history(history)
    return mean(var(mean(organisms, [EnvironmentType(), PopulationType()]), [EnvironmentType()]))
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
@testset "DMD" begin
    # Create some linear data
    A = [0.9 -0.2; 0.0 0.2]
    y = [[10.; -10.]]
    for i in 1:10
        push!(y, A*y[end])
    end
    X = hcat(y...)
    @test_throws AssertionError ExactDMD(X[:, 1:end-2], dt = -1.0)
    estimator = ExactDMD(X[:,1:end-2])
    @test isstable(estimator)
    @test estimator.Ã ≈ A
    @test eigvals(estimator) ≈ eigvals(A)
    @test eigvecs(estimator) ≈ eigvecs(A)
    @test_nowarn dynamics(estimator)
    @test_throws AssertionError dynamics(estimator, discrete = false)
    @test_nowarn update!(estimator, X[:, end-1], X[:,end])
end

@testset "EDMD" begin
    # Test for linear system
    function linear_sys(u, p, t)
        x = -0.9*u[1]
        y = -0.3*u[2]
        return [x;y]
    end

    u0 = [π; 1.0]
    tspan = (0.0, 20.0)
    prob = DiscreteProblem(linear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())

    @variables u[1:2]
    h = [1u[1]; 1u[2]; sin(u[1]); cos(u[1]); u[1]*u[2]]
    basis = Basis(h, u)

    @test_throws AssertionError ExtendedDMD(sol[:,:], basis, dt = -1.0)
    estimator = ExtendedDMD(sol[:,:], basis)
    @test basis == estimator.basis
    basis_2 = reduce_basis(estimator, threshold = 1e-5)
    @test size(basis_2)[1] < size(basis)[1]
    estimator_2 = ExtendedDMD(sol[:,:], basis_2)
    p1 = DiscreteProblem(dynamics(estimator), u0, tspan, [])
    s1 = solve(p1,FunctionMap())
    p2 = DiscreteProblem(dynamics(estimator_2), u0, tspan, [])
    s2 = solve(p2,FunctionMap())
    p3 = DiscreteProblem(linear_dynamics(estimator_2), estimator_2(u0), tspan, [])
    s3 = solve(p3,FunctionMap())
    @test sol[:,:] ≈ s1[:,:]
    @test sol[:,:] ≈ s2[:,:]
    @test estimator_2.basis(sol[:,:])≈ s3[:,:]
    @test eigvals(estimator_2) ≈ [-0.9; -0.3]

    # Test for nonlinear system
    function nonlinear_sys(du, u, p, t)
        du[1] = sin(u[1])
        du[2] = -0.3*u[2] -0.9*u[1]
    end

    prob = DiscreteProblem(nonlinear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())
    estimator = ExtendedDMD(sol[:,:], basis)
    p4 = DiscreteProblem(dynamics(estimator), u0, tspan, [])
    s4 = solve(p4,FunctionMap())
    @test sol[:,:] ≈ s4[:,:]
end

@testset "DMDc" begin
    # Define measurements from unstable system with known control input
    X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007]
    U = [-4 -2 -1 -0.5]
    B = Float32[1; 0]

    # But with a little more knowledge
    sys = DMDc(X, U, B = B)
    @test isa(get_dynamics(sys), ExactDMD)
    @test sys.koopman.Ã ≈[1.5 0; 0 0.1]
    @test get_input_map(sys) ≈ [1.0; 0.0]
    @test !isstable(sys)
    @test_nowarn eigen(sys)

    # Check the solution of an unforced and forced system against each other
    dudt_ = dynamics(sys)
    prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
    sol_unforced = solve(prob,  FunctionMap())

    dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1])
    prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
    sol = solve(prob, FunctionMap())

    @test all(abs.(diff(sol[1,:])) .< 1e-5)
    @test sol[2,:] ≈ sol_unforced[2,:]
end


text:
using HTTP
using JSON
using DataFrames

"""
For obtaining the address in the format to use with the API key `address()` can be used.

Example:

julia> address("1600 Amphitheatre Parkway, Mountain View, CA")

"""

function address(street_address::Vector{String})
    street_address = replace.(street_address, r"\s*(,)?\s+" => "+")
    return(street_address)
end

"""
        `julia> geocode(street_address::String, api_key::String, method::String)`

JuliaMaps has `geocode()` for obtaining the coordingates of locations. It will produce `longitude` and `latitude` of location in `Dict`.

    Usage:
    `julia> geocode("full address", "api key", "method")`

    Note: Currently, `geocode()` can call from two methods: "google" and "osm".

If there are more than one addresses, it'd be reueired to use broadcasting operator with `geocode()`.

    Example:
    `julia> geocode.(address_list, "api","method")`

Method:
- google: Calls data from google maps. In this case, it'd be reqired to use google maps api key.
- osm: Calls data from open street maps, using MapQuest api. Visit: https://developer.mapquest.com/documentation/
"""

function geocode(street_address::String, api_key::String, method::String)
    street_address = replace.(street_address, r"\s*(,)?\s+" => "+")
    if method == "google"
        url = "https://maps.googleapis.com/maps/api/geocode/json?address=$street_address&key=$api_key"
        response = HTTP.get(url)
        gem = geomgoogle(JSON.parse(String(response.body)))
    elseif method == "osm"
        url = "http://www.mapquestapi.com/geocoding/v1/address?key=$api_key&location=$street_address"
        response = HTTP.get(url)
        gem = geomOSM(JSON.parse(String(response.body)))
    else
        gem = error("Please use method correctly. Currently JuliaMaps.jl only supports two methods: google and osm")
    end
    return(gem)
end

"""
If there were more than one addresses in geocoding, while using `geocode()`, all the longitudes and latitudes can be printed in a data frame using `geometry(geocod_output)`.

    Usage:
    `julia> geometry(geocode_output)`

Note: This function can only be used if there were more than one addresses in `geocode()`.
"""

function geometry(geocode_output)
    lat = getindex.(geocode_output, "lat")
    lng = getindex.(geocode_output, "lng")
    d   = DataFrame(Longitude = lng, Latitude = lat)
    return(d)
end


text:
using ReactionCoordinates
using Test

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


text:
module ExampleRajith
using ForwardDiff

include("extra_file.jl")

export my_f, derivative_of_my_f

end


text:
module TestMacro

include("preamble.jl")

@testset begin
    @test (@! push!(empty!([0, 1]), 2, 3)) == [2, 3]
    @test (@! push!(empty!((0, 1)), 2, 3)) == (2, 3)
end

@testset ".=" begin
    @test !@isdefined y
    @test begin
        x = SVector(1, 2)
        @! y .= x .* 2
    end === SVector(2, 4)

    @test !@isdefined y
    @test begin
        x = SVector(1, 2)
        y = SVector(0, 0)
        @! y .= x .* 2
    end === SVector(2, 4)

    let x = [1, 2]
        y = [0, 0]
        @test (@! y .= x .* 2)::Vector{Int} == [2, 4]
    end

    let x = [1, 2]
        y = SVector(0, 0)
        @test (@! y .= x .* 2)::Vector{Int} == [2, 4]
    end
end

@testset "@." begin
    @test !@isdefined y
    @test begin
        x = SVector(1, 2)
        @! @. y = x * 2
    end === SVector(2, 4)

    @test !@isdefined y
    @test begin
        x = SVector(1, 2)
        y = SVector(0, 0)
        @! @. y = x * 2
    end === SVector(2, 4)

    let x = [1, 2]
        y = [0, 0]
        @test (@! @. y = x * 2)::Vector{Int} == [2, 4]
    end

    let x = [1, 2]
        y = SVector(0, 0)
        @test (@! @. y = x * 2)::Vector{Int} == [2, 4]
    end
end

end  # module


text:
using DataStructures

#. HASHTABLES

#.. EasyTable
# Building upon Julia's Dict

mutable struct EasyTable{T, I <: Unsigned}
    val::Vector{T}
    rev::Dict{T, I}
end

function easytable(val, ind_type = UInt32)
    # tbl_val = collect(val)
    tbl_val = val
    tbl_rev = Dict{eltype(val), ind_type}(broadcast(x -> reverse(x), enumerate(tbl_val)))
    EasyTable{eltype(val), ind_type}(tbl_val, tbl_rev)
end

ind_type(table::EasyTable{T, I}) where {T, I} = I
Base.getindex(table::EasyTable{T}, i) where T = table.val[i]
Base.length(table::EasyTable) = length(table.val)

@inline find(table::EasyTable{T, I}, v::T) where {T, I} = get(table.rev, v, zero(I))
@inline function findorpush!(table::EasyTable{T, I}, v::T) where {T, I}
    n = find(table, v)
    if n > 0
        return n
    else
        push!(table.val, v)
        n = I(length(table.val))
        table.rev[v] = n
        return n
    end
end

function indexpolynomial(tbl::EasyTable{M, I}, p::Polynomial{M, T}) where {M, I, T}
    Polynomial{I, T}([findorpush!(tbl, m) for m in p.mo], p.co)
end

function unindexpolynomial(tbl::EasyTable{M, I}, p::Polynomial{I, T}) where {M, I, T}
    Polynomial{M, T}([tbl[i] for i in p.mo], p.co)
end

#.. MonomialHashTable
# Mostly for fun, may not be a good candidate for GB computations

mutable struct MonomialHashTable{N, E, I <: Unsigned, B <: Unsigned}
    val::Vector{Monomial{N, E}}
    bitmasks::Vector{B}
    idx::Vector{I}

    max_powers::MVector{N, E}
    min_powers::MVector{N, E}
    bitmask_powers::Dict{Int, Vector{E}}
    mask::I
    size::Int
    maxloop::Int
    hits::Int
    totaldivs::Int

    function MonomialHashTable{N, E, I, B}(; deg_bound = 0) where {N, E, I, B}
        tsize = 16
        @assert 2*tsize <= typemax(I)
        nbits = ndigits(typemax(B), base = 2)
        masks_per_var = even_partition(nbits, N)
        min_powers = SVector{N, E}(zeros(E, N))
        max_powers = SVector{N, E}(repeat([E(deg_bound)], N))
        bitmask_powers = Dict([(i, zeros(E, masks_per_var[i])) for i in 1:N])
        [bitmask_powers[i] = even_between(min_powers[i], max_powers[i], length(bitmask_powers[i]))
         for i in 1:N]
        new(Monomial{N, E}[], B[], zeros(Int, 2*tsize),
            max_powers, min_powers, bitmask_powers, tsize-1, 0, 4, 0, 0)
    end
end

Base.@propagate_inbounds Base.getindex(table::MonomialHashTable, i) = table.val[i]
Base.length(table::MonomialHashTable) = length(table.val)

index1(table::MonomialHashTable{N, E, I}, v::Monomial{N, E}) where {N, E, I} = 2*( (hash(v) % I)&table.mask ) + 1
index2(table::MonomialHashTable{N, E, I}, v::Monomial{N, E}) where {N, E, I} = 2*( ((hash(v) >> 32) % I)&table.mask ) + 2

function find(table::MonomialHashTable{N, E, I}, v::Monomial{N, E}) where {N, E, I}
    @inbounds begin
        n = table.idx[index1(table, v)]
        if n > 0 && table[n] == v
            return n
        end

        n = table.idx[index2(table, v)]
        if n > 0 && table[n] == v
            return n
        end
        return zero(I)
    end
end

function rehash!(table::MonomialHashTable{N, E, I}) where {N, E, I}
    #@warn "REHASH $(length(table)) elems, $(length(table.idx)) index table"

    table.idx = zeros(I, 2*length(table.idx))
    table.maxloop += 2
    table.mask = 2*table.mask + 1

    for n in 1:length(table.val)
        register!(table, I(n))
    end
end

function register!(table::MonomialHashTable{N, E, I}, n::I) where {N, E, I}
    #@warn "register $n"
    @inbounds for _ in 1:table.maxloop
        i1 = index1(table, table[n])
        #@warn "push $n in i1=$i1"
        n, table.idx[i1] = table.idx[i1], n
        if iszero(n)
            return
        end

        i2 = index2(table, table[n])
        #@warn "push $n in i2=$i2"
        n, table.idx[i2] = table.idx[i2], n
        if iszero(n)
            return
        end
    end

    rehash!(table)
end

function findorpush!(table::MonomialHashTable{N, E, I, B}, v::Monomial{N, E}) where {N, E, I, B}
    @inbounds begin
        n = table.idx[index1(table, v)]
        if n > 0 && table[n] == v
            return n
        end

        n = table.idx[index2(table, v)]
        if n > 0 && table[n] == v
            return n
        end
    end

    #@warn "Insert $v"
    push!(table.val, v)
    @inbounds for i in eachindex(v.exponents)
        e = v.exponents[i]
        if e > table.max_powers[i]
            table.max_powers[i] = e
        end
    end
    push!(table.bitmasks, bitmask(v, table.bitmask_powers, masktype = B))
    table.size += 1
    n = I(table.size)

    register!(table, n)
    return n
end

function remask!(table::MonomialHashTable{N, E, I, B}) where {N, E, I, B}

    [table.bitmask_powers[i] = even_between(table.min_powers[i], table.max_powers[i], length(table.bitmask_powers[i]))
     for i in 1:N]
    table.bitmasks = broadcast(v -> bitmask(v, table.bitmask_powers), table.val)
end

#. INDEXED MONOMIALS

mutable struct IxMonomialΓ{I<:Unsigned, N, E, B, MΓ<:MonomialContext{Monomial{N, E}}}<:MonomialContext{I}
    ctx::MΓ
    table::MonomialHashTable{N, E, I, B}
end

# need to change this?
const IxPolynomialΓ{I, T, moΓ <: IxMonomialΓ{I}, coΓ} = PolynomialΓ{I, T, moΓ, coΓ}

function ixmonomialctx(moctx=nothing; indices=UInt32, mask_type=UInt32, remask_after=1, deg_bound = 0, kwargs...)
    if isnothing(moctx)
        moctx = monomialctx(;kwargs...)
    end
    idxmoctx = IxMonomialΓ{indices, params(moctx)..., mask_type, typeof(moctx)}(moctx, MonomialHashTable{params(moctx)..., indices, mask_type}(deg_bound = deg_bound))
    idxmoctx
end

(ctx::IxMonomialΓ{I, N, E, B})(x::Monomial{N, E}) where {I, N, E, B} = findorpush!(ctx.table, ctx.ctx(x))
(ctx::IxMonomialΓ{I})(i::I) where I = i
(ctx::IxMonomialΓ)(x::AA.MPolyElem) = findorpush!(ctx.table, ctx.ctx(x))

Base.@propagate_inbounds Base.getindex(ctx::IxMonomialΓ{I}, i::I) where I = ctx.table[i]

function lt(ix::IxMonomialΓ{I}, a::I, b::I) where I
    lt(ix.ctx, ix[a], ix[b])
end

#.. Monomial interface

pretty_print(ctx::IxMonomialΓ{I}, m::I) where I = "$(ctx[m])"

valtype(ctx::IxMonomialΓ) = eltype(ctx.ctx)
variables(ctx::IxMonomialΓ) = variables(ctx.ctx)
nvars(ctx::IxMonomialΓ) = nvars(ctx.ctx)
exponenttype(ctx::IxMonomialΓ) = exponenttype(ctx.ctx)

termorder(ctx::IxMonomialΓ) = termorder(ctx.ctx)
exponents(ctx::IxMonomialΓ{I}, i::I) where I = exponents(ctx.ctx, ctx[i])
degree(ctx::IxMonomialΓ{I}, i::I) where I = sum(exponents(ctx, i))

Base.one(ctx::IxMonomialΓ) = ctx(one(valtype(ctx)))

mul(ctx::IxMonomialΓ{I}, i::I, j::I) where I = ctx(mul(ctx.ctx, ctx[i], ctx[j]))
div(ctx::IxMonomialΓ{I}, i::I, j::I) where I = ctx(div(ctx.ctx, ctx[i], ctx[j]))
lcm(ctx::IxMonomialΓ{I}, i::I, j::I) where I = ctx(lcm(ctx.ctx, ctx[i], ctx[j]))
gcd(ctx::IxMonomialΓ{I}, i::I, j::I) where I = ctx(gcd(ctx.ctx, ctx[i], ctx[j]))

function divides(ctx::IxMonomialΓ{I, N},
                 i::I,
                 j::I) where {I, N}

    # iszero(ctx.table.bitmasks[i]  & (~ ctx.table.bitmasks[j])) && divides(ctx.ctx, ctx[i], ctx[j])
    ctx.table.totaldivs += 1
    if iszero(ctx.table.bitmasks[i]  & (~ ctx.table.bitmasks[j]))
        return divides(ctx.ctx, ctx[i], ctx[j])
    end
    ctx.table.hits += 1
    false
end


text:
using DiffPrecTest
using JLD

resPath = ARGS[1]
outPrefix = ARGS[2]

# resPath = /scratch/mkolar/diffTest/sim1/res_2_2_$(rep).jld

NUM_REP = 1000
results = Array{Any}(undef, 9, NUM_REP)
for rep=1:NUM_REP
    global results
    fname = "$(resPath)_$(rep).jld"

    # try
      file = jldopen(fname, "r")
      res = read(file, "est")
      close(file)

      results[:, rep] = deepcopy(res)
    # catch
    #     @show fname
    # end
end

@save "$(outPrefix).jld" results


text:
using Lowranksvd
using Documenter

DocMeta.setdocmeta!(Lowranksvd, :DocTestSetup, :(using Lowranksvd); recursive=true)

makedocs(;
    modules=[Lowranksvd],
    authors="jieli-matrix <li_j20@fudan.edu.cn> and contributors",
    repo="https://github.com/jieli-matrix/Lowranksvd.jl/blob/{commit}{path}#{line}",
    sitename="Lowranksvd.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jieli-matrix.github.io/Lowranksvd.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jieli-matrix/Lowranksvd.jl",
)


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

using EmpiricalDistributions
using Test

using Random, LinearAlgebra
using Distributions, StatsBase, ArraysOfArrays
using Adapt

@testset "mv_binned_dist" begin
    μ = [1.23, -0.67]
    Σ = [0.45 0.32; 0.32 0.76]' * [0.45 0.32; 0.32 0.76]
    true_dist = MvNormal(μ, Σ)
    h = Histogram((μ[1]-10Σ[1]:Σ[1]/10:μ[1]+10Σ[1], μ[2]-10Σ[4]:Σ[4]/10:μ[2]+10Σ[4]))
    n = 10^7
    r = rand(true_dist, n)
    append!(h, (r[1, :], r[2, :]))

    @test @inferred(MvBinnedDist(h)) isa MvBinnedDist{<:Real,2}
    d = MvBinnedDist(h)
    @test @inferred(MvBinnedDist{Float32}(h)) isa MvBinnedDist{Float32,2}
    d_32 = MvBinnedDist{Float32}(h)
    @test @inferred(MvBinnedDist{Float64}(h)) isa MvBinnedDist{Float64,2}
    d_64 = MvBinnedDist{Float64}(h)

    @inferred(adapt(EmpiricalDistributions._AdaptToArray(), d))._edges[1] isa Array

    @test @inferred(convert(Histogram, d)) == normalize(h)

    @test @inferred(length(d)) == 2
    @test @inferred(size(d)) == (2,)
    @test @inferred(eltype(d)) == Float64

    @test all(isapprox.(mean(true_dist), @inferred(mean(d)), atol = 0.01))
    @test all(isapprox.(mode(true_dist), @inferred(mode(d)), atol = 0.2))
    @test all(isapprox.(var(true_dist), @inferred(var(d)), atol = 0.01))
    @test all(isapprox.(cov(true_dist), @inferred(cov(d)), atol = 0.01))

    r = zeros(2, 10^6)
    @test @inferred(rand!(d, r)) === r

    fit_dist = fit(MvNormal, r)
    @test all(isapprox.(μ, fit_dist.μ,     atol = 0.01))
    @test all(isapprox.(Σ, fit_dist.Σ.mat, atol = 0.01))

    @test @inferred(pdf(d, r[:,1])) isa Real
    @test isapprox(pdf.(Ref(true_dist), nestedview(r)), pdf.(Ref(d), nestedview(r)), rtol = 0.1)
    @test @inferred(pdf(d, [1000, 1000])) == 0
    @test @inferred(pdf(d, [-1000, -1000])) == 0

    @test @inferred(logpdf(d, r[:,1])) == log(pdf(d, r[:,1]))
    @test @inferred(logpdf(d, [1000, 1000])) == -Inf
    @test @inferred(logpdf(d, [-1000, -1000])) == -Inf

    @test @inferred(rand(d)) isa Vector{Float64}
    @test size(rand(d)) == (length(d),)
    @test @inferred(rand(d_32)) isa Vector{Float32}
    @test @inferred(rand(d_64)) isa Vector{Float64}

    @test @inferred(rand(d, 10)) isa Matrix{Float64}
    @test size(rand(d, 10)) == (length(d), 10)
end


text:
module TestElastic
using Test, GaussianProcesses

@testset "elastic" begin
    N = 2; nobs = 20;
    for m in [MeanLin(rand(N)), MeanZero(), MeanConst(rand())]
        for k in [Mat52Iso(rand(), rand()), SEIso(rand(), rand()),
                  Mat32Ard(rand(N), rand()), RQArd(rand(N), rand(), rand()),
                  LinIso(rand()), LinArd(rand(N))]
            x0 = Array{Float64, 2}(undef, N, 0); y0 = Array{Float64, 1}(undef, 0);
            x = rand(N, nobs); y = rand(nobs);
            gp1 = GPE(x, y, m, k, -2.)
            gp2 = ElasticGPE(x0, y0, m, k, -2.)
            append!(gp2, x[:, 1:7], y[1:7])
            append!(gp2, x[:, 8], y[8])
            append!(gp2, x[:, 9:20], y[9:20])
            @test gp1.cK.chol.U ≈ view(gp2.cK.chol).U
            xtest = rand(N, 2)
            mu1, sig1 = predict_y(gp1, xtest)
            mu2, sig2 = predict_y(gp2, xtest)
            @test mu1 ≈ mu2
            @test sig1 ≈ sig2
            @test gp1.mll ≈ gp2.mll
            @test isapprox(gp1.alpha, gp2.alpha, atol = 1e-4)
            optimize!(gp1)
            optimize!(gp2)
            @test gp1.mll ≈ gp2.mll
            @test isapprox(gp1.alpha, gp2.alpha, atol = 1e-4)
        end
    end
    gp = ElasticGPE(N, kernel = SEIso(0., 0.), capacity = 50)
    @test gp.cK.chol.capacity == 50
    @test gp.cK.mat.m.capacity == gp.data.R.capacity == (50, 50)
end

end


text:
using Pkg
Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
using MLCourse, CSV, MLJ, DataFrames, MLJLinearModels, Random, Plots
include("./Data_Processing.jl")

#Load the data
drop = generate(option = "drop", std = "false", valid = "false", test = "true");
med = generate(option = "med", std = "false", valid = "false", test = "true");

#Test if the machine has a greater AUC on the drop data set or the med data set
mach_drop = machine(LogisticClassifier(), select(drop.train[:,:], Not(:precipitation_nextday)), drop.train.precipitation_nextday) |> fit!;
mach_med = machine(LogisticClassifier(), select(med.train[:,:], Not(:precipitation_nextday)), med.train.precipitation_nextday) |> fit!;
losses(mach_drop, select(drop.test[:,:], Not(:precipitation_nextday)), drop.test.precipitation_nextday) # 0.79
losses(mach_med, select(med.test[:,:], Not(:precipitation_nextday)), med.test.precipitation_nextday) # 0.88
# AUC_drop = 0.786, AUC_med = 0.880

#Function to tune the model with L1 and L2 regularization with log scaled lambda values
function TunedModel_L2(data)
	Random.seed!(2711)
	model = LogisticClassifier(penalty = :l2)
	self_tuning_model = TunedModel(model = model,
	                         resampling = CV(nfolds = 10),
	                         tuning = Grid(goal = 100),
	                         range = range(model, :lambda,
									    lower = 1e1, upper = 1e3, scale = :log),
	                         measure = auc)
    self_tuning_mach = machine(self_tuning_model,
							select(data[:,:], Not(:precipitation_nextday)),
							data.precipitation_nextday) |> fit!;
	self_tuning_mach
end

function TunedModel_L1(data)
	Random.seed!(2711)
	model = LogisticClassifier(penalty = :l1)
	self_tuning_model = TunedModel(model = model,
	                         resampling = CV(nfolds = 10),
	                         tuning = Grid(goal = 100),
	                         range = range(model, :lambda,
									    lower = 1e-1, upper = 1e2, scale = :log),
	                         measure = auc)
    self_tuning_mach = machine(self_tuning_model,
							select(data[:,:], Not(:precipitation_nextday)),
							data.precipitation_nextday) |> fit!;
	self_tuning_mach
end

#Load the standardized data
drop_std = generate(option = "drop", std = "true", valid = "false", test = "true");
med_std = generate(option = "med", std = "true", valid = "false", test = "true");

#Test de best model on the drop and med data set with L1 or L2 regularization
rep_drop_L2 = report(TunedModel_L2(drop_std.train));
scatter(reshape(rep_drop_L2.plotting.parameter_values, :), rep_drop_L2.plotting.measurements, xlabel = "Lamba", ylabel = "AUC")
rep_drop_L1 = report(TunedModel_L1(drop_std.train));
scatter(reshape(rep_drop_L1.plotting.parameter_values, :), rep_drop_L1.plotting.measurements, xlabel = "Lamba", ylabel = "AUC")
rep_med_L2 = report(TunedModel_L2(med_std.train));
scatter(reshape(rep_med_L2.plotting.parameter_values, :), rep_med_L2.plotting.measurements, xlabel = "Lamba", ylabel = "AUC")
rep_med_L1 = report(TunedModel_L1(med_std.train));
scatter(reshape(rep_med_L1.plotting.parameter_values, :), rep_med_L1.plotting.measurements, xlabel = "Lamba", ylabel = "AUC")

#Test the best models over the test set with the best hyper-parameters to pick the best one
best_mach_drop_L2 = machine(LogisticClassifier(penalty = :l2, lambda = rep_drop_L2.best_model.lambda), select(drop_std.train[:,:], Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)|> fit!
best_mach_drop_L1 = machine(LogisticClassifier(penalty = :l1, lambda = rep_drop_L1.best_model.lambda), select(drop_std.train[:,:], Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)|> fit!
best_mach_med_L2 = machine(LogisticClassifier(penalty = :l2, lambda = rep_med_L2.best_model.lambda), select(med_std.train[:,:], Not(:precipitation_nextday)), med_std.train.precipitation_nextday)|> fit!
best_mach_med_L1 = machine(LogisticClassifier(penalty = :l1, lambda = rep_med_L1.best_model.lambda), select(med_std.train[:,:], Not(:precipitation_nextday)), med_std.train.precipitation_nextday)|> fit!
losses(best_mach_drop_L2, select(drop_std.test, Not(:precipitation_nextday)), drop_std.test.precipitation_nextday)
losses(best_mach_drop_L1, select(drop_std.test, Not(:precipitation_nextday)), drop_std.test.precipitation_nextday)
losses(best_mach_med_L2, select(med_std.test, Not(:precipitation_nextday)), med_std.test.precipitation_nextday)
losses(best_mach_med_L1, select(med_std.test, Not(:precipitation_nextday)), med_std.test.precipitation_nextday)
# AUC_DROP_L2 = 0.917, AUC_DROP_L1 = 0.914, AUC_MED_L2 = 0.914, AUC_MED_L1 = 0.912

#Write in the submission file with the best machine trained on all the data
train, test = generate(option = "drop", std = "true", valid = "false", test = "false");
best_mach = machine(LogisticClassifier(penalty = :l2, lambda = rep_drop_L2.best_model.lambda), select(train[:,:], Not(:precipitation_nextday)), train.precipitation_nextday)|> fit!;
sample = CSV.read(joinpath(@__DIR__, "data", "sample_submission.csv"), DataFrame);
pred = pdf.(predict(best_mach, test), true);
sample.precipitation_nextday = pred;
CSV.write(joinpath(@__DIR__, "results", "Linear_submission.csv"), sample);

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
using Documenter, InteractBase, Literate
src = joinpath(@__DIR__, "src")
Literate.markdown(joinpath(src, "tutorial.jl"), src, codefence = "```julia" => "```")

makedocs(
    format = :html,
    sitename = "InteractBase",
    authors = "Pietro Vertechi",
    pages = [
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "API reference" => "api_reference.md",
    ]
)

deploydocs(
    repo = "github.com/piever/InteractBase.jl.git",
    target = "build",
    julia  = "0.6",
    osname = "linux",
    deps   = nothing,
    make   = nothing
)


text:
module PadeApproximation

using ..Polynomials
indeterminate = Polynomials.indeterminate

using ..PolyCompat
export Pade, padeval

#=
Pade approximation

Note: This can be moved to polynomials/Polynomial.jl after Poly type is removed
=#

export Pade

"""
    Pade(::Polynomial, m::Integer, n::Integer)
    Pade(::Polynomial, ::Polynomial)

Return Pade approximation of polynomial.

# References
[Pade Approximant](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant)
"""
struct Pade{T <: Number,S <: Number}
    p::Union{Poly{T}, Polynomial{T}}
    q::Union{Poly{S}, Polynomial{S}}
    var::Symbol
    function Pade{T,S}(p::Union{Poly{T}, Polynomial{T}}, q::Union{Poly{S}, Polynomial{S}}) where {T,S}
        if indeterminate(p) != indeterminate(q) error("Polynomials must have same variable") end
        new{T,S}(p, q, indeterminate(p))
    end
end

Pade(p::Polynomial{T}, q::Polynomial{S}) where {T <: Number,S <: Number} = Pade{T,S}(p, q)

function Pade(c::Polynomial{T}, m::Integer, n::Integer) where {T}
    m + n < length(c) || error("m + n must be less than the length of the Polynomial")
    rold = Polynomial([zeros(T, m + n + 1);one(T)], indeterminate(c))
    rnew = Polynomial(c[0:m + n], indeterminate(c))
    uold = Polynomial([one(T)], indeterminate(c))
    vold = Polynomial([zero(T)], indeterminate(c))
    unew, vnew = vold, uold
    @inbounds for i = 1:n
        temp0, temp1, temp2 = rnew, unew, vnew
        q, rnew = divrem(rold, rnew)
        unew, vnew = uold - q * unew, vold - q * vnew
        rold, uold, vold = temp0, temp1, temp2

    end
    if vnew[0] == 0
        d = gcd(rnew, vnew)
        rnew = rnew ÷ d
        vnew = vnew ÷ d
    end
    Pade(rnew / vnew[0], vnew / vnew[0])
end

Pade(p::Poly{T}, q::Poly{S}) where {T <: Number,S <: Number} = Pade{T,S}(p, q)

function Pade(c::Poly{T}, m::Integer, n::Integer) where {T}
    m + n < length(c) || error("m + n must be less than the length of the polynomial")
    rold = Poly([zeros(T, m + n + 1);one(T)], indeterminate(c))
    rnew = Poly(c[0:m + n], indeterminate(c))
    uold = Poly([one(T)], indeterminate(c))
    vold = Poly([zero(T)], indeterminate(c))
    unew, vnew = vold, uold
    @inbounds for i = 1:n
        temp0, temp1, temp2 = rnew, unew, vnew
        q, rnew = divrem(rold, rnew)
        unew, vnew = uold - q * unew, vold - q * vnew
        rold, uold, vold = temp0, temp1, temp2

    end
    if vnew[0] == 0
        d = gcd(rnew, vnew)
        rnew = rnew ÷ d
        vnew = vnew ÷ d
    end
    Pade(rnew / vnew[0], vnew / vnew[0])
end

"""
    (::Pade)(x)

Evaluate the Pade approximant at the given point.

# Examples
```jldoctest pade
julia> using Polynomials, Polynomials.PolyCompat, SpecialFunctions

julia> p = Polynomial(@.(1 // BigInt(gamma(1:17))));

julia> pade = Pade(p, 8, 8);

julia> pade(1.0) ≈ exp(1.0)
true

```
"""
(PQ::Pade)(x) = PQ.p(x) / PQ.q(x)

## Compat
padeval(PQ::Pade, x::Number) = PQ(x)
padeval(PQ::Pade, x) = PQ.(x)

end


text:
using MultivariateStats
using Base.Test

import MultivariateStats: qnormalize!

srand(34568)

dx = 5
dy = 6
p = 3

# CCA with zero means

X = rand(dx, 100)
Y = rand(dy, 100)

Px = qr(randn(dx, dx))[1][:, 1:p]
Py = qr(randn(dy, dy))[1][:, 1:p]

M = CCA(Float64[], Float64[], Px, Py, [0.8, 0.6, 0.4])

@test xindim(M) == dx
@test yindim(M) == dy
@test xmean(M) == zeros(dx)
@test ymean(M) == zeros(dy)
@test xprojection(M) == Px
@test yprojection(M) == Py
@test correlations(M) == [0.8, 0.6, 0.4]

@test_approx_eq xtransform(M, X) Px'X
@test_approx_eq ytransform(M, Y) Py'Y

## CCA with nonzero means

ux = randn(dx)
uy = randn(dy)

M = CCA(ux, uy, Px, Py, [0.8, 0.6, 0.4])

@test xindim(M) == dx
@test yindim(M) == dy
@test xmean(M) == ux
@test ymean(M) == uy
@test xprojection(M) == Px
@test yprojection(M) == Py
@test correlations(M) == [0.8, 0.6, 0.4]

@test_approx_eq xtransform(M, X) Px' * (X .- ux)
@test_approx_eq ytransform(M, Y) Py' * (Y .- uy)

## prepare data

n = 1000
dg = 10
G = randn(dg, n)

X = randn(dx, dg) * G + 0.2 * randn(dx, n)
Y = randn(dy, dg) * G + 0.2 * randn(dy, n)
xm = vec(mean(X, 2))
ym = vec(mean(Y, 2))
Zx = X .- xm
Zy = Y .- ym

Cxx = cov(X; vardim=2)
Cyy = cov(Y; vardim=2)
Cxy = cov(X, Y; vardim=2)
Cyx = Cxy'

## ccacov

# X ~ Y
M = fit(CCA, X, Y; method=:cov, outdim=p)
Px = xprojection(M)
Py = yprojection(M)
rho = correlations(M)
@test xindim(M) == dx
@test yindim(M) == dy
@test outdim(M) == p
@test xmean(M) == xm
@test ymean(M) == ym
@test issorted(rho; rev=true)

@test_approx_eq Px' * Cxx * Px eye(p)
@test_approx_eq Py' * Cyy * Py eye(p)
@test_approx_eq Cxy * (Cyy \ Cyx) * Px scale(Cxx * Px, rho.^2)
@test_approx_eq Cyx * (Cxx \ Cxy) * Py scale(Cyy * Py, rho.^2)
@test_approx_eq Px qnormalize!(Cxx \ (Cxy * Py), Cxx)
@test_approx_eq Py qnormalize!(Cyy \ (Cyx * Px), Cyy)

# Y ~ X
M = fit(CCA, Y, X; method=:cov, outdim=p)
Py = xprojection(M)
Px = yprojection(M)
rho = correlations(M)
@test xindim(M) == dy
@test yindim(M) == dx
@test outdim(M) == p
@test xmean(M) == ym
@test ymean(M) == xm
@test issorted(rho; rev=true)

@test_approx_eq Px' * Cxx * Px eye(p)
@test_approx_eq Py' * Cyy * Py eye(p)
@test_approx_eq Cxy * (Cyy \ Cyx) * Px scale(Cxx * Px, rho.^2)
@test_approx_eq Cyx * (Cxx \ Cxy) * Py scale(Cyy * Py, rho.^2)
@test_approx_eq Px qnormalize!(Cxx \ (Cxy * Py), Cxx)
@test_approx_eq Py qnormalize!(Cyy \ (Cyx * Px), Cyy)

## ccasvd

# n > d
M = fit(CCA, X, Y; method=:svd, outdim=p)
Px = xprojection(M)
Py = yprojection(M)
rho = correlations(M)
@test xindim(M) == dx
@test yindim(M) == dy
@test outdim(M) == p
@test xmean(M) == xm
@test ymean(M) == ym
@test issorted(rho; rev=true)

@test_approx_eq Px' * Cxx * Px eye(p)
@test_approx_eq Py' * Cyy * Py eye(p)
@test_approx_eq Cxy * (Cyy \ Cyx) * Px scale(Cxx * Px, rho.^2)
@test_approx_eq Cyx * (Cxx \ Cxy) * Py scale(Cyy * Py, rho.^2)
@test_approx_eq Px qnormalize!(Cxx \ (Cxy * Py), Cxx)
@test_approx_eq Py qnormalize!(Cyy \ (Cyx * Px), Cyy)


text:
# A `Results` object stores information about the results of benchmarking an
# expression.
#
# Fields:
#
#     precompiled::Bool: During benchmarking, did we ensure that the
#         benchmarkable function was precompiled? We do this for all
#         functions that can be executed at least twice without exceeding
#         the time budget specified by the user. Otherwise, we only measure
#         the expression's execution once and store this flag to indicate that
#         our single measurement potentially includes compilation time.
#
#     multiple_samples::Bool: During benchmarking, did we gather more than one
#         sample? If so, we will attempt to report results that acknowledge
#         the variability in our sampled observations. If not, we will only
#         report a single point estimate of the expression's performance
#         without any measure of our uncertainty in that point estimate.
#
#     search_performed::Bool: During benchmarking, did we perform a geometric
#         search to determine the minimum number of times we must evaluate the
#         expression being benchmarked before an individual sample can be
#         considered an unbiased estimate of the expression's performance? If
#         so, downstream analyses should use the slope of the linear regression
#         model, `elapsed_time ~ 1 + evaluations`, as their estimate of the
#         time it takes to evaluate the expression once. If not, we know that
#         `evaluation[i] == 1` for all `i`.
#
#     samples::Samples: A record of all samples that were recorded during
#         benchmarking.
#
#     time_used::Float64: The time (in nanoseconds) that was consumed by the
#         benchmarking process.

immutable Results
    precompiled::Bool
    multiple_samples::Bool
    search_performed::Bool
    samples::Samples
    time_used::Float64
end

# A `SummaryStatistics` object stores the results of a statistic analysis of
# a `Results` object. The precise analysis strategy employed depends on the
# structure of the `Results` object:
#
#     (1) If only a single sample of a single evaluation was recorded, the
#     analysis reports only point estimates.
#
#     (2) If multiple samples of a single evaluation were recorded, the
#     analysis reports point estimates and CI's determined by straight
#     summary statistic calculations.
#
#     (3) If a geometric search was performed to generate samples that
#     represent multiple evaluations, an OLS regression is fit that
#     estimates the model `elapsed_time ~ 1 + evaluations`. The slope
#     of the `evaluations` term is treated as the best estimate of the
#     elapsed time of a single evaluation.
#
# For both strategies (2) and (3), we try to make up for a lack of IID
# samples by using 6-sigma CI's instead of the traditional 2-sigma CI'
# reported in most applied statistical work.
#
# In order to estimate GC time, we assume that the relationship betweeen GC
# time and total time is constant with respect to the number of evaluations.
# As such, we avoid using an OLS fit for estimating GC time.
#
# We also assume that the ratio, `bytes_allocated / evaluations` is a
# constant that only exhibits upward-biased noise. As such, we take the
# minimum value of this ratio to determine the memory allocation behavior of
# an expression.

immutable SummaryStatistics
    n::Int
    n_evaluations::Int
    elapsed_time_lower::Nullable{Float64}
    elapsed_time_center::Float64
    elapsed_time_upper::Nullable{Float64}
    gc_proportion_lower::Nullable{Float64}
    gc_proportion_center::Float64
    gc_proportion_upper::Nullable{Float64}
    bytes_allocated::Int
    allocations::Int
    r²::Nullable{Float64}

    function SummaryStatistics(r::Results)
        s = r.samples
        n = length(s.evaluations)
        n_evaluations = convert(Int, sum(s.evaluations))
        if !r.search_performed
            if !r.multiple_samples
                @assert n == 1
                @assert all(s.evaluations .== 1.0)
                m = s.elapsed_times[1]
                gc_proportion = s.gc_times[1] / s.elapsed_times[1]
                elapsed_time_center = m
                elapsed_time_lower = Nullable{Float64}()
                elapsed_time_upper = Nullable{Float64}()
                r² = Nullable{Float64}()
                gc_proportion_center = 100.0 * gc_proportion
                gc_proportion_lower = Nullable{Float64}()
                gc_proportion_upper = Nullable{Float64}()
            else
                @assert all(s.evaluations .== 1.0)
                m = mean(s.elapsed_times)
                sem = std(s.elapsed_times) / sqrt(n)
                gc_proportion = mean(s.gc_times ./ s.elapsed_times)
                gc_proportion_sem = std(s.gc_times ./ s.elapsed_times) / sqrt(n)
                r² = Nullable{Float64}()
                elapsed_time_center = m
                elapsed_time_lower = max(0.0, m - 6.0 * sem)
                elapsed_time_upper = m + 6.0 * sem
                gc_proportion_center = 100.0 * gc_proportion
                gc_proportion_lower = Nullable{Float64}(
                    max(
                        0.0,
                        gc_proportion_center - 6.0 * 100 * gc_proportion_sem
                    )
                )
                gc_proportion_upper = Nullable{Float64}(
                    min(
                        100.0,
                        gc_proportion_center + 6.0 * 100 * gc_proportion_sem
                    )
                )
            end
        else
            a, b, ols_r² = ols(s.evaluations, s.elapsed_times)
            sem = sem_ols(s.evaluations, s.elapsed_times)
            gc_proportion = mean(s.gc_times ./ s.elapsed_times)
            gc_proportion_sem = std(s.gc_times ./ s.elapsed_times) / sqrt(n)
            r² = Nullable{Float64}(ols_r²)
            elapsed_time_center = b
            elapsed_time_lower = max(0.0, b - 6.0 * sem)
            elapsed_time_upper = b + 6.0 * sem
            gc_proportion_center = 100.0 * gc_proportion
            gc_proportion_lower = Nullable{Float64}(
                max(
                    0.0,
                    gc_proportion_center - 6.0 * 100 * gc_proportion_sem
                )
            )
            gc_proportion_upper = Nullable{Float64}(
                min(
                    100.0,
                    gc_proportion_center + 6.0 * 100 * gc_proportion_sem
                )
            )
        end

        i = indmin(s.bytes_allocated ./ s.evaluations)

        bytes_allocated = fld(
            s.bytes_allocated[i],
            convert(UInt, s.evaluations[i])
        )
        allocations = fld(
            s.allocations[i],
            convert(UInt, s.evaluations[i])
        )

        new(
            n,
            n_evaluations,
            elapsed_time_lower,
            elapsed_time_center,
            elapsed_time_upper,
            gc_proportion_lower,
            gc_proportion_center,
            gc_proportion_upper,
            bytes_allocated,
            allocations,
            r²,
        )
    end
end

function pretty_time_string(t)
    if t < 1_000.0
        @sprintf("%.2f ns", t)
    elseif t < 1_000_000.0
        @sprintf("%.2f μs", t / 1_000.0)
    elseif t < 1_000_000_000.0
        @sprintf("%.2f ms", t / 1_000_000.0)
    else # if t < 1_000_000_000_000.0
        @sprintf("%.2f s", t / 1_000_000_000.0)
    end
end

function pretty_memory_string(b)
    if b < 1_024.0
        @sprintf("%.2f bytes", b)
    elseif b < 1_024.0^2
        @sprintf("%.2f kb", b / 1_024.0)
    elseif b < 1_024.0^3
        @sprintf("%.2f mb", b / 1_024.0^2)
    else # if b < 1_024.0^4
        @sprintf("%.2f gb", b / 1_024.0^3)
    end
end

# Pretty-print information about the results of benchmarking an expression.
#
# Arguments:
#
#     io::IO: An `IO` object to be written to.
#
#     r::Results: The `Results` object that we want to print to `io`.
#

function Base.show(io::IO, r::Results)
    stats = SummaryStatistics(r)

    max_length = 24
    @printf(io, "================ Benchmark Results ========================\n")

    if !r.precompiled
        @printf(io, "Warning: function may not have been precompiled\n")
    end
    if isnull(stats.elapsed_time_lower) || isnull(stats.elapsed_time_upper)
        @printf(
            io,
            "%s: %s\n",
            lpad("Time per evaluation", max_length),
            pretty_time_string(stats.elapsed_time_center),
        )
    else
        @printf(
            io,
            "%s: %s [%s, %s]\n",
            lpad("Time per evaluation", max_length),
            pretty_time_string(stats.elapsed_time_center),
            pretty_time_string(get(stats.elapsed_time_lower)),
            pretty_time_string(get(stats.elapsed_time_upper)),
        )
    end
    if isnull(stats.elapsed_time_lower) || isnull(stats.elapsed_time_upper)
        @printf(
            io,
            "%s: %.2f%%\n",
            lpad("Proportion of time in GC", max_length),
            stats.gc_proportion_center
        )
    else
        @printf(
            io,
            "%s: %.2f%% [%.2f%%, %.2f%%]\n",
            lpad("Proportion of time in GC", max_length),
            stats.gc_proportion_center,
            get(stats.gc_proportion_lower),
            get(stats.gc_proportion_upper),
        )
    end
    @printf(
        io,
        "%s: %s\n",
        lpad("Memory allocated", max_length),
        pretty_memory_string(stats.bytes_allocated),
    )
    @printf(
        io,
        "%s: %d allocations\n",
        lpad("Number of allocations", max_length),
        stats.allocations,
    )
    @printf(
        io,
        "%s: %d\n",
        lpad("Number of samples", max_length),
        stats.n
    )
    @printf(
        io,
        "%s: %d\n",
        lpad("Number of evaluations", max_length),
        stats.n_evaluations
    )
    if r.search_performed
        @printf(
            io,
            "%s: %.3f\n",
            lpad("R² of OLS model", max_length),
            get(stats.r², NaN),
        )
    end
    @printf(
        io,
        "%s: %.2f s\n",
        lpad("Time spent benchmarking", max_length),
        r.time_used,
    )
    # @printf(
    #     io,
    #     "%s: %s\n",
    #     lpad("Precompiled", max_length),
    #     string(r.precompiled)
    # )
    # @printf(
    #     io,
    #     "%s: %s\n",
    #     lpad("Multiple samples", max_length),
    #     string(r.multiple_samples),
    # )
    # @printf(
    #     io,
    #     "%s: %s",
    #     lpad("Search performed", max_length),
    #     string(r.search_performed),
    # )
end


text:
# Module-related types and auxiliary functions

export
    CuModule, CuModuleFile

include(joinpath("module", "jit.jl"))

const CuModule_t = Ptr{Cvoid}

"""
    CuModule(data, options::Dict{CUjit_option,Any})
    CuModuleFile(path, options::Dict{CUjit_option,Any})

Create a CUDA module from a data, or a file containing data. The data may be PTX code, a
CUBIN, or a FATBIN.

The `options` is an optional dictionary of JIT options and their respective value.
"""
mutable struct CuModule
    handle::CuModule_t
    ctx::CuContext

    function CuModule(data, options::Dict{CUjit_option,Any}=Dict{CUjit_option,Any}())
        handle_ref = Ref{CuModule_t}()

        options[ERROR_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
        @debug begin
            options[INFO_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
            options[LOG_VERBOSE] = true
            "JIT compiling code" # FIXME: remove this useless message
        end
        optionKeys, optionVals = encode(options)

        try
            @apicall(:cuModuleLoadDataEx,
                     (Ptr{CuModule_t}, Ptr{Cchar}, Cuint, Ptr{CUjit_option}, Ptr{Ptr{Cvoid}}),
                     handle_ref, data, length(optionKeys), optionKeys, optionVals)
        catch err
            (err == ERROR_NO_BINARY_FOR_GPU || err == ERROR_INVALID_IMAGE || err == ERROR_INVALID_PTX) || rethrow(err)
            options = decode(optionKeys, optionVals)
            rethrow(CuError(err.code, options[ERROR_LOG_BUFFER]))
        end

        @debug begin
            options = decode(optionKeys, optionVals)
            if isempty(options[INFO_LOG_BUFFER])
                """JIT info log is empty"""
            else
                """JIT info log:
                   $(options[INFO_LOG_BUFFER])"""
            end
        end

        ctx = CuCurrentContext()
        obj = new(handle_ref[], ctx)
        finalizer(unsafe_unload!, obj)
        return obj
    end
end

function unsafe_unload!(mod::CuModule)
    if isvalid(mod.ctx)
        @apicall(:cuModuleUnload, (CuModule_t,), mod)
    end
end

Base.unsafe_convert(::Type{CuModule_t}, mod::CuModule) = mod.handle

Base.:(==)(a::CuModule, b::CuModule) = a.handle == b.handle
Base.hash(mod::CuModule, h::UInt) = hash(mod.handle, h)

CuModuleFile(path) = CuModule(read(path, String))

include(joinpath("module", "function.jl"))
include(joinpath("module", "global.jl"))
include(joinpath("module", "linker.jl"))