text:
export HPolygon, add_point!, npoints, RandomHPolygon, sides
export polygon_check

"""
`HPolygon()` creates a new polygon (with no points).

`HPolygon(list)` creates a polygon whose points are specified in `list`.

See: `add_point!`
"""
struct HPolygon <: HObject
    plist::Array{HPoint,1}
    attr::Dict{Symbol,Any}
    function HPolygon()
        X = new(HPoint[],Dict{Symbol,Any}())
        set_color(X)
        set_thickness(X)
        set_line_style(X)
        return X
    end
end

"""
`add_point!(X::HPolygon, P::HPoint)` adds the point `P`
as the last point of the polygon `X`
"""
add_point!(X::HPolygon, P::HPoint) = push!(X.plist,P)

"""
`endpoints(X::HPolygon)` returns the list of vertices (in order)
of the polygon.
"""
endpoints(X::HPolygon) = deepcopy(X.plist)

"""
`sides(P:::HPolygon/HTriangle)` returns a list of the line segments
that are the sides of the polygon.
"""
function sides(X::HPolygon)::Array{HSegment,1}
    n = npoints(X)
    if n < 2
        return HSegment[]
    end

    result = Array{HSegment,1}(undef,n)
    for k=1:n-1
        a = X.plist[k]
        b = X.plist[k+1]
        result[k] = a+b
    end
    result[end] = X.plist[end] + X.plist[1]
    return result
end

sides(T::HTriangle) = sides(HPolygon(T))

"""
`npoints(X::HPolygon)` returns the number of points on the polygon.
"""
npoints(X::HPolygon) = length(X.plist)

function HPolygon(pts::Array{HPoint,1})
    X = HPolygon()
    for p in pts
        add_point!(X,p)
    end
    return X
end

HPolygon(pts...) = HPolygon(collect(pts))
HPolygon(X::HPolygon) = HPolygon(X.plist)  # copy constructor

"""
`polygon_check(X::HPolygon,quiet=true)` checks that the polygon is nondegenerate.
Possible degeneracies are:
+ Repeated vertices
+ Fewer than three distinct vertices
+ Angles that are either 0 degrees or 180 degrees
If `quiet` is `false`, then a reason for the failed check is printed.
"""
function polygon_check(X::HPolygon, quiet::Bool=true)::Bool

    # Check that the endpoints are all distinct
    n = npoints(X)

    C = HContainer(X.plist...)
    if length(C) != n
        quiet || println("The polygon has repeated vertices")
        return false
    end

    if length(C) < 3
        quiet || println("The polygon is degenerate (fewer than 3 distinct vertices)")
        return false
    end

    angs = angles(X)
    zero_angs = angs .< (THRESHOLD * eps(1.0))
    if any(zero_angs)
        quiet || println("The polygon has 0-degree angles")
        return false
    end

    big_angs = angs .> (pi - THRESHOLD*eps(1.0))
    if any(big_angs)
        quiet || println("The polygon has 180-degree angles")
        return false
    end

    return true
end

function HPolygon(T::HTriangle)
    a,b,c = endpoints(T)
    return HPolygon(a,b,c)
end

function HTriangle(X::HPolygon)
    @assert npoints(X)==3 "Can only convert a 3-point HPolygon into an HTriangle"
    return HTriangle(X.plist...)
end

"""
`RandomHPolygon(n::Int,simple::Bool=false)` create a new `HPolygon` with
`n` points chosen at random. With `simple` set to `true`, return a polygon
that does not self-intersect.
"""
function RandomHPolygon(n::Int, simple::Bool=false)
    @assert n>=0 "Number of vertices must be nonnegative"
    if simple && n>3
        P = RandomHPolygon(n)
        while !is_simple(P)
            P = RandomHPolygon(n)
        end
        return P
    end

    pts = [ RandomHPoint() for j=1:n ]
    return HPolygon(pts)
end

"""
`angles(P::HPolygon)` returns a list of the angles at the vertices of `P`.

+ The results are always in the interval `[0,pi]`.
+ The order of the angles is the order of the vertices in `P.plist`.
"""
function angles(P::HPolygon)::Array{Float64,1}
    n = npoints(P)
    result = zeros(Float64,n)
    if n < 3
        return result
    end

    # first angle
    result[1] = angle(P.plist[end],P.plist[1],P.plist[2])

    for j=2:n-1
        result[j] = angle(P.plist[j-1], P.plist[j], P.plist[j+1])
    end

    result[n] = angle(P.plist[n-1],P.plist[n],P.plist[1])

    return result
end

function perimeter(P::HPolygon)
    n = npoints(P)
    if n < 2
        return 0.0
    end
    result = dist(P.plist[1],P.plist[end])
    if n==2
        return 2*result
    end
    for j=1:n-1
        result += dist(P.plist[j],P.plist[j+1])
    end
    return result
end

function show(io::IO,X::HPolygon)
    print(io,"HPolygon with $(npoints(X)) points")
end

# require 0 <= k < n
"""
`_cycle(A,k)` returns a `k`-step shift of `A`. We require
`k` to be in the interval `[0,n-1]` where `n=length(A)`.
No checking is done.
"""
function _cycle(A::Array{T,1}, k::Int) where T
    n = length(A)
    B = Array{T,1}(undef,n)
    for t = 1:n-k
        @inbounds B[t] = A[t+k]
    end
    for t = 1:k
        @inbounds B[n-k+t] = A[t]
    end
    return B
end

"""
`_cyclic_equal(A,B)` checks if some cyclic shift of one list
equals the other.
"""
function _cyclic_equal(A::Array{S,1}, B::Array{T,1}) where {S,T}
    n = length(A)
    if length(B) != n
        return false
    end
    for s=0:n-1
        if A == _cycle(B,s)
            return true
        end
    end
    return false
end

(==)(X::HPolygon, Y::HPolygon) = _cyclic_equal(X.plist, Y.plist) ||
    _cyclic_equal(X.plist,reverse(Y.plist))

"""
`is_simple(X::HPolygon)` determines if the polygon edges do not self-intersect.
Be sure the polygon is legit using `polygon_check` first.
"""
function is_simple(X::HPolygon)
    n = npoints(X)
    if n <= 3
        return true   # triangles are fine
    end
    slist = sides(X)
    for i=1:n-2
        for j=i+2:n
            if i==1 && j==n
                continue
            end
            if meet_check(slist[i], slist[j])
                return false
            end
        end
    end
    return true
end


text:
module ExampleHouseholdGetter
#
# This module contains code to fetch some test households from CSV files.
# It differs from the main FRSHouseholdGetter in that it doesn't bother with weighting, and households
# to be accessed by names rather than indexes.
#
using DataFrames
using CSV

using ScottishTaxBenefitModel
using .Definitions
using .ModelHousehold: Household
using .HouseholdFromFrame: load_hhld_from_frame
using .RunSettings: Settings

export  initialise, get_household

EXAMPLE_HOUSEHOLDS = Dict{String,Household}()

KEYMAP = Vector{AbstractString}()

"""
return number of households available
"""
function initialise(
    settings       :: Settings
    ;
    # fixme move these to settings
    household_name :: AbstractString = "example_households",
    people_name    :: AbstractString = "example_people",
     ) :: Vector{AbstractString}

    global KEYMAP
    global EXAMPLE_HOUSEHOLDS
    KEYMAP = Vector{AbstractString}()
    hh_dataset = CSV.File("$(MODEL_DATA_DIR)/$(household_name).tab", delim='\t' ) |> DataFrame
    people_dataset = CSV.File("$(MODEL_DATA_DIR)/$(people_name).tab", delim='\t' ) |> DataFrame
    npeople = size( people_dataset)[1]
    nhhlds = size( hh_dataset )[1]
    for hseq in 1:nhhlds
        hhf = hh_dataset[hseq,:]
        push!( KEYMAP, hhf.name )
        println( "loading $(hhf.name)")
        EXAMPLE_HOUSEHOLDS[hhf.name] = load_hhld_from_frame(
            hseq, hhf, people_dataset, ExampleSource, settings )
    end
    return KEYMAP
end

function example_names()
    return KEYMAP
end

function get_household( pos :: Integer ) :: Household
    key = KEYMAP[pos]
    return EXAMPLE_HOUSEHOLDS[key]
end

function get_household( name :: AbstractString ) :: Household
    # global EXAMPLE_HOUSEHOLDS
    return EXAMPLE_HOUSEHOLDS[name]
end

end

ext: jl

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
"""
```julia
    Conjugate( ϵ::T, scaled::Bool = true ) where {T <: Number}
```

Terminate the Fast Gradient Method when the conjugate residual value is below `ϵ`.
The value of `ϵ` is divided by the number of variables in the problem when `scaled`
is true, producing a scaled termination criteria.

This criteria is based on the results in §6.4 of [1].

[1] S. Richter, ‘Computational complexity certification of gradient methods for real-time model predictive control’, Doctoral Thesis, ETH Zurich, 2012.

See also: [`Gradient`](@ref) and [`Best`](@ref)
"""
mutable struct Conjugate{T} <: AbstractStopCondition where {T <: Number}
    ϵ::T

    scale::Bool   # Should a scaled termination criteria be used
    n::Int        # Number of variables in the problem
    ϵₛ::T         # Scaled termination criteria

    function Conjugate( ϵ::T; scaled::Bool = true ) where {T <: Number}
        return new{T}( ϵ, scaled, 1, ϵ )
    end
end

function init!( cond::Conjugate, n, L, μ )
    cond.n = n

    if( cond.scale )
        cond.ϵₛ = cond.ϵ * 1/n
    else
        cond.ϵₛ = cond.ϵ
    end
end

function (cond::Conjugate)( state::FGMState )
    return compute_conjugate( state ) < cond.ϵₛ
end

function compute_conjugate( state::FGMState )
    return abs( state.xₙ'*state.∇xₙ + norm( state.∇xₙ, 1 ) )
end

"""
```julia
    Gradient( ϵ::T, scaled::Bool = true ) where {T <: Number}
```

Terminate the Fast Gradient Method when the gradient value is below `ϵ`.
The value of `ϵ` is divided by the number of variables in the problem when `scaled`
is true, producing a scaled termination criteria.

This criteria is based on the results in §6.4 of [1].

[1] S. Richter, ‘Computational complexity certification of gradient methods for real-time model predictive control’, Doctoral Thesis, ETH Zurich, 2012.

See also: [`Conjugate`](@ref) and [`Best`](@ref)
"""
mutable struct Gradient{T} <: AbstractStopCondition where {T <: Number}
    ϵ::T
    L::T
    μ::T
    c::T

    scale::Bool   # Should a scaled termination criteria be used
    n::Int        # Number of variables in the problem
    ϵₛ::T         # Scaled termination criteria

    function Gradient( ϵ::T; scaled::Bool = true ) where {T <: Number}
        return new{T}( ϵ, scaled, 0.0, 0.0, scaled, 1, ϵ )
    end
end

function init!( cond::Gradient, n, L, μ )
    cond.μ = μ     #
    cond.L = L     #
    cond.c = 0.5 * ( 1/cond.μ - 1/cond.L )  # Coefficient of the computation

    cond.n = n

    if( cond.scale )
        cond.ϵₛ = cond.ϵ * 1/n
    else
        cond.ϵₛ = cond.ϵ
    end
end

function (cond::Gradient)( state::FGMState )
    return abs( cond.c * norm( cond.L* ( state.y - state.xₙ ) )^2 ) < cond.ϵₛ
end

"""
```julia
    Best( ϵ₁::T, ϵ₂:T = ϵ₁; scaled::Bool = true ) where {T <: Number}
```

Terminate the Fast Gradient Method when either the gradient value is below `ϵ₁` or the
conjugate residual value is below `ϵ₂`. `ϵ₁` and `ϵ₂` are divided by the number of variables
in the problem to produce a scaled termination criteria when `scaled` is true.

This criteria is based on the results in §6.4 of [1].

[1] S. Richter, ‘Computational complexity certification of gradient methods for real-time model predictive control’, Doctoral Thesis, ETH Zurich, 2012.

See also: [`Gradient`](@ref) and [`Conjugate`](@ref)
"""
struct Best <: AbstractStopCondition
    grad::Gradient
    conj::Conjugate

    function Best( ϵ₁::T, ϵ₂::T = ϵ₁; scaled::Bool = true ) where {T <: Number}
        return new( Gradient( ϵ₁; scaled = scaled ), Conjugate( ϵ₂; scaled = scaled ) )
    end
end

function init!( cond::Best, n, L, μ )
    init!( cond.grad, n, L, μ )
    init!( cond.conj, n, L, μ )
end

function (cond::Best)( state::FGMState )
    return cond.grad( state ) || cond.conj( state )
end


text:
## --- String parsing functions

    @test parsedlm("1,2,3\n4,5,6\n7,8,9\n", ',', Float64) == reshape(1:9,3,3)'
    @test parsedlm("1,2,3,4\n5,6,7,8\n9,10,11,12\n13,14,15,16", ',', Int64) == reshape(1:16,4,4)'

    A = delim_string_function(x -> delim_string_parse(x, ',', Float32),
        "1,2,3,4\n5,6,7,8\n9,10,11,12\n13,14,15,16", '\n', Array{Float32,1})
    @test isa(A, Array{Array{Float32,1},1})
    @test all([A[i][j] == (i-1)*4 + j for i=1:4, j=1:4])

    A = delim_string_function(x -> delim_string_parse(x, ',', Int64, merge=true, undefval=0),
        "\n1,2,3,,4\n5,6,,7,8\n9,10,,,,11,12\n\n\n13,14,15,16", '\n', Array{Int64,1}, merge=true)
    @test all([A[i][j] == (i-1)*4 + j for i=1:4, j=1:4])

## --- Elementify/unelementify functions

    elements = ["U" "Lv" "Te" "O" "Er" "W" "Re" "j" "asdf" "Zr" "Al" "S" "K" "V" "N" "Ga" "I"]
    data = vcat(elements, hcat(rand(1000, length(elements)-1), string.(rand("abcdefghijklmnopqrstuvwxyz0123456789",1000))))
    datatuple = elementify(data,importas=:Tuple)::NamedTuple
    datadict = elementify(data,importas=:Dict)::Dict

    @test isa(datatuple, NamedTuple)
    @test unelementify(datatuple) == data
    @test isa(datadict, Dict)
    @test unelementify(datadict) == data

    # Test adding or averaging option for numeric elements
    addtest = ["a" "b" "a";1 2 3]
    avg = elementify(addtest)
    add = elementify(addtest, sumduplicates=true)
    @test avg["elements"] == avg["elements"]
    @test avg["a"] == 2
    @test add["a"] == 4

## --- Import / export functions

    @test exportdataset(datatuple, "tupledataset.csv", ',') == nothing
    @test importdataset("tupledataset.csv", ',', importas=:Tuple) == datatuple

    @test exportdataset(datadict, datadict["elements"], "dictdataset.csv", ',') == nothing
    @test importdataset("dictdataset.csv", ',', importas=:Dict, mindefinedcolumns=2) == datadict

## --  Normalization functions

    dataarray = rand(1000, length(elements))
    data = vcat(elements, dataarray)
    datatuple = elementify(data,importas=:Tuple)::NamedTuple
    datadict = elementify(data,importas=:Dict)::Dict

    renormalize!(dataarray, total=100)
    @test nansum(dataarray) ≈ 100
    renormalize!(dataarray, dim=1, total=100)
    @test all(nansum(dataarray, dims=1) .≈ 100)
    renormalize!(dataarray, dim=2, total=100)
    @test all(nansum(dataarray, dims=2) .≈ 100)

    # Renormalization functions on NamedTuple-based dataset
    datatuple = elementify(unelementify(datatuple, findnumeric=true), importas=:Tuple)
    renormalize!(datatuple, total=100)
    @test all(sum(unelementify(datatuple, floatout=true),dims=2) .≈ 100)

    # Renormalization functions on Dict-based dataset
    datadict = elementify(unelementify(datadict, findnumeric=true), importas=:Dict)
    renormalize!(datadict, datadict["elements"], total=100.)
    @test all(sum(unelementify(datadict, floatout=true),dims=2) .≈ 100)

    # Internal standardization functions
    @test isnan(StatGeochem.floatify("asdf"))
    @test StatGeochem.floatify("12345", Float64) === 12345.0
    @test StatGeochem.floatify("12345", Float32) === 12345f0
    @test StatGeochem.floatify(12345, Float64) === 12345.0
    @test StatGeochem.floatify(12345, Float32) === 12345f0

    @test isa(StatGeochem._columnformat(["asdf","qwer","zxcv"], false), Array{String,1})
    @test isa(StatGeochem._columnformat([1f0, 2f0, 3f0], false), Array{Float32,1})
    @test isa(StatGeochem._columnformat([1., 2., 3.], false), Array{Float64,1})
    @test isa(StatGeochem._columnformat([0x01,0x02,0x03], false), Array{UInt8,1})
    @test isa(StatGeochem._columnformat([1,2,3], false), Array{Int64,1})
    @test all(StatGeochem._columnformat([0x01,2,"3"], false) .=== [0x01,2,"3"])
    @test StatGeochem._columnformat([0x01,2,"3"], true) == [1,2,3]

## --- Concatenating and merging datasets

    d2 = concatenatedatasets(datadict, datadict)
    @test isa(d2, Dict)

    d2array = unelementify(d2, floatout=true)
    @test isa(d2array, Array{Float64,2})
    @test size(d2array) == (2000, length(datadict["elements"]))

## ---


text:
#=

Thin-veener over the VISA shared library.
See VPP-4.3.2 document for details.

=#

############################ Types #############################################

#Vi datatypes
#Cribbed from VPP-4.3.2 section 3.1 table and/or visa.h
#It's most likely we don't actually need all of these but they're easy to
#generate with some metaprogramming

for typePair = [("UInt32", UInt32),
				("Int32", Int32),
				("UInt64", UInt64),
				("Int64", Int64),
				("UInt16", UInt16),
				("Int16", Int16),
				("UInt8", UInt8),
				("Int8", Int8),
				("Addr", Void),
				("Char", Int8),
				("Byte", UInt8),
				("Boolean", UInt16),
				("Real32", Float32),
				("Real64", Float64),
				("Status", Int32),
				("Version", UInt32),
				("Object", UInt32),
				("Session", UInt32)
				]

	viTypeName = symbol("Vi"*typePair[1])
	viConsructorName = symbol("vi"*typePair[1])
	viPTypeName = symbol("ViP"*typePair[1])
	viATypeName = symbol("ViA"*typePair[1])
	@eval begin
		typealias $viTypeName $typePair[2]
		$viConsructorName(x) = convert($viTypeName, x)
		typealias $viPTypeName Ptr{$viTypeName}
		typealias $viATypeName Array{$viTypeName, 1}
	end
end

for typePair = [("Buf", "PByte"),
				("String", "PChar"),
				("Rsrc", "String")
				]
	viTypeName = symbol("Vi"*typePair[1])
	viPTypeName = symbol("ViP"*typePair[1])
	viATypeName = symbol("ViA"*typePair[1])

	mappedViType = symbol("Vi"*typePair[2])

	@eval begin
		typealias $viTypeName $mappedViType
		typealias $viPTypeName $mappedViType
		typealias $viATypeName Array{$viTypeName, 1}
	end
end

typealias ViEvent ViObject
typealias ViPEvent Ptr{ViEvent}
typealias ViFindList ViObject
typealias ViPFindList Ptr{ViFindList}
typealias ViString ViPChar
typealias ViRsrc ViString
typealias ViBuf ViPByte;
typealias ViAccessMode ViUInt32
typealias ViAttr ViUInt32
typealias ViEventType ViUInt32
typealias ViEventFilter ViUInt32

########################## Constants ###########################################

# Completion and Error Codes ----------------------------------------------*/
include("codes.jl")

#Atributes and other definitions
include("constants.jl")

######################### Functions ############################################

#Helper macro to make VISA call and check the status for an error
macro check_status(viCall)
	return quote
		status = $viCall
		if status < VI_SUCCESS
			errMsg = codes[status]
			error("VISA C call failed with status $(errMsg[1]): $(errMsg[2])")
		end
		status
	end
end

#- Resource Manager Functions and Operations -------------------------------#
function viOpenDefaultRM()
	rm = ViSession[0]
	@check_status ccall((:viOpenDefaultRM, libvisa), ViStatus, (ViPSession,), pointer(rm))
	rm[1]
end

function viFindRsrc(sesn::ViSession, expr::AbstractString)
	returnCount = ViUInt32[0]
	findList = ViFindList[0]
	desc = Array(ViChar, VI_FIND_BUFLEN)
	@check_status ccall((:viFindRsrc, libvisa), ViStatus,
						(ViSession, ViString, ViPFindList, ViPUInt32, ViPChar),
						sesn, expr, findList, returnCount, desc)

	#Create the array of instrument strings and push them on
	instrStrs = ASCIIString[bytestring(convert(Ptr{UInt8}, pointer(desc)))]
	while (returnCount[1] > 1)
		@check_status ccall((:viFindNext, libvisa), ViStatus,
						(ViFindList, ViPChar), findList[1], desc)
		returnCount[1] -= 1
		push!(instrStrs, bytestring(convert(Ptr{UInt8}, pointer(desc))))
	end

	instrStrs
end

# ViStatus _VI_FUNC  viParseRsrc     (ViSession rmSesn, ViRsrc rsrcName,
#                                     ViPUInt16 intfType, ViPUInt16 intfNum);
# ViStatus _VI_FUNC  viParseRsrcEx   (ViSession rmSesn, ViRsrc rsrcName, ViPUInt16 intfType,
#                                     ViPUInt16 intfNum, ViChar _VI_FAR rsrcClass[],
#                                     ViChar _VI_FAR expandedUnaliasedName[],
#                                     ViChar _VI_FAR aliasIfExists[]);

function viOpen(sesn::ViSession, name::ASCIIString; mode::ViAccessMode=VI_NO_LOCK, timeout::ViUInt32=VI_TMO_IMMEDIATE)
	#Pointer for the instrument handle
	instrHandle = ViSession[0]
	@check_status ccall((:viOpen, libvisa), ViStatus,
						(ViSession, ViRsrc, ViAccessMode, ViUInt32, ViPSession),
						sesn, name, mode, timeout, instrHandle)
	instrHandle[1]
end

function viClose(viObj::ViObject)
	@check_status ccall((:viClose, libvisa), ViStatus, (ViObject,), viObj)
end

# #- Resource Template Operations --------------------------------------------*/

function viSetAttribute(viObj::ViObject, attrName::ViAttr, attrValue::ViAttrState)
	@check_status ccall((:viSetAttribute, libvisa), ViStatus,
						(ViObject, ViAttr, ViAttrState),
						viObj, attrName, attrValue)
end

function viGetAttribute(viObj::ViObject, attrName::ViAttr)
	value = ViAttrState[0]
	@check_status ccall((:viGetAttribute, libvisa), ViStatus,
						(ViObject, ViAttr, Ptr{Void}),
						viObj, attrName, value)
	value[]
end

# ViStatus _VI_FUNC  viStatusDesc    (ViObject vi, ViStatus status, ViChar _VI_FAR desc[]);
# ViStatus _VI_FUNC  viTerminate     (ViObject vi, ViUInt16 degree, ViJobId jobId);

# ViStatus _VI_FUNC  viLock          (ViSession vi, ViAccessMode lockType, ViUInt32 timeout,
#                                     ViKeyId requestedKey, ViChar _VI_FAR accessKey[]);
# ViStatus _VI_FUNC  viUnlock        (ViSession vi);

function viEnableEvent(instrHandle::ViSession, eventType::Integer,
					   mechanism::Integer)
	@check_status ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16, ViEventFilter),
						 instrHandle, eventType, mechanism, 0)
end

function viDisableEvent(instrHandle::ViSession, eventType::Integer,
					   mechanism::Integer)
	@check_status ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16),
						 instrHandle, eventType, mechanism)
end

function viDiscardEvents(instrHandle::ViSession, eventType::ViEventType,
					   mechanism::UInt16)
	@check_status ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16),
						 instrHandle, eventType, mechanism)
end

function viWaitOnEvent(instrHandle::ViSession, eventType::ViEventType, timeout::UInt32 = VI_TMO_INFINITE)
	outType = Array(ViEventType)
	outEvent = Array(ViEvent)
	@check_status ccall((:viWaitOnEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt32, Ptr{ViEventType}, Ptr{ViEvent}),
						 instrHandle, eventType, timeout, outType, outEvent)
	(outType[], outEvent[])
end

# ViStatus _VI_FUNC  viWaitOnEvent   (ViSession vi, ViEventType inEventType, ViUInt32 timeout,
#                                     ViPEventType outEventType, ViPEvent outContext);

# ViStatus _VI_FUNC  viDisableEvent  (ViSession vi, ViEventType eventType, ViUInt16 mechanism);
# ViStatus _VI_FUNC  viDiscardEvents (ViSession vi, ViEventType eventType, ViUInt16 mechanism);
# ViStatus _VI_FUNC  viWaitOnEvent   (ViSession vi, ViEventType inEventType, ViUInt32 timeout,
#                                     ViPEventType outEventType, ViPEvent outContext);
# ViStatus _VI_FUNC  viInstallHandler(ViSession vi, ViEventType eventType, ViHndlr handler,
#                                     ViAddr userHandle);
# ViStatus _VI_FUNC  viUninstallHandler(ViSession vi, ViEventType eventType, ViHndlr handler,
#                                       ViAddr userHandle);

#- Basic I/O Operations ----------------------------------------------------#

function viWrite(instrHandle::ViSession, data::Union{ASCIIString, Vector{UInt8}})
	bytesWritten = ViUInt32[0]
	@check_status ccall((:viWrite, libvisa), ViStatus,
						(ViSession, ViBuf, ViUInt32, ViPUInt32),
						instrHandle, data, length(data), bytesWritten)
	bytesWritten[1]
end

function viRead!(instrHandle::ViSession, buffer::Array{UInt8})
	bytesRead = ViUInt32[0]
	status = @check_status ccall((:viRead, libvisa), ViStatus,
						(ViSession, ViBuf, ViUInt32, ViPUInt32),
						instrHandle, buffer, sizeof(buffer), bytesRead)
	return (status != VI_SUCCESS_MAX_CNT, bytesRead[])
end

function viRead(instrHandle::ViSession; bufSize::UInt32=0x00000400)
	buf = Array(UInt8, bufSize)
	(done, bytesRead) = viRead!(instrHandle, buf)
	buf[1:bytesRead]
end

function readavailable(instrHandle::ViSession)
	ret = IOBuffer()
	buf = Array(UInt8, 0x400)
	while true
		(done, bytesRead) = viRead!(instrHandle, buf)
		write(ret,buf[1:bytesRead])
		if done
			break
		end
	end
	takebuf_array(ret)
end

# ViStatus _VI_FUNC  viReadAsync     (ViSession vi, ViPBuf buf, ViUInt32 cnt, ViPJobId  jobId);
# ViStatus _VI_FUNC  viReadToFile    (ViSession vi, ViConstString filename, ViUInt32 cnt,
#                                     ViPUInt32 retCnt);
# ViStatus _VI_FUNC  viWriteAsync    (ViSession vi, ViBuf  buf, ViUInt32 cnt, ViPJobId  jobId);
# ViStatus _VI_FUNC  viWriteFromFile (ViSession vi, ViConstString filename, ViUInt32 cnt,
#                                     ViPUInt32 retCnt);
# ViStatus _VI_FUNC  viAssertTrigger (ViSession vi, ViUInt16 protocol);
# ViStatus _VI_FUNC  viReadSTB       (ViSession vi, ViPUInt16 status);
# ViStatus _VI_FUNC  viClear         (ViSession vi);


text:
# This section defines utility functions/structs that are used throughout Baytes, but do not need to be extended
############################################################################################

#########################################
"""
$(TYPEDEF)
Abstract super type for Tuning Algorithms.
"""
abstract type AbstractTune end

#########################################
"""
$(TYPEDEF)
Lower level configuration for individual Algorithms.
"""
abstract type AbstractConfiguration end

#########################################
function init end

"""
    $(FUNCTIONNAME)
Inplace version of [`init`](@ref).
"""
function init! end

"""
    $(FUNCTIONNAME)
Update struct and return new type.
"""
function update end

"""
    $(FUNCTIONNAME)
Update fields of existing struct.
"""
function update! end

############################################################################################
# Export
export
    AbstractTune,
    AbstractDiagnostics,
    AbstractConfiguration,

    init,
    init!,
    update,
    update!


text:
struct CMMatrix{T, C} <: AbstractMatrix{T} # M_r
    cc::Vector{C} # vector of conjugacy classes to fix the order
    r::Int # the index of conjugacy class
    m::Matrix{T} # cache of class coefficients

    function CMMatrix(cc::A, r::Int, T::Type=Int) where {C, A<:AbstractVector{C}}
        M = -ones(T, length(cc), length(cc))
        new{T, C}(cc, r, M)
    end
end

Base.size(M::CMMatrix) = size(M.m)
Base.IndexStyle(::Type{<:CMMatrix}) = IndexCartesian()

function Base.getindex(M::CMMatrix, s::Integer, t::Integer)
    if isone(-M.m[s,t])
        out = one(first(first(M.cc)))

        M.m[s,:] .= 0
        r = M.r

        for g in M.cc[r]
            for h in M.cc[s]
                out = PermutationGroups.mul!(out, g, h)
                for t in 1:size(M, 2)
                    if out == first(M.cc[t])
                        M.m[s, t] += 1
                        break
                    end
                end
            end
        end

    end
    return M.m[s,t]
end


text:
"""
Update prior in `Param` either using modm or modi

* `priori` : prior model on inversion grid
"""
function update_prior!(pa; priori=nothing, priorm=nothing)

	if(priori===nothing) # use current modi as prior!
		if(priorm===nothing)
			copyto!(pa.priori, pa.modi)
		else
			Models.interp_spray!(priorm, pa.priori, :interp)
			copyto!(pa.priori, pa.modi)
		end
	else
		!(isapprox(priori,pa.modi)) && error("priori model has to be on igrid")
		copyto!(pa.priori, priori)
	end
	# update pa.mx.prior
	Models.Seismic_get!(pa.mx.prior,pa.priori,pa.parameterization)
end


text:
@testset "approximate conditioning" begin

end


text:
using Chain
using DataFrames
using Test

using Stonks: SchemaValidationError, to_dataframe
using Stonks.Stores: FileStore, load, save
using Stonks.Stores:
  generate_partition_pairs,
  filter_files,
  list_nested_files,
  list_partition_nesting,
  validate_schema
using Stonks.Models: AssetPrice, AssetInfo

include("test_utils.jl")

@testset "Store functionality" begin
  symbols = ["AAPL", "MSFT"]
  dest = joinpath(@__DIR__, "data/test_stonks")
  prices = to_dataframe(test_price_data())
  info = to_dataframe(test_info_data())
  ds = FileStore{AssetPrice}(; path=dest, ids=["symbol"])
  dsp = FileStore{AssetPrice}(; path=dest, ids=["symbol"], partitions=["symbol"])

  @testset "Schema validation" begin
    @test validate_schema(AssetPrice, prices)
    @test isa(validate_schema(AssetInfo, prices), SchemaValidationError)
  end

  @testset "Store writes - valid schema with no partitions" begin
    save(ds, prices)
    @test readdir(dest) == ["data.csv"]
    @test filesize("$dest/data.csv") > 0
    rm(dest; recursive=true, force=true)
  end

  @testset "Store writes - valid schema with partitions" begin
    save(dsp, prices)
    @test readdir(dest) == map(s -> "symbol=$s", symbols)
    @test minimum(map(s -> filesize("$dest/symbol=$s/data.csv"), symbols)) > 0
    rm(dest; recursive=true, force=true)
  end

  @testset "Store writes - invalid schema" begin
    @test isa(save(dsp, info), SchemaValidationError)
  end

  @testset "Store kwargs constructors" begin
    # ids which are not members of AssetPrice => throw error
    try_ds = try
      FileStore{AssetPrice}(; path=dest, ids=["foo", "bar"])
    catch err
      err
    end
    @test typeof(try_ds) === DomainError
    # partitions which are not members of AssetPrice => throw error
    try_ds = try
      FileStore{AssetPrice}(; path=dest, ids=["symbol"], partitions=["foo", "bar"])
    catch err
      err
    end
    @test typeof(try_ds) === DomainError
    # more than 2 ids should throw a DomainError
    try_ds = try
      FileStore{AssetPrice}(; path=dest, ids=["symbol", "date", "close"])
    catch err
      err
    end
    @test typeof(try_ds) === DomainError
  end

  @testset "Default arguments" begin
    ds = FileStore{AssetPrice}(; path=dest, ids=[:symbol])
    @test ds.time_column == "date"
    @test ds.format == "csv"
    @test ds.partitions == []
    ds = FileStore{AssetInfo}(; path=dest, ids=[:symbol])
    @test ismissing(ds.time_column)
  end

  @testset "Partition pairs" begin
    partitions = generate_partition_pairs(prices, dsp.partitions)
    @test length(partitions) == 2
    @test sort(map(x -> x.key, partitions)) == map(s -> "symbol=$s", symbols)
  end

  @testset "Store read" begin
    save(dsp, prices)
    data = load(dsp)
    @test readdir(dest) == map(s -> "symbol=$s", symbols)
    @test nrow(data) == nrow(prices)
    get_type_param(::Vector{S}) where {S} = S
    df_types = [(name=col, type=get_type_param(data[:, col])) for col in names(data)]
    @test map(x -> x.name, df_types) == [String(x) for x in fieldnames(AssetPrice)]
    @test validate_schema(AssetPrice, data) == true
    rm(dsp.path; recursive=true, force=true)
  end

  @testset "Initialize a datastore" begin
    isdir(dest) && rm(dest; recursive=true, force=true)
    ds = FileStore{AssetPrice}(; path=dest, ids=["symbol"])
    df = load(ds)
    @test nrow(df) == 0
    @test validate_schema(AssetPrice, df)
    rm(dest; recursive=true, force=true)
  end

  @testset "Initialize a datastore with 1 partition" begin
    ds = FileStore{AssetPrice}(; path=dest, ids=["symbol"], partitions=["symbol"])
    df = load(ds)
    @test nrow(df) == 0
    @test validate_schema(AssetPrice, df)
    @test readdir("$dest/symbol=__init__") == ["data.csv"]
    rm(dest; recursive=true, force=true)
  end

  @testset "Initialize a datastore with multiple partitions" begin
    ds = FileStore{AssetPrice}(;
      path=dest, ids=["symbol"], partitions=["symbol", "date"]
    )
    df = load(ds)
    @test nrow(df) == 0
    @test validate_schema(AssetPrice, df)
    @test list_partition_nesting(ds.path) == ds.partitions
    @test readdir("$dest/symbol=__init__/date=__init__") == ["data.csv"]
    rm(dest; recursive=true, force=true)
  end

  @testset "Filter files with a datastore with one partition" begin
    ds = FileStore{AssetPrice}(; path=dest, ids=["symbol"], partitions=["symbol"])
    ref_date, symb = (last_workday(), ["AAPL", "MSFT", "IBM", "GOOG"])
    prices = to_dataframe(fake_price_data(7, ref_date, symb))
    save(ds, prices)
    all_files = list_nested_files(ds.path)
    p_files = filter_files(ds, Dict("symbol" => symbols[1:2]))
    @test length(p_files) < length(all_files)
    @test all(map(x -> !contains(x, "symbol=GOOG"), p_files))

    err = filter_files(ds, Dict("foo" => ["bar", "baz"]))
    @test typeof(err) <: Exception
    rm(dest; recursive=true, force=true)
  end

  @testset "Filter files with a datastore with multiple partitions" begin
    ds = FileStore{AssetPrice}(;
      path=dest, ids=["symbol"], partitions=["symbol", "date"]
    )
    ref_date, symb = (last_workday(), ["AAPL", "MSFT", "IBM", "GOOG"])
    prices = to_dataframe(fake_price_data(7, ref_date, symb))
    save(ds, prices)
    all_files = list_nested_files(ds.path)
    p_files = @chain begin
      map(i -> Dates.format(ref_date - Day(i), "Y-mm-dd"), 1:3)
      filter_files(ds, Dict("symbol" => symbols[1:2], "date" => _))
    end
    @test length(p_files) < length(all_files)
    @test all(map(x -> !contains(x, "symbol=GOOG"), p_files))

    symbols_new = vcat(symb[1:2], ["FOO", "BAR"])
    files_new = filter_files(ds, Dict("symbol" => symbols_new))
    @test length(files_new) < length(all_files) && length(files_new) > length(p_files)
    rm(dest; recursive=true, force=true)
  end

  @testset "Load a datastore using partition pruning" begin
    ds = FileStore{AssetPrice}(; path=dest, ids=["symbol"], partitions=["symbol"])
    ref_date, symb = (last_workday(), ["AAPL", "MSFT", "IBM", "GOOG"])
    prices = to_dataframe(fake_price_data(7, ref_date, symb))
    save(ds, prices)
    df = load(ds, Dict("symbol" => ["AAPL", "MSFT"]))
    @test sort(unique(df[:, :symbol])) == ["AAPL", "MSFT"]
    data = load(ds, Dict("foo" => ["BAR", "BAZ"]))
    @test typeof(data) <: Exception
    rm(dest; recursive=true, force=true)
  end

  isdir(dest) && rm(dest; recursive=true, force=true)
end


text:
#=
day1:
- Julia version: 1.5.2
- Author: Paul.Mealor
- Date: 2020-12-01
=#

using AoC

step(pos, line::AbstractString) = step(pos, parseLine(line, String, Int)...)
step(pos, direction, distance) = step(pos, Val(Symbol(direction)), distance)
step(pos, ::Val{:forward},distance) = pos + [1,0,0]*distance + [0,1,0]*pos[3]*distance
step(pos, ::Val{:up}, distance) = pos + [0,0,1]*distance
step(pos, ::Val{:down},distance) = pos + [0,0,-1] * distance

function exec(lines)
    pos=[0,0,0]
    for line in lines
        pos = step(pos, line)
    end
    return pos[1:2]
end

show(exampleLines(2,1   ) |> ll -> @time exec(ll))
show(lines(2) |> ll -> @time exec(ll))
show(exampleLines(2,1) |> ll -> @time *(exec(ll)...))
show(lines(2) |> ll -> @time *(exec(ll)...))


text:
module AxisSets

using AutoHashEquals
using AxisKeys
using FeatureTransforms
using Impute
using NamedDims
using OrderedCollections
using ReadOnlyArrays

using Impute: DeclareMissings, Filter, Imputor, Validator

export KeyedDataset

# There's a few places calling `only` is convenient, even for older Julia releases
if VERSION < v"1.4"
    function _only(x)
        if isempty(x)
            throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
        elseif length(x) > 1
            throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
        else
            first(x)
        end
    end
else
    _only(x) = only(x)
end

include("flatten.jl")
include("patterns.jl")
include("dataset.jl")
include("indexing.jl")
include("functions.jl")
include("impute.jl")
include("featuretransforms.jl")
include("utils.jl")

end


text:
export Lehmer, random_sample

# MINSTD
const M_BIT_LEN = 31
const M_PLUS_1 = UInt64(2) ^ M_BIT_LEN
const LCG_M = M_PLUS_1 - UInt64(1)  # Mersenne prime, so smaller numbers are always coprime to it
const LCG_A = UInt64(7) ^ 5

"Lehmer RNG, a special case of Linear congruential generator where c = 0."
mutable struct Lehmer
    x::UInt32

    Lehmer(seed::UInt32) = new(seed)
end

function Base.iterate(l::Lehmer, state=l)
    # Fast computes `(a * x) % m` when `m = (2 ^ e - 1)`:
    #
    #     r ⧋ (a * x) % m  s.t. ∃ k: (a * x) = k * m + r
    #
    #     ∃ k: a * x = k * (2 ^ e - 1) + r  ⟹
    #     ∃ k: a * x = k * (2 ^ e) + (r - k)
    #
    # Sometimes, `r - k` is negative, so we have
    #
    #     ∃ k: a * x = (k - 1) * (2 ^ e) + (2 ^ e + r - (k - 1) - 1)
    #
    # Since `r - k << (2 ^ e)`,
    #
    #     k = (a * x) ÷ (2 ^ e)
    #     r = (a * x) % (2 ^ e) + k

    ax = LCG_A * l.x
    k = ax >> 31
    r = (ax & LCG_M) + k
    rr = r + UInt64(1)
    if !iszero(rr & M_PLUS_1)
        r = rr & LCG_M
    end
    l.x = r

    (r, l)
end

UINT32_MAX = typemax(UInt32)

"Randomly samples `k` indices without replacement from 1 to n by Fisher–Yates shuffle."
function random_sample(lehmer::Lehmer, n::UInt32, k::UInt32)::Vector{UInt32}
    samples = collect(UInt32(1):n)

    # Random values from RNG which are greater than `skip_threshold` will be skipped. For dealing
    # with modulo bias.
    for i in 1:k
        n_choices = n - i + 1
        remainder = UINT32_MAX % n_choices
        skip_threshold = UINT32_MAX - (remainder + 1) % n_choices
        j = i + first(Iterators.filter(x -> x <= skip_threshold, lehmer)) % n_choices
        samples[i], samples[j] = samples[j], samples[i]
    end

    samples[1:k]
end

# For testing
function next_lehmer(x::UInt32)::UInt32
    UInt32((LCG_A * x) % LCG_M)
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# Julia-specific extensions to the LSP

@dict_readable struct VersionedTextDocumentPositionParams
    textDocument::TextDocumentIdentifier
    version::Int
    position::Position
end

include("messagedefs.jl")


text:
function paramwalk!(err::LegHemodynamics.Errors,pbar::Vector{Float64},θ::Vector{Float64})
    for i = 1:length(θ)
        θ[i] = rand(Distributions.TruncatedNormal(err.a*θ[i]+(1-err.a)*pbar[i],err.h*
            err.pdev[i],err.lb[i],err.ub[i]));
    end
    return θ
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using Test
using DataWrangler

@testset "normalize" begin

    # normalize! Vector
    x = [1.,2,3,4,5]
    normalize!(x)
    @test x ≈ [-1.2649110640673518, -0.6324555320336759, 0.0, 0.6324555320336759, 1.2649110640673518]

    x = [1.,2,3,4,5]
    normalize!(x; method = "z-score")
    @test x ≈ [-1.2649110640673518, -0.6324555320336759, 0.0, 0.6324555320336759, 1.2649110640673518]

    x = [1.,2,3,4,5]
    normalize!(x; method = "min-max")
    @test x == [0.0, 0.25, 0.5, 0.75, 1.0]

    x = [1.,2,3,4,5]
    normalize!(x; method = "softmax")
    @test x ≈ [0.011656230956039607, 0.03168492079612427, 0.0861285444362687, 0.23412165725273662, 0.6364086465588308]

    x = [1.,2,3,4,5]
    normalize!(x; method = "sigmoid")
    @test x ≈ [0.7310585786300049, 0.8807970779778823, 0.9525741268224334, 0.9820137900379085, 0.9933071490757153]

    x = [1. 2 3; 1 2 3; 1 2 3]
    normalize!(x;  method = "min-max")
    @test x ≈ [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0]

    # normalize Vector
    x = [1,2,3,4,5]
    @test normalize(x) ≈ [-1.2649110640673518, -0.6324555320336759, 0.0, 0.6324555320336759, 1.2649110640673518]

    @test normalize(x; method = "z-score") ≈ [-1.2649110640673518, -0.6324555320336759, 0.0, 0.6324555320336759, 1.2649110640673518]

    @test normalize(x; method = "min-max") == [0.0, 0.25, 0.5, 0.75, 1.0]

    @test normalize(x; method = "softmax") ≈ [0.011656230956039607, 0.03168492079612427, 0.0861285444362687, 0.23412165725273662, 0.6364086465588308]

    @test normalize(x; method = "sigmoid") ≈ [0.7310585786300049, 0.8807970779778823, 0.9525741268224334, 0.9820137900379085, 0.9933071490757153]

    x = [1 2 3; 1 2 3; 1 2 3]
    @test normalize(x; method = "min-max") ≈ [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0]

end


text:
using Swalbe, Test, Statistics, Parameters, Random, LazySets, BSON, FileIO, DataFrames

@testset "Swalbe.jl" begin
    println("Testing allocations")
    include("initialize.jl")
    println("Testing initial states")
    include("initialvalues.jl")
    println("Testing finite differences")
    include("differences.jl")
    println("Testing film pressure")
    include("pressure.jl")
    println("Testing collision operator")
    include("collide.jl")
    println("Testing equilibrium distribution")
    include("equilibrium.jl")
    println("Testing forces")
    include("forcing.jl")
    println("Testing moment calculation")
    include("moments.jl")
    println("Testing measurements")
    include("measures.jl")
    println("Testing with relevant simulations")
    include("simulate.jl")
end


text:
using LombScargle, Distances

export IrregularLombScargle

"""
    IrregularLombScargle(t; tol = 1, n_total = 100000, n_acc = 50000, q = 1)

[`IrregularLombScargle`](@ref) surrogates for unevenly sampled time series
with supporting time steps `t`, generated using the simulated annealing algorithm
described in [^SchreiberSchmitz1999].

[`IrregularLombScargle`](@ref) surrogates (given enough iterations and a low enough
tolerance) preserve the periodogram and the amplitude
distribution of the original signal. For time series with equidistant time steps,
surrogates generated by this method result in surrogates similar to those produced
by the [`IAAFT`](@ref) method.

This algorithm starts with a random permutation of the original data. Then it iteratively
approaches the power spectrum of the original data by swapping two randomly selected values
in the surrogate data if the Minkowski distance of order `q` between the power spectrum of
the surrogate data and the original data is less than before. The iteration procedure ends
when the relative deviation between the periodograms is less than `tol` or when `n_total`
number of tries or `n_acc` number of actual swaps is reached.

[^SchmitzSchreiber1999]: A.Schmitz T.Schreiber (1999). "Testing for nonlinearity in unevenly sampled time series" [Phys. Rev E](https://journaIrregularLombScargle.aps.org/pre/pdf/10.1103/PhysRevE.59.4044)
"""
struct IrregularLombScargle{T<:AbstractVector, S<:Real} <: Surrogate
    t::T
    tol::S
    n_total::Int
    n_acc::Int
    q::Int
end
IrregularLombScargle(t; tol = 1.0, n_total = 100000, n_acc = 50000, q = 1) = IrregularLombScargle(t, tol, n_total, n_acc, q)

function surrogenerator(x, method::IrregularLombScargle, rng = Random.default_rng())
    # Plan Lombscargle periodogram. Use default flags.
    lsplan = LombScargle.plan(method.t, x, fit_mean = false)

    # Compute initial periodogram.
    x_ls = lombscargle(lsplan)

    # We have to copy the power vector here, because we are reusing `lsplan` later on
    xpower = copy(x_ls.power)

    # Use Minkowski distance of order q
    dist = Distances.Minkowski(method.q)

    # When re-computing the Lomb-Scargle periodogram, we will use the
    # `_periodogram!` method, which re-uses the lsplan with a shuffled
    # time vector. This is the same as shuffling the signal, so the
    # surrogate starts out as a shuffled version of `t`.
    s = shuffle(rng, method.t)

    # Initialize a candidate surrogate.
    candidate = zero(s)

    init = (
        lsplan = lsplan,
        xpower = xpower,
        n = length(x),
        dist = dist,
        candidate = candidate,
    )

    return SurrogateGenerator(method, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:IrregularLombScargle})()
    lsplan, xpower, n, dist, candidate = sg.init
    t, tol = sg.method.t, sg.method.tol
    x, s, rng = sg.x, sg.s, sg.rng

    # Power spectrum for the randomly shuffled signal.
    spower = LombScargle._periodogram!(lsplan.P, s, lsplan)

    # Compare power spectra for original (`xpower`) and randomly shuffled signal (`spower`).
    lossold = Distances.evaluate(dist, xpower, spower)

    range = 1:n
    i = j = 0
    while i < sg.method.n_total && j < sg.method.n_acc
        if mod(i, 2000) == 0
            @info "iterations: $i, swaps: $j, loss: $lossold"
        end

        # Initially, the new surrogate is identical to the existing surrogate.
        copy!(candidate, s)

        # Swap two random points and re-compute power spectrum for the candidate.
        swap_elements!(candidate, s, rng, range)

        # This is the bottleneck.
        spower = LombScargle._periodogram!(lsplan.P, candidate, lsplan)

        # If spectra are more similar after the swap, accept the new
        # surrogate. Otherwise, do a new swap.
        lossnew = evaluate(dist, xpower, spower)

        if lossnew < lossold
            lossnew <= tol && break
            s = copy(candidate)
            lossold = lossnew
            j += 1
        end
        i += 1
    end
    @info "Terminated simulated annealing process after $i iterations and $j swaps. Loss: $lossold"

    # Use the permutation of the time vector to permute the signal vector
    # This gives us the inverse permutation from t to perm
    perm = sortperm(sortperm(s))

    # Check if this worked as expected
    @assert t[perm] == s

    # Re-scale back to original time series values.
    return x[perm]
end

function swap_elements!(candidate, s, rng, range)
    k = sample(rng, range)
    l = sample(rng, range)
    candidate[k] = s[l]
    candidate[l] = s[k]
end


text:
using LinearAlgebra
using Statistics

x = [ 2, -1, 2];

# norm of x is written in Julia as norm(x)
norm(x)

sqrt(x'*x)

sqrt(sum(x.^2))

x = randn(10); y = randn(10)

lhs = norm(x+y)

rhs = norm(x) + norm(y)

# rms value of a vector is rms(x) = norm(x)/sqrt(length(x))
rms(x) = norm(x) / sqrt(length(x))

t = 0:0.01:1; # list of times

x = cos.(8*t) - 2*sin.(11*t)

rms(x)

using Plots

plot(t,x)

plot!(t, avg(x)*ones(length(x)))

plot!(t, (avg(x)+rms(x))*ones(length(x)), color = :green)

plot!(t, (avg(x)-rms(x))*ones(length(x)), color = :green)

plot!(legend = false)

# Define Chebyshev bound function
# The Chebyshev inequality states that the number of entries of an
# n-vector x that have absolute value at least a is no more than
# ∥x∥2/a2 = nrms(x)2/a2. If this number is, say, 12.15,
# we can conclude that no more that 12 entries have absolute value
# at least a, since the number of entries is an integer.

cheb_bound(x,a) = floor(norm(x)^2/a);
a = 1.5;
cheb_bound(x,a)
# Number of entries of x with |x_i| >= a
sum(abs.(x) .>= a)

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
function laser_over_sensors(laser::PlaneWave, sensors)
    direction, s, Δ = laser.direction, laser.s, laser.Δ

    E₀ = estimate_E₀(laser)
    atomic_component = select_sensor_axes(sensors, direction)

    Threads.@threads for n in 1:Nsensors #
        @inbounds E[n] = E₀ * cis(k₀ * atomic_component[n])
    end
    return E
end

function laser_over_sensors(laser::Gaussian3D, sensors)
    ω₀, s, Δ = laser.ω₀, laser.s, laser.Δ
    Nsensors = length(sensors)

    E = zeros(ComplexF64, Nsensors)
    E₀ = estimate_E₀(laser)

    Threads.@threads for n in 1:Nsensors #
        oneSensor = get_one_sensor(sensors, n)
        @inbounds E[n] = _core_Gaussian(oneSensor, E₀, ω₀)
    end
    return E
end

### --------------- PLANE WAVES ---------------
laser_over_atoms(laser::PlaneWave, atoms) = apply_plane_wave_field(laser, atoms)

function apply_plane_wave_field(laser, atoms)
    direction, s, Δ = laser.direction, laser.s, laser.Δ

    E₀ = estimate_E₀(laser)
    atomic_component = select_atoms_axes(atoms, direction)
    return E₀ .* exp.(+im * k₀ * atomic_component)
end

estimate_E₀(laser) = √(laser.s * (1 + 4(laser.Δ / Γ)^2) / 2)
### --------------- GAUSSIANS ---------------
laser_over_atoms(laser::Gaussian3D, atoms) = apply_gaussian3D_field(laser, atoms)

"""
    Gaussian Beam along z direction
"""
function apply_gaussian3D_field(laser, atoms)
    ω₀, s, Δ = laser.ω₀, laser.s, laser.Δ

    E = zeros(ComplexF64, atoms.N)
    E₀ = estimate_E₀(laser)

    Threads.@threads for n in 1:(atoms.N) #
        oneAtom = get_one_atom(atoms, n)
        @inbounds E[n] = _core_Gaussian(oneAtom, E₀, ω₀)
    end
    return E
end

function _core_Gaussian(oneAtom, E₀, ω₀)
    x, y, z = oneAtom[1], oneAtom[2], oneAtom[3]

    ## This formula is stable for z==0
    # Ref: Eq 3.11 from "CHAPTER 3. PROPAGATION AND FOCUSING OF OPTICAL FIELDS"
    denominator_factor = 1 .+ 2im .* z / (k₀ * ω₀^2)
    Eᵢ = E₀ .* exp.(+im * k₀ * z)
    Eᵢ = Eᵢ .* exp.(-(x .^ 2 + y .^ 2) ./ (denominator_factor .* ω₀^2))
    Eᵢ = Eᵢ ./ denominator_factor
    return Eᵢ
end

function estimate_waist(atoms::Cube)
    return size(atoms) / 4
end
function estimate_waist(atoms::Sphere)
    return size(atoms) / 4
end


text:
import Base: @deprecate

@deprecate DataFrame(t::Type, nrows::Integer, ncols::Integer) DataFrame([Vector{t}(undef, nrows) for i in 1:ncols])

@deprecate DataFrame(column_eltypes::AbstractVector{<:Type},
                     nrows::Integer) DataFrame(column_eltypes, Symbol.('x' .* string.(1:length(column_eltypes))), nrows)

function DataFrame(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol},
                   categorical::AbstractVector{Bool}, nrows::Integer;
                   makeunique::Bool=false)::DataFrame where T<:Type
    Base.depwarn("`DataFrame` constructor with `categorical` positional argument is deprecated. " *
                 "Instead use `DataFrame(columns, names)` constructor.",
                 :DataFrame)
    updated_types = convert(Vector{Type}, column_eltypes)
    if length(categorical) != length(column_eltypes)
        throw(DimensionMismatch("arguments column_eltypes and categorical must have the same length " *
                                "(got $(length(column_eltypes)) and $(length(categorical)))"))
    end
    for i in eachindex(categorical)
        categorical[i] || continue
        elty = CategoricalValue{nonmissingtype(updated_types[i]),
                                CategoricalArrays.DefaultRefType}
        if updated_types[i] >: Missing
            updated_types[i] = Union{elty, Missing}
        else
            updated_types[i] = elty
        end
    end
    return DataFrame(updated_types, cnames, nrows, makeunique=makeunique)
end

import Base: show
@deprecate show(io::IO, df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol, summary::Bool) show(io, df, allcols=allcols, rowlabel=rowlabel, summary=summary)
@deprecate show(io::IO, df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol) show(io, df, allcols=allcols, rowlabel=rowlabel)
@deprecate show(io::IO, df::AbstractDataFrame, allcols::Bool) show(io, df, allcols=allcols)
@deprecate show(df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol, summary::Bool) show(df, allcols=allcols, rowlabel=rowlabel, summary=summary)
@deprecate show(df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol) show(df, allcols=allcols, rowlabel=rowlabel)
@deprecate show(df::AbstractDataFrame, allcols::Bool) show(df, allcols=allcols)
@deprecate showall(io::IO, df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol, summary::Bool) show(io, df, allrows=true, allcols=allcols, rowlabel=rowlabel, summary=summary)
@deprecate showall(io::IO, df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol) show(io, df, allrows=true, allcols=allcols, rowlabel=rowlabel)
@deprecate showall(io::IO, df::AbstractDataFrame, allcols::Bool = true) show(io, df, allrows=true, allcols=allcols)
@deprecate showall(df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol, summary::Bool) show(df, allrows=true, allcols=allcols, rowlabel=rowlabel, summary=summary)
@deprecate showall(df::AbstractDataFrame, allcols::Bool, rowlabel::Symbol) show(df, allrows=true, allcols=allcols, rowlabel=rowlabel)
@deprecate showall(df::AbstractDataFrame, allcols::Bool = true) show(df, allrows=true, allcols=allcols)
@deprecate showall(io::IO, dfvec::AbstractVector{T}) where {T <: AbstractDataFrame} foreach(df->show(io, df, allrows=true, allcols=true), dfvec)
@deprecate showall(dfvec::AbstractVector{T}) where {T <: AbstractDataFrame} foreach(df->show(df, allrows=true, allcols=true), dfvec)
@deprecate showall(io::IO, df::GroupedDataFrame) show(io, df, allgroups=true)
@deprecate showall(df::GroupedDataFrame) show(df, allgroups=true)

@deprecate colwise(f, d::AbstractDataFrame) [f(col) for col in eachcol(d)]
@deprecate colwise(fns::Union{AbstractVector, Tuple}, d::AbstractDataFrame) [f(col) for f in fns, col in eachcol(d)]
@deprecate colwise(f, gd::GroupedDataFrame) [[f(col) for col in eachcol(d)] for d in gd]
@deprecate colwise(fns::Union{AbstractVector, Tuple}, gd::GroupedDataFrame) [[f(col) for f in fns, col in eachcol(d)] for d in gd]

@deprecate deletecols!(df::DataFrame, inds) select!(df, Not(inds))
@deprecate deletecols(df::DataFrame, inds; copycols::Bool=true) select(df, Not(inds), copycols=copycols)

import Base: getindex
@deprecate getindex(df::DataFrame, col_ind::ColumnIndex) df[!, col_ind]
@deprecate getindex(df::DataFrame, col_inds::Union{AbstractVector, Regex, Not}) df[:, col_inds]
@deprecate getindex(df::DataFrame, ::Colon) df[:, :]
@deprecate getindex(sdf::SubDataFrame, colind::ColumnIndex) sdf[!, colind]
@deprecate getindex(sdf::SubDataFrame, colinds::Union{AbstractVector, Regex, Not}) sdf[!, colinds]
@deprecate getindex(sdf::SubDataFrame, ::Colon) sdf[!, :]

import Base: view
@deprecate view(adf::AbstractDataFrame, colind::ColumnIndex) view(adf, :, colind)
@deprecate view(adf::AbstractDataFrame, colinds) view(adf, :, colinds)

import Base: setindex!
@deprecate setindex!(df::DataFrame, x::Nothing, col_ind::Int) select!(df, Not(col_ind))
@deprecate setindex!(sdf::SubDataFrame, val::Any, colinds::Any) (sdf[:, colinds] = val; sdf)
@deprecate setindex!(df::DataFrame, v::AbstractVector, col_ind::ColumnIndex) (df[!, col_ind] = v; df)

# df[SingleColumnIndex] = Single Item (EXPANDS TO NROW(df) if NCOL(df) > 0)
function Base.setindex!(df::DataFrame, v, col_ind::ColumnIndex)
    if haskey(index(df), col_ind)
        df[:, col_ind] .= v
        Base.depwarn("Implicit broadcasting to an existing column in DataFrame assignment is deprecated." *
                     "Use an explicit broadcast with `df[:, col_ind] .= v`", :setindex!)
    else
        if ncol(df) == 0
            df[!, col_ind] = [v]
            Base.depwarn("Implicit broadcasting to a new column in DataFrame assignment is deprecated." *
                         "Use `df[!, col_ind] = [v]` when `df` has zero columns", :setindex!)
        else
            df[!, col_ind] .= v
            Base.depwarn("Implicit broadcasting to a new column in DataFrame assignment is deprecated." *
                         "Use `df[!, col_ind] .= v`  when `df` has some columns", :setindex!)
        end
    end
    return df
end

# df[MultiColumnIndex] = DataFrame
function Base.setindex!(df::DataFrame, new_df::DataFrame, col_inds::AbstractVector{Bool})
    setindex!(df, new_df, findall(col_inds))
end
@deprecate setindex!(df::DataFrame, new_df::DataFrame,
                     col_inds::AbstractVector{<:ColumnIndex}) foreach(((j, colind),) -> (df[!, colind] = new_df[!, j]),
                                                                      enumerate(col_inds))

# df[MultiColumnIndex] = AbstractVector (REPEATED FOR EACH COLUMN)
@deprecate setindex!(df::DataFrame, v::AbstractVector,
                     col_inds::AbstractVector{Bool}) (foreach(c -> (df[!, c] = copy(v)), findall(col_inds)); df)
@deprecate setindex!(df::DataFrame, v::AbstractVector,
                     col_inds::AbstractVector{<:ColumnIndex}) (foreach(c -> (df[!, c] = copy(v)), col_inds); df)

# df[MultiColumnIndex] = Single Item (REPEATED FOR EACH COLUMN; EXPANDS TO NROW(df) if NCOL(df) > 0)
function Base.setindex!(df::DataFrame,
                        val::Any,
                        col_inds::AbstractVector{Bool})
    setindex!(df, val, findall(col_inds))
end
function Base.setindex!(df::DataFrame, val::Any, col_inds::AbstractVector{<:ColumnIndex})
    for col_ind in col_inds
        df[col_ind] = val
    end
    Base.depwarn("implicit broadcasting in setindex! is deprecated; " *
                 "use `df[:, col_inds] .= Ref(v)` broadcasting assignment to change the columns in place", :setindex!)
    return df
end

# df[:] = AbstractVector or Single Item
function Base.setindex!(df::DataFrame, v, ::Colon)
    df[1:size(df, 2)] = v
    Base.depwarn("`df[:] = v` syntax is deprecated; " *
                 "use `df[:, :] .= Ref(v)` broadcasting assignment to change the columns in place", :setindex!)
    df
end

# df[SingleRowIndex, MultiColumnIndex] = 1-Row DataFrame
@deprecate setindex!(df::DataFrame, new_df::DataFrame, row_ind::Integer,
                     col_inds::AbstractVector{Bool}) (foreach(((i, c),) -> (df[row_ind, c] = new_df[1, i]),
                                                              enumerate(findall(col_inds))); df)
@deprecate setindex!(df::DataFrame, new_df::DataFrame, row_ind::Integer,
                     col_inds::AbstractVector{<:ColumnIndex}) (foreach(((i, c),) -> (df[row_ind, c] = new_df[1, i]),
                                                                       enumerate(col_inds)); df)

# df[MultiRowIndex, SingleColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                        v::Any,
                        row_inds::AbstractVector,
                        col_ind::ColumnIndex)
    insert_multiple_entries!(df, v, row_inds, col_ind)
    Base.depwarn("implicit broadcasting in setindex! is deprecated; " *
                 "use `df[row_inds, col_ind] .= Ref(v)` broadcasting assignment to change the column in place", :setindex!)
    return df
end

# df[MultiRowIndex, MultiColumnIndex] = AbstractVector
function Base.setindex!(df::DataFrame,
                        v::AbstractVector,
                        row_inds::AbstractVector,
                        col_inds::AbstractVector)
    col_inds_norm = index(df)[col_inds]
    for col_ind in col_inds_norm
        insert_multiple_entries!(df, v, row_inds, col_ind)
    end
    Base.depwarn("implicit vector broadcasting in setindex! is deprecated; " *
                 "use `df[row_inds, col_inds] .= v` broadcasting assignment to change the columns in place", :setindex!)
    return df
end

# df[MultiRowIndex, MultiColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                        v::Any,
                        row_inds::AbstractVector,
                        col_inds::AbstractVector)
    col_inds_norm = index(df)[col_inds]
    for col_ind in col_inds_norm
        insert_multiple_entries!(df, v, row_inds, col_ind)
    end
    Base.depwarn("implicit vector broadcasting in setindex! is deprecated; " *
                 "use `df[row_inds, col_inds] .= Ref(v)` broadcasting assignment to change the columns in place", :setindex!)
    return df
end

# df[:, :] = ...
function Base.setindex!(df::DataFrame, v, ::Colon, ::Colon)
    df[1:size(df, 1), 1:size(df, 2)] = v
    Base.depwarn("implicit vector broadcasting in setindex! is deprecated; " *
        "use `df[:, col_inds] .= Ref(v)` broadcasting assignment to change the columns in place", :setindex!)
    df
end

# df[Any, :] = ...
function Base.setindex!(df::DataFrame, v, row_inds, ::Colon)
    df[row_inds, 1:size(df, 2)] = v
    Base.depwarn("implicit vector broadcasting in setindex! is deprecated; " *
                 "use `df[row_inds, col_inds] .= Ref(v)` broadcasting assignment", :setindex!)
    df
end

# df[:, Any] = ...
function Base.setindex!(df::DataFrame, v, ::Colon, col_inds)
    df[col_inds] = v
    Base.depwarn("implicit vector broadcasting in setindex! is deprecated; " *
                 "use `df[:, col_inds] .= Ref(v)` broadcasting assignment to change the columns in place", :setindex!)
    df
end

function Base.dotview(df::AbstractDataFrame, col::ColumnIndex)
    Base.depwarn("in broadcasted assignment use `df[:, col]` instead of `df[col]`", :dotview)
    @view df[:, col]
end

import Base: setproperty!
@deprecate setproperty!(df::DataFrame, col_ind::Symbol, v) (df[!, col_ind] .= v)
@deprecate setproperty!(df::SubDataFrame, col_ind::Symbol, v) (df[:, col_ind] .= v)

@deprecate eltypes(df::AbstractDataFrame) eltype.(eachcol(df))

@deprecate permutecols!(df::DataFrame, p::AbstractVector) select!(df, p)
@deprecate names!(x::Index, nms::Vector{Symbol};
                  makeunique::Bool=false) rename!(x, nms, makeunique=makeunique)
@deprecate names!(df::AbstractDataFrame, vals::Vector{Symbol};
                  makeunique::Bool=false) rename!(df, vals, makeunique=makeunique)

import DataAPI: describe
@deprecate describe(io::IO, df::AbstractDataFrame, stats::Union{Symbol, Pair{Symbol}}...;
                    cols=:) describe(df, stats..., cols=cols)

@deprecate stackdf(args...; kwargs...) stack(args...; kwargs..., view=true)
@deprecate meltdf(args...; kwargs...) melt(args...; kwargs..., view=true)

@deprecate melt(df::AbstractDataFrame, id_vars;
                variable_name::Symbol=:variable, value_name::Symbol=:value,
                view::Bool=false) stack(df, Not(id_vars); variable_name=variable_name,
                                        value_name=value_name, view=view)

@deprecate melt(df::AbstractDataFrame, id_vars, measure_vars;
                variable_name::Symbol=:variable, value_name::Symbol=:value,
                view::Bool=false) stack(df, measure_vars, id_vars; variable_name=variable_name,
                                        value_name=value_name, view=view)
@deprecate melt(df::AbstractDataFrame; variable_name::Symbol=:variable, value_name::Symbol=:value,
                view::Bool=false) stack(df; variable_name=variable_name, value_name=value_name,
                                        view=view)

import Base: lastindex
@deprecate lastindex(df::AbstractDataFrame) ncol(df)

function insertcols!(df::DataFrame, col_ind::Int; makeunique::Bool=false, name_col...)
    Base.depwarn("inserting colums using a keyword argument is deprecated, use " *
                 "`insertcols!(df, col_ind, :$(keys(name_col)[1]) => $(name_col[1]), " *
                 "makeunique=$makeunique)` instead. ", :insertcols!)
    length(name_col) == 1 || throw(ArgumentError("one and only one column must be provided"))
    return insertcols!(df, col_ind, makeunique=makeunique, keys(name_col)[1] => name_col[1])
end

function Base.join(df1::AbstractDataFrame, df2::AbstractDataFrame;
                   on::Union{<:OnType, AbstractVector} = Symbol[],
                   kind::Symbol = :inner, makeunique::Bool=false,
                   indicator::Union{Nothing, Symbol} = nothing,
                   validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false))
    if kind == :inner
        Base.depwarn("inner joining data frames using join is deprecated, use " *
                     "`innerjoin(df1, df2, on=$on, makeunique=$makeunique, validate=$validate)` instead", :join)
        if isnothing(indicator)
            innerjoin(df1, df2, on=on, makeunique=makeunique, validate=validate)
        else
            throw(ArgumentError("`indicator` keyword argument is not supported for `innerjoin`"))
        end
    elseif kind == :left
        Base.depwarn("left joining data frames using join is deprecated, use " *
                     "`leftjoin(df1, df2, on=$on, makeunique=$makeunique, indicator=$indicator, validate=$validate)` instead", :join)
        leftjoin(df1, df2, on=on, makeunique=makeunique, indicator=indicator, validate=validate)
    elseif kind == :right
        Base.depwarn("right joining data frames using join is deprecated, use " *
                     "`rightjoin(df1, df2, on=$on, makeunique=$makeunique, indicator=$indicator, validate=$validate)` instead", :join)
        rightjoin(df1, df2, on=on, makeunique=makeunique, indicator=indicator, validate=validate)
    elseif kind == :outer
        Base.depwarn("outer joining data frames using join is deprecated, use " *
                     "`outerjoin(df1, df2, on=$on, makeunique=$makeunique, indicator=$indicator, validate=$validate)` instead", :join)
        outerjoin(df1, df2, on=on, makeunique=makeunique, indicator=indicator, validate=validate)
    elseif kind == :semi
        Base.depwarn("semi joining data frames using join is deprecated, use " *
                     "`semijoin(df1, df2, on=$on, makeunique=$makeunique, validate=$validate)` instead", :join)
        semijoin(df1, df2, on=on, makeunique=makeunique, validate=validate)
    elseif kind == :anti
        Base.depwarn("anti joining data frames using join is deprecated, use " *
                     "`antijoin(df1, df2, on=$on, makeunique=$makeunique, validate=$validate)` instead", :join)
        antijoin(df1, df2, on=on, makeunique=makeunique, validate=validate)
    elseif kind == :cross
        Base.depwarn("cross joining data frames using join is deprecated, use " *
                     "`crossjoin(df1, df2, makeunique=$makeunique)` instead", :join)
        crossjoin(df1, df2, makeunique=makeunique)
    else
        throw(ArgumentError("Unknown kind of join requested: $kind"))
    end
end

function Base.join(df1::AbstractDataFrame, df2::AbstractDataFrame,
                   dfs::AbstractDataFrame...;
                   on::Union{Symbol, AbstractVector{Symbol}} = Symbol[],
                   kind::Symbol = :inner, makeunique::Bool=false,
                   validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false))
    if kind == :inner
        Base.depwarn("inner joining data frames using join is deprecated, use " *
                     "`innerjoin(df1, df2, dfs..., on=$on, makeunique=$makeunique, validate=$validate)` instead", :join)
        innerjoin(df1, df2, dfs..., on=on, makeunique=makeunique, validate=validate)
    elseif kind == :outer
        Base.depwarn("outer joining data frames using join is deprecated, use " *
                     "`outerjoin(df1, df2, dfs..., on=$on, makeunique=$makeunique, validate=$validate)` instead", :join)
        outerjoin(df1, df2, dfs..., on=on, makeunique=makeunique, validate=validate)
    elseif kind == :cross
        Base.depwarn("cross joining data frames using join is deprecated, use " *
                     "`crossjoin(df1, df2, dfs..., makeunique=$makeunique)` instead", :join)
        crossjoin(df1, df2, dfs..., makeunique=makeunique)
    else
        throw(ArgumentError("Only inner, outer, and cross joins are supported when " *
                            "joining more than two data frames"))
    end
end

@deprecate eachcol(df::AbstractDataFrame, names::Bool) names ? collect(pairs(eachcol(df))) : collect(eachcol(df))

@deprecate groupvars(gd::GroupedDataFrame) groupcols(gd)

export aggregate

function aggregate(d::AbstractDataFrame, f::Any; sort::Bool=false)
    df = combine(d, names(d) .=> [f])
    if sort
        Base.depwarn("`aggregate(d, f, sort=true)` is deprecated. " *
                     "Instead use `sort!(combine(d, names(d) .=> f))`.", :aggregate)
        sort!(df)
    else
        Base.depwarn("`aggregate(d, f)` is deprecated. " *
                     "Instead use `combine(d, names(d) .=> f)`.", :aggregate)
    end
    return df
end

function aggregate(d::AbstractDataFrame, fs::AbstractVector; sort::Bool=false)
    df = hcat([combine(d, names(d) .=> [f]) for f in fs]..., makeunique=true)
    if sort
        Base.depwarn("`aggregate(d, fs, sort=true)` is deprecated. Instead" *
                     " use `sort!(combine(d, [names(d) .=> f for f in fs]...))` " *
                     "if functions in `fs` have unique names.", :aggregate)
        sort!(df)
    else
        Base.depwarn("`aggregate(d, fs)` is deprecated. Instead use " *
                     "`combine(d, [names(d) .=> f for f in fs]...)` if functions " *
                     "in `fs` have unique names.", :aggregate)
    end
    return df
end

function aggregate(gd::GroupedDataFrame, f::Any; sort::Bool=false)
    df = combine(gd, valuecols(gd) .=> [f])
    if sort
        Base.depwarn("`aggregate(gd, f, sort=true)` is deprecated. Instead use" *
                     " `df = combine(gd, valuecols(gd) .=> f); " *
                     "sort!(df, Not(groupcols(gd)))`.",
                     :aggregate)
        sort!(df, names(df, Not(groupcols(gd))))
    else
        Base.depwarn("`aggregate(gd, f)` is deprecated. Instead use" *
                     " `combine(gd, valuecols(gd) .=> f)`",
                     :aggregate)
    end
    return df
end

function aggregate(gd::GroupedDataFrame, fs::AbstractVector; sort::Bool=false)
    df = hcat([combine(gd, valuecols(gd) .=> [f], keepkeys=i==1) for (i, f) in enumerate(fs)]...,
                       makeunique=true)
    if sort
        Base.depwarn("`aggregate(gd, fs, sort=true)` is deprecated. Instead use " *
                     "`df = combine(gd, [names(gd) .=> f for f in fs]...); " *
                     "sort!(df, names(df, Not(groupcols(gd))))`" *
                     " if functions in `fs` have unique names.", :aggregate)
        sort!(df, Not(groupcols(gd)))
    else
        Base.depwarn("`aggregate(gd, fs)` is deprecated. Instead" *
                     " use `combine(gd, [names(gd) .=> f for f in fs]...)`" *
                     " if functions in `fs` have unique names.", :aggregate)
    end
    return df
end

function aggregate(d::AbstractDataFrame, cols, f::Any;
                   sort::Bool=false, skipmissing::Bool=false)
    Base.depwarn("`aggregate(d, cols, f, sort=$sort, skipmissing=$skipmissing)` " *
                 "is deprecated. Instead use " *
                 "combine(groupby(d, cols, sort=$sort, skipmissing=$skipmissing)," *
                 " names(d, Not(cols)) .=> f)`",
                 :aggregate)
    gd = groupby(d, cols, sort=sort, skipmissing=skipmissing)
    df = combine(gd, valuecols(gd) .=> [f])
    return df
end

function aggregate(d::AbstractDataFrame, cols, fs::AbstractVector;
                   sort::Bool=false, skipmissing::Bool=false)
    Base.depwarn("`aggregate(d, cols, fs, sort=$sort, skipmissing=$skipmissing)` " *
                 " is deprecated. Instead use " *
                 "combine(groupby(d, cols, sort=$sort, skipmissing=$skipmissing), "*
                 "[names(d, Not(cols)) .=> f for f in fs]...)`" *
                 " if functions in `fs` have unique names.", :aggregate)
    gd = groupby(d, cols, sort=sort, skipmissing=skipmissing)
    df = hcat([combine(gd, valuecols(gd) .=> [f], keepkeys=i==1) for (i, f) in enumerate(fs)]...,
                       makeunique=true)
    return df
end

@deprecate deleterows!(df::DataFrame, inds) delete!(df, inds)

@deprecate by(f::Base.Callable, d::AbstractDataFrame, cols::Any;
    sort::Bool=false, skipmissing::Bool=false) combine(f,
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(f::Pair{<:ColumnIndex}, d::AbstractDataFrame, cols::Any;
    sort::Bool=false, skipmissing::Bool=false) combine(f,
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(f::Pair, d::AbstractDataFrame, cols::Any;
    sort::Bool=false, skipmissing::Bool=false) combine(AsTable(first(f)) => last(f),
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(d::AbstractDataFrame, cols::Any, f::Base.Callable;
    sort::Bool=false, skipmissing::Bool=false) combine(f,
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(d::AbstractDataFrame, cols::Any, f::Pair{<:ColumnIndex};
    sort::Bool=false, skipmissing::Bool=false) combine(f,
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(d::AbstractDataFrame, cols::Any, f::Pair;
    sort::Bool=false, skipmissing::Bool=false) combine(AsTable(first(f)) => last(f),
    groupby(d, cols, sort=sort, skipmissing=skipmissing))
@deprecate by(d::AbstractDataFrame, cols::Any, f::Pair{<:ColumnIndex}...;
    sort::Bool=false, skipmissing::Bool=false) combine(groupby(d, cols, sort=sort, skipmissing=skipmissing), f...)
@deprecate by(d::AbstractDataFrame, cols::Any, f::Pair...;
    sort::Bool=false, skipmissing::Bool=false) combine(groupby(d, cols, sort=sort, skipmissing=skipmissing),
    [(col isa ColumnIndex ? col : AsTable(col)) => fun for (col, fun) in f]...)
@deprecate by(d::AbstractDataFrame, cols::Any;
    sort::Bool=false, skipmissing::Bool=false, f...) combine(groupby(d, cols,
    sort=sort, skipmissing=skipmissing),
    [(in_col isa ColumnIndex ? in_col : AsTable(in_col)) => fun => out_col for (out_col, (in_col, fun)) in f]...)

import Base: map
@deprecate map(f::Base.Callable, gd::GroupedDataFrame) combine(f, gd, ungroup=false)
@deprecate map(f::Pair{<:ColumnIndex}, gd::GroupedDataFrame) combine(f, gd, ungroup=false)
@deprecate map(f::Pair, gd::GroupedDataFrame) combine(AsTable(first(f)) => last(f), gd, ungroup=false)


text:
@inline function Base.iterate(graph::PharmGraph)
    iteration = Base.iterate(MetaGraphs.vertices(graph.g))
    vertex_integer, next_state = iteration
    class = _vertex_integer_to_class(graph, vertex_integer)
    result = (class, next_state)
    return result
end

@inline function Base.iterate(graph::PharmGraph, state)
    iteration = Base.iterate(MetaGraphs.vertices(graph.g), state)
    iteration === nothing && return nothing
    vertex_integer, next_state = iteration
    class = _vertex_integer_to_class(graph, vertex_integer)
    result = (class, next_state)
    return result
end

@inline function Base.IteratorSize(::Type{PharmGraph})
    return Base.HasLength()
end

@inline function Base.IteratorEltype(::Type{PharmGraph})
    return Base.HasEltype()
end

@inline function Base.eltype(graph::PharmGraph)
    return PharmClass
end

@inline function Base.length(graph::PharmGraph)
    return MetaGraphs.nv(graph.g)
end


text:
function nlp_cvx_501_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - n dimentional nlp
    # Variants
    #   010 - binding constraint (all variables non-zero)
    #   011 - extra sqrt

    function nd_shpere(n = 2)
        model = Model(optimizer)

        @variable(model, vars[1:n])

        @objective(model, Min, sum(-x for x in vars))
        @NLconstraint(model, sum(x^2 for x in vars) <= 1.0)

        optimize!(model)

        #println(getobjectivevalue(model))

        check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
        check_objective(model, -n / sqrt(n), tol = objective_tol)
        return check_solution(
            vars,
            [1 / sqrt(n) for i in vars],
            tol = primal_tol,
        )
    end

    for n in 1:20
        nd_shpere(n)
    end
end


text:
module Bosonic
using SparseArrays
using LinearAlgebra

include("operators_fixed.jl")
include("states_fixed.jl")
include("correlations.jl")

export basis_m, bdb, bbd, Op_fixed, dim, nume, le, basis, b, bd, state_form
export State_fixed, st, ope, nume, dim, basis, typ, rhosp
end


text: eye(N) = Matrix{Float64}(I, N, N)


text:
# Test multihypo computation assembly

using Distributions
using DocStringExtensions
using Base: Test

using IncrementalInference

@testset "test IncrementalInference.assembleHypothesesElements! without multihypothesis..." begin

s2_1_gt1 = 1:2
s2_1_gt2 = (1:20,Int[])
s2_1_gt3 = [(1,1:2); (2,Int[])]
s2_1_gt4 = Int[]
s2_1 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 1, 2 )
@test sum(s2_1_gt1 .- s2_1[1]) == 0
@test sum(s2_1_gt2[1] .- s2_1[2][1]) == 0
@test sum(s2_1_gt2[2] .- s2_1[2][2]) == 0
@test s2_1_gt3[1][1] == s2_1[3][1][1]
@test sum(s2_1_gt3[1][2] .- s2_1[3][1][2]) == 0
@test s2_1_gt3[2][1] == s2_1[3][2][1]
@test sum(s2_1_gt3[2][2] .- s2_1[3][2][2]) == 0
@test sum(s2_1_gt4 .- s2_1[4])  == 0

s2_2_gt1 = 1:2
s2_2_gt2 = (1:20,Int[])
s2_2_gt3 = [(1,1:2); (2,Int[])]
s2_2_gt4 = Int[]

s2_2 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 2, 2 )

@test sum(s2_2_gt1 .- s2_2[1]) == 0
@test sum(s2_2_gt2[1] .- s2_2[2][1]) == 0
@test sum(s2_2_gt2[2] .- s2_2[2][2]) == 0
@test s2_2_gt3[1][1] == s2_2[3][1][1]
@test sum(s2_2_gt3[1][2] .- s2_2[3][1][2]) == 0
@test s2_2_gt3[2][1] == s2_2[3][2][1]
@test sum(s2_2_gt3[2][2] .- s2_2[3][2][2]) == 0
@test sum(s2_2_gt4 .- s2_2[4])  == 0

end

@testset "test IncrementalInference.assembleHypothesesElements! with bi-modality..." begin

s3_1_gt1 = [1]
s3_1_gt2 = (0,5,5,20)
s3_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3])]
s3_1_gt4 = 20

s3_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 20, 1, 3 )

@test sum(s3_1_gt1 - s3_1[1]) == 0
@test sum(s3_1_gt2[1] .- s3_1[2][1]) == 0
@test length(s3_1[2][2]) > s3_1_gt2[2]
@test length(s3_1[2][3]) > s3_1_gt2[3]
@test length(s3_1[2][2]) + length(s3_1[2][3]) == s3_1_gt2[4]
@test s3_1_gt3[1][1] == s3_1[3][1][1]
@test s3_1_gt3[2][1] == s3_1[3][2][1]
@test s3_1_gt3[3][1] == s3_1[3][3][1]
@test sum(s3_1_gt3[1][2] .- s3_1[3][1][2]) == 0
@test sum(s3_1_gt3[2][2] .- s3_1[3][2][2]) == 0
@test sum(s3_1_gt3[3][2] .- s3_1[3][3][2]) == 0

@test sum(s3_1[4] .== 2) > s3_1_gt2[2]
@test sum(s3_1[4] .== 3) > s3_1_gt2[3]

@test sum( [1:20;][s3_1[4] .== 2] .== s3_1[2][2] ) == length(s3_1[2][2])
@test sum( [1:20;][s3_1[4] .== 3] .== s3_1[2][3] ) == length(s3_1[2][3])
@test length(s3_1[4]) == s3_1_gt4

s3_2_gt1 = [1]
s3_2_gt2 = (0,5,5,20)
s3_2_gt3 = [(1,Int[1;2]); (2,Int[2;3]); (3,Int[2;3])]
s3_2_gt4 = 20

s3_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 20, 2, 3 )

@test sum(s3_2_gt1 - s3_2[1]) == 0
@test sum(s3_2_gt2[1] .- s3_2[2][1]) == 0
@test length(s3_2[2][2]) > s3_2_gt2[2]
@test length(s3_2[2][3]) > s3_2_gt2[3]
@test length(s3_2[2][2]) + length(s3_2[2][3]) == s3_2_gt2[4]
@test s3_2_gt3[1][1] == s3_2[3][1][1]
@test s3_2_gt3[2][1] == s3_2[3][2][1]
@test s3_2_gt3[3][1] == s3_2[3][3][1]
@test sum(s3_2_gt3[1][2] .- s3_2[3][1][2]) == 0
@test sum(s3_2_gt3[2][2] .- s3_2[3][2][2]) == 0
@test sum(s3_2_gt3[3][2] .- s3_2[3][3][2]) == 0

@test sum(s3_2[4] .== 2) > s3_2_gt2[2]
@test sum(s3_2[4] .== 3) > s3_2_gt2[3]

@test sum( [1:20;][s3_2[4] .== 2] .== s3_2[2][2] ) == length(s3_2[2][2])
@test sum( [1:20;][s3_2[4] .== 3] .== s3_2[2][3] ) == length(s3_2[2][3])
@test length(s3_2[4]) == s3_2_gt4

s3_3_gt1 = [1]
s3_3_gt2 = (0,5,5,20)
s3_3_gt3 = [(1,Int[1;3]); (2,Int[2;3]); (3,Int[2;3])]
s3_3_gt4 = 20

s3_3 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 20, 3, 3 )

@test sum(s3_3_gt1 - s3_3[1]) == 0
@test sum(s3_3_gt2[1] .- s3_3[2][1]) == 0
@test length(s3_3[2][2]) > s3_3_gt2[2]
@test length(s3_3[2][3]) > s3_3_gt2[3]
@test length(s3_3[2][2]) + length(s3_3[2][3]) == s3_3_gt2[4]
@test s3_3_gt3[1][1] == s3_3[3][1][1]
@test s3_3_gt3[2][1] == s3_3[3][2][1]
@test s3_3_gt3[3][1] == s3_3[3][3][1]
@test sum(s3_3_gt3[1][2] .- s3_3[3][1][2]) == 0
@test sum(s3_3_gt3[2][2] .- s3_3[3][2][2]) == 0
@test sum(s3_3_gt3[3][2] .- s3_3[3][3][2]) == 0

@test sum(s3_3[4] .== 2) > s3_3_gt2[2]
@test sum(s3_3[4] .== 3) > s3_3_gt2[3]

@test sum( [1:20;][s3_3[4] .== 2] .== s3_3[2][2] ) == length(s3_3[2][2])
@test sum( [1:20;][s3_3[4] .== 3] .== s3_3[2][3] ) == length(s3_3[2][3])
@test length(s3_3[4]) == s3_3_gt4

end

@testset "test IncrementalInference.assembleHypothesesElements! with tri-modality..." begin

# s3_4_gt1 = [1]
# s3_4_gt2 = (0,5,5,20)
# s3_4_gt3 = [(1,Int[1;3]); (2,Int[2;3]); (3,Int[2;3])]
# s3_4_gt4 = 20

s4_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.3;0.3;0.4]), 20, 2, 4 )

end


text:
# Contains evaluation functions
# ----------------------------------
"""
	 Evaluates the FIS

	eval_fis(fis, input_values,defuzz_method = "WTAV")

	 Parameters
	 ----------
	 `fis` is the inference system to evaluate
	 `input_values` is a Vector of inputs
	 `defuzz_method` is the method for defuzzification, see defuzz function definition

"""
function eval_fis(fis::FISMamdani, input_values::Vector{<:AbstractFloat},defuzz_method = "WTAV")

	firing_strengths = AbstractFloat[]
	for rule in fis.rules
		tmp_strengths = AbstractFloat[]
		for i in 1:length(rule.input_mf_names)
			push!(tmp_strengths, fis.input_mfs_dicts[i][rule.input_mf_names[i]].eval(input_values[i]))
		end
			push!(firing_strengths, firing(tmp_strengths,rule.firing_method ))
	end
	defuzz(firing_strengths, fis, defuzz_method)

end

"""
	 Evaluates the FIS

	eval_fis(fis, input_values)

	 Parameters
	 ----------
	 `fis` is the inference system to evaluate
	 `input_values` is a Vector of inputs

"""
function eval_fis(fis::FISSugeno, input_values::Vector{<:AbstractFloat})

	firing_strengths = AbstractFloat[]
	for rule in fis.rules
		tmp_strengths = AbstractFloat[]
		for i in 1:length(rule.input_mf_names)
			push!(tmp_strengths, fis.input_mfs_dicts[i][rule.input_mf_names[i]].eval(input_values[i]))
		end
			push!(firing_strengths, firing(tmp_strengths,rule.firing_method ))
	end

	push!(input_values, 1)
	n_firing_strengths = firing_strengths / sum(firing_strengths)
	out = 0.0
	for i = 1:length(fis.rules)
		out += n_firing_strengths[i] * (fis.rules[i].output_mf' * input_values)[1]
	end
	pop!(input_values)
	return out

end

"""
	 Defuzzifies the output using the given firing strengths

	defuzz(firing_strengths, rules, output_mfs_dict, defuzz_method)

	 Parameters
	 ----------
	 `firing_strengths` is a Vector of firing strengths
	 		one for each output membership function
	 `rules` is a Vector of Rule
	 `output_mfs_dict` is a Dict of output membership functions
	 `defuzz_method` is the method for defuzzification
	 		"MOM" - Mean of Maximum
	 		"WTAV" - Weighted Average
			"COG" - Center of Gravity
"""

include("Utils.jl")

function defuzz(firing_strengths::Vector{AbstractFloat}, fis::FISMamdani, defuzz_method::AbstractString)
	rules = fis.rules
	output_mfs_dict = fis.output_mfs_dict
	domain = fis.output_domain
	if defuzz_method == "MOM"
		max_firing_index = argmax(firing_strengths)
		max_fired_mf_name = rules[max_firing_index].output_mf
		output_mfs_dict[max_fired_mf_name].mean_at(maximum(firing_strengths))
	elseif defuzz_method == "WTAV"
		mean_vec = AbstractFloat[]
		for i in 1:length(rules)
			push!(mean_vec, output_mfs_dict[rules[i].output_mf].mean_at(firing_strengths[i]))
		end
		(mean_vec' * firing_strengths)[1] / sum(firing_strengths)
	elseif defuzz_method == "COG"
		cutted = []
		for i in 1:length(rules)
			push!(cutted, output_mfs_dict[rules[i].output_mf].cut(firing_strengths[i]))
		end
		agg(x) = maximum(map(y -> y.eval(x), cutted))
		agg_aux(x) = x*agg(x)
		n = sum(map(y -> y.get_n(1e-5), cutted))
		return trapz(agg_aux, domain[1], domain[2], n)/trapz(agg, domain[1], domain[2], n)
	end

end

"""
	 Docs goes here
"""
function firing(tmp_strengths::Vector{<:AbstractFloat},firing_method::AbstractString)
	if firing_method == "MIN"
		return	minimum_value(tmp_strengths)
	elseif firing_method == "A-PROD"
		return algebraic_product(tmp_strengths)
	elseif firing_method == "B-DIF"
		return bounded_difference(tmp_strengths)
	elseif firing_method == "D-PROD"
		return drastic_product(tmp_strengths)
	elseif firing_method == "E-PROD"
		return einstein_product(tmp_strengths)
	elseif firing_method == "H-PROD"
		return hamacher_product(tmp_strengths)
	elseif firing_method == "MAX"
		return	maximum_value(tmp_strengths)
	elseif firing_method == "A-SUM"
		return algebraic_sum(tmp_strengths)
	elseif firing_method == "B-SUM"
		return bounded_sum(tmp_strengths)
	elseif firing_method == "D-SUM"
		return drastic_sum(tmp_strengths)
	elseif firing_method == "E-SUM"
		return einstein_sum(tmp_strengths)
	elseif firing_method == "H-SUM"
		return hamacher_sum(tmp_strengths)
	end
end


text:
"""
A lightweitht julia wrapper for [WORLD](https://github.com/mmorise/World),
a high-quality speech analysis, manipulation and synthesis system.
WORLD provides a way to decompose a speech signal into:

- Fundamental frequency (F0)
- spectral envelope
- aperiodicity

and re-synthesize a speech signal from these paramters. Please see the project
page for more details on the WORLD.

!!! note
    WORLD.jl is based on a fork of WORLD ([r9y9/World-cmake](https://github.com/r9y9/World-cmake)).

[https://github.com/r9y9/WORLD.jl](https://github.com/r9y9/WORLD.jl)

## Usage

In the following examples, suppose `x::Vector{Float64}` is a input monoral
speech signal like:

![](assets/x.png)

### F0 estimation

#### Harvest

```julia
opt = HarvestOption(71.0, 800.0, period)
f0, timeaxis = harvest(x, fs, opt)
```

![](assets/f0_by_harvest.png)

#### Dio

```julia
opt = DioOption(f0floor=71.0, f0ceil=800.0, channels_in_octave=2.0,
        period=period, speed=1)
f0, timeaxis = dio(x, fs, opt)
```

![](assets/f0_by_dio.png)

#### StoneMask

```julia
f0 = stonemask(x, fs, timeaxis, f0)
```

![](assets/f0_refinement.png)

### Spectral envelope estimation by CheapTrick

```julia
spectrogram = cheaptrick(x, fs, timeaxis, f0)
```

![](assets/envelope_by_cheaptrick.png)

### Aperiodicity ratio estimation by D4C

```julia
aperiodicity = d4c(x, fs, timeaxis, f0)
```

![](assets/aperiodicity_by_d4c.png)

### Synthesis

```julia
y = synthesis(f0, spectrogram, aperiodicity, period, fs, length(x))
```

![](assets/synthesis.png)

### Compact speech parameterization

Raw spectrum envelope and aperiodicity spectrum are relatively high dimentional
(offen more than 513 or 1025) so one might want to get more compact
representation. To handle this situation, WORLD provides coding/decoding APIs
for spectrum envelope and aperiodicity. Additionally, WORLD.jl provides
conversions from spectrum envelope to mel-cepstrum and vice versa. You can choose
any of coding/decoding APIs depends on your purpose.

#### spectrum envelope to mel-cepstrum

```julia
mc = sp2mc(spectrogram, order, α) # e.g. order=40, α=0.41
```

where `order` is the order of mel-cepstrum (except for 0th) and α is a frequency
warping parameter.

![](assets/melcepstrum.png)

#### mel-cepstrum to spectrum envelope

```julia
approximate_spectrogram = mc2sp(mc, α, get_fftsize_for_cheaptrick(fs))
```

![](assets/envelope_reconstructed_from_melcepstrum.png)

#### Code aperiodicity

```julia
coded_aperiodicity = code_aperiodicity(aperiodicity, fs)
```

![](assets/coded_aperiodicity.png)

#### Decode aperiodicity

```julia
decoded_aperiodicity = decode_aperiodicity(coded_aperiodicity, fs)
```

![](assets/decoded_aperiodicity.png)

For the complete code of visualizations shown above, please check
[the IJulia notebook](http://nbviewer.jupyter.org/github/r9y9/WORLD.jl/blob/master/docs/src/assets/WORLD-demo.ipynb).

## Exports

$(EXPORTS)
"""
module WORLD

using DocStringExtensions

export DioOption, HarvestOption, CheapTrickOption, D4COption, dio, harvest,
    stonemask, cheaptrick, d4c, synthesis, get_fftsize_for_cheaptrick,
    interp1!, interp1, sp2mc, mc2sp,
    get_number_of_aperiodicities, code_aperiodicity, decode_aperiodicity,
    code_spectral_envelope, decode_spectral_envelope

### Binary dependency loading ###
deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("WORLD not properly installed. Please run Pkg.build(\"WORLD\")")
end

### Check binary dependency compatibility ###
# This should prevent segfaults when you are using WORLD on system paths (e.g.
# /usr/local/lib/libworld.so) and accidentaly change the verion of WORLD.
try
    # function GetWORLDVersion was added in v0.2.1-2
    versionstr = unsafe_string(convert(Ptr{UInt8},
        ccall((:GetWORLDVersion, libworld), Ptr{Cchar}, ())))
    global const version = VersionNumber(versionstr)
catch e
    try
        ccall((:GetFFTSizeForStar, libworld), Cint, (Cint,), 48000)
        global const version = v"0.1.4"
    catch e
        global const version = v"0.2.0"
    end
end

# v0.2.1-4 breaks API of CheapTrick and D4C, and it's hard to support both old
# and new WORLD library, so throws an error when older version of WORLD is
# detected.
if version < v"0.3.0"
    error("""Binary dependency incompatibility

          WORLD version 0.3.0 or later must be required ($(version) is detected).
          If you have WORLD installed on your system paths, please get the latest
          stable WORLD and install it, and then re-build WORLD.jl with:

          julia>  Pkg.build(\"WORLD\")

          If not, simpily re-building WORLD.jl with the above command. This
          should install the proper version of WORLD.

          If you want to use WORLD.jl with older WORLD, please downgrade WORLD.jl
          for your need. E.g.

          julia> cd(Pkg.dir(\"WORLD\"))
          shell> git checkout v0.1.3 # that supprts v0.2.0 to v0.2.1-3.
         """)
end

for name in ["bridge", "codec", "mcep", "deprecated"]
    include(string(name, ".jl"))
end

end # module WORLD


text:
function tarai(x::Int64, y::Int64, z::Int64)
  if x <= y
    return y
  else
    return tarai(tarai(x - 1, y, z), tarai(y - 1, z, x), tarai(z - 1, x, y))
  end
end

tarai2(x, y, z) = if x <= y return y else return tarai(tarai(x-1, y, z), tarai(y-1, z, x), tarai(z-1, x, y)) end

println(tarai(12, 6, 0))
println(tarai2(12, 6, 0))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
using Oceananigans.Architectures: arch_array, architecture

Lz(c::AbstractDataField) = c.grid.Lz
Nz(c::AbstractDataField) = c.grid.Nz
Δz(g::RegularRectilinearGrid, i::Int) = g.Δz

function integrate_range(c, cgrid, i₁::Int, i₂::Int)
    total = 0
    for i = i₁:i₂
        @inbounds total += c[i] * Δz(cgrid, i)
    end
    return total
end

function integral(c::AbstractDataField, z₋, z₊=0)

    @assert z₊ > c.grid.zF[1] "Integration region lies outside the domain."
    @assert z₊ > z₋ "Invalid integration range: upper limit greater than lower limit."

    # Find region bounded by the face ≤ z₊ and the face ≤ z₁
    i₁ = searchsortedfirst(c.grid.zF, z₋) - 1
    i₂ = searchsortedfirst(c.grid.zF, z₊) - 1

    if i₂ ≠ i₁
        # Calculate interior integral, recalling that the
        # top interior cell has index i₂-2.
        total = integrate_range(c, c.grid, i₁+1, i₂-1)

        # Add contribution to integral from fractional bottom part,
        # if that region is a part of the grid.
        if i₁ > 0
            total += c[i₁] * (c.grid.zF[i₁+1] - z₋)
        end

        # Add contribution to integral from fractional top part
        total += c[i₂] * (z₊ - c.grid.zF[i₂])
    else
        total = c[i₁] * (z₊ - z₋)
    end

    return total
end

function integral(c, cgrid, z₋, z₊=0)

    @assert z₊ > cgrid.zF[1] "Integration region lies outside the domain."
    @assert z₊ > z₋ "Invalid integration range: upper limit greater than lower limit."

    # Find region bounded by the face ≤ z₊ and the face ≤ z₁
    i₁ = searchsortedfirst(cgrid.zF, z₋) - 1
    i₂ = searchsortedfirst(cgrid.zF, z₊) - 1

    if i₂ ≠ i₁
        # Calculate interior integral, recalling that the
        # top interior cell has index i₂-2.
        total = integrate_range(c, cgrid, i₁+1, i₂-1)

        # Add contribution to integral from fractional bottom part,
        # if that region is a part of the grid.
        if i₁ > 0
            total += c[i₁] * (cgrid.zF[i₁+1] - z₋)
        end

        # Add contribution to integral from fractional top part
        total += c[i₂] * (z₊ - cgrid.zF[i₂])
    else
        total = c[i₁] * (z₊ - z₋)
    end

    return total
end

# Set interior of field `c` to values of `data`
function set!(c::AbstractField, data::AbstractArray)

    arch = architecture(c)

    # Reshape `data` to the size of `c`'s interior
    d = arch_array(arch, reshape(data, size(c)))

    # Sets the interior of field `c` to values of `data`
    c .= d

end

horizontal_size(grid) = (grid.Nx, grid.Ny)
extent(grid) = (grid.Lx, grid.Ly, grid.Lz)

# Set two fields to one another... some shenanigans
#
_set_similar_fields!(c::AbstractField{Ac, G}, d::AbstractField{Ad, G}) where {Ac, Ad, G} =
    c.data .= convert(typeof(c.data), d.data)

function interp_and_set!(c1::AbstractField{A1, G1}, c2::AbstractField{A2, G2}) where {A1, A2, G1, G2}

    grid1 = c1.grid
    grid2 = c2.grid

    @assert extent(grid1) == extent(grid2) "Physical domains differ between the two fields."

    for j in 1:grid1.Ny, i in 1:grid1.Nx
        for k in 1:grid1.Nz
            @inbounds c1[i,j,k] = integral(c2[i,j,:], grid2, grid1.zF[k], grid1.zF[k+1]) / Δz(grid1, i)
        end
    end

    return nothing
end

"""
    set!(c::AbstractField{Ac, G}, d::AbstractField{Ad, G}) where {Ac, Ad, G}

Set the data of field `c` to the data of field `d`, adjusted to field `c`'s grid.

The columns are assumed to be independent and thus the fields must have the same
horizontal resolution. This implementation does not accommodate 3D grids with
dependent columns.
"""
function set!(c::AbstractField{Ac, G}, d::AbstractField{Ad, G}) where {Ac, Ad, G}

    s1 = horizontal_size(c.grid)
    s2 = horizontal_size(d.grid)
    @assert s1 == s2 "Field grids have a different number of columns."

    if s1 != (1, 1)
        @assert c.grid isa OneDimensionalEnsembleGrid && d.grid isa OneDimensionalEnsembleGrid "Field has dependent columns."
    end

    if Lz(c) == Lz(d) && Nz(c) == Nz(d)
        return _set_similar_fields!(c, d)
    else
        return interp_and_set!(c, d)
    end

end


text:
using MechanicalSketch
import MechanicalSketch: °, background, sethue, O, PT, EM, finish
import MechanicalSketch: color_from_palette, color_with_lumin
import MechanicalSketch: m, °, kN, mm, PALETTE
import MechanicalSketch: m, m², m³, s, N, kN, kPa, g, kg
import MechanicalSketch: arrow, drawcart, circle, foil_draw
import MechanicalSketch: line, showpower, empty_figure
import MechanicalSketch: rope_pos_tension, setline, kW

let
BACKCOLOR = PALETTE[8];
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_13.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
end
restart()

#  Physical constants

"Rope shape coefficient at Reynold's number range. It is
applied to the velocity components in a plane normal to the section"
Cs = 1.4

"Density air"
ρ_air = 1.322kg/m³

#  Scale-general constants

"Rope minimum breaking stress"
σ_MB = 41.2kN / (π / 4  * (6mm)^2 ) |> N/mm^2
# or
σ_MB = kN / (π / 4  * (6mm)^2 ) |> N/mm^2
# or
σ_MB = 5768kN / (π / 4  * (96mm)^2 ) |> N/mm^2
# or
σ_MB = 407.1kN / (π / 4  * (20mm)^2 ) |> N/mm^2

"Wing loading, dimensioning, from Boeing 737"
p_wing_dim = 690kg*g/m² |> kPa

# Model

"Rope safety factor"
Sf_s = 1.75

"Cart position"
p = (3m, -9m)

"Wind speed scalar"
v_w = 10m/s

"Wind velocity vector"
w = (0m/s, v_w)

"Cart velocity"
v_cart = -10m/s

"Generated wind speed vector, same for kite and cart since the model is quasi static"
w_generated = (-v_cart, 0m/s)

"Relative wind velocity vector"
w_rel = w + w_generated

"Relative wind velocity"
v = hypot(w_rel)

"Relative wind direction to screen x-axis"
α_w_rel = atan(w_rel[2], w_rel[1]) |> °

"Diameter rope"
diameter_rope = 6mm

"Wing area, kite"
A_k = σ_MB * π / 4 * diameter_rope^2 / Sf_s / p_wing_dim |> m²

"Aspect ratio kite"
AR = 7

"Average chord length, kite"
l_c = A_k / AR / m

"Optimal angle of attack to wind (positive clockwise by convention, defined by chord, not by zero lift)"
α_a = 8°

"Chord angle to global x, positive around z"
α_chord = α_w_rel - α_a

"Lift coefficient at optimal angle of attack, considering aspect ratio"
Cl = 1.25

"Lift to drag, optimal angle"
Cl_Cd = 14

"Foil pressure on wing projected on chord-span area"
pl = 0.5 * Cl * ρ_air * v^2 |> kPa

"Foil shear stress on wing projected on chord-span area, including body drag"
pd = pl / Cl_Cd

"Foil force component in chord coordinates"
Lift, Drag = A_k .* (pl, pd) .|> kN

"Foil force vector in chord coordinates"
Fk = (Drag, Lift)

"Change the diameter of rope from the sensible values"
diameter_rope = 20 * diameter_rope |> m

"Length of rope"
Ls = 15m

"Number of rope sections (one less than rope nodes)"
Ns = 100

"Rope positions from foil origo, oriented in global axes"
psf, Ts = rope_pos_tension(Ls, Ns, diameter_rope, Cs, ρ_air, v, α_w_rel, Fk, α_chord);

"Foil position"
p_foil = p - psf[end]

"Rope positions, global axes"
ps =  map(p-> p + p_foil, psf);

"Tangent at rope end"
Δ = ps[end-1] - ps[end]

"Unit vector, rope end"
eΔ = Δ ./ hypot(Δ)

"Force vector, rope end"
Fc = eΔ .* Ts[end]

"Braking force scalar, cart, positive is actual braking"
F_brake = -Fc[1]

"Braking power, cart"
power = -F_brake * v_cart |> kW

"Horizontal force, foil"
Fk_x =  -Fk[2] * sin(α_chord) + Fk[1] *cos(α_chord)

"Working power, kite"
power_kite = Fk_x * v_cart |> kW

"Horizontal force, rope"
Fs =  -Fk_x - F_brake |> N

"Braking power, rope"
power_rope = -Fs  * v_cart |> kW

###########################
# Drawing only
###########################

#"Cart position"
circle(O + p, 0.1m, :stroke)
drawcart(p = O +p)

#"Generated wind speed vector, same for kite and cart since the model is quasi static"
sethue(color_with_lumin(color_from_palette("green"), 85))
arrow(O + p + 0.5w_generated + 2w,
    w_generated,
    backgroundcolor = BACKCOLOR,
    labellength = true)

#"Relative wind velocity"
sethue(sethue(PALETTE[1]))
arrow(O + p + 0.5w_rel , w_rel,
    backgroundcolor = BACKCOLOR,
    labellength = true)

# Draw the foil
sethue(color_from_palette("blue"))
foil_draw(O + p_foil, α = α_chord, rel_scale = 1, l = l_c, backgroundcolor = BACKCOLOR);

# And the force on foil
sethue(color_from_palette("indigo"))
arrow(O + p_foil, Fk, α = α_chord,
    backgroundcolor = BACKCOLOR,
    labellength = true)

# Draw the relative wind in front of the kite
sethue(PALETTE[1])
arrow(O + p_foil - 1.5w_rel , w_rel,
    backgroundcolor = BACKCOLOR,
    labellength = false)

#"Rope positions, global axes"
setline(diameter_rope)
for i = 1:Ns
    line(O + ps[i], O + ps[i+1], :stroke)
end
setline(0.5PT)

#"Force vector, rope end"
sethue(PALETTE[5])
arrow(O + p, Fc,
    backgroundcolor = BACKCOLOR,
    labellength = true)

#"Braking force scalar, cart, positive is actual braking"
sethue(PALETTE[6])
arrow(O + p, F_brake,
    backgroundcolor = BACKCOLOR,
    labellength = true)

#"Braking power, cart"
sethue(color_from_palette("yellow"))
showpower(O + p + (10EM, 0), power; backgroundcolor = BACKCOLOR)

#"Horizontal force, foil"
sethue(PALETTE[6])
arrow(O + p_foil, Fk_x,
    backgroundcolor = BACKCOLOR,
    labellength = true)

#"Working power, kite"
sethue(color_from_palette("yellow"))
showpower(O + p_foil + (10EM, 0), power_kite; backgroundcolor = BACKCOLOR)

#"Horizontal force, rope"
sethue(color_with_lumin(PALETTE[6], 10))
arrow(O + 0.5p_foil + 0.5p, Fs,
    backgroundcolor = BACKCOLOR,
    labellength = true)

#"Braking power, rope"
sethue(color_from_palette("yellow"))
showpower(O + 0.5p_foil + 0.5p + (10EM, 0), power_rope; backgroundcolor = BACKCOLOR)

finish()
nothing
end


text:
using ACEbase
using Test

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