text:
process(tea, args...;kw...) = tea(args...;kw...)

function Base.put!(p::PotID, flavor)
    try
        put!(p[], flavor)
    catch e
        # TODO add test
        if e isa PotNotRegisteredError
            rethrow(e)
        else
            @error e
            boil(p)
            put!(p, flavor)
        end
    end
end

#####

struct CallMsg{A}
    args::A
    kw
    promise
end

is_prioritized(::CallMsg{<:Tuple{<:AbstractSysMsg}}) = true

function (p::PotID)(args...;kw...)
    promise = Promise(whereis(p))  # !!! the result should reside in the same place
    put!(p, CallMsg(args, kw.data, promise))
    promise
end

# ??? non specialized tea?
function process(tea, msg::CallMsg)
    try
        res = process(tea, msg.args...;msg.kw...)
        put!(msg.promise, res)
    catch err
        ce = CapturedException(err, catch_backtrace())
        put!(msg.promise, Failure(ce))
        rethrow(err)
    end
end

#####

function Base.:(|>)(x, p::PotID)
    put!(p, x)
    nothing
end

#####

struct GetPropMsg
    prop::Symbol
end

Base.getproperty(p::PotID, prop::Symbol) = p(GetPropMsg(prop))

process(tea, msg::GetPropMsg) = getproperty(tea, msg.prop)

##### SysMsg

struct Exit
end

const EXIT = Exit()

"""
    CloseWhenIdleMsg(t::Int)

Signal a Pot and its children to close the channel and release claimed resources if the Pot has been idle for `t` seconds.
"""
struct CloseWhenIdleMsg <: AbstractSysMsg
    t::Int
end

function process(tea, msg::CloseWhenIdleMsg)
    t_idle = (now() - _self().last_update) / Millisecond(1_000)
    if t_idle >= msg.t && isempty(_self().ch)
        for c in children()
            msg |> c
        end
        EXIT
    end
end

#####

"""
Close the active channel and remove the registered `Pot`.
"""
struct RemoveMsg <: AbstractSysMsg
end

Base.rm(p::PotID) = p(RemoveMsg())

function process(tea, msg::RemoveMsg)
    # !!! note the order
    for c in children()
        c(msg)[]
    end
    unregister(self())
    close(_self().ch)
    EXIT
end

#####

struct ResumeMsg <: AbstractSysMsg
end

process(tea, ::ResumeMsg) = nothing

#####

struct RestartMsg <: AbstractSysMsg
end

const RESTART = RestartMsg()

process(tea, ::RestartMsg) = RESTART

struct PreRestartMsg <: AbstractSysMsg
end

process(tea, ::PreRestartMsg) = nothing

struct PostRestartMsg <: AbstractSysMsg
end

process(tea, ::PostRestartMsg) = nothing

process(tea, ::Failure) = RESTART


text:
using Barcode
using Test
using FileIO
using Plots # it is in runtests.jl just to generate the image for README
            # it will be removed from tests once Documentation is ready

@testset "Encoding" begin
    @testset "Code128a" begin
        let code = Barcode.encode("A", :code128a)
            @test code == [
                "START A"
                "A"
                "CHECKSUM"
                "STOP"
            ]
        end
        let code = Barcode.encode("\x02A\x03", :code128a)
            @test code == [
                "START A"
                "STX"
                "A"
                "ETX"
                "CHECKSUM"
                "STOP"
            ]
        end
    end

    @testset "Code128b" begin
        let code = Barcode.encode("A", :code128b)
            @test code == [
                "START B"
                "A"
                "CHECKSUM"
                "STOP"
            ]
        end
        let code = Barcode.encode("aBc", :code128b)
            @test code == [
                "START B"
                "a"
                "B"
                "c"
                "CHECKSUM"
                "STOP"
            ]
        end
    end

    @testset "Code128b" begin
        let code = Barcode.encode("A", :code128b)
            @test code == [
                "START B"
                "A"
                "CHECKSUM"
                "STOP"
            ]
        end
        let code = Barcode.encode("aBc", :code128b)
            @test code == [
                "START B"
                "a"
                "B"
                "c"
                "CHECKSUM"
                "STOP"
            ]
        end
    end

    @testset "Code128c" begin
        let code = Barcode.encode("00", :code128c)
            @test code == [
                "START C"
                "00"
                "CHECKSUM"
                "STOP"
            ]
        end
        let code = Barcode.encode("012345", :code128c)
            @test code == [
                "START C"
                "01"
                "23"
                "45"
                "CHECKSUM"
                "STOP"
            ]
        end
    end

    @testset "Code128" begin
        let code = Barcode.encode("CSE370", :code128)
            @test code == [
                "START B"
                "C"
                "S"
                "E"
                "3"
                "7"
                "0"
                "CHECKSUM"
                "STOP"
            ]
        end

        let code = Barcode.encode("\x02Aa\tA0902a93892\x03000a\x04z", :code128)
            @test code == [
                "START A"
                "STX"
                "A"
                "SHIFT B"
                "a"
                "HT"
                "A"
                "CODE C"
                "09"
                "02"
                "CODE B"
                "a"
                "9"
                "CODE C"
                "38"
                "92"
                "CODE A"
                "ETX"
                "0"
                "0"
                "0"
                "SHIFT B"
                "a"
                "EOT"
                "CODE B"
                "z"
                "CHECKSUM"
                "STOP"
            ]
        end
    end
end

@testset "Decoding" begin
    for msg in (
        "ABC",
        "abc",
        "0123",
        "\x02ABC\x03",
        "\x02ABcZ\x03",
        "\x02AB012345\x03",
        "A b!\t0012\nZ z@\t0013\nAz `\t9999\naA \x7f\t1357"
    )
        code = Barcode.encode(msg, :code128)
        @test msg == Barcode.decode(code, :code128)
    end
end

@testset "Patterns" begin

    @testset "code128 subtypes" begin
        let pattern = Barcode.barcode_pattern("A", :code128a)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010000100" * # START A
                "10100011000" * # A
                "10100011000" * # checksum 33 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone

        end

        let pattern = Barcode.barcode_pattern("a", :code128b)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010010000" * # START B
                "10010110000" * # a
                "10010000110" * # checksum 66 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end

        let pattern = Barcode.barcode_pattern("00", :code128c)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010011100" * # START C
                "11011001100" * # 00
                "11001100110" * # checksum 2 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end
    end

    @testset "code128 auto" begin
        let pattern = Barcode.barcode_pattern("A", :code128)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010010000" * # START B
                "10100011000" * # A
                "10001011000" * # checksum 34 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end

        let pattern = Barcode.barcode_pattern("a", :code128)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010010000" * # START B
                "10010110000" * # a
                "10010000110" * # checksum 66 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end

        let pattern = Barcode.barcode_pattern("00", :code128)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010011100" * # START C
                "11011001100" * # 00
                "11001100110" * # checksum 2 pattern
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end

        let pattern = Barcode.barcode_pattern("\x02A\tB\x07\x03", :code128)
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010000100" * # START A
                "10010000110" * # \x02 = STX = Start of Text
                "10100011000" * # A
                "10000110100" * # \t = Horizontal Tab
                "10001011000" * # B
                "10011010000" * # \x07 = BEL = Bell
                "10000101100" * # \x03 = ETX = End of Text
                "10001100100" * # CHECKSUM
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end
    end

    @testset "mixed subtypes" begin

        let pattern = Barcode.barcode_pattern(
                [
                    "START A"
                    "A"
                    "B"
                    "SHIFT B"
                    "a"
                    "A"
                    "CODE C"
                    "00"
                    "CHECKSUM"
                    "STOP"
                ],
                :code128
            )
            @test pattern ==
                "00000000000" * # Quiet zone
                "11010000100" * # "START A" (103)
                "10100011000" * # "A" (33)
                "10001011000" * # "B" (34)
                "11110100010" * # "SHIFT B" (98)
                "10010110000" * # "a" (65)
                "10100011000" * # "A" (33)
                "10111011110" * # "CODE C" (99)
                "11011001100" * # "00" (0)
                "11000010010" * # CHECKSUM (75)
                # ( 103 + 1 * 33 + 2 * 34 + 3 * 98 + 4 * 65 + 5 * 33 + 6 * 99 + 7 * 0 ) % 103
                "11000111010" * # STOP
                "11" * # END
                "00000000000" # Quiet zone
        end

        let pattern = Barcode.barcode_pattern(
                [
                    "START C"
                    "FNC 1"
                    "42"
                    "18"
                    "40"
                    "20"
                    "50"
                    "CODE A"
                    "0"
                    "CHECKSUM"
                    "STOP"
                ],
                :code128,
            )

        end
    end
end

@testset "Depattern" begin
    for msg in (
        "ABC",
        "abc",
        "0123",
        "\x02ABC\x03",
        "\x02ABcZ\x03",
        "\x02AB012345\x03",
        "A b!\t0012\nZ z@\t0013\nAz `\t9999\naA \x7f\t1357"
    )
        code = Barcode.encode(msg, :code128)
        pattern = Barcode.barcode_pattern(code, :code128)
        @test msg == Barcode.barcode_decode(pattern, :code128)
    end
end

@testset "Images" begin
    @testset "save Images.jl" begin
        let zip_code = "12.345-678"
            zip_code = replace(zip_code, r"\s|\.|\-" => "")
            pattern = Barcode.barcode_pattern(zip_code, :code128)
            img = Barcode.barcode_img(pattern)
            @test FileIO.save("../img/zipcode_img.png", img) === nothing
        end

        let pattern = Barcode.barcode_pattern("\x02abc1234\x03", :code128)
            img = Barcode.barcode_img(pattern)
            @test FileIO.save("../img/abc1234.png", img) === nothing
        end
    end

    @testset "Plot barcode" begin
        let zip_code = "12.345-678"
            zip_code = replace(zip_code, r"\s|\.|\-" => "")
            pattern = Barcode.barcode_pattern(zip_code, :code128)
            x, w = Barcode.barcode_positions(pattern)
            Plots.plot(
                [x'; x' + w'], ones(2, length(x)), color = :black, fill = true,
                xlims = (1, length(pattern)),  ylims = (0, 1), border = :none,
                legend = nothing, size = (400, 150)
            )
            @test Plots.savefig("../img/zipcode_plot.png") === nothing
        end
    end
end

ext: jl


text:
const Instruction = Tuple{Char,Int}

function vfold(M::Matrix{Bool}, x::Int) :: Matrix{Bool}
  sY = size(M)[2]
  for i in 1:min(x-1,sY-x)
    M[:,x-i] = M[:,x-i] .| M[:,x+i]
  end
  return deepcopy(M[:,1:x-1])
end

function hfold(M::Matrix{Bool}, y::Int) :: Matrix{Bool}
  sX = size(M)[1]
  for i in 1:min(y-1,sX-y)
    M[y-i,:] = M[y-i,:] .| M[y+i,:]
  end
  return deepcopy(M[1:y-1,:])
end

fold(M::Matrix{Bool}, I::Instruction) = (I[1]=='x') ? vfold(M,I[2]) : hfold(M,I[2])

parseDot(line::String) = CartesianIndex(map(s->parse(Int,s),split(line,','))...)

function parseInstruction(line::String) :: Instruction
  s = split(line,' ')[3]
  return (first(s), parse(Int,s[3:end]))
end

function parseLine(line::String)
  if isempty(line)
    return
  elseif first(line) == 'f'
    parseInstruction(line)
  else
    parseDot(line)
  end
end

function matrix(dots::Vector{CartesianIndex{2}}) :: Matrix{Bool}
  sX = maximum(map(x->x[1],dots))+1
  sY = maximum(map(x->x[2],dots))+1
  M :: Matrix{Bool} = zeros(sX,sY)
  map(x->M[x] = true, dots)
  return deepcopy(M)
end

function readInput()
  f = open("./inputs/day13.txt")

  data = map( line -> parseLine(line), eachline(f))
  dots::Vector{CartesianIndex{2}} = filter(x->typeof(x) == CartesianIndex{2}, data)
  instructions::Vector{Instruction} = filter(x->typeof(x) == Instruction, data)
  return dots, instructions
end

function test()
  M :: Matrix{Bool} = zeros(5,5)
  M[1,5] = true
  M[2,4] = true
  M[3,2] = true
  M[4,1] = true
  M[5,1] = true
  M[5,4] = true
  M[5,5] = true
  display(M)
  M2 = hfold(M,2)
  display(M2)
  vfold(M2,3)
end

test()

dots, instructions = readInput()
dots = map(p -> p + CartesianIndex(1,1), dots)
instructions = map(x -> (x[1],x[2]+1), instructions)
M = matrix(dots)

M1 = deepcopy(M)

MA = [M1]
for f in instructions
  push!(MA,fold(deepcopy(MA[end]),f))
end

display(MA[end])


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, DataFrames, Plots
	using Distributions, LinearAlgebra, LaTeXStrings
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.32"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	N = 10^4

	SigY = [ 6 4 ;
			 4 9]
	muY = [15 ;
		   20]
	A = cholesky(SigY).L

	rngGens = [()->rand(Normal()),
			   ()->rand(Uniform(-sqrt(3),sqrt(3))),
			   ()->rand(Exponential())-1]

	rv(rg) = A*[rg(),rg()] + muY

	data = [[rv(r) for _ in 1:N] for r in rngGens]
end

# ╔═╡ e752604c-2211-11eb-36e4-9d0ca870afa5
stats(data) = begin
    data1, data2 = first.(data),last.(data)
    (
		round(mean(data1),digits=2),
		round(mean(data2),digits=2),
    	round(var(data1),digits=2),
		round(var(data2),digits=2),
        round(cov(data1,data2),digits=2)
	)
end

# ╔═╡ e753223c-2211-11eb-1b21-77985295723c
begin
	df = DataFrame()
	for d in data
		r = stats(d)
		append!(df, DataFrame(:Mean1 => r[1], :Mean2 => r[2], :Var1 => r[3],
				:Var2 => r[4], :Cov => r[5]))
	end
	df
end

# ╔═╡ e75f86ee-2211-11eb-24b6-3921f08d4bfd
begin
	scatter(first.(data[1]), last.(data[1]), c=:blue, ms=1, msw=0, label="Normal")
	scatter!(first.(data[2]), last.(data[2]), c=:red, ms=1, msw=0, label="Uniform")
	scatter!(first.(data[3]),last.(data[3]),c=:green, ms=1,msw=0,label="Exponential",
		xlims=(0,40), ylims=(0,40), legend=:bottomright, ratio=:equal,
		xlabel=L"X_1", ylabel=L"X_2")
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.32"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═e752604c-2211-11eb-36e4-9d0ca870afa5
# ╠═e753223c-2211-11eb-1b21-77985295723c
# ╠═e75f86ee-2211-11eb-24b6-3921f08d4bfd
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
"""
    Code for convenience functions for calculating CCFs for an AbstractChunkOfSpectrum (see RvSpectMLBase)
Author: Eric Ford
Created: August 2020
"""

"""  `calc_ccf_chunk!( cccf_out, chunk, ccf_plan )`
Convenience function to compute CCF for one chunk of spectrum, evaluated using mask_shape and line list from ccf plan

# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Returns:
- `ccf_out`
"""
function calc_ccf_chunk!(ccf_out::AbstractArray{T1,1}, chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T2} = chunk.var,
                 assume_sorted::Bool = false ) where { T1<:Real, T2<:Real, PlanT<:AbstractCCFPlan }
  @assert assume_sorted || issorted( plan.line_list.λ )
  @assert length(ccf_out) == calc_length_ccf_v_grid(plan)
  ccf_1D!(ccf_out, chunk.λ, chunk.flux, plan; assume_sorted=true)
  return ccf_out
end

"""  `calc_ccf_and_var_chunk!( chunk, ccf_plan )`
Convenience function to compute CCF and variance of each "CCF pixel" for one chunk of spectrum, evaluated using mask_shape and line list from `ccf_plan`.
# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `ccf_var_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `var`:  `AbstractArray` with variance to use for each pixel (overides value in chunk)
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Returns Named Tuple with:
- `ccf_out`:
- `ccf_var_out`:
"""
function calc_ccf_and_var_chunk!(ccf_out::AbstractArray{T1,1}, ccf_var_out::AbstractArray{T2,1},
                chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T3} = chunk.var, ccf_var_scale::Real = 1.0,
                 assume_sorted::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, PlanT<:AbstractCCFPlan }
  @assert assume_sorted || issorted( plan.line_list.λ )
  @assert length(ccf_out) == calc_length_ccf_v_grid(plan)
  ccf_1D!(ccf_out, ccf_var_out, chunk.λ, chunk.flux, var, plan; ccf_var_scale=ccf_var_scale, assume_sorted=true)
  return (ccf=ccf_out, ccf_var=ccf_var_out)
end

"""  `calc_ccf_and_covar_chunk!( chunk, ccf_plan )`
Convenience function to compute CCF and covariance of each pair of "CCF pixels" for one chunk of spectrum, evaluated using mask_shape and line list from `ccf_plan`.
# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `ccf_covar_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `var`:  `AbstractArray` with variance to use for each pixel (overides value in chunk)
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Returns Named Tuple with:
- `ccf_out`:
- `ccf_covar_out`:
"""
function calc_ccf_and_covar_chunk!(ccf_out::AbstractArray{T1,1}, ccf_covar_out::AbstractArray{T2,2},
                chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T3} = chunk.var, ccf_var_scale::Real = 1.0,
                 assume_sorted::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, PlanT<:AbstractCCFPlan }
  @assert assume_sorted || issorted( plan.line_list.λ )
  @assert length(ccf_out) == calc_length_ccf_v_grid(plan)
  ccf_1D!(ccf_out, ccf_covar_out, chunk.λ, chunk.flux, var, plan; ccf_var_scale=ccf_var_scale, assume_sorted=true)
  return (ccf=ccf_out, ccf_covar=ccf_covar_out)
end

"""  `calc_ccf_chunk( chunk, ccf_plan )`
Convenience function to compute CCF for one chunk of spectrum.
# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
- `calc_ccf_var`:  if true also computes estimate of variance for each value of ccf
# Returns:
- CCF for one chunk of spectrum, evaluated using mask_shape and line list from ccf plan
"""
function calc_ccf_chunk(chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T} = chunk.var, ccf_var_scale::Real = 1.0, Δfwhm::Real = 0,
                 assume_sorted::Bool = false, calc_ccf_var::Bool = false ) where { T<:Real, PlanT<:AbstractCCFPlan }
  if Δfwhm > 0
     this_plan_for_chunk = copy(plan)
     increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
  else
     this_plan_for_chunk = plan
  end
  len_v_grid = calc_length_ccf_v_grid(plan)
  ccf_out = zeros(len_v_grid)
  if calc_ccf_var
    ccf_var_out = zeros(len_v_grid)
    return calc_ccf_and_var_chunk!(ccf_out, ccf_var_out, chunk, this_plan_for_chunk, var=var, ccf_var_scale=ccf_var_scale, assume_sorted=assume_sorted )
  else
    return calc_ccf_chunk!(ccf_out, chunk, this_plan_for_chunk, var=var, assume_sorted=assume_sorted )
  end
end

"""  `calc_ccf_and_var_chunk( chunk, ccf_plan )`
Convenience function to compute CCF and variance of each "CCF pixel" for one chunk of spectrum, evaluated using mask_shape and line list from `ccf_plan`.
# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `ccf_var_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `var`:  `AbstractArray` with variance to use for each pixel (overides value in chunk)
`- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Returns Named Tuple with:
- `ccf_out`:
- `ccf_var_out`:
"""
function calc_ccf_and_var_chunk(chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T} = chunk.var, ccf_var_scale::Real = 1.0, Δfwhm::Real = 0,
                 assume_sorted::Bool = false ) where { T<:Real, PlanT<:AbstractCCFPlan }
  @assert assume_sorted || issorted( plan.line_list.λ )
  if Δfwhm > 0
     this_plan_for_chunk = copy(plan)
     increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
  else
     this_plan_for_chunk = plan
  end
  len_v_grid = calc_length_ccf_v_grid(plan)
  ccf_out = zeros(len_v_grid)
  ccf_var_out = zeros(len_v_grid)
  calc_ccf_and_var_chunk!(ccf_out, ccf_var_out, chunk, this_plan_for_chunk, var=var, ccf_var_scale=ccf_var_scale, assume_sorted=true )
  return (ccf=ccf_out, ccf_var=ccf_var_out)
end

"""  `calc_ccf_and_covar_chunk( chunk, ccf_plan )`
Convenience function to compute CCF and covariance of each pair of "CCF pixels" for one chunk of spectrum, evaluated using mask_shape and line list from `ccf_plan`.
# Inputs:
- `ccf_out`:  `AbstractArray` to store output
- `ccf_covar_out`:  `AbstractArray` to store output
- `chunk`: ChunkOfSpectrum to compute CCF for
- `ccf_plan`: for now, just a BasicCCFPlan that provides line_list, mask_shape and other parameters for calculating CCF
# Optional Arguments:
- `var`:  `AbstractArray` with variance to use for each pixel (overides value in chunk)
`- `assume_sorted`:  if true, skips checking the line_list is sorted by wavelength
# Returns Named Tuple with:
- `ccf_out`:
- `ccf_covar_out`:
"""
function calc_ccf_and_covar_chunk(chunk::AbstractChunkOfSpectrum, plan::PlanT = BasicCCFPlan()
                 ; var::AbstractVector{T} = chunk.var, ccf_var_scale::Real = 1.0, Δfwhm::Real = 0,
                 assume_sorted::Bool = false ) where { T<:Real, PlanT<:AbstractCCFPlan }
  @assert assume_sorted || issorted( plan.line_list.λ )
  if Δfwhm > 0
     this_plan_for_chunk = copy(plan)
     increase_mask_fwhm!(this_plan_for_chunk,Δfwhm)
  else
     this_plan_for_chunk = plan
  end
  len_v_grid = calc_length_ccf_v_grid(plan)
  ccf_out = zeros(len_v_grid)
  ccf_covar_out = zeros(len_v_grid,len_v_grid)
  calc_ccf_and_covar_chunk!(ccf_out, ccf_covar_out, chunk, this_plan_for_chunk, var=var, ccf_var_scale=ccf_var_scale, assume_sorted=true )
  return (ccf=ccf_out, ccf_covar=ccf_covar_out)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
module TestResampling

# using Revise
using Test
using MLJ
using MLJBase
import MLJModels
import StatsBase
using Distributed
import ComputationalResources: CPU1, CPUProcesses, CPUThreads
import Random.seed!
seed!(1234)

include("foobarmodel.jl")
@everywhere workers() include("foobarmodel.jl")

macro testset_accelerated(name::String, var, ex)
    testset_accelerated(name, var, ex)
end
macro testset_accelerated(name::String, var, opts::Expr, ex)
    testset_accelerated(name, var, ex; eval(opts)...)
end
function testset_accelerated(name::String, var, ex; exclude=[])
    final_ex = quote
        $var = CPU1()
        @testset $name $ex
    end
    resources = Any[CPUProcesses()]
    @static if VERSION >= v"1.3.0-DEV.573"
        push!(resources, CPUThreads())
    end
    for res in resources
        if any(x->typeof(res)<:x, exclude)
            push!(final_ex.args, quote
                $var = $res
                @testset $(name*" (accelerated with $(typeof(res).name))") begin
                    @test_broken false
                end
            end)
        else
            push!(final_ex.args, quote
                $var = $res
                @testset $(name*" (accelerated with $(typeof(res).name))") $ex
            end)
        end
    end
    return esc(final_ex)
end

@test CV(nfolds=6) == CV(nfolds=6)
@test CV(nfolds=5) != CV(nfolds=6)
@test MLJ.train_test_pairs(CV(), 1:10) !=
     MLJ.train_test_pairs(CV(shuffle=true), 1:10)
@test MLJ.train_test_pairs(Holdout(), 1:10) !=
     MLJ.train_test_pairs(Holdout(shuffle=true), 1:10)

@testset "checking measure/model compatibility" begin
    model = ConstantRegressor()
    y = rand(4)
    override=false
    @test MLJ._check_measure(:junk, :junk, :junk, :junk, true) == nothing
    @test_throws(ArgumentError,
                  MLJ._check_measure(model, rms, y, predict, override))
    @test MLJ._check_measure(model, rms, y, predict_mean, override) == nothing
    @test MLJ._check_measure(model, rms, y, predict_median, override) == nothing
    y=categorical(collect("abc"))
    @test_throws(ArgumentError,
                 MLJ._check_measure(model, rms, y, predict_median, override))
    model = ConstantClassifier()
    @test_throws(ArgumentError,
       MLJ._check_measure(model, misclassification_rate, y, predict, override))
    @test MLJ._check_measure(model, misclassification_rate, y,
                            predict_mode, override) == nothing
    model = MLJModels.DeterministicConstantClassifier()
    @test_throws ArgumentError MLJ._check_measure(model, cross_entropy, y,
                            predict, override)
end

@testset_accelerated "folds specified" accel (exclude=[CPUProcesses],) begin
    x1 = ones(10)
    x2 = ones(10)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    my_rms(yhat, y) = sqrt(mean((yhat -y).^2))
    my_mav(yhat, y) = abs.(yhat - y)
    MLJ.reports_each_observation(::typeof(my_mav)) = true

    resampling = [(3:10, 1:2),
                  ([1, 2, 5, 6, 7, 8, 9, 10], 3:4),
                  ([1, 2, 3, 4, 7, 8, 9, 10], 5:6),
                  ([1, 2, 3, 4, 5, 6, 9, 10], 7:8),
                  (1:8, 9:10)]

    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)

    # check detection of incompatible measure (cross_entropy):
    @test_throws ArgumentError evaluate!(mach, resampling=resampling,
                                         measure=[cross_entropy, rmslp1], acceleration=accel)
    result = evaluate!(mach, resampling=resampling,
                       measure=[my_rms, my_mav, rmslp1], acceleration=accel)
    v = [1/2, 3/4, 1/2, 3/4, 1/2]
    @test result.per_fold[1] ≈ v
    @test result.per_fold[2] ≈ v
    @test result.per_fold[3][1] ≈ abs(log(2) - log(2.5))
    @test ismissing(result.per_observation[1])
    @test result.per_observation[2][1] ≈ [1/2, 1/2]
    @test result.per_observation[2][2] ≈ [3/4, 3/4]
    @test result.measurement[1] ≈ mean(v)
    @test result.measurement[2] ≈ mean(v)
end

@testset_accelerated "holdout" accel begin
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0]

    @test MLJBase.show_as_constructed(Holdout)
    holdout = Holdout(fraction_train=0.75)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    result = evaluate!(mach, resampling=holdout,
                       measure=[rms, rmslp1], acceleration=accel)
    result = evaluate!(mach, verbosity=0, resampling=holdout, acceleration=accel)
    result.measurement[1] ≈ 2/3

    # test direct evaluation of a model + data:
    result = evaluate(model, X, y, verbosity=0,
                      resampling=holdout, measure=rms)
    @test result.measurement[1] ≈ 2/3

    X = (x=rand(100),)
    y = rand(100)
    mach = machine(model, X, y)
    evaluate!(mach, verbosity=0,
              resampling=Holdout(shuffle=true, rng=123), acceleration=accel)
    e1 = evaluate!(mach, verbosity=0,
                   resampling=Holdout(shuffle=true),
                   acceleration=accel).measurement[1]
    @test e1 != evaluate!(mach, verbosity=0,
                          resampling=Holdout(), acceleration=accel).measurement[1]
end

@testset_accelerated "cv" accel begin
    x1 = ones(10)
    x2 = ones(10)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    @test MLJBase.show_as_constructed(CV)
    cv=CV(nfolds=5)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    result = evaluate!(mach, resampling=cv, measure=[rms, rmslp1],
                       acceleration=accel)
    @test result.per_fold[1] ≈ [1/2, 3/4, 1/2, 3/4, 1/2]

    shuffled = evaluate!(mach, resampling=CV(shuffle=true),
                          acceleration=accel) # using rms default
    @test shuffled.measurement[1] != result.measurement[1]
end

@testset "stratified_cv" begin

    # check in explicit example:
    y = categorical(['c', 'a', 'b', 'a', 'c', 'x',
                 'c', 'a', 'a', 'b', 'b', 'b', 'b', 'b'])
    rows = [14, 13, 12, 11, 10, 9, 8, 7, 5, 4, 3, 2, 1]
    @test y[rows] == collect("bbbbbaaccabac")
    scv = StratifiedCV(nfolds=3)
    pairs = MLJ.train_test_pairs(scv, rows, nothing, y)
    @test pairs == [([12, 11, 10, 8, 5, 4, 3, 2, 1], [14, 13, 9, 7]),
                    ([14, 13, 10, 9, 7, 4, 3, 2, 1], [12, 11, 8, 5]),
                    ([14, 13, 12, 11, 9, 8, 7, 5], [10, 4, 3, 2, 1])]
    scv_random = StratifiedCV(nfolds=3, shuffle=true)
    pairs_random = MLJ.train_test_pairs(scv_random, rows, nothing, y)
    @test pairs != pairs_random

    # wrong target type throws error:
    @test_throws Exception MLJ.train_test_pairs(scv, rows, nothing, get.(y))

    # too many folds throws error:
    @test_throws Exception MLJ.train_test_pairs(StratifiedCV(nfolds=4),
                                                rows, nothing, y)

    # check class distribution is preserved in a larger randomized example:
    N = 30
    y = shuffle(vcat(fill(:a, N), fill(:b, 2N),
                        fill(:c, 3N), fill(:d, 4N))) |> categorical;
    d = fit(UnivariateFinite, y)
    pairs = MLJ.train_test_pairs(scv, 1:10N, nothing, y)
    folds = vcat(first.(pairs), last.(pairs))
    @test all([fit(UnivariateFinite, y[fold]) ≈ d for fold in folds])

end

@testset_accelerated "sample weights in evaluation" accel begin

    # cv:
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 2.0, 3.0, 1.0]
    w = 1:4
    cv=CV(nfolds=2)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    e = evaluate!(mach, resampling=cv, measure=l1,
                  weights=w, verbosity=0, acceleration=accel).measurement[1]
    @test e ≈ (1/3 + 13/14)/2
end

@testset_accelerated "resampler as machine" accel (exclude=[CPUProcesses],) begin
    N = 50
    X = (x1=rand(N), x2=rand(N), x3=rand(N))
    y = X.x1 -2X.x2 + 0.05*rand(N)

    ridge_model = FooBarRegressor(lambda=20.0)
    holdout = Holdout(fraction_train=0.75)
    resampler = Resampler(resampling=holdout, model=ridge_model, measure=mav)
    resampling_machine = machine(resampler, X, y)
    @test_logs((:info, r"^Training"), fit!(resampling_machine))
    e1=evaluate(resampling_machine).measurement[1]
    mach = machine(ridge_model, X, y)
    @test e1 ≈  evaluate!(mach, resampling=holdout,
                          measure=mav, verbosity=0,
                          acceleration=accel).measurement[1]
    ridge_model.lambda=1.0
    fit!(resampling_machine, verbosity=2)
    e2=evaluate(resampling_machine).measurement[1]
    @test e1 != e2
    resampler.weights = rand(N)
    fit!(resampling_machine, verbosity=0)
    e3=evaluate(resampling_machine).measurement[1]
    @test e3 != e2

    @test MLJBase.package_name(Resampler) == "MLJ"
    @test MLJBase.is_wrapper(Resampler)
    rnd = randn(5)
    @test evaluate(resampler, rnd) === rnd
end

struct DummyResamplingStrategy <: MLJ.ResamplingStrategy end

@testset_accelerated "custom strategy using resampling depending on X, y" accel begin
    function MLJ.train_test_pairs(resampling::DummyResamplingStrategy,
                              rows, X, y)
        train = filter(rows) do j
            y[j] == y[1]
        end
        test = setdiff(rows, train)
        return [(train, test),]
    end

    X = (x = rand(8), )
    y = categorical([:x, :y, :x, :x, :y, :x, :x, :y])
    @test MLJ.train_test_pairs(DummyResamplingStrategy(), 2:6, X, y) ==
        [([3, 4, 6], [2, 5]),]

    e = evaluate(ConstantClassifier(), X, y,
                 measure=misclassification_rate,
                 resampling=DummyResamplingStrategy(),
                 operation=predict_mode,
                 acceleration=accel)
    @test e.measurement[1] ≈ 1.0
end

@testset_accelerated "sample weights in training and evaluation" accel begin
    yraw = ["Perry", "Antonia", "Perry", "Antonia", "Skater"]
    X = (x=rand(5),)
    y = categorical(yraw)
    w = [1, 10, 1, 10, 5]

    # without weights:
    mach = machine(ConstantClassifier(), X, y)
    e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                  operation=predict_mode, measure=misclassification_rate,
                  acceleration=accel)
    @test e.measurement[1] ≈ 1.0

    # with weights in training and evaluation:
    mach = machine(ConstantClassifier(), X, y, w)
    e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                  operation=predict_mode, measure=misclassification_rate,
                  acceleration=accel)
    @test e.measurement[1] ≈ 1/3

    # with weights in training but overriden in evaluation:
    e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                  operation=predict_mode, measure=misclassification_rate,
                  weights = fill(1, 5), acceleration=accel)
    @test e.measurement[1] ≈ 1/2

    @test_throws(DimensionMismatch,
                 evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                           operation=predict_mode,
                           measure=misclassification_rate,
                           weights = fill(1, 100), acceleration=accel))

    @test_throws(ArgumentError,
                 evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                           operation=predict_mode,
                           measure=misclassification_rate,
                           weights = fill('a', 5), acceleration=accel))

    # resampling on a subset of all rows:
    model = @load KNNClassifier

    N = 200
    X = (x = rand(3N), );
    y = categorical(rand("abcd", 3N));
    w = rand(3N);
    rows = StatsBase.sample(1:3N, 2N, replace=false);
    Xsmall = selectrows(X, rows);
    ysmall = selectrows(y, rows);
    wsmall = selectrows(w, rows);

    mach1 = machine(model, Xsmall, ysmall, wsmall)
    e1 = evaluate!(mach1, resampling=CV(),
                   measure=misclassification_rate,
                   operation=predict_mode, acceleration=accel)

    mach2 = machine(model, X, y, w)
    e2 = evaluate!(mach2, resampling=CV(),
                   measure=misclassification_rate,
                   operation=predict_mode,
                   rows=rows, acceleration=accel)

    @test e1.per_fold ≈ e2.per_fold

    # resampler as machine with evaluation weights not specified:
    resampler = Resampler(model=model, resampling=CV();
                          measure=misclassification_rate,
                          operation=predict_mode)
    resampling_machine = machine(resampler, X, y, w)
    fit!(resampling_machine)
    e1 = evaluate(resampling_machine).measurement[1]
    mach = machine(model, X, y, w)
    e2 = evaluate!(mach, resampling=CV();
                   measure=misclassification_rate,
                   operation=predict_mode, acceleration=accel).measurement[1]
    @test e1 ≈ e2

    # resampler as machine with evaluation weights specified:
    weval = rand(3N);
    resampler = Resampler(model=model, resampling=CV();
                          measure=misclassification_rate,
                          operation=predict_mode,
                          weights=weval, acceleration=accel)
    resampling_machine = machine(resampler, X, y, w)
    fit!(resampling_machine)
    e1 = evaluate(resampling_machine).measurement[1]
    mach = machine(model, X, y, w)
    e2 = evaluate!(mach, resampling=CV();
                   measure=misclassification_rate,
                   operation=predict_mode,
                   weights=weval, acceleration=accel).measurement[1]
    @test e1 ≈ e2

end

end
true


text:
using ITensors

#
# DMRG calculation of the extended Hubbard model
# ground state wavefunction, and spin densities
#

let
  N = 20
  Npart = 10
  t1 = 1.0
  t2 = 0.2
  U  = 1.0
  V1 = 0.5

  sites = siteinds("Electron", N;
                   conserve_qns = true)

  ampo = AutoMPO()
  for b=1:N-1
    ampo += -t1, "Cdagup", b,    "Cup", b+1
    ampo += -t1, "Cdagup", b+1,  "Cup", b
    ampo += -t1, "Cdagdn", b,    "Cdn", b+1
    ampo += -t1, "Cdagdn", b+1,  "Cdn", b
    ampo +=  V1,   "Ntot", b,   "Ntot", b+1
  end
  for b=1:N-2
    ampo += -t2, "Cdagup", b,   "Cup", b+2
    ampo += -t2, "Cdagup", b+2, "Cup", b
    ampo += -t2, "Cdagdn", b,   "Cdn", b+2
    ampo += -t2, "Cdagdn", b+2, "Cdn", b
  end
  for i=1:N
    ampo += U, "Nupdn", i
  end
  H = MPO(ampo,sites)

  sweeps = Sweeps(6)
  maxdim!(sweeps,50,100,200,400,800,800)
  cutoff!(sweeps,1E-12)
  @show sweeps

  state = ["Emp" for n=1:N]
  p = Npart
  for i=N:-1:1
    if p > i
      println("Doubly occupying site $i")
      state[i] = "UpDn"
      p -= 2
    elseif p > 0
      println("Singly occupying site $i")
      state[i] = (isodd(i) ? "Up" : "Dn")
      p -= 1
    end
  end
  # Initialize wavefunction to be bond
  # dimension 10 random MPS with number
  # of particles the same as `state`
  psi0 = randomMPS(sites, state, 10)

  # Check total number of particles:
  @show flux(psi0)

  # Start DMRG calculation:
  energy, psi = dmrg(H, psi0, sweeps)

  upd = fill(0.0, N)
  dnd = fill(0.0, N)
  for j=1:N
    orthogonalize!(psi,j)
    psidag_j = dag(prime(psi[j], "Site"))
    upd[j] = scalar(psidag_j * op(sites, "Nup", j) * psi[j])
    dnd[j] = scalar(psidag_j * op(sites, "Ndn", j) * psi[j])
  end

  println("Up Density:")
  for j=1:N
    println("$j $(upd[j])")
  end
  println()

  println("Dn Density:")
  for j=1:N
    println("$j $(dnd[j])")
  end
  println()

  println("Total Density:")
  for j=1:N
    println("$j $(upd[j]+dnd[j])")
  end
  println()

  println("\nGround State Energy = $energy")

end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_commit(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)

	# Commitment state for each resource in each time step
	commit = zeros(G,T)
	for i in inputs["COMMIT"]
		commit[i,:] = value.(EP[:vCOMMIT])[i,:]
	end
	dfCommit = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone])
	dfCommit = hcat(dfCommit, DataFrame(commit, :auto))
	auxNew_Names=[Symbol("Resource");Symbol("Zone");[Symbol("t$t") for t in 1:T]]
	rename!(dfCommit,auxNew_Names)
	CSV.write(joinpath(path, "commit.csv"), dftranspose(dfCommit, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# inflation.jl - Estructura básica para cómputo de estimadores de inflación interanual

"""
    abstract type InflationFunction <: Function

Tipo abstracto para representar las funciones de inflación que operan sobre
[`CountryStructure`](@ref) y [`VarCPIBase`](@ref). Permiten computar la medida
de ritmo inflacionario interanual, el índice de precios dado por la metodología
y las variaciones intermensuales del índice de precios.
"""
abstract type InflationFunction <: Function end

# La mayoría de funciones de inflación devuelven una medida
"""
    num_measures(::InflationFunction)

Devuelve la cantidad de medidas devueltas por la función de inflación. Las
funciones de [`EnsembleFunction`](@ref) pueden computar varias medidas de
inflación simultáneamente.
"""
num_measures(::InflationFunction) = 1

# Funciones de inflación deben extender el método measure_name para devolver su nombre
"""
    measure_name(inflfn::InflationFunction)

Este método permite obtener el nombre convencional de una medida de inflación.
"""
measure_name(inflfn::InflationFunction) =
    error("Se debe extender el método `measure_name` para la función de inflación " * string(nameof(inflfn)))

# Funciones de inflación pueden extender opcionalmente un método measure_tag
# para indicar en archivos de resultados. Si no, se tomará por defecto el nombre
# del tipo de la función de inflación
"""
    measure_tag(inflfn::InflationFunction)

Obtiene una etiqueta de la medida de inflación. Se puede utilizar para guardar
como parámetro en archivos de resultados de evaluación.
"""
measure_tag(inflfn::InflationFunction) = string(nameof(inflfn))

# Método para obtener parámetros
"""
    params(inflfn::InflationFunction)
Método para obtener parámetros de la función de inflación. Devuelve una tupla
con el conjunto de parámetros utilizado por la función de inflación `inflfn`.
Este método debe redefinirse en las nuevas medidas de inflación si estas están
parametrizadas.
"""
params(inflfn::InflationFunction) = getproperty.(Ref(inflfn), propertynames(inflfn))

## Tipos para resultados, utilizados para el despacho de métodos
"Tipo abstracto para manejar el despacho de las funciones de inflación"
abstract type CPIResult end
"Tipo concreto único para obtener el índice de una función de inflación"
struct CPIIndex <: CPIResult end
"""
Tipo concreto único para obtener la variación intermensual de una función de inflación
"""
struct CPIVarInterm <: CPIResult end

## Esquema general de cómputo (programación genérica con tipos abstractos):
# - La función sobre CountryStructure devuelve la inflación interanual sobre
#   todas las bases que componen
# - Esta llama a la función de inflación que recibe `CPIIndex`.
# - Y esta a su vez, llama a la función de inflación que recibe `CPIVarInterm`.
# - Finalmente, se llama al método que opera sobre objetos `VarCPIBase` para
#   obtener la variación intermensual resumen. De tal manera que la mayoría de
#   funciones solamente requieren definir su operación sobre los contenedores
#   `VarCPIBase` y devolver una variación intermensual resumen.

function (inflfn::InflationFunction)(cs::CountryStructure)
    cpi_index = inflfn(cs, CPIIndex())
    varinteran(cpi_index)
end

function (inflfn::InflationFunction)(cs::CountryStructure, ::CPIIndex)
    v_interm = inflfn(cs, CPIVarInterm())
    capitalize!(v_interm, 100) # v_interm -> cpi_index
    v_interm
end

function (inflfn::InflationFunction)(cs::CountryStructure, ::CPIVarInterm)
    # Acá se llama a inflfn(base), en donde base es de tipo VarCPIBase. Esta
    # es la función que debe definirse para cualquier medida de inflación.
    v_interm = mapfoldl(inflfn, vcat, cs.base)
    v_interm
end

# Funciones de inflación deben extender el método que opera sobre VarCPIBase
# para computar la variación intermensual resumen
function (inflfn::InflationFunction)(::VarCPIBase)
    error("Se debe extender un método para computar la variación intermensual resumen de la medida de inflación " * string(nameof(inflfn)))
end


text:
__precompile__()

module EasyPhys

export plt,
       Fitter,
       set!,
       set_data!,
       free!,
       fix!,
       guess!,
       apply_mask!,
       ignore_outliers!,
       fit!,
       parameter_covariance,
       studentized_residuals,
       reduced_χ²,
       apply_f,
       plot,
       plot!

using Compat
import Compat.String

include("utilities.jl")
include("parameters.jl")
include("fitter.jl")
include("plot.jl")

end


text:
export Parameter
export make_parameter, share_parameter

type Parameter <: AbstractParameter
  name          :: AbstractString
  blob          :: Blob
  gradient      :: Blob
  initializer   :: Initializer
  regularizer   :: Regularizer
  constraint    :: Constraint
  learning_rate :: AbstractFloat # relative learning rate

  rc            :: RefCounter
end

# construct a new parameter
Parameter(name,blob,gradient,initializer,regularizer,constraint,lr) =
    Parameter(name, blob, gradient, initializer, regularizer, constraint, lr, RefCounter(1))

function make_parameter{N}(backend::Backend, name::AbstractString,
    data_type::Type, dims::NTuple{N,Int},
    init::Initializer, regu::Regularizer, cons::Constraint, lr::AbstractFloat)

  blob = make_blob(backend, data_type, dims)
  grad = make_blob(backend, data_type, dims)
  rc   = RefCounter(1)

  Parameter(name, blob, grad, init, regu, cons, lr, rc)
end

# make a shared parameter
function share_parameter(backend::Backend, param::Parameter)
  blob  = param.blob
  grad  = make_blob(backend, eltype(blob), size(blob))
  init  = NullInitializer()
  regu  = param.regularizer
  cons  = param.constraint
  lr    = param.learning_rate
  rc    = ref(param.rc)

  Parameter(param.name, blob, grad, init, regu, cons, lr, rc)
end

function destroy(param::Parameter)
  destroy(param.gradient)
  if dec(param.rc) == 0
    destroy(param.blob)
  end
end

type SparseParam <: AbstractParameter
  name          :: AbstractString
  blob          :: Blob
  gradient      :: Blob
  gradient_cols :: Blob
  initializer   :: Initializer
  learning_rate :: AbstractFloat # relative learning rate

  rc            :: RefCounter
end

function make_sparse_param{N}(backend::Backend, name::AbstractString,
    data_type::Type, dims::NTuple{N,Int}, batch_size::Int,
    init::Initializer, lr::AbstractFloat)

  blob = make_blob(backend, data_type, dims)
  grad = make_blob(backend, data_type, dims)
  rc   = RefCounter(1)

  Parameter(name, blob, grad, init, regu, cons, lr, rc)
end

function destroy(param::SparseParam)
  destroy(param.gradient)
  destroy(param.gradient_cols)
  if dec(param.rc) == 0
    destroy(param.blob)
  end
end


text:
"""
    struct Gaussian

# Fields
- `μ`: Mean.
- `σ`: Standard deviation.
"""
struct Gaussian
    μ
    σ
end

sample(d::Gaussian; num_samples::Integer) =
    d.μ .+ d.σ .* randn_gpu(data_eltype(d.μ), size(d.μ)..., num_samples)

Statistics.mean(d::Gaussian) = d.μ
Statistics.std(d::Gaussian) = d.σ

logpdf(d::Gaussian, x) = gaussian_logpdf(x, d.μ, d.σ)
kl(p::Gaussian, q::Gaussian) = kl(p.μ, p.σ, q.μ, q.σ)

Base.map(f, d::Gaussian) = Gaussian(f(d.μ), f(d.σ))

"""
    struct Dirac

The log-pdf of a `Dirac` is defined as `0` everywhere, and the KL divergence between
any two `Dirac`s is also defined as `0`.

# Fields
- `x`: Position.
"""

struct Dirac
    x
end

sample(d::Dirac; num_samples::Integer) = d.x

Statistics.mean(d::Dirac) = d.x
Statistics.std(d::Dirac) = 0

logpdf(d::Dirac, x) = 0
kl(p::Dirac, q::Dirac) = 0

Base.map(f, d::Dirac) = Dirac(f(d.x))

# We need to help `sum` to work with scalars and the keyword argument `dims`.
Base.sum(x::Real; dims=nothing) = x

"""
    struct MultivariateGaussian

# Fields
- `μ`: Mean.
- `Σ`: Covariance matrix.
"""
struct MultivariateGaussian
    μ
    Σ
end

function sample(d::MultivariateGaussian; num_samples::Integer)
    n, num_channels, batch_size = size(d.μ)

    μ = d.μ
    Σ = d.Σ

    samples_channels = []
    for c = 1:num_channels
        samples_batches = []
        for b = 1:batch_size
            ε = cholesky(Σ[:, :, c, b]).U' * randn_gpu(data_eltype(d.μ), n, num_samples)
            push!(samples_batches, μ[:, c:c, b:b] .+ reshape(ε, n, 1, 1, num_samples))
        end
        push!(samples_channels, cat(samples_batches..., dims=3))
    end
    sample = cat(samples_channels..., dims=2)

    return sample
end

Statistics.mean(d::MultivariateGaussian) = d.μ
function Statistics.std(d::MultivariateGaussian)
    Σ = d.Σ
    n = size(Σ, 1)
    return sqrt.(cat([Σ[i, i:i, :, :] for i = 1:n]..., dims=1))
end
Statistics.var(d::MultivariateGaussian) = d.Σ

function logpdf(d::MultivariateGaussian, x)
    _, num_channels, batch_size = size(x)

    μ = d.μ
    Σ = d.Σ

    logpdfs_channels = []
    for c = 1:num_channels
        logpdfs_batches = []
        for b = 1:batch_size
            push!(logpdfs_batches, gaussian_logpdf(x[:, c, b], μ[:, c, b], Σ[:, :, c, b]))
        end
        push!(logpdfs_channels, cat(logpdfs_batches..., dims=3))
    end
    logpdfs = cat(logpdfs_channels..., dims=2)

    return logpdfs
end

Base.map(f, d::MultivariateGaussian) = Gaussian(f(d.μ), f(d.Σ))


text:
#!/usr/local/bin/julia
# author: zhwang
# date: 2022-03-01
# description: 检查待提交文件中的信息，函数名等，并生成函数信息文件 --- 已有函数信息 ---

function get_info(body::AbstractString)
    ## 获取代码内容
    code_reg = r"## @lc code=start\n([\s\S]*)\n## @lc code=end"
    code = match(code_reg, body)
    isnothing(code) && throw("未识别到代码内容")
    code = code.captures[1]
    ## 获取函数内容
    func_reg1 = r"^function (\w+!?\([\w:{}, ]*\))"
    func_reg2 = r"^(\w+!?\([\w:{}, ]*\))[\w: ]*="
    funcs = String[]
    for line in split(code, '\n')
        (func = match(func_reg1, line)) ≢ nothing && push!(funcs, func.captures[1])
        (func = match(func_reg2, line)) ≢ nothing && push!(funcs, func.captures[1])
    end
    ## 获取信息部分
    author = match(r"# author: ([\w ]+)", body).captures[1]
    date = match(r"# date: ([\w-]+)", body).captures[1]
    any(isempty, [date, author]) && throw("未检索到信息！")
    (code, unique!(funcs)), (author, date)
end
### 检查和准备阶段 ###
cd(Base.source_dir()[1:end-5])
current_files = readdir()
!all(file in current_files for file in ["src", "tools"]) && throw("Script file misplacing: could not find source path `src/`")

problem_path = "submit/problems/"
test_path = "submit/test/"
problems, tests = readdir.([problem_path, test_path])
isempty(problems) && throw("问题集为空，请检查！")
problems != tests && throw("问题集和测试集不匹配，请检查")
sort!(problems; by=i -> parse(Int, split(i, '.')[1]))
num = length(problems)

### 获取信息 ###
codes, authors, dates, funcs = String[], String[], String[], Vector{String}[]
for file in problems
    # println("正在处理文件 $file")
    content = read(open(problem_path * file, "r"), String)
    (code, func), (author, date) = get_info(content)
    isempty(func) && throw("未识别到函数：$file")
    push!.([authors, dates, codes], [author, date, code])
    append!(funcs, [func])
end
println("识别完毕，一共 $num 个文件")

### 检查信息 ###
a = findfirst(i->authors[i]!="zhwang", 1:num)
!isnothing(a) && throw("作者信息错误：$(file[a])")
a = findfirst(i->!occursin(r"2022-\d\d-\d\d", dates[i]), 1:num)
!isnothing(a) && throw("日期未修改：$(file[a])")
println("日期和作者名称检查完毕")

### 检查函数信息 ###
# 1. 检查函数派发是否重复信息
# 2. 输出函数名
func_name(func::String) = replace(func, r"(?<name>\w*!?)\(.*\)" => s"\g<name>")
func_method(func::String) = replace(func, r"(?<sym>[(, ]\w+)::[\w{}]+" => s"\g<sym>")

output = "submit/local_functions.md"
io = open(output, "w")
println(io, "# 基本信息\n## 新提交的函数名\n")
println(io, "| 问题编号 | 函数名 |\n| :---: | :---: |")

for (i, file) in enumerate(problems)
    ind = split(file, '.')[1]
    func = unique(func_name.(funcs[i]))
    for f in func
        println(io, "| $ind | `$f` |")
    end
end

println(io, "\n## 函数方法\n")
println(io, "| 问题编号 | 函数方法 |\n| :---: | :---: |")

for (i, file) in enumerate(problems)
    ind = split(file, '.')[1]
    func = unique(func_method.(funcs[i]))
    for f in func
        println(io, "| $ind | `$f` |")
    end
end

close(io)
println("写入完毕，文件保存在 $output 中")

# ## 获取全局信息
# println("检查主仓库的函数信息")
# main_problem_path = "../LeetCode.jl/src/problems/"
# main_problems = readdir(main_problem_path)
# isempty(main_problems) && throw("未能读取到主仓库文件")
# deleteat!(main_problems, length(main_problems))
# sort!(main_problems; by=i -> parse(Int, split(i, '.')[1]))
# total_num = length(main_problems)

# main_funcs = Vector{String}[]
# fails = String[]
# for file in main_problems
#     # println("正在处理文件 $file")
#     content = read(open(main_problem_path * file, "r"), String)
#     (_, func), _ = get_info(content)
#     if isempty(func)
#         println("未识别到函数：$file")
#         push!(fails, file)
#     else
#         append!(main_funcs, [func])
#     end
# end
# println("识别完毕，$(total_num - length(fails)) 个文件成功识别，$(length(fails)) 个文件识别失败")

# for i in fails
#     print(stdout, i)
# end
# for (i, file) in enumerate(main_problems)
#     file in fails && deleteat!(main_problems, i)
# end


text:
using JuMP, Ipopt

@testset "dat export import consistency check" begin
    # min x - y
    # st  x + x^2 + x*y + y^2 <= 1
    #     -2 <= x, y <= 2
    # solution: x+y = -1/3
    # optimal objective -1-4/sqrt(3)

    qcqp = Problem()

    x = MPC.Variable("x", Complex)
    y = MPC.Variable("y", Real)
    b = MPC.Variable("b", Bool)

    set_objective!(qcqp, x - y)

    add_constraint!(qcqp, "ctr1", (x + Base.power_by_squaring(x,2) + x*y + Base.power_by_squaring(y,2)) << 1)
    add_constraint!(qcqp, "ctr2", (Base.power_by_squaring(x,3) + (1.0+3im)*conj(Base.power_by_squaring(x,2))*x*y + Base.power_by_squaring(y,2)) << (1 + 2im))
    add_constraint!(qcqp, "ctr3", -3im << (5*b + 5im*y) << 1)
    add_constraint!(qcqp, "x_bounds", -2 << x << 2)
    add_constraint!(qcqp, "y_bounds", -2 << y << 2)

    point = Point([x], [3])
    dat_exportpath = joinpath(pwd(),"tmp_datexport")
    # dat_exportpath = joinpath(Pkg.dir("MathProgComplex"), "tmp_datexport")
    mkpath(splitdir(dat_exportpath)[1])

    export_to_dat(qcqp, dat_exportpath, filename="POP.dat", point=point)

    qcqp2, point2 = import_from_dat(joinpath(dat_exportpath, "POP.dat"))

    rm(joinpath(dat_exportpath, "POP.dat"))

    @test qcqp.variables == qcqp2.variables
    @test qcqp.objective == qcqp2.objective

    @test Set(keys(qcqp.constraints)) == Set(keys(qcqp2.constraints))

    for (ctrname, ctr) in qcqp.constraints
        @test qcqp.constraints[ctrname].ub == qcqp2.constraints[ctrname].ub
        @test qcqp.constraints[ctrname].lb == qcqp2.constraints[ctrname].lb
        @test qcqp.constraints[ctrname].p == qcqp2.constraints[ctrname].p
    end

    @test point == point2

    rm(dat_exportpath, recursive=true)
end

@testset "WB5.dat import and Ipopt solve" begin
    instancepath = joinpath(pwd(), "test", "instances")
    # instancepath = joinpath(Pkg.dir("MathProgComplex"), "test", "instances")
    WB5_cplx, initpt = import_from_dat(joinpath(instancepath, "WB5.dat"))

    WB5 = pb_cplx2real(WB5_cplx)

    mysolver = Ipopt.Optimizer
    m, variables_jump, ctr_jump, ctr_exp = get_JuMP_cartesian_model(WB5, mysolver)
    optimize!(m)

    sol = get_JuMP_solution(m, variables_jump, WB5)

    @test get_objective(WB5, sol) ≈ 1146.478 atol=1e-3
end


text:
function verified_read_data!(S::GphysData, fmt::String, fpat::Union{String, Array{String,1}};
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  strict  ::Bool    = true,
  swap    ::Bool    = false,              # do byte swap?
  v       ::Int64   = KW.v,               # verbosity level
  vl      ::Bool    = false,              # verbose logging
  allow_empty::Bool = false
  )

  read_data!(S, fmt, fpat,
    full    = full,
    cf      = cf,
    jst     = jst,
    nx_add  = nx_add,
    nx_new  = nx_new,
    strict  = strict,
    swap    = swap,
    v       = v,
    vl      = vl
    )

  basic_checks(S, allow_empty = allow_empty)
  return nothing
end

function verified_read_data(fmt::String, fpat::Union{String, Array{String,1}};
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  strict  ::Bool    = true,
  swap    ::Bool    = false,              # do byte swap?
  v       ::Int64   = KW.v,               # verbosity level
  vl      ::Bool    = false,              # verbose logging
  allow_empty::Bool = false
  )

  S = read_data(fmt, fpat,
    full    = full,
    cf      = cf,
    jst     = jst,
    nx_add  = nx_add,
    nx_new  = nx_new,
    strict  = strict,
    swap    = swap,
    v       = v,
    vl      = vl
    )

  basic_checks(S, allow_empty = allow_empty)
  return S
end


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

function _reshape_v(target::AnyDensityLike, v::Any)
    shape = varshape(convert(AbstractDensity, target))
    reshape_variate(shape, v)
end

function _reshape_vs(target::AnyDensityLike, vs::AbstractVector)
    shape = varshape(convert(AbstractDensity, target))
    reshape_variates(shape, vs)
end

_reshape_rand_output(x::Any) = x
_reshape_rand_output(x::AbstractMatrix) = nestedview(x)

_rand_v(rng::AbstractRNG, src::Distribution) = varshape(src)(rand(rng, bat_sampler(unshaped(src))))
_rand_v(rng::AbstractRNG, src::DistLikeDensity) = varshape(src)(rand(rng, bat_sampler(unshaped(src))))
_rand_v(rng::AbstractRNG, src::Distribution, n::Integer) = _reshape_rand_output(rand(rng, bat_sampler(src), n))
_rand_v(rng::AbstractRNG, src::DistLikeDensity, n::Integer) = _reshape_rand_output(rand(rng, bat_sampler(src), n))

function _rand_v(rng::AbstractRNG, src::AnyIIDSampleable)
    _rand_v(rng, convert(DistLikeDensity, src))
end

function _rand_v(rng::AbstractRNG, src::AnyIIDSampleable, n::Integer)
    _rand_v(rng, convert(DistLikeDensity, src), n)
end

_rand_v(rng::AbstractRNG, src::DensitySampleVector) =
    first(_rand_v(rng, src, 1))

function _rand_v(rng::AbstractRNG, src::DensitySampleVector, n::Integer)
    orig_idxs = eachindex(src)
    weights = FrequencyWeights(float(src.weight))
    resampled_idxs = sample(orig_idxs, weights, n, replace=false, ordered=false)
    samples = src[resampled_idxs]
end

function _rand_v_for_target(rng::AbstractRNG, target::AnySampleable, src::Any)
    _reshape_v(target, _rand_v(rng, src))
end

function _rand_v_for_target(rng::AbstractRNG, target::AnySampleable, src::Any, n::Integer)
    _reshape_vs(target, _rand_v(rng, src, n))
end

_get_first(x::ShapedAsNTArray) = view(x, firstindex(x))
_get_first(x) = first(x)

function _rand_v_for_target(rng::AbstractRNG, target::AnySampleable, src::DensitySampleVector)
    _get_first(_rand_v_for_target(rng, target, src, 1))
end

function _rand_v_for_target(rng::AbstractRNG, target::AnySampleable, src::DensitySampleVector, n::Integer)
    bat_sample(rng, src, RandResampling(nsamples = n)).result.v
end

"""
    struct InitFromTarget <: InitvalAlgorithm

Generates initial values for sampling, optimization, etc. by direct i.i.d.
sampling a suitable component of that target density (e.g. it's prior)
that supports it.

* If the target is supports direct i.i.d. sampling, e.g. because it is a
  distribution, initial values are sampled directly from the target.

* If the target is a posterior density, initial values are sampled from the
  prior (or the prior's prior if the prior is a posterior itself, etc.).

* If the target is a sampled density, initial values are (re-)sampled from
  the available samples.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct InitFromTarget <: InitvalAlgorithm end
export InitFromTarget

function get_initsrc_from_target end

get_initsrc_from_target(target::AnyIIDSampleable) = target
get_initsrc_from_target(target::RenormalizedDensity{<:DistributionDensity}) = bat_sampler(target)

get_initsrc_from_target(target::AbstractPosteriorDensity) = get_initsrc_from_target(getprior(target))

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::InitFromTarget)
    (result = _rand_v_for_target(rng, target, get_initsrc_from_target(target)),)
end

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::InitFromTarget)
    (result = _rand_v_for_target(rng, target, get_initsrc_from_target(target), n),)
end

function bat_initval_impl(rng::AbstractRNG, target::ReshapedDensity, algorithm::InitFromTarget)
    v_orig = bat_initval_impl(rng, parent(target), algorithm).result
    v = varshape(target)(unshaped(v_orig))
    (result = v,)
end

function bat_initval_impl(rng::AbstractRNG, target::ReshapedDensity, n::Integer, algorithm::InitFromTarget)
    v_orig = bat_initval_impl(rng, parent(target), n, algorithm).result
    v = varshape(target).(unshaped.(v_orig))
    (result = v,)
end

function bat_initval_impl(rng::AbstractRNG, target::TransformedDensity, algorithm::InitFromTarget)
    v_orig = bat_initval_impl(rng, target.orig, algorithm).result
    v = target.trafo(v_orig)
    (result = v,)
end

function bat_initval_impl(rng::AbstractRNG, target::TransformedDensity, n::Integer, algorithm::InitFromTarget)
    v = bat_initval_impl(rng, target.orig, n, algorithm).result
    # ToDo: Improve implementation
    trg_vs = valshape(target.trafo(first(v)))
    unshaped_v = unshaped.(v)
    unshaped_v .= unshaped.(target.trafo.(v))
    (result = trg_vs.(unshaped_v),)
end

"""
    struct InitFromSamples <: InitvalAlgorithm

Generates initial values for sampling, optimization, etc. by direct sampling
from a given i.i.d. sampleable source.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct InitFromSamples{SV<:DensitySampleVector} <: InitvalAlgorithm
    samples::SV
end
export InitFromSamples

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::InitFromSamples)
    (result = _rand_v(rng, target, algorithm.samples),)
end

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::InitFromSamples)
    (result = _rand_v(rng, target, algorithm.samples, n),)
end

"""
    struct InitFromIID <: InitvalAlgorithm

Generates initial values for sampling, optimization, etc. by random
resampling from a given set of samples.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct InitFromIID{D<:AnyIIDSampleable} <: InitvalAlgorithm
    src::D
end
export InitFromIID

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::InitFromIID)
    (result = _rand_v(rng, target, algorithm.src),)
end

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::InitFromIID)
    (result = _rand_v(rng, target, algorithm.src, n),)
end

"""
    struct ExplicitInit <: InitvalAlgorithm

Uses initial values from a given vector of one or more values/variates. The
values are used in the order they appear in the vector, not randomly.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
struct ExplicitInit{V<:AbstractVector} <: InitvalAlgorithm
    xs::V
end
export ExplicitInit

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::ExplicitInit)
    (result = first(algorithm.xs),)
end

function bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::ExplicitInit)
    xs = algorithm.xs
    idxs = eachindex(xs)
    (result = xs[idxs[1:n]],)
end

function apply_trafo_to_init(trafo::AbstractVariateTransform, initalg::ExplicitInit)
    xs_tr = trafo.(initalg.xs)
    ExplicitInit(xs_tr)
end


text:
# Types and methods for counting mutations
# ========================================
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

include("site_types.jl")

typealias FourBitAlphs Union{DNAAlphabet{4},RNAAlphabet{4}}
typealias TwoBitAlphs Union{DNAAlphabet{2},RNAAlphabet{2}}
typealias NucAlphs Union{DNAAlphabet, RNAAlphabet}

# Includes for naieve site counting.
include("naive/is_site.jl")
include("naive/count_sites_naive.jl")

"""
    count_sites{T<:Site}(::Type{T}, a::BioSequence, b::BioSequence)

Mutations are counted using the `count_sites` method.

If the concrete site case type inherits from NoPairDel, then the result will
be a simple integer value of the number of counts of site.

If the concrete site case type inherits from PairDel, then the result
will be a tuple. As a result, counting methods for such types also count and
return the number of sites which could not be unambiguously determined.
This second count is important for some downstream purposes, for example
evolutionary/genetic distance computations in which pairwise deletion of
ambiguous sites is necessary.
"""
function count_sites{T<:Site}(::Type{T}, a::BioSequence, b::BioSequence)
    return count_sites_naive(T, a, b)
end

"""
    count_sites{T<:NoPairDel}(::Type{T}, a::EachWindowIterator, b::EachWindowIterator)

A method of the count_sites function for counting sites of SiteType `T` between two aligned
sequences.

Crucially, this method accepts two EachWindowIterators as parameters `a` and `b`
instead of two `BioSequence`s. As a result, the two sequences are iterated over
with a sliding window, and a site count computed for each window.

This function returns two vectors, one with the ranges
"""
function count_sites{T<:Site}(::Type{T}, a::EachWindowIterator, b::EachWindowIterator)
    itr = zip(a, b)
    results = Vector{IntervalValue{Int, Int}}(length(itr))
    i = 1
    for (wina, winb) in itr
        interval = wina[1]
        results[i] = IntervalValue(first(interval), last(interval), count_sites(T, wina[2], winb[2]))
        i += 1
    end
    return results
end

function count_sites{T<:Mutation}(::Type{T}, a::EachWindowIterator, b::EachWindowIterator)
    itr = zip(a, b)
    len = length(itr)
    results = Vector{IntervalValue{Int}}(len)
    undetermined = Vector{IntervalValue{Int}}(len)
    i = 1
    for (wina, winb) in itr
        interval = wina[1]
        counts = count_sites(T, wina[2], winb[2])
        results[i] = IntervalValue(first(interval), last(interval), counts[1])
        undetermined[i] = IntervalValue(first(interval), last(interval), counts[2])
        i += 1
    end
    return results, undetermined
end

"Count the number of sites between two nucleotide sequences using a sliding window."
function count_sites{T<:Site}(::Type{T}, a::BioSequence, b::BioSequence, width::Int, step::Int)
    awindows = eachwindow(a, width, step)
    bwindows = eachwindow(b, width, step)
    return count_sites(T, awindows, bwindows)
end

"""
Count the number of mutations between nucleotide sequences in a pairwise manner.
"""
function count_sites{T<:Site,A<:Alphabet}(::Type{T}, sequences::Array{BioSequence{A}})
    len = length(sequences)
    counts = PairwiseListMatrix(Int, len, false)
    @inbounds for i in 1:len, j in (i + 1):len
        counts[i, j] = count_sites(T, sequences[i], sequences[j])
    end
    return counts
end

"""
Count the number of mutations between nucleotide sequences in a pairwise manner,
using a sliding window of a given width and step.
"""
function count_sites{T<:Site,A<:Alphabet}(::Type{T}, sequences::Array{BioSequence{A}}, width::Int, step::Int)
    len = length(sequences)
    counts = PairwiseListMatrix(Vector{IntervalValue{Int}}, len, false)
    @inbounds for i in 1:len, j in (i + 1):len
        counts[i, j] = count_sites(T, sequences[i], sequences[j], width, step)
    end
    return counts
end

"""
Count the number of mutations between nucleotide sequences in a pairwise manner.
"""
function count_sites{T<:Mutation,A<:Alphabet}(::Type{T}, sequences::Array{BioSequence{A}})
    len = length(sequences)
    counts = PairwiseListMatrix(Int, len, false)
    undetermined = PairwiseListMatrix(Int, len, false)
    @inbounds for i in 1:len, j in (i + 1):len
        counts[i, j], undetermined[i, j] = count_sites(T, sequences[i], sequences[j])
    end
    return counts, undetermined
end

"""
Count the number of mutations between nucleotide sequences in a pairwise manner,
using a sliding window of a given width and step.
"""
function count_sites{T<:Mutation,A<:Alphabet}(::Type{T}, sequences::Array{BioSequence{A}}, width::Int, step::Int)
    len = length(sequences)
    counts = PairwiseListMatrix(Vector{IntervalValue{Int}}, len, false)
    undetermined = PairwiseListMatrix(Vector{IntervalValue{Int}}, len, false)
    @inbounds for i in 1:len, j in (i + 1):len
        counts[i, j], undetermined[i, j] = count_sites(T, sequences[i], sequences[j], width, step)
    end
    return counts, undetermined
end

"Count the number of sites of a given SiteType `T` in a biological sequence."
function count_sites{T<:Site}(::Type{T}, seq::BioSequence)
    return count_sites_naive(T, seq)
end


text:
module CommTest

using Comm
using Base.Test

function remote1(i, s, a)
    @test Comm.myproc() == (Comm.nprocs() < 2 ? 1 : 2)
    @test i == 1
    @test s == "a"
    @test a == [1.0]
    rexec(remote2, Comm.nprocs())
end

function remote2()
    @test Comm.myproc() == Comm.nprocs()
    rexec(1) do
        remote3()
    end
end

function remote3()
    @test Comm.myproc() == 1
    global DONE = true
end

function inc()
    # rexec(() -> global COUNTER += 1, 1)
    rexec(inc1, 1)
end
function inc1()
    global COUNTER += 1
end

function main()
    # local function
    let
        @test Comm.myproc() == 1
        global DONE = false
        rexec(mod1(2, Comm.nprocs())) do
            remote1(1, "a", [1.0])
        end
        while !DONE yield() end
    end

    let
        @test Comm.myproc() == 1
        global DONE = false
        i=1
        rexec(()->(s="a"; remote1(i, s, [1.0])), mod1(2, Comm.nprocs()))
        while !DONE yield() end
    end

    let
        @test Comm.myproc() == 1
        global DONE = false
        i=1
        s="a"
        f() = remote1(i, s, [1.0])
        rexec(f, mod1(2, Comm.nprocs()))
        while !DONE yield() end
    end

    let
        global COUNTER = 0
        rexec_everywhere(inc)
        while COUNTER < Comm.nprocs() yield() end
    end

    # # local macro
    # let
    #     while !DONE yield() end
    #     global DONE = false
    #     @rexec mod1(2, Comm.nprocs()) remote1(1, "a", [1.0])
    #     while COUNTER < Comm.nprocs() yield() end
    # end
    #
    # let
    #     while !DONE yield() end
    #     global DONE = false
    #     i=1
    #     @rexec mod1(2, Comm.nprocs()) (s="a"; remote1(i, s, [1.0]))
    #     while COUNTER < Comm.nprocs() yield() end
    # end
    #
    # let
    #     global COUNTER = 0
    #     @rexec_everywhere inc()
    #     while COUNTER < Comm.nprocs() yield() end
    # end
end

end


text:
# ================================
# Hybrid flowpipe
# ================================

"""
    HybridFlowpipe{N, D, FT<:AbstractFlowpipe, VOA<:VectorOfArray{N, D, Vector{FT}}} <: AbstractFlowpipe

Type that wraps a vector of flowpipes of the same type, such that they are
contiguous in time.

### Fields

- `Fk`  -- vector of flowpipes
- `ext` -- (optional, default: empty) dictionary for extensions

### Notes

The evaluation functions (in time) for this type do not assume that the flowpipes are contiguous in time.
That is, the final time of the `i`-th flowpipe does not match the start time of the `i+1`-th flowpipe.
"""
struct HybridFlowpipe{N, RT<:AbstractReachSet{N}, FT<:AbstractFlowpipe} <: AbstractFlowpipe
    Fk::VectorOfArray{RT, 2, Vector{FT}}
    ext::Dict{Symbol, Any}
end

function HybridFlowpipe(Fk::Vector{FT}) where {N, RT<:AbstractReachSet{N}, FT<:Flowpipe{N, RT}}
    voa = VectorOfArray{RT, 2, Vector{FT}}(Fk)
    ext = Dict{Symbol, Any}()
    return HybridFlowpipe{N, RT, FT}(voa, ext)
end

function HybridFlowpipe(Fk::Vector{FT}, ext::Dict{Symbol, Any}) where {N, RT<:AbstractReachSet{N}, FT<:Flowpipe{N, RT}}
    voa = VectorOfArray{RT, 2, Vector{FT}}(Fk)
    return HybridFlowpipe{N, RT, FT}(voa, ext)
end

function HybridFlowpipe(Fk::Vector{SFT}) where {N, RT, FT<:Flowpipe{N, RT}, NT<:Number, SFT<:ShiftedFlowpipe{FT, NT}}
    voa = VectorOfArray{RT, 2, Vector{SFT}}(Fk)
    ext = Dict{Symbol, Any}()
    return HybridFlowpipe{N, RT, SFT}(voa, ext)
end

# interface functions
array(fp::HybridFlowpipe) = fp.Fk
flowpipe(fp::HybridFlowpipe) = fp
numtype(::HybridFlowpipe{N}) where {N} = N
setrep(::Type{HybridFlowpipe{N, RT, FT}}) where {N, RT, FT} = setrep(RT)
setrep(::HybridFlowpipe{N, RT, FT}) where {N, RT, FT} = setrep(RT)
rsetrep(::Type{HybridFlowpipe{N, RT, FT}}) where {N, RT, FT} = RT
numrsets(fp::HybridFlowpipe) = mapreduce(length, +, fp)

# indexing: fp[j, i] returning the j-th reach-set of the i-th flowpipe
Base.getindex(fp::HybridFlowpipe, I::Int...) = getindex(fp.Fk, I...)

function tspan(fp::HybridFlowpipe)
    ti = minimum(tstart, fp)
    tf = maximum(tend, fp)
    return TimeInterval(ti, tf)
end

function Base.similar(fp::HybridFlowpipe{N, RT, FT}) where {N, RT, FT}
    return HybridFlowpipe(Vector{FT}())
end

function overapproximate(fp::HybridFlowpipe, args...)
    return HybridFlowpipe([overapproximate(F, args...) for F in fp], fp.ext)
end

function project(fp::HybridFlowpipe, args...)
    return [project(F, args...) for F in fp]
end

# LazySets interface
function LazySets.ρ(d::AbstractVector, fp::HybridFlowpipe)
    return maximum(ρ(d, F) for F in array(fp))
end

#function LazySets.σ(d::AbstractVector, fp::HybridFlowpipe)
#    error("not implemented")
#end

Base.:⊆(F::HybridFlowpipe, X::LazySet) = all(fp ⊆ X for fp in F)
Base.:⊆(F::HybridFlowpipe, Y::AbstractLazyReachSet) = all(fp ⊆ set(Y) for fp in F)

# evaluation for scalars
function (fp::HybridFlowpipe{N, RT})(t::Number) where {N, RT<:AbstractReachSet{N}}
    if t ∉ tspan(fp)
        throw(ArgumentError("time $t does not belong to the time span, " *
                            "$(tspan(fp)), of the given flowpipe"))
    end
    vec = Vector{RT}()
    for (k, fk) in enumerate(array(fp)) # loop over flowpipes
        if t ∈ tspan(fk)
            for Ri in fk  # loop over reach-sets for this flowpipe
                if t ∈ tspan(Ri)
                    push!(vec, Ri)
                end
            end
        end
    end
    return vec
end

# evaluation for time intervals
function (fp::HybridFlowpipe{N, RT})(dt::TimeInterval) where {N, RT<:AbstractReachSet{N}}
    if !(dt ⊆ tspan(fp)) # ⊊ doesn't work see, IntervalArithmetic#409
        throw(ArgumentError("time interval $dt does not belong to the time span, " *
                            "$(tspan(fp)), of the given flowpipe"))
    end
    vec = Vector{RT}()
    for (k, fk) in enumerate(array(fp)) # loop over flowpipes
        dtint = dt ∩ tspan(fk)
        if !isempty(dtint)
            for R in fk(dtint)  # loop over reach-sets for this flowpipe
                push!(vec, R)
            end
        end
    end
    return vec
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Sundials, ParameterizedFunctions, Plots, ODE, ODEInterfaceDiffEq, ODEInterface, LSODA
gr()
using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

van = @ode_def begin
  dy = μ*((1-x^2)*y - x)
  dx = 1*y
end μ

prob = ODEProblem(van,[0;2.],(0.0,6.3),1e6)
abstols = 1.0 ./ 10.0 .^ (5:9)
reltols = 1.0 ./ 10.0 .^ (2:6)

sol = solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

plot(sol,ylim=[-4;4])

plot(sol)

#sol = solve(prob,ode23s()); println("Total ODE.jl steps: $(length(sol))")
#using GeometricIntegratorsDiffEq
#try
#    sol = solve(prob,GIRadIIA3(),dt=1/1000)
#catch e
#    println(e)
#end

sol = solve(prob,ARKODE(),abstol=1e-4,reltol=1e-2);

sol = solve(prob,ARKODE(nonlinear_convergence_coefficient = 1e-6),abstol=1e-4,reltol=1e-1);

sol = solve(prob,ARKODE(order=3),abstol=1e-4,reltol=1e-1);

sol = solve(prob,ARKODE(nonlinear_convergence_coefficient = 1e-6,order=3),abstol=1e-4,reltol=1e-1);

sol = solve(prob,ARKODE(order=5,nonlinear_convergence_coefficient = 1e-3),abstol=1e-4,reltol=1e-1);

sol = solve(prob,ARKODE(order=5,nonlinear_convergence_coefficient = 1e-4),abstol=1e-4,reltol=1e-1);

setups = [
          #Dict(:alg=>ROCK2())    #Unstable
          #Dict(:alg=>ROCK4())    #needs more iterations
          #Dict(:alg=>ESERK5()),
          ]

setups = [
  #Dict(:alg=>Hairer4()),
  #Dict(:alg=>Hairer42()),
  #Dict(:alg=>Cash4()),
]

sol = solve(prob,EXPRB53s3(),dt=2.0^(-8));
sol = solve(prob,EPIRK4s3B(),dt=2.0^(-8));
sol = solve(prob,EPIRK5P2(),dt=2.0^(-8));

abstols = 1.0 ./ 10.0 .^ (4:7)
reltols = 1.0 ./ 10.0 .^ (1:4)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ddebdf()),
          Dict(:alg=>rodas()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas3()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>rodas()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>radau()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>ROS34PW1a()),
          ]
gr()
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Kvaerno3()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(nonlinear_convergence_coefficient = 1e-6)),
          Dict(:alg=>SDIRK2()),
          Dict(:alg=>radau())]
names = ["Rosenbrock23" "Kvaerno3" "KenCarp4" "TRBDF2" "KenCarp3" "ARKODE" "SDIRK2" "radau"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5,nonlinear_convergence_coefficient = 1e-4)),
          Dict(:alg=>ARKODE(nonlinear_convergence_coefficient = 1e-6)),
          Dict(:alg=>ARKODE(nonlinear_convergence_coefficient = 1e-6,order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ImplicitEulerExtrapolation()),
          #Dict(:alg=>ImplicitDeuflhardExtrapolation()), # Diverges
          #Dict(:alg=>ImplicitHairerWannerExtrapolation()), # Diverges
          Dict(:alg=>ABDF2()),
          #Dict(:alg=>QNDF()), # ???
          #Dict(:alg=>Exprb43()), # Diverges
          Dict(:alg=>Exprb32()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (4:7)
reltols = 1.0 ./ 10.0 .^ (1:4)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ddebdf()),
          Dict(:alg=>rodas()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      error_estimator=:l2,appxsol=test_sol,maxiters=Int(1e5),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas3()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>rodas()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>radau()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>ROS34PW1a()),
          ]
gr()
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimator=:l2,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23(),:dense=>false),
          Dict(:alg=>Kvaerno3(),:dense=>false),
          Dict(:alg=>KenCarp4(),:dense=>false),
          Dict(:alg=>TRBDF2(),:dense=>false),
          Dict(:alg=>KenCarp3(),:dense=>false),
          Dict(:alg=>SDIRK2(),:dense=>false),
          Dict(:alg=>radau())]
names = ["Rosenbrock23" "Kvaerno3" "KenCarp4" "TRBDF2" "KenCarp3" "SDIRK2" "radau"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,appxsol=test_sol,maxiters=Int(1e5),error_estimator=:l2,seconds=5)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ImplicitEulerExtrapolation()),
          #Dict(:alg=>ImplicitDeuflhardExtrapolation()), # Diverges
          #Dict(:alg=>ImplicitHairerWannerExtrapolation()), # Diverges
          Dict(:alg=>ABDF2()),
          #Dict(:alg=>QNDF()), # ???
          #Dict(:alg=>Exprb43()), # Diverges
          Dict(:alg=>Exprb32()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimator=:l2,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:11)
reltols = 1.0 ./ 10.0 .^ (4:8)
setups = [Dict(:alg=>Rodas3()),
          Dict(:alg=>GRK4A()),
          Dict(:alg=>Rodas4P()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>Rodas5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e6),seconds=5)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:11)
reltols = 1.0 ./ 10.0 .^ (4:8)
setups = [Dict(:alg=>Rodas3()),
          Dict(:alg=>Kvaerno4()),
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>ARKODE()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau()),
          Dict(:alg=>Rodas5())]
names = ["Rodas3" "Kvaerno4" "Kvaerno5" "CVODE_BDF" "KenCarp4" "KenCarp5" "ARKODE" "Rodas4" "radau" "Rodas5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e6),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rodas3()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau()),
          Dict(:alg=>Rodas5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e6),seconds=5)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:11)
reltols = 1.0 ./ 10.0 .^ (4:8)
setups = [Dict(:alg=>Rodas3()),
          Dict(:alg=>GRK4A()),
          Dict(:alg=>Rodas4P()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>Rodas5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:l2,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e6),seconds=5)
plot(wp)

setups = [Dict(:alg=>Rodas3()),
          Dict(:alg=>Kvaerno4()),
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau()),
          Dict(:alg=>Rodas5())]
names = ["Rodas3" "Kvaerno4" "Kvaerno5" "CVODE_BDF" "KenCarp4" "KenCarp5" "Rodas4" "radau" "Rodas5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,appxsol=test_sol,maxiters=Int(1e6),error_estimate=:l2,seconds=5)
plot(wp)

setups = [Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau()),
          Dict(:alg=>Rodas5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e6),error_estimate=:l2,seconds=5)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


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

"""
    Chain(p1, p2, ..., pn)

A polygonal chain from a sequence of points `p1`, `p2`, ..., `pn`.
See https://en.wikipedia.org/wiki/Polygonal_chain.
"""
struct Chain{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polytope{1,Dim,T}
  vertices::V
end

# specialize Polytope outer constructor to use standard vector
Chain(vertices::Vararg{<:Point}) = Chain(collect(vertices))

Chain(vertices::CircularVector) =
  Chain([collect(vertices); vertices[begin]])

nvertices(c::Chain) = length(c.vertices) - isclosed(c)

"""
    vertices(chain)

Return the vertices of the `chain`. In the case that the
chain is closed, the returned vector is circular and can
be indexed with arbitrarily negative or positive indices.
"""
function vertices(c::Chain)
  if isclosed(c)
    vs = @view c.vertices[begin:end-1]
    CircularVector(vs)
  else
    c.vertices
  end
end

"""
    segments(chain)

Return the segments linking consecutive points of the `chain`.
"""
function segments(c::Chain)
  vs = c.vertices
  n = length(vs)
  (Segment(view(vs, [i,i+1])) for i in 1:n-1)
end

"""
    isclosed(chain)

Tells whether or not the chain is closed.

A closed chain is also known as a ring.
"""
isclosed(c::Chain) = first(c.vertices) == last(c.vertices)

"""
   issimple(chain)

Tells whether or not the `chain` is simple.

A chain is simple when all its segments only
intersect at end points.
"""
function issimple(c::Chain)
  vs = c.vertices
  ss = collect(segments(c))
  for i in 1:length(ss)
    for j in i+1:length(ss)
      I = intersecttype(ss[i], ss[j])
      if !(I isa CornerTouchingSegments || I isa NonIntersectingSegments)
        return false
      end
    end
  end
  true
end

"""
    windingnumber(point, chain)

Winding number of `point` with respect to the `chain`.

## References

* Balbes, R. and Siegel, J. 1990. [A robust method for calculating
  the simplicity and orientation of planar polygons]
  (https://www.sciencedirect.com/science/article/abs/pii/0167839691900198)
"""
function windingnumber(p::Point, c::Chain)
  vₒ, vs = p, c.vertices
  sum(∠(vs[i], vₒ, vs[i+1]) for i in 1:length(vs)-1)
end

abstract type OrientationMethod end

struct WindingOrientation <: OrientationMethod end

struct TriangleOrientation <: OrientationMethod end

"""
    orientation(chain, [method])

Returns the orientation of the `chain` as either
counter-clockwise (CCW) or clockwise (CW).

Optionally, specify the orientation `method`:

* `WindingOrientation()` - Balbes, R. and Siegel, J. 1990.
* `TriangleOrientation()` - Held, M. 1998.

Default method is `WindingOrientation()`.

## References

* Balbes, R. and Siegel, J. 1990. [A robust method for calculating
  the simplicity and orientation of planar polygons]
  (https://www.sciencedirect.com/science/article/abs/pii/0167839691900198)
* Held, M. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
"""
orientation(c::Chain) = orientation(c, WindingOrientation())

function orientation(c::Chain{Dim,T}, ::WindingOrientation) where {Dim,T}
  # pick any segment
  x1, x2 = c.vertices[1:2]
  x̄ = centroid(Segment(x1, x2))
  w = windingnumber(x̄, c) - ∠(x1, x̄, x2)
  isapprox(w, T(π), atol=atol(T)) ? :CCW : :CW
end

function orientation(c::Chain{Dim,T}, ::TriangleOrientation) where {Dim,T}
  v = vertices(c)
  Δ(i) = signarea(v[1], v[i], v[i+1])
  a = mapreduce(Δ, +, 2:length(v)-1)
  a ≥ zero(T) ? :CCW : :CW
end

"""
    unique(chain)

Return a new `chain` without duplicate vertices.
"""
Base.unique(c::Chain) = unique!(deepcopy(c))

"""
    unique!(chain)

Remove duplicate vertices in the `chain`.
"""
function Base.unique!(c::Chain)
  # sort vertices lexicographically
  verts = c.vertices
  perms = sortperm(coordinates.(verts))

  # remove true duplicates
  keep = Int[]
  sorted = @view verts[perms]
  for i in 1:length(sorted)-1
    if sorted[i] != sorted[i+1]
      # save index in the original vector
      push!(keep, perms[i])
    end
  end
  push!(keep, last(perms))

  # preserve chain order
  sort!(keep)

  # update vertices in place
  copy!(verts, verts[keep])

  c
end

"""
    close!(chain)

Close the `chain`, i.e. repeat the first vertex
at the end of the vertex list.
"""
function close!(c::Chain)
  push!(c.vertices, first(c.vertices))
  c
end

"""
    open!(chain)

Open the `chain`, i.e. remove the last vertex.
"""
function open!(c::Chain)
  pop!(c.vertices)
  c
end

"""
    reverse!(chain)

Reverse the `chain` vertices in place.
"""
function Base.reverse!(c::Chain)
  reverse!(c.vertices)
  c
end

"""
    reverse(chain)

Reverse the `chain` vertices.
"""
Base.reverse(c::Chain) = reverse!(deepcopy(c))

"""
    angles(chain)

Return angles `∠(vᵢ-₁, vᵢ, vᵢ+₁)` at all vertices
`vᵢ` of the `chain`. If the chain is open, the first
and last vertices have no angles. Positive angles
represent a CCW rotation whereas negative angles
represent a CW rotation. In either case, the
absolute value of the angles returned is never
greater than `π`.
"""
function angles(c::Chain)
  θs = map(2:length(c.vertices)-1) do i
    p1 = c.vertices[i-1]
    p2 = c.vertices[i  ]
    p3 = c.vertices[i+1]
    ∠(p1, p2, p3)
  end

  if isclosed(c)
    p1 = c.vertices[end-1]
    p2 = c.vertices[1]
    p3 = c.vertices[2]
    pushfirst!(θs, ∠(p1, p2, p3))
  end

  θs
end

"""
    innerangles(chain)

Return inner angles of the *closed* `chain`. Inner
angles are always positive, and unlike `angles`
they can be greater than `π`.
"""
function innerangles(c::Chain{Dim,T}) where {Dim,T}
  @assert isclosed(c) "Inner angles only defined for closed chains"

  # correct sign of angles in case orientation is CW
  θs = orientation(c) == :CW ? -angles(c) : angles(c)

  [θ > 0 ? 2*T(π) - θ : -θ for θ in θs]
end

Base.view(c::Chain, inds) = Chain(view(vertices(c), inds))

function Base.show(io::IO, c::Chain)
  N = length(c.vertices)
  print(io, "$N-chain")
end

function Base.show(io::IO, ::MIME"text/plain", c::Chain{Dim,T}) where {Dim,T}
  v = c.vertices
  N = length(v)
  I, J = N > 10 ? (5, N-4) : (N, N+1)
  lines = [["  └─$(v[i])" for i in 1:I]
           (N > 10 ? ["  ⋮"] : [])
           ["  └─$(v[i])" for i in J:N]]
  println(io, "$N-chain{$Dim,$T}")
  print(io, join(lines, "\n"))
end


text:
using Documenter, TextAnalysis

load_dir(x) = map(file -> joinpath("lib", x, file), readdir(joinpath(Base.source_dir(), "src", "lib", x)))

makedocs(
   modules = [TextAnalysis],
   clean = false,
   format = [:html],#, :latex],
   sitename = "TextAnalysis",
   pages = Any[
       "Home" => "index.md",
       "Documents" => "documents.md",
       "Corpus" => "corpus.md",
       "Features" => "features.md",
       "Semantic Analysis" => "semantic.md"
   ],
   assets = ["assets/custom.css", "assets/custom.js"]
)

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