text:
@testset "IO HDF5 PointCloud" begin
    filename = joinpath(tempdir(), "test_pc.h5")

    try
        N = 1000
        pos = rand(SVector{3,Float32}, N)
        int = rand(UInt16, N)
        class = zeros(Int, N)
        colors = rand(RGB{N0f8}, N)

        pc = Table(position = pos, intensity = int, classification = class, color = colors)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename; spacing = 0.1)
        @test pc2.position == pc.position
        @test pc2.intensity == pc.intensity
        @test pc2.classification == pc.classification
        @test pc2.color == pc.color
        @test pc2.position isa AcceleratedArray

        colors2 = rand(RGB{N0f16}, N)
        pc = Table(position = pos, intensity = int, classification = class, color = colors2)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename)
        @test pc2.position == pc.position
        @test pc2.intensity == pc.intensity
        @test pc2.classification == pc.classification
        @test pc2.color == pc.color
        @test !(pc2.position isa AcceleratedArray)

        pc = empty(pc)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename)
        @test pc2 == pc
    finally
        rm(filename, force=true)
    end
end

@testset "IO HDF5 PointCloud - GeoRepo2 format" begin
    filename = joinpath(tempdir(), "test_pc2.h5")

    try
        N = 1000
        pos = rand(SVector{3,Float32}, N)
        int = rand(UInt16, N)
        cluster = zeros(Int, N)
        colors = rand(RGB{N0f8}, N)

        pc = Table(position = pos, intensity = int, clusterid = cluster, color = colors)
        save_pointcloud(filename, pc; format = "XYZIUrgb")
        pc2 = load_pointcloud(filename; spacing = 0.1, format = "XYZIUrgb")
        @test pc2.position == pc.position
        @test pc2.intensity == pc.intensity
        @test pc2.clusterid == pc.clusterid
        @test pc2.color == pc.color
        @test pc2.position isa AcceleratedArray

        colors2 = rand(RGB{N0f16}, N)
        pc = Table(position = pos, intensity = int, clusterid = cluster, color = colors2)
        save_pointcloud(filename, pc; format = "XYZIUrgb")
        pc2 = load_pointcloud(filename; format = "XYZIUrgb")
        @test pc2.position == pc.position
        @test pc2.intensity == pc.intensity
        @test pc2.clusterid == pc.clusterid
        @test pc2.color == convert(Vector{RGB{N0f8}}, pc.color) # color is saved in 8-bit format
        @test !(pc2.position isa AcceleratedArray)

        pc = empty(pc)
        save_pointcloud(filename, pc; format = "XYZIUrgb")
        pc2 = load_pointcloud(filename; format = "XYZIUrgb")
        @test pc2 == pc
    finally
        rm(filename, force=true)
    end
end

@testset "IO LAS PointCloud" begin
    filename = joinpath(tempdir(), "test_pc.las")

    try
        N = 1000
        pos = rand(SVector{3,Float32}, N)
        int = rand(UInt16, N)
        class = zeros(Int, N)
        colors = rand(RGB{N0f8}, N)

        pc = Table(position = pos, intensity = int, classification = class, color = colors)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename; spacing = 0.1)
        @test all(((p1, p2),) -> isapprox(p1, p2; atol = 1e-3), zip(pc.position, pc2.position))
        @test pc2.intensity == pc.intensity
        @test pc2.classification == pc.classification
        @test pc2.color == pc.color
        @test pc2.position isa AcceleratedArray

        colors2 = rand(RGB{N0f16}, N)
        pc = Table(position = pos, intensity = int, classification = class, color = colors2)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename)
        @test all(((p1, p2),) -> isapprox(p1, p2; atol = 1e-3), zip(pc.position, pc2.position))
        @test pc2.intensity == pc.intensity
        @test pc2.classification == pc.classification
        @test pc2.color == pc.color
        @test !(pc2.position isa AcceleratedArray)

        pc = empty(pc)
        save_pointcloud(filename, pc)
        pc2 = load_pointcloud(filename)
        @test pc2.position == pc.position
        @test pc2.intensity == pc.intensity
        @test pc2.classification == pc.classification
        @test_broken pc2.color == pc.color # empty point cloud defaults to point type 0, should fix at some point
    finally
        rm(filename, force=true)
    end
end


text:
using YAML
using DataStructures
using FileIO

# Functions for parsing and writing dsl files
export parse_dsl, parse_dsl_files, generate_dsl_files

# Define compatible DSL file version number, which must match when parsed.
const DSL_VERSION = VersionNumber("0.4.0")

###############################################################################
# DSL Parsing functions
###############################################################################
"""
    write_version_header(io::IO)

Function to write version head to DSL file with current version constant
"""
function write_version_header(io::IO)
    write(io, "# version: $(DSL_VERSION)\n")
end

"""
    check_compatible_version_dsl(line::String)

Checks if version is defined in line and checks compatibility with VERSION_DSL
"""
function check_compatible_version_dsl(line::String)
    exists_version_dsl = startswith(strip(line), '#') && occursin("version:", line)
    is_compatible::Bool = true

    if exists_version_dsl
        version_dsl = strip(last(split(line,"version:")))
        version_dsl = VersionNumber(version_dsl)

        # Simple logic enforcing matching versions, which can be extended
        is_compatible = version_dsl == DSL_VERSION
    else
        is_compatible = false
        version_dsl = nothing
    end

    return is_compatible, version_dsl
end

"""
    parse_command(line::String)

Parse a DSL command
"""
function parse_command(line::AbstractString)
    m = match(r"^([a-z]*)", line)
    command = nothing
    if m !== nothing
        type = m.captures[1]
        command = if type == "load" LoadCommand(line)
        elseif type == "view" command = ViewCommand(line)
        elseif type == "ncon" command = ContractCommand(line)
        elseif type == "save" command = SaveCommand(line)
        elseif type == "output" command = OutputCommand(line)
        elseif type == "reshape" command = ReshapeCommand(line)
        else
            error("$(type) command has not been implemented yet")
        end
    end
    command
end

"""
    parse_dsl(buffer::Vector{String})

Parse a list of DSL commands and generate a CommandList for execution
"""
function parse_dsl(buffer::Vector{<:AbstractString})
    line = string(strip(first(buffer)))
    is_compatible, version_dsl = check_compatible_version_dsl(line)
    if !is_compatible
        throw(ArgumentError("DSL version not compatible:\n\t'$version_dsl', expected '$DSL_VERSION'"))
    end

    # find index of first line that doesn't start with "#"
    cmd_idx = findfirst(x -> x[1] != '#', buffer)
    metadata_str = join(map(x -> replace(x, r"^# " => ""), buffer[2:cmd_idx-1]), "\n")
    metadata = length(metadata_str) > 0 ? YAML.load(metadata_str) : OrderedDict()

    cmds = Vector{AbstractCommand}()
    for line in buffer[cmd_idx:end]
        line_command = string(first(split(line, '#')))
        if !isempty(line_command)
            push!(cmds, parse_command(line_command))
        end
    end
    build_tree(cmds), metadata
end

"""
    parse_dsl(filename::String)

Read a DSL file and generate a CommandList for execution
"""
function parse_dsl(filename::String)
    return open(filename) do file
        parse_dsl(readlines(file))
    end
end

"""
    parse_dsl_files(dsl_file::String, data_file::String)

Read a DSL file and tensors file and return compute tree and meta data
"""
function parse_dsl_files(dsl_file::String, data_file::String)
    root_node, metadata = parse_dsl(dsl_file)
    @assert splitext(data_file)[end] == ".jld2" "Data file should have \".jld2\" suffix"
    tensors = Dict(Symbol(x) => y for (x, y) in pairs(load(data_file)))
    ComputeGraph(root_node, tensors), metadata
end

###############################################################################
# DSL writing functions
###############################################################################

"""
    Base.write(io::IO, Union{ComputeNode, ComputeGraph}; metadata=nothing)

Write the compute tree to an ascii file. Optionally write meta_data
if present
"""
function Base.write(io::IO, cn::ComputeNode; metadata=nothing)
    # first we write the version header
    write_version_header(io)

    if metadata !== nothing
        # prepend each link of yaml output with "# "
        yml_str = YAML.write(metadata)
        yml_str = join(["# " * x for x in split(yml_str, "\n")], "\n") * "\n"
        write(io, yml_str)
    end

    for each in PostOrderDFS(cn)
        write(io, each.op)
    end
end

"""
    generate_dsl_files(compute_tree::ComputeGraph,
                       prefix::String;
                       force::Bool=true,
                       metadata=nothing)

Function to create dsl and data files to contract the given tensor network circuit
with the plan provided
"""
function generate_dsl_files(compute_tree::ComputeGraph,
                            prefix::String;
                            force::Bool=true,
                            metadata=nothing)

    dsl_filename = "$(prefix).qx"
    data_filename = "$(prefix).jld2"

    @assert force || !isfile(dsl_filename) "Error $(dsl_filename) already exists"
    @assert force || !isfile(data_filename) "Error $(data_filename) already exists"

    open(dsl_filename, "w") do dsl_io
        write(dsl_io, compute_tree.root; metadata=metadata)
    end

    save(data_filename, Dict(String(x) => y for (x, y) in pairs(compute_tree.tensors)))
    nothing
end

ext: jl


text:
using ColormapsOcean
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
cmaps=[:algae,:amp,:balance,:curl,:deep,:delta,:dense,:gray,:haline,:ice,:matter,:oxy,:phase,:solar ,:speed,:tempo,:thermal,:turbid]
sz=[size(cmocean[cm],1) for cm in cmaps]
@test size(sz, 1)==1
@test all([z==256 || z==512 for z in sz])


text:
#
#--------------------------------------#
export meshplt
#--------------------------------------#
function meshplt(x,y,u;a=45,b=60)
    p = plot(x,y,u,legend=false,c=:grays,camera=(a,b))
    p = plot!(x',y',u',legend=false,c=:grays,camera=(a,b))
    return p
end
#--------------------------------------#
function meshplt(u,msh::Mesh;a=45,b=60)
    p = meshplt(msh.x,msh.y,u;a=a,b=b)
    return p
end
#--------------------------------------#
function meshplt(f::Field,msh::Mesh;a=45,b=60)
    p = meshplt(f.u,msh;a=a,b=b)
    return p
end
#--------------------------------------#


text:
module Insolation

using Dates, DelimitedFiles, Interpolations
using CLIMAParameters
using CLIMAParameters.Planet
using CLIMAParameters: AbstractParameterSet
const APS = AbstractParameterSet

export orbital_params

mock_t_range = 1.0:1.0:10.0;
mock_x_data = rand(10);
const e_spline_etp = Ref(CubicSplineInterpolation(mock_t_range, mock_x_data, extrapolation_bc = NaN))
const γ_spline_etp = Ref(CubicSplineInterpolation(mock_t_range, mock_x_data, extrapolation_bc = NaN))
const ϖ_spline_etp = Ref(CubicSplineInterpolation(mock_t_range, mock_x_data, extrapolation_bc = NaN))
function __init__()
    datapath = joinpath(@__DIR__, "../src/data/", "INSOL.LA2004.BTL.csv");
    x, _ = readdlm(datapath, ',', Float64, header=true);
    t_range = x[1,1]*1e3 : 1e3 : x[end,1]*1e3; # array of every 1 kyr to range of years
    e_spline_etp[] = CubicSplineInterpolation(t_range, x[:,2], extrapolation_bc = NaN);
    γ_spline_etp[] = CubicSplineInterpolation(t_range, x[:,3], extrapolation_bc = NaN);
    ϖ_spline_etp[] = CubicSplineInterpolation(t_range, x[:,4], extrapolation_bc = NaN);
end
function e_spline(args...)
    return e_spline_etp[](args...)
end
function γ_spline(args...)
    return γ_spline_etp[](args...)
end
function ϖ_spline(args...)
    return ϖ_spline_etp[](args...)
end

"""
    orbital_params(dt::FT) where {FT <: Real}

This function returns the orbital parameters (ϖ, γ, e) at dt (years) since J2000 epoch.
Data are read from file and interpolation functions are created in __init__() method.
The functions are stored as global variables that are used inside Insolation.jl.
The parameters vary due to Milankovitch cycles.
Orbital parameters from the Laskar 2004 paper are in the "src/data/INSOL.LA2004.BTL.csv" file.
"""
function orbital_params(dt::FT) where {FT <: Real}
    return ϖ_spline(dt), γ_spline(dt), e_spline(dt)
end

include("ZenithAngleCalc.jl")
include("InsolationCalc.jl")

end # module


text:
using Unicode
function bacs()
  bulls = cows = turns = 0
  a = Any[]
  while length(unique(a)) < 4
    push!(a,rand('0':'9'))
    end
  my_guess = unique(a)
  print("\nBULLS and COWS\n===============")
  prompt = "\nEnter four distinct digits or <return> to quit"
  println(prompt)
  while (bulls != 4)
    print("Guess> ");
    s = chomp(readline())
    if (length(s) == 0)
      print("My guess was "); [print(my_guess[i]) for i=1:4]
      return
    end
    guess = collect(s)
    if !(length(unique(guess)) == length(guess) == 4 && all(Unicode.isdigit,guess))
      println(prompt)
      continue
    end
    bulls = sum(map(==, guess, my_guess))
    cows = length(intersect(guess,my_guess)) - bulls
    println("$bulls bulls and $cows cows!")
    turns += 1
  end
  println("\nYou guessed my number in $turns turns.")
end

BULLS and COWS
===============
Enter four distinct digits or <return> to quit
Guess> 1234
0 bulls and 1 cows!
Guess> 5678
0 bulls and 1 cows!
Guess> 1590
2 bulls and 0 cows!
Guess> 2690
2 bulls and 0 cows!
Guess> 3790
2 bulls and 0 cows!
Guess> 4890
2 bulls and 2 cows!
Guess> 8490
4 bulls and 0 cows!

You guessed my number in 7 turns.


text:
using Mimi

@defcomp RegionAggregatorSum begin
    countries = Dimension()
    regions = Dimension()

    region_mapping::Dict{Symbol,Vector{Symbol}} = Parameter()
    input = Parameter(index=[time,inputregions])
    output = Parameter(index=[time,outputregions])
end


text:
"""
    stack(df::AbstractDataFrame[, measure_vars[, id_vars] ];
          variable_name=:variable, value_name=:value,
          view::Bool=false, variable_eltype::Type=String)

Stack a data frame `df`, i.e. convert it from wide to long format.

Return the long-format `DataFrame` with: columns for each of the `id_vars`,
column `value_name` (`:value` by default)
holding the values of the stacked columns (`measure_vars`), and
column `variable_name` (`:variable` by default) a vector holding
the name of the corresponding `measure_vars` variable.

If `view=true` then return a stacked view of a data frame (long format).
The result is a view because the columns are special `AbstractVectors`
that return views into the original data frame.

# Arguments
- `df` : the AbstractDataFrame to be stacked
- `measure_vars` : the columns to be stacked (the measurement variables),
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  If neither `measure_vars` or `id_vars` are given, `measure_vars`
  defaults to all floating point columns.
- `id_vars` : the identifier columns that are repeated during stacking,
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  Defaults to all variables that are not `measure_vars`
- `variable_name` : the name (`Symbol` or string) of the new stacked column that
  shall hold the names of each of `measure_vars`
- `value_name` : the name (`Symbol` or string) of the new stacked column containing
  the values from each of `measure_vars`
- `view` : whether the stacked data frame should be a view rather than contain
  freshly allocated vectors.
- `variable_eltype` : determines the element type of column `variable_name`.
  By default a `PooledArray{String}` is created.
  If `variable_eltype=Symbol` a `PooledVector{Symbol}` is created,
  and if `variable_eltype=CategoricalValue{String}`
  a `CategoricalArray{String}` is produced (call `using CategoricalArrays` first if needed)
  Passing any other type `T` will produce a `PooledVector{T}` column
  as long as it supports conversion from `String`.
  When `view=true`, a `RepeatedVector{T}` is produced.

# Examples
```jldoctest
julia> df = DataFrame(a=repeat(1:3, inner=2),
                      b=repeat(1:2, inner=3),
                      c=repeat(1:1, inner=6),
                      d=repeat(1:6, inner=1),
                      e=string.('a':'f'))
6×5 DataFrame
 Row │ a      b      c      d      e
     │ Int64  Int64  Int64  Int64  String
─────┼────────────────────────────────────
   1 │     1      1      1      1  a
   2 │     1      1      1      2  b
   3 │     2      1      1      3  c
   4 │     2      2      1      4  d
   5 │     3      2      1      5  e
   6 │     3      2      1      6  f

julia> stack(df, [:c, :d])
12×5 DataFrame
 Row │ a      b      e       variable  value
     │ Int64  Int64  String  String    Int64
─────┼───────────────────────────────────────
   1 │     1      1  a       c             1
   2 │     1      1  b       c             1
   3 │     2      1  c       c             1
   4 │     2      2  d       c             1
   5 │     3      2  e       c             1
   6 │     3      2  f       c             1
   7 │     1      1  a       d             1
   8 │     1      1  b       d             2
   9 │     2      1  c       d             3
  10 │     2      2  d       d             4
  11 │     3      2  e       d             5
  12 │     3      2  f       d             6

julia> stack(df, [:c, :d], [:a])
12×3 DataFrame
 Row │ a      variable  value
     │ Int64  String    Int64
─────┼────────────────────────
   1 │     1  c             1
   2 │     1  c             1
   3 │     2  c             1
   4 │     2  c             1
   5 │     3  c             1
   6 │     3  c             1
   7 │     1  d             1
   8 │     1  d             2
   9 │     2  d             3
  10 │     2  d             4
  11 │     3  d             5
  12 │     3  d             6

julia> stack(df, Not([:a, :b, :e]))
12×5 DataFrame
 Row │ a      b      e       variable  value
     │ Int64  Int64  String  String    Int64
─────┼───────────────────────────────────────
   1 │     1      1  a       c             1
   2 │     1      1  b       c             1
   3 │     2      1  c       c             1
   4 │     2      2  d       c             1
   5 │     3      2  e       c             1
   6 │     3      2  f       c             1
   7 │     1      1  a       d             1
   8 │     1      1  b       d             2
   9 │     2      1  c       d             3
  10 │     2      2  d       d             4
  11 │     3      2  e       d             5
  12 │     3      2  f       d             6

julia> stack(df, Not([:a, :b, :e]), variable_name=:somemeasure)
12×5 DataFrame
 Row │ a      b      e       somemeasure  value
     │ Int64  Int64  String  String       Int64
─────┼──────────────────────────────────────────
   1 │     1      1  a       c                1
   2 │     1      1  b       c                1
   3 │     2      1  c       c                1
   4 │     2      2  d       c                1
   5 │     3      2  e       c                1
   6 │     3      2  f       c                1
   7 │     1      1  a       d                1
   8 │     1      1  b       d                2
   9 │     2      1  c       d                3
  10 │     2      2  d       d                4
  11 │     3      2  e       d                5
  12 │     3      2  f       d                6
```
"""
function stack(df::AbstractDataFrame,
               measure_vars = findall(col -> eltype(col) <: Union{AbstractFloat, Missing},
                                      eachcol(df)),
               id_vars = Not(measure_vars);
               variable_name::SymbolOrString=:variable,
               value_name::SymbolOrString=:value, view::Bool=false,
               variable_eltype::Type=String)
    variable_name_s = Symbol(variable_name)
    value_name_s = Symbol(value_name)
    # getindex from index returns either Int or AbstractVector{Int}
    mv_tmp = index(df)[measure_vars]
    ints_measure_vars = mv_tmp isa Int ? [mv_tmp] : mv_tmp
    idv_tmp = index(df)[id_vars]
    ints_id_vars = idv_tmp isa Int ? [idv_tmp] : idv_tmp
    if view
        return _stackview(df, ints_measure_vars, ints_id_vars,
                          variable_name=variable_name_s,
                          value_name=value_name_s,
                          variable_eltype=variable_eltype)
    end
    N = length(ints_measure_vars)
    cnames = _names(df)[ints_id_vars]
    push!(cnames, variable_name_s)
    push!(cnames, value_name_s)
    if variable_eltype === Symbol
        catnms = PooledArray(_names(df)[ints_measure_vars])
    elseif variable_eltype === String
        catnms = PooledArray(names(df, ints_measure_vars))
    else
        # this covers CategoricalArray{String} in particular
        # (note that copyto! inserts levels in their order of appearance)
        nms = names(df, ints_measure_vars)
        simnms = similar(nms, variable_eltype)
        catnms = simnms isa Vector ? PooledArray(catnms) : simnms
        copyto!(catnms, nms)
    end
    return DataFrame(AbstractVector[[repeat(df[!, c], outer=N) for c in ints_id_vars]..., # id_var columns
                                    repeat(catnms, inner=nrow(df)),                       # variable
                                    vcat([df[!, c] for c in ints_measure_vars]...)],      # value
                     cnames, copycols=false)
end

function _stackview(df::AbstractDataFrame, measure_vars::AbstractVector{Int},
                    id_vars::AbstractVector{Int}; variable_name::Symbol,
                    value_name::Symbol, variable_eltype::Type)
    N = length(measure_vars)
    cnames = _names(df)[id_vars]
    push!(cnames, variable_name)
    push!(cnames, value_name)
    if variable_eltype === Symbol
        catnms = _names(df)[measure_vars]
    elseif variable_eltype === String
        catnms = names(df, measure_vars)
    else
        # this covers CategoricalArray{String} in particular,
        # as copyto! inserts levels in their order of appearance
        nms = names(df, measure_vars)
        catnms = copyto!(similar(nms, variable_eltype), nms)
    end
    return DataFrame(AbstractVector[[RepeatedVector(df[!, c], 1, N) for c in id_vars]..., # id_var columns
                                    RepeatedVector(catnms, nrow(df), 1),                  # variable
                                    StackedVector(Any[df[!, c] for c in measure_vars])],  # value
                     cnames, copycols=false)
end

"""
    unstack(df::AbstractDataFrame, rowkeys, colkey, value; renamecols::Function=identity,
            allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing)
    unstack(df::AbstractDataFrame, colkey, value; renamecols::Function=identity,
            allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing)
    unstack(df::AbstractDataFrame; renamecols::Function=identity,
            allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing)

Unstack data frame `df`, i.e. convert it from long to wide format.

Row and column keys will be ordered in the order of their first appearance.

# Positional arguments
- `df` : the AbstractDataFrame to be unstacked
- `rowkeys` : the columns with a unique key for each row, if not given,
  find a key by grouping on anything not a `colkey` or `value`.
  Can be any column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  If `rowkeys` contains no columns all rows are assumed to have the same key.
- `colkey` : the column ($COLUMNINDEX_STR) holding the column names in wide format,
  defaults to `:variable`
- `value` : the value column ($COLUMNINDEX_STR), defaults to `:value`

# Keyword arguments

- `renamecols`: a function called on each unique value in `colkey`; it must return
  the name of the column to be created (typically as a string or a `Symbol`).
  Duplicates in resulting names when converted to `Symbol` are not allowed.
  By default no transformation is performed.
- `allowmissing`: if `false` (the default) then an error will be thrown if `colkey`
  contains `missing` values; if `true` then a column referring to `missing` value
  will be created.
- `allowduplicates`: if `false` (the default) then an error an error will be thrown
  if combination of `rowkeys` and `colkey` contains duplicate entries; if `true`
  then  then the last encountered `value` will be retained.
- `fill`: missing row/column combinations are filled with this value. The default
  is `missing`. If the `value` column is a `CategoricalVector` and `fill`
  is not `missing` then in order to keep unstacked value columns also
  `CategoricalVector` the `fill` must be passed as `CategoricalValue`

# Examples

```jldoctest
julia> wide = DataFrame(id=1:6,
                        a=repeat(1:3, inner=2),
                        b=repeat(1.0:2.0, inner=3),
                        c=repeat(1.0:1.0, inner=6),
                        d=repeat(1.0:3.0, inner=2))
6×5 DataFrame
 Row │ id     a      b        c        d
     │ Int64  Int64  Float64  Float64  Float64
─────┼─────────────────────────────────────────
   1 │     1      1      1.0      1.0      1.0
   2 │     2      1      1.0      1.0      1.0
   3 │     3      2      1.0      1.0      2.0
   4 │     4      2      2.0      1.0      2.0
   5 │     5      3      2.0      1.0      3.0
   6 │     6      3      2.0      1.0      3.0

julia> long = stack(wide)
18×4 DataFrame
 Row │ id     a      variable  value
     │ Int64  Int64  String    Float64
─────┼─────────────────────────────────
   1 │     1      1  b             1.0
   2 │     2      1  b             1.0
   3 │     3      2  b             1.0
   4 │     4      2  b             2.0
   5 │     5      3  b             2.0
   6 │     6      3  b             2.0
   7 │     1      1  c             1.0
   8 │     2      1  c             1.0
  ⋮  │   ⋮      ⋮       ⋮         ⋮
  12 │     6      3  c             1.0
  13 │     1      1  d             1.0
  14 │     2      1  d             1.0
  15 │     3      2  d             2.0
  16 │     4      2  d             2.0
  17 │     5      3  d             3.0
  18 │     6      3  d             3.0
                         3 rows omitted

julia> unstack(long)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :variable, :value)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :id, :variable, :value)
6×4 DataFrame
 Row │ id     b         c         d
     │ Int64  Float64?  Float64?  Float64?
─────┼─────────────────────────────────────
   1 │     1       1.0       1.0       1.0
   2 │     2       1.0       1.0       1.0
   3 │     3       1.0       1.0       2.0
   4 │     4       2.0       1.0       2.0
   5 │     5       2.0       1.0       3.0
   6 │     6       2.0       1.0       3.0

julia> unstack(long, [:id, :a], :variable, :value)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :id, :variable, :value, renamecols=x->Symbol(:_, x))
6×4 DataFrame
 Row │ id     _b        _c        _d
     │ Int64  Float64?  Float64?  Float64?
─────┼─────────────────────────────────────
   1 │     1       1.0       1.0       1.0
   2 │     2       1.0       1.0       1.0
   3 │     3       1.0       1.0       2.0
   4 │     4       2.0       1.0       2.0
   5 │     5       2.0       1.0       3.0
   6 │     6       2.0       1.0       3.0

julia> df = DataFrame(id=["1", "1", "2"],
                      variable=["Var1", "Var2", "Var1"],
                      value=[1, 2, 3])
3×3 DataFrame
 Row │ id      variable  value
     │ String  String    Int64
─────┼─────────────────────────
   1 │ 1       Var1          1
   2 │ 1       Var2          2
   3 │ 2       Var1          3

julia> unstack(df, :variable, :value, fill=0)
2×3 DataFrame
 Row │ id      Var1   Var2
     │ String  Int64  Int64
─────┼──────────────────────
   1 │ 1           1      2
   2 │ 2           3      0
```
Note that there are some differences between the widened results above.
"""
function unstack(df::AbstractDataFrame, rowkeys, colkey::ColumnIndex,
                 value::ColumnIndex; renamecols::Function=identity,
                 allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing)
    rowkey_ints = vcat(index(df)[rowkeys])
    @assert rowkey_ints isa AbstractVector{Int}
    g_rowkey = groupby(df, rowkey_ints)
    g_colkey = groupby(df, colkey)
    valuecol = df[!, value]
    return _unstack(df, rowkey_ints, index(df)[colkey], g_colkey,
                    valuecol, g_rowkey, renamecols, allowmissing, allowduplicates, fill)
end

function unstack(df::AbstractDataFrame, colkey::ColumnIndex, value::ColumnIndex;
                 renamecols::Function=identity,
                 allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing)
    colkey_int = index(df)[colkey]
    value_int = index(df)[value]
    return unstack(df, Not(colkey_int, value_int), colkey_int, value_int,
            renamecols=renamecols, allowmissing=allowmissing,
            allowduplicates=allowduplicates, fill=fill)
end

unstack(df::AbstractDataFrame; renamecols::Function=identity,
        allowmissing::Bool=false, allowduplicates::Bool=false, fill=missing) =
    unstack(df, :variable, :value, renamecols=renamecols, allowmissing=allowmissing,
            allowduplicates=allowduplicates, fill=fill)

# we take into account the fact that idx, starts and ends are computed lazily
# so we rather directly reference the gdf.groups
# this function is tailor made for unstack so it does assume that no groups were
# dropped (i.e. gdf.groups does not contain 0 entries)
function find_group_row(gdf::GroupedDataFrame)
    rows = zeros(Int, length(gdf))
    isempty(rows) && return rows

    filled = 0
    i = 1
    groups = gdf.groups
    while filled < length(gdf)
        group = groups[i]
        if rows[group] == 0
            rows[group] = i
            filled += 1
        end
        i += 1
    end
    return rows # return row index of first occurrence of each group in gdf.groups
end

function _unstack(df::AbstractDataFrame, rowkeys::AbstractVector{Int},
                  colkey::Int, g_colkey::GroupedDataFrame,
                  valuecol::AbstractVector, g_rowkey::GroupedDataFrame,
                  renamecols::Function,
                  allowmissing::Bool, allowduplicates::Bool, fill)
    rowref = g_rowkey.groups
    row_group_row_idxs = find_group_row(g_rowkey)
    Nrow = length(g_rowkey)

    @assert groupcols(g_colkey) == _names(df)[colkey:colkey]
    colref = g_colkey.groups
    Ncol = length(g_colkey)
    col_group_row_idxs = find_group_row(g_colkey)
    colref_map = df[col_group_row_idxs, colkey]
    if any(ismissing, colref_map) && !allowmissing
        throw(ArgumentError("Missing value in variable :$(_names(df)[colkey]). " *
                            "Pass `allowmissing=true` to skip missings."))
    end
    unstacked_val = [fill!(similar(valuecol,
                                   promote_type(eltype(valuecol), typeof(fill)),
                                   Nrow),
                           fill) for _ in 1:Ncol]

    mask_filled = falses(Nrow, Ncol)

    @assert length(rowref) == length(colref) == length(valuecol)
    for (k, (row_id, col_id, val)) in enumerate(zip(rowref, colref, valuecol))
        if !allowduplicates && mask_filled[row_id, col_id]
            throw(ArgumentError("Duplicate entries in unstack at row $k for key "*
                                "$(tuple((df[k, s] for s in rowkeys)...)) and variable $(colref_map[col_id]). " *
                                "Pass allowduplicates=true to allow them."))
        end
        unstacked_val[col_id][row_id] = val
        mask_filled[row_id, col_id] = true
    end

    # note that Symbol(renamecols(x)) must produce unique column names
    # and names between df1 and df2 must be unique
    df1 = df[row_group_row_idxs, g_rowkey.cols]
    df2 = DataFrame(unstacked_val, Symbol[Symbol(renamecols(x)) for x in colref_map],
                    copycols=false)

    @assert length(col_group_row_idxs) == ncol(df2)
    # avoid reordering when col_group_row_idxs was already ordered
    if !issorted(col_group_row_idxs)
        df2 = df2[!, sortperm(col_group_row_idxs)]
    end

    res_df = hcat(df1, df2, copycols=false)

    @assert length(row_group_row_idxs) == nrow(res_df)
    # avoid reordering when col_group_row_idxs was already ordered
    if !issorted(row_group_row_idxs)
        res_df = res_df[sortperm(row_group_row_idxs), :]
    end

    return res_df
end

"""
    StackedVector <: AbstractVector

An `AbstractVector` that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

# Constructor
```julia
StackedVector(d::AbstractVector)
```

# Arguments
- `d...` : one or more AbstractVectors

# Examples
```julia
StackedVector(Any[[1, 2], [9, 10], [11, 12]])  # [1, 2, 9, 10, 11, 12]
```
"""
struct StackedVector{T} <: AbstractVector{T}
    components::Vector{Any}
end

StackedVector(d::AbstractVector) =
    StackedVector{promote_type(map(eltype, d)...)}(d)

function Base.getindex(v::StackedVector{T}, i::Int)::T where T
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    return v.components[j][k]
end

Base.IndexStyle(::Type{StackedVector}) = Base.IndexLinear()
Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.eltype(v::Type{StackedVector{T}}) where {T} = T
Base.similar(v::StackedVector, T::Type, dims::Union{Integer, AbstractUnitRange}...) =
    similar(v.components[1], T, dims...)

"""
    RepeatedVector{T} <: AbstractVector{T}

An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

# Constructor
```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

# Arguments
- `parent` : the AbstractVector that's repeated
- `inner` : the numer of times each element is repeated
- `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

# Examples
```julia
RepeatedVector([1, 2], 3, 1)   # [1, 1, 1, 2, 2, 2]
RepeatedVector([1, 2], 1, 3)   # [1, 2, 1, 2, 1, 2]
RepeatedVector([1, 2], 2, 2)   # [1, 2, 1, 2, 1, 2, 1, 2]
```
"""
struct RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

Base.parent(v::RepeatedVector) = v.parent

function Base.getindex(v::RepeatedVector, i::Int)
    N = length(parent(v))
    idx = Base.fld1(mod1(i, v.inner*N), v.inner)
    parent(v)[idx]
end

Base.IndexStyle(::Type{<:RepeatedVector}) = Base.IndexLinear()
Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(parent(v))
Base.eltype(v::Type{RepeatedVector{T}}) where {T} = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(parent(v)), v.inner, v.outer)
Base.similar(v::RepeatedVector, T::Type, dims::Dims) = similar(parent(v), T, dims)
Base.unique(v::RepeatedVector) = unique(parent(v))

Base.transpose(::AbstractDataFrame, args...; kwargs...) =
    throw(ArgumentError("`transpose` not defined for `AbstractDataFrame`s. Try `permutedims` instead"))

"""
    permutedims(df::AbstractDataFrame, src_namescol::Union{Int, Symbol, AbstractString},
                [dest_namescol::Union{Symbol, AbstractString}];
                makeunique::Bool=false)

Turn `df` on its side such that rows become columns
and values in the column indexed by `src_namescol` become the names of new columns.
In the resulting `DataFrame`, column names of `df` will become the first column
with name specified by `dest_namescol`.

# Arguments
- `df` : the `AbstractDataFrame`
- `src_namescol` : the column that will become the new header.
  This column's element type must be `AbstractString` or `Symbol`.
- `dest_namescol` : the name of the first column in the returned `DataFrame`.
  Defaults to the same name as `src_namescol`.
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found; if `true`, duplicate names will be suffixed
  with `_i` (`i` starting at 1 for the first duplicate).

Note: The element types of columns in resulting `DataFrame`
(other than the first column, which always has element type `String`)
will depend on the element types of _all_ input columns
based on the result of `promote_type`.
That is, if the source data frame contains `Int` and `Float64` columns,
resulting columns will have element type `Float64`. If the source has
`Int` and `String` columns, resulting columns will have element type `Any`.

# Examples

```jldoctest
julia> df1 = DataFrame(a=["x", "y"], b=[1.0, 2.0], c=[3, 4], d=[true, false])
2×4 DataFrame
 Row │ a       b        c      d
     │ String  Float64  Int64  Bool
─────┼───────────────────────────────
   1 │ x           1.0      3   true
   2 │ y           2.0      4  false

julia> permutedims(df1, 1) # note the column types
3×3 DataFrame
 Row │ a       x        y
     │ String  Float64  Float64
─────┼──────────────────────────
   1 │ b           1.0      2.0
   2 │ c           3.0      4.0
   3 │ d           1.0      0.0

julia> df2 = DataFrame(a=["x", "y"], b=[1, "two"], c=[3, 4], d=[true, false])
2×4 DataFrame
 Row │ a       b    c      d
     │ String  Any  Int64  Bool
─────┼───────────────────────────
   1 │ x       1        3   true
   2 │ y       two      4  false

julia> permutedims(df2, 1, "different_name")
3×3 DataFrame
 Row │ different_name  x     y
     │ String          Any   Any
─────┼─────────────────────────────
   1 │ b               1     two
   2 │ c               3     4
   3 │ d               true  false
```
"""
function Base.permutedims(df::AbstractDataFrame, src_namescol::ColumnIndex,
                          dest_namescol::Union{Symbol, AbstractString};
                          makeunique::Bool=false)

    if src_namescol isa Integer
        1 <= src_namescol <= ncol(df) || throw(BoundsError(index(df), src_namescol))
    end
    eltype(df[!, src_namescol]) <: SymbolOrString ||
        throw(ArgumentError("src_namescol must have eltype `Symbol` or `<:AbstractString`"))

    df_notsrc = df[!, Not(src_namescol)]
    df_permuted = DataFrame(dest_namescol => names(df_notsrc))

    if ncol(df_notsrc) == 0
        df_tmp = DataFrame(AbstractVector[[] for _ in 1:nrow(df)], df[!, src_namescol],
                           makeunique=makeunique, copycols=false)
    else
        m = permutedims(Matrix(df_notsrc))
        df_tmp = rename!(DataFrame(Tables.table(m)), df[!, src_namescol], makeunique=makeunique)
    end
    return hcat!(df_permuted, df_tmp, makeunique=makeunique, copycols=false)
end

function Base.permutedims(df::AbstractDataFrame, src_namescol::ColumnIndex;
                          makeunique::Bool=false)
    if src_namescol isa Integer
        1 <= src_namescol <= ncol(df) || throw(BoundsError(index(df), src_namescol))
        dest_namescol = _names(df)[src_namescol]
    else
        dest_namescol = src_namescol
    end
    return permutedims(df, src_namescol, dest_namescol; makeunique=makeunique)
end


text:
import Gaugefields:staggered_U

struct StaggeredFermion_4D_wing{NC} <: AbstractFermionfields_4D{NC}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NDW::Int64
    NG::Int64 #size of the Gamma matrix. In Staggered fermion, this is one.
    NV::Int64
    f::Array{ComplexF64,6}
    Dirac_operator::String

    function StaggeredFermion_4D_wing(NC,NX,NY,NZ,NT)
        NG = 1
        NDW = 1
        NV = NC*NX*NY*NZ*NT*NG
        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW,NG)
        Dirac_operator = "Staggered"
        return new{NC}(NC,NX,NY,NZ,NT,NDW,NG,NV,f,Dirac_operator)
    end
end

function Base.size(x::StaggeredFermion_4D_wing{NC}) where NC
    return (x.NC,x.NX,x.NY,x.NZ,x.NT,x.NG)
    #return (x.NV,)
end

function Base.length(x::StaggeredFermion_4D_wing{NC}) where {NC}
    return NC*x.NX*x.NY*x.NZ*x.NT*x.NG
end

function Base.similar(x::T) where T <: StaggeredFermion_4D_wing
    return StaggeredFermion_4D_wing(x.NC,x.NX,x.NY,x.NZ,x.NT)
end

function Dx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1},boundarycondition) where  {T <: StaggeredFermion_4D_wing,G <:AbstractGaugefields}
    #temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    #clear!(temp)
    set_wing_fermion!(x,boundarycondition)
    clear_fermion!(xout)
    for ν=1:4
        xplus = shift_fermion(x,ν)
        Us = staggered_U(U[ν],ν)
        mul!(temp1,Us,xplus)

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)
        Uminus_s = staggered_U(Uminus,ν)
        mul!(temp2,Uminus_s',xminus)

        add_fermion!(xout,0.5,temp1,-0.5,temp2)

        #fermion_shift!(temp1,U,ν,x)
        #fermion_shift!(temp2,U,-ν,x)
        #add!(xout,0.5,temp1,-0.5,temp2)

    end

    set_wing_fermion!(xout,boundarycondition)

    return
end

function clear_fermion!(x::StaggeredFermion_4D_wing{NC},evensite) where NC
    ibush = ifelse(evensite,0,1)
    for it=1:x.NT
        for iz=1:x.NZ
            for iy=1:x.NY
                xran =1+(1+ibush+iy+iz+it)%2:2:x.NX
                for ix in xran
                    @simd for ic=1:NC
                        x[ic,ix,iy,iz,it,1] = 0
                    end
                end
            end
        end
    end
    return
end

function set_wing_fermion!(a::StaggeredFermion_4D_wing{NC},boundarycondition) where NC
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX

    #!  X-direction
    for ialpha=1:1
        for it=1:NT
            for iz = 1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,0,iy,iz,it,ialpha] = boundarycondition[1]*a[k,NX,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:1
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,NX+1,iy,iz,it,ialpha] = boundarycondition[1]*a[k,1,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    #Y-direction
    for ialpha = 1:1
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,0,iz,it,ialpha] = boundarycondition[2]*a[k,ix,NY,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:1
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,NY+1,iz,it,ialpha] = boundarycondition[2]*a[k,ix,1,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:1
        # Z-direction
        for it=1:NT
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,0,it,ialpha] = boundarycondition[3]*a[k,ix,iy,NZ,it,ialpha]
                        a[k,ix,iy,NZ+1,it,ialpha] = boundarycondition[3]*a[k,ix,iy,1,it,ialpha]

                    end
                end
            end
        end

        #T-direction
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,iz,0,ialpha] = boundarycondition[4]*a[k,ix,iy,iz,NT,ialpha]
                        a[k,ix,iy,iz,NT+1,ialpha] = boundarycondition[4]*a[k,ix,iy,iz,1,ialpha]
                    end
                end
            end
        end

    end

end


text:
type FDMMesh{X,MIN,MAX}
  Δxs::X
  mins::MIN
  maxs::MAX
  grids::Tuple
  dims::Int
  gridSize::Tuple
  square::Bool
end

function FDMMesh(Δxs::AbstractArray;mins=[0;0],maxs=[1;1],buildMesh=true)
  if length(mins)!=length(maxs) || length(Δxs)!=length(mins) error("DimensionMismatch") end
  dims = length(mins)
  if buildMesh && dims == 2
    grids = meshgrid(mins[1]:Δxs[1]:maxs[1],mins[2]:Δxs[2]:maxs[2])
    gridSize = size(grids[1])
  elseif buildMesh && dims == 3
    grids = meshgrid(mins[1]:Δxs[1]:maxs[1],mins[2]:Δxs[2]:maxs[2],mins[3]:Δxs[3]:maxs[3])
    gridSize = size(grids[1])
  else
    grids = nothing
    gridSize = (maxs-mins)./Δxs
  end
  square = minimum(map((x)->(x == gridSize[1]),gridSize))
  FDMMesh(Δxs,mins,maxs,grids,dims,gridSize,square)
end
FDMMesh(Δx::Number;mins=[0;0],maxs=[1;1],buildMesh=true) = FDMMesh(Δx*ones(mins),mins=mins,maxs=maxs,buildMesh=buildMesh)

"""
`size(mesh::FDMMesh)`

Returns `gridSize`.
"""
Base.size(mesh::FDMMesh) = mesh.gridSize

type StokesProblem{F1,F2,G,Ug,Vg,UA,VA,PA}
  f₁::F1
  f₂::F2
  g::G
  ugD::Ug
  vgD::Vg
  uanalytic::UA
  vanalytic::VA
  panalytic::PA
  trueknown::Bool
end

StokesProblem(f₁,f₂,g,uanalytic,vanalytic,panalytic) = StokesProblem(f₁,f₂,g,uanalytic,vanalytic,uanalytic,vanalytic,panalytic,true)
StokesProblem(f₁,f₂,g,ugD,vgD) = StokesProblem(f₁,f₂,g,ugD,vgD,nothing,nothing,nothing,false)


text:
module DailyTreasuryYieldCurve

using Reexport
using EzXML
import HTTP
@reexport using Dates
using Missings
@reexport using DataFrames
using Interpolations

export getyieldcurves
export RateInterpolator, createRateInterpolator

const DATAFEED = "https://data.treasury.gov/feed.svc/DailyTreasuryYieldCurveRateData"
const DATAFEEDREAL = "http://data.treasury.gov/feed.svc/DailyTreasuryRealYieldCurveRateData"

# these fields are provided by the feed (as of 20200413)
# note that we must use :m1 because :1m is not a valid symbol
const COLNAMES = [:id,:date,:m1,:m2,:m3,:m6,:y1,:y2,:y3,:y5,:y7,:y10,:y20,:y30,:y30dup]
const COLNAMESREAL = [:id,:date,:y5,:y7,:y10,:y20,:y30]

"""
    getyieldcurves(;real::Bool=false,begdt::Date=Date(1990,1,2),enddt::Date=today())

Download the whole published history of daily US Treasury yield curves from the official data feed.

Optionally, pass a filename if you have already downloaded the data, eg. `getyieldcurves(fn::AbstractString)`.

By default, gets the nominal yield curve. Pass `realrates=true` to get the real yield curve.  The nominal series starts in 1990, while the real series starts in 2003.

Returns a `DataFrame`.
"""
function getyieldcurves(;realrates::Bool=false,begdt::Date=Date(1990,1,2),enddt::Date=today())
    startyr=max(year(begdt),1990)-1
    endyr=min(year(enddt),year(today()))+1
    datefilter = string("?\$filter=","year(NEW_DATE)%20gt%20",startyr,"%20and%20","year(NEW_DATE)%20lt%20",endyr)
    if realrates
        h = HTTP.get(string(DATAFEEDREAL,datefilter))
    else
        h = HTTP.get(string(DATAFEED,datefilter))
    end
    thexml = parsexml(String(h.body))
    df =_parseyieldcurves(thexml,realrates)
    filter!(row -> begdt <= row.date <= enddt,df)
    return df
end

function getyieldcurves(fn::AbstractString;realrates::Bool=false,begdt::Date=Date(1990,1,2),enddt::Date=today())
    thexml = readxml(fn)
    df =_parseyieldcurves(thexml,realrates)
    filter!(row -> begdt <= row.date <= enddt,df)
    return df
end

"""
    _parseyieldcurves(thexml,realrates)

Parser function for Treasury yield curve data. (unexported)
"""
function _parseyieldcurves(thexml,realrates)
    if realrates
        return _parserealcurves(thexml)
    else
        return _parsenominalcurves(thexml)
    end
end

function _parsenominalcurves(thexml)
    r = thexml.root
    cnod=findall("/*/*/*[position()=7]",r)
    @assert length(cnod) >= 1 "Rate XML parse error likely."

    stripsplitstrip(s) = strip.(split(strip(s),"\n"))

    strarr = map(stripsplitstrip,nodecontent.(cnod))
    df=DataFrame(Tuple.(strarr))
    @assert nrow(df) >= 1 "Rate XML parse error likely."

    rename!(df,COLNAMES)
    select!(df,Not([:id; :y30dup]))
    for c in COLNAMES[2:end-1]
        if c == :date
            df[!,c] = Date.(parse.(DateTime,df[!,c]))
        else
            df[!,c] = passmissing(x->parse(Float64,x)).(replace(df[!,c],""=>missing))
        end
    end
    sort!(df,:date)
    return df
end

function _parserealcurves(thexml)
    r = thexml.root
    cnod=findall("/*/*/*[position()=7]",r)
    @assert length(cnod) >= 1 "Rate XML parse error likely."

    nodestr = nodecontent.(cnod)

    function splitandfill(s)
        spl = split(s)
        L = length(spl)
        L == length(COLNAMESREAL) && return spl

        strvec = Vector{eltype(spl)}(undef,length(COLNAMESREAL))
        for i in 1:length(COLNAMESREAL)
            if i <= L
                strvec[i] = spl[i]
            else
                strvec[i] = ""
            end
        end
        return strvec
    end

    strarr = splitandfill.(nodestr)
    df=DataFrame(Tuple.(strarr))
    @assert nrow(df) >= 1 "Rate XML parse error likely."

    rename!(df,COLNAMESREAL)
    select!(df,Not(:id))
    for c in COLNAMESREAL[2:end]
        if c == :date
            df[!,c] = Date.(parse.(DateTime,df[!,c]))
        else
            df[!,c] = passmissing(x->parse(Float64,x)).(replace(df[!,c],""=>missing))
        end
    end
    sort!(df,:date)
    return df
end

# include additional files
include("interp.jl")

end # module


text:
using PyPlot
using Optim

function get_ln_delta(N, d_VC, epsilon)
    return log(4) + d_VC * log(2N) -(epsilon^2)*N/8
end;

d_VC = 10;
epsilon = 0.05;
N_arr = [400000. 420000. 440000. 460000. 480000.];
for N in N_arr
    delta = exp(get_ln_delta(N, d_VC, epsilon));
    println("N = $N, delta = $delta")
end

function get_error(N, d_VC, epsilon, target_delta)
    delta = exp(get_ln_delta(N, d_VC, epsilon));
    return (delta - target_delta)^2;
end;

result = optimize(N -> get_error(N, 10, 0.05, 0.05), 400000, 500000)

N = result.minimum;
err = result.f_minimum;
if result.converged == true
    println("Solution:  N = $N  (error = $err)")
end

function get_vc_bound(N,d_VC,delta)
    inner = log(4) + d_VC * log(2N) - log(delta);
    return sqrt(8inner/N);
end;

function get_rademacher_penalty(N,d_VC,delta)
    inner = 2(log(2N) + d_VC * log(N))/N
    return sqrt(inner) + sqrt(-2log(delta)/N) + 1/N;
end;

function get_parrondo_and_van_den_broek(N,d_VC,delta)
    pvdb_err(epsilon,N,d_VC,delta) = (epsilon - sqrt((2epsilon + log(6) + d_VC*log(2N) - log(delta))/N))^2;
    result = optimize(epsilon -> pvdb_err(epsilon,N,d_VC,delta), 0.0, 10.0);
    if result.converged==true
        return result.minimum
    else
        return NaN
    end;
end;

function get_devroye(N,d_VC,delta)
    devroye_err(epsilon,N,d_VC,delta) = (epsilon - sqrt((4epsilon*(1+epsilon) + log(4) + d_VC*log(N^2) - log(delta))/(2N)))^2;
    result = optimize(epsilon -> devroye_err(epsilon,N,d_VC,delta), 0.0, 20.0);
    if result.converged==true
        return result.minimum
    else
        return NaN
    end;
end;

d_VC = 50;
delta = 0.05;

#N_arr = logspace(0,4,11);  #from 1 to 10k
N_arr = [5, 10, 50, 100, 500, 1000, 5000, 10000];
epsilon_vc = [ get_vc_bound(N,d_VC,delta) for N in N_arr ];
epsilon_rp = [ get_rademacher_penalty(N,d_VC,delta) for N in N_arr ];
epsilon_pv = [ get_parrondo_and_van_den_broek(N,d_VC,delta) for N in N_arr ];
epsilon_dv = [ get_devroye(N,d_VC,delta) for N in N_arr ];

function plot_epsilons(N_arr,epsilon_arr,label_arr)
    fig = figure(figsize=(7,4))
    ax = fig[:add_subplot](1,1,1)
    for i in 1:size(label_arr)[1]
        #ax[:loglog](N_arr,epsilon_arr[:,i],label=label_arr[i]);
        loglog(N_arr,epsilon_arr[:,i],"x-",label=label_arr[i]);
    end;
    #ax[:legend](loc="best",frameon=false);
    legend(loc="best",frameon=false,fontsize=10)
    grid(true);
end;

epsilon_arr = [ epsilon_vc epsilon_rp epsilon_pv epsilon_dv ];
label_arr = [ "VC", "Rademacher", "Parrondo & Van den Broek", "Devroye" ];
plot_epsilons(N_arr,epsilon_arr,label_arr);

println("N = $(N_arr[end])")
println("ϵ (VC)                       = $(epsilon_vc[end])")
println("ϵ (Rademacher)               = $(epsilon_rp[end])")
println("ϵ (Parrondo & Van den Broek) = $(epsilon_pv[end])")
println("ϵ (Devroye)                  = $(epsilon_dv[end])")

println("N = $(N_arr[1])")
println("ϵ (VC)                       = $(epsilon_vc[1])")
println("ϵ (Rademacher)               = $(epsilon_rp[1])")
println("ϵ (Parrondo & Van den Broek) = $(epsilon_pv[1])")
println("ϵ (Devroye)                  = $(epsilon_dv[1])")

function generate_dataset(;n_points=200)
    x = linspace(-1,1,n_points);
    y = sin(π*x);
    x_sample = rand(2);
    y_sample = sin(π*x_sample);
    return x,y,x_sample,y_sample;
end;

function fit_a(x_sample,y_sample)
    sum = x_sample[1]*y_sample[1]+x_sample[2]*y_sample[2]
    return sum/(x_sample[1]^2+x_sample[2]^2);
end;

function predict_y_a(x,a)
    return a .* x;
end;

function plot_xy(x,y,x_sample,y_sample,y_pred)
    figure(figsize=(7,5));
    plot(x,y,"m-",label=L"$\sin(\pi x)$");
    plot(x,y_pred,"g--",label="predicted");
    plot(x_sample,y_sample,"bo",label="sample data");
    grid(true);
    legend(loc="best",frameon=false);
end;
x,y,x_sample,y_sample = generate_dataset();
y_pred = predict_y_a(x,fit_a(x_sample,y_sample));
plot_xy(x,y,x_sample,y_sample,y_pred);

function get_avg_a(n_trials; n_points=200)
    a_arr = Float64[];
    for i in 1:n_trials
        x_sample = rand(2);
        y_sample = sin(π*x_sample);
        a = fit_a(x_sample,y_sample);
        push!(a_arr,a);
    end
    return mean(a_arr);
end;
a = get_avg_a(1e6);    #1.423, 1.425, 1.424, 1.427
println("<g(x)> = $(round(a,2))x");

function get_bias_variance(avg_a, n_trials; n_points=200)
    bias  = Float64[];
    var   = Float64[];
    for i in 1:n_trials
        x_sample = rand(2);
        y_sample = sin(π*x_sample);
        a = fit_a(x_sample,y_sample);
        g_bar_x = avg_a*x_sample;
        push!(bias,mean((g_bar_x - y_sample).^2));
        push!(var,mean((g_bar_x - a .* x_sample).^2));
    end
    return bias, var;
end;

bias,var = get_bias_variance(a, 1e7);
println("<bias> = $(mean(bias))");
println("<var>  = $(mean(var))");

max_stars_repo_name: nathanielng/machine-learning

max_issues_repo_name: nathanielng/machine-learning

max_forks_repo_name: nathanielng/machine-learning


text:
using CSV, JLD2, DataStructures, Statistics, Plots
include("src/SchoolOutbreak.jl")
@load "data/anonymizedstudents.jld2"

function temporalpatterns(students)
    strata=(student->getfield.(Ref(student),(:schoolID,:gradeID,:classID)).|>first).(filter(x->x.isinfected,students))
    onsets=getfield.(students,:onset)
    stindex=unique(strata)|>sort
    OrderedDict(stindex.=>[onsets[findall(x->st==x,strata)]|>sort for st in stindex])
end
@time onsetbystratum=temporalpatterns(anonymizedstudents)

gr(fontfamily="Helvetica",foreground_color_legend = nothing,background_color_legend = nothing, titlefontsize=9,legendfontsize=9,guidefontsize=10,grid=false, tick_direction=:out)
function onsetbylevel(onsetbystratum, level,origin=mean)
    ukeys=unique(getindex.(keys(onsetbystratum),Ref(1:level)))
    vectors=(vcat((filter(x->x[1][1:level]==ukey,onsetbystratum)|>values)...) for ukey in ukeys)
    filter(x->x!=0,vcat(((v->v.-origin(v)).(vectors))...))
end
devonsets=onsetbylevel.(Ref(onsetbystratum),0:3)
plot(((v,c)->histogram(v,xlim=(-100,100),color=c,normalize=true)).(devonsets,[6,9,:skyblue,1])...,size=(550,550),layout=(4,1),titles=["all" "within-school" "within-grade" "within-class"],xlabel=[fill("",1,3) "deviation of onset dates from mean"],labels="SD = ".*string.(round.(std.(devonsets)',digits=1)))

## Clustering simulation
using JLD2
@load "Matsumoto_studentsdata.jld2" students
duringbreak=(d2di(Date("2014-12-27")):d2di(Date("2015-1-7"))).+1 # 27 Dec - 7 Jan was the final day of school in Matsumoto, shifted one day to account for presymptomatic infection
@time simstudents=SchoolOutbreak.Students.([filter(x->x.stratum[1]==schoolid,students) for schoolid in (1:29)[Not([3,5,21])]],Ref(parameters),[Dict(duringbreak=>duringbreak#=,duringclosure=duringclosure[x]=#) for x in (1:29)[Not([3,5,21])]]);
chain1=CSV.read("output/Matsumoto_MCMCsamples.csv",DataFrame)

function parametersfromchain(chain, iter)
    (β=SchoolOutbreak.posterior(chain,:βs,iter).samples[:,[1;1:4]]|>vec,
    γ=SchoolOutbreak.posterior(chain,:γs,iter).samples[:,[1;1:3;3]]|>vec,
    δ=SchoolOutbreak.posterior(chain,:δs,iter).samples[:,[1;1:3;3]]|>vec,
    rcom=SchoolOutbreak.posterior(chain,:rcom,iter).samples|>first|>fill,
    suscoef=SchoolOutbreak.posterior(chain,:suscoef,iter).samples|>vec,
    infcoef=SchoolOutbreak.posterior(chain,:infcoef,iter).samples|>vec)
end
function parametersfromchain(chain::DataFrame, iter)
    chainmat=Matrix(chain)
    (β=(chainmat[iter,[2;2:5]]|>vec) .|>exp,
    γ=chainmat[iter,[6;6:8;8]]|>vec,
    δ=chainmat[iter,[9;9:11;11]]|>vec,
    rcom=((chainmat[iter,1])|>exp).*1,
    suscoef=chainmat[iter,12:15]|>vec,
    infcoef=chainmat[iter,16:20]|>vec)
end
@time updateparameters=parametersfromchain.(Ref(chain1),(1:50).*20);
Random.seed!(2021)
@time outbreaks=SchoolOutbreak.simulateoutbreaks!.(Ref(simstudents), updateparameters, Ref(1:250),initcases=0);
incidences=SchoolOutbreak.incidence(outbreaks);

# plot group-level deviation of onset dates
function onset(outbreak,students)
    onsetpairs=((ob,st)->reduce(vcat,(Tuple.(getindex.(getfield.(st.members[ids],:stratum),Ref(1:3))).=>date) for (date, ids) in enumerate(ob))).(outbreak,students) # Pairs: stratum => onsetdate
    onsetpairs=vcat(onsetpairs...)
    ukeys=unique(first.(onsetpairs))|>sort
    Dict(ukeys.=>[filter(x->first(x)==ukey,onsetpairs).|>last for ukey in ukeys])
end
onsets=[onset(outbreaks[n],simstudents) for n in 1:50]
gr(fontfamily="Helvetica",foreground_color_legend = nothing,background_color_legend = nothing, titlefontsize=9,legendfontsize=9,guidefontsize=10,grid=false, tick_direction=:out)
#=joinbylevel(v::Array{<:Array},level)=level==0 ? [joinbylevel(vcat(v...),0)] : joinbylevel.(v,level-1)
joinbylevel(v::Array{<:Number},level)=v
function onsetbylevel(simonsets, level,origin=mean)
    vectors=joinbylevel(simonsets,level)
    vcat(((v->v.-origin(v)).(vectors))...)
    end=#
function onsetbylevel(onsetbystratum, level,origin=mean)
    ukeys=unique(getindex.(keys(onsetbystratum),Ref(1:level)))
    vectors=[vcat((filter(x->x[1][1:level]==ukey,onsetbystratum)|>values)...) for ukey in ukeys]
    vcat(((v->v.-origin(v)).(filter(x->length(x)>1,vectors)))...)
end
devonsets=vcat.((onsetbylevel.(Ref(onsets[n]),0:3) for n in (1:50).*1)...)
Plots.plot(((v,c)->histogram(v,xlim=(-100,100),color=c,normalize=true)).(devonsets,[6,9,:skyblue,1])...,size=(550,550),layout=(4,1),titles=["all" "within-school" "within-grade" "within-class"],xlabel=[fill("",1,3) "deviation of onset dates from mean"],labels="SD = ".*string.(round.(std.(devonsets)',digits=1)))


text:
# Copyright (c) 2021, Rolfe Power

module Pat

using DocStringExtensions

# Constant Definitions
include("constants.jl")

# Dimensions
include("dimensions.jl")

# Celestial Bodies
include("bodies.jl")

# Two Body
include("twobody/twobody.jl")

end


text:
module RegService

using Sockets
using Distributed
using Pkg
using Logging

import ..RegEdit: register
using ..Messaging

const config = Dict{String,Any}()

include("management.jl")

"""
    service(zsock::ReplySocket)
Start a server that calls `RegEdit.register` on incoming
registration requests from `zsock`. Registration requests
must be serialized `RegEdit.RegisterParams` objects. Close
the socket to halt the server.

Parameters:
- `zsock::ReplySocket`

Returns:
nothing
"""
function service(zsock::ReplySocket)
    while true
        ret = recvsend(zsock) do regp
            if isempty(regp.gitconfig)
                haskey(config, "user") && (regp.gitconfig["user.name"] = config["user"])
                haskey(config, "email") && (regp.gitconfig["user.email"] = config["email"])
            end
            register(regp)
        end
        ret || return
    end
    nothing
end

function main()
    if isempty(ARGS)
        println("Usage: julia -e 'using Registrator; Registrator.RegService.main()' <configuration>")
        return
    end

    merge!(config, Pkg.TOML.parsefile(ARGS[1])["regservice"])
    global_logger(SimpleLogger(stdout, get_log_level(config["log_level"])))

    zsock = ReplySocket(get(config, "port", 5555))

    @info("Starting registration service...")
    t = @async status_monitor(config["stop_file"], zsock)
    service(zsock)
    wait(t)

    # The !stopped! part is grep'ed by restart.sh, do not change
    @warn("Registration service !stopped!")
end

end    # module


text:
struct PrimHeapEntry{T<:Real}
    edge::Edge
    dist::T
end

isless(e1::PrimHeapEntry, e2::PrimHeapEntry) = e1.dist < e2.dist

"""
    prim_mst(g, distmx=weights(g))

Return a vector of edges representing the minimum spanning tree of a connected, undirected graph `g` with optional
distance matrix `distmx` using [Prim's algorithm](https://en.wikipedia.org/wiki/Prim%27s_algorithm).
Return a vector of edges.
"""
function prim_mst end
@traitfn function prim_mst(
    g::::(!IsDirected),
    distmx::AbstractMatrix = weights(g)
    )
    pq = Vector{PrimHeapEntry}()
    mst = Vector{Edge}()
    marked = zeros(Bool, nv(g))

    sizehint!(pq, ne(g))
    sizehint!(mst, ne(g))
    visit!(g, 1, marked, pq, distmx)

    while !isempty(pq)
        heap_entry = heappop!(pq)
        v = src(heap_entry.edge)
        w = dst(heap_entry.edge)

        marked[v] && marked[w] && continue
        push!(mst, heap_entry.edge)

        !marked[v] && visit!(g, v, marked, pq, distmx)
        !marked[w] && visit!(g, w, marked, pq, distmx)
    end

    return mst
end

"""
    visit!(g, v, marked, pq, distmx)

Mark the vertex `v` of graph `g` true in the array `marked` and enter all its
edges into priority queue `pq` with its `distmx` values as a PrimHeapEntry.
"""
function visit!(
    g::AbstractGraph,
    v::Integer,
    marked::AbstractVector{Bool},
    pq::AbstractVector,
    distmx::AbstractMatrix
)
    marked[v] = true
    for w in out_neighbors(g, v)
        if !marked[w]
            x = min(v, w)
            y = max(v, w)
            heappush!(pq, PrimHeapEntry(Edge(x, y), distmx[x, y]))
        end
    end
end


text:
using MyPackagenew
using Test

my_funct(2,3)

der(2,1)

@testset "MyPackagenew.jl" begin
@test my_funct(2,1) == 5
@test my_funct(2,2) == 6
@test my_funct(2,2) == 6
end


text:
using Base.Test
using SparseRF
srand(16)

n,m = 10^3, 5 ;
features = rand(n,m);
weights = rand(-1:1,m);
labels = _int(features * weights);

maxdepth = 3
model = build_tree(labels, features, 0, maxdepth)
@test depth(model) == maxdepth

println("\n##### nfoldCV Classification Tree #####")
accuracy = nfoldCV_tree(labels, features, 0.9, 3)
@test mean(accuracy) > 0.7

println("\n##### nfoldCV Classification Forest #####")
accuracy = nfoldCV_forest(labels, features, 2, 10, 3)
@test mean(accuracy) > 0.7

println("\n##### nfoldCV Adaboosted Stumps #####")
accuracy = nfoldCV_stumps(labels, features, 7, 3)
@test mean(accuracy) > 0.5


text:
include("lehmer.jl")
include("soliton.jl")
include("codec.jl")

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


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

name = "Qt5Svg"
version = v"5.15.3"

# Collection of sources required to build qt5
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtsvg-everywhere-opensource-src-$version.tar.xz",
                  "3adc41dfcc67bbe3b8ff553bdac30ee75e270745536a58e54cdb741fa0505d89"),
]

script = raw"""
cd $WORKSPACE/srcdir
mkdir build
cd build/

qtsrcdir=`ls -d ../qtsvg-*`

qmake $qtsrcdir
make -j${nproc}

# It appends the root path here with the qmake dir, which is /workspace/artifacts/blah...
make INSTALL_ROOT=$PWD/dummyroot install
rsync -ua dummyroot/workspace/*/artifacts/*/* $prefix

sed -i "s?^prefix=.*?prefix=$prefix?" $prefix/lib/pkgconfig/Qt5Svg.pc

install_license $qtsrcdir/LICENSE.LGPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms_linux = [
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "freebsd"),
    Platform("powerpc64le", "linux"; libc="glibc"),
]
platforms_linux = expand_cxxstring_abis(platforms_linux)
platforms_win = expand_cxxstring_abis([Platform("x86_64", "windows"), Platform("i686", "windows")])
platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt5Svg", "libQt5Svg", "QtSvg"], :libqt5svg),
]

products_macos = [
    FrameworkProduct("QtSvg", :libqt5svg),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Qt5Base_jll"),
]

include("../../fancy_toys.jl")

julia_compat = "1.6"

if any(should_build_platform.(triplet.(platforms_linux)))
    build_tarballs(ARGS, name, version, sources, script, platforms_linux, products, dependencies; preferred_gcc_version = v"7", julia_compat)
end
if any(should_build_platform.(triplet.(platforms_win)))
    build_tarballs(ARGS, name, version, sources, script, platforms_win, products, dependencies; preferred_gcc_version = v"8", julia_compat)
end
if any(should_build_platform.(triplet.(platforms_macos)))
    build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"7", julia_compat)
end

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
export iauTpors
"""
In the tangent plane projection, given the rectangular coordinates
of a star and its spherical coordinates, determine the spherical
coordinates of the tangent point.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   xi,eta     double  rectangular coordinates of star image (Note 2)
   a,b        double  star's spherical coordinates (Note 3)

Returned:
   *a01,*b01  double  tangent point's spherical coordinates, Soln. 1
   *a02,*b02  double  tangent point's spherical coordinates, Soln. 2

Returned (function value):
               int     number of solutions:
                     0 = no solutions returned (Note 5)
                     1 = only the first solution is useful (Note 6)
                     2 = both solutions are useful (Note 6)

Notes:

   1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

   2. The eta axis points due north in the adopted coordinate system.
   If the spherical coordinates are observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the spherical coordinates are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

   3. All angular arguments are in radians.

   4. The angles a01 and a02 are returned in the range 0-2pi.  The
   angles b01 and b02 are returned in the range +/-pi, but in the
   usual, non-pole-crossing, case, the range is +/-pi/2.

   5. Cases where there is no solution can arise only near the poles.
   For example, it is clearly impossible for a star at the pole
   itself to have a non-zero xi value, and hence it is meaningless
   to ask where the tangent point would have to be to bring about
   this combination of xi and dec.

   6. Also near the poles, cases can arise where there are two useful
   solutions.  The return value indicates whether the second of the
   two solutions returned is useful;  1 indicates only one useful
   solution, the usual case.

   7. The basis of the algorithm is to solve the spherical triangle PSC,
   where P is the north celestial pole, S is the star and C is the
   tangent point.  The spherical coordinates of the tangent point are
   [a0,b0];  writing rho^2 = (xi^2+eta^2) and r^2 = (1+rho^2), side c
   is then (pi/2-b), side p is sqrt(xi^2+eta^2) and side s (to be
   found) is (pi/2-b0).  Angle C is given by sin(C) = xi/rho and
   cos(C) = eta/rho.  Angle P (to be found) is the longitude
   difference between star and tangent point (a-a0).

   8. This function is a member of the following set:

      spherical      vector         solve for

      iauTpxes      iauTpxev         xi,eta
      iauTpsts      iauTpstv          star
      > iauTpors <    iauTporv         origin

Called:
   iauAnp       normalize angle into range 0 to 2pi

References:

   Calabretta M.R. & Greisen, E.W., 2002, "Representations of
   celestial coordinates in FITS", Astron.Astrophys. 395, 1077

   Green, R.M., "Spherical Astronomy", Cambridge University Press,
   1987, Chapter 13.

This revision:   2018 January 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# int iauTpors(double xi, double eta, double a, double b,
#              double *a01, double *b01, double *a02, double *b02)

function iauTpors(xi::Real, eta::Real, a::Real, b::Real)
   ref_a01 = Ref{Float64}(0.0)
   ref_b01 = Ref{Float64}(0.0)
   ref_a02 = Ref{Float64}(0.0)
   ref_b02 = Ref{Float64}(0.0)

   status = ccall((:iauTpors, libsofa_c), Cint,
       (Cdouble, Cdouble, Cdouble, Cdouble,
       Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
       convert(Float64, xi),
       convert(Float64, eta),
       convert(Float64, a),
       convert(Float64, b),
       ref_a01, ref_b01, ref_a02, ref_b02)

   return status, ref_a01[], ref_b01[], ref_a02[], ref_b02[]
end


text:
abstract type AbstractChargeCloud end

struct PointCharge{T <: SSDFloat, S} <: AbstractChargeCloud
    charge::T
    pos::AbstractCoordinatePoint{T, 3, S}
end

struct NBodyChargeCloud{T <: SSDFloat, S} <: AbstractChargeCloud
    # To be done...
end


text:
module AODictionary

"""
I initially forked this from https://gist.github.com/bzinberg/fb3981dd2251fb25cdbcd4f90b8b9a72
but that makes a gist, which doesn't really play nicely with Julia's packaging system
https://gist.github.com/lawless-m/4812e1aef4b4782302c2866761f2f3e0
So now it's here.
"""

export AODict

"""
Append-only dictionary that:

* Preserves insertion order
* Supports concise and efficient linear indexing into the sequence of `(key,
  value)` pairs
* Supports concise and efficient linear indexing into the sequence of values.

Example:

```jldoctest
>>> d = AODict(:a => 2, :b => 4, :c => 6);
>>> d[:b]
4
>>> d.seq[2]
:b => 4
>>> d.seqvals[2]
4
>>> d.seqvals[2:end]
2-element Array{Int64,1}:
 4
 6
```
"""
struct AODict{K, V} <: AbstractDict{K, V}
  dict::Dict{K, V}
  keys::Vector{K}
  index::Dict{K, Int}
end

function AODict{K, V}(kvs) where {K, V}
  ks = collect(map(first, kvs))
  dict = Dict{K, V}(kvs)
  idx = Dict([k=>i for (i,k) in enumerate(ks)])
  return AODict{K, V}(dict, ks, idx)
end

AODict{K, V}(kvs::Pair...) where {K, V} = AODict{K, V}(kvs)
AODict(kvs::Pair{K, V}...) where {K, V} = AODict{K, V}(kvs)
AODict(kvs) = AODict{Any, Any}(kvs)

function Base.iterate(d::AODict)
  s = Base.iterate(d.keys)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (k => d.dict[k], state)
end

function Base.iterate(d::AODict, state)
  s = Base.iterate(d.keys, state)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (k => d.dict[k], state)
end

Base.length(d::AODict) = Base.length(d.keys)

Base.getindex(d::AODict{K}, key::K) where K = d.dict[key]

function Base.setindex!(d::AODict{K, V}, val::V, key::K) where {K, V}
  @assert key ∉ keys(d.dict)
  d.dict[key] = val
  push!(d.keys, key)
  d.index[key] = length(d.keys)
  return val
end

# Custom overload of `getproperty` so that we can write `d.seq` and `d.seqvals`
function Base.getproperty(d::AODict, prop::Symbol)
  if prop == :seq
    return _SequentialAOOD(d)
  elseif prop == :seqvals
    return _SequentialAOODVals(d)
  end
  return Base.getfield(d, prop)
end

Base.propertynames(d::AODict) = [Base.fieldnames(d)..., :seq, :seqvals]

Base.in(v, d::AODict) = v in d.keys

"""
Read-only, array-like indexable accessor for the `(key, value)` pairs of an
`AODict`.
"""
struct _SequentialAOOD
  d::AODict
end

function Base.iterate(seq::_SequentialAOOD)
  s = Base.iterate(seq.d.keys)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (k => seq.d.dict[k], state)
end

function Base.iterate(seq::_SequentialAOOD, state)
  s = Base.iterate(seq.d.keys, state)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (k => seq.d.dict[k], state)
end

Base.length(seq::_SequentialAOOD) = Base.length(seq.d)

function Base.getindex(seq::_SequentialAOOD, i::Integer)
  k = seq.d.keys[i]
  return k => seq.d.dict[k]
end

function Base.getindex(seq::_SequentialAOOD, r::AbstractRange)
  ks = seq.d.keys[r]
  return [k => seq.d.dict[k] for k in ks]
end

Base.firstindex(seq::_SequentialAOOD) = Base.firstindex(seq.d.keys)
Base.lastindex(seq::_SequentialAOOD) = Base.lastindex(seq.d.keys)

"""
Read-only, array-like indexable accessor for the values of an
`AODict`.
"""
struct _SequentialAOODVals
  d::AODict
end

function Base.iterate(seq::_SequentialAOODVals)
  s = Base.iterate(seq.d.keys)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (seq.d.dict[k], state)
end

function Base.iterate(seq::_SequentialAOODVals, state)
  s = Base.iterate(seq.d.keys, state)
  if isnothing(s); return nothing; end
  (k, state) = s
  return (seq.d.dict[k], state)
end

Base.length(seq::_SequentialAOODVals) = Base.length(seq.d)

function Base.getindex(seq::_SequentialAOODVals, i::Integer)
  k = seq.d.keys[i]
  return seq.d.dict[k]
end

function Base.getindex(seq::_SequentialAOODVals, r::AbstractRange)
  ks = seq.d.keys[r]
  return [seq.d.dict[k] for k in ks]
end

Base.firstindex(seq::_SequentialAOODVals) = Base.firstindex(seq.d.keys)
Base.lastindex(seq::_SequentialAOODVals) = Base.lastindex(seq.d.keys)

###
end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
res = ŷ .- y
stem(res)

savefig("assets/literate/ISL-lab-3-res.svg") # hide


text:
using MesocosmMetabolism
using Test

@testset "MesocosmMetabolism.jl" begin
    @test 3 == 3
end

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
doc"""
    Beta(α,β)

The *Beta distribution* has probability density function

$f(x; \alpha, \beta) = \frac{1}{B(\alpha, \beta)}
 x^{\alpha - 1} (1 - x)^{\beta - 1}, \quad x \in [0, 1]$

The Beta distribution is related to the [`Gamma`](:func:`Gamma`) distribution via the
property that if $X \sim \operatorname{Gamma}(\alpha)$ and $Y \sim \operatorname{Gamma}
(\beta)$ independently, then $X / (X + Y) \sim \operatorname{Beta}(\alpha, \beta)$.

```julia
Beta()        # equivalent to Beta(1, 1)
Beta(a)       # equivalent to Beta(a, a)
Beta(a, b)    # Beta distribution with shape parameters a and b

params(d)     # Get the parameters, i.e. (a, b)
```

External links

* [Beta distribution on Wikipedia](http://en.wikipedia.org/wiki/Beta_distribution)

"""

immutable Beta{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    β::T

    function Beta(α::T, β::T)
        @check_args(Beta, α > zero(α) && β > zero(β))
        new(α, β)
    end
end

Beta{T<:Real}(α::T, β::T) = Beta{T}(α, β)
Beta(α::Real, β::Real) = Beta(promote(α, β)...)
Beta(α::Integer, β::Integer) = Beta(Float64(α), Float64(β))
Beta(α::Real) = Beta(α, α)
Beta() = Beta(1, 1)

@distr_support Beta 0.0 1.0

#### Conversions
function convert{T<:Real}(::Type{Beta{T}}, α::Real, β::Real)
    Beta(T(α), T(β))
end
function convert{T <: Real, S <: Real}(::Type{Beta{T}}, d::Beta{S})
    Beta(T(d.α), T(d.β))
end

#### Parameters

params(d::Beta) = (d.α, d.β)
@inline partype{T<:Real}(d::Beta{T}) = T

#### Statistics

mean(d::Beta) = ((α, β) = params(d); α / (α + β))

function mode(d::Beta)
    (α, β) = params(d)
    (α > 1 && β > 1) || error("mode is defined only when α > 1 and β > 1.")
    return (α - 1) / (α + β - 2)
end

modes(d::Beta) = [mode(d)]

function var(d::Beta)
    (α, β) = params(d)
    s = α + β
    return (α * β) / (abs2(s) * (s + 1))
end

meanlogx(d::Beta) = ((α, β) = params(d); digamma(α) - digamma(α + β))

varlogx(d::Beta) = ((α, β) = params(d); trigamma(α) - trigamma(α + β))
stdlogx(d::Beta) = sqrt(varlogx(d))

function skewness(d::Beta)
    (α, β) = params(d)
    if α == β
        return zero(α)
    else
        s = α + β
        (2(β - α) * sqrt(s + 1)) / ((s + 2) * sqrt(α * β))
    end
end

function kurtosis(d::Beta)
    α, β = params(d)
    s = α + β
    p = α * β
    6(abs2(α - β) * (s + 1) - p * (s + 2)) / (p * (s + 2) * (s + 3))
end

function entropy(d::Beta)
    α, β = params(d)
    s = α + β
    lbeta(α, β) - (α - 1) * digamma(α) - (β - 1) * digamma(β) +
        (s - 2) * digamma(s)
end

#### Evaluation

@_delegate_statsfuns Beta beta α β

gradlogpdf{T<:Real}(d::Beta{T}, x::Real) =
    ((α, β) = params(d); 0 <= x <= 1 ? (α - 1) / x - (β - 1) / (1 - x) : zero(T))

#### Sampling

rand(d::Beta) = StatsFuns.RFunctions.betarand(d.α, d.β)

#### Fit model

# TODO: add MLE method (should be similar to Dirichlet)

# This is a moment-matching method (not MLE)
#
function fit{T<:Real}(::Type{Beta}, x::AbstractArray{T})
    x_bar = mean(x)
    v_bar = varm(x, x_bar)
    α = x_bar * (((x_bar * (1 - x_bar)) / v_bar) - 1)
    β = (1 - x_bar) * (((x_bar * (1 - x_bar)) / v_bar) - 1)
    Beta(α, β)
end


text:
parse_currency(currency) = Unitful.uparse(currency, unit_context=UnitfulAssets)

pyname(o::PyObject) = o.__class__.__name__

function check_pytype(x, o::PyObject)
    if string(nameof(x)) != pyname(o)
        error("Type mismatch: $(string(nameof(x))) != $(pyname(o))")
    end
    return nothing
end


text:
#
# Correctness Tests
#

using Compat: uninitialized
using Test, Dates, Printf
using DataArrays

my_tests = ["abstractarray.jl",
            "booleans.jl",
            "constructors.jl",
            "containers.jl",
            "conversions.jl",
            "data.jl",
            "dataarray.jl",
            "datamatrix.jl",
            "indexing.jl",
            "linalg.jl",
            "nas.jl",
            "operators.jl",
            "reduce.jl",
            "reducedim.jl",
            "broadcast.jl",
            "padding.jl",
            "pooleddataarray.jl",
            "extras.jl",
            "sort.jl",
            "statistics.jl",
            "literals.jl",
            "newtests/dataarray.jl",
            "newtests/datamatrix.jl",
            "newtests/datavector.jl"]

for test in my_tests
    include(test)
end

@testset "Ambiguities" begin
    @test_broken isempty(detect_ambiguities(DataArrays, Base, Core))
end


text:
#### Path

# First just do it for the case where there we assume
# tainted gotoifnots do not go in a loop!
# TODO: write a thing to detect this! (overdub predicates only in tainted ifs)
# implement snapshotting function state as an optimization for branch exploration
mutable struct Path
    path::BitVector
    cursor::Int
end

Path() = Path([], 1)

function increment!(bitvec)
    for i=1:length(bitvec)
        if bitvec[i] === true
            bitvec[i] = false
        else
            bitvec[i] = true
            break
        end
    end
end

function reset!(p::Path)
    p.cursor=1
    increment!(p.path)
    nothing
end

function alldone(p::Path) # must be called at the end of the function!
    all(identity, p.path)
end

function this_here_predicate!(p::Path)
    if p.cursor > length(p.path)
        push!(p.path, false)
    else
        p.path[p.cursor]
    end
    val = p.path[p.cursor]
    p.cursor+=1
    val
end

alldone(c::SparsityContext) = alldone(c.metadata[2])
reset!(c::SparsityContext) = reset!(c.metadata[2])
this_here_predicate!(c::SparsityContext) = this_here_predicate!(c.metadata[2])

#=
julia> p=Path()
Path(Bool[], 1)

julia> alldone(p) # must be called at the end of a full run
true

julia> this_here_predicate!(p)
false

julia> alldone(p) # must be called at the end of a full run
false

julia> this_here_predicate!(p)
false

julia> p
Path(Bool[false, false], 3)

julia> alldone(p)
false

julia> reset!(p)

julia> p
Path(Bool[true, false], 1)

julia> this_here_predicate!(p)
true

julia> this_here_predicate!(p)
false

julia> alldone(p)
false

julia> reset!(p)

julia> p
Path(Bool[false, true], 1)

julia> this_here_predicate!(p)
false

julia> this_here_predicate!(p)
true

julia> reset!(p)

julia> this_here_predicate!(p)
true

julia> this_here_predicate!(p)
true

julia> alldone(p)
true
=#


text:
using Dates

lo, hi = 2008, 2121
xmas = collect(Date(lo, 12, 25):Year(1):Date(hi, 12, 25))
filter!(xmas) do dt
    dayofweek(dt) == Dates.Sunday
end

println("Years from $lo to $hi having Christmas on Sunday: ")
foreach(println, year.(xmas))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ###############
# # recreate manual-traitdispatch.jl: f2

# @traitfn ttf2{X,Y<:Integer; D1{Y}, D4{X,Y}}(x::X,y::Y) = x + sin(y)
# @traitfn ttf2{S,T<:Integer; D1{S}, D1{T}  }(s::S,t::T) = sin(s) - sin(t)
# @traitfn ttf2{X,Y<:AbstractFloat; D1{X}, D1{Y}  }(x::X,y::Y) = cos(x) - cos(y)

# @test ttf2(4,5)==(4 + sin(5))
# @test ttf2(4,5.) == cos(4)-cos(5.)

# @test ttf2(MTT1(5),4) == sin(5)-sin(4)
# @test ttf2(MTT2(5),3) == sin(3)+5

# @test length(methods(f2))==length(methods(ttf2))

# note that some of below functions are also tested in perf/perf.jl.
# So a update here might need an update there.

################
@traitfn ft1{X,Y; Eq{X,Y}}(x::X,y::Y) = x==y
ff1(x,y) = x==y
@test ff1(4,5)==ft1(4,5)

#################
@traitfn function ft2{X,Y; Arith{X,Y}}(x::X,y::Y)
    out = zero(promote(x,y)[1])
    for xe in 1:round(Int,x)
        out += xe + y
    end
    out
end
function ff2{X,Y}(x::X,y::Y)
    out = zero(promote(x,y)[1])
    for xe in 1:round(Int,x)
        out += xe + y
    end
    out
end

@test ff2(7.3,5.)==ft2(7.3,5.)

##################

@traitfn ft3{X1,X2; IsIterable{X1}}(x::X1,y::X2) = (out=X2[]; for i in x; push!(out,i+y) end; out)
ff3{Y}(x,y::Y) = Y[i+y for i in x]

@test ff3([1:10;],5)==ft3([1:10;],5)
@test ff3(BigFloat[1:10;],BigInt(5))==ft3(BigFloat[1:10;],BigInt(5))

###############
@traitdef MyTr{X,Y} begin
    foobar(X,Y) -> Bool
end
type A
    a::Int
end
foobar(a::A, b::A) = a.a==b.a
@test istrait(MyTr{A,A})  # true
@test istrait(MyTr{Int,Int})==false

# make a function which dispatches on traits:
@traitfn ft111{X,Y; Cmp{X,Y}}(x::X,y::Y)  = x>y ? 5 : 6
@traitfn ft111{X,Y; MyTr{X,Y}}(x::X,y::Y) = foobar(x,y) ? -99 : -999

@test ft111(4,5)==6
@test ft111(A(5), A(6))==-999

@test_throws TraitException ft111("asdf", 5)
foobar(a::AbstractString, b::Int) = length(a)==b
@test ft111("asdf", 4)==-99

## dispatch using subtraits
@traitdef MyTr2{X,Y} <: MyTr{X,Y} begin
    bar(X,Y) -> Bool
end

@traitfn gt1{X,Y; MyTr2{X,Y}}(x::X,y::Y)  = "MyTr2"
@traitfn gt1{X,Y; MyTr{X,Y}}(x::X,y::Y)   = "MyTr"

type B1
    a::Int
end
foobar(a::B1, b::B1) = a.a==b.a
type B2
    a::Int
end
foobar(a::B2, b::B2) = a.a==b.a
bar(a::B2, b::B2) = a.a==b.a
@test istrait(MyTr{B1,B1})  # true
@test istrait(MyTr2{B1,B1})==false
@test istrait(MyTr2{B2,B2})

@test gt1(B1(1), B1(1))=="MyTr"
@test gt1(B2(1), B2(1))=="MyTr2"

##########
## adding trait methods to existing functions:
@traitdef MyTr7{X} begin
    bobo(X) -> AbstractString
end

abstract AMTyp7
type MTyp71 <: AMTyp7
end
type MTyp72 <: AMTyp7
end
# add to above trait
@traitimpl MyTr7{MTyp71} begin
    bobo(X::MTyp71) = "Yeah"
end

import Base.sin
@traitfn sin{X; MyTr7{X}}(x::X) = bobo(x)
@test sin(MTyp71())=="Yeah"
@test_throws TraitException sin(MTyp72())

@traitfn Base.cos{X; MyTr7{X}}(x::X) = bobo(x)
@test cos(MTyp71())=="Yeah"
@test_throws TraitException cos(MTyp72())

######
# Ambiguities
######
@traitdef TrTr1{X} begin
    len1(X)
end

@traitfn tf7465{X<:Integer,Y; TrTr1{X}}(x::X,y::Y) = x-y
@traitfn tf7465{X<:Integer,Y; TrTr1{Y}}(x::X,y::Y) = x+y

@traitimpl TrTr1{Int} begin
    len1(x::Int) = 5
end

# now ambiguous:
@test_throws Traits.TraitException tf7465(5,6)
# resolve
@traitfn tf7465{X<:Integer,Y; TrTr1{X}, TrTr1{Y}}(x::X,y::Y) = x*y
@test tf7465(5,6)==5*6

@traitdef TrTr22{X,Y} <: TrTr1{X}, TrTr1{Y} begin
    # empty
end

@traitfn tf7465{X<:Integer,Y; TrTr22{X,Y}}(x::X,y::Y) = x*y*1000
# errors again because ambigours again
@test_throws Traits.TraitException tf7465(5,6)

## single argument ambiguities
####
@traitdef TrTr2{X} begin
    len2(X)
end

@traitfn tttf{X; TrTr1{X}}(x::X) = len1(x)
@traitfn tttf{X; TrTr2{X}}(x::X) = len2(x)

@test tttf(5)==len1(5)  # works as TrTr2{Int} is not a trait

@traitimpl TrTr1{Float64} begin
    len1(x::Float64) = 5
end
@traitimpl TrTr2{Float64} begin
    len2(x::Float64) = 55
end

@test_throws Traits.TraitException tttf(5.)  # errors

@traitfn tttf{X; TrTr1{X}, TrTr2{X}}(x::X) = len2(x) # pick len2 over len1
@test tttf(5.)==len2(5.)

# # this can also be resolved by using a sub-trait, essentially equivalent:
# @traitdef STrTr{X} <: TrTr1{X}, TrTr2{X} begin
#     # empty
# end
# # no need to implement it for Float64: it's already a member:
# @test istrait(STrTr{Float64})

## test resolution of single argument ambiguities with subtraits
####
# came up here: https://github.com/mauro3/Traits.jl/pull/4#issuecomment-69742371
@traitdef U1{X} begin
    u1(X)
end
@traitdef U2{X} <: U1{X} begin
    u2(X)
end

@traitdef UU1{X} <: U1{X} begin
    uu1(X)
end
@traitdef UU2{X} <: U2{X} begin
    uu2(X)
end

@traitimpl U1{Int} begin
    u1(x::Int) = 1
end
@traitimpl U2{Int} begin
    u2(x::Int) = 2
end
@traitimpl UU1{Int} begin
    uu1(x::Int) = 11
end
@traitimpl UU2{Int} begin
    uu2(x::Int) = 12
end

@traitfn tttf238{X; UU1{X}}(x::X) = "this should loose"
@traitfn tttf238{X; UU2{X}}(x::X) = "this should win"

# This test in traitdispatch.jl should probably pass:
if dispatch_bug1
    @test_throws Traits.TraitException tttf238(5)=="this should win"
else
    @test tttf238(5)=="this should win"
end

# however if U2 were a subtrait of something else but U1 then dispatch should be ambiguous:
@traitdef V1{X} begin
    v1(X)
end
@traitdef Other1{X} begin
    other1(X)
end

@traitdef V2{X} <: Other1{X} begin
    v2(X)
end

@traitdef VV1{X} <: V1{X} begin
    vv1(X)
end
@traitdef VV2{X} <: V2{X} begin
    vv2(X)
end

@traitimpl V1{Int} begin
    v1(x::Int) = 1
end
@traitimpl Other1{Int} begin
    other1(x::Int) = -1
end
@traitimpl V2{Int} begin
    v2(x::Int) = 2
end
@traitimpl VV1{Int} begin
    vv1(x::Int) = 11
end
@traitimpl VV2{Int} begin
    vv2(x::Int) = 12
end

@traitfn tttf240{X; VV1{X}}(x::X) = "neither should win 1"
@traitfn tttf240{X; VV2{X}}(x::X) = "neither should win 2"
@test_throws Traits.TraitException tttf240(5)


text:
using Libdl

const _DEPS_FILE = joinpath(dirname(@__FILE__), "deps.jl")
if isfile(_DEPS_FILE)
    rm(_DEPS_FILE)
end

function write_depsfile(path)
    open(_DEPS_FILE, "w") do f
        println(f, "const libcplex = \"$(escape_string(path))\"")
    end
end

function library_name(v)
    cpx_prefix = Sys.iswindows() ? "" : "lib"
    return "$(cpx_prefix)cplex$(v).$(Libdl.dlext)"
end

function possible_path(cplex_studio_path::AbstractString)
    if Sys.iswindows()
        return escape_string("C:\\Program Files\\IBM\\ILOG\\$cplex_studio_path\\cplex\\bin\\x64_win64\\")
    elseif Sys.isapple()
        return "/Applications/$cplex_studio_path/cplex/bin/x86-64_osx/"
    else
        return "/opt/$cplex_studio_path/cplex/bin/x86-64_linux/"
    end
end

function get_error_message_if_not_found()
    return """
    Unable to install CPLEX.jl.

    The versions of CPLEX supported by CPLEX.jl are:

    * 12.10
    * 20.1

    You must download and install one of these versions separately.

    You should set the `CPLEX_STUDIO_BINARIES` environment variable to point to
    the install location then try again. For example (updating the path to the
    correct location if needed):

    ```
    ENV["CPLEX_STUDIO_BINARIES"] = "$(possible_path("CPLEX_Studio201"))"
    import Pkg
    Pkg.add("CPLEX")
    Pkg.build("CPLEX")
    ```

    See the CPLEX.jl README at https://github.com/jump-dev/CPLEX.jl for further
    instructions.
    """
end

function try_local_installation()
    # Find the path to the CPLEX executable.
    cplex_path = try
        @static if Sys.isapple() || Sys.isunix()
            dirname(strip(read(`which cplex`, String)))
        elseif Sys.iswindows()
            dirname(strip(read(`where cplex`, String)))
        end
    catch
        nothing
    end

    # Iterate through a series of places where CPLEX could be found: either in
    # the path (directly the callable library or the CPLEX executable) or from
    # an environment variable.
    cpxvers = [
        "1210", "12100",
        "201", "2010", "20100",
    ]
    base_env = "CPLEX_STUDIO_BINARIES"

    libnames = String["cplex"]
    for v in reverse(cpxvers)
        name = library_name(v)
        push!(libnames, name)
        if cplex_path !== nothing
            push!(libnames, joinpath(cplex_path, name))
        end
        for env in [base_env, base_env * v]
            if !haskey(ENV, env)
                continue
            end
            for d in split(ENV[env], ';')
                push!(libnames, joinpath(d, name))
            end
        end

        for product in ["CPLEX_Studio$v", "CPLEX_Enterprise_Server$v/CPLEX_Studio"]
            guessed_file = joinpath(possible_path(product), name)
            if isfile(guessed_file)
                push!(libnames, guessed_file)
            end
        end
    end

    # Perform the actual search in the potential places.
    for l in libnames
        d = Libdl.dlopen_e(l)
        if d == C_NULL
            continue
        end
        write_depsfile(Libdl.dlpath(d))
        @info("Using CPLEX found in location `$(l)`")
        return
    end

    error(get_error_message_if_not_found())
end

function try_ci_installation()
    CPLEX_VERSION = ENV["CPLEX_VERSION"]
    url = ENV["SECRET_CPLEX_URL_" * CPLEX_VERSION]
    local_filename = joinpath(@__DIR__, "libcplex.so")
    download(url, local_filename)
    write_depsfile(local_filename)
end

if get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
    # We need to be able to install and load this package without error for
    # Julia's registry AutoMerge to work. Just write a fake libcplex path.
    write_depsfile("julia_registryci_automerge")
elseif get(ENV, "SECRET_CPLEX_URL_12100", "") != ""
    try_ci_installation()
else
    try_local_installation()
end


text:
using Test
using SynthControl, TreatmentPanels, DataFrames, Dates

@testset "SimpleSCM" begin

  outcome_A = collect(range(1.0; stop = 3.0, length = 10))
  outcome_B = collect(range(2.0; stop = 1.0, length = 10))
  outcome_C = 0.7 .* outcome_A .+ 0.3 .* outcome_B
  outcome_C[8:end] .+= 1.0

  df = DataFrame(id = repeat(["A", "B", "C"], inner = 10),
    time = repeat(Date(2000):Year(1):Date(2009), 3),
    value = [outcome_A; outcome_B; outcome_C])

  tp = BalancedPanel(df, "C" => Date(2007);
    outcome_var = :value, id_var = :id, t_var = :time)

  # Test that BalancedPanel works
  @test tp isa BalancedPanel{SingleUnitTreatment{Continuous}}

  s = SimpleSCM(tp)

  @test !(isfitted(s))

  fit!(s)

  @test isfitted(s)

end

@testset "Synthetic Differences-in-Differences" begin

  bp = load_germany_panel()

  @test bp isa BalancedPanel{SingleUnitTreatment{Continuous}}

  sdid_model = SyntheticDiD(bp)

  fit!(sdid_model)

  @test only(sdid_model.τ̂) ≈ -3121.95 atol=0.1
end

@testset "Data sets" begin

  df = load_germany()
  @test only(df[df.country .== "USA" .&& df.year .== 1961, :gdp]) == 2929

  bp = load_germany_panel()
  @test typeof(bp) <: BalancedPanel

  df = load_basque()
  @test only(df[df.regionno .== 1 .&& df.year .== 1957, :gdpcap]) == 2.60361314875437

  bp = load_basque_panel()
  @test typeof(bp) <: BalancedPanel

  df = load_brexit()
  @test only(df[df.country .== "Australia" .&& df.quarter .== Date(2009, 1, 1), :realgdp]) == 1.04

  bp = load_brexit_panel()
  @test typeof(bp) <: BalancedPanel

end

max_stars_repo_name: nilshg/SynthControl

max_issues_repo_name: nilshg/SynthControl

max_forks_repo_name: nilshg/SynthControl


text:
using Random
using Base.Meta: isexpr

const rseed = Ref(Random.GLOBAL_RNG)  # to get new random directories (see julia #24445)
if isempty(methods(Random.seed!, Tuple{typeof(rseed[])}))
    # Julia 1.3-rc1 doesn't have this, fixed in https://github.com/JuliaLang/julia/pull/32961
    Random.seed!(rng::typeof(rseed[])) = Random.seed!(rng, nothing)
end
function randtmp()
    Random.seed!(rseed[])
    dirname = joinpath(tempdir(), randstring(10))
    rseed[] = Random.GLOBAL_RNG
    return dirname
end

const to_remove = String[]

function newtestdir()
    testdir = randtmp()
    mkdir(testdir)
    push!(to_remove, testdir)
    push!(LOAD_PATH, testdir)
    return testdir
end

@static if Sys.isapple()
    const mtimedelay = 2.1  # so the defining files are old enough not to trigger mtime criterion
else
    const mtimedelay = 0.1
end

yry() = (sleep(mtimedelay); revise(); sleep(mtimedelay))

function collectexprs(rex::Revise.RelocatableExpr)
    items = []
    for item in Revise.LineSkippingIterator(rex.ex.args)
        push!(items, isa(item, Expr) ? Revise.RelocatableExpr(item) : item)
    end
    items
end

function get_docstring(obj)
    while !isa(obj, AbstractString)
        fn = fieldnames(typeof(obj))
        if :content ∈ fn
            obj = obj.content[1]
        elseif :code ∈ fn
            obj = obj.code
        else
            error("unknown object ", obj)
        end
    end
    return obj
end

function get_code(f, typ)
    # Julia 1.5 introduces ":code_coverage_effect" exprs
    ci = code_typed(f, typ)[1].first
    code = copy(ci.code)
    while !isempty(code) && isexpr(code[1], :code_coverage_effect)
        popfirst!(code)
    end
    return code
end

function do_test(name)
    runtest = isempty(ARGS) || name in ARGS
    # Sometimes we get "no output received for 10 minutes" on CI,
    # to debug this it may be useful to know what test is being run.
    runtest && haskey(ENV, "CI") && println("Starting test ", name)
    return runtest
end

if !isempty(ARGS) && "REVISE_TESTS_WATCH_FILES" ∈ ARGS
    Revise.watching_files[] = true
    println("Running tests with `Revise.watching_files[] = true`")
    idx = findall(isequal("REVISE_TESTS_WATCH_FILES"), ARGS)
    deleteat!(ARGS, idx)
end


text:
#!/usr/bin/env julia

# Copyright 2022 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import assignment8
using Test

@testset "assignment8.jl" begin
    datafile = realpath(dirname(@__FILE__)*"/../data/data.csv")
    strain, stress = assignment8.convert_to_true_stress_and_strain(datafile)
    @test all(strain[100:102] .≈ [0.0005020476229990764, 0.0005009517727055945, 0.0005108676748850077])
    @test all(stress[512:515] .≈ [47382.09298647613, 47360.04462664843, 47414.967086017925, 47423.4567510546])
    @test assignment8.compute_toughness(datafile) ≈ 70836.1434834513
end