text:
include("liftedfrombindeps.jl")

import WinReg

# define current version:
mosekver = open(joinpath(dirname(@__FILE__),"..","MOSEKVER")) do f readline(f) end
mskvmajor,mskvminor = split(mosekver,".")

mskplatform,distroext =
  if Sys.ARCH == :i386 || Sys.ARCH == :i686
    if     Sys.islinux()   "linux32x86",  ".tar.bz2"
    elseif Sys.iswindows() "win32x86",  ".zip"
    else   error("Platform not supported")
    end
  elseif Sys.ARCH == :x86_64
    if     Sys.islinux()   "linux64x86",".tar.bz2"
    elseif Sys.isapple()   "osx64x86",  ".tar.bz2"
    elseif Sys.iswindows() "win64x86",  ".zip"
    else   error("Platform not supported")
    end
  elseif Sys.ARCH == :aarch64
    if     Sys.islinux()   "linuxaarch64",".tar.bz2"
    elseif Sys.isapple()   "osxaarch64",  ".tar.bz2"
    else   error("Platform not supported")
    end
  else
    error("Platform not supported")
  end

bindepsdir = dirname(@__FILE__)

function getregistrykey(base::UInt32,path::String,key::String)
    WinReg.querykey(base,path,key)
end

function hasregistrykey(base::UInt32,path::String,key::String)
    if ! Sys.iswindows()
        false
    else
        try
            getregistrykey(base,path,key)
	    true
        catch err
            false
        end
    end
end

function findlibs(path::AbstractString,mskvmajor::AbstractString,mskvminor::AbstractString)
    moseklib =
        if Sys.ARCH == :i386 || Sys.ARCH == :i686
            if     Sys.iswindows() "mosek$(mskvmajor)_$(mskvminor).dll"
            else   error("Platform not supported")
            end
        elseif Sys.ARCH == :x86_64
            if     Sys.islinux()   "libmosek64.so.$(mskvmajor).$(mskvminor)"
            elseif Sys.isapple()   "libmosek64.$(mskvmajor).$(mskvminor).dylib"
            elseif Sys.iswindows() "mosek64_$(mskvmajor)_$(mskvminor).dll"
            else   error("Platform not supported")
            end
        else
            error("Platform not supported")
        end

    let moseklibpath = joinpath(path,moseklib)

        if ! isfile(moseklibpath)
            error("Library '$moseklib' not found ($path)")
        else
            moseklibpath
        end
    end
end

function versionFromBindir(bindir ::AbstractString)
    try
        mosekbin = if Sys.iswindows() "mosek.exe" else "mosek" end
        txt = read(`$bindir/$mosekbin`,String)
        @info(txt)
        m = match(r"\s*MOSEK Version ([0-9]+\.[0-9]+\.[0-9])",txt)
        if m == nothing
            return nothing
        else
            return m.captures[1]
        end
    catch
        return nothing
    end
end

function bindirIsCurrentVersion(bindir)
    ver = versionFromBindir(bindir)
    if ver != nothing
        @info("Got version: $ver, expected version: $mskvmajor.$mskvminor")
        ver = split(ver,".")
        return ver[1] == mskvmajor && ver[2] == mskvminor
    else
        return false
    end
end

function collect_output(cmd::Cmd)
    out = Pipe()
    err = Pipe()

    process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))

    stdout = @async (read(out))
    stderr = @async (read(err))

    wait(process)
    close(out.in)
    close(err.in)
    (
        process.exitcode,
        fetch(stdout),
        fetch(stderr)
    )
end

# Detect previous installation method:
#   "internal" -> the MOSEK distro was downloaded and installed by the installer
#   "external" -> the MOSEK distro was detected in some other location and not downloaded by the installer
# If the previous installation was internal, we will not attempt to locate a MOSEK installation
#
# This ensures that you can rebuild an internal installation without
# the installer checking for other installations.
instmethod =
    try
        open(joinpath(bindepsdir,"inst_method"),"r") do f
            strip(read(f,String))
        end
    catch
        nothing
    end

curmosekbindir =
    try
        open(joinpath(bindepsdir,"mosekbindir"),"r") do f
            strip(read(f,String))
        end
    catch
        nothing
    end
forcedownload = haskey(ENV,"MOSEKJL_FORCE_DOWNLOAD")

mskbindir =
# 1. If MOSEKBINDIR we use that path (and if it is not valid we produce an error message)
    if  ! forcedownload && haskey(ENV,"MOSEKBINDIR")
        mosekbindir = ENV["MOSEKBINDIR"]

        if ! bindirIsCurrentVersion(mosekbindir)
            error("MOSEKBINDIR ($mosekbindir) does not point to a MOSEK bin directory")
        end
        instmethod = "external"

        mosekbindir
# 2a. If last build used a user-specified MOSEK, we check if that is still valid and use that again
    elseif ! forcedownload &&
        instmethod == "external" &&
        curmosekbindir != nothing &&
        bindirIsCurrentVersion(curmosekbindir)

        curmosekbindir
# 2b. Otherwise, look in the UNIX default installation path
    elseif ! forcedownload &&
        haskey(ENV,"HOME") &&
        bindirIsCurrentVersion(joinpath(ENV["HOME"],"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin"))

        instmethod = "external"

        joinpath(ENV["HOME"],"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin")
# 2c. Windows default home
    elseif ! forcedownload &&
        haskey(ENV,"HOMEDRIVE") &&
        haskey(ENV,"HOMEPATH") &&
        bindirIsCurrentVersion(joinpath(string(ENV["HOMEDRIVE"],ENV["HOMEPATH"]),"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin"))

        home = string(ENV["HOMEDRIVE"],ENV["HOMEPATH"])

        instmethod = "external"

        joinpath(home,"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin")
# 2d. Window global installation
    elseif ! forcedownload && hasregistrykey(WinReg.HKEY_LOCAL_MACHINE,"SOFTWARE\\mosek$mskvmajor$mskvminor","InstallDir")
        instmethod = "external"
        p = getregistrykey(WinReg.HKEY_LOCAL_MACHINE,"SOFTWARE\\MOSEK$mskvmajor$mskvminor","InstallDir")
	# Fix for broken  registry entry:
	p = replace(p,"\\\\" => "\\")
        joinpath(p,"tools","platform",mskplatform,"bin")
# 3. Otherwise, fetch the MOSEK distro and unpack it
    else
        srcdir   = joinpath(bindepsdir,"src")
        dldir    = joinpath(bindepsdir,"downloads")
        archname = "mosektools$(mskplatform)$(distroext)"
        hosturl  = "https://www.mosek.com/downloads/default_dns.txt"

        mkpath(dldir)

        downloadhost = "download.mosek.com"
        #dlcmd = download_cmd(hosturl, joinpath(dldir,"downloadhostname"))
        #@info("Download command: $dlcmd")
        #(res,stdout,stderr) = collect_output(dlcmd)
        #downloadhost =
        #    if res != 0
        #        @error(String(stderr))
        #        error("Failed to get MOSEK download host")
        #    else
        #        open(joinpath(dldir,"downloadhostname"),"r") do f
        #            strip(read(f,String))
        #        end
        #    end

        verurl   = "https://$downloadhost/stable/$mskvmajor.$mskvminor/version"

        cur_version =
            if isfile(joinpath(bindepsdir,"version"))
                open(joinpath(bindepsdir,"version"),"r") do f
                    cur_version = strip(read(f,String))
                end
            else
                nothing
            end

        @info("Get latest MOSEK version ($verurl)")
        success(download_cmd(verurl, joinpath(dldir,"new_version"))) || error("Failed to get MOSEK version")

        new_version =
            open(joinpath(dldir,"new_version"),"r") do f
                strip(read(f,String))
            end
        if cur_version === nothing
            @info("Latest MOSEK version = $new_version, nothing currently installed")
        else
            @info("Latest MOSEK version = $new_version, currently installed = $cur_version")
        end

        instmethod = "internal"

        if  cur_version == nothing ||
            cur_version != new_version ||
            !bindirIsCurrentVersion(joinpath(srcdir,"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin"))

            verarr = split(new_version,'.')

            archurl = "https://$downloadhost/stable/$(new_version)/$archname"
            @info("Download MOSEK distro ($archurl)")

            basename,ext,sndext = splittarpath(archname)

            success(download_cmd(archurl, joinpath(dldir,archname))) || error("Failed to download MOSEK distro")

            mkpath(srcdir)
            @info("Unpack MOSEK distro ($dldir/$archname -> $srcdir)")
            success(unpack_cmd(joinpath(dldir,archname),srcdir, ext, sndext)) || error("Failed to unpack MOSEK distro")

            @info("MOSEK installation complete.")
            joinpath(srcdir,"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin")
        else
            @info("Update not necessary")
            joinpath(srcdir,"mosek","$mskvmajor.$mskvminor","tools","platform",mskplatform,"bin")
        end
    end

mskbindir = replace(mskbindir,"\\" => "/")

version = versionFromBindir(mskbindir)
@info("mskbindir = $mskbindir")

if version == nothing
    error("MOSEK package is broken")
else
    open(joinpath(bindepsdir,"version"),"w") do f
        write(f,version)
    end
end

verarr = split(version,'.')

libmosekpath= findlibs(mskbindir,verarr[1],verarr[2])

if instmethod == "external"
    @info("""Found MOSEK $version at $mskbindir""")
end

open(joinpath(bindepsdir,"inst_method"),"w") do f
    write(f,instmethod)
end
open(joinpath(bindepsdir,"mosekbindir"),"w") do f
    write(f,mskbindir)
end

libmosekpath = escape_string(normpath(libmosekpath))

open(joinpath(bindepsdir,"deps.jl"),"w") do f
    write(f,"""
# This is an auto-generated file; do not edit
# Macro to load a library
import Libdl
macro checked_lib(libname, path)
    (Libdl.dlopen_e(path) == C_NULL) && error("Unable to load \\n\\n\$libname (\$path)\\n\\nPlease re-run Pkg.build(package), and restart Julia.")
    quote const \$(esc(libname)) = \$path end
end

# Load dependencies
@checked_lib libmosek      "$libmosekpath"
""")

end # open


text:
#handle case with 1 link index at edges.
function fix_autoMPO1(W::ITensor)::ITensor
    ils=filterinds(W,"Link")
    iss=filterinds(W,"Site")
    @assert length(ils)==1
    @assert length(iss)==2
    il=ils[1] #link index
    Dw=dim(il)
    p=collect(1:Dw)
    p[2],p[Dw]=p[Dw],p[2]
    W1=ITensor(il,iss...)
    for js in eachindval(iss)
        for jl in eachindval(il)
            Sl=slice(W,jl)
            assign!(W1,Sl,il=>p[jl.second])
        end
    end
    return W1
end

function fix_autoMPO2(W::ITensor)::ITensor
    ils=filterinds(W,"Link")
    if length(ils)==1
        return fix_autoMPO1(W)
    end
    iss=filterinds(W,"Site")
    @assert length(ils)==2
    @assert length(iss)==2
    d,n,r,c=parse_links(W)
    Dw1,Dw2=dim(r),dim(c)
    #
    #  set up perm arrays to swap row and col 2 with N
    #
    pr=collect(1:Dw1)
    pc=collect(1:Dw2)
    pr[2],pr[Dw1]=pr[Dw1],pr[2]
    pc[2],pc[Dw2]=pc[Dw2],pc[2]
    W1=ITensor(r,c,iss...)
    for jr in eachindval(r)
        for jc in eachindval(c)
            sl=slice(W,jr,jc)
            assign!(W1,sl,r=>pr[jr.second],c=>pc[jc.second])
        end
    end
    return W1
end

@doc """
    fix_autoMPO!(H::MPO)

Convert AutoMPO output into lower regular form.  If `Dw` x `Dw'` are the dimensions of
the Link indices for any site then this routine simply does 2 swaps:
1. Swap row 2 with row `Dw`
2. Swap column 2 with column `Dw'`

# Arguments
- `H::MPO` : MPO to be characterized.

MPO is modified in place into lower regular form

# Examples

```julia
julia> using ITensors
julia> using ITensorMPOCompression
julia> N=5
julia> sites = siteinds("S=1/2",5);
julia> ampo = OpSum();
julia> for j=1:N-1
julia>     add!(ampo,1.0, "Sz", j, "Sz", j+1);
julia>     add!(ampo,0.5, "S+", j, "S-", j+1);
julia>     add!(ampo,0.5, "S-", j, "S+", j+1);
julia> end
julia> H=MPO(ampo,sites);
julia> pprint(H[2])
I 0 0 0 0
0 I S S S
S 0 0 0 0
S 0 0 0 0
S 0 0 0 0

julia> fix_autoMPO!(H)
julia> pprint(H[2])
I 0 0 0 0
S 0 0 0 0
S 0 0 0 0
S 0 0 0 0
0 S S S I
```
"""
function fix_autoMPO!(H::MPO)
    N=length(H)
    for n in 1:N
        H[n]=fix_autoMPO2(H[n])
    end
end

make_Heisenberg_AutoMPO(sites,NNN::Int64,hz::Float64,ul::reg_form,J::Float64=1.0)::MPO =
    make_Heisenberg_AutoMPO(sites,NNN,hz,J)

@doc """
    make_Heisenberg_AutoMPO(sites,[NNN=1[,hz=0.0[,J=1.0]]])

Use `ITensor.autoMPO` to build up a Heisenberg model Hamiltonian with up to `NNN` neighbour
interactions.  The interactions are hard coded to decay like J/(i-j) between sites `i` and `j`.
The MPO is returned in lower regular form.

# Arguments
- `sites` : Site set defining the lattice of sites.
- `NNN::Int64` : Number of nearest neighbour interactions to include in `H`
- `hz::Float64=0.0` : External magnetic field in `z` direction.
- `J::Float64=1.0` : Nearest neighbour interaction strength.

"""
function make_Heisenberg_AutoMPO(sites,NNN::Int64=1,hz::Float64=0.0,J::Float64=1.0)::MPO
    N=length(sites)
    @assert(N>=NNN)
    ampo = OpSum()
    for j=1:N
        add!(ampo, hz   ,"Sz", j)
    end
    for dj=1:NNN
        f=J/dj
        for j=1:N-dj
            add!(ampo, f    ,"Sz", j, "Sz", j+dj)
            add!(ampo, f*0.5,"S+", j, "S-", j+dj)
            add!(ampo, f*0.5,"S-", j, "S+", j+dj)
        end
    end
    mpo=MPO(ampo,sites)
    fix_autoMPO!(mpo) #swap row[2]<->row[Dw] and col[2]<->col[Dw]
    return mpo
end

make_transIsing_AutoMPO(sites,NNN::Int64,hx::Float64,ul::reg_form,J::Float64=1.0)::MPO =
     make_transIsing_AutoMPO(sites,NNN,hx,J)

@doc """
    make\\_transIsing\\_AutoMPO(sites,NNN[,hx=0.0[,J=1.0]])

 Use `ITensor.autoMPO` to build up a transverse Ising model Hamiltonian with up to `NNN` neighbour
 interactions.  The interactions are hard coded to decay like J/(i-j) between sites `i` and `j`.
 The MPO is returned in lower regular form.

 # Arguments
 - `sites` : Site set defining the lattice of sites.
 - `NNN::Int64` : Number of nearest neighbour interactions to include in `H`
 - `hx::Float64=0.0` : External magnetic field in `x` direction.
 - `J::Float64=1.0` : Nearest neighbour interaction strength.

 """
function make_transIsing_AutoMPO(sites,NNN::Int64=1,hx::Float64=0.0,J::Float64=1.0)::MPO
    do_field = hx!=0.0
    N=length(sites)
    @assert(N>NNN)
    ampo = OpSum()
    if do_field
        for j=1:N
            add!(ampo, hx   ,"Sx", j)
        end
    end
    for dj=1:NNN
        f=J/dj
        for j=1:N-dj
            add!(ampo, f    ,"Sz", j, "Sz", j+dj)
        end
    end
    mpo=MPO(ampo,sites)
    fix_autoMPO!(mpo) #swap row[2]<->row[Dw] and col[2]<->col[Dw]
    return mpo
end

@doc """
    make_transIsing_MPO(sites[,NNN=1[,hx=0.0[,ul=lower[,J=1.0]]]])

Directly coded build up of a transverse Ising model Hamiltonian with up to `NNN` neighbour
interactions.  The interactions are hard coded to decay like J/(i-j) between sites `i` and `j`.

# Arguments
- `sites` : Site set defining the lattice of sites.
- `NNN::Int64=1` : Number of nearest neighbour interactions to include in `H`
- `hx::Float64=0.0` : External magnetic field in `x` direction.
- `ul::reg_form=lower` : build H with `lower` or `upper` regular form.
- `J::Float64=1.0` : Nearest neighbour interaction strength.

"""
function make_transIsing_MPO(sites,NNN::Int64=1,hx::Float64=0.0,ul::reg_form=lower,J::Float64=1.0;kwargs...)::MPO
    pbc::Bool=get(kwargs,:pbc,false)
    mpo=MPO(sites) #make and MPO only to get the indices
    prev_link=Nothing
    for n in 1:length(sites)
        mpo[n]=make_transIsing_op(mpo[n],prev_link,NNN,J,hx,ul)
        prev_link=filterinds(mpo[n],tags="Link,l=$n")[1]
    end
    if pbc
        return mpo #leave end Ws as full matrices, i.e. periodic boundary conditions.
    else
        return ITensorMPOCompression.to_openbc(mpo) #contract with l* and *r at the edges.
    end
end

#  It turns out that the trans-Ising model
#  Dw = 2+Sum(i,i=1..NNN) =2 + NNN*(NNN+1)/2
#
function transIsing_Dw(NNN::Int64)::Int64
    return 2+NNN*(NNN+1)/2
end

function make_Ising_index(Dw::Int64,tags::String,use_qn::Bool,dir)
    if (use_qn)
        ind=Index(QN("Sz",0)=>Dw;dir=dir,tags=tags)
    else
        ind=Index(Dw,tags)
    end
    return ind
end

# NNN = Number of Nearest Neighbours, for example
#    NNN=1 corresponds to nearest neighbour
#    NNN=2 corresponds to nearest and next nearest neighbour
function make_transIsing_op(Wref::ITensor,prev_link,NNN::Int64,J::Float64,hx::Float64=0.0,ul::reg_form=lower)::ITensor
    @assert NNN>=1
    do_field = hx!=0.0
    Dw::Int64=transIsing_Dw(NNN)
    is=filterinds(Wref,tags="Site")[1] #get any site index for generating operators
    use_qn=hasqns(is)
    d,n,r,c=parse_links(Wref)
    if tags(r)==TagSet("")
        r=make_Ising_index(Dw,"Link,l=$(n-1)",use_qn,ITensors.In)
    else
        r=prev_link
    end
    if tags(c)==TagSet("")
        c=make_Ising_index(Dw,"Link,l=$n",use_qn,ITensors.Out)
    else
        c=redim(c,Dw)
    end
    iblock=1;

    W=ITensor(r,dag(c),is,dag(is'))
    Id=op(is,"Id")
    Sz=op(is,"Sz")
    if do_field
        Sx=op(is,"Sx")
    end
    assign!(W,Id,r=>1 ,c=>1 )
    assign!(W,Id,r=>Dw,c=>Dw)

    if ul==lower
        if do_field
            assign!(W ,hx*Sx,r=>Dw,c=>1); #add field term
        end
        #very hard to explain this without a diagram.
        for iNN in 1:NNN
            assign!(W,Sz,r=>iblock+1,c=>1)
            for jNN in 1:iNN-1
                assign!(W,Id,r=>iblock+1+jNN,c=>iblock+jNN)
            end
            Jn=J/(iNN) #interactions need to decay with distance in order for H to extensive
            assign!(W,Jn*Sz,r=>Dw,c=>iblock+iNN)
            iblock+=iNN
        end
    else
        if do_field
            assign!(W,hx*Sx,r=>1,c=>Dw ); #add field term
        end
        #very hard to explain this without a diagram.
        for iNN in 1:NNN
            assign!(W,Sz,r=>1,c=>iblock+1)
            for jNN in 1:iNN-1
                assign!(W,Id,r=>iblock+jNN,c=>iblock+1+jNN)
            end
            Jn=J/(iNN) #interactions need to decay with distance in order for H to extensive
            assign!(W,Jn*Sz,r=>iblock+iNN,c=>Dw)
            iblock+=iNN
        end
    end
    return W
end

function to_openbc(mpo::MPO)::MPO
    N=length(mpo)
    l,r=get_lr(mpo)
    mpo[1]=l*mpo[1]
    mpo[N]=mpo[N]*r
    @assert length(filterinds(inds(mpo[1]),tags="Link"))==1
    @assert length(filterinds(inds(mpo[N]),tags="Link"))==1
    return mpo
end

function get_lr(mpo::MPO)::Tuple{ITensor,ITensor}
    ul::reg_form = is_lower_regular_form(mpo,1e-14) ? lower : upper

    N=length(mpo)
    W1=mpo[1]
    llink=filterinds(inds(W1),tags="l=0")[1]
    l=ITensor(0.0,dag(llink))

    WN=mpo[N]
    rlink=filterinds(inds(WN),tags="l=$N")[1]
    r=ITensor(0.0,dag(rlink))
    if ul==lower
        l[llink=>dim(llink)]=1.0
        r[rlink=>1]=1.0
    else
        l[llink=>1]=1.0
        r[rlink=>dim(rlink)]=1.0
    end

    return l,r
end

function fast_GS(H::MPO,sites)::Tuple{Float64,MPS}
    psi0  = randomMPS(sites,length(H))
    sweeps = Sweeps(5)
    setmaxdim!(sweeps, 2,4,8,16,32)
    setcutoff!(sweeps, 1E-10)
    E,psi= dmrg(H,psi0, sweeps;outputlevel=0)
    return E,psi
end

ext: jl


text:
function calculate_velocities(
    m::AbstractMetricParams{T},
    init_pos,
    α_generator,
    β::T,
) where {T}
    [map_impact_parameters(m, init_pos, α, β) for α in α_generator]
end

function calculate_velocities(
    m::AbstractMetricParams{T},
    init_pos,
    α_genetator,
    β_generator,
) where {T}
    [map_impact_parameters(m, init_pos, α, β) for α in α_genetator, β in β_generator]
end

"""
In-place specialisation
"""
function calculate_velocities!(
    vs,
    m::AbstractMetricParams{T},
    init_pos,
    α_generator,
    β::T,
) where {T}
    for (i, α) in enumerate(α_generator)
        vs[i] = map_impact_parameters(m, init_pos, α, β)
    end
end

# have to use a slight 0.001 offset to avoid integrating α=0.0 geodesics
x_to_α(X, x_mid, fov_factor) = (X + 0.001 - x_mid) / fov_factor
y_to_β(Y, y_mid, fov_factor) = (Y - y_mid) / fov_factor


text:
type Ledger
    accounts::Set{Account}
    entries::SortedMultiDict{Date,Entry}
    lastid::Int64
end

Ledger() = Ledger(
    Set{Account}(),
    SortedMultiDict(Date[], Entry[]),
    0)

nextid!(l::Ledger) = (l.lastid += 1)
function push!(l::Ledger, ent::Entry)
    id!(ent, nextid!(l))
    insert!(l.entries, date(ent), ent)
    l
end
function register!(l::Ledger, acc::Account)
    id!(acc, nextid!(l))
    push!(l.accounts, acc)
    acc
end

# Make account type functions.
for acctype in [:Asset, :Liability, :Equity, :Revenue, :Expense, :Trading]
    # Get symbols.
    lc = lowercase(string(acctype))
    asym = Expr(:quote, Symbol(lc))  # :asset
    mutfn = Symbol(lc, "!")  # :asset!

    # Make constructor (Asset).
    @eval ($acctype)(name) = Account(-1, $asym, Compat.UTF8String(name))

    # Define asset! method.
    @eval begin
        ($mutfn)(l::Ledger, name::AbstractString) =
            register!(l, ($acctype)(name))
    end
end

function transfer!(
    sys::Ledger,
    date::Date,
    credit::Account,
    debit::Account,
    amount::Monetary,
    description::AbstractString = "")
    ent = Entry(date, Compat.UTF8String(description), [
        Split(debit, amount),
        Split(credit, -amount)
    ])
    push!(sys, ent)
    ent
end

# access functions
function debit(ledger::Ledger, account::Account, cur::Symbol)
    netdebit = zero(Monetary{cur})
    for (date, tx) in ledger.entries
        for split in tx.splits
            if split.account ≡ account && currency(split) == cur
                netdebit += debit(split)
            end
        end
    end
    netdebit
end

immutable Balance
    account::Account
    netdebit::Monetary
end
currency(b::Balance) = currency(b.netdebit)
debit(b::Balance) = b.amount

function balances(ledger::Ledger)
    ret = Dict{Account, Basket}()
    for (date, tx) in ledger.entries
        for split in tx.splits
            symb = currency(split)
            key = split.account
            if haskey(ret, key)
                ret[key][symb] += debit(split)
            else
                ret[key] = Basket(debit(split))
            end
        end
    end
    ret
end


text:
using JuMP
using MathOptInterface

const MOI = MathOptInterface

struct Solution
    lsol::Array{Float64,2}
    πsol::Vector{Float64}
end

struct Result
    status::MOI.TerminationStatusCode
    sol::Solution
    value::Float64
end

struct LP <: PipeDimensioningSolver
    lpsolver
end

function make_model(inst::Instance, topo::Topology, optimizer)
    nnodes = length(inst.nodes)
    length(topo.nodes) == nnodes ||
        throw(ArgumentError("Steiner nodes not allowed"))
    arcs, narcs = topo.arcs, length(topo.arcs)
    ndiams = length(inst.diameters)

    q = uniq_flow(inst, topo)

    model = JuMP.Model(optimizer)

    # squared pressure variables at nodes
    lb = [b.lb^2 for b in inst.pressure]
    ub = [b.ub^2 for b in inst.pressure]
    @variable(model, lb[v] <= π[v=1:nnodes] <= ub[v])

    # relative length of pipe segments with specific diameter
    @variable(model, 0.0 <= l[1:narcs, 1:ndiams] <= 1.0)

    # convex combination of pipe segments
    @constraint(model, convexcomb[a=1:narcs], sum(l[a,i] for i=1:ndiams) == 1.0)

    # pressure loss constraint, with variable factor from diameter choice
    L = pipelengths(topo)
    D5 = [diam.value^5 for diam in inst.diameters]
    C = inst.ploss_coeff * L .* q .* abs.(q)
    @constraint(model, ploss[a=1:narcs], π[arcs[a].tail] - π[arcs[a].head]
                == sum(C[a] / D5[i] * l[a,i] for i=1:ndiams))

    # minimize total construction cost
    cost = [diam.cost for diam in inst.diameters]
    @objective(model, Min, sum(cost[i] * L[a] * l[a,i] for a=1:narcs for i=1:ndiams))

    model, π, l
end

function PipeLayout.optimize(inst::Instance, topo::Topology, solver::LP)
    model, π, l = make_model(inst, topo, solver.lpsolver)
    JuMP.optimize!(model)
    status = JuMP.termination_status(model)
    objval = status == MOI.OPTIMAL ? JuMP.objective_value(model) : Inf
    sol = Solution(JuMP.value.(l), JuMP.value.(π))
    Result(status, sol, objval)
end


text:
@doc raw"""
    LogEuclideanMetric <: RiemannianMetric

The LogEuclidean Metric consists of the Euclidean metric applied to all elements after mapping them
into the Lie Algebra, i.e. performing a matrix logarithm beforehand.
"""
struct LogEuclideanMetric <: RiemannianMetric end

@doc raw"""
    distance(M::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric}, p, q)

Compute the distance on the [`SymmetricPositiveDefinite`](@ref) manifold between
`p` and `q` as a [`MetricManifold`](@ref) with [`LogEuclideanMetric`](@ref).
The formula reads

```math
    d_{\mathcal P(n)}(p,q) = \lVert \operatorname{Log} p - \operatorname{Log} q \rVert_{\mathrm{F}}
```

where $\operatorname{Log}$ denotes the matrix logarithm and
$\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
"""
function distance(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogEuclideanMetric},
    p,
    q,
) where {N}
    return norm(log(Symmetric(p)) - log(Symmetric(q)))
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Particle Cloud Simulation Example
#
#

#
#
#

#
# Here we simulate the trajectory of a particle cloud in a randomly generated flow field which is a crude representation for an ocean eddy field.
#
# As an exercise left to the user, directions are provided e.g. to modify the size of the domain or the initial position of the particle cloud.
#
# The demonstrated methods that are generally applicable for research purposes (e.g. with flow field derived from models or data). The notebook also demoes simple methods to create plots and movies.
#
# - 1. setup the software and initialize example
# - 2. simulate trajectories & plot results
# - 3. exercises (user)

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### For More Documentation
#
# On Julia :
# - <https://julialang.org>
# - <https://docs.julialang.org>
#
# On this notebook :
# - <https://docs.juliadiffeq.org/latest>
# - <https://en.wikipedia.org/wiki/Displacement_(vector)>
# - <https://juliaclimate.github.io/IndividualDisplacements.jl/dev>
# - <https://juliaclimate.github.io/MeshArrays.jl/dev>

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 1.1 Import Software
# -

using OrdinaryDiffEq, Plots, ColorSchemes
using IndividualDisplacements, MeshArrays
using Statistics, Random, DataFrames
p=dirname(pathof(MeshArrays)); include(joinpath(p,"../examples/Demos.jl"))
include("helper_functions.jl")

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ## 1.2 Setup Problem
# -

# Put grid variables in a dictionary.

np=24
Γ=SetupPeriodicDomain(np);

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# Derive flow field (`𝑃["u0"]` etc) from randomly generated ϕ streamfunction
# -

𝑃,ϕ=SetupRandomFlow(Γ);
#𝑃["t1"]=100.0

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 1.3 Initial Positions
#
# Set the initial positions to a region of the gridded domain. These will provide the initial condition `u0` for the differential equation solver.
#
# _Skipped at first in presentation mode_

# + {"slideshow": {"slide_type": "skip"}}
x0=np*(0.25:0.02:0.75)
y0=np*(0.25:0.02:0.75)
#x0=np*(0.2:0.005:0.25)
#y0=np*(0.6:0.005:0.65)

#reformat data into u0 initial condition
x0=vec(x0)*ones(1,length(y0))
y0=ones(size(x0,1),1)*transpose(vec(y0))
u0=transpose([x0[:] y0[:] ones(size(x0[:]))]);

# + {"slideshow": {"slide_type": "subslide"}}
x=vec(Γ["XC"][1][:,1])
y=vec(Γ["YC"][1][1,:])
χ=transpose(ϕ[1])

plt=contourf(x,y,χ,linewidth = 0.1,clims=(-0.4,0.4),
             title="Flow Stream (color) & Initial Positions (dots)")
scatter!(plt,u0[1,:],u0[2,:],c=:white,markersize=3.0,leg=false)

#pth=tempdir()*"/"
#savefig(plt,pth*"RandomFlow.png")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 2.1 Compute Trajectories
#
# As done in the `SolidBodyRotation` example earlier, we simply use the `OrdinaryDiffEq` package to solve for the particle trajectories.
# -

𝑇 = (𝑃["t0"],𝑃["t1"])
prob = ODEProblem(⬡!,u0,𝑇,𝑃)
sol = solve(prob,Tsit5(),reltol=1e-5,abstol=1e-5)
size(sol)

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ## 2.2 Post-Process Output
#
# Reformat the `sol` output from the `solve` function into a `DataFrame` for practical use, plotting, etc.
#
# _Code cell is skipped at first in presentation mode_

# + {"slideshow": {"slide_type": "skip"}}
#x,y axes etc
x=sol[1,:,:]
y=sol[2,:,:]
fIndex=sol[3,:,:]
ID=collect(1:size(sol,2))*ones(1,size(sol,3))

#put together as a DataFrame
df = DataFrame(ID=Int.(ID[:]), x=mod.(x[:],Ref(np)), y=mod.(y[:],Ref(np)), fIndex=fIndex[:]);

# + {"slideshow": {"slide_type": "skip"}}
#add a time axis
nf=size(u0,2)
nt=size(df,1)/nf
t=[ceil(i/nf)-1 for i in 1:nt*nf]

df[!,:t]=(𝑃["t1"]-𝑃["t0"])/t[end].*t;

#df

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 2.3 Plot Results
#
# - define `scatter_subset` convenience function
# - generate animation using `scatter_subset`

# + {"slideshow": {"slide_type": "subslide"}}
function scatter_subset(df,t)
    dt=5.0
    df_t = df[ (df.t.>t-dt).&(df.t.<=t) , :]

    x=vec(Γ["XC"][1][:,1])
    y=vec(Γ["YC"][1][1,:])
    χ=transpose(ϕ[1])
    contourf(x,y,χ,c = :blues,linewidth = 0.1,clims=(-0.4,0.4),colorbar=false)

    scatter!(df_t.x,df_t.y,markersize=2.0,c=:red,
    xlims=(0,np),ylims=(0,np),leg=:none,marker = (:circle, stroke(0)))
end
# + {"slideshow": {"slide_type": "subslide"}}
#@gif for t in 0:2.0:maximum(df[!,:t])
#   scatter_subset(df,t)
#end

anim = @animate for t in 0:2.0:maximum(df[!,:t])
   scatter_subset(df,t)
end
pth=tempdir()*"/"
gif(anim, pth*"RandomFlow.gif", fps = 15)
# -

# ## 3. Exercises
#
#
# ### Try uncomment some of these lines:
#
# - `#df` (see the [package docs](https://juliadata.github.io/DataFrames.jl/stable/))
# - the three lines below at once, or one at a time
#     - `#x0=np*(0.2:0.005:0.25)`
#     - `#y0=np*(0.6:0.005:0.65)`
#     - `𝑃["t1"]=100`
#
# _Tip: select `run all below` from the `Cell` drop down menu to avoid replacing the random flow field every time._
#
# ### Try the following exercise (after uncommenting the above lines):
#
# Can you find a choice of `x0` & `y0` that reduces convergence in one direction? Or increases convergence in one direction?
#
# - more `convergence` ~ less relative `separation` ~ reduced xx,yy below
# - less `convergence` ~ more relative `separation` ~ increased xx,yy below
#
# _See code & plot below_

# + {"slideshow": {"slide_type": "skip"}}
#Same DataFrame as eefore but with original x,y
df2=deepcopy(df)
df2.x[:]=sol[1,:,:]
df2.y[:]=sol[2,:,:];

# +
#Compute relative separation estimates
tt=0:2.0:maximum(df[!,:t])
xx=zeros(length(tt)); yy=similar(xx)
dt=5.0
for t in 1:length(tt)
    df_t = df2[ (df2.t.>tt[t]-dt).&(df2.t.<=tt[t]) , :]
    xx[t]=std(df_t.x)
    yy[t]=std(df_t.y)
end

#Plot relative separation estimates
plot(tt,xx,label="std(x)",title="relative separation estimates")
plot!(tt,yy,label="std(y)")
# -

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
#========== making ForwardDiff work with LoopVectorization ==========#

using .LoopVectorization
using Core: VecElement

using .ForwardDiff
using .ForwardDiff: Dual, Partials, partials

#=
# using Tullio.LoopVectorization: LoopVectorization, SVec, vconvert, SIMDPirates

s1 = SVec{2,Float64}(5.5, 6.6) # SVec{2,Float64}<5.5, 6.6>
# dump(s1)
# SVec{2,Float64}
#   data: Tuple{VecElement{Float64},VecElement{Float64}}
#     1: VecElement{Float64}
#       value: Float64 5.5
#     2: VecElement{Float64}
#       value: Float64 6.6
s1[2]
s1 |> typeof |> parentmodule # VectorizationBase

# @inline svec(tup::NTuple{N,T}) where {N,T} = SVec{N,T}(tup...)

d1 = Dual(1.23, (4,0,0))
typeof(d1) # Dual{Nothing,Float64,3}
# dump(d1)
# Dual{Nothing,Float64,2}
#   value: Float64 1.23
#   partials: Partials{2,Float64}
#     values: Tuple{Float64,Float64}
#       1: Float64 4.0
#       2: Float64 0.0
#       3: Float64 0.0
d1.partials # Partials{3,Float64}
d1.partials[1]

partials(d1, 1)
# @inline val(d::Dual) = d.value

=#

ForwardDiff.can_dual(::Type{<:SVec}) = true

@inline function Base.:+(x::Dual{Z,T,D}, sv::SVec{N,S}) where {Z,T<:Number,D,N,S}
    y = x.value + sv
    ps = ntuple(d -> x.partials.values[d] + zero(sv), Val(D))
    TS = SVec{N,promote_type(T,S)}
    Dual{Z,TS,D}(y, Partials{D,TS}(ps))
end
@inline function Base.:+(sv::SVec{N,S}, x::Dual{Z,T,D}) where {Z,T<:Number,D,N,S}
    y = x.value + sv
    ps = ntuple(d -> x.partials.values[d] + zero(sv), Val(D))
    TS = SVec{N,promote_type(T,S)}
    Dual{Z,TS,D}(y, Partials{D,TS}(ps))
end

@inline function Base.:*(x::Dual{Z,SVec{N,T},D}, sv::SVec{N,S}) where {Z,T,D,N,S}
    y = x.value * sv
    ps = ntuple(d -> x.partials.values[d] * sv, Val(D))
    TS = SVec{N,promote_type(T,S)}
    Dual{Z,typeof(y),D}(y, Partials{D,typeof(y)}(ps))
end
@inline function Base.:*(sv::SVec{N,S}, x::Dual{Z,SVec{N,T},D}) where {Z,T,D,N,S}
    y = sv * x.value
    ps = ntuple(d -> sv * x.partials.values[d], Val(D))
    TS = SVec{N,promote_type(T,S)}
    Dual{Z,TS,D}(y, Partials{D,TS}(ps))
end

@inline function Base.:*(p::Partials{D,SVec{N,T}}, sv::SVec{N,S}) where {T,D,N,S}
    TS = SVec{N,promote_type(T,S)}
    Partials{D,TS}(ntuple(d -> p.values[d] * sv, Val(D)))
end
@inline function Base.:*(sv::SVec{N,S}, p::Partials{D,SVec{N,T}}) where {T,D,N,S}
    TS = SVec{N,promote_type(T,S)}
    Partials{D,TS}(ntuple(d -> sv * p.values[d], Val(D)))
end

#========== the end ==========#


text:
module ReverseDiffX

export NotTracked

using MacroTools, LinearAlgebra, ..ReverseDiff, StaticArrays
using Base.Broadcast: BroadcastStyle, ArrayStyle, Broadcasted, broadcasted
using ForwardDiff: ForwardDiff, Dual
using ..ReverseDiff: SpecialInstruction, value, value!, deriv, track, record!, tape, unseed!, @grad, TrackedReal, TrackedVector, TrackedMatrix, TrackedArray
using ..DistributionsAD: DistributionsAD

const TrackedVecOrMat{V,D} = Union{TrackedVector{V,D},TrackedMatrix{V,D}}

import SpecialFunctions, NaNMath
import ..DistributionsAD: turing_chol, symm_turing_chol, _mv_categorical_logpdf
import Base.Broadcast: materialize
import StatsFuns: logsumexp
import ZygoteRules

const RDBroadcasted{F, T} = Broadcasted{<:Any, <:Any, F, T}

using Distributions, PDMats
import Distributions: logpdf,
                      Gamma,
                      MvNormal,
                      MvLogNormal,
                      Dirichlet,
                      Wishart,
                      InverseWishart,
                      PoissonBinomial,
                      isprobvec

using ..DistributionsAD: TuringMvNormal,
                         TuringMvLogNormal,
                         TuringWishart,
                         TuringInverseWishart,
                         TuringDirichlet,
                         TuringPoissonBinomial,
                         TuringScalMvNormal,
                         TuringDiagMvNormal,
                         TuringDenseMvNormal

include("reversediffx.jl")

function PoissonBinomial(p::TrackedArray{<:Real}; check_args=true)
    return TuringPoissonBinomial(p; check_args = check_args)
end

Gamma(α::TrackedReal, θ::Real; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::Real, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::TrackedReal, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
pgamma(α, θ; check_args=true) = Gamma(promote(α, θ)..., check_args = check_args)
Gamma(α::T; check_args=true) where {T <: TrackedReal} = Gamma(α, one(T), check_args = check_args)
function Gamma(α::T, θ::T; check_args=true) where {T <: TrackedReal}
    check_args && Distributions.@check_args(Gamma, α > zero(α) && θ > zero(θ))
    return Gamma{T}(α, θ)
end

# Work around to stop TrackedReal of Inf and -Inf from producing NaN in the derivative
function Base.minimum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * minimum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(minimum)(d.ρ))
    end
end
function Base.maximum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * maximum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(maximum)(d.ρ))
    end
end

for T in (:TrackedVector, :TrackedMatrix)
    @eval begin
        function logpdf(d::MvNormal{<:Any, <:PDMats.ScalMat}, x::$T)
            logpdf(TuringScalMvNormal(d.μ, d.Σ.value), x)
        end
        function logpdf(d::MvNormal{<:Any, <:PDMats.PDiagMat}, x::$T)
            logpdf(TuringDiagMvNormal(d.μ, d.Σ.diag), x)
        end
        function logpdf(d::MvNormal{<:Any, <:PDMats.PDMat}, x::$T)
            logpdf(TuringDenseMvNormal(d.μ, d.Σ.chol), x)
        end

        function logpdf(d::MvLogNormal{<:Any, <:PDMats.ScalMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringScalMvNormal(d.normal.μ, d.normal.Σ.value)), x)
        end
        function logpdf(d::MvLogNormal{<:Any, <:PDMats.PDiagMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringDiagMvNormal(d.normal.μ, d.normal.Σ.diag)), x)
        end
        function logpdf(d::MvLogNormal{<:Any, <:PDMats.PDMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringDenseMvNormal(d.normal.μ, d.normal.Σ.chol)), x)
        end
    end
end

# zero mean, dense covariance
MvNormal(A::TrackedMatrix) = TuringMvNormal(A)

# zero mean, diagonal covariance
MvNormal(σ::TrackedVector) = TuringMvNormal(σ)

# dense mean, dense covariance
MvNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)

# dense mean, diagonal covariance
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvNormal(m, D)
end

# dense mean, diagonal covariance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)

# dense mean, constant variance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Real) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)

# dense mean, constant variance
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvNormal(m, A)
end

# zero mean,, constant variance
MvNormal(d::Int, σ::TrackedReal) = TuringMvNormal(d, σ)

# zero mean, dense covariance
MvLogNormal(A::TrackedMatrix) = TuringMvLogNormal(TuringMvNormal(A))

# zero mean, diagonal covariance
MvLogNormal(σ::TrackedVector) = TuringMvLogNormal(TuringMvNormal(σ))

# dense mean, dense covariance
MvLogNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))

# dense mean, diagonal covariance
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end

# dense mean, diagonal covariance
MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::TrackedVector{<:Real}, σ::Real)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end

# zero mean,, constant variance
MvLogNormal(d::Int, σ::TrackedReal) = TuringMvLogNormal(TuringMvNormal(d, σ))

Dirichlet(alpha::TrackedVector) = TuringDirichlet(alpha)
Dirichlet(d::Integer, alpha::TrackedReal) = TuringDirichlet(d, alpha)
for func_header in [
    :(simplex_logpdf(alpha::TrackedVector, lmnB::Real, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::TrackedVector)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::TrackedVector)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::Real, x::TrackedVector)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::TrackedReal, x::TrackedVector)),

    :(simplex_logpdf(alpha::TrackedVector, lmnB::Real, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::TrackedMatrix)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::TrackedMatrix)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::Real, x::TrackedMatrix)),
    :(simplex_logpdf(alpha::TrackedVector, lmnB::TrackedReal, x::TrackedMatrix)),
]
    @eval $func_header = track(simplex_logpdf, alpha, lmnB, x)
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractVector)
    simplex_logpdf(value(alpha), value(lmnB), value(x)), Δ -> begin
        (Δ .* log.(value(x)), -Δ, Δ .* (value(alpha) .- 1))
    end
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractMatrix)
    simplex_logpdf(value(alpha), value(lmnB), value(x)), Δ -> begin
        (log.(value(x)) * Δ, -sum(Δ), repeat(value(alpha) .- 1, 1, size(x, 2)) * Diagonal(Δ))
    end
end

function logpdf(d::MatrixBeta, X::AbstractArray{<:TrackedMatrix{<:Real}})
    return map(x -> logpdf(d, x), X)
end

Wishart(df::TrackedReal, S::Matrix{<:Real}) = TuringWishart(df, S)
Wishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringWishart(df, S)
Wishart(df::Real, S::TrackedMatrix) = TuringWishart(df, S)
Wishart(df::TrackedReal, S::TrackedMatrix) = TuringWishart(df, S)

InverseWishart(df::TrackedReal, S::Matrix{<:Real}) = TuringInverseWishart(df, S)
InverseWishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringInverseWishart(df, S)
InverseWishart(df::Real, S::TrackedMatrix) = TuringInverseWishart(df, S)
InverseWishart(df::TrackedReal, S::TrackedMatrix) = TuringInverseWishart(df, S)

function logpdf(d::Wishart, X::TrackedMatrix)
    return logpdf(TuringWishart(d), X)
end
function logpdf(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringWishart(d), X)
end

function logpdf(d::InverseWishart, X::TrackedMatrix)
    return logpdf(TuringInverseWishart(d), X)
end
function logpdf(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringInverseWishart(d), X)
end

# isprobvec

function isprobvec(p::TrackedArray{<:Real})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end
function isprobvec(p::SubArray{<:TrackedReal, 1, <:TrackedArray{<:Real}})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end

_mv_categorical_logpdf(ps::TrackedMatrix, x) = track(_mv_categorical_logpdf, ps, x)
@grad function _mv_categorical_logpdf(ps, x)
    ps_data = value(ps)
    probs = view(ps_data, x, :)
    ps_grad = zero(ps_data)
    sum(log, probs), Δ -> begin
        ps_grad .= 0
        ps_grad[x,:] .= Δ ./ probs
        return (ps_grad, nothing)
    end
end

end


text:
"""
    atmosphere_isa(height)
Calculate temperature, pressure, density and sound velocity for the given geopotential
height according to International Standard Atmosphere 1976.

# References
1.U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976

| Layer | h (m) | p (Pa)  | T (K)  | ``α`` (K/m) |
| ----- | ----- | ------- | ------ | ----------- |
| 0     | 0     | 101325  | 288.15 | -0.0065     |
| 1     | 11000 | 22632.1 | 216.65 | 0           |
| 2     | 20000 | 5474.89 | 216.65 | 0.001       |
| 3     | 32000 | 868.019 | 228.65 | 0.0028      |
| 4     | 47000 | 110.906 | 270.65 | 0           |
| 5     | 51000 | 66.9389 | 270.65 | -0.0028     |
| 6     | 71000 | 3.95642 | 214.65 | -0.002      |

Source: [https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere](https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere)

"""
function atmosphere_isa(height)

    if 0 <= height < 11000  # Troposphere
        α = -0.0065  # K/m
        T0 = 288.15      # K
        p0 = 101325.0    # Pa

        T = T0 + α * height
        p = p0 * (T0 / T) ^ (gD / (R_AIR * α))

    elseif 11000 <= height < 20000  # Tropopause
        T = 216.65    # K
        p0 = 22632.1  # Pa
        h0 = 11000    # m

        p = p0 * exp(-gD * (height - h0) / (R_AIR * T))

    elseif 20000 <= height < 32000  # Stratosphere 1
        α = 0.001    # K/m
        T0 = 216.65      # K
        p0 = 5474.89     # Pa
        h0 = 20000       # m

        T = T0 + α * (height - h0)
        p = p0 * (T0 / T) ^ (gD / (R_AIR * α))

    elseif 32000 <= height < 47000  # Stratosphere 2
        α = 0.0028   # K/m
        T0 = 228.65      # K
        p0 = 868.019     # Pa
        h0 = 32000       # m

        T = T0 + α * (height - h0)
        p = p0 * (T0 / T) ^ (gD / (R_AIR * α))

    elseif 47000 <= height < 51000  # Stratopause
        T = 270.65    # K
        p0 = 110.906  # Pa
        h0 = 47000    # m
        h0 = 47000    # m

        p = p0 * exp(-gD * (height - h0) / (R_AIR * T))

    elseif 51000 <= height < 71000  # Mesosphere 1
        α = -0.0028  # K/m
        T0 = 270.65      # K
        p0 = 66.9389     # Pa
        h0 = 51000       # m

        T = T0 + α * (height - h0)
        p = p0 * (T0 / T) ^ (gD / (R_AIR * α))

    elseif 71000 <= height <= 84500  # Mesosphere 2
        α = -0.002   # K/m
        T0 = 214.65      # K
        p0 = 3.95642     # Pa
        h0 = 71000       # m

        T = T0 + α * (height - h0)
        p = p0 * (T0 / T) ^ (gD / (R_AIR * α))

    else
        throw(DomainError(height, "height must be between 0 m and 84500 m"))
    end

    rho = p / (R_AIR * T)
    a   = sqrt(γ_AIR * R_AIR * T)

    return [T, p, rho, a]
end


text:
using PropertyUtils
using Test

#-----------------------------------------------------------------------------# setup
mutable struct A
    a
end
Base.getproperty(a::A, x::Symbol) = "hello"

struct B
    b
end

#-----------------------------------------------------------------------------# tests
@testset "PropertyUtils" begin

@testset "@with" begin
    result = @with (x=1, y=2) x + y
    @test result == 3
end

@testset "fields" begin
    a = A(1)
    @test a.a == "hello"
    @test fields(a).a == 1
    fields(a).a = 2
    @test fields(a).a == 2
end

@testset "indexes" begin
    d = Dict(:x => 1, :y => 2)
    @test indexes(d).x == 1
    indexes(d).z = 3
    @test d[:z] == 3
end

@testset "joinprops" begin
    j = joinprops(fields(A(1)), B(2))
    result = @with j a + b
    @test result == 3
    j.a = 2
    result2 = @with j a + b
    @test result2 == 4
end

end


text:
using StatsBase

@testset "presmoothing" begin
    # df = 1
    smtab = presmoothing(sgtabX, Equate.fml₁)
    m = smtab.table.scale'smtab.table.prob
    @test m ≈ mean(smtab.raw)
    # df = 2
    smtab = presmoothing(sgtabX, Equate.fml₂)
    m = smtab.table.scale'smtab.table.prob
    s = ((smtab.table.scale .- m).^2)'smtab.table.prob
    @test s ≈ var(smtab.raw; corrected = false)
    # df = 3
    smtab = presmoothing(sgtabX, Equate.fml₃)
    m = smtab.table.scale'smtab.table.prob
    s = ((smtab.table.scale .- m).^2)'smtab.table.prob
    sk = ((smtab.table.scale .- m).^3)'smtab.table.prob / sqrt(s)^3
    # @test sk ≈ skewness(smtab.raw)
    @test round(sk; digits = 5) ≈ round(skewness(smtab.raw); digits = 5)
    # df = 4
    smtab = presmoothing(sgtabX, Equate.fml₄)
    m = smtab.table.scale'smtab.table.prob
    s = ((smtab.table.scale .- m).^2)'smtab.table.prob
    ku = ((smtab.table.scale .- m).^4)'smtab.table.prob / s^2 -3
    @test ku ≈ kurtosis(smtab.raw)
end

@testset "KernelSmoothing (Gaussian)" begin
    fml = LogLinearFormula(6)
    smtab = presmoothing(sgtabX, fml)
    bbw = EstBandwidth(smtab)
    @test bbw.ls_success
    ktab = KernelSmoothing(smtab; hX = exp(bbw.minimizer[1]))
    @test typeof(ktab) == Equate.KernelFreqTab
end

@testset "Examine smoothing" begin
    fit = presmoothing(sgtabX, 6)
    @test fit isa DataFrame
end


text:
######################################################################
# Packet Construction & Destructuring methods
# -----
# Licensed under Apache License 2.0
module PacketFactories
import ExtraFun: Ident
using ..Pia
using ..Packets
import ..Packets: packetid, classid

# ===== Factory Methods =====
export packetbytes
"""`packetbytes(classid::PacketClass, packetid, args...; kwargs...)`
Construct packet of class `classid` and type `packetid`, passing individually required `args` and `kwargs` as necessary.
While this is a factory method, the underlying construction methods may require specific arguments."""
packetbytes(packetid::Enum, args...; kwargs...) = packetbytes(classid(packetid), packetid, args...; kwargs...)
function packetbytes(classid::PacketClass, packetid::Enum, args...; kwargs...)
    buffer = IOBuffer()
    write(buffer, UInt8(classid))
    write(buffer, UInt8(packetid))
    packetbytes(Ident{classid}(), Ident{packetid}(), buffer, args...; kwargs...)
    return take!(buffer)
end

export packetcontent
"""`packetcontent(bytes)`
Destructure packet contents from the supplied packet data. Throws respective errors if not possible."""
function packetcontent(bytes; without_id::Bool = false)
    ids, bytes = bytes[1:2], bytes[3:length(bytes)]

    cls = PacketClass(ids[1])
    pkt = packetid(Ident{cls}(), ids[2])
    ret = packetcontent(Ident{cls}(), Ident{pkt}(), bytes)
    if without_id
        return ret
    else
        return (cls, pkt), ret
    end
end

export peekpacket
"""`peekpacket(bytes)`
Peek the packet class & packet ID from the packet header."""
function peekpacket(bytes)
    cls = PacketClass(bytes[1])
    pkt = packetid(cls, bytes[2])
    return cls, pkt
end

# ===== Connectivity Class Packets =====
# --- Handshake Start ---
function packetbytes(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeStart}, buffer::IOBuffer, appname::AbstractString, version::VersionNumber, handshake_id::UInt8)
    if !isempty(version.prerelease) || !isempty(version.build)
        throw(ArgumentError("Version must only consist of major.minor.patch components"))
    end
    if length(appname) > typemax(UInt8)
        throw(ArgumentError("Application name too long ($(typemax(UInt8)) characters)"))
    end

    appname_str = string(appname)

    write(buffer, UInt8(length(appname_str)))
    write(buffer, appname_str)
    write(buffer, UInt8.([version.major, version.minor, version.patch]))
    write(buffer, handshake_id)
end

function packetcontent(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeStart}, bytes)
    buffer = IOBuffer(bytes)
    appname_len = read(buffer, UInt8)
    appname     = read(buffer, appname_len)
    major, minor, patch = read(buffer, 3)
    handshake_id = read(buffer, UInt8)
    @assert isempty(buffer) "unexpected excess packet data"
    return unsafe_string(pointer(appname)), VersionNumber(major, minor, patch), handshake_id
end

# --- Handshake Response ---
function packetbytes(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeResponse}, buffer::IOBuffer, handshake_id::UInt8)
    write(buffer, handshake_id)
end

function packetcontent(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeResponse}, bytes)
    buffer = IOBuffer(bytes)
    handshake_incr_id = read(buffer, UInt8)
    @assert isempty(buffer) "unexpected excess packet data"
    return handshake_incr_id
end

# --- Handshake Confirm ---
function packetbytes(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeConfirm}, buffer::IOBuffer)
    nothing
end

function packetcontent(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.HandshakeConfirm}, bytes)
    nothing
end

# --- Goodbye ---
function packetbytes(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.Goodbye}, buffer::IOBuffer)
    nothing
end

function packetcontent(::Ident{Packets.ConnectivityClass}, ::Ident{Packets.Goodbye}, bytes)
    nothing
end

# ===== Discovery Class Packets =====
# --- Query Peer ---
function packetbytes(::Ident{Packets.DiscoveryClass}, ::Ident{Packets.QueryPeer}, buffer::IOBuffer, address::Address{S}) where S
    @assert S isa Integer "invalid address length parameter"
    write(buffer, hton(UInt32(S)))
    write(buffer, address)
end

function packetcontent(::Ident{Packets.DiscoveryClass}, ::Ident{Packets.QueryPeer}, bytes)
    S = ntoh(read(buffer, UInt32))
    addr = read(buffer, Address{S})
    return addr
end

# ===== Utility Class Packets =====
# --- Wakeup Loop ---
function packetbytes(::Ident{Packets.UtilityClass}, ::Ident{Packets.Wakeup}, buffer::IOBuffer)
    nothing
end

function packetcontent(::Ident{Packets.UtilityClass}, ::Ident{Packets.Wakeup}, bytes)
    nothing
end

# ===== Data Class Packets =====
# --- AppData ---
function packetbytes(::Ident{Packets.DataClass}, ::Ident{Packets.AppData}, buffer::IOBuffer, data)
    write(buffer, data)
end

function packetcontent(::Ident{Packets.DataClass}, ::Ident{Packets.AppData}, bytes)
    return bytes
end

end # module PacketFactories


text:
# This example demonstrated how to calculate the bubble diagram of free electrons using the Monte Carlo module

using QuantumStatistics, LinearAlgebra, Random, Printf, BenchmarkTools, InteractiveUtils, Parameters
# using ProfileView

const Steps = 1e8

include("parameter.jl")

@with_kw struct Para
    n::Int = 0 # external Matsubara frequency
    Qsize::Int = 16
    extQ::Vector{SVector{3,Float64}} = [@SVector [q, 0.0, 0.0] for q in LinRange(0.0, 3.0 * kF, Qsize)]
end

function integrand(config)
    if config.curr != 1
        error("impossible")
    end
    para = config.para

    T, K, Ext = config.var[1], config.var[2], config.var[3]
    k = K[1]
    Tin, Tout = 0.0, T[1]
    extidx = Ext[1]
    q = para.extQ[extidx] # external momentum
    kq = k + q
    τ = (Tout - Tin) / β
    ω1 = (dot(k, k) - kF^2) / (2me) * β
    g1 = Spectral.kernelFermiT(τ, ω1)
    ω2 = (dot(kq, kq) - kF^2) / (2me) * β
    g2 = Spectral.kernelFermiT(-τ, ω2)
    phase = 1.0 / (2π)^3
    return g1 * g2 * spin * phase * cos(2π * para.n * τ)
end

function measure(config)
    obs = config.observable
    factor = 1.0 / config.reweight[config.curr]
    extidx = config.var[3][1]
    weight = integrand(config)
    obs[extidx] += weight / abs(weight) * factor
end

function run(steps)

    para = Para()
    @unpack extQ, Qsize = para

    T = MonteCarlo.Tau(β, β / 2.0)
    K = MonteCarlo.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
    Ext = MonteCarlo.Discrete(1, length(extQ)) # external variable is specified

    dof = [[1, 1, 1],] # degrees of freedom of the normalization diagram and the bubble
    obs = zeros(Float64, Qsize) # observable for the normalization diagram and the bubble

    config = MonteCarlo.Configuration(steps, (T, K, Ext), dof, obs; para=para)
    avg, std = MonteCarlo.sample(config, integrand, measure; print=0, Nblock=16)
    # @profview MonteCarlo.sample(config, integrand, measure; print=0, Nblock=1)
    # sleep(100)

    if isnothing(avg) == false
        @unpack n, extQ = Para()

        for (idx, q) in enumerate(extQ)
            q = q[1]
            p, err = TwoPoint.LindhardΩnFiniteTemperature(dim, q, n, kF, β, me, spin)
            @printf("%10.6f  %10.6f ± %10.6f  %10.6f ± %10.6f\n", q / kF, avg[idx], std[idx], p, err)
        end

    end
end

run(Steps)
# @time run(Steps)


text:
using BarycentricLagrange
using Compat
using Compat.Test
using Compat.LinearAlgebra

@testset "BarycentricLagrange" begin

    function dotest(F)
        # Chebyshev source nodes
        n = 100
        x = cos.((0:n)*pi/n)
        w = bclag_interp_weights(x)
        # Equispaced evluation nodes
        xx = collect(-1:0.01:1)

        f = F.(x)
        ref = F.(xx)

        @test norm(bclag_interp_eval(x, f, xx)-ref, Inf) <  2e-15
        @test norm(bclag_interp_eval(x, f, xx, w)-ref, Inf) <  2e-15

        M1 = bclag_interp_matrix(x, xx)
        @test norm(M1*f-ref, Inf) <  2e-15

        M2 = bclag_interp_matrix(x, xx, w)
        @test norm(M2*f-ref, Inf) <  2e-15
    end

    @testset "Real-valued" begin
        F = x -> cos(x)
        dotest(F)
    end

    @testset "Complex-valued" begin
        F = x -> exp(1im*x)
        dotest(F)
    end

end


text:
# -*- mode: julia; -*-
module WAV
export wavread, wavwrite, wavplay, WAVE_FORMAT_PCM, WAVE_FORMAT_IEEE_FLOAT, WAVE_FORMAT_ALAW, WAVE_FORMAT_MULAW
import Base.unbox, Base.box

@linux? include("wavplay.jl") : (wavplay(args...) = warn("wavplay is not currently implemented on $OS_NAME"))

# The WAV specification states that numbers are written to disk in little endian form.
write_le(stream::IO, value::Uint8) = write(stream, value)

function write_le(stream::IO, value::Uint16)
    write(stream, uint8((value & 0x00ff)     ))
    write(stream, uint8((value & 0xff00) >> 8))
end

function write_le(stream::IO, value::Uint32)
    write(stream, uint8((value & 0x000000ff)      ))
    write(stream, uint8((value & 0x0000ff00) >>  8))
    write(stream, uint8((value & 0x00ff0000) >> 16))
    write(stream, uint8((value & 0xff000000) >> 24))
end

function write_le(stream::IO, value::Uint64)
    write(stream, uint8((value & 0x00000000000000ff)      ))
    write(stream, uint8((value & 0x000000000000ff00) >>  8))
    write(stream, uint8((value & 0x0000000000ff0000) >> 16))
    write(stream, uint8((value & 0x00000000ff000000) >> 24))
    write(stream, uint8((value & 0x000000ff00000000) >> 32))
    write(stream, uint8((value & 0x0000ff0000000000) >> 40))
    write(stream, uint8((value & 0x00ff000000000000) >> 48))
    write(stream, uint8((value & 0xff00000000000000) >> 56))
end

write_le(stream::IO, value::Int16) = write_le(stream, uint16(value))
write_le(stream::IO, value::Int32) = write_le(stream, uint32(value))
write_le(stream::IO, value::Int64) = write_le(stream, uint64(value))
write_le(stream::IO, value::Float32) = write_le(stream, box(Uint32, unbox(Float32, value)))
write_le(stream::IO, value::Float64) = write_le(stream, box(Uint64, unbox(Float64, value)))

read_le(stream::IO, x::Type{Uint8}) = read(stream, x)

function read_le(stream::IO, ::Type{Uint16})
    bytes = uint16(read(stream, Uint8, 2))
    (bytes[2] << 8) | bytes[1]
end

function read_le(stream::IO, ::Type{Uint32})
    bytes = uint32(read(stream, Uint8, 4))
    (bytes[4] << 24) | (bytes[3] << 16) | (bytes[2] << 8) | bytes[1]
end

function read_le(stream::IO, ::Type{Uint64})
    bytes = uint64(read(stream, Uint8, 8))
    (bytes[8] << 56) | (bytes[7] << 48) | (bytes[6] << 40) | (bytes[5] << 32) | (bytes[4] << 24) | (bytes[3] << 16) | (bytes[2] << 8) | bytes[1]
end

read_le(stream::IO, ::Type{Int16}) = int16(read_le(stream, Uint16))
read_le(stream::IO, ::Type{Int32}) = int32(read_le(stream, Uint32))
read_le(stream::IO, ::Type{Int64}) = int64(read_le(stream, Uint64))
read_le(stream::IO, ::Type{Float32}) = box(Float32, unbox(Uint32, read_le(stream, Uint32)))
read_le(stream::IO, ::Type{Float64}) = box(Float64, unbox(Uint64, read_le(stream, Uint64)))

# Required WAV Chunk; The format chunk describes how the waveform data is stored
type WAVFormat
    compression_code::Uint16
    nchannels::Uint16
    sample_rate::Uint32
    bps::Uint32 # average bytes per second
    block_align::Uint16
    nbits::Uint16
    extra_bytes::Array{Uint8}

    data_length::Uint32
end
WAVFormat() = WAVFormat(uint16(0), uint16(0), uint32(0), uint32(0), uint16(0), uint16(0), Array(Uint8), uint32(0))
WAVFormat(comp, chan, fs, bytes, ba, nbits) = WAVFormat(comp, chan, fs, bytes, ba, nbits, Array(Uint8), uint32(0))

const WAVE_FORMAT_PCM        = 0x0001 # PCM
const WAVE_FORMAT_IEEE_FLOAT = 0x0003 # IEEE float
const WAVE_FORMAT_ALAW       = 0x0006 # A-Law
const WAVE_FORMAT_MULAW      = 0x0007 # Mu-Law
const WAVE_FORMAT_EXTENSIBLE = 0xfffe # Extension!

# used by WAVE_FORMAT_EXTENSIBLE
type WAVFormatExtension
    valid_bits_per_sample::Uint16
    channel_mask::Uint32
    sub_format::Array{Uint8} # 16 byte GUID
end
WAVFormatExtension() = WAVFormatExtension(uint16(0), uint32(0), b"")

# DEFINE_GUIDSTRUCT("00000001-0000-0010-8000-00aa00389b71", KSDATAFORMAT_SUBTYPE_PCM);
const KSDATAFORMAT_SUBTYPE_PCM = [
0x01, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                  ]
# DEFINE_GUIDSTRUCT("00000003-0000-0010-8000-00aa00389b71", KSDATAFORMAT_SUBTYPE_IEEE_FLOAT);
const KSDATAFORMAT_SUBTYPE_IEEE_FLOAT = [
0x03, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                         ]
# DEFINE_GUID(KSDATAFORMAT_SUBTYPE_MULAW, 0x00000007, 0x0000, 0x0010, 0x80, 0x00, 0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71);
const KSDATAFORMAT_SUBTYPE_MULAW = [
0x07, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                    ]
#DEFINE_GUID(KSDATAFORMAT_SUBTYPE_ALAW, 0x00000006, 0x0000, 0x0010, 0x80, 0x00, 0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71);
const KSDATAFORMAT_SUBTYPE_ALAW = [
0x06, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                   ]

function WAVFormatExtension(bytes::Array{Uint8})
    if length(bytes) != 22
        error("There are not the right number of bytes for the WAVFormat Extension")
    end
    # split bytes into valid_bits_per_sample, channel_mask, and sub_format
    valid_bits_per_sample = (uint16(bytes[2]) << 8) | uint16(bytes[1])
    channel_mask = (uint32(bytes[6]) << 24) | (uint32(bytes[5]) << 16) | (uint32(bytes[4]) << 8) | uint32(bytes[3])
    sub_format = bytes[7:end]
    return WAVFormatExtension(valid_bits_per_sample, channel_mask, sub_format)
end

function read_header(io::IO)
    # check if the given file has a valid RIFF header
    riff = read(io, Uint8, 4)
    if riff !=  b"RIFF"
        error("$filename is not a valid WAV file: The RIFF header is invalid")
    end

    chunk_size = read_le(io, Uint32)

    # check if this is a WAV file
    format = read(io, Uint8, 4)
    if format != b"WAVE"
        error("$filename is not a valid WAV file: the format is not WAVE")
    end
    return chunk_size
end

function write_header(io::IO, fmt::WAVFormat)
    write(io, b"RIFF") # RIFF header
    write_le(io, uint32(36 + fmt.data_length)) # chunk_size
    write(io, b"WAVE")
end

function write_header(io::IO, fmt::WAVFormat, ext_fmt::WAVFormatExtension)
    write(io, b"RIFF") # RIFF header
    write_le(io, uint32(60 + fmt.data_length)) # chunk_size
    write(io, b"WAVE")
end

function read_format(io::IO, chunk_size::Uint32)
    # can I read in all of the fields at once?
    orig_chunk_size = int(chunk_size)
    if chunk_size < 16
        error("The WAVE Format chunk must be at least 16 bytes")
    end
    format = WAVFormat(read_le(io, Uint16), # Compression Code
                       read_le(io, Uint16), # Number of Channels
                       read_le(io, Uint32), # Sample Rate
                       read_le(io, Uint32), # bytes per second
                       read_le(io, Uint16), # block align
                       read_le(io, Uint16)) # bits per sample
    chunk_size -= 16
    if chunk_size > 0
        # TODO add error checking for size mismatches
        extra_bytes = read_le(io, Uint16)
        format.extra_bytes = read(io, Uint8, extra_bytes)
    end
    return format
end

function write_format(io::IO, fmt::WAVFormat, ext_length::Integer)
    # write the fmt subchunk header
    write(io, b"fmt ")
    write_le(io, uint32(16 + ext_length)) # subchunk length; 16 is size of base format chunk

    write_le(io, fmt.compression_code) # audio format (Uint16)
    write_le(io, fmt.nchannels) # number of channels (Uint16)
    write_le(io, fmt.sample_rate) # sample rate (Uint32)
    write_le(io, fmt.bps) # byte rate (Uint32)
    write_le(io, fmt.block_align) # byte align (Uint16)
    write_le(io, fmt.nbits) # number of bits per sample (UInt16)
end
write_format(io::IO, fmt::WAVFormat) = write_format(io, fmt, 0)

function write_format(io::IO, fmt::WAVFormat, ext::WAVFormatExtension)
    write_format(io, fmt, 24) # 24 is the added length needed to encode the extension
    write_le(io, uint16(22))
    write_le(io, ext.valid_bits_per_sample)
    write_le(io, ext.channel_mask)
    @assert length(ext.sub_format) == 16
    write(io, ext.sub_format)
end

function pcm_container_type(nbits::Unsigned)
    if nbits > 32
        return Int64
    elseif nbits > 16
        return Int32
    elseif nbits > 8
        return Int16
    end
    return  Uint8
end

ieee_float_container_type(nbits::Unsigned) = (nbits == 32 ? Float32 : (nbits == 64 ? Float64 : error("$nbits bits is not supported for WAVE_FORMAT_IEEE_FLOAT.")))

function read_pcm_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange::Range1)
    samples = Array(pcm_container_type(fmt.nbits), length(subrange), fmt.nchannels)
    if isempty(subrange)
        return samples
    end
    const nbytes = iceil(fmt.nbits / 8)
    const bitshift::Array{Uint} = linspace(0, 64, 9)
    const mask = unsigned(1) << (fmt.nbits - 1)
    const signextend_mask = ~unsigned(0) << fmt.nbits
    skip(io, uint((first(subrange) - 1) * nbytes * fmt.nchannels))
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            raw_sample = read(io, Uint8, nbytes)
            my_sample = uint64(0)
            for k = 1:nbytes
                my_sample |= uint64(raw_sample[k]) << bitshift[k]
            end
            my_sample >>= nbytes * 8 - fmt.nbits
            # sign extend negative values
            if fmt.nbits > 8 && (my_sample & mask > 0)
                my_sample |= signextend_mask
            end
            samples[i, j] = convert(eltype(samples), my_sample)
        end
    end
    samples
end

# "subrange" is a placeholder for None in this method
function read_pcm_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange)
    const nblocks = uint(chunk_size / fmt.block_align) # each block stores fmt.nchannels channels
    if nblocks == 0
        return Array(pcm_container_type(fmt.nbits), 0, fmt.nchannels)
    end
    read_pcm_samples(io, chunk_size, fmt, 1:nblocks)
end

function read_ieee_float_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange::Range1)
    const floatType = ieee_float_container_type(fmt.nbits)
    if isempty(subrange)
        return Array(floatType, nblocks, fmt.nchannels)
    end
    const nblocks = length(subrange)
    samples = Array(floatType, nblocks, fmt.nchannels)
    skip(io, uint((first(subrange) - 1) * (fmt.nbits / 8) * fmt.nchannels))
    for i = 1:nblocks
        for j = 1:fmt.nchannels
            samples[i, j] = read_le(io, floatType)
        end
    end
    samples
end

function read_companded_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange, table::Array{Int16})
    const nblocks = uint(chunk_size / fmt.block_align) # each block stores fmt.nchannels channels
    if nblocks == 0
        return Array(Int16, 0, fmt.nchannels)
    end
    read_companded_samples(io, chunk_size, fmt, 1:nblocks, table)
end

function read_companded_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange::Range1, table::Array{Int16})
    if isempty(subrange)
        return Array(Int16, 0, fmt.nchannels)
    end
    const nblocks = length(subrange)
    samples = Array(Int16, nblocks, fmt.nchannels)
    skip(io, uint((first(subrange) - 1) * fmt.nchannels))
    for i = 1:nblocks
        for j = 1:fmt.nchannels
            # add one to value from blocks because A-law stores values from 0 to 255.
            const compressedByte::Uint8 = clamp(read_le(io, Uint8), 0, 255)
            # Julia indexing is 1-based; I need a value from 1 to 256
            samples[i, j] = table[compressedByte + 1]
        end
    end
    return samples
end

function read_mulaw_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange)
    # Quantized μ-law algorithm -- Use a look up table to convert
    # From Wikipedia, ITU-T Recommendation G.711 and G.191 specify the following intervals:
    #
    # ---------------------------------------+--------------------------------
    #  14 bit Binary Linear input code       | 8 bit Compressed code
    # ---------------------------------------+--------------------------------
    # +8158 to +4063 in 16 intervals of 256  |  0x80 + interval number
    # +4062 to +2015 in 16 intervals of 128  |  0x90 + interval number
    # +2014 to +991 in 16 intervals of 64    |  0xA0 + interval number
    # +990 to +479 in 16 intervals of 32     |  0xB0 + interval number
    # +478 to +223 in 16 intervals of 16     |  0xC0 + interval number
    # +222 to +95 in 16 intervals of 8       |  0xD0 + interval number
    # +94 to +31 in 16 intervals of 4        |  0xE0 + interval number
    # +30 to +1 in 15 intervals of 2         |  0xF0 + interval number
    # 0                                      |  0xFF
    # −1                                     |  0x7F
    # −31 to −2 in 15 intervals of 2         |  0x70 + interval number
    # −95 to −32 in 16 intervals of 4        |  0x60 + interval number
    # −223 to −96 in 16 intervals of 8       |  0x50 + interval number
    # −479 to −224 in 16 intervals of 16     |  0x40 + interval number
    # −991 to −480 in 16 intervals of 32     |  0x30 + interval number
    # −2015 to −992 in 16 intervals of 64    |  0x20 + interval number
    # −4063 to −2016 in 16 intervals of 128  |  0x10 + interval number
    # −8159 to −4064 in 16 intervals of 256  |  0x00 + interval number
    # ---------------------------------------+--------------------------------
    const MuLawDecompressTable::Array{Int16} =
    [
    -32124,-31100,-30076,-29052,-28028,-27004,-25980,-24956,
    -23932,-22908,-21884,-20860,-19836,-18812,-17788,-16764,
    -15996,-15484,-14972,-14460,-13948,-13436,-12924,-12412,
    -11900,-11388,-10876,-10364, -9852, -9340, -8828, -8316,
    -7932, -7676, -7420, -7164, -6908, -6652, -6396, -6140,
    -5884, -5628, -5372, -5116, -4860, -4604, -4348, -4092,
    -3900, -3772, -3644, -3516, -3388, -3260, -3132, -3004,
    -2876, -2748, -2620, -2492, -2364, -2236, -2108, -1980,
    -1884, -1820, -1756, -1692, -1628, -1564, -1500, -1436,
    -1372, -1308, -1244, -1180, -1116, -1052,  -988,  -924,
    -876,  -844,  -812,  -780,  -748,  -716,  -684,  -652,
    -620,  -588,  -556,  -524,  -492,  -460,  -428,  -396,
    -372,  -356,  -340,  -324,  -308,  -292,  -276,  -260,
    -244,  -228,  -212,  -196,  -180,  -164,  -148,  -132,
    -120,  -112,  -104,   -96,   -88,   -80,   -72,   -64,
    -56,   -48,   -40,   -32,   -24,   -16,    -8,     -1,
    32124, 31100, 30076, 29052, 28028, 27004, 25980, 24956,
    23932, 22908, 21884, 20860, 19836, 18812, 17788, 16764,
    15996, 15484, 14972, 14460, 13948, 13436, 12924, 12412,
    11900, 11388, 10876, 10364,  9852,  9340,  8828,  8316,
    7932,  7676,  7420,  7164,  6908,  6652,  6396,  6140,
    5884,  5628,  5372,  5116,  4860,  4604,  4348,  4092,
    3900,  3772,  3644,  3516,  3388,  3260,  3132,  3004,
    2876,  2748,  2620,  2492,  2364,  2236,  2108,  1980,
    1884,  1820,  1756,  1692,  1628,  1564,  1500,  1436,
    1372,  1308,  1244,  1180,  1116,  1052,   988,   924,
    876,   844,   812,   780,   748,   716,   684,   652,
    620,   588,   556,   524,   492,   460,   428,   396,
    372,   356,   340,   324,   308,   292,   276,   260,
    244,   228,   212,   196,   180,   164,   148,   132,
    120,   112,   104,    96,    88,    80,    72,    64,
    56,    48,    40,    32,    24,    16,     8,     0
     ]
    @assert length(MuLawDecompressTable) == 256
    return read_companded_samples(io, chunk_size, fmt, subrange, MuLawDecompressTable)
end

function read_alaw_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange)
    # Quantized A-law algorithm -- Use a look up table to convert
    const ALawDecompressTable::Array{Int16} =
    [
    -5504, -5248, -6016, -5760, -4480, -4224, -4992, -4736,
    -7552, -7296, -8064, -7808, -6528, -6272, -7040, -6784,
    -2752, -2624, -3008, -2880, -2240, -2112, -2496, -2368,
    -3776, -3648, -4032, -3904, -3264, -3136, -3520, -3392,
    -22016,-20992,-24064,-23040,-17920,-16896,-19968,-18944,
    -30208,-29184,-32256,-31232,-26112,-25088,-28160,-27136,
    -11008,-10496,-12032,-11520,-8960, -8448, -9984, -9472,
    -15104,-14592,-16128,-15616,-13056,-12544,-14080,-13568,
    -344,  -328,  -376,  -360,  -280,  -264,  -312,  -296,
    -472,  -456,  -504,  -488,  -408,  -392,  -440,  -424,
    -88,   -72,   -120,  -104,  -24,   -8,    -56,   -40,
    -216,  -200,  -248,  -232,  -152,  -136,  -184,  -168,
    -1376, -1312, -1504, -1440, -1120, -1056, -1248, -1184,
    -1888, -1824, -2016, -1952, -1632, -1568, -1760, -1696,
    -688,  -656,  -752,  -720,  -560,  -528,  -624,  -592,
    -944,  -912,  -1008, -976,  -816,  -784,  -880,  -848,
    5504,  5248,  6016,  5760,  4480,  4224,  4992,  4736,
    7552,  7296,  8064,  7808,  6528,  6272,  7040,  6784,
    2752,  2624,  3008,  2880,  2240,  2112,  2496,  2368,
    3776,  3648,  4032,  3904,  3264,  3136,  3520,  3392,
    22016, 20992, 24064, 23040, 17920, 16896, 19968, 18944,
    30208, 29184, 32256, 31232, 26112, 25088, 28160, 27136,
    11008, 10496, 12032, 11520, 8960,  8448,  9984,  9472,
    15104, 14592, 16128, 15616, 13056, 12544, 14080, 13568,
    344,   328,   376,   360,   280,   264,   312,   296,
    472,   456,   504,   488,   408,   392,   440,   424,
    88,    72,   120,   104,    24,     8,    56,    40,
    216,   200,   248,   232,   152,   136,   184,   168,
    1376,  1312,  1504,  1440,  1120,  1056,  1248,  1184,
    1888,  1824,  2016,  1952,  1632,  1568,  1760,  1696,
    688,   656,   752,   720,   560,   528,   624,   592,
    944,   912,  1008,   976,   816,   784,   880,   848
     ]
    @assert length(ALawDecompressTable) == 256
    return read_companded_samples(io, chunk_size, fmt, subrange, ALawDecompressTable)
end

function compress_sample_mulaw(sample::Int16)
    const MuLawCompressTable::Array{Uint8} =
    [
    0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
    4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
     ]
    @assert length(MuLawCompressTable) == 256
    const cBias::Int16 = 0x84
    const cClip::Int16 = 32635

    const sign = (sample >>> 8) & 0x80
    if sign != 0
        sample = -sample
    end
    if sample > cClip
        sample = cClip
    end
    sample = sample + cBias
    const exponent::Uint8 = MuLawCompressTable[(sample >>> 7) + 1]
    const mantissa = (sample >> (exponent+3)) & 0x0F
    const compressedByte = ~ (sign | (exponent << 4) | mantissa)
    return compressedByte
end

function compress_sample_alaw(sample::Int16)
    const ALawCompressTable::Array{Uint8} =
    [
    1,1,2,2,3,3,3,3,
    4,4,4,4,4,4,4,4,
    5,5,5,5,5,5,5,5,
    5,5,5,5,5,5,5,5,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7
     ]
    @assert length(ALawCompressTable) == 128
    const cBias::Int16 = 0x84
    const cClip::Int16 = 32635
    const sign::Int = ((~sample >>> 8) & 0x80)
    if sign == 0
        sample = -sample
    end
    if sample > cClip
        sample = cClip
    end
    compressedByte::Uint8 = 0
    if sample >= 256
        const exponent::Int = ALawCompressTable[((sample >>> 8) & 0x7f) + 1]
        const mantissa::Int = (sample >>> (exponent + 3) ) & 0x0f
        compressedByte = uint8((exponent << 4) | mantissa)
    else
        compressedByte = uint8(sample >>> 4)
    end
    compressedByte $= (sign $ 0x55)
end

function write_companded_samples{T<:Integer}(io::IO, fmt::WAVFormat, samples::Array{T}, compander::Function)
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            const compressedByte::Uint8 = compander(samples[i, j])
            write_le(io, compressedByte)
        end
    end
end

function write_companded_samples{T<:FloatingPoint}(io::IO, fmt::WAVFormat, samples::Array{T}, compander::Function)
    samples = convert(Array{Int16}, round(samples * typemax(Int16)))
    write_companded_samples(io, fmt, samples, compander)
end

function read_ieee_float_samples(io::IO, chunk_size::Unsigned, fmt::WAVFormat, subrange)
    const nblocks = uint(chunk_size / fmt.block_align) # each block stores fmt.nchannels channels
    if nblocks == 0
        return Array(ieee_float_container_type(fmt.nbits), 0, fmt.nchannels)
    end
    read_ieee_float_samples(io, chunk_size, fmt, 1:nblocks)
end

# PCM data is two's-complement except for resolutions of 1-8 bits, which are represented as offset binary.

# support every bit width from 1 to 8 bits
convert_pcm_to_double(samples::Array{Uint8}, nbits::Integer) = convert(Array{Float64}, samples) ./ (2^nbits - 1) .* 2.0 .- 1.0
convert_pcm_to_double(samples::Array{Int8}, nbits::Integer) = error("WAV files use offset binary for less than 9 bits")
# support every bit width from 9 to 64 bits
convert_pcm_to_double{T<:Signed}(samples::Array{T}, nbits::Integer) = convert(Array{Float64}, samples) / (2^(nbits - 1) - 1)

function read_data(io::IO, chunk_size::Uint32, fmt::WAVFormat, format::String, subrange::Any)
    # "format" is the format of values, while "fmt" is the WAV file level format
    samples = None
    convert_to_double = x -> convert(Array{Float64}, x)
    if fmt.compression_code == WAVE_FORMAT_EXTENSIBLE
        ext_fmt = WAVFormatExtension(fmt.extra_bytes)
        if ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_PCM
            fmt.nbits = ext_fmt.valid_bits_per_sample
            samples = read_pcm_samples(io, chunk_size, fmt, subrange)
            convert_to_double = x -> convert_pcm_to_double(x, fmt.nbits)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_IEEE_FLOAT
            fmt.nbits = ext_fmt.valid_bits_per_sample
            samples = read_ieee_float_samples(io, chunk_size, fmt, subrange)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_ALAW
            fmt.nbits = 8
            samples = read_alaw_samples(io, chunk_size, fmt, subrange)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_MULAW
            fmt.nbits = 8
            samples = read_mulaw_samples(io, chunk_size, fmt, subrange)
        else
            error("$ext_fmt -- WAVE_FORMAT_EXTENSIBLE Not done yet!")
        end
    elseif fmt.compression_code == WAVE_FORMAT_PCM
        samples = read_pcm_samples(io, chunk_size, fmt, subrange)
        convert_to_double = x -> convert_pcm_to_double(x, fmt.nbits)
    elseif fmt.compression_code == WAVE_FORMAT_IEEE_FLOAT
        samples = read_ieee_float_samples(io, chunk_size, fmt, subrange)
    elseif fmt.compression_code == WAVE_FORMAT_MULAW
        samples = read_mulaw_samples(io, chunk_size, fmt, subrange)
        convert_to_double = x -> convert_pcm_to_double(x, 16)
    elseif fmt.compression_code == WAVE_FORMAT_ALAW
        samples = read_alaw_samples(io, chunk_size, fmt, subrange)
        convert_to_double = x -> convert_pcm_to_double(x, 16)
    else
        error("$(fmt.compression_code) is an unsupported compression code!")
    end
    if format == "double"
        samples = convert_to_double(samples)
    end
    samples
end

function write_pcm_samples{T<:Integer}(io::IO, fmt::WAVFormat, samples::Array{T})
    # number of bytes per sample
    const nbytes = iceil(fmt.nbits / 8)
    const bitshift::Array{Uint} = linspace(0, 64, 9)
    const minval = fmt.nbits > 8 ? -2^(fmt.nbits - 1) : -2^(fmt.nbits)
    const maxval = fmt.nbits > 8 ? 2^(fmt.nbits - 1) - 1 : 2^(fmt.nbits) - 1
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            my_sample = clamp(samples[i, j], minval, maxval)
            # shift my_sample into the N most significant bits
            my_sample <<= nbytes * 8 - fmt.nbits
            mask = uint64(0xff)
            for k = 1:nbytes
                write_le(io, uint8((my_sample & mask) >> bitshift[k]))
                mask <<= 8
            end
        end
    end
end

function write_pcm_samples{T<:FloatingPoint}(io::IO, fmt::WAVFormat, samples::Array{T})
    # Scale the floating point values to the PCM range
    if fmt.nbits > 8
        # two's complement
        samples = convert(Array{pcm_container_type(fmt.nbits)}, round(samples * (2^(fmt.nbits - 1) - 1)))
    else
        # offset binary
        samples = convert(Array{Uint8}, round((samples .+ 1.0) / 2.0 * (2^fmt.nbits - 1)))
    end
    return write_pcm_samples(io, fmt, samples)
end

function write_ieee_float_samples{T<:FloatingPoint}(io::IO, fmt::WAVFormat, samples::Array{T})
    const floatType = ieee_float_container_type(fmt.nbits)
    samples = convert(Array{floatType}, samples)
    const minval = convert(floatType, -1.0)
    const maxval = convert(floatType, 1.0)
    # Interleave the channel samples before writing to the stream.
    for i = 1:size(samples, 1) # for each sample
        for j = 1:size(samples, 2) # for each channel
            write_le(io, clamp(samples[i, j], minval, maxval))
        end
    end
end

function write_data(io::IO, fmt::WAVFormat, ext_fmt::WAVFormatExtension, samples::Array)
    if fmt.compression_code == WAVE_FORMAT_EXTENSIBLE
        if ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_PCM
            fmt.nbits = ext_fmt.valid_bits_per_sample
            return write_pcm_samples(io, fmt, samples)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_IEEE_FLOAT
            fmt.nbits = ext_fmt.valid_bits_per_sample
            return write_ieee_float_samples(io, fmt, samples)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_ALAW
            fmt.nbits = 8
            return write_companded_samples(io, fmt, samples, compress_sample_alaw)
        elseif ext_fmt.sub_format == KSDATAFORMAT_SUBTYPE_MULAW
            fmt.nbits = 8
            return write_companded_samples(io, fmt, samples, compress_sample_mulaw)
        else
            error("$ext_fmt -- WAVE_FORMAT_EXTENSIBLE Not done yet!")
        end
    elseif fmt.compression_code == WAVE_FORMAT_PCM
        return write_pcm_samples(io, fmt, samples)
    elseif fmt.compression_code == WAVE_FORMAT_IEEE_FLOAT
        return write_ieee_float_samples(io, fmt, samples)
    elseif fmt.compression_code == WAVE_FORMAT_MULAW
        return write_companded_samples(io, fmt, samples, compress_sample_mulaw)
    elseif fmt.compression_code == WAVE_FORMAT_ALAW
        return write_companded_samples(io, fmt, samples, compress_sample_alaw)
    else
        error("$(fmt.compression_code) is an unsupported compression code.")
    end
end

make_range(subrange) = subrange
make_range(subrange::Number) = 1:convert(Int, subrange)
make_range(subrange::Range1) = convert(Range1{Int}, subrange)

function wavread(io::IO; subrange=None, format="double")
    chunk_size = read_header(io)
    fmt = WAVFormat()
    samples = Array(Float64)

    # Note: This assumes that the format chunk is written in the file before the data chunk. The
    # specification does not require this assumption, but most real files are written that way.

    # Subtract the size of the format field from chunk_size; now it holds the size
    # of all the sub-chunks
    chunk_size -= 4
    while chunk_size > 0
        # Read subchunk ID and size
        subchunk_id = read(io, Uint8, 4)
        subchunk_size = read_le(io, Uint32)
        chunk_size -= 8 + subchunk_size
        # check the subchunk ID
        if subchunk_id == b"fmt "
            fmt = read_format(io, subchunk_size)
        elseif subchunk_id == b"data"
            if format == "size"
                return int(subchunk_size / fmt.block_align), int(fmt.nchannels)
            end
            samples = read_data(io, subchunk_size, fmt, format, make_range(subrange))
        else
            # return unknown sub-chunks?
            # Note: Ignoring unknown sub chunks for now
            skip(io, subchunk_size)
        end
    end
    return samples, fmt.sample_rate, fmt.nbits, None
end

function wavread(filename::String; subrange=Any, format="double")
    io = open(filename, "r")
    finalizer(io, close)
    const result = wavread(io, subrange=subrange, format=format)
    close(io)
    return result
end

# These are the MATLAB compatible signatures
wavread(filename::String, fmt::String) = wavread(filename, format=fmt)
wavread(filename::String, N::Int) = wavread(filename, subrange=N)
wavread(filename::String, N::Range1) = wavread(filename, subrange=N)
wavread(filename::String, N::Int, fmt::String) = wavread(filename, subrange=N, format=fmt)
wavread(filename::String, N::Range1, fmt::String) = wavread(filename, subrange=N, format=fmt)

get_default_compression{T<:Integer}(samples::Array{T}) = WAVE_FORMAT_PCM
get_default_compression{T<:FloatingPoint}(samples::Array{T}) = WAVE_FORMAT_IEEE_FLOAT
get_default_pcm_precision(samples::Array{Uint8}) = 8
get_default_pcm_precision(samples::Array{Int16}) = 16
get_default_pcm_precision(samples) = 24

function get_default_precision(samples, compression)
    if compression == WAVE_FORMAT_ALAW || compression == WAVE_FORMAT_MULAW
        return 8
    elseif compression == WAVE_FORMAT_IEEE_FLOAT
        return 32
    end
    get_default_pcm_precision(samples)
end

function wavwrite(samples::Array, io::IO; Fs=8000, nbits=0, compression=0)
    if compression == 0
        compression = get_default_compression(samples)
    elseif compression == WAVE_FORMAT_ALAW || compression == WAVE_FORMAT_MULAW
        nbits = 8
    end
    if nbits == 0
        nbits = get_default_precision(samples, compression)
    end
    fmt = WAVFormat()
    fmt.compression_code = compression
    fmt.nchannels = size(samples, 2)
    fmt.sample_rate = Fs
    fmt.nbits = iceil(nbits / 8) * 8
    fmt.block_align = fmt.nbits / 8 * fmt.nchannels
    fmt.bps = fmt.sample_rate * fmt.block_align
    fmt.data_length = size(samples, 1) * fmt.block_align

    ext = WAVFormatExtension()
    if fmt.nchannels > 2 || fmt.nbits > 16 || fmt.nbits != nbits
        fmt.compression_code = WAVE_FORMAT_EXTENSIBLE
        ext.valid_bits_per_sample = nbits
        ext.channel_mask = 0
        if compression == WAVE_FORMAT_PCM
            ext.sub_format = KSDATAFORMAT_SUBTYPE_PCM
        elseif compression == WAVE_FORMAT_IEEE_FLOAT
            ext.sub_format = KSDATAFORMAT_SUBTYPE_IEEE_FLOAT
        elseif compression == WAVE_FORMAT_ALAW
            ext.sub_format = KSDATAFORMAT_SUBTYPE_ALAW
        elseif compression == WAVE_FORMAT_MULAW
            ext.sub_format = KSDATAFORMAT_SUBTYPE_MULAW
        else
            error("Unsupported extension sub format: $compression")
        end
        write_header(io, fmt, ext)
        write_format(io, fmt, ext)
    else
        write_header(io, fmt)
        write_format(io, fmt)
    end

    # write the data subchunk header
    write(io, b"data")
    write_le(io, fmt.data_length) # Uint32
    write_data(io, fmt, ext, samples)
end

function wavwrite(samples::Array, filename::String; Fs=8000, nbits=0, compression=0)
    io = open(filename, "w")
    finalizer(io, close)
    const result = wavwrite(samples, io, Fs=Fs, nbits=nbits, compression=compression)
    close(io)
    return result
end

wavwrite(y::Array, f::Real, filename::String) = wavwrite(y, filename, Fs=f)
wavwrite(y::Array, f::Real, N::Real, filename::String) = wavwrite(y, filename, Fs=f, nbits=N)

# support for writing native arrays...
wavwrite{T<:Integer}(y::Array{T}, io::IO) = wavwrite(y, io, nbits=sizeof(T)*8)
wavwrite{T<:Integer}(y::Array{T}, filename::String) = wavwrite(y, filename, nbits=sizeof(T)*8)
wavwrite(y::Array{Int32}, io::IO) = wavwrite(y, io, nbits=24)
wavwrite(y::Array{Int32}, filename::String) = wavwrite(y, filename, nbits=24)
wavwrite{T<:FloatingPoint}(y::Array{T}, io::IO) = wavwrite(y, io, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)
wavwrite{T<:FloatingPoint}(y::Array{T}, filename::String) = wavwrite(y, filename, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)

end # module


text:
"""
    DataFunction

A type wrapping a data generating function.
"""
immutable DataFunction{T<:Number,D<:Number,A<:AbstractArray} <: AbstractData
    index::Vector
    values::A
    modifyvalues::Bool
    generator::Function
    args::Tuple
    DataFunction(i::AbstractVector{T},d::AbstractArray{D},w::Bool,g::Function,args...) = new(collect(i),d,w,g,args)
end

"""
    data(:function!,index::AbstractVector,out::AbstractArray,f!::Function,[,...args])

Create a `DataFunction` object which encapsulates an in-place data generating function.

The `:function!` variant expects a pre-allocated output array `out`,
and a data generating function which takes `out` as its first argument.

# Arguments
* `index::AbstractVector`: an index for the rows of the result (e.g., a time index)
* `out::AbstractArray`: array with one row per data point and one column per variable
* `f!::Function`: a data generating function taking as first argument `out`
* `args...`: additional arguments of the data generating function
---
"""
function data(::Type{Val{:function!}},i::AbstractVector,d::AbstractArray,g::Function,args...)
    @assert length(i) == size(d,1)
    g(d,args...)
    DataFunction{eltype(i),eltype(d),typeof(d)}(i,d,true,g,args...)
end

"""
    data(:function,index::AbstractVector,f::Function[,args...])

Create a `DataFunction` object which encapsulates a data generating function.

The `:function` variant can be used if the data generating function
cannot generate its result in a pre-allocated output array.
The result should contain one row per data point and one column per variable.

# Arguments
* `index::AbstractVector`: an index for the rows of the result (e.g., a time index)
* `f::Function`: a function returning an array of values as a result
* `args...`: additional arguments of the data generating function
---
"""
function data(::Type{Val{:function}},i::AbstractVector,g::Function,args...)
    d = g(args...)
    @assert length(i) == size(d,1)
    DataFunction{eltype(i),eltype(d),typeof(d)}(i,d,false,g,args...)
end

@inline _generate(::Type{Val{true}},d::DataFunction) = d.generator(d.values,d.args...)
@inline _generate(::Type{Val{false}},d::DataFunction) = copy!(d.values,d.generator(d.args...))

numvalues(d::DataFunction) = size(d.values,1)
numvars(d::DataFunction) = size(d.values,2)
eltype(d::DataFunction) = eltype(d.values)
generate!(d::DataFunction) = (_generate(Val{d.modifyvalues},d) ; d)
dataindex(d::DataFunction) = d.index
datavalues(d::DataFunction) = d.values
datatypename(d::DataFunction) = "DataFunction"


text:
using LoggingProfiler
using Test

@testset "LoggingProfiler.jl" begin
	function foo(x, y)
	  z =  x * y
	  z + sin(y)
	end

	@testset "whitelisting" begin
		ci = @code_lowered foo(1.0, 1.0)
		@test map(LoggingProfiler.timable, ci.code) == Bool[0, 0, 1, 1, 0]
		whitelist!(LoggingProfiler.timable_list, :sin)
		@test map(LoggingProfiler.timable, ci.code) == Bool[0, 0, 1, 0, 0]
		whitelist!(LoggingProfiler.timable_list, Core.GlobalRef(Main, :+))
		@test map(LoggingProfiler.timable, ci.code) == Bool[0, 0, 1, 1, 0]
	end

	LoggingProfiler.clear!()
	@record foo(1.0, 1.0)
	calls = LoggingProfiler.to[1];
	@test  calls.i[] == 364

	@test LoggingProfiler.iscallalist(calls, 1) == false
	@test LoggingProfiler.iscallalist(calls, 2) == true
	@test LoggingProfiler.iscallalist(calls, 3) == false

	LoggingProfiler.clear!()
	@record foo(1.0, 1.0)
	events = LoggingProfiler.tape2structure()
	LoggingProfiler._visualize("/tmp/test.html", events[1])
end

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler