text:
using Clang.cindex
using Clang.wrap_c
using Compat

POLYLIB_INCLUDE = "./polylib-5.22.5/include"

petsc_header = [joinpath(POLYLIB_INCLUDE, "polylib/polylib64.h")]

# Set up include paths
clang_includes = ASCIIString[]
push!(clang_includes, "/usr/lib/llvm-3.4/lib/clang/3.4/include")
push!(clang_includes, POLYLIB_INCLUDE)

# Clang arguments
clang_extraargs = ["-v"]
# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

# Callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
    return startswith(dirname(cursor_header), POLYLIB_INCLUDE)
end

lib_file(hdr::ASCIIString) = "polylib"
output_file(hdr::ASCIIString) = "PolyLib.jl"

const wc = wrap_c.init(;
                        headers = petsc_header,
                        output_file = "libpolylib64_h.jl",
                        common_file = "libpolylib64_common.jl",
                        clang_includes      = clang_includes,
                        clang_args          = clang_extraargs,
                        header_wrapped      = wrap_header,
                        header_library      = lib_file,
                        header_outputfile   = output_file,
                        clang_diagnostics = true)

run(wc)


text:
module BookTemplate

using Reexport: @reexport
@reexport begin
using Books:
    build_all,
    gen
using DataFrames:
    DataFrame,
    filter!,
    filter,
    select!,
    select
end # @reexport

include("data.jl")
export example_dataframe

"""
    build()

This function is called during CI.
"""
function build()
    println("Building Example Book")
    # To avoid publishing broken websites.
    fail_on_error = true
    gen(; fail_on_error)
    build_all(; fail_on_error)
end

end # module

ext: jl

max_stars_repo_name: tbeason/InvestmentsClass

max_issues_repo_name: tbeason/InvestmentsClass

max_forks_repo_name: tbeason/InvestmentsClass


text:
module RawFile

export saveraw,readraw,rawsize,appendraw
export RawFile,RawFileIter,start,done,next

token = "RAWF"
version = 2
endtoken = "FIN"

mutable struct RawHeader
    version::UInt8
    eltype::Type
    sizes::Vector{Int64}
end

import Base.write

function write(f::IO,x::RawHeader)
    write(f,token)
    write(f,x.version)
    typename = string(x.eltype)
    write(f,UInt8(length(typename)))
    write(f,typename)
    write(f,UInt8(length(x.sizes)))
    write(f,x.sizes)
    return nothing
end

function readheader(f::IO)
    tok = String(read(f,length(token)))
    tok != token && error("Invalid token in RawFile")
    ver = read(f,UInt8)
    ver != version && warn("Incorrect RawFile version ($ver, current version $version)")
    lentypename = read(f,UInt8)
    typename = String(read(f,lentypename))
    eltype = eval(Meta.parse(typename))
    nd = read(f,UInt8)
    sizes = read!(f,Array{Int64}(undef,nd))
    RawHeader(ver,eltype,sizes)
end

"""
    saveraw{T<:Number,V}(a::AbstractArray{T,V},f::IO)

Save `a` to the open IO stream `f`.
"""
function saveraw(a::AbstractArray{T,V},f::IO) where {T<:Number,V}
    header = RawHeader(version,eltype(a),collect(size(a)))
    write(f,header)
    write(f,a)
    write(f,endtoken)
    return nothing
end

"""
    saveraw{T<:Number,V}(a::AbstractArray{T,V},fname::String)

Save `a` to the file `fname`.
"""
function saveraw(a::AbstractArray{T,V},fname::String) where {T<:Number,V}
    open(fname,"w") do f
        saveraw(a,f)
    end
    return nothing
end

"""
    appendraw{T<:Number,V}(a::AbstractArray{T,V},fname::String)

Append the `AbstractArray` `a` to the file `fname`, along last dimension. Requires `a` to be the same
`Type` and shape (excluding last dimension)
"""
function appendraw(a::AbstractArray{T,V},fname::String) where {T<:Number,V}
    open(fname,"r+") do f
        h = readheader(f)
        h.eltype != eltype(a) && error("Trying to append RawFile with different data Type")
        Tuple(h.sizes[1:end-1]) != size(a)[1:end-1] && error("Trying to append RawFile with Array of different shape")
        h.sizes[end] += size(a)[end]
        seekstart(f)
        write(f,h)
        seekend(f)
        seek(f,position(f)-length(endtoken))
        write(f,a)
        write(f,endtoken)
    end
    return nothing
end

appendraw(a::T,fname::String) where {T<:Number} = appendraw([a],fname)

"""
    readraw(f::IO)

Read the open IO stream `f` and return a reconstructed `Array` with the proper `Type` and size.
"""
function readraw(f::IO)
    h = readheader(f)
    d = read!(f,Array{h.eltype}(undef,Tuple(h.sizes)))
    endtok = String(read(f,length(endtoken)))
    endtok != endtoken && error("Invalid end of RawFile")
    return d
end

"""
    readraw(fname::String)

Read the file `fname` and return a reconstructed `Array` with the proper `Type` and size.
"""
function readraw(fname::String)
    open(fname) do f
        return readraw(f)
    end
end

"""
    rawsize(f::IO)

Read size from IO stream `f`, returns `Tuple`
"""
function rawsize(f::IO)
    h = readheader(f)
    return Tuple(h.sizes)
end

"""
    rawsize(fname::String)

Read size from file `fname`, returns `Tuple`
"""
function rawsize(fname::String)
    open(fname) do f
        return rawsize(f)
    end
end

mutable struct PartialRaw
    f::IO
    eltype::Type
    sizes::Array{Int64,1}
    total::Int
end

"""
    saveraw(func::Function,fname::String)

Save `Array`s progressively produced in `func`. Each partial `Array` needs to have
the same dimensions (excluding the last dimension), and are concatenated along the
last dimension.

# Examples
```julia-repl
julia> saveraw("test.raw") do f
           for i=1:10
               write(f,rand(100,10,5))
           end
       end

julia> rawsize("test.raw")
(100, 10, 50)
```
"""
function saveraw(func::Function,fname::String)
    open(fname,"w") do f
        p = PartialRaw(f,Int,[],0)
        func(p)
        write(f,endtoken)
        seekstart(f)
        p.sizes[end] = p.total
        write(f,RawHeader(version,p.eltype,p.sizes))
    end
    return nothing
end

function write(p::PartialRaw,d::AbstractArray{T,V}) where {T<:Number,V}
    if length(p.sizes)==0
        # This is the first Array we've seen
        p.sizes = collect(size(d))
        p.eltype = eltype(d)
        write(p.f,RawHeader(version,p.eltype,p.sizes))
    else
        Tuple(p.sizes[1:end-1]) != size(d)[1:end-1]  && error("Cannot write partial RawFile, all dimensions other than the last must be the same")
        p.eltype != eltype(d) && error("Cannot write partial RawFile, all data must be the same type")
    end
    write(p.f,d)
    p.total += size(d)[end]
    return nothing
end

"""
    RawFileIter(fname::String,num_batch::Int)

Creates an interator that will iterate through the file `fname`, returning `Array`s
with the last dimension `num_batch` (or less if EOF is reached).

# Examples
```julia-repl
julia> for d in RawFileIter("test.raw",20)
                  @info(size(d))
              end
[ Info: (100, 10, 20)
[ Info: (100, 10, 20)
[ Info: (100, 10, 10)
```
"""
mutable struct RawFileIter
    fname::String
    num_batch::Int
end

mutable struct RawFileState
    i::Int
    total_length::Int
    num_batch::Int
    batch_step::Int
    batch_size::Array{Int,1}
    f::IO
    eltype::Type
end

import Base.iterate,Base.length

function finalize(s::RawFileState)
    close(s.f)
end

function iterate(r::RawFileIter)
    f = open(r.fname)
    h = readheader(f)
    batch_step = reduce(*,h.sizes[1:end-1])
    total_length = batch_step*h.sizes[end]
    batch_size = copy(h.sizes)
    i = 0
    s = RawFileState(i,total_length,r.num_batch,batch_step,batch_size,f,h.eltype)
    finalizer(finalize,s)

    return iterate(r,s)
end

function iterate(r::RawFileIter,state::RawFileState)
    if state.i < state.total_length
        this_len = Int(min(state.num_batch,(state.total_length-state.i)/state.batch_step))
        state.batch_size[end] = this_len
        d = read!(state.f,Array{state.eltype}(undef,Tuple(state.batch_size)))
        state.i += state.batch_step*this_len
        return (d,state)
    else
        endtok = String(read(state.f,length(endtoken)))
        endtok != endtoken && error("Invalid end of RawFile")
        return nothing
    end
end

"""
    readraw(func::Function,fname::String,batch::Int)

Convenience function for iteratively reading a `RawFile` in chunks. Calls `RawFileIter`
and runs `func` on each chunk.

# Examples
```julia-repl
julia> readraw("test.raw",20) do c
           @info(size(c))
       end
[ Info: (100, 10, 20)
[ Info: (100, 10, 20)
[ Info: (100, 10, 10)
```
"""
function readraw(func::Function,fname::String,batch::Int)
    for d in RawFileIter(fname,batch)
        func(d)
    end
end

end


text:
# Greetings! from path ~/.julia/config/startup.jl

############################################
## 1.Set parameters
############################################
# processes and threads parameters
cores = 2                                                           # Int64; total number of processes. cores number should be less than time-chunck number corresponded with "select_time".
threads = 1                                                         # Int64; total number of threads. threads number should be less than N, N is the number of sta-pairs in each time-chunck.
flag = true                                                         # true of false; print computing time for debugging purpose.
precompile = false                                                  # true or false; precompile or not.

# absolute path parameters
FFTDIR = expanduser("/Users/yf/3.Project/2-1.Velocity_Change/FFT")  # dir where FFT data are stored.
CCFDIR = expanduser("./CORR_fft")                                   # dir to store CC data.

# select time-chunck folder for cross-correlate
select_time = false                                                 # true or false; whether specify time-chunck folder path.
TIMETXT = expanduser("./time_cc.txt")                               # time-chunck folder path.

# select station for cross
select_sta = "sr"                                                   # "no", "sta" or "sr"; "sta" cross-correlate only works for selected net-sta pairs base on STATXT; "sr" cross-correlate only works for specified source and receiver based on SRTXT.
STATXT = expanduser("./sta_cc.txt")                                 # station infomation including: network-name and station-name.
SRTXT = expanduser("./sr_cc.txt")                                   # source and receiver information including: source-1 receiver-1 r-2 r-3...\nsource-2 r-1 r-2 r-3...
channel_regular = "*CHN.*"                                          # channel(regular expression) of FFT-data; we will use 'glob(comp,FFTDIR)' to select data in 'FFTDIR' path.
comp = "CHN"                                                        # component name corresponds to "channel_regular"; we will use 'load_fft(file,comp)' to read FFTData in 'FFTDIR' path, and "comp" is the group name of jld2 file(HDF5 in julia).

# some control parameters
cc_method = "CC"                                                    # "CC" and "PCC"; "CC" for pure cross correlation, or "PCC" for Phase cross-correlation, see [Ventosa et al., 2019].
auto_corr = true                                                    # true or false; perform auto-correlation or not.
cross_corr = true                                                   # true or false; perform cross-correlation or not.

# cross-correlation parameters
maxlag = 20.0                                                       # Float64; lags of cross-correlation to save (sec).
substack = true                                                     # true or false; sub-stack one time-chunck result of cross-correlation or not.
substack_method = ["mean","pws"]                                           # "mean", "pws", and "robust"; Vector{String}, you can choose one or more like ["mean","pws","robust"].

############################################
# 2.run NoiseCC.FFT_CC
############################################
# a.using and import
# using Distributed
# NoiseCC_path = "/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/"
# push!(LOAD_PATH, NoiseCC_path)
# @everywhere using NoiseCC # or @everywhere import NoiseCC

# b.include
using Distributed
addprocs(cores, topology = :master_worker, exeflags = "-t $(threads)")
@everywhere include(
    "/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/NoiseCC.jl",
)

NoiseCC.CC(
    cores,
    threads,
    flag,
    precompile,
    FFTDIR,
    CCFDIR,
    select_time,
    TIMETXT,
    select_sta,
    STATXT,
    SRTXT,
    channel_regular,
    comp,
    cc_method,
    auto_corr,
    cross_corr,
    maxlag,
    substack,
    substack_method,
)

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
using MultiJuMP, JuMP
using Test
#using AmplNLWriter
using Ipopt

# facts("Utopia and Nadir points") do
#     # Create model
#     m = MultiModel()
#     @variable(m, -1 <= x <= 1)
#     @NLexpression(m, f1, x^2)
#     @NLexpression(m, f2, (x-0.5)^2)

#     multim = getMultiData(m)
#     multim.f1 = f1
#     multim.f2 = f2

#     @fact solve(m, method = :WS) --> :Optimal
#     @fact multim.utopia --> roughly([0.0, 0.0], 1e-5)
#     @fact multim.nadir --> roughly([0.25, 0.25], 1e-5)
# end

@testset "NBI optimisation" begin
    m = MultiModel(solver = IpoptSolver(print_level=0))
    @variable(m, x[i=1:5])
    @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = getMultiData(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    solve(m, method = :NBI)

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0, 7.16982, 4.48665, 2.08008]
    f2true = [2.13057, -4.01115, -2.78066, -1.45458, 0.0513689]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "WS optimisation" begin
    m = MultiModel(solver = IpoptSolver(print_level=0))
    @variable(m, x[i=1:5])
    @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = getMultiData(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    solve(m, method = :WS)

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0, 10.0, 2.88267, 0.726877]
    f2true = [2.13057, -4.01115, -4.01115, -0.510963, 1.48613]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "EPS optimisation" begin
    m = MultiModel(solver = IpoptSolver(print_level=0))
    @variable(m, x[i=1:5])
        @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = getMultiData(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    solve(m, method = :EPS)

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0, 7.63877, 5.27754, 2.91631]
    f2true = [2.13057, -4.01115, -2.99356, -1.86918, -0.532781]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "TODO: Initial value test" begin
    # TODO: Let f1 have two minima and
    # check that initial value gives correct minima
    # TODO: Add anonymous variables
    @test 1 == 1
end

@testset "TODO: NBI Inequality test" begin
    # TODO: Test that the inequality-extension of NBI works
    # We can use the nc_paper_2.jl example
    @test 1 == 1
end


text:
# Packages required for testing
using Test
using Random
using Logging

# Package Under Test
using JuliaPackageTemplate

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Fix randomness during tests
Random.seed!(0)

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test isapprox(a,b, rtol=rtol, atol=atol)
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        @test isapprox(a, y, rtol=rtol, atol=atol)
    end
end

# Define package tests
@time @testset "JuliaPackageTemplate Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "JuliaPackageTemplate.YourSubmodule" begin
        include(joinpath(testdir, "test_submodule.jl"))
    end
    @time @testset "JuliaPackageTemplate.SingleSatellite" begin
        include(joinpath(testdir, "test_rubberducks.jl"))
    end
end


text:
using ProjectedGPE
using Base.Test

# Run tests
tic()
include("modeorthonorm.jl")
@testset "Hermite basis tests" begin include("testHermites.jl") end
toc()


text:
import YAML
config = YAML.load(open("config.yaml"))

# Read the array from a text file

# Instead of using a Latin Hypercube, write out parameters in several dimensions, starting from a central point and then doing ranges. 5 points in each sample.

nus = linspace(345.7956, 345.7975, 5)
masses = linspace(0.8, 1.3, 5)
r_cs = linspace(300, 600, 5)
T_10s = linspace(60, 130, 5)
qs = linspace(0.3, 0.7, 5)
logM_gass = linspace(-3, -1.5, 5)
ksis = linspace(0.2, 0.5, 5)
incls = linspace(35, 55, 5)
PAs = linspace(145, 160, 5)

np = 9

total = 5 * np

pars = Array(Float64, (np, total))

parranges = Float64[nus masses r_cs T_10s qs logM_gass ksis incls PAs]

println(parranges)

centers = Int64[3, 3, 3, 3, 3, 3, 3, 3]

colinds = Int64[i for i=1:np]
println(colinds)

for col=1:np
  for row=1:5
    k = (col - 1) * 5 + row

    # insert!(copy(centers)), index, item)
    ind = copy(centers)
    insert!(ind, col, row)

    parsub = Array(Float64, np)

    for (l, pair) in enumerate(zip(ind, colinds))
      parsub[l] = parranges[pair...]
    end

    pars[:,k] = parsub

  end
end

println(size(pars'))
println(pars')

# Now, write this parameter list to disk
writedlm("parameters_star.dat", pars')

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
#Numbers for 87Rb
#23/07/2017

function Rb87numbers(nn,ll,jj)
    #Low level parameters
    D1trans = 1.56051536 #D1 line transition energy in eV
    D2trans = 1.58999085 #D2 line transition energy in eV
    n4l2j5_2en = 2.40116159
    n4l2j3_2en = 2.40121692
    n6l0j1_2en = 2.49759249
    n6l1j1_2en = 2.94203794
    n6l1j3_2en = 2.95165365
    n5l0j1_2en = 0
    ionlim = 4.1771270 #Rb groundstate energy in eV
    #Other parameters
    ZZ = 37 #Nucler charge of Rubidium
    spin = 1/2 #Electron spin
    #Quantum defects from Li et al. (2003) and Han et al. (2006)
    if ll==0
        delta_0 = 3.1311804
        delta_2 = 0.1784
    elseif ll==1&&jj==1/2
        delta_0 = 2.65448849
        delta_2 = 0.2900
    elseif ll==1&&jj==3/2
        delta_0 = 2.6416737
        delta_2 = 0.2950
    elseif ll==2&&jj==3/2
        delta_0 = 1.34809171
        delta_2 = -0.60286
    elseif ll==2&&jj==5/2
        delta_0 = 1.34646572
        delta_2 = -0.59600
    elseif ll==3&&jj==5/2
        delta_0 = 0.0165192
        delta_2 = -0.085
    elseif ll==3&&jj==7/2
        delta_0 = 0.0165437
        delta_2 = -0.086
    else
        delta_0 = 0
        delta_2 = 0
    end
    #Approximation for quantum efect assuming n > 20
    delta_nlj = delta_0 + delta_2/(nn-delta_0)^2

    #Now get the energy eigenvalue
    if nn==5&&ll==0
        eneigval = -ionlim/atomenergy
    elseif nn==5&&ll==1&&jj==1/2
        eneigval = (D1trans-ionlim)/atomenergy
    elseif nn==5&&ll==1&&jj==3/2
        eneigval = (D2trans-ionlim)/atomenergy
    elseif nn==4&&ll==2&&jj==5/2
        eneigval = (n4l2j5_2en-ionlim)/atomenergy
    elseif nn==4&&ll==2&&jj==3/2
        eneigval = (n4l2j3_2en-ionlim)/atomenergy
    elseif nn==6&&ll==0
        eneigval = (n6l0j1_2en-ionlim)/atomenergy
    elseif nn==6&&ll==1&&jj==1/2
        eneigval = (n6l1j1_2en-ionlim)/atomenergy
    elseif nn==6&&ll==1&&jj==3/2
        eneigval = (n6l1j3_2en-ionlim)/atomenergy
    else
        eneigval = -Ry/(2*(nn-delta_nlj)^2)
    end

    #Modified Coulomb potential model parameters from Marinescu et al. (1994)
    alpha_c = 9.0760
    if ll==0
        a_1 = 3.69628474
        a_2 = 1.64915255
        a_3 = -9.86069196
        a_4 = 0.19579987
        r_c = 1.66242117
    elseif ll==1
        a_1 = 4.44088978
        a_2 = 1.92828831
        a_3 = -16.79597770
        a_4 = -0.81633314
        r_c = 1.50195124
    elseif 11==2
        a_1 = 3.78717363
        a_2 = 1.57027864
        a_3 = -11.65588970
        a_4 = 0.52942835
        r_c = 4.86851938
    else
        a_1 = 2.39848933
        a_2 = 1.76810544
        a_3 = -12.07106780
        a_4 = 0.77256589
        r_c = 4.79831327
    end

    params = [alpha_c,a_1,a_2,a_3,a_4,r_c]

    return ZZ, spin, eneigval, params, delta_nlj
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using ExploreFunction, Test, Random

# Only basic tests for now

@testset "generatePerturbation" begin
    @testset "Power" begin
        rng = MersenneTwister(1234)
        perturb = generatePerturbation(rng,1000000,10,1)
        @test 10 ≈ sqrt(sum(abs2,perturb)/1000000) atol=0.1
    end
    @testset "length" begin
        perturb = generatePerturbation(100,1,1)
        @test length(perturb) == 100
    end
end

# Should add more tests here at some point
@testset "MCMCTrajectory" begin
    @testset "Dimensions" begin
        x₀ = [1, 2]
        fₓ(x) = x
        traject, ~ = MCMCTrajectory(fₓ,x₀,1,100,1,10)
        @test size(traject) == (4,100)
    end
end

@testset "minimumdistance" begin
    A = zeros(1,3)
    B = zeros(1,4)
    A[1,:] = [1; 3; 4]
    B[1,:] = [5; 6; 7; 8]
    @test minimumdistance(A,B) == 1
    @test minimumdistance(B,A) == 1
    @test minimumdistance(A,A) == 0
    B[1,:] = [9; 6; 7; 8]
    @test minimumdistance(A,B) == 2
end

@testset "detectTransitions" begin
    disc_test = zeros(100)
    threshold = 0.5
    pathlength = 21
    disc_test[10] = 1
    disc_test[30] = 1

    unstable, stable = detectTransitions(disc_test,threshold,pathlength)

    @test unstable == [20, 40]
    stable_true = zeros(Int,98)

    for i in 1:98
        if i<10
            stable_true[i] = i + 10
        elseif i<29
            stable_true[i] = i + 11
        else
            stable_true[i] = i + 12
        end
    end
    @test stable == stable_true
end

@testset "GenerateConstraints" begin
    directions = randn(3,5)
    gradients = randn(3,5)

    A1_true = directions[:,1:2]
    A2_true = -directions[:,4:5]
    b1_true = zeros(2)
    b2_true = zeros(2)
    b1_true[1] = A1_true[:,1]' * gradients[:,1]
    b1_true[2] = A1_true[:,2]' * gradients[:,2]
    b2_true[1] = A2_true[:,1]' * gradients[:,4]
    b2_true[2] = A2_true[:,2]' * gradients[:,5]

    A1, A2, b1, b2 = ExploreFunction.GenerateConstraints(gradients, directions)
    @test A1 == A1_true'
    @test A2 == A2_true'
    @test b1 == b1_true
    @test b2 == b2_true
end

@testset "SolveMinProg" begin
    A1 = [1/sqrt(2)  1/sqrt(2);
          1/sqrt(2) -1/sqrt(2)]
    A2 = -A1;
    b1 = ones(2)
    b2 = ones(2)
    true_min_dist = 2sqrt(2)
    @test true_min_dist ≈ ExploreFunction.SolveMinProb(A1,A2,b1,b2, 0, false) atol=0.01
end

@testset "MinimumHoleSize" begin
    gradients = ones(1,5)
    gradients[1,3] = 0
    gradients[1,4:5] .= -1
    directions = ones(1,5)

    true_size = 2
    @test true_size ≈ MinimumHoleSize(gradients,directions, 0) atol=0.01
end

@testset "FormNetwork" begin
    dataset = 0.01randn(2,90)
    dataset[1,11:20] += range(0,1,length=10)
    dataset[1,21:end] .+= 1
    dataset[2,31:40] += range(0,1,length=10)
    dataset[2,41:end] .+= 1
    dataset[1,51:60] -= range(0,1,length=10)
    dataset[1,61:end] .-= 1
    dataset[2,71:80] -= range(0,1,length=10)
    dataset[2,81:end] .-= 1

    transitions = zeros(40)
    transitions[1:10] = 11:20
    transitions[11:20] = 31:40
    transitions[21:30] = 51:60
    transitions[31:40] = 71:80

    threshold = 0.2

    edges, mean = FormNetwork(dataset,transitions,threshold)
    # Sort the edges lexicographical order
    sort!(edges,by= x -> 10*minimum(x) + maximum(x))
    edges_comp = zeros(Int,2,length(edges))
    for i in 1:length(edges)
        if edges[i][1] < edges[i][2]
            edges_comp[1,i] = edges[i][1]
            edges_comp[2,i] = edges[i][2]
        else
            edges_comp[2,i] = edges[i][1]
            edges_comp[1,i] = edges[i][2]
        end
    end

    true_means = zeros(size(mean))
    true_means[:,1] = [0,0]
    true_means[:,2] = [1,0]
    true_means[:,3] = [1,1]
    true_means[:,4] = [0,1]

    true_edges = [  1 1 2 3;
                    2 4 3 4]

    @test mean ≈ true_means atol = 0.1
    @test true_edges == edges_comp
end


text:
using AlphaZero
using ProgressMeter
using Statistics: mean

# include("../games/tictactoe/main.jl")
# using .Tictactoe: Game

include("../games/connect-four/main.jl")
using .ConnectFour: Game

baseline = MinMax.Player{Game}(depth=5, amplify_rewards=true, τ=0.2)
mcts = MCTS.Env{Game}(MCTS.RolloutOracle{Game}())
mcts = MctsPlayer(mcts, niters=1000, τ=ConstSchedule(0.5))
player = TwoPlayers(mcts, baseline)

num_games = 200
bar = Progress(num_games)
rewards = map(1:num_games) do _
  trace = play_game(player)
  next!(bar)
  return total_reward(trace)
end

println("Average reward: $(mean(rewards))")


text:
using Base: @_inline_meta

function findstate(V::Vector{Ket{T,D}},s::Union{Tuple,Vector}) where {T,D}
    if length(s) != D || any(s .>= dims(V[1]))
        throw(ArgumentError("state $s incompatible with system dimensions $(dims(V[1]))"))
    end
    ti = tensored_sub2ind(dims(V[1]),s)
    j, jval = 1, abs2(V[1][ti]) # initial values
    for i = 2:length(V)
        ival = abs2(V[i][ti])
        if ival > jval; j, jval = i, ival; end
    end
    return j
end

# the below function is useful to find states at resonance
function findstate(V::Vector{Ket{T,D}},s1::Union{Tuple,Vector},s2::Union{Tuple,Vector}) where {T,D}
    if length(s1) != D || length(s2) != D || any(s1 .>= dims(V[1])) || any(s2 .>= dims(V[1]))
        throw(ArgumentError("state $s incompatible with system dimensions $(dims(V[1]))"))
    end
    ti1 = tensored_sub2ind(dims(V[1]),s1)
    ti2 = tensored_sub2ind(dims(V[1]),s2)
    j1, j1val = 1, zero(abs2(V[1][ti1])) # initial values
    j2, j2val = 1, j1val
    for i = 1:length(V)
        i1val = abs2(V[i][ti1]); i2val = abs2(V[i][ti2])
        if i1val > j1val # found bigger value for state 1
            if abs2(V[j1][ti2]) > j2val # check if current index is better for state 2
                j2, j2val = j1, abs2(V[j1][ti2])
            end
            j1, j1val = i, i1val # now do the swap
            continue
        end
        if i2val > j2val # found bigger value for state 2
            if abs2(V[j2][ti1]) > j1val # check if current index is better for state 1
                j1, j1val = j2, abs2(V[j2][ti1])
            end
            j2, j2val = i, i2val # now do the swap
            continue
        end
    end
    return j1, j2
end

# index in a tensored system (WARNING: subscripts must be 0-indexed, but returns 1-indexed index)
tensored_sub2ind(dims::Dims{N}, inds::Vector{<:Integer}) where {N} =
    (checkbounds(inds,N); @_inline_meta; _tensored_sub2ind(dims, inds, 1, 1))
tensored_sub2ind(dims::Dims{N}, inds::Dims{N}) where {N} =
    (@_inline_meta; _tensored_sub2ind(dims, inds, 1, 1))
_tensored_sub2ind(dims::Dims{N}, inds, D, ind) where {N} =
    (@_inline_meta; _tensored_sub2ind(front(dims), inds, D*dims[N], ind+inds[N]*D))
_tensored_sub2ind(dims::Tuple{Int}, inds, D, ind) = (ind + inds[1]*D)

# give the tensored subscripts corresponding to the index (0-indexed)
tensored_ind2sub(dims::Dims, ind::Integer) =
    (@_inline_meta; next=ind÷dims[end]; (tensored_ind2sub(front(dims), next)..., ind-dims[end]*next))
tensored_ind2sub(dims::Dims{1}, ind::Integer) =
    (@_inline_meta; (ind,))


text:
# This file is a part of the ExamplePackage.

module ExamplePackage
include("sort_numbers.jl")
end


text:
#############################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# JuMP implementation of the cannery problem
# Dantzig, Linear Programming and Extensions,
# Princeton University Press, Princeton, NJ, 1963.
#
# Author: Louis Luangkesorn
# Date January 30, 2015
#############################################################################

using JuMP

function PrintSolution(status, plants, markets, ship)
    println("RESULTS:")
    if status == :Optimal
      for i = 1:length(plants)
        for j in 1:length(markets)
          println("  $(plants[i]) $(markets[j]) = $(getvalue(ship[i,j]))")
        end
      end
    else
        println("  No solution")
    end
    println("")
end

function solveCannery(plants, markets, capacity, demand, distance, freight)
  numplants = length(plants)
  nummarkets = length(markets)
  cannery = Model()

  @variable(cannery, ship[1:numplants, 1:nummarkets] >= 0)

  # Ship no more than plant capacity
  @constraint(cannery, xyconstr[i=1:numplants],
                   sum{ship[i,j], j=1:nummarkets}<=capacity[i])

  # Ship at least market demand
  @constraint(cannery, xyconstr[j=1:nummarkets],
               sum{ship[i,j], i=1:numplants}>=demand[j])

  # Minimize transporatation cost
  @objective(cannery, Min,
              sum{distance[i,j]* freight * ship[i,j], i=1:numplants, j=1:nummarkets})

  solution = solve(cannery)
  if solution == :Optimal
    result = ship
    PrintSolution(solution, plants, markets, ship)
  else
    result = solution
    print("No solution")
  end
  result
end

# PARAMETERS

plants = ["Seattle", "San-Diego"]
markets = ["New-York", "Chicago", "Topeka"]

# capacity and demand in cases
capacitycases = [350, 600]
demandcases = [300, 300, 300]

# distance in thousand miles
distanceKmiles = [2.5 1.7 1.8;
            2.5 1.8 1.4]

# cost per case per thousand miles
freightcost = 90

solveCannery(plants, markets, capacitycases, demandcases, distanceKmiles, freightcost)


text:
function direct_neighbors(graph, node)
    neighbors = graph[node]
    if node ∉ neighbors
        pushfirst!(neighbors, node)
    end
    return neighbors
end

function direct_neighbors(graph::Matrix{Bool}, node)
    nodes = collect(1:size(graph)[1])
    neighbors = nodes[graph[node,:]]
    if node ∉ neighbors
        pushfirst!(neighbors, node)
    end
    return neighbors
end

function reachable(graph, node)
    return dfs([], graph, node)
end

function dfs(visited, graph, node)
    if node ∉ visited
        push!(visited, node)
        for neighbor in direct_neighbors(graph, node)
            dfs(visited, graph, neighbor)
        end
    end
    return visited
end

function conn_comp(graph)
    comps = Array{Int}[]
    for node in collect(1:size(graph)[1])
        if length(comps)<1
            push!(comps, reachable(graph, node))
        elseif  sum(map(x->node in x, comps)) < 1  #if node has not be found/added to a connection yet
            push!(comps, reachable(graph, node)) #then add it
        end
    end
    return comps
end


text:
# cartesian and homogeneous orientation determinants
using DoubleFloats
#  iszero, isnonzero, isone                 #  value == 0, value != 0, value == 1
#  ispositive, isnegative,                  #  value >  0, value <  0

## I don't want Float64 to widen to BigFloat, but to Double64
widemult(a::Float64, b::Real) = Double64(a)*Double64(b)
widemult(a::Union{Float32,Int}, b::Float64) = Double64(a)*Double64(b)
widemult(a, b) = widemul(a,b)

##

" det2
	determinants using wider multiplication
"
det2(p,q, a,b) = @inline @inbounds muladd(p[a],q[b], -widemult(p[b],q[a])) # 2 mult, add, neg
det3(p,q,r, a,b,c) = muladd(p[a], det2(q,r,b,c), # 9 mult, 5 add, 3 neg,
					 muladd(p[b], det2(r,q,a,c),
						    p[c]*det2(q,r,a,b)))
det4ck(p,q,r,s, a,b,c,d) = muladd(p[a], det3(q,r,s, b,c,d), # 40 mult, 23 add, 12 neg
						 muladd(p[b], det3(q,r,s, a,d,c),
						 muladd(p[c], det3(q,r,s, a,b,d),
								p[d]*det3(q,r,s, a,c,b))))

# 2x2 minors is better
det4(p,q,r,s, a,b,c,d) = Double64(det2(p,q, a,b))*Double64(det2(r,s, c,d)) - # 30 mult, 17 add, 12 neg
						 Double64(det2(p,q, a,c))*Double64(det2(r,s, b,d)) +
						 Double64(det2(p,q, a,d))*Double64(det2(r,s, b,c)) +
						 Double64(det2(p,q, b,c))*Double64(det2(r,s, a,d)) -
						 Double64(det2(p,q, b,d))*Double64(det2(r,s, a,c)) +
						 Double64(det2(p,q, c,d))*Double64(det2(r,s, a,b))

Plucker3L(p,q) = (q.-p, det2(p,q, 1,2), det2(p,q, 1,3), det2(p,q, 1,4), det2(p,q, 2,3), det2(p,q, 2,4), det(p,q, 3,4))

det4dbl(p,q,r,s, a,b,c,d) = muladd(det2(p,q, a,b), det2(r,s, c,d), # 30 mult, 17 add, 12 neg
						 muladd(det2(p,q, c,a), det2(r,s, b,d),
						 muladd(det2(p,q, a,d), det2(r,s, b,c),
						 muladd(det2(p,q, b,c), det2(r,s, a,d),
						 muladd(det2(p,q, d,b), det2(r,s, a,c), det2(p,q, c,d))*det2(r,s, a,b)) ) ) )

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
## Generic and Fallbacks
"""
    discount(yc, to)
    discount(yc, from,to)

The discount factor for the yield curve `yc` for times `from` through `to`.
"""
discount(yc::YieldCurve, time) = exp(-yc.zero(time) * time)
discount(yc::AbstractYield, from, to) = discount(yc, to) / discount(yc, from)

"""
    forward(yc, from, to, CompoundingFrequency=Periodic(1))

The forward `Rate` implied by the yield curve `yc` between times `from` and `to`.
"""
function forward(yc::AbstractYield, from, to)
    return forward(yc, from, to, Periodic(1))
end

function forward(yc::AbstractYield, from, to, cf::CompoundingFrequency)
    r = Periodic((accumulation(yc, to) / accumulation(yc, from))^(1 / (to - from)) - 1, 1)
    return convert(cf, r)
end

function forward(yc::AbstractYield, from)
    to = from + 1
    return forward(yc, from, to)
end

"""
    accumulation(yc, from, to)

The accumulation factor for the yield curve `yc` for times `from` through `to`.
"""
function accumulation(yc::AbstractYield, time)
    return 1 ./ discount(yc, time)
end

function accumulation(yc::AbstractYield, from, to)
    return 1 ./ discount(yc, from, to)
end


text:
# https://discourse.julialang.org/t/generate-random-value-from-a-given-function-out-of-box/5793/3
using Waterfall

include(joinpath(WATERFALL_DIR,"src","includes.jl"))

DATA_DIR = joinpath(WATERFALL_DIR,"data","pvrd")
df = CSV.read(joinpath(DATA_DIR,"pvrd2-investment-metrics.csv"), DataFrames.DataFrame)
df[!,:Label] .= string.(titlecase.(df[:,:Process])," (",df[:,:Step],")")
df[!,:Units] .= "Efficiency (%)"

# Define global arguments.
distribution=:normal
fuzziness=(0.01,0.3)
kwargs = (label=:Label, distribution=distribution, fuzziness=fuzziness)

# Define highlights to mark.
prob = [0.75,0.25]
highlights=["mean"; [("quantile",p) for p in prob]]

for samples in [1,5,10,50]

    pdata = Plot(df; samples=samples, ylabel=:Units, kwargs...)

    # Iterate over plot type.
    for T in [Horizontal,Parallel,Vertical,Violin]
        p = Plot{T}(pdata)

        # Show different combinations of labeling.
        # for hh in [[], [1], 1:length(highlights)]
        for hh in [[]]
            (samples==1 && !isempty(hh)) && continue

            f = filename(p, highlights[hh]; distribution=distribution)
            Printf.@printf("\nPlotting and saving figure to %s", f)

            Luxor.@png begin
                Luxor.fontsize(14)
                Luxor.setmatrix([1 0 0 1 LEFT_BORDER TOP_BORDER])

                draw(p; distribution=distribution, samples=samples)
                # _draw_highlight(pdata, highlights[hh])

            end WIDTH+LEFT_BORDER+RIGHT_BORDER HEIGHT+TOP_BORDER+BOTTOM_BORDER f
        end
    end
end


text:
function weighted_moment_distance(wm::BlockBootstrapWeightMatrix, sims)
    observed_moments = select_moments(wm.obs)

    n_replications = size(sims, 2)
    G = zeros(n_replications, 5)

    sim_moments_matrix = get_summary_stats(sims, wm.obs)
    G = sim_moments_matrix - repeat(observed_moments', n_replications, 1)
    G = mean(G, dims = 1)
    return (G * wm.weights * G')[]
end


text:
import DataFrames: DataFrame, DataFrameRow, nrow, ncol, mapcols, missings
import Random: shuffle
import Statistics: mean, median, std
import CUDA: CuVector, CuMatrix
import Tables

export num_examples, num_features,
    example, feature_values,
    iscomplete, isweighted, isfpdata, isbinarydata,
    num_chunks, chunks, eltype,
    weigh_samples, split_sample_weights, get_weights,
    shuffle_examples, batch, batch_size, isbatched,
    threshold, soften, marginal_prob, bagging_dataset,
    random_sample,
    to_gpu, to_cpu, isgpu, same_device,
    ll_per_example, bits_per_pixel,
    make_missing_mcar, impute

##############################
# Generic DataFrame Extensions
##############################

import Base: eltype # extend

"Find a type that can capture all column values"
eltype(df::DataFrame) = reduce(typejoin, eltype.(eachcol(df)))

import Base: _msk_end

# lift BitVector function to DataFrame of BitVectors
@inline _msk_end(df::DataFrame) = _msk_end(nrow(df))

##############################
# GPU-CPU memory helper function
##############################

"Move data to the GPU"
to_gpu(d::Union{CuBitVector,CuArray}) = d
to_gpu(m::Array) = CuArray(m)
to_gpu(v::BitVector) =
    AbstractBitVector(to_gpu(chunks(v)), length(v))
to_gpu(df::DataFrame) = isgpu(df) ? df : mapcols(to_gpu, df)
to_gpu(df::Vector{DataFrame}) = map(df) do d
    to_gpu(d)
end

"Move data to the CPU"
to_cpu(m::Union{BitVector,Array}) = m
to_cpu(m::CuArray) = Array(m)
to_cpu(cv::CuBitVector) =
    AbstractBitVector(to_cpu(cv.chunks), length(cv))
to_cpu(df::DataFrame) = mapcols(to_cpu, df)
to_cpu(df::Vector{DataFrame}) = map(df) do d
    to_cpu(d)
end

"Check whether data resides on the GPU"
isgpu(::Union{Array, BitArray}) = false
isgpu(::Union{CuArray, CuBitVector}) = true
isgpu(df::DataFrame) = all(isgpu, eachcol(df))
isgpu(df::Vector{DataFrame}) = all(isgpu, df)

"Ensure that `x` resides on the same device as `data`"
same_device(x, data) = isgpu(data) ? to_gpu(x) : to_cpu(x)

##############################
# Dataset Utilities
##############################

# Basic data-related utilities.
# These methods assume the following possible representations.
# 1/ Complete data
#  - `Bool` data:
#      * `DataFrame` of `(Cu)BitVector`
#  - `Float` data:
#      * `DataFrame` of `(Cu)Vector{Float}`
# 2/ Incomplete data:
#  - `Bool` data:
#      * `DataFrame` of `Vector{Union{Bool,Missing}}`
#      * `DataFrame` of `CuVector{Union{Missing, Bool}}`
#  - `Float` data:
#      * `DataFrame` of `Vector{Union{Float,Missing}}`
#      * `DataFrame` of `CuVector{Union{Missing, Float}}`

"""
    num_examples(df::DataFrame)

Number of examples in data
"""
num_examples(df::DataFrame) = nrow(df)
num_examples(df::Vector{DataFrame}) = mapreduce(num_examples, +, df)

"""
    num_features(df::DataFrame)

Number of features in the data
"""
num_features(df::DataFrame) = isweighted(df) ? ncol(df) - 1 : ncol(df)
num_features(df::Vector{DataFrame}) = num_features(df[1])

"Get the ith example"
example(d::DataFrame,i) = vec(convert(Array, collect(d[i,:])))

"Get the ith feature values"
feature_values(df::DataFrame, i) = df[!,i]

"Is the data complete (no missing values)?"
iscomplete(data::DataFrame) =
    all(iscomplete_col, eachcol_unweighted(data))
iscomplete(data::Vector{DataFrame}) = all(iscomplete, data)

"Is the data column complete (no missing values)?"
iscomplete_col(::AbstractVector{Bool}) = true
iscomplete_col(::AbstractVector{<:Int}) = true
iscomplete_col(::AbstractVector{<:AbstractFloat}) = true
iscomplete_col(x::AbstractVector{Union{Bool,Missing}}) = !any(ismissing, x)
iscomplete_col(x::AbstractVector{Union{<:AbstractFloat,Missing}}) = !any(ismissing, x)

"Is the dataset binary?"
isbinarydata(df::DataFrame) = all(eachcol_unweighted(df)) do col
    nonmissingtype(eltype(col)) <: Union{Bool}
end
isbinarydata(df::Vector{DataFrame}) =
    all(isbinarydata, df)

"Is the dataset consisting of floating point data?"
isfpdata(df::DataFrame) = all(eachcol_unweighted(df)) do col
    nonmissingtype(eltype(col)) <: AbstractFloat
end
isfpdata(df::Vector{DataFrame}) =
    all(isfpdata, df)

"For binary complete data, how many `UInt64` bit strings are needed to store one feature?"
num_chunks(d::DataFrame) = num_chunks(feature_values(d,1))
num_chunks(v::AbstractBitVector) = length(v.chunks)

"Retrieve chunks of bit vector"
chunks(v::AbstractBitVector) = v.chunks
chunks(df::DataFrame, i) = chunks(feature_values(df, i))

# WEIGHTED DATASETS

"Iterate over columns, excluding the sample weight column"
eachcol_unweighted(data::DataFrame) =
    isweighted(data) ? eachcol(data)[1:end-1] : eachcol(data)

"Create a weighted copy of the data set"
weigh_samples(df::DataFrame, weights) = begin
    @assert length(weights) == num_examples(df)
    df2 = copy(df)
    df2.weight = weights
    df2
end

"Split a weighted dataset into unweighted dataset
 and its corresponding weights."
split_sample_weights(df::DataFrame) = begin
    @assert isweighted(df) "`df` is not weighted."
    df[!, 1:end-1], df[!, end]
end

split_sample_weights(batches::Vector{DataFrame}) = begin
    @assert isweighted(batches) "`batches` is not weighted."
    samples = map(b -> b[!, 1:end-1], batches)
    weights = get_weights(batches)
    (samples, weights)
end

"Get the weights from a weighted dataset."
get_weights(df::DataFrame) = begin
    @assert isweighted(df) "`df` is not weighted."
    df[!, end]
end
get_weights(df::Vector{DataFrame}) = begin
    map(get_weights, df)
end

"Is the dataset weighted?"
isweighted(df::DataFrame) =
    names(df)[end] == "weight"
isweighted(df::Vector{DataFrame}) =
    all(d -> isweighted(d), df)

# DATA TRANSFORMATIONS

"""
    shuffle_examples(df::DataFrame)

Shuffle the examples in the data
"""
shuffle_examples(data) = data[shuffle(axes(data, 1)), :]

"Create mini-batches"
function batch(data::DataFrame, batchsize=1024; shuffle::Bool = true)
    from_gpu = isgpu(data)
    # move data to CPU to do minibatching
    from_gpu && (data = to_cpu(data))

    shuffle && (data = shuffle_examples(data))

    n_examples = num_examples(data)
    data = map(1:batchsize:n_examples) do start_index
        stop_index = min(start_index + batchsize - 1, n_examples)
        data[start_index:stop_index, :]
    end

    # move data back to GPU if that's where it came from
    from_gpu ? to_gpu(data) : data
end

"Is the dataset batched?"
isbatched(::DataFrame) = false
isbatched(::Vector{DataFrame}) = true

"Batch size of the dataset"
batch_size(data::Vector{DataFrame}) = num_examples(data[1])

"Randomly draw samples from the dataset with replacement"
function random_sample(data::DataFrame, num_samples=1, weighted=isweighted(data))
    # TODO add keyword to fix the random seed
    example_idxs = if weighted
        data, weights = split_sample_weights(data)
        cum_weights = cumsum(weights)
        total_weight = cum_weights[end]
        map(1:num_samples) do idx
            randval = rand() * total_weight
            i, j, k = 1, num_examples(data), 0
            #  do binary search for sample
            while i < j - 1
                k = (i + j) ÷ 2
                if randval < cum_weights[k]
                    j = k
                else
                    i = k
                end
            end
            i == j ? i : (randval >= cum_weights[i] ? j : i)
        end
    else
        map(1:num_samples) do idx
            rand(1:num_examples(data))
        end
    end
    data[example_idxs, :]
end

"Returns an array of DataFrames where each DataFrame is randomly sampled from the original dataset `data`"
function bagging_dataset(data::DataFrame; num_bags::Integer = 1, frac_examples::AbstractFloat = 1.0,
                         batch_size::Integer = 0)
    map(1:num_bags) do dataset_idx
        sample = random_sample(data, floor(UInt32, num_examples(data) * frac_examples))
        batch_size == 0 ? sample : batch(sample, batch_size)
    end
end

"Threshold a numeric dataset making it binary"
threshold(train, valid, test) = threshold(train, valid, test, 0.05) # default threshold offset (used for MNIST)

function threshold(train::DataFrame, valid, test, offset)
    @assert isfpdata(train) "DataFrame to be thresholded contains non-numeric columns."
    train = Tables.matrix(train)
    valid = issomething(valid) ? Tables.matrix(valid) : nothing
    test = issomething(test) ? Tables.matrix(test) : nothing
    means = mean(train, dims=1)
    stds = std(train, dims=1)
    threshold_value = means .+ offset .* stds
    train = DataFrame(BitArray(train .> threshold_value), :auto)
    valid = issomething(valid) ? DataFrame(BitArray(valid .> threshold_value), :auto) : nothing
    test = issomething(test) ? DataFrame(BitArray(test .> threshold_value), :auto) : nothing
    return train, valid, test
end

"Turn binary data into floating point data close to 0 and 1."
soften(data::DataFrame, softness=0.05; scale_by_marginal=true, precision=Float32) = begin
    @assert !isfpdata(data) "'soften' does not support floating point data."
    n_col = ncol(data)
    data_weighted = isweighted(data)

    marginals = marginal_prob(data; precision = precision)
    vs = AbstractVector[]
    col_idx = 1
    for v in eachcol(data)
        if data_weighted && col_idx == n_col
            fv = v
        elseif scale_by_marginal
            fv = v .* precision(1.0 - softness) .+ (precision(softness) * marginals[col_idx])
        else
            fv = v .* precision(1.0 - 2.0 * softness) .+ precision(softness)
        end
        push!(vs, fv === v ? copy(fv) : fv)

        col_idx += 1
    end
    DataFrame(vs, names(data), copycols = false)
end

"Compute the marginal prob of each feature in a binary dataset."
marginal_prob(data; precision=Float32) = begin
    @assert !isfpdata(data) "'marginal_prob' does not support floating point data."
    n_examples = num_examples(data)
    map(1:num_features(data)) do idx
        precision(sum(coalesce.(feature_values(data, idx), 0.5))) / n_examples
    end
end

# LIKELIHOOD HELPERS

"Normalize the given log-likelihood by the number of examples in `data`"
ll_per_example(ll, data) = ll / num_examples(data)

"Normalize the given log-likelihood as bits per pixel in `data`"
bits_per_pixel(ll, data) = -(ll_per_example(ll, data)  / num_features(data)) / log(2)

"Computer the per-example log-likelihood of a fully factorized ML model on Bool data"
function fully_factorized_log_likelihood(data; pseudocount=0)
    @assert isbinarydata(data) && iscomplete(data) "This method requires complete binary data"
    m = Tables.matrix(data)
    counts = sum(m, dims=1)
    smoothed_counts = counts .+ pseudocount/2.0
    log_estimates = log.(smoothed_counts ./ (num_examples(data) + pseudocount))
    ll = sum(log_estimates .* counts)
    ll += sum(log.(1 .- exp.(log_estimates)) .* (num_examples(data) .- counts))
    ll_per_example(ll, data)
end

##############################
# Imputations & Missing value generation
##############################

"""
    make_missing_mcar(d::DataFrame; keep_prob::Float64=0.8)

Returns a copy of dataframe with making some features missing as MCAR, with
`keep_prob` as probability of keeping each feature.
"""
function make_missing_mcar(d::DataFrame; keep_prob::Float64=0.8)
    m = missings(eltype(d), num_examples(d), num_features(d))
    flag = rand(num_examples(d), num_features(d)) .<= keep_prob
    m[flag] .= Matrix(d)[flag]
    DataFrame(m, names(d))
end;

"""
Return a copy of Imputed values of X  (potentially statistics from another DataFrame)

For example, to impute using same DataFrame:

    impute(X; method=:median)

If you want to use another DataFrame to provide imputation statistics:

    impute(test_x, train_x; method=:mean)

Supported methods are `:median`, `:mean`, `:one`, `:zero`
"""
function impute(X::DataFrame; method=:median)
    impute(X, X; method=method)
end
function impute(X::DataFrame, train::DataFrame; method=:median)
    type = typeintersect(eltype(X), eltype(train))
    @assert type !== Union{}

    if typeintersect(type, Bool) == Bool
        type = Bool
    elseif typeintersect(type, AbstractFloat) <: AbstractFloat
        type = typeintersect(type, AbstractFloat)
    end
    @assert !(type isa Union) "Type error $(type)"

    if method == :median
        impute_function = median
    elseif method == :mean
        impute_function = mean
    elseif method == :one
        impute_function = (x -> one(type))
    elseif method == :zero
        impute_function = (x -> zero(type))
    else
        throw("Unsupported imputation type $(method)")
    end

    X_impute = deepcopy(X)
    for feature = 1:size(X)[2]
        mask_train = ismissing.(train[:, feature])
        mask_x     = ismissing.(X[:, feature])

        cur_impute = impute_function(train[:, feature][.!(mask_train)] )

        if type == Bool
            X_impute[mask_x, feature] .= Bool(cur_impute .>= 0.5)
        else
            X_impute[mask_x, feature] .= type(cur_impute)
        end
    end

    # For Bool return BitArray instead
    if type == Bool
        return DataFrame(BitArray(Tables.matrix(X_impute)), names(X_impute))
    else
        return X_impute::DataFrame
    end
end


text:
@testset "Day2.jl" begin
    # Write your tests here.
    @test solve(Day2_1, "$(@__DIR__)/assets/day2_test1.txt") == 2

    @test solve(Day2_2, "$(@__DIR__)/assets/day2_test1.txt") == 1
end


text:
export interaction_matrix, PowerLaw

abstract type Interaction end

"""
    interaction_matrix(interaction, geometry, positions)
    interaction_matrix(interaction, distances)

Compute the interaction matrix from given distances or geometry and positions.
The latter is only supported for anistrope interactions as these do only depend on the distance.
"""
function interaction_matrix end

abstract type AnisotropeInteraction <: Interaction end

interaction_matrix(interaction::AnisotropeInteraction, geometry::Geometry, positions) = interaction_matrix(interaction, distance_matrix(geometry, positions))

interaction_matrix(interaction::AnisotropeInteraction, distances) = _interaction_strength.(Ref(interaction), distances)

"""
    _interaction_strength(<:AnistropeInteraction, distance)
Compute interaction strength given the distance.
"""
function _interaction_strength end

"""
    PowerLaw{T}

Interaction strength scales proportional to distance to the power -α.

# Fields
- `α::T`: Exponent
"""
struct PowerLaw <: AnisotropeInteraction
    α::Float64
    PowerLaw(α) = new(convert(Float64, α))
end

_interaction_strength(int::PowerLaw, distance::Float64) = float(distance > 0 ? distance^-int.α : 0)


text:
module BasiTest
using SimplePNGs
using Test

include("common.jl")
using .TestCommon: load_json

pl(name) = SimplePNGs.load(joinpath("PngSuite", name*".png"))

@testset "Basic format test files (Adam-7 interlaced)" begin
    @testset "black & white" begin
        img1 = load_json("basn0g01")
        img2 = pl("basi0g01")
        @test img1 == img2
    end

    @testset "2 bit (4 level) grayscale" begin
        img1 = load_json("basn0g02")
        img2 = pl("basi0g02")
        @test img1 == img2
    end

    @testset "4 bit (16 level) grayscale" begin
        img1 = load_json("basn0g04")
        img2 = pl("basi0g04")
        @test img1 == img2
    end

    @testset "8 bit (256 level) grayscale" begin
        img1 = load_json("basn0g08")
        img2 = pl("basi0g08")
        @test img1 == img2
    end

    @testset "16 bit (64k level) grayscale" begin
        img1 = load_json("basn0g16")
        img2 = pl("basi0g16")
        @test img1 == img2
    end

    @testset "3x8 bits rgb color" begin
        img1 = load_json("basn2c08")
        img2 = pl("basi2c08")
        @test img1 == img2
    end

    @testset "3x16 bits rgb color" begin
        img1 = load_json("basn2c16")
        img2 = pl("basi2c16")
        @test img1 == img2
    end

    @testset "1 bit (2 color) paletted" begin
        img1 = load_json("basn3p01")
        img2 = pl("basi3p01")
        @test img1 == img2
    end

    @testset "2 bit (4 color) paletted" begin
        img1 = load_json("basn3p02")
        img2 = pl("basi3p02")
        @test img1 == img2
    end

    @testset "4 bit (16 color) paletted" begin
        img1 = load_json("basn3p04")
        img2 = pl("basi3p04")
        @test img1 == img2
    end

    @testset "8 bit (256 color) paletted" begin
        img1 = load_json("basn3p08")
        img2 = pl("basi3p08")
        @test img1 == img2
    end

    @testset "8 bit grayscale + 8 bit alpha-channel" begin
        img1 = load_json("basn4a08")
        img2 = pl("basi4a08")
        @test img1 == img2
    end

    @testset "16 bit grayscale + 16 bit alpha-channel" begin
        img1 = load_json("basn4a16")
        img2 = pl("basi4a16")
        @test img1 == img2
    end

    @testset "3x8 bits rgb color + 8 bit alpha-channel" begin
        img1 = load_json("basn6a08")
        img2 = pl("basi6a08")
        @test img1 == img2
    end

    @testset "3x16 bits rgb color + 16 bit alpha-channel" begin
        img1 = load_json("basn6a16")
        img2 = pl("basi6a16")
        @test img1 == img2
    end
end

end # module


text:
# For debugging
function markit(pcpt, cpt)
    circle(pcpt, r = 10mm)
    line(pcpt, cpt + (-0.2EM, 0), :stroke)
    circle(cpt + (-0.2EM, 0), r = 20mm)
end

# Plots default
default(titlefont = (20, "times"), legendfontsize = 26,
        guidefont = (30, :darkgreen), tickfont = (24, :orange),
        framestyle = :origin, minorgrid = true,
        legend = :topleft, linewidth = 4,
        bottom_margin = 18px, left_margin = 18px)


text:
__precompile__(true)

"""
A highly optimized Business Days calculator written in Julia language. Also known as Working Days calculator.

Website: https://github.com/felipenoris/BusinessDays.jl
"""
module BusinessDays

using Base.Dates

# types.jl
export
    HolidayCalendar,
    CompositeHolidayCalendar

# isholiday.jl
export
    isholiday

# bdays.jl : main functions for BusinessDays module
export
    isweekend,
    isbday,
    tobday,
    advancebdays,
    bdays,
    listholidays,
    listbdays

include("dateutils.jl")
include("types.jl")
include("isholiday.jl")
include("bdayscache.jl")
include("bdays.jl")
include("bdaysvecfun.jl")
include("composite.jl")
include("query.jl")

end # module BusinessDays


text:
using LIKWID
using Base.Threads
using Test
using Random

println("Before: ", LIKWID.get_processor_ids())
@test typeof(LIKWID.get_processor_id()) == Int
@test typeof(LIKWID.pinprocess(0)) == Bool
@test typeof(LIKWID.pinthread(0)) == Bool

topo = LIKWID.get_cpu_topology()
ncores = topo.numSockets * topo.numCoresPerSocket
cpus_firstN = 0:nthreads()-1
cpus_firstN_shuffeled = shuffle(cpus_firstN)
cpus_rand = shuffle(0:ncores-1)[1:nthreads()]

LIKWID.pinthreads(cpus_firstN)
println("First N: ", LIKWID.get_processor_ids())
@test LIKWID.get_processor_ids() == cpus_firstN

LIKWID.pinthreads(cpus_firstN_shuffeled)
println("First N (shuffled): ", LIKWID.get_processor_ids())
@test LIKWID.get_processor_ids() == cpus_firstN_shuffeled

LIKWID.pinthreads(cpus_rand)
println("Rand: ", LIKWID.get_processor_ids())
@test LIKWID.get_processor_ids() == cpus_rand


text:
type _Lead_Lag
  C :: Array{Float64,1}
  T :: Array{Float64,1}
  C_0 :: Float64
end

function _Lead_Lag(x:: Array{Float64,1}, y :: Array{Float64,1},  Tidx :: Array{Int64,1} )
  T = Tidx +0.0
  C = fill(NaN,length(Tidx))
  xa = x - mean(x)
  ya = y - mean(y)
  for (ii, tt) in enumerate(Tidx)
    x1 = xa[1+tt:end]
    y2 = ya[1:end-tt]
    C[ii] = mean(x1.*y2)
  end
  C_0 = mean(xa.*ya)
  return _Lead_Lag(C,T,C_0)
end

type _Power_Density
  PD :: Array{Float64,1}
  T :: Array{Float64,1}
  Freq :: Array{Float64,1}
  PD_0 :: Float64
end

function _Power_Density( x:: Array{Float64,1},y :: Array{Float64,1}, Tidx :: Array{Int64,1})
  return _Power_Density( _Lead_Lag(x,y,Tidx))
end

function _Power_Density( LL :: _Lead_Lag)
  T = LL.T
  PD = fill(NaN, length(T))
  C = LL.C
  Freq = (1:length(LL.T))/maximum(LL.T) * 0.5
  for (ii, ff) in enumerate(Freq)
    PD[ii] = sum( cos( T*2*pi*ff ) .* C .* 2) + LL.C_0
  end

  PD_0 = 2*sum(C) + LL.C_0
  return _Power_Density(PD, T,Freq, PD_0)
end

export _Lead_Lag, _Power_Density


text:
using Test
using CUDAnative, CuArrays
using CUDAatomics

if CUDAnative.functional()
    function kern_atomicadd( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicadd(a,b[i], 2)
        return nothing
    end
    d_a = CuArray(zeros(Float32, 2))
    d_b = CuArray(Float32.(collect(1:1024)))
    @cuda threads=32 blocks=32 kern_atomicadd(d_a, d_b)
    @test abs(Array(d_a)[2] - sum(Array(d_b))) < 1.0f-7

    function kern_atomicsub( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicsub(a,b[i])
        return nothing
    end
    d_a = CuArray(zeros(Int32, 1))
    d_b = CuArray(Int32.(collect(1:1024)))
    @cuda threads=32 blocks=32 kern_atomicsub(d_a, d_b)
    @test Array(d_a)[1] == -524800

    function kern_atomicmin( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicmin(a,b[i])
        return nothing
    end
    d_a = CuArray(Array{Int32}([1025]))
    d_b = CuArray(Int32.(collect(1:1024)))
    @cuda threads=32 blocks=32 kern_atomicmin(d_a, d_b)
    @test Array(d_a)[1] == 1

    function kern_atomicmax( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicmax(a,b[i])
        return nothing
    end
    d_a = CuArray(Array{Int32}([1]))
    d_b = CuArray(Int32.(collect(1:1024)))
    @cuda threads=32 blocks=32 kern_atomicmax(d_a, d_b)
    @test Array(d_a)[1] == 1024

    function kern_atomicinc( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicinc(a,b[i])
        return nothing
    end
    d_a = CuArray(Array{UInt32}([UInt32(0)]))
    d_b = CuArray(repeat(Array{UInt32}([UInt32(1024)]), 1024))
    @cuda threads=32 blocks=32 kern_atomicinc(d_a, d_b)
    @test Array(d_a)[1] == 0x00000400

    function kern_atomicdec( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicdec(a,b[i])
        return nothing
    end
    d_a = CuArray(Array{UInt32}([UInt32(1025)]))
    d_b = CuArray(repeat(Array{UInt32}([UInt32(1025)]), 1024))
    @cuda threads=32 blocks=32 kern_atomicdec(d_a, d_b)
    @test Array(d_a)[1] == 0x00000001

    function kern_atomicand( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicand(a,b[i])
        return nothing
    end
    d_a = CuArray(Array{UInt32}([UInt32(1389)]))
    d_b = CuArray(repeat(Array{UInt32}([UInt32(1023)]), 1024))
    @cuda threads=32 blocks=32 kern_atomicand(d_a, d_b)
    @test Array(d_a)[1] == 0x0000016d

    function kern_atomicexch( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomicexch(a,b[i])
        return nothing
    end
    d_a = CuArray(zeros(Float32, 1))
    d_b = CuArray(Float32.(collect(1:1024)))
    @cuda threads=32 blocks=32 kern_atomicexch(d_a, d_b)
    @test findfirst( Array(d_b).==Array(d_a) ) < 1025

    function kern_atomiccas( a, b )
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        atomiccas(a, i, b[i])
        return nothing
    end
    d_a = CuArray(Array{Int32}([17]))
    d_b = CuArray(Array{Int32}(collect(1025:2048)))
    @cuda threads=32 blocks=32 kern_atomiccas(d_a, d_b)
    @test findfirst(Array(d_b).==Array(d_a)) == 17
end


text:
# Adapted directly from Radford Neal's R code

# Alternative interface using Density type
hmc_sampler(x::Array{Float64}, dd::DifferentiableDensity, epsilon::Float64, L::Int64) = hmc_sampler(x, dd.f, dd.gradient, epsilon, L)
# Some default values for epsilon and L
hmc_sampler(x::Array{Float64}, dd::DifferentiableDensity, gx::Float64) = hmc_sampler(x,dd,.01,20)

function hmc_sampler(x::Array{Float64}, dd::DifferentiableDensity, params::Array)
  L = convert(Int64, params[2])
  hmc_sampler(x,dd,params[1],L)  # TODO Use a different data structure for params
end

function hmc_sampler(current_q::Array{Float64}, U::Function, grad_U::Function, epsilon::Float64, L::Int64, )

  q = current_q
  p = randn(length(q))'  # independent standard normal variates (row vector)
  current_p = p

  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U(q) / 2

  # Alternate full steps for position and momentum
  for i in 1:L
    # Make a full step for the position
    q = q + epsilon * p
    # Make a full step for the momentum, except at end of trajectory
    if i!=L
        p = p - epsilon * grad_U(q)
    end
  end

  # Make a half step for momentum at the end.
  p = p - epsilon * grad_U(q) / 2

  # Negate momentum at end of trajectory to make the proposal symmetric
  p = -p

  # Evaluate potential and kinetic energies at start and end of trajectory
  current_U = U(current_q)
  current_K = sum(current_p.^2) / 2
  proposed_U = U(q)
  proposed_K = sum(p.^2) / 2

  # Accept or reject the state at end of trajectory, returning either
  # the position at the end of the trajectory or the initial position
  if rand() < exp(current_U-proposed_U+current_K-proposed_K)
    return q, proposed_U  # accept
  else
    return current_q, current_U  # reject
  end
end

function bounded_hmc_sampler(x::Array{Float64}, dd::DifferentiableDensity, opts::Options)
  bounded_hmc_sampler(x,dd.f,dd.gradient,opts)
end

function bounded_hmc_sampler(current_q::Array{Float64}, U::Function, grad_U::Function, opts::Options)

  @defaults opts stepsize=.01 numsteps=50 bounded=[] lower_bounds=[] upper_bounds=[]

  epsilon = stepsize
  L = numsteps

  q = current_q
  p = randn(length(q))'  # independent standard normal variates (row vector)
  current_p = p

  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U(q) / 2

  # Alternate full steps for position and momentum
  for i in 1:L
    # Make a full step for the position
    q = q + epsilon * p
    # Make a full step for the momentum, except at end of trajectory
    if i!=L
        p = p - epsilon * grad_U(q)
    end

    # If out of bounds, reflect
    for j in bounded
        if q[j] < lower_bounds[j]
            q[j] = lower_bounds[j] + (lower_bounds[j] - q[j])
            p[j] = -p[j]
        end
        if q[j] > upper_bounds[j]
            q[j] = upper_bounds[j] - (q[j] - upper_bounds[j])
            p[j] = -p[j]
        end
    end
  end

  # Make a half step for momentum at the end.
  p = p - epsilon * grad_U(q) / 2

  # Negate momentum at end of trajectory to make the proposal symmetric
  p = -p

  # Evaluate potential and kinetic energies at start and end of trajectory
  current_U = U(current_q)
  current_K = sum(current_p.^2) / 2
  proposed_U = U(q)
  proposed_K = sum(p.^2) / 2

  @check_used opts
  # Accept or reject the state at end of trajectory, returning either
  # the position at the end of the trajectory or the initial position
  if rand() < exp(current_U-proposed_U+current_K-proposed_K)
    return q, proposed_U  # accept
  else
    return current_q, current_U  # reject
  end
end


text:
# Test ArbFloats.jl

using Test
using ArbFloats

function package_directory(pkgName::String)
    pkgdir = Base.find_package(pkgName)
    nothing == pkgdir && throw(ErrorException(string("Package $pkgName not found.")))
    return abspath(joinpath(split(pkgdir, pkgName)[1], pkgName))
end

d = joinpath(package_directory("ArbFloats"),"test")
test_files = [
                "TestAsType.jl",
                "TestAsNumber.jl",
                "TestAsInterval.jl"
              ]

#= prepare to test =#

# include(joinpath(d,"prepareToTest.jl"))

#= begin tests =#

# println("Linting ...")
# using Lint
# @test isempty(lintpkg( "MyPackage", returnMsgs=true))
# println("Done.")

#println("Testing ...")
for f in test_files
#    println(f)
   # include(f)
end
#println("Done testing.")

# end testing ArbFloats


text:
module FF

export ff

import Base: /, *, +, -, promote_rule
import BitIntegers: UBI, @uint256_str, UInt256

const _P_ = uint256"21888242871839275222246405745257275088548364400416034343698204186575808495617"

struct ff <: Integer
    v::UInt256
    ff(v::T) where T <: Integer = begin
        t = UInt256(mod(v, _P_))
        return (t >= 0) ? new(t) : new(t + _P_)
    end
    ff(v::ff) = v
end

@generated promote_rule(::Type{ff}, ::Type{<:UBI}) = ff
(::Type{T})(x::ff) where T <: UBI = T(x.v)

(+)(a::ff, b::ff) = ff(a.v + b.v)
(-)(a::ff, b::ff)  = ff(_P_ + a.v - b.v)
(-)(a::ff)  = ff(_P_ - a.v)
(*)(a::ff, b::ff) = ff(widemul(a.v, b.v) % _P_)
(/)(a::ff, b::ff) = ff(widemul(a.v, invmod(b.v, _P_)) % _P_)

iszero(x::ff) = iszero(x.v)
end #module


text: (b"[encoded data, 8192 characters]", (128, 4))

max_stars_repo_name: joshua-maros/julia-blue-noise

max_issues_repo_name: joshua-maros/julia-blue-noise

max_forks_repo_name: joshua-maros/julia-blue-noise


text:
# ---
# title: 332. Reconstruct Itinerary
# id: problem332
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Graph
# link: <https://leetcode.com/problems/reconstruct-itinerary/description/>
# hidden: true
# ---
#
# Given a list of airline tickets represented by pairs of departure and arrival
# airports `[from, to]`, reconstruct the itinerary in order. All of the tickets
# belong to a man who departs from `JFK`. Thus, the itinerary must begin with
# `JFK`.
#
# **Note:**
#
#   1. If there are multiple valid itineraries, you should return the itinerary that has the smallest lexical order when read as a single string. For example, the itinerary `["JFK", "LGA"]` has a smaller lexical order than `["JFK", "LGB"]`.
#   2. All airports are represented by three capital letters (IATA code).
#   3. You may assume all tickets form at least one valid itinerary.
#   4. One must use all the tickets once and only once.
#
# **Example 1:**
#
#
#
#     Input:[["MUC", "LHR"], ["JFK", "MUC"], ["SFO", "SJC"], ["LHR", "SFO"]]
#     Output:["JFK", "MUC", "LHR", "SFO", "SJC"]
#
#
# **Example 2:**
#
#
#
#     Input:[["JFK","SFO"],["JFK","ATL"],["SFO","ATL"],["ATL","JFK"],["ATL","SFO"]]
#     Output:["JFK","ATL","JFK","SFO","ATL","SFO"]
#     Explanation: Another possible reconstruction is ["JFK","SFO","ATL","JFK","ATL","SFO"].
#                  But it is larger in lexical order.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using MARY_fRG.Basics

@testset "等边六边型" begin
    hexa1 = EqHexagon(Point2D(0., 0.), 1)
    println(hexa1)
    hexa2 = EqHexagon(Point2D(-2., 3.), 1)
    @test abs(area(hexa2) - 3*sqrt(3) / 2) < 1e-6
end


text:
"""
    structure(items)

Recursively format expressions as lists: [operator, operands...]

# Notes:
list trees (aka S-expressions) are convenient to work with because modifying it with a recursive function is easy.
Also, WebAssembly supports text format written in S-expressions so if all required functionality like
overloading, phinodes and such was built into WebAssembly one could essentially do a one liner:
webassemblytext = translate(structure(code_typed(somejuliafunction))

# Details:
- Expressions dont always have the operator in head, sometimes its in args[1] and head is just :call
- pi et al are refs, so if eval(item) is a number just use the number instead of the ref
- Const can hold anything inside it, use the value instead of the container
"""
structure(items::Array) = structure.(items)
structure(item) = item
structure(item::Const) = structure(item.val) #item.val can be anything, including another Expr, so convert that aswell
function structure(item::Expr)
    if item.head == :(call)
        head = item.args[1]
        args = item.args[2:end]
    else
        head = item.head
        args = item.args
    end
    return [head; structure(args)]
end

function structure(item::GlobalRef)
    evaluatedref = Base.eval(Evalscope, item)
    # evaulate refs to constants into constants, also turn π, ℯ et al. to float rather than irrational
    if typeof(evaluatedref) <: Number
        if typeof(evaluatedref) <: AbstractFloat || typeof(evaluatedref) <: Irrational
            return AbstractFloat(evaluatedref)
        else
            return evaluatedref
        end
    else
        return item
    end
end

"""
    restructure(ci::CodeInfo, i::Integer, ssa::Array, items::Array)

Restructure items for more straightforward translation.

# Details
for example, rewriting expressions like this:
- [mul,a,b,c,d] => [mul,d,[mul,c,[mul,a,b]]]
- [ifselse, cond, a, b] => [select, a, b, cond]
"""
restructure(ci::CodeInfo, i::Integer, ssa::Array, item) = item

#=
function restructure(ci::CodeInfo, i::Integer, ssa::Array, item)
    #use this function to find potential stuff that needs specialized restructuring
    println("restructure, default item: ",item)
    println("restructure, default typeof(item): ",typeof(item))
    return item
end
=#
restructure(ci::CodeInfo, i::Integer, ssa::Array, item::GotoNode) = Any[Expr(:(goto), item.label)]
restructure(ci::CodeInfo, i::Integer, ssa::Array, item::GotoIfNot) = [Expr(:(gotoif), item.dest), ["i32.eqz"; item.cond]]
restructure(ci::CodeInfo, i::Integer, ssa::Array, item::ReturnNode) = Any[Symbol("return"), item.val]

function restructure(ci::CodeInfo, i::Integer, ssa::Array, items::Array)
    if length(items) > 3 && hasname(items[1], keys(floatops))
        # expand N-ary representation ([mul,a,b,c,d] => [mul,d,[mul,c,[mul,a,b]]])
        expanded = items[1:3]
        for i = 4:length(items)
            expanded = [items[1], items[i], expanded]
        end

        return expanded

    elseif hasname(items[1], :(getfield))
        fieldnumber = items[3]
        if isa(items[2], SlotNumber)
            id = items[2].id
            name = fieldnumber == 1 ? "$(ci.slotnames[id])" : "$(ci.slotnames[id])$(fieldnumber)" # get value or index
            return "(local.get \$$name)"
        elseif isa(items[2], SSAValue)
            if isa(ssa[items[2].id], Array)
                # getfield refers to some expression, Julia seems to handle this by a some macro that.. does something
                error("getfield with expression as target, unsure what to return. If tuple assigment like a,b = myfunc(). Solve by doing ab = myfunc() followed by ab[1] or ab[2] when needed. getfield is completely cost free.")
                return ssa[items[2].id][2]
            else
                # getfield refers to an iterator tuple
                target = items[2].id
                id = ssa[target].id
                name = fieldnumber == 1 ? "_$(id)" : "_$(id)i" # get value or index
                # maybe return nothing if fieldnumber == 2.. this makes for x in v work...
                # because the reference that will use up this line is replaced in specializediterate
                # fieldnumber == 2 && return nothing
                return "(local.get \$$name)"
            end
        else
            return "unknown_getfield_on_$(items[2])"
        end

    elseif hasname(items[1], :(iterate))
        return specializediterate(ci, i, ssa, items)
    elseif hasname(items[1], :(ifelse))
        # rewrite ifelse as select [ifselse, cond, a, b] => [select, a, b, cond]
        return [items[1],items[3],items[4],items[2]]
    elseif hasname(items[1], :(:))
        # return items[2:end]
        return nothing
    else
        return restructure.((ci,), (i,), (ssa,), items)
    end
end

function specializediterate(ci, i, ssa, items)
    # single arg => initial value, initial index
    # two args => next value, next index
    ssaref = items[2]
    target = items[2].id
    iterator = ssa[target]
    iteratortype = ci.ssavaluetypes[target]
    if iteratortype <: Array
        if length(items) == 2
            return [GlobalRef(Evalscope, :(iteratearray_init)); ssaref]
        else
            # return [GlobalRef(Evalscope, :(iteratearray)); ssaref; items[3]]
            id = ssa[ssa[items[3].id][2].id].id
            iteratorindex = "(local.get \$_$(id)i)"
            return [GlobalRef(Evalscope, :(iteratearray)); ssaref; iteratorindex]
        end
    # for range functions, pass iteratorvalue instead of iteratorindex as iteratorindex
    elseif iteratortype <: UnitRange
        if length(items) == 2
            return [GlobalRef(Evalscope, :(iterateunitrange_init)); iterator[2:end]]
        else
            # return [GlobalRef(Evalscope, :(iterateunitrange)); iterator[2:end]; items[3]]
            # println("UnitRange, ssa[items[3].id]: ", ssa[items[3].id])
            id = ssa[ssa[items[3].id][2].id].id # walk along some refs to find which slotname it is
            iteratorval = "(local.get \$_$(id))"
            # ssa[items[3].id] = nothing
            return [GlobalRef(Evalscope, :(iterateunitrange)); iterator[2:end]; iteratorval]
        end
    elseif iteratortype <: StepRange
        if length(items) == 2
            return [GlobalRef(Evalscope, :(iteratesteprange_init)); iterator[2:end]]
        else
            # return [GlobalRef(Evalscope, :(iteratesteprange)); iterator[2:end]; items[3]]
            id = ssa[ssa[items[3].id][2].id].id
            iteratorval = "(local.get \$_$(id))"
            return [GlobalRef(Evalscope, :(iteratesteprange)); iterator[2:end]; iteratorval]
        end
    elseif iteratortype <: StepRangeLen
        if length(items) == 2
            return [GlobalRef(Evalscope, :(iteratesteprangelen_init)); iterator[2:end]]
        else
            id = ssa[ssa[items[3].id][2].id].id
            iteratorval = "(local.get \$_$(id))"
            return [GlobalRef(Evalscope, :(iteratesteprangelen)); iterator[2:end]; iteratorval]
        end
    else
        if length(items) == 2
            return [GlobalRef(Evalscope, :(iteratecollection_init)); ssaref]
        else
            return [GlobalRef(Evalscope, :(iteratecollection)); ssaref; items[3]]
        end
    end
end


text:
# 0.18 deprecations. Remove before releasing 0.19
@deprecate path(t::Trie, str::AbstractString) partial_path(t::Trie, str::AbstractString)
@deprecate find_root find_root!
@deprecate top first
@deprecate reverse_iter Iterators.reverse
# Deprecations from #700
Base.@deprecate_binding DisjointSets DisjointSet
Base.@deprecate_binding IntDisjointSets IntDisjointSet
@deprecate DisjointSets(xs...) DisjointSet(xs)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export calcdJdu_CS, calcObjectiveFn, obj_zero, calcVV, calcdRdA_CS

"""
  calculates dJdu using the complex step method
"""
function calcdJdu_CS(mesh::AbstractMesh{Tmsh}, sbp::AbstractSBP, eqn::AbstractSolutionData{Tsol}, opts, h, t) where {Tmsh, Tsol}

  # complex step it
  pert = complex(0, 1e-20)

  dJdu = zeros(Tsol, length(eqn.q_vec))

  for i = 1:length(eqn.q_vec)
    eqn.q_vec[i] += pert

    evalResidual(mesh, sbp, eqn, opts, t)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

    J_arr = calcObjectiveFn(mesh, sbp, eqn, opts, h, t)
    J = J_arr[1]
    # println("=== in dJdu_CS: typeof(J_arr): ", typeof(J_arr))
    # println("=== in dJdu_CS: typeof(J): ", typeof(J))
    # println("=== in dJdu_CS: typeof(integrand_deriv): ", typeof(integrand_deriv))
    # println("=== in dJdu_CS: typeof(pert): ", typeof(pert))
    dJdu[i] = imag(J)/norm(pert)
    eqn.q_vec[i] -= pert
  end

  evalResidual(mesh, sbp, eqn, opts, t)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

  return dJdu

end

"""
  calculates dJdu using the finite difference method
"""
function calcdJdu_FD(mesh::AbstractMesh{Tmsh}, sbp::AbstractSBP, eqn::AbstractSolutionData{Tsol}, opts, h, t) where {Tmsh, Tsol}

  pert = 1e-6

  dJdu = zeros(Tsol, length(eqn.q_vec))

  # unperturbed
  J_unpert_arr = calcObjectiveFn(mesh, sbp, eqn, opts, h, t)
  J_unpert = J_unpert_arr[1]

  for i = 1:length(eqn.q_vec)

    eqn.q_vec[i] += pert

    evalResidual(mesh, sbp, eqn, opts, t)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

    J_arr = calcObjectiveFn(mesh, sbp, eqn, opts, h, t)
    J = J_arr[1]

    dJdu[i] = (J - J_unpert)/norm(pert)
    eqn.q_vec[i] -= pert

  end

  evalResidual(mesh, sbp, eqn, opts, t)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

  return dJdu

end

"""
  calculates the time integral of the objective function,
    storing it in a file each time as well as returning it
"""
function calcTimeIntObjFn(mesh, sbp, eqn, opts, h, t; isDeriv=false, isFirstOrLastStep=false)

  J_total = readdlm(filename)

  if isFirstOrLastStep
    fac = 1.0
  else
    fac = 2.0
  end

  J_i = calcObjectiveFn(mesh, sbp, eqn, opts, h, t; isDeriv)
  J_total += (h*0.5)*fac*J_i

  filename = "J_total.dat"
  writedlm(filename, J_total)

  return J_total

end

"""
  calculates the objective function

  J = int(u^2)

  or optionally (if isDeriv is passed as true), calculates dJdu = int(2u)
"""
function calcObjectiveFn(mesh::AbstractMesh{Tmsh}, sbp::AbstractSBP, eqn::AbstractSolutionData{Tsol}, opts, h, t; isDeriv=false) where {Tmsh, Tsol}

  array1DTo3D(mesh, sbp, eqn, opts, eqn.q, eqn.q_vec)
  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  # TODO: get functional edges in a non BS way
  # 1 corresponds to right side of square domain
  #   Numbering is CCW, with 0 on bottom
  functional_edges = 1
  nDof = 1

  if isDeriv == false
    J = zeros(Tsol, 1)
  else
    dJdu = zeros(Tsol, mesh.numNodesPerElement, mesh.numEl)
  end

  # mesh.bndry_geo_nums:
  #   [1]
  #   [0, 2, 3]

  for itr = 1:length(functional_edges)                  # 1:1
    g_edge_number = functional_edges[itr]               # 1
    itr2 = 0
    for itr2 = 1:mesh.numBC                             # 1:2
      if findfirst(mesh.bndry_geo_nums[itr2], g_edge_number) > 0
        break
      end
    end

    # itr2: 1
    # mesh.bndry_offsets:   [1; 3; 9]
    # mesh.bndryfaces:
    #   Boundary element, face = 7, 1
    #   Boundary element, face = 8, 1
    #   Boundary element, face = 1, 3
    #   Boundary element, face = 1, 1
    #   Boundary element, face = 2, 3
    #   Boundary element, face = 4, 1
    #   Boundary element, face = 6, 2
    #   Boundary element, face = 8, 2

    start_index = mesh.bndry_offsets[itr2]        # start_index = 1
    end_index = mesh.bndry_offsets[itr2+1]        # end_index = 3
    idx_range = start_index:(end_index-1)         # idx_range = 1:2
    bndry_facenums = sview(mesh.bndryfaces, idx_range)

    # bndry_facenums:
    #   Boundary element, face = 7, 1
    #   Boundary element, face = 8, 1

    nfaces = length(bndry_facenums)               # nfaces = 2

    # for scalar field version
    integrand = zeros(Tsol, mesh.sbpface.numnodes, nfaces)     # size: (2, 2)

    if isDeriv == false                 # calculates J = int(u^2)

      for i = 1:nfaces                # 1:2
        bndry_i = bndry_facenums[i]
        global_facenum = idx_range[i]   # i of 1:2

        for j = 1:mesh.sbpface.numnodes       # 1:2
          q = eqn.q_bndry[:, j, global_facenum]

          integrand[j, i] = q[1]*q[1]
          # integrand[j, i] = q[1]
          # println(" ;;;;;;;== q[1]: ", q[1])      # using this, did verify that q[1] = 3 when eqn_poly.q_vec[dof_ix] = 3

          # 20170622: need to scale integrand!
          # TODO: get dimensionality in a non BS way
          dimensionality = 2
          #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          # NOTE as of merge with full code, 20170913:
          #   mesh.jac_bndry[j, i] no longer works with the new curvilinear code.
          #   You need to figure out how to fix it if you want to run this case.
          #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          scaling_factor = 1/(mesh.jac_bndry[j, i]^*(1/dimensionality))
          integrand[j, i] = integrand[j, i]*scaling_factor

        end   # end of loop: j = 1:mesh.sbpfacenumnodes

        # use integratefunctional, not boundaryintegrate: why?

      end   # end of loop: i = 1:nfaces
      J = integratefunctional!(mesh.sbpface, mesh.bndryfaces[idx_range], integrand)

    else                                # calculates dJdu = deriv(int(u^2)) = 2*u

      # for u_ix = 1:length(eqn.q_vec)
        for i = 1:nfaces                # 1:2
          bndry_i = bndry_facenums[i]
          global_facenum = idx_range[i]   # i of 1:2

          for j = 1:mesh.sbpface.numnodes       # 1:2
            q = eqn.q_bndry[:, j, global_facenum]

            integrand[j, i] = 2*q[1]

            # 20170622: need to scale integrand!
            # TODO: get dimensionality in a non BS way
            dimensionality = 2

            #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            # NOTE as of merge with full code, 20170913:
            #   mesh.jac_bndry[j, i] no longer works with the new curvilinear code.
            #   You need to figure out how to fix it if you want to run this case.
            #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            scaling_factor = 1/(mesh.jac_bndry[j, i]^*(1/dimensionality))
            integrand[j, i] = integrand[j, i]*scaling_factor

            # println(" ``````````````````````` j: $j  i: $i  q: ", q[1], "  integrand[j,i]: ", integrand[j,i])

          end   # end of loop: j = 1:mesh.sbpfacenumnodes

          # dJdu[u_ix] = mesh.sbpface.wface[j] * integrand[j,i]
          # use integratefunctional, not boundaryintegrate: why?
        end   # end of loop: i = 1:nfaces
      # end   # end of loop: u_ix = 1:length(eqn.q_vec)

        # dJdu = integratefunctional!(mesh.sbpface, mesh.bndryfaces[idx_range], integrand)
        # dJdu = integratefunctional!(mesh.sbpface, mesh.bndryfaces[idx_range], integrand)

        boundaryintegrate!(mesh.sbpface, mesh.bndryfaces[idx_range], integrand, dJdu)

    end   # end of isDeriv check around loops

  end   # end of loop: itr = 1:length(functional_edges)

  # 20170830 NOTE NOTE: this section is for setting the objective function to a constant (and dJdu to 0)
  #   the result of this is that the adjoint should be zero.
  #   Verified that the adjoint is zero when this is done. I'm leaving it commented out here for posterity,
  #   until the full adjoint works.
  #=
  if isDeriv == false
    println("type of J: ", typeof(J))
    println("size of J: ", size(J))
    J = 4.0
    return J
  else
    dJdu = zeros(dJdu)
    return dJdu
  end
  =#

  if isDeriv == false
    return J
  else
    return dJdu
  end

end

"""
calculates v, which is dudA, for the advection adjoint test
"""
function calcVV(mesh, sbp, eqn, opts, t)
  # eqn (or adj) is unaltered here. In fact, the only thing used
  #   from eqn is eqn.params.omega and eqn.params.sin_amplitude

  v = zeros(eqn.q_vec)

  for dof_ix = 1:mesh.numDof
    # st: subscript_tuple
    st = ind2sub((mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl), dof_ix)
    dofnum = st[1]    # dof number in this node
    nodenum = st[2]   # node number in this element
    elnum = st[3]     # element number

    coords_this_dof = getindex(mesh.coords, :, nodenum, elnum)

    x = coords_this_dof[1]

    omega = eqn.params.omega
    A = eqn.params.sin_amplitude

    # v is dudA. since u = A*sin(-x + omega*t), dudA = sin(-x + omega*t)
    v[dof_ix] = sin(-x + omega*t)

  end

  return v

end

"""
calculates dRdA for the advection adjoint test, using complex step
"""
function calcdRdA_CS(mesh, sbp, eqn, opts, i, t)

  # println(" calculating dRdA using CS")

  # complex step perturbation
  pert = complex(0, 1e-20)

  println(" ~~~~~~~~~~~~~~~~~~~~ eltype(mesh.jac): ", eltype(mesh.jac))

  # println(" <<< dRdA_CS: eqn.params.sin_amplitude: ", eqn.params.sin_amplitude)

  eqn_temp = eqn_deepcopy(mesh, sbp, eqn, opts)
  check_q_qvec_consistency(mesh, sbp, eqn_temp, opts)

  eqn_temp.params.sin_amplitude += pert
  # eqn_temp.params.omega += pert
  # println(" <<< dRdA_CS: eqn_temp.params.sin_amplitude: ", eqn_temp.params.sin_amplitude)

  # Somewhere here, we should see that all elements on the inflow should have a non zero imag portion in res.
  #   because of upwinding, outflow bndry shouldn't have it
  #   because of normal vector w/ adv vel, top & bot shouldn't have it.

  # But we don't for el 1.

  evalResidual(mesh, sbp, eqn_temp, opts, t)
  array3DTo1D(mesh, sbp, eqn_temp, opts, eqn_temp.res, eqn_temp.res_vec)

  dRdA = imag(eqn_temp.res_vec)/norm(pert)
  # println("dRdA: ", dRdA)

  eqn_temp.params.sin_amplitude -= pert
  # eqn_temp.params.omega -= pert

  return dRdA
end

"""
calculates dRdA for the advection adjoint test, using finite difference
"""
function calcdRdA_FD(mesh, sbp, eqn, opts, i, t)

  # println(" calculating dRdA using FD")

  # finite difference perturbation
  pert = 1e-6

  # println(" <<< dRdA_FD: eqn.params.sin_amplitude: ", eqn.params.sin_amplitude)

  # filename_eqn = string("dRdA_FD_eqn_i-",i,".dat")
  # filename_eqn_temp = string("dRdA_FD_eqn_temp_i-",i,".dat")
  # writedlm(filename_eqn, eqn.q_vec)

  eqn_temp = eqn_deepcopy(mesh, sbp, eqn, opts)
  check_q_qvec_consistency(mesh, sbp, eqn_temp, opts)
  # verified: eqn_temp.q_vec & eqn.q_vec values match
  if norm(eqn_temp.q_vec - eqn.q_vec) > 1e-4
    error("Problem in calcdRdA_FD's eqn_deepcopy")
  end

  # writedlm(filename_eqn_temp, eqn_temp.q_vec)
  println("     ~~~~~~~~~~~~~~~ in calcdRdA_FD ~~~~~~~~~~~~~~")
  println("       eqn_temp.params.sin_amplitude = ", eqn_temp.params.sin_amplitude)
  println("       eqn.params.sin_amplitude = ", eqn.params.sin_amplitude)
  println("       ... perturbing eqn_temp.params.sin_amplitude")
  eqn_temp.params.sin_amplitude += pert
  println("       eqn_temp.params.sin_amplitude = ", eqn_temp.params.sin_amplitude)
  println("       eqn.params.sin_amplitude = ", eqn.params.sin_amplitude)
  # eqn_temp.params.omega += pert
  # println(" <<< dRdA_FD: eqn_temp.params.sin_amplitude: ", eqn_temp.params.sin_amplitude)

  # println(" <<< dRdA_FD: using forward difference")

  evalResidual(mesh, sbp, eqn_temp, opts, t)
  array3DTo1D(mesh, sbp, eqn_temp, opts, eqn_temp.res, eqn_temp.res_vec)

  eqn_temp.params.sin_amplitude -= pert
  # eqn_temp.params.omega -= pert

  # this is required! without it, eqn.res_vec apparently has stale data.
  #   without this => difference between CS & FD is ~1e10
  #   TODO: why is this required
  evalResidual(mesh, sbp, eqn, opts, t)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  # TODO TODO 20170908

  # println(" <<< dRdA_FD: norm(eqn_temp.res_vec): ", norm(eqn_temp.res_vec))
  # println(" <<< dRdA_FD: norm(eqn.res_vec): ", norm(eqn.res_vec))

  # Ready to finish FD approx.
  #   f(x+h) is eqn_temp.res_vec
  #   f(x) is eqn.res_vec

  dRdA = (eqn_temp.res_vec - eqn.res_vec)/norm(pert)
  #println("dRdA: ", dRdA)

  #=
  #----------------
  # central difference
  eqn_temp_fwd = eqn_deepcopy(mesh, sbp, eqn, opts)
  eqn_temp_fwd.params.sin_amplitude += pert

  eqn_temp_bck = eqn_deepcopy(mesh, sbp, eqn, opts)
  eqn_temp_bck.params.sin_amplitude -= pert

  println(" <<< dRdA_FD: using central difference")

  evalResidual(mesh, sbp, eqn_temp_fwd, opts)
  array3DTo1D(mesh, sbp, eqn_temp_fwd, opts, eqn_temp_fwd.res, eqn_temp_fwd.res_vec)
  evalResidual(mesh, sbp, eqn_temp_bck, opts)
  array3DTo1D(mesh, sbp, eqn_temp_bck, opts, eqn_temp_bck.res, eqn_temp_bck.res_vec)

  eqn_temp_fwd.params.sin_amplitude -= pert
  eqn_temp_bck.params.sin_amplitude += pert

  dRdA = (eqn_temp_fwd.res_vec - eqn_temp_bck.res_vec)/(2*norm(pert))
  =#

  return dRdA
end

"""
  obj_zero

  Inputs: none
  Outputs: 0

  Zero-valued objective function
"""
function obj_zero()
  return 0.0
end


text:
using IdentityMatrix
using Test
using LinearAlgebra
using FillArrays
using MethodInSrc

@testset  "norms" begin
    for ncols in (0, 1, 3, 10)
        for T in (Float64, Int)
            m = Eye{T}(ncols)
            md = Matrix(m)
            for op in (opnorm, norm, x -> norm(x, 1), x -> norm(x, 2), x -> norm(x, Inf), x -> norm(x, -Inf), x -> norm(x, 3//2))
                @test op(m) == op(md)
            end
            @test @isinsrc opnorm(m)
        end
    end
end

@testset "Eye" begin
    for ncols in (1, 2, 5)
        IM = Eye(ncols)
        IMd = Matrix(IM)
        @test isa(imag(IM), Diagonal)
        for f in (imag, isposdef, sum, prod, first, last,
                  minimum, maximum, extrema, triu, triu!, tril, tril!, inv,
                  diag, det, logdet, sqrt)
            @test @insrc(f(IM)) == f(copy(IMd))
        end
        for f in (iszero, )
            @test @ninsrc(f(IM)) == f(copy(IMd))
        end
        for p in -5:5
            @test IM == @insrc IM^p
        end
        @test isa(copy(IM), Diagonal)
        @test isa(@insrc(Matrix(IM)), Matrix)
    end
end

@testset "mul and div" begin
    for ncols in (1, 3 ,10)
        IM = Eye(ncols)
        @test IM === @insrc IM * IM
        @test IM === @insrc IM / IM
        @test IM === @insrc IM \ IM

        rm = rand(ncols, ncols)
        @test rm === @insrc IM * rm
        @test rm === @insrc rm * IM
        @test rm === @insrc rm / IM
        @test inv(rm) == @insrc IM / rm
        @test inv(rm) == @insrc rm \ IM
        @test rm == @insrc IM \ rm
    end
end

@testset "algebra with AbstractFill" begin
    for d in (1, 2, 3, 10)
       for  T1 in (Int, Float64, ComplexF64)
           for T2 in (Int, Float64, ComplexF64)
               @test @insrc(Eye{T1}(d) * Zeros{T2}(d, d)) == Zeros{typeof(one(T1)*one(T2))}(d, d)
               @test @insrc(Eye{T1}(d) * Zeros{T2}(d)) == Zeros{typeof(one(T1)*one(T2))}(d)
               @test @insrc(Eye{T1}(d) * Ones{T2}(d, d)) == Ones{typeof(one(T1)*one(T2))}(d, d)
               @test @insrc(Eye{T1}(d) * Ones{T2}(d)) == Ones{typeof(one(T1)*one(T2))}(d)
               for val in (1, 2)
                   @test @insrc(Eye{T1}(d) * Fill{T2}(val, d, d)) == Fill{typeof(one(T1)*one(T2))}(val, d, d)
                   @test @insrc(Eye{T1}(d) * Fill{T2}(val, d)) == Fill{typeof(one(T1)*one(T2))}(val, d)
               end
           end
       end
    end
end

@testset "any all iszero isone" begin
    for T in (Int, Float64, ComplexF64)
        d1 = 1
        m = Eye{T}(d1)
        @test ! @inmodule FillArrays any(iszero, m)
        @test ! @ninsrc all(iszero, m)
        @test @ninsrc any(isone, m)
        @test @ninsrc all(isone, m)

        onem = Ones{T}(d1, d1)
        @test @ninsrc isone(onem)
        @test ! iszero(onem)

        zerom = Zeros{T}(d1, d1)
        @test ! isone(zerom)
        @test  iszero(zerom)

        fillm0 = Fill(T(0), d1, d1)
        @test ! isone(fillm0)
        @test   iszero(fillm0)

        fillm1 = Fill(T(1), d1, d1)
        @test isone(fillm1)
        @test ! iszero(fillm1)

        fillm2 = Fill(T(2), d1, d1)
        @test ! isone(fillm2)
        @test ! iszero(fillm2)

        for d in (2, 3)
            m = Eye{T}(d)
            @test any(iszero, m)
            @test ! all(iszero, m)
            @test any(isone, m)
            @test ! all(isone, m)

            m1 = Ones{T}(d, d)
            @test ! isone(m1)
            @test ! iszero(m1)
            @test all(isone, m1)
            @test ! all(iszero, m1)

            m2 = Zeros{T}(d, d)
            @test ! isone(m2)
            @test iszero(m2)
            @test ! all(isone, m2)
            @test  all(iszero, m2)

            m3 = Fill(T(2), d, d)
            @test ! isone(m3)
            @test ! iszero(m3)
            @test ! all(isone, m3)
            @test ! all(iszero, m3)
        end
    end
end

@testset "reduction types" begin
    for ncols in (1, 3, 10)
        for T in (Int, Float64, Int32, BigInt) # Rational{Int})
            for fop in (sum, prod)
                mi = Eye{T}(ncols)
                midense = Matrix(mi)
                r = @insrc fop(mi)
                rdense = fop(midense)
                @test r == rdense
                @test typeof(r) == eltype(mi)
            end
        end
    end
end

@testset "transform" begin
    m = Eye(3)
    for f in (triu, triu!, tril, tril!, inv)
        @test m == @insrc f(m) # FIXME: why does === fail ?
    end
end

@testset "kron Eye" begin
    for ncols1 in (1, 2, 3, 4)
        for ncols2 in (1, 2, 3, 4)
            im = Eye(ncols1)
            imd = Matrix(im)
            md = rand(ncols2, ncols2)
            @test kron(imd, md) == @insrc kron(im, md)
            @test kron(md, imd) == @insrc kron(md, im)
            @test isone(@isinsrc(kron(im, im)))
            @test isa(kron(im, im), typeof(im))
            @test size(@insrc(kron(im, zeros(0, 0)))) == (0,0)
            @test size(kron(zeros(0, 0), im)) == (0,0)
        end
    end
end

@testset "kron Diagonal" begin
    Ntests = 10
    for m in (1, 2, 3, 4)
        for m1 in (1, 2, 3, 4, 7)
            for n1 in (1, 2, 3, 4, 7)
                for icount in 1:Ntests
                    md = Diagonal(rand(m))
                    mdd = Matrix(md)
                    m2 = rand(m1, n1)
                    @test kron(md, m2) == kron(mdd, m2)
                    @test kron(m2, md) == kron(m2, mdd)
                end
            end
        end
    end
end

@testset "eigen" begin
    for ncols in (0, 1, 3, 5)
        im = Eye(ncols)
        md = Matrix(im)
        ri = @insrc eigen(im)
        rd = eigen(md)
        @test isapprox(ri.values, rd.values)
        @test isapprox(ri.vectors, rd.vectors)
        @test isapprox(eigvals(md), @insrc(eigvals(im)))
        @test isapprox(eigvecs(im), eigvecs(md))
        @test isa(ri.vectors, typeof(im))
    end
end

@testset "one oneunit zero" begin
    for ncols in (0, 1, 3, 10)
        for T in (Float64, Int)
            m = Eye{T}(ncols)
            @test one(m) === m
            @test oneunit(m) === m
            @test zero(m) == Diagonal(Zeros{eltype(m)}(size(m, 1)))
        end
    end
end

@testset "Fill" begin
    for ncols in (1, 3, 10)
        for T in (Float64,)
            for val in (1, 2, 2.5, 10)
                fm = Diagonal(Fill{T}(val, ncols))
                fmd = Matrix(fm)
                @test inv(fm) == inv(fmd)
            end
        end
    end
end


text:
#####
##### Profiling Self-play
#####

# using Revise; Revise.includet("scripts/profile/self_play.jl")
# profile_self_play()

using AlphaZero
using Setfield

using Profile
using ProfileSVG
# using ProfileView
# import StatProfilerHTML

import CUDA

function profile_self_play(
  exp::Experiment = Examples.experiments["connect-four"];
  num_games=512,
  num_workers=256,
  batch_size=128,
  num_filters=64)

  exp = @set exp.netparams.num_filters = num_filters
  exp = @set exp.params.self_play.sim.num_workers = num_workers
  exp = @set exp.params.self_play.sim.batch_size = batch_size
  exp = @set exp.params.self_play.sim.num_games = num_games

  session = Session(exp, autosave=false, dir="sessions/profile-backprop-$(exp.name)")
  env = session.env
  UI.Log.section(session.logger, 1, "Profiling data generation")
  CUDA.@time AlphaZero.self_play_step!(env, session)
  return
end

function full_profile()
  profile_self_play() # Compilation
  profile_self_play()
  profile_self_play() # Double-checking
end

function flame_graph()
  profile_self_play() # Compilation
  Profile.init(n=10^8, delay=0.01)
  Profile.clear()
  @profile profile_self_play()
  ProfileSVG.save("self-play-profile.svg")
  # StatProfilerHTML.statprofilehtml()
  # ProfileView.view(C=true)
end

full_profile()

# Num cores experiments
# 1: 203s / 13% GC
# 2: 156s / 16% GC
# 3: 131s / 19% GC
# 4: 128s / 20% GC
# 6: 125s / 20% GC

# Knet vs Flux
# 1. Flux, 64 filters:  92 samples/s  (20% in GC)
# 2. Flux, 128 filters: 52 samples/s
# 3. Knet, 64 filters:  82 samples/s  (35% in GC)
# 4. Knet, 128 filters: 42 samples/s


text:
export Comment

"""
    Comment(text::AbstractString)

`Comment` objects `show` a FastHenry comment.
"""
struct Comment <: Element
  text :: String
end