text:
using Base.Cartesian
using BitBasis
using Combinatorics
#using Debugger
using SparseArrays
using Arpack

#spacial_dims=Int8(2)
#num_spin=2;
#num_species=2;
#L_x=2;
#L_y=1;
#num_sites = L_x*L_y
#num_n=num_spin*num_species*num_sites;
#N=2^num_n;
#t=1.0;
#U=1.0;
#mu=1.0;

## TO DO - assign types to all variables and functions
function occupation_ed(specs,nums,temp::Float64)

    N=length(nums)
    occ=zeros(Float64,length(nums[1]))
    z=0.0
    for i in 1:N
        num=nums[i]
        bolzman=sum(exp.(-(1/temp).*specs[i]))
        occ.+=bolzman.*num
        z+=bolzman
    end
    return occ./z
end

function spectrum(spacial_dims::NTuple{DIMS,Int64},num_species::Int64, t::Float64, U::Float64,μ::Float64) where {DIMS}
    num_spin=2
    N_max=prod(spacial_dims)
    nt=0:N_max
    n=[]
    for i in 1:num_spin*num_species
        push!(n,nt)
    end
    en_spec=[]
    n_vecs=collect(Iterators.product(n...))[:]
    for n_vec in n_vecs
        ham=hamiltonian_sub(spacial_dims,n_vec,t,U,μ)
        E=eigen(ham)
        push!(en_spec,E.values)
    end
    return en_spec , n_vecs
end

function hamiltonian_sub(spacial_dims::NTuple{DIMS,Int64},Num, t::Float64, U::Float64,μ::Float64) where {DIMS}

   states =  build_states(spacial_dims,Num)
   N=length(states)
   state_num=Int64[]
   for c in 1:N#state in states
    state_c = states[c]
        push!(state_num,packbits(state_c[:]))
   end

   H_sub=zeros(Float64,N,N)#spzeros(Float64,N,N)
   for i in 1:N
    state_i=states[i]
        state_sp , co = hamiltonian_operator(state_i,t,U,μ)
        for j in 1:length(state_sp)
            state_j=state_sp[j]
            index_j = findall(x->x==packbits(state_j[:]),state_num)
            H_sub[i,index_j...]+=co[j]
        end
    end
    return H_sub
end

function build_states(dims::NTuple{DIMS,Int64},Num::NTuple{DIMS2,Int64}) where {DIMS} where {DIMS2}
    #Num is a vector of  number operators for the different types of electrons (classified by their specie and spin)
    #dims is a vector of the spacial dimensions.
    # returen all states in the sub-space of the Fock space
    num_sites=prod(dims)
    states=[] # TO DO -set type to integer
    for Nsp in Num
        st = [ones(Int8,Nsp);zeros(Int8,num_sites-Nsp)]
        bit_combs=unique(permutations(st) )
        num_comb = [] # TO DO -set type to integer
        for bit_comb in bit_combs
            push!(num_comb,packbits(bit_comb)  )
        end
        push!(states,num_comb)
    end
    states_comb = collect(Iterators.product(states... ))[:]
    states_bin=[]
    for sc in states_comb
        push!(states_bin,bin_states(sc,dims,Int8(length(Num)/2)))
    end
       return states_bin
end

function bin_states(state::NTuple{DIMS2,Int64},dims::NTuple{DIMS,Int64},Nspecies::Int8) where {DIMS} where {DIMS2}
    Nspin=2
     bs=zeros(Int8,dims...,Nspecies,Nspin)
     for ns in 1:Nspecies
         for np in 1:Nspin
            bs[ntuple(k->:,length(dims))...,ns,np]=reshape(bitarray(state[ns+(np-1)*Nspecies],prod(dims)),dims)
         end
    end
    return bs
end

function hamiltonian_operator(state::Array{Int8,DIMS2} ,  t::Float64,U::Float64,μ::Float64) where {DIMS2}

    dims = size(state)
    spacial_dims = length(dims)-2
    state_sp = Array{Int8,length(dims)}[]
    co=Float64[]

    diag_term = U*sum( sum(state[ntuple(k->:,spacial_dims+1)...,1].-state[ntuple(k->:,spacial_dims+1)...,2] ,dims=spacial_dims+1 ).^2 ) #interaction term
    diag_term = -diag_term - μ*sum(state) #chemical potential term
    push!(state_sp,state)
    push!(co,diag_term)

    #hopping term
    occupations = findall(x->x==1, state) # indices which are  occupied
    for occupation in occupations
        for dir in 1:length(spacial_dims)
            for lr in 1:2
                index_i = occupation
                index_j = hop(index_i,dir,lr,dims)
                state_f , co_temp = hopping_operatr(state,index_i,index_j)
                if co_temp != 0
                    push!(state_sp,state_f)
                    push!(co,-t*co_temp)
                end
            end
        end
    end
    return state_sp , co
end

function create(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    # creation operator at (x,y,species,spin)
    new_state = copy(state)
    ex = 0  # denotes if the state was annihilated by the creation operator.
    if state[index] == 0
        new_state[index] = 1
        ex = check_fermion_comm(state ,index)#1
    end
        return new_state , ex
end

function annihilate(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    # annihilation operator at (x,y,spin,species)
    new_state = copy(state)
    ex = 0  # denotes if the state was annihilated by the annihilation operator.
    if state[index] == 1
        new_state[index] = 0
        ex = check_fermion_comm(state ,index)#1
    end
        return new_state , ex
end

function check_fermion_comm(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    state_l=LinearIndices(state)
    index_l=state_l[index]
    sign = sum(state[:][index_l+1:end])
    ex = (-1)^(sign-1)#(-1)^(sign-1)
    return ex
end

function hopping_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #hopping from  (x_i,y_i,species_i,spin_i) to  (x_j,y_j,species_f,spin_j)
    state_temp , ex1 = annihilate(state,index_i)
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = create(state_temp, index_j)
    return state_f , ex1*ex2
end

function hop(index::CartesianIndex{DIMS},dir::Int64,lr::Int64,dims::NTuple{DIMS,Int64}) where {DIMS}
    # update index
    if (lr==1)
      hop_index= index[dir]==dims[dir] ? 1 : index[dir]+1
    else
      hop_index= index[dir]==1 ? dims[dir] : index[dir]-1
    end
    # generate a new CartesianIndex with updated index
    CartesianIndex(Base.setindex(Tuple(index), hop_index, dir))

end

function build_states(dims::NTuple{DIMS,Int64}) where {DIMS}#(N,L_x,L_y,num_species,num_spin)
    # returen all states in the Fock space,  each state is reshaped to a an array of dims=(L_x,L_y,num_spin,num_species)
    num_n=prod(dims)
    N=2^num_n
    states=[];
    for c in 0:N-1
        state = dec2bin(c,num_n)
        push!(states,reshape(state , dims...))
    end
    return states
end


text:
# ---
# title: 423. Reconstruct Original Digits from English
# id: problem423
# author: Indigo
# date: 2021-06-20
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/reconstruct-original-digits-from-english/description/>
# hidden: true
# ---
#
# Given a **non-empty** string containing an out-of-order English representation
# of digits `0-9`, output the digits in ascending order.
#
# **Note:**
#
#   1. Input contains only lowercase English letters.
#   2. Input is guaranteed to be valid and can be transformed to its original digits. That means invalid inputs such as "abc" or "zerone" are not permitted.
#   3. Input length is less than 50,000.
#
# **Example 1:**
#
#
#
#     Input: "owoztneoer"
#
#     Output: "012"
#
#
# **Example 2:**
#
#
#
#     Input: "fviefuro"
#
#     Output: "45"
#
#
#
## @lc code=start
using LeetCode

function original_digits(s::String)
    cnt = counter(s)
    n0 = cnt['z']
    n2 = cnt['w']
    n8 = cnt['g']
    n6 = cnt['x']
    n3 = cnt['t'] - n2 - n8
    n4 = cnt['r'] - n3 - n0
    n7 = cnt['s'] - n6
    n1 = cnt['o'] - n4 - n2 - n0
    n5 = cnt['v'] - n7
    n9 = cnt['i'] - n8 - n6 - n5
    ns = (n0, n1, n2, n3, n4, n5, n6, n7, n8, n9)
    return join(string(i - 1)^n for (i, n) in enumerate(ns))
end
## @lc code=end

ext: jl


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

type PropertyKind
    value::lib.chfl_property_kind

"""
The possible types of Properties are:

- ``Chemfiles.PROPERTY_BOOL`` for storing bools
- ``Chemfiles.PROPERTY_DOUBLE`` for storing doubles
- ``Chemfiles.PROPERTY_STRING`` for storing strings
- ``Chemfiles.PROPERTY_VECTOR3D`` for storing vectors
"""
    function PropertyKind(value)
        value = lib.chfl_property_kind(value)
        if value in [lib.CHFL_PROPERTY_BOOL, lib.CHFL_PROPERTY_DOUBLE, lib.CHFL_PROPERTY_STRING, lib.CHFL_PROPERTY_VECTOR3D]
            return new(value)
        else
            throw(ChemfilesError("Invalid value for conversion to PropertyType: $value"))
        end
    end
end

Base.:(==)(x::PropertyKind, y::PropertyKind) = x.value == y.value

const PROPERTY_BOOL = PropertyKind(lib.CHFL_PROPERTY_BOOL)
const PROPERTY_DOUBLE = PropertyKind(lib.CHFL_PROPERTY_DOUBLE)
const PROPERTY_STRING = PropertyKind(lib.CHFL_PROPERTY_STRING)
const PROPERTY_VECTOR3D = PropertyKind(lib.CHFL_PROPERTY_VECTOR3D)

"""
Create a ``Bool`` ``Property``.
"""
function Property(value::Bool)
    return Property(lib.chfl_property_bool(convert(UInt8,value)))
end

"""
Create a ``Float64`` ``Property``.
"""
function Property(value::Float64)
    return Property(lib.chfl_property_double(value))
end

"""
Create a ``String`` ``Property``.
"""
function Property(value::String)
    return Property(lib.chfl_property_string(pointer(value)))
end

"""
Create a ``Vector`` ``Property``.
"""
function Property(value::Vector{Float64})
    return Property(lib.chfl_property_vector3d(value))
end

"""
Obtain the kind of property.
"""
function kind(property::Property)
    result = Ref{lib.chfl_property_kind}(0)
    _check(
        lib.chfl_property_get_kind(property.handle, result)
    )
    return PropertyKind(result[])
end

"""
Obtain the value stored by a property.
"""
function Base.get(property::Property)
    property_kind = kind(property)
    if property_kind == PROPERTY_BOOL
        result = Ref{UInt8}(0)
        _check(
            lib.chfl_property_get_bool(property.handle, result)
        )
        return convert(Bool, result[])
    elseif property_kind == PROPERTY_DOUBLE
        result = Ref{Cdouble}(0)
        _check(
            lib.chfl_property_get_double(property.handle, result)
        )
        return result[]
    elseif property_kind == PROPERTY_STRING
        return _call_with_growing_buffer(
            (buffer, size) -> _check(lib.chfl_property_get_string(property.handle, buffer, size))
        )
    elseif property_kind == PROPERTY_VECTOR3D
        result = Float64[0, 0, 0]
        _check(
            lib.chfl_property_get_vector3d(property.handle, result)
        )
        return result
    else
        throw( ChemfilesError("Invalid kind of property property $property_kind") )
    end

end

"""
Free the allocated memory for the ``Property`` object.
"""
function _free(property::Property)
    lib.chfl_property_free(property.handle)
end


text:
#
# widgets.jl -
#
# Implement Tk (and TTk) widgets
#

"""
    @TkWidget cls cmd pfx

create widget class `cls` based on Tk command `cmd` and using prefix `pfx` for
automatically defined widget names.  If `pfx` starts with a dot, a toplevel
widget class is assumed.  For now, `cmd` and `pfx` must be string literals.

"""
macro TkWidget(_cls, cmd, pfx)

    cls = esc(_cls)
    isa(cmd, String) || error("command must be a string literal")
    isa(pfx, String) || error("prefix must be a string literal")

    pfx[1] == '.' ? quote

        struct $cls <: TkRootWidget
            interp::TclInterp
            path::String
            obj::TclObj{$cls}
            function $cls(interp::TclInterp, name::Name=autoname($pfx);
                          kwds...)
                # It is sufficient to ensure that Tk is loaded for root widgets
                # because other widgets must have a parent.
                tkstart(interp)
                path = __widgetpath(__string(name))
                obj = __createwidget($cls, interp, $cmd, path; kwds...)
                new(interp, path, obj)
            end
        end

        $cls(name::Name = autoname($pfx); kwds...) =
            $cls(getinterp(), name; kwds...)

        (w::$cls)(args...; kwds...) = Tcl.eval(w, args...; kwds...)

    end : quote

        struct $cls <: TkWidget
            parent::TkWidget
            interp::TclInterp
            path::String
            obj::TclObj{$cls}
            function $cls(parent::TkWidget, child::Name=autoname($pfx);
                          kwds...)
                interp = getinterp(parent)
                path = __widgetpath(parent, __string(child))
                obj = __createwidget($cls, interp, $cmd, path; kwds...)
                new(parent, interp, path, obj)
            end
        end

        (w::$cls)(args...; kwds...) = Tcl.eval(w, args...; kwds...)

    end
end

@TkWidget TkToplevel      "::toplevel"          ".top"
@TkWidget TkMenu          "::menu"              ".mnu"

@TkWidget TtkButton       "::ttk::button"       "btn"
@TkWidget TtkCheckbutton  "::ttk::checkbutton"  "cbt"
@TkWidget TtkCombobox     "::ttk::combobox"     "cbx"
@TkWidget TtkEntry        "::ttk::entry"        "ent"
@TkWidget TtkFrame        "::ttk::frame"        "frm"
@TkWidget TtkLabel        "::ttk::label"        "lab"
@TkWidget TtkLabelframe   "::ttk::labelframe"   "lfr"
@TkWidget TtkMenubutton   "::ttk::menubutton"   "mbt"
@TkWidget TtkNotebook     "::ttk::notebook"     "nbk"
@TkWidget TtkPanedwindow  "::ttk::panedwindow"  "pwn"
@TkWidget TtkProgressbar  "::ttk::progressbar"  "pgb"
@TkWidget TtkRadiobutton  "::ttk::radiobutton"  "rbt"
@TkWidget TtkScale        "::ttk::scale"        "scl"
@TkWidget TtkScrollbar    "::ttk::scrollbar"    "sbr"
@TkWidget TtkSeparator    "::ttk::separator"    "sep"
@TkWidget TtkSizegrip     "::ttk::sizegrip"     "szg"
@TkWidget TtkSpinbox      "::ttk::spinbox"      "sbx"
@TkWidget TtkTreeview     "::ttk::treeview"     "trv"
@TkWidget TkButton        "::button"            "btn"
@TkWidget TkCanvas        "::canvas"            "cnv"
@TkWidget TkCheckbutton   "::checkbutton"       "cbt"
@TkWidget TkEntry         "::entry"             "ent"
@TkWidget TkFrame         "::frame"             "frm"
@TkWidget TkLabel         "::label"             "lab"
@TkWidget TkLabelframe    "::labelframe"        "lfr"
@TkWidget TkListbox       "::listbox"           "lbx"
@TkWidget TkMenubutton    "::menubutton"        "mbt"
@TkWidget TkMessage       "::message"           "msg"
@TkWidget TkPanedwindow   "::panedwindow"       "pwn"
@TkWidget TkRadiobutton   "::radiobutton"       "rbt"
@TkWidget TkScale         "::scale"             "scl"
@TkWidget TkScrollbar     "::scrollbar"         "sbr"
@TkWidget TkSpinbox       "::spinbox"           "sbx"
@TkWidget TkText          "::text"              "txt"

# Private method called to check/build the path of a child widget.
function __widgetpath(parent::TkWidget, child::String) :: String
    if findfirst(isequal('.'), child) === nothing
        Tcl.error("illegal window name \"$(child)\"")
    end
    parentpath = getpath(parent)
    return (parentpath == "." ? "." : parentpath*".")*child
end

# Private method called to check the path of a root widget.
function __widgetpath(path::String) :: String
    if path[1] != '.'
        Tcl.error("root window name must start with a dot")
    end
    if findnext(isequal('.'), child, 2) === nothing
        Tcl.error("illegal root window name \"$(path)\"")
    end
    return path
end

# Private method called to create a widget.
function __createwidget(::Type{T}, interp::TclInterp,
                        cmd::String, path::String;
                        kwds...) where {T<:Union{TkWidget,TkRootWidget}}
    local objptr :: TclObjPtr
    if interp(Bool, "winfo exists", path)
        # If widget already exists, it will be simply re-used, so we just apply
        # configuration options if any.  FIXME: there must be a way to check
        # the correctness of the widget class.
        if length(kwds) > 0
            exec(path, "configure"; kwds...)
        end
        objptr = __newobj(path)
    else
        # Widget does not already exists, create it with configuration options.
        if exec(TclStatus, interp, cmd, path; kwds...) != TCL_OK
            Tcl.error(interp)
        end
        objptr = Tcl_GetObjResult(interp.ptr)
    end
    return TclObj{T}(objptr)
end

"""
    TkToplevel([interp], ".")

yields the toplevel Tk window for interpreter `interp` (or for the initial
interpreter if this argument is omitted).  This also takes care of loading Tk
extension in the interpreter and starting the event loop.

To create a new toplevel window:

    TkToplevel([interp,] path; kwds...)

""" TkToplevel

getinterp(w::TkWidget) = w.interp
getpath(w::TkWidget) = w.path
getparent(w::TkWidget) = w.parent
getparent(::TkRootWidget) = nothing
TclObj(w::TkWidget) = w.obj
AtomicType(::Type{<:TkWidget}) = Atomic()
__objptr(w::TkWidget) = __objptr(w.obj)

getpath(root::TkWidget, args::AbstractString...) =
    getpath(getpath(root), args...)

getpath(arg0::AbstractString, args::AbstractString...) =
   join(((arg0 == "." ? "" : arg0), args...), '.')

Tcl.eval(w::TkWidget, args...) =
    Tcl.eval(getinterp(w), w.obj, args...)

exec(w::TkWidget, args...; kwds...) =
    exec(getinterp(w), w.obj, args...; kwds...)

"""
If Tk package is not yet loaded in interpreter `interp` (or in the initial
interpreter if this argument is missing), then:

```julia
tkstart([interp]) -> interp
```

will load Tk package and start the event loop.  The returned value is the
interpreter into which Tk has been started.  Note that this method also takes
care of withdrawing the root window "." to avoid its destruction as this would
terminate the Tcl application.

"""
function tkstart(interp::TclInterp = getinterp()) :: TclInterp
    if ! interp(Bool, "info exists tk_version")
        local status::TclStatus
        status = interp(TclStatus, "package require Tk")
        status == TCL_OK && (status = interp(TclStatus, "package require Ttk"))
        status == TCL_OK && (status = interp(TclStatus, "wm withdraw ."))
        status == TCL_OK || Tcl.error(interp)
        resume()
    end
    return interp
end

"""
    Tcl.configure(w)

yields all the options of Tk widget `w`, while:

    Tcl.configure(w, opt1=val1, opt2=val2)

change some options of widget `w`.  Options names may be specified as `String`
or `Symbol`.  Another way to change the settings is:

    w[opt1] = val1
    w[opt2] = val2

"""
configure(w::TkWidget; kwds...) = exec(w, "configure"; kwds...)

"""

    Tcl.cget(w, opt)

yields the value of the option `opt` for widget `w`.  Option `opt` may be
specified as a `String` or as a `Symbol`.  Another way to obtain an option
value is:

    w[opt]

"""
cget(w::TkWidget, opt::Name) = exec(w, "cget", "-"*string(opt))

Base.getindex(w::TkWidget, key::Name) = cget(w, key)
Base.setindex!(w::TkWidget, value, key::Name) = begin
    exec(w, "configure", "-"*string(key), value)
    w
end

"""
    Tcl.grid(args...; kwds...)
    Tcl.pack(args...; kwds...)
    Tcl.place(args...; kwds...)

communicate with one of the Tk geometry manager.  One of the arguments must be
an instance of `TkWidget`.  For instance (assuming `top` is some frame or
toplevel widget):

    using Tcl
    Tcl.pack(TkButton(top, "b1"; text="Send message",
                      command = (args...) -> println("message sent!")),
             side = "bottom")

"""
function grid end
@doc @doc(grid) pack
@doc @doc(grid) place

for cmd in (:grid, :pack, :place)
    @eval begin
        function $cmd(args...; kwds...)
            for arg in args
                if isa(arg, TkWidget)
                    interp = getinterp(arg)
                    return exec(interp, $(string(cmd)), args...; kwds...)
                end
            end
            Tcl.error("missing a widget argument")
        end
    end
end

# Base.bind is overloaded because it already exists for sockets, but there
# should be no conflicts.
"""
    bind(w, ...)

binds events to widget `w` or yields bindings for widget `w`.  With a single
argument

    bind(w)

yields binded sequences for widget `w`; while

    bind(w, seq)

yields the specific bindings for the sequence of events `seq` and

    bind(w, seq, script)

arranges to invoke `script` whenever any event of the sequence `seq` occurs for
widget `w`.  For instance:

    bind(w, "<ButtonPress>", "+puts click")

To deal with class bindings, the Tcl interpreter may be provided (otherwise the
initial interpreter will be used):

    bind([interp,] classname, args...)

where `classname` is the name of the widget class (a string or a symbol).

"""
Base.bind(arg0::TkWidget, args...) = bind(getinterp(arg0), arg0, args...)
Base.bind(arg0::Name, args...) = bind(getinterp(), arg0, args...)
Base.bind(interp::TclInterp, args...) = exec(interp, "bind", args...)


text:
x = range(0, stop=1, length=100)
y = x.^2
fill_between(x, y, x; color="orchid", alpha=0.2)
fill_between!(x, 0.5, 0.6; from=0.4, to=0.9, alpha=0.2)
style = (color="darkslategray", ls="--", lw=0.01)
hline(0.5; style...)
hline(0.6; style...)
vline(0.4; style...)
vline(0.9; style...)


text:
println("---")
using FactCheck, CompressedVolumes, FunctionalData

shouldtest(f, a) = length(ARGS) == 0 || in(a, ARGS) ? facts(f, a) : nothing
shouldtestcontext(f, a) = length(ARGS) < 2 || a == ARGS[2] ? facts(f, a) : nothing

shouldtest("overall") do
end

siz(a) = reshape([size(a)...], (3,1))

shouldtest("init") do

    # shouldtestcontext("fromarray") do
		# v = zeros(512,512,10)
		# for i in 1:size(v,3) v[:,:,i] = CompressedVolumes.rawimread("/Users/rene/temp/in/$i.png") end
		# v ./= maximum(v);
		# @time a = encode!(CompressedVolume(v, qp = 0))
		# @time r = decode(a)
		# @fact size(r) --> size(v)
		# @fact r[1:10] --> v[1:10]
		# @fact mean(r) --> mean(v)
    # end

	shouldtestcontext("fromsize") do
		a = CompressedVolume((5,4,3), qp = 0)
		c = decode(encode!(a))
		@fact siz(c) --> a.siz
		@fact c --> zeros(5,4,3)

		b = encode((5,4,3), qp = 0, defaultvalue = 1)
   		d = decode(b)
		@fact siz(d) --> b.siz
		@fact d --> ones(5,4,3)
	end
end

shouldtest("getindex") do
	m = 5; n = 4; o = 3
	v = zeros(m,n,o)
	v[:] = 1:length(v)
	# @show v
	a = encode(v, qp = 0)
	b = a[1:m,1:n,1:o]
	@fact b --> v

	a2 = encode(v, qp = 0, blocksiz = [1;1;1])
	b2 = a2[1:m,1:n,1:o]
	@fact b2 --> v

	a3 = encode(v, qp = 0, blocksiz = [2;3;2])
	b3 = a3[1:m,1:n,1:o]
	@fact b3 --> v
end

shouldtest("setindex!") do
	m = 5; n = 4; o = 3
	v = zeros(m,n,o)
	v[:] = 1:length(v)
	# @show v
	a = CompressedVolume(v)
	a[1,1,1] = 10.
	@fact a[1,1,1] --> 10.
	@fact a[1:2,1,1] --> reshape([10., 2.], 2,1,1)

	a = CompressedVolume((m,n,o))
	a[1,1,1] = 10.
	@fact a[1,1,1] --> 10.
	@fact a[1:2,1,1] --> reshape([10., 0.], 2,1,1)

	a[1,1:4,1] = 10.
	@fact a[1,3,1] --> 10.
	@fact a[1:2,2,1] --> reshape([10., 0.], 2,1,1)
end

shouldtest("compact!") do
	v = rand(1:100,3,4,2)
	a = encode(v, qp = 0)
	b = decode(a)
	compact!(a)
	c = decode(a)
	@fact maximum(abs(c-v))<0.5 --> true
end


text:
using CrystalNets
using Test, Random
using PeriodicGraphs
using StaticArrays
import Base.Threads

const cifs = joinpath(@__DIR__, "cif")
const crystalnetsdir = normpath(@__DIR__, "..")

@testset "Archive" begin
    Threads.@threads for (genome, id) in collect(CrystalNets.CRYSTAL_NETS_ARCHIVE)
        (id == "sxt" || id == "llw-z") && continue # special case for these known unstable nets
        @test reckognize_topology(topological_genome(PeriodicGraph(genome))) == id
    end
end

@testset "Module" begin
    Threads.@threads for target in
        ["pcu", "afy, AFY", "apc, APC", "bam", "bcf", "cdp", "cnd", "ecb", "fiv",
         "ftd", "ftj", "ins", "kgt", "mot", "moz", "muh", "pbz", "qom", "sig",
         "sma", "sod-f", "sod-h", "utj", "utp"]
       @info "Testing $target"
       graph = PeriodicGraph(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE[target])
       n = PeriodicGraphs.nv(graph)
       for k in 1:50
           r = randperm(n)
           offsets = [SVector{3,Int}([rand(-3:3) for _ in 1:3]) for _ in 1:n]
           graph = swap_axes!(offset_representatives!(graph[r], offsets), randperm(3))
           @test reckognize_topology(topological_genome(graph)) == target
       end
   end

   @test_broken reckognize_topology(topological_genome(PeriodicGraph(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE["sxt"]))) == "sxt"
   @test_broken reckognize_topology(topological_genome(PeriodicGraph(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE["llw-z"]))) == "llw-z"

   @info """The following warnings about guessing bonds are expected."""
   @test reckognize_topology(topological_genome(CrystalNet(parse_chemfile(joinpath(cifs, "Moganite.cif"))))) == "mog"
end

@testset "Executable" begin
    safestdout = deepcopy(stdout)
    safeARGS = deepcopy(ARGS)
    piperead, _ = redirect_stdout()
    empty!(ARGS)

    push!(ARGS, "-g", "3   1 2  0 0 0   1 2  0 0 1   1 2  0 1 0   1 2  1 0 0")
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "dia"

    empty!(ARGS)
    path = joinpath(cifs, "ABW.cif")
    push!(ARGS, path)
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "sra, ABW"

    empty!(ARGS)
    path = joinpath(cifs, "ABW.cif")
    push!(ARGS, "-a", CrystalNets.arc_location*"rcsr.arc", path)
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "sra"

    empty!(ARGS)
    path = joinpath(cifs, "RRO.cif")
    push!(ARGS, path)
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "RRO"

    empty!(ARGS)
    path = joinpath(cifs, "RRO.cif")
    push!(ARGS, "-a", CrystalNets.arc_location*"rcsr.arc", path)
    @test CrystalNets.julia_main() == 1
    @test readuntil(piperead, '\n') == "UNKNOWN"

    empty!(ARGS)
    path = joinpath(cifs, "HKUST-1.cif")
    push!(ARGS, "-c", "mof", path)
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "tbo"

    empty!(ARGS)
    path = joinpath(cifs, "Diamond.cif")
    push!(ARGS, "-c", "atom", path)
    @test CrystalNets.julia_main() == 0
    @test readuntil(piperead, '\n') == "dia"

    empty!(ARGS)
    push!(ARGS, "--help")
    @test CrystalNets.julia_main() == 0
    @test startswith(readuntil(piperead, '\n'), "usage: CrystalNets")
    @test occursin("CRYSTAL_FILE", readuntil(piperead, '\n'))
    @test occursin("Form B", readuntil(piperead, '\n'))
    @test occursin("Form C", readuntil(piperead, '\n'))
    @test isempty(readuntil(piperead, '\n'))
    @test isempty(readuntil(piperead, '\n'))
    @test readuntil(piperead, '\n') == "Automatic reckognition of crystal net topologies."

    redirect_stdout(safestdout)
    empty!(ARGS)
    append!(ARGS, safeARGS)
end


text:
# Probability computations

"""
	$(SIGNATURES)

Validate a probability matrix. Elements are expected to sum to 1.
"""
function validate_prob_matrix(probM :: AbstractMatrix{T1}) where T1 <: AbstractFloat
    isValid = all_at_least(probM, 0.0)  &&  all_at_most(probM, 1.0);
    isValid = isValid  &&  isapprox(sum(probM), 1.0);
    return isValid
end

"""
	$(SIGNATURES)

Validate a probability vector. Elements are expected to sum to 1.
"""
validate_prob_vector(probV :: AbstractVector{T1}) where T1 <: AbstractFloat =
    validate_prob_matrix(probV');

"""
	$(SIGNATURES)

Conditional probabilities of rows or columns given a probability matrix.
"""
prob_j(m :: AbstractMatrix{T1}) where T1 <: AbstractFloat = vec(sum(m, dims = 2));
prob_j(m :: AbstractMatrix{T1}, j :: Integer) where T1 <: AbstractFloat =
    sum(m[j,:]);
prob_j(m :: AbstractMatrix{T1}, j) where T1 <: AbstractFloat =
    sum(m[j, :], dims = 2);

"""
	$(SIGNATURES)

Conditional probabilities of rows or columns given a probability matrix.
"""
prob_k(m :: AbstractMatrix{T1}) where T1 <: AbstractFloat = vec(sum(m, dims = 1));
prob_k(m :: AbstractMatrix{T1}, k :: Integer) where T1 <: AbstractFloat =
    sum(m[:, k]);
prob_k(m :: AbstractMatrix{T1}, k) where T1 <: AbstractFloat =
    vec(sum(m[:, k], dims = 1));

"""
    $(SIGNATURES)

Conditional probability (j | k).
"""
prob_j_k(m :: AbstractMatrix, j, k :: Integer) =
    m[j, k] ./ prob_k(m, k);

"""
    $(SIGNATURES)

Conditional probability (k | j).
"""
prob_k_j(m :: AbstractMatrix, k, j :: Integer) =
    m[j, k] ./ prob_j(m, j);

## -----------  Expected values

"""
	$(SIGNATURES)

Compute expected value of `x`, given row index `j`. Based on matrix of values `xM[j,k]` and matrix of probabilities (sum to 1) `prM[j, k]`.
"""
ev_given_j(x :: AbstractMatrix, prM :: AbstractMatrix, j :: Integer) =
    sum([(prob_k_j(prM, k, j) * x[j, k])  for k = 1 : size(x, 2) ]);
ev_given_j(x :: AbstractMatrix, prM :: AbstractMatrix) =
    [ev_given_j(x, prM, j)  for j = 1 : size(x, 1)];

"""
	$(SIGNATURES)

Compute expected value of `x`, given column index `k`. Based on matrix of values `xM[j,k]` and matrix of probabilities (sum to 1) `prM[j, k]`.
"""
ev_given_k(x :: AbstractMatrix, prM :: AbstractMatrix, k :: Integer) =
    sum([(prob_j_k(prM, j, k) * x[j, k])  for j = 1 : size(x, 1) ]);
ev_given_k(x :: AbstractMatrix, prM :: AbstractMatrix) =
    [ev_given_k(x, prM, k)  for k = 1 : size(x, 2)];

    """
    $(SIGNATURES)

Given an array of probabilities: ensure that all are in [0, 1].
Error if bounds violation larger than rounding errors.
Make sure that sum does not exceed an upper bound.
"""
function scale_prob_array!(m :: AbstractArray{F1}; maxSum :: F1 = one(F1)) where F1 <: AbstractFloat

    fSmall = F1(.0000001);
    pSum = sum(m);
    @assert (pSum < one(F1) + fSmall)  "Sum too large: $pSum"
    @assert all_at_least(m, -fSmall)  "Negative probabilities"
    @assert all_at_most(m, one(F1) + fSmall)  "Probabilities above 1"

    if any(x -> x < zero(F1), m)
        m[findall(x -> x < zero(F1), m)] .= zero(F1);
    end
    if pSum > maxSum
        m .*= ((maxSum - fSmall) / pSum);
    end
end

# ---------


text:
y= [2, 5, 12]

x = [12, 10]
println(y)
println(x)

println("Results without code")
println(x)

y = 1:5
println(y)

a = "Don't print me"
println(a)

println("No markup for results.")

println(collect(0:10:1000))

println(collect(0:10:1000))

println(collect(0:10:1000))


text:
#
#  Split - Apply - Combine operations
#

##############################################################################
##
## GroupedDataTable...
##
##############################################################################

"""
The result of a `groupby` operation on an AbstractDataTable; a
view into the AbstractDataTable grouped by rows.

Not meant to be constructed directly, see `groupby`.
"""
type GroupedDataTable
    parent::AbstractDataTable
    cols::Vector         # columns used for sorting
    idx::Vector{Int}     # indexing vector when sorted by the given columns
    starts::Vector{Int}  # starts of groups
    ends::Vector{Int}    # ends of groups
end

#
# Split
#

"""
A view of an AbstractDataTable split into row groups

```julia
groupby(d::AbstractDataTable, cols)
groupby(cols)
```

### Arguments

* `d` : an AbstractDataTable to split (optional, see [Returns](#returns))
* `cols` : data table columns to group by

### Returns

* `::GroupedDataTable` : a grouped view into `d`
* `::Function`: a function `x -> groupby(x, cols)` (if `d` is not specified)

### Details

An iterator over a `GroupedDataTable` returns a `SubDataTable` view
for each grouping into `d`. A `GroupedDataTable` also supports
indexing by groups and `map`.

See the following for additional split-apply-combine operations:

* `by` : split-apply-combine using functions
* `aggregate` : split-apply-combine; applies functions in the form of a cross product
* `combine` : combine (obviously)
* `colwise` : apply a function to each column in an AbstractDataTable or GroupedDataTable

Piping methods `|>` are also provided.

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
gd = groupby(dt, :a)
gd[1]
last(gd)
vcat([g[:b] for g in gd]...)
for g in gd
    println(g)
end
map(d -> mean(dropnull(d[:c])), gd)   # returns a GroupApplied object
combine(map(d -> mean(dropnull(d[:c])), gd))
dt |> groupby(:a) |> [sum, length]
dt |> groupby([:a, :b]) |> [sum, length]
```

"""
function groupby{T}(dt::AbstractDataTable, cols::Vector{T}; sort::Bool = false)
    sdt = dt[cols]
    dt_groups = group_rows(sdt)
    # sort the groups
    if sort
        group_perm = sortperm(view(sdt, dt_groups.rperm[dt_groups.starts]))
        permute!(dt_groups.starts, group_perm)
        Base.permute!!(dt_groups.stops, group_perm)
    end
    GroupedDataTable(dt, cols, dt_groups.rperm,
                     dt_groups.starts, dt_groups.stops)
end
groupby(d::AbstractDataTable, cols; sort::Bool = false) = groupby(d, [cols], sort = sort)

# add a function curry
groupby{T}(cols::Vector{T}; sort::Bool = false) = x -> groupby(x, cols, sort = sort)
groupby(cols; sort::Bool = false) = x -> groupby(x, cols, sort = sort)

Base.start(gd::GroupedDataTable) = 1
Base.next(gd::GroupedDataTable, state::Int) =
    (view(gd.parent, gd.idx[gd.starts[state]:gd.ends[state]]),
     state + 1)
Base.done(gd::GroupedDataTable, state::Int) = state > length(gd.starts)
Base.length(gd::GroupedDataTable) = length(gd.starts)
Base.endof(gd::GroupedDataTable) = length(gd.starts)
Base.first(gd::GroupedDataTable) = gd[1]
Base.last(gd::GroupedDataTable) = gd[end]

Base.getindex(gd::GroupedDataTable, idx::Int) =
    view(gd.parent, gd.idx[gd.starts[idx]:gd.ends[idx]])
Base.getindex(gd::GroupedDataTable, I::AbstractArray{Bool}) =
    GroupedDataTable(gd.parent, gd.cols, gd.idx, gd.starts[I], gd.ends[I])

Base.names(gd::GroupedDataTable) = names(gd.parent)
_names(gd::GroupedDataTable) = _names(gd.parent)

##############################################################################
##
## GroupApplied...
##    the result of a split-apply operation
##    TODOs:
##      - better name?
##      - ref
##      - keys, vals
##      - length
##      - start, next, done -- should this return (k,v) or just v?
##      - make it a real associative type? Is there a need to look up key columns?
##
##############################################################################

"""
The result of a `map` operation on a GroupedDataTable; mainly for use
with `combine`

Not meant to be constructed directly, see `groupby` abnd
`combine`. Minimal support is provided for this type. `map` is
provided for a GroupApplied object.

"""
immutable GroupApplied{T<:AbstractDataTable}
    gd::GroupedDataTable
    vals::Vector{T}

    @compat function (::Type{GroupApplied})(gd::GroupedDataTable, vals::Vector)
        length(gd) == length(vals) ||
            throw(DimensionMismatch("GroupApplied requires keys and vals be of equal length (got $(length(gd)) and $(length(vals)))."))
        new{eltype(vals)}(gd, vals)
    end
end

#
# Apply / map
#

# map() sweeps along groups
function Base.map(f::Function, gd::GroupedDataTable)
    GroupApplied(gd, [wrap(f(dt)) for dt in gd])
end
function Base.map(f::Function, ga::GroupApplied)
    GroupApplied(ga.gd, [wrap(f(dt)) for dt in ga.vals])
end

wrap(dt::AbstractDataTable) = dt
wrap(A::Matrix) = convert(DataTable, A)
wrap(s::Any) = DataTable(x1 = s)

"""
Combine a GroupApplied object (rudimentary)

```julia
combine(ga::GroupApplied)
```

### Arguments

* `ga` : a GroupApplied

### Returns

* `::DataTable`

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
combine(map(d -> mean(dropnull(d[:c])), gd))
```

"""
function combine(ga::GroupApplied)
    gd, vals = ga.gd, ga.vals
    valscat = vcat(vals...)
    idx = Vector{Int}(size(valscat, 1))
    j = 0
    @inbounds for (start, val) in zip(gd.starts, vals)
        n = size(val, 1)
        idx[j + (1:n)] = gd.idx[start]
        j += n
    end
    hcat!(gd.parent[idx, gd.cols], valscat)
end

"""
Apply a function to each column in an AbstractDataTable or
GroupedDataTable

```julia
colwise(f::Function, d)
colwise(d)
```

### Arguments

* `f` : a function or vector of functions
* `d` : an AbstractDataTable of GroupedDataTable

If `d` is not provided, a curried version of groupby is given.

### Returns

* various, depending on the call

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
colwise(sum, dt)
colwise([sum, lenth], dt)
colwise((minimum, maximum), dt)
colwise(sum, groupby(dt, :a))
```

"""
function colwise(f, d::AbstractDataTable)
    x = [f(d[i]) for i in 1:ncol(d)]
    if eltype(x) <: Nullable
        return NullableArray(x)
    else
        return x
    end
end
# apply several functions to each column in a DataTable
function colwise(fns::Union{AbstractVector, Tuple}, d::AbstractDataTable)
    x = [f(d[i]) for f in fns, i in 1:ncol(d)]
    if eltype(x) <: Nullable
        return NullableArray(x)
    else
        return x
    end
end
colwise(f, gd::GroupedDataTable) = [colwise(f, g) for g in gd]
colwise(f) = x -> colwise(f, x)

"""
Split-apply-combine in one step; apply `f` to each grouping in `d`
based on columns `col`

```julia
by(d::AbstractDataTable, cols, f::Function; sort::Bool = false)
by(f::Function, d::AbstractDataTable, cols; sort::Bool = false)
```

### Arguments

* `d` : an AbstractDataTable
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `f` : a function to be applied to groups; expects each argument to
  be an AbstractDataTable
* `sort`: sort row groups (no sorting by default)

`f` can return a value, a vector, or a DataTable. For a value or
vector, these are merged into a column along with the `cols` keys. For
a DataTable, `cols` are combined along columns with the resulting
DataTable. Returning a DataTable is the clearest because it allows
column labeling.

A method is defined with `f` as the first argument, so do-block
notation can be used.

`by(d, cols, f)` is equivalent to `combine(map(f, groupby(d, cols)))`.

### Returns

* `::DataTable`

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
by(dt, :a, d -> sum(d[:c]))
by(dt, :a, d -> 2 * dropnull(d[:c]))
by(dt, :a, d -> DataTable(c_sum = sum(d[:c]), c_mean = mean(dropnull(d[:c]))))
by(dt, :a, d -> DataTable(c = d[:c], c_mean = mean(dropnull(d[:c]))))
by(dt, [:a, :b]) do d
    DataTable(m = mean(dropnull(d[:c])), v = var(dropnull(d[:c])))
end
```

"""
by(d::AbstractDataTable, cols, f::Function; sort::Bool = false) =
    combine(map(f, groupby(d, cols, sort = sort)))
by(f::Function, d::AbstractDataTable, cols; sort::Bool = false) =
    by(d, cols, f, sort = sort)

#
# Aggregate convenience functions
#

# Applies a set of functions over a DataTable, in the from of a cross-product
"""
Split-apply-combine that applies a set of functions over columns of an
AbstractDataTable or GroupedDataTable

```julia
aggregate(d::AbstractDataTable, cols, fs)
aggregate(gd::GroupedDataTable, fs)
```

### Arguments

* `d` : an AbstractDataTable
* `gd` : a GroupedDataTable
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `fs` : a function or vector of functions to be applied to vectors
  within groups; expects each argument to be a column vector

Each `fs` should return a value or vector. All returns must be the
same length.

### Returns

* `::DataTable`

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
aggregate(dt, :a, sum)
aggregate(dt, :a, [sum, x->mean(dropnull(x))])
aggregate(groupby(dt, :a), [sum, x->mean(dropnull(x))])
dt |> groupby(:a) |> [sum, x->mean(dropnull(x))]   # equivalent
```

"""
aggregate(d::AbstractDataTable, fs::Function; sort::Bool=false) = aggregate(d, [fs], sort=sort)
function aggregate{T<:Function}(d::AbstractDataTable, fs::Vector{T}; sort::Bool=false)
    headers = _makeheaders(fs, _names(d))
    _aggregate(d, fs, headers, sort)
end

# Applies aggregate to non-key cols of each SubDataTable of a GroupedDataTable
aggregate(gd::GroupedDataTable, f::Function; sort::Bool=false) = aggregate(gd, [f], sort=sort)
function aggregate{T<:Function}(gd::GroupedDataTable, fs::Vector{T}; sort::Bool=false)
    headers = _makeheaders(fs, setdiff(_names(gd), gd.cols))
    res = combine(map(x -> _aggregate(without(x, gd.cols), fs, headers), gd))
    sort && sort!(res, cols=headers)
    res
end

(|>)(gd::GroupedDataTable, fs::Function) = aggregate(gd, fs)
(|>){T<:Function}(gd::GroupedDataTable, fs::Vector{T}) = aggregate(gd, fs)

# Groups DataTable by cols before applying aggregate
function aggregate{S<:ColumnIndex, T <:Function}(d::AbstractDataTable,
                                                 cols::Union{S, AbstractVector{S}},
                                                 fs::Union{T, Vector{T}};
                                                 sort::Bool=false)
    aggregate(groupby(d, cols, sort=sort), fs)
end

function _makeheaders{T<:Function}(fs::Vector{T}, cn::Vector{Symbol})
    fnames = _fnames(fs) # see other/utils.jl
    [Symbol(colname,'_',fname) for fname in fnames for colname in cn]
end

function _aggregate{T<:Function}(d::AbstractDataTable, fs::Vector{T}, headers::Vector{Symbol}, sort::Bool=false)
    res = DataTable(Any[vcat(f(d[i])) for f in fs for i in 1:size(d, 2)], headers)
    sort && sort!(res, cols=headers)
    res
end


text:
using LinearAlgebra
sq(x) = x*x

"Fits a straight line through a set of points, `y = a₁ + a₂ * x`"
function linear_fit(x, y)

    sx = sum(x)
    sy = sum(y)

    m = length(x)

    sx2 = zero(sx)
    sy2 = zero(sy)
    sxy = zero(sx*sy)

    for i = 1:m
        sx2 += x[i]*x[i]
        sy2 += y[i]*y[i]
        sxy += x[i]*y[i]
    end

    a0 = (sx2*sy - sxy*sx) / ( m*sx2 - sx*sx )
    a1 = (m*sxy - sx*sy) / (m*sx2 - sx*sx)

    return [a0, a1]
end

"Fits a log function through a set of points: `y = a₁+ a₂*log(x)`"
log_fit(x, y) = linear_fit(log.(x), y)

"Fits a power law through a set of points: `y = a₁*x^a₂`"
function power_fit(x, y)
    fit = linear_fit(log.(x), log.(y))
    [exp(fit[1]), fit[2]]
end

"Fits an `exp` through a set of points: `y = a₁*exp(a₂*x)`"
function exp_fit(x, y)
    fit = linear_fit(x, log.(y))
    [exp(fit[1]), fit[2]]
end

"""
Create Vandermonde matrix for simple polynomial fit
"""
function vandermondepoly(x, y, n)
    m = length(x)
    A = zeros(eltype(y), m, n+1)
    A[:,1] .= 1.0

    for i = 1:n
        for k in 1:m
            A[k, i+1] = A[k, i] * x[k]
        end
    end
    return A
end

"""
    Calculates the coefficients of a linear model using Least Squares

Given a Vandermonde matrix, this function uses the QR factorization to compute the Least Squares
fit of a linear model. The code probably could be optimized.

"""
function fit_linear_model(Av, y)
    m, n = size(Av)
    F = qr(Av)

    d = (F.Q * Matrix(I, m, n))' * y
    R = UpperTriangular(F.R)

    coefs = R\d

    return coefs

end

"""
Fits a polynomial of degree `n` through a set of points.

Simple algorithm, doesn't use orthogonal polynomials or any such thing
and therefore unconditioned matrices are possible. Use it only for low
degree polynomials.

This function returns a the coefficients of the polynomial.
"""
function poly_fit(x, y, n)

    A = vandermondepoly(x, y, n)
    coefs = fit_linear_model(A, y)
    return coefs
end

"""
High Level interface for fitting straight lines
"""
struct LinearFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    LinearFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
LinearFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = LinearFit{T}(linear_fit(x, y))

"High Level interface for fitting log laws"
struct LogFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    LogFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
LogFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = LogFit{T}(log_fit(x, y))

"High Level interface for fitting power laws"
struct PowerFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    PowerFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
PowerFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = PowerFit{T}(power_fit(x, y))

"High Level interface for fitting exp laws"
struct ExpFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    ExpFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
ExpFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = ExpFit{T}(exp_fit(x, y))

"""
# Generic interface for curve fitting.

The same function `curve_fit` can be used to fit the data depending on fit type,
shich is specified in the first parameter. This function returns an object that
can be used to estimate the value of the fitting model using function `apply_fit`.

## A few examples:

 * `f = curve_fit(LinearFit, x, y)`
 * `f = curve_fit(Poly, x, y, n)`
"""
curve_fit(::Type{T}, x, y) where {T<:LeastSquares} = T(x, y)
curve_fit(::Type{T}, x, y, args...) where {T<:LeastSquares} = T(x, y, args...)
curve_fit(::Type{Poly}, x, y, n=1) = Poly(poly_fit(x, y, n))

(f::LinearFit)(x) = f.coefs[1] + f.coefs[2] *x
(f::PowerFit)(x) = f.coefs[1] * x ^ f.coefs[2]
(f::LogFit)(x) = f.coefs[1] + f.coefs[2] * log(x)
(f::ExpFit)(x) = f.coefs[1] * exp(f.coefs[2] * x)


text:
using CuArrays, Test

include("util.jl")

using Random
Random.seed!(1)

using CUDAnative
using CUDAdrv

# GPUArrays has a testsuite that isn't part of the main package.
# Include it directly.
import GPUArrays
gpuarrays = pathof(GPUArrays)
gpuarrays_root = dirname(dirname(gpuarrays))
include(joinpath(gpuarrays_root, "test", "testsuite.jl"))

testf(f, xs...; kwargs...) = TestSuite.compare(f, CuArray, xs...; kwargs...)

import CuArrays: allowscalar, @allowscalar
allowscalar(false)

@testset "CuArrays" begin

@testset "GPUArrays test suite" begin
  TestSuite.test(CuArray)
end

include("base.jl")
include("blas.jl")
include("rand.jl")
include("fft.jl")
include("sparse.jl")
include("solver.jl")
include("sparse_solver.jl")
include("dnn.jl")
include("tensor.jl")
include("forwarddiff.jl")

CuArrays.memory_status()
CuArrays.pool_timings()
CuArrays.alloc_timings()
CuArrays.reset_timers!()

end


text:
module MesaX
include("structs.jl")
include("parser.jl")
import Random, JSON

greet() = print("Hello World!")
greet_alien() = print("Hello ", Random.randstring(8))

catfood = ["foo", "bar", "baz"]

moose = (cat = "33", dog = 33)

export greet, greet_alien

##
#     q = Meta.parse("""function foo(i::Array{Int64})::Int64
#     sum::Int64 = 0
#     for x in i
#         sum += x
#     end
#     sum
# end """)
#     dogs = Meta.parse("""ba = rand(Int64, 50000000)""")
#     z = Meta.parse("""foo(ba)""")
#     zp = Meta.parse("""sum(ba)""")
#     Base.eval(q)
#     Base.eval(dogs)
#     @time Base.eval(z)
#     @time Base.eval(zp)
##

end # module

max_stars_repo_name: mesa-x/spreadsheet

max_issues_repo_name: mesa-x/spreadsheet

max_forks_repo_name: mesa-x/spreadsheet


text:
__precompile__(false)
module MicroBench

# The main reason for this module's existence
using CachedArrays

# stdlib
using Dates
using LinearAlgebra
using Random
using Statistics

# "internal" dependencies
using MaxLFSR
using KernelBenchmarks  # for fast initialization of arrays.
using SystemSnoop
using CounterTools
using MattDaemon
using ExperimentsDB

# external deps
using DataStructures
using ProgressMeter
using StructArrays

# debug mode enables much smaller kernels that don't exceed the size of the Local
const DEBUG = false

function __init__()
    if DEBUG
        @warn """
        Running in debug mode!!

        Make sure you aren't trying to record data!!
        """
    end
end

# Will throw an error if this is not defined
const IS_2LM = parse(Bool, ENV["JULIA_IS_2LM"])
const IS_1LM = !IS_2LM

# More Prefetch definitions
maybesuper(x...) = CachedArrays.maybesuper(x...)

function LinearAlgebra.mul!(Y::AbstractCachedArray, A::AbstractCachedArray, B::AbstractCachedArray)
    CachedArrays.shallowfetch!(Y)
    invoke(
        LinearAlgebra.mul!,
        Tuple{maybesuper(Y), maybesuper(A), maybesuper(B)},
        CachedArrays.unlock(Y),
        A,
        B,
    )
    return nothing
end

include("counters.jl")
include("1d/top_1d.jl")
include("2d/top_2d.jl")

end # module


text:
# ---
# title: 937. Reorder Data in Log Files
# id: problem937
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/reorder-data-in-log-files/description/>
# hidden: true
# ---
#
# You have an array of `logs`.  Each log is a space delimited string of words.
#
# For each log, the first word in each log is an alphanumeric _identifier_.
# Then, either:
#
#   * Each word after the identifier will consist only of lowercase letters, or;
#   * Each word after the identifier will consist only of digits.
#
# We will call these two varieties of logs _letter-logs_ and _digit-logs_.  It
# is guaranteed that each log has at least one word after its identifier.
#
# Reorder the logs so that all of the letter-logs come before any digit-log.
# The letter-logs are ordered lexicographically ignoring identifier, with the
# identifier used in case of ties.  The digit-logs should be put in their
# original order.
#
# Return the final order of the logs.
#
#
#
# **Example 1:**
#
#
#
#     Input: logs = ["dig1 8 1 5 1","let1 art can","dig2 3 6","let2 own kit dig","let3 art zero"]
#     Output: ["let1 art can","let3 art zero","let2 own kit dig","dig1 8 1 5 1","dig2 3 6"]
#
#
#
#
# **Constraints:**
#
#   1. `0 <= logs.length <= 100`
#   2. `3 <= logs[i].length <= 100`
#   3. `logs[i]` is guaranteed to have an identifier, and a word after the identifier.
#
#
## @lc code=start
using LeetCode

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


text:
module TestLift

using Maybe
using Test

include("utils.jl")

@testset "identity" begin
    @test (@? nothing) === Some(nothing)
    @test (@? identity(nothing)) === nothing
    @test (@? Some(nothing)) === Some(nothing)
    @test (@? Some(Some(nothing))) === Some(nothing)
    @test (@? Some(Some(Some(nothing)))) === Some(nothing)
    @test (@? Some(Some(Some(nothing)))) === Some(nothing)
    @test (@? $Some($Some($Some(nothing)))) === Some(Some(Some(Some(nothing))))
end

@testset "return returns nothing" begin
    @? function f1()
        return
    end
    f2 = @? function ()
        return
    end
    f3 = @? () -> return
    @testset for f in map(something, [f1, f2, f3])
        @test f() === nothing
    end
end

@testset "return nothing returns something" begin
    @? function f1()
        return nothing
    end
    @? function f2()
    end
    f3 = @? function ()
        return nothing
    end
    f4 = @? function ()
    end
    @? f5() = nothing
    @? f6() = begin end
    f7 = @? () -> nothing
    f8 = @? () -> begin end
    @testset for f in map(something, [f1, f2, f3, f4, f5, f6, f7, f8])
        @test f() === Some(nothing)
    end
end
# Since `fᵢ` and `fᵢ₊₁` (`i` is odd) above should work identically,
# `return nothing` should be lifted to `return Some(nothing)`.

@testset "\$return(x)" begin
    f1() = @? $return(1)
    @? function f2()
        $return(1)
    end
    f3 = @? function ()
        $return(1)
    end
    @? f4() = $return(1)
    f5 = @? () -> $return(1)
    @testset for f in map(something, [f1, f2, f3, f4, f5])
        @test f() === 1
    end
end

@testset "keyword arguments" begin
    f(a; b = 1, c = 2) = (a, b, c)
    f() = nothing
    @test (@? f(1; b = f())) === nothing
    @test (@? f(1, b = f())) === nothing
    @test (@? f(b = f(), 1)) === nothing
    @test (@? f(1; b = f(1; b = nothing)[2])) === Some((1, nothing, 2))
    @test (@? f(1, b = f(1, b = nothing)[2])) === Some((1, nothing, 2))
    @test (@? f(b = f(b = nothing, 1)[2], 1)) === Some((1, nothing, 2))
end

@testset "splatting (positional arguments)" begin
    t(args...) = args
    f() = nothing
    g() = (1, 2)
    h() = Some((1, 2))
    @test (@? t(0, f()...)) === nothing
    @test (@? t(0, g()...)) === Some((0, 1, 2))
    @test (@? t(0, h()...)) === Some((0, 1, 2))
    x = 0
    @test (@? t(x, f()...)) === nothing
    @test (@? t(x, g()...)) === Some((0, 1, 2))
    @test (@? t(x, h()...)) === Some((0, 1, 2))
    @test (@? t(f()..., x)) === nothing
    @test (@? t(g()..., x)) === Some((1, 2, 0))
    @test (@? t(h()..., x)) === Some((1, 2, 0))
    xname = :x
    @test (@? t(xname => x, f()...)) === nothing
    @test (@? t(xname => x, g()...)) === Some((xname => 0, 1, 2))
    @test (@? t(xname => x, h()...)) === Some((xname => 0, 1, 2))
end

@testset "splatting (keyword arguments)" begin
    t(; kwargs...) = values(kwargs)
    f() = nothing
    g() = (a = 1, b = 2)
    h() = Some((a = 1, b = 2))
    @test (@? t(; x = 0, f()...)) === nothing
    @test (@? t(; x = 0, g()...)) === Some((x = 0, a = 1, b = 2))
    @test (@? t(; x = 0, h()...)) === Some((x = 0, a = 1, b = 2))
    x = 0
    @test (@? t(; x, y = 1, f()...)) === nothing
    @test (@? t(; x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? t(; x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    xname = :x
    @test (@? t(; xname => x, y = 1, f()...)) === nothing
    @test (@? t(; xname => x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? t(; xname => x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
end

@testset "splatting (mixed)" begin
    t(args...; kwargs...) = (args, values(kwargs))
    f() = nothing
    gt() = (1, 2)
    ht() = Some((1, 2))
    gk() = (a = 1, b = 2)
    hk() = Some((a = 1, b = 2))
    @test (@? t(x = 0, f()...)) === nothing
    @test (@? t(x = 0, gk()...)) === Some(t(x = 0, gk()...))
    @test (@? t(x = 0, hk()...)) === Some(t(x = 0, gk()...))
    @test (@? t(0, f()...; f()...)) === nothing
    @test (@? t(0, gt()...; gk()...)) === Some(t(0, gt()...; gk()...))
    @test (@? t(0, ht()...; hk()...)) === Some(t(0, gt()...; gk()...))
    x = 0
    @test (@? t(x, y = 1, f()...)) === nothing
    @test (@? t(x, y = 1, gk()...)) === Some(t(x, y = 1, gk()...))
    @test (@? t(x, y = 1, hk()...)) === Some(t(x, y = 1, gk()...))
    xname = :x
    @test (@? t(xname => x, y = 1, f()...)) === nothing
    @test (@? t(xname => x, y = 1, gk()...)) === Some(t(xname => x, y = 1, gk()...))
    @test (@? t(xname => x, y = 1, hk()...)) === Some(t(xname => x, y = 1, gk()...))
end

@testset "splatting (tuple)" begin
    f() = nothing
    g() = (1, 2)
    h() = Some((1, 2))
    @test (@? (0, f()...)) === nothing
    @test (@? (0, g()...)) === Some((0, 1, 2))
    @test (@? (0, h()...)) === Some((0, 1, 2))
    @test (@? tuple(0, f()...)) === nothing
    @test (@? tuple(0, g()...)) === Some((0, 1, 2))
    @test (@? tuple(0, h()...)) === Some((0, 1, 2))
    x = 0
    @test (@? (x, f()...)) === nothing
    @test (@? (x, g()...)) === Some((0, 1, 2))
    @test (@? (x, h()...)) === Some((0, 1, 2))
    @test (@? tuple(x, f()...)) === nothing
    @test (@? tuple(x, g()...)) === Some((0, 1, 2))
    @test (@? tuple(x, h()...)) === Some((0, 1, 2))
    xname = :x
    @test (@? (xname => x, f()...)) === nothing
    @test (@? (xname => x, g()...)) === Some((:x => 0, 1, 2))
    @test (@? (xname => x, h()...)) === Some((:x => 0, 1, 2))
    @test (@? tuple(xname => x, f()...)) === nothing
    @test (@? tuple(xname => x, g()...)) === Some((:x => 0, 1, 2))
    @test (@? tuple(xname => x, h()...)) === Some((:x => 0, 1, 2))
end

@testset "splatting (named tuple)" begin
    f() = nothing
    g() = (a = 1, b = 2)
    h() = Some((a = 1, b = 2))
    @test (@? (x = 0, f()...)) === nothing
    @test (@? (x = 0, g()...)) === Some((x = 0, a = 1, b = 2))
    @test (@? (x = 0, h()...)) === Some((x = 0, a = 1, b = 2))
    @test (@? (; x = 0, f()...)) === nothing
    @test (@? (; x = 0, g()...)) === Some((x = 0, a = 1, b = 2))
    @test (@? (; x = 0, h()...)) === Some((x = 0, a = 1, b = 2))
    x = 0
    @test (@? (x, y = 1, f()...)) === nothing
    @test (@? (x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (; x, y = 1, f()...)) === nothing
    @test (@? (; x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (; x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    xname = :x
    @test (@? (xname => x, y = 1, f()...)) === nothing
    @test (@? (xname => x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (xname => x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (; xname => x, y = 1, f()...)) === nothing
    @test (@? (; xname => x, y = 1, g()...)) === Some((x = 0, y = 1, a = 1, b = 2))
    @test (@? (; xname => x, y = 1, h()...)) === Some((x = 0, y = 1, a = 1, b = 2))
end

struct Indexable end
Base.getindex(::Indexable, args...; kwargs...) = (args, values(kwargs))

@testset "Indexable" begin
    A = Indexable()
    ent = NamedTuple()
    @test A[1] === ((1,), ent)
    @test A[1, 2] === ((1, 2), ent)
    @test A[1, a = 2] === ((1,), (a = 2,))
end

@testset "getindex" begin
    A = Indexable()
    @test (@? A[1]) === Some(A[1])
    @test (@? A[1, 2]) === Some(A[1, 2])
    @test (@? A[1, a = 2]) === Some(A[1, a = 2])

    vov = [[11, 12], [21, 22, 23]]
    @test (@? vov[1][1]) === Some(vov[1][1])
    @test (@? vov[2][3]) === Some(vov[2][3])
    @test (@? vov[-1][1]) === nothing
    @test (@? vov[5][1]) === nothing
    @test (@? vov[1][-1]) === nothing
    @test (@? vov[1][5]) === nothing
end

@testset "getproperty" begin
    x = (a = (b = 1, c = nothing),)
    @test (@? x.a.b) === Some(1)
    @test (@? x.a.c) === Some(nothing)
    @test (@? x.a.c.d) === nothing
    @test (@? x.a.d) === nothing
    @test (@? x.d.e) === nothing
    nthg = nothing
    @test (@? nthg.a) === nothing
    n() = nothing
    @test (@? n().a) === nothing
end

@testset "destructuring bind" begin
    n() = nothing
    @test (@? begin
        a, b = n()
        a
    end) === nothing
    @test (@? begin
        a, b = (1, 2)
        a
    end) === Some(1)
end

@testset "control flows" begin
    t(_...) = true
    f(_...) = false
    n() = nothing
    @test (@? t() ? 1 : 2) === Some(1)
    @test (@? f() ? 1 : 2) === Some(2)
    @test (@? n() === nothing ? 1 : 2) === nothing
    @test (@? n() !== nothing ? 1 : 2) === nothing
    @test (@? $n() === nothing ? 1 : 2) === Some(1)
    @test (@? $n() !== nothing ? 1 : 2) === Some(2)
    @test (@? t() && t() ? 1 : 2) === Some(1)
    @test (@? t() && f() ? 1 : 2) === Some(2)
    @test (@? t() || t() ? 1 : 2) === Some(1)
    @test (@? t() || f() ? 1 : 2) === Some(1)
    @test (@? t() & t() ? 1 : 2) === Some(1)
    @test (@? t() & f() ? 1 : 2) === Some(2)
    @test (@? t() | t() ? 1 : 2) === Some(1)
    @test (@? t() | f() ? 1 : 2) === Some(1)
    @test (@? 1 < 2 < 3 ? 1 : 2) === Some(1)
    @test (@? 1 < 20 < 3 ? 1 : 2) === Some(2)
    @test (@? n() < 2 < 3 ? 1 : 2) === nothing
    @test (@? 1 < n() < 3 ? 1 : 2) === nothing
    @test (@? 1 < 2 < n() ? 1 : 2) === nothing
    @test (@? begin
        x = 2
        if t()
            x -= 1
        end
        x
    end) === Some(1)
    @test (@? begin
        x = 2
        if f()
            x -= 1
        end
        x
    end) === Some(2)
    @test (@? begin
        x = 2
        if t(n())
            x -= 1
        end
        x
    end) === nothing
    @test begin
        @? begin
            ys = []
            for x in 1:3
                push!(ys, x)
            end
            ys
        end
    end |> something |> ==(1:3)
    @test begin
        xs = [1, nothing, 2]
        @? begin
            ys = []
            for x in xs
                push!(ys, x)
            end
            ys
        end
    end === nothing
end

@testset "do" begin
    d = Dict(:a => 1)
    @test (@? get(d, :a) do; end) === Some(1)
    @test (@? get(d, :b) do; 2 end) === Some(2)
    @test (@? get(d, :b) do; return nothing end) === Some(nothing)
    @test (@? get(d, :b) do; end) === Some(nothing)
    @test (@? $get(d, :a) do; end) === Some(Some(1))
    @test (@? $get(d, :b) do; 2 end) === Some(Some(2))
    @test (@? $get(d, :b) do; return nothing end) === Some(Some(nothing))
    @test (@? $get(d, :b) do; end) === Some(Some(nothing))
end

@testset "unused" begin
    n() = nothing
    @test (@? begin
        n()
        1
    end) === nothing
end

@testset "literals" begin
    @test_broken (@? [1, nothing, 2]) |> something |> ==([1, Some(nothing), 2])
    @test_broken (@? [1; nothing; 2]) |> something |> ==([1; Some(nothing); 2])
    @test_broken (@? [1 2; nothing 4; 5 6]) |> something |> ==([1 2; Some(nothing) 4; 5 6])
    @test (@? [1, identity(nothing), 2]) === nothing
    @test (@? [1; identity(nothing); 2]) === nothing
    @test (@? [1 2; identity(nothing) 4; 5 6]) === nothing
end

@testset "type assertions" begin
    # Is this even possible with the current framework?
    d = Dict{Symbol,Any}(:a => Dict{Symbol,Any}(:b => 1))
    @test (@? begin
        a::Dict{Symbol,Any} = d[:a]
        b::Integer = a[:b]
        b + 1
    end) === Some(2)
end

@testset "dot-call" begin
    err = ArgumentError("dot-call is not supported yet")
    @test_throws err @callmacro @? a .+ b
    @test_throws err @callmacro @? f.(a, b)
    @test_throws err @callmacro @? g(f.(a, b))
    @test_throws err @callmacro @? f.(a, b) + c
    @test_throws err @callmacro @? f(a .+ b)
end

end  # module


text:
export DCTOp, \

mutable struct DCTOp{T,F1<:FuncOrNothing,F2<:FuncOrNothing,F3<:FuncOrNothing} <: AbstractLinearOperator{T,F1,F2,F3}
  nrow   :: Int
  ncol   :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod   :: F1 # apply the operator to a vector
  tprod  :: F2 # apply the transpose operator to a vector
  ctprod :: F3 # apply the transpose conjugate operator to a vector
end

"""
  DCTOp(T::Type, shape::Tuple, dcttype=2)

returns a `DCTOp <: AbstractLinearOperator` which performs a DCT on a given input array.

# Arguments:
* `T::Type`       - type of the array to transform
* `shape::Tuple`  - size of the array to transform
* `dcttype`       - type of DCT (currently `2` and `4` are supported)
"""
function DCTOp(T::Type, shape::Tuple, dcttype=2)
  if dcttype == 2
    return DCTOp{T,Function,Nothing,Function}(prod(shape), prod(shape), true, false
            , x->vec((dct(reshape(x,shape))))
            , nothing
            , y->vec((idct(reshape(y,shape)))))
  elseif dcttype == 4
    factor = sqrt(1.0/(prod(shape)* 2^length(shape)) )
    return DCTOp{T,Function,Function,Nothing}(prod(shape), prod(shape), true, false
            , x->vec((FFTW.r2r(reshape(x,shape),FFTW.REDFT11).*factor))
            , x->vec((FFTW.r2r(reshape(x,shape),FFTW.REDFT11).*factor))
            , nothing )
  else
    error("DCT type $(dcttype) not supported")
  end
end

\(A::DCTOp, x::Vector) = adjoint(A) * x


text:
function glplot(arg1, style = :default;
        screen = viewing_screen(), camera = :perspective,
        kw_args...
    )
    robj = visualize(arg1, style; kw_args...)
    _view(robj, screen, camera = camera)
    register_plot!(robj, screen)
    robj
end

function register_plot!(robj::Vector, screen = viewing_screen(); create_gizmo = false)
    vcat(map(robj) do elem
        register_plot!(elem, screen, create_gizmo = create_gizmo)
    end...)
end
function register_plot!(robj::Context, screen = viewing_screen(); create_gizmo = false)
    register_plot!(robj.children, screen, create_gizmo = create_gizmo)
end

function left_clicked(button::Set{Int})
    return GLAbstraction.singlepressed(button, GLFW.MOUSE_BUTTON_LEFT)
end

function to_clip(p, width, proj_view)
    clip = proj_view * Vec4f0(p, 1f0)
    clip = clip ./ clip[4]
    nv = (clip + 1f0) * 0.5f0
    Point2f0(nv[1], nv[2]) .* Point2f0(width-1)
end

function gizmo(w, is_selected)
    colors = Array(RGBA{Float32}, 6)
    positions = Array(Point3f0, 6)
    for axis in 1:3
        c = RGBA{Float32}(ntuple(i->i==axis ? 1f0 : 0f0, 3)...)
        v = unit(Vec3f0, axis)
        positions[2(axis-1)+1] = Point3f0(0)
        positions[2(axis-1)+2] = v
        colors[2(axis-1)+1] = RGBA{Float32}(0,0,0,1)
        colors[2(axis-1)+2] = c
    end
    robj = visualize(
        positions, :linesegment,
        color = colors,
        prerender = ()->glDisable(GL_DEPTH_TEST), # draw over other items
        postrender = ()->glEnable(GL_DEPTH_TEST),
        thickness = 10f0
    ).children[]
    _view(robj, camera = :perspective)
    @materialize mouseposition, mouse_buttons_pressed, buttons_pressed = w.inputs
    m2id = mouse2id(w)
    show_gizmo = map(buttons_pressed) do bp
        GLAbstraction.singlepressed(bp, GLFW.KEY_M) && value(is_selected)
    end
    set_arg!(robj, :visible, show_gizmo)
    start_drag = foldp((false, 0), m2id) do v0, id
        ison = id.id == robj.id && id.index > 0 && id.index <= 6 && value(show_gizmo)
        if ison
            return (true, id.index)
        else
            return (false, v0[2])
        end
    end
    left_pressed = const_lift(GLAbstraction.pressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    view = w.cameras[:perspective].projectionview
    ddiff = GLAbstraction.dragged_diff(mouseposition, left_pressed, map(first, start_drag))

    translation = foldp(Vec3f0(0), ddiff) do v0, diff
        idx = value(start_drag)[2]
        pv = value(view)
        if idx != 0
            i = ceil(Int, idx/2)
            uv = unit(Vec3f0, i)
            wh = widths(w)
            nv2d = to_clip(uv, wh, pv)
            origin = to_clip(Vec3f0(0), wh, pv)
            nv2d = nv2d-origin
            mv = -Vec2f0(diff).*10f0
            return v0 + (uv * (dot(mv, normalize(nv2d)) ./ (norm(wh))))
        end
        v0
    end
    model = map(translationmatrix, translation)
    set_arg!(robj, :model, model)
    model
end

"""
There's a lot of noise if we just go through all parameters of a `RenderObject`.
This function filters out the internal values
"""
function is_editable(k, v_v)
    v = value(v_v)
    !(
        k == :objectid ||
        k == :is_fully_opaque ||
        k == :instances ||
        k == Symbol("position.multiplicator") ||
        k == Symbol("position.dims") ||
        k == Symbol("spatialorder") ||
        k == Symbol("resolution") ||
        k == Symbol("fxaa") ||
        k == Symbol("light") ||
        k == Symbol("light") ||
        k == Symbol("doc_string") ||
        k == Symbol("faces") ||
        k == Symbol("image") ||
        k == Symbol("vertices") ||
        k == Symbol("texturecoordinates") ||
        k == Symbol("ranges") ||
        k == Symbol("model") ||
        k == :visible ||
        startswith(string(k), "boundingbox") ||
        (k == Symbol("color") && isa(v, AbstractArray)) ||
        k in fieldnames(PerspectiveCamera) ||
        k == :instances ||
        isa(v, Symbol) ||
        isa(v, Void) ||
        isa(v, NativeMesh) ||
        isa(v, Int) ||
        isa(v, Int32) ||
        isa(v, UInt32) ||
        isa(v, UInt) ||
        (isa(v, FixedVector) && eltype(v) <: Integer)
    )
end

function to_edit_dict(robj)
    filter(collect(robj.uniforms)) do kv
        is_editable(kv[1], kv[2])
    end
end

function register_plot!(
        robj::RenderObject, screen = viewing_screen();
        create_gizmo = false
    )
    left_gap = 1.5mm
    visible_button, visible_s = visible_toggle()
    set_arg!(robj, :visible, visible_s)
    delete_button, del_signal = button(imload("delete.png"), edit_screen())
    edit_button, no_edit_signal = edit_toggle()

    item_height = Signal(0)
    not_del_signal = droprepeats(foldp(false, del_signal) do deleted, to_delete
        deleted && return deleted
        if to_delete
            push!(item_height, 0)
            delete!(screen, robj)
        end
        return to_delete
    end)

    edit_signal = map(!, no_edit_signal)
    item_screen = widget_screen!(delete=not_del_signal)
    edititemarea = map(edit_item_area, item_screen.area, item_height, Signal(left_gap))
    edit_item_screen = Screen(edit_screen(), area = edititemarea)
    preserve(foldp((false, value(item_height)), edit_signal) do v0, edit
        if edit
            if !v0[1] # only do this at the first time
                dict = to_edit_dict(robj)
                nh = if !isempty(dict)
                    vis, signal_dict = extract_edit_menu(
                        dict,
                        edit_item_screen,
                        edit_signal,
                    )
                    for (k, v) in signal_dict
                        robj.uniforms[k] = v
                    end
                    _view(vis, edit_item_screen, camera = :fixed_pixel)
                    new_heights = widths(value(boundingbox(vis)))[2]
                    ceil(Int, new_heights)
                else
                    0
                end
                push!(item_height, nh)
                return true, nh
            else
                push!(item_height, v0[2])
            end
        else
            push!(item_height, 0)
        end
        return v0
    end)
    @materialize buttons_pressed, mouse_buttons_pressed, mouseinside = screen.inputs
    dc = doubleclick(mouse_buttons_pressed, 0.1)
    selected = foldp(false, mouse2id(screen), dc, mouseinside) do v0, id, dc, mi
        dc && mi && return id.id == robj.id # if double clicked and inside and hovers over robj
        return v0
    end
    color = map(selected) do s
        s ? RGBA{Float32}(0.9, 0.99, 1, 1) : RGBA{Float32}(1, 1, 1, 1)
    end
    select_icon = visualize(SimpleRectangle(0,0,1.5mm, 10mm), color=color)
    add_widget!(select_icon, visible_button, delete_button, edit_button, window=item_screen)
    if create_gizmo
        model = gizmo(screen, selected)
        set_arg!(robj, :model, model)
    end
    [del_signal]
end


text:
__precompile__()
module NamedTuples

export @NT, NamedTuple, setindex, delete

abstract type NamedTuple end

Base.keys( t::NamedTuple ) = fieldnames( t )
Base.values( t::NamedTuple ) = [ getfield( t, i ) for i in 1:nfields( t ) ]
Base.haskey( t::NamedTuple, k ) = k in keys(t)
Base.length( t::NamedTuple ) = nfields( t )
# Iteration
Base.start( t::NamedTuple ) = 1
Base.done( t::NamedTuple, iter ) = iter > nfields( t )
Base.next( t::NamedTuple, iter ) = ( getfield( t, iter ), iter + 1 )
Base.endof( t::NamedTuple ) = length( t )
Base.last( t::NamedTuple ) = t[end]
function Base.show( io::IO, t::NamedTuple )
    print(io, "(")
    first = true
    for (k,v) in zip(keys(t),values(t))
        !first && print(io, ", ")
        print(io, k, " = "); show(io, v)
        first = false
    end
    print(io, ")")
end
# Make this indexable so that it works like a Tuple
Base.getindex( t::NamedTuple, i::Int ) = getfield( t, i )
# We also support indexing by symbol
Base.getindex( t::NamedTuple, i::Symbol ) = getfield( t, i )
Base.getindex( t::NamedTuple, i::Symbol, default ) = get( t, i, default )
# This is a linear lookup...
Base.get( t::NamedTuple, i::Symbol, default ) = i in keys(t) ? t[i] : default
# Deep compare
import Base: ==

@generated function ==( lhs::NamedTuple, rhs::NamedTuple)
    if !isequal(fieldnames(lhs), fieldnames(rhs)) || lhs !== rhs
        return false
    end

    q = quote end

    for i in 1:length( fieldnames(lhs) )
        push!(q.args, :(lhs[$(i)] == rhs[$(i)] || return false))
    end

    return q
end

# Deep hash
@generated function Base.hash(nt::NamedTuple, hs::UInt64)
    q = quote
        h = 17
    end

    for i in 1:length(fieldnames(nt))
        push!(q.args, :(h = h * 23 + hash(nt[$(i)], hs)))
    end

    return q
end

# Helper type, for transforming parse tree to NameTuple definition
struct ParseNode{T} end

function trans( ::Type{ParseNode{:Symbol}}, expr::Expr)
    (expr.args[1],nothing,nothing)
end

function trans( ::Union{Type{ParseNode{:(=)}},Type{ParseNode{:kw}}}, expr::Expr)
    (sym, typ ) = trans( expr.args[1])
    return (sym, typ, expr.args[2] )
end

function trans( ::Type{ParseNode{:call}}, expr::Expr)
    if expr.args[1] == :(=>)
        Base.depwarn("\"=>\" syntax for NamedTuple construction is deprecated, use \"=\" instead.", Symbol("@NT"))
        (sym, typ ) = trans( expr.args[1])
        return (sym, typ, expr.args[2] )
    end
    return (nothing, nothing, expr)
end

# Allow unary
function trans( ::Type{ParseNode{:(::)}}, expr::Expr)
    if( length( expr.args ) > 1 )
        return ( expr.args[1], expr.args[2], nothing)
    else
        return ( nothing, expr.args[1], nothing)
    end
end

function trans( expr::Expr )
    trans( ParseNode{expr.head}, expr)
end

function trans( sym::Symbol )
    return (sym, nothing, nothing)
end

function trans( ::Type{ParseNode{:quote}}, expr::Expr )
    return trans( expr.args[1] )
end

# Literal nodes
function trans{T}( lit::T )
    return (nothing, nothing, lit)
end

function trans{T}( ::Type{ParseNode{T}}, expr::Expr)
    return (nothing, nothing, expr)
end

function gen_namedtuple_ctor_body(n::Int, args)
    types = [ :(typeof($x)) for x in args ]
    cnvt = [ :(convert(fieldtype(TT,$n),$(args[n]))) for n = 1:n ]
    if n == 0
        texpr = :T
    else
        texpr = :(NT{$(types...)})
    end
    if isless(Base.VERSION, v"0.6.0-")
        tcond = :(NT === NT.name.primary)
    else
        tcond = :(isa(NT,UnionAll))
    end
    quote
        if $tcond
            TT = $texpr
        else
            TT = NT
        end
        if nfields(TT) !== $n
            throw(ArgumentError("wrong number of arguments to named tuple constructor"))
        end
        $(Expr(:new, :TT, cnvt...))
    end
end

# constructor for all NamedTuples
@generated function (::Type{NT}){NT<:NamedTuple}(args...)
    n = length(args)
    aexprs = [ :(args[$i]) for i = 1:n ]
    return gen_namedtuple_ctor_body(n, aexprs)
end

# specialized for certain argument counts
for n = 0:5
    args = [ Symbol("x$n") for n = 1:n ]
    @eval function (::Type{NT}){NT<:NamedTuple}($(args...))
        $(gen_namedtuple_ctor_body(n, args))
    end
end

# Create a NameTuple type, if a type with these field names has not already
# been constructed.
# TODO: to make modules containing named tuples precompile-able, change `= NamedTuples` to `= current_module()`
function create_namedtuple_type(fields::Vector{Symbol}, mod::Module = NamedTuples)
    escaped_fieldnames = [replace(string(i), "_", "__") for i in fields]
    name = Symbol( string( "_NT_", join( escaped_fieldnames, "_")) )
    if !isdefined(mod, name)
        len = length( fields )
        types = [Symbol("T$n") for n in 1:len]
        tfields = [ Expr(:(::), Symbol( fields[n] ), Symbol( "T$n") ) for n in 1:len ]
        def = Expr(:type, false, Expr( :(<:), Expr( :curly, name, types... ), GlobalRef(NamedTuples, :NamedTuple) ),
                   Expr(:block, tfields...,
                        Expr(:tuple)))  # suppress default constructors
        eval(mod, def)
    end
    return getfield(mod, name)
end

#
# Given a symbol list create the NamedTuple
#
@doc doc"Given a symbol vector create the `NamedTuple`" ->
function make_tuple( syms::Vector{Symbol} )
    return create_namedtuple_type( syms )
end

#
# Given an expression vector create the NamedTuple
#
@doc doc"Given an expression vector create the `NamedTuple`" ->
function make_tuple( exprs::Vector)
    len    = length( exprs )
    fields = Array{Symbol}(len)
    values = Array{Any}(len)
    typs   = Array{Any}(len)

    # Are we directly constructing the type, if so all values must be
    # supplied by the caller, we use this state to ensure this
    construct = false
    # handle the case where this is defining a datatype
    for i in 1:len
        expr = exprs[i]
        ( sym, typ, val ) = trans( expr )
        if( construct == true && val == nothing || ( i > 1 && construct == false && val != nothing ))
            error( "Invalid tuple, all values must be specified during construction @ ($expr)")
        end
        construct  = val != nothing
        fields[i]  = sym != nothing?sym:Symbol( "_$(i)_")
        typs[i] = typ
        # On construction ensure that the types are consitent with the declared types, if applicable
        values[i]  = ( typ != nothing && construct)? Expr( :call, :convert, typ, val ) : val
    end

    ty = create_namedtuple_type( fields )

    # Either call the constructor with the supplied values or return the type
    if( !construct )
        if len == 0
            return ty
        end
        return Expr( :curly, ty, typs... )
    else
        return Expr( :call, ty, values ... )
    end
end

@doc doc"""
Syntax

    @NT( a, b )                 -> Defines a tuple with a and b as members
    @NT( a::Int64, b::Float64 ) -> Defines a tuple with the specific arg types as members
    @NT( a = 1, b = "hello")  -> Defines and constructs a tuple with the specifed members and values
    @NT( a, b )( 1, "hello")    -> Is equivalent to the above definition
    @NT( a::Int64 )( 2.0 )      -> Calls `convert( Int64, 2.0 )` on construction and sets `a`

NamedTuples may be used anywhere you would use a regular Tuple, this includes method definition and return values.

    module Test
    using NamedTuples

    function foo( y )
        a = 1
        x = 3
        return  @NT( a = 1, b = "world", c = "hello", d=a/x, y = a/y  )
    end
    function bar( nt::@NT( a::Int64, c::ASCIIString ))
        return repeat( nt.c, nt.a )
    end

    end

    Test.foo( 1 ) # Returns a NamedTuple of 5 elements
    Test.bar( @NT( a= 2, c="hello")) # Returns `hellohello`
""" ->
macro NT( expr... )
    return esc(make_tuple( collect( expr )))
end

# Helper function for 0.4 compat
if VERSION < v"0.5.0"
    getfieldname( t, i ) = fieldnames(t)[i]
else
    getfieldname( t, i ) = fieldname( t, i )
end

@inline function Base.map(f, nt::NamedTuple, nts::NamedTuple...)
    # this method makes sure we don't define a map(f) method
    _map(f, nt, nts...)
end

@generated function _map(f, nts::NamedTuple...)
    fields = fieldnames(nts[1])
    for x in nts[2:end]
        if !isequal(fieldnames(x), fields)
            throw(ArgumentError("All NamedTuple inputs to map must have the same fields in the same order"))
        end
    end
    N = nfields(nts[1])
    M = length(nts)

    # This type will already exist if this function may be called
    NT = create_namedtuple_type(fields, moduleof(nts[1]))
    args = Expr[:(f($(Expr[:(getfield(nts[$i], $j)) for i = 1:M]...))) for j = 1:N]
    quote
        $NT($(args...))
    end
end

@generated function Base.isless(t1::NamedTuple, t2::NamedTuple)
    if !isequal(fieldnames(t1), fieldnames(t2))
        throw(ArgumentError("NamedTuple inputs to isless must have the same fields in the same order"))
    end
    quote
        Base.@nexprs $(nfields(t1)) i -> begin
            a_i = getfield(t1, i)
            b_i = getfield(t2, i)
            if !isequal(a_i, b_i)
                return isless(a_i, b_i)
            end
        end
        return false
    end
end

function Base.getindex( t::NamedTuple, rng::AbstractVector )
    names = unique( Symbol[ isa(i,Symbol) ? i : getfieldname(typeof(t),i) for i in rng ] )
    ty = create_namedtuple_type( names )
    ty([ getfield( t, i ) for i in names ]...)
end

@doc doc"""
Merge two NamedTuples favoring the lhs
Order is preserved lhs names come first.
This copies the underlying data.
""" ->
function Base.merge( lhs::NamedTuple, rhs::NamedTuple )
    nms = unique( vcat( fieldnames( lhs ), fieldnames( rhs )) )
    ty = create_namedtuple_type( nms )
    # FIXME should handle the type only case
    vals = [ haskey( lhs, nm ) ? lhs[nm] : rhs[nm] for nm in nms ]
    ty(vals...)
end

@doc doc"""
Create a new NamedTuple with the new value set on it, either overwriting
the old value or appending a new value.
This copies the underlying data.
""" ->
function setindex{V}( t::NamedTuple, key::Symbol, val::V)
    nt = create_namedtuple_type( [key] )( val )
    return merge( t, nt )
end

@doc doc"""
Create a new NamedTuple with the specified element removed.
""" ->
function delete( t::NamedTuple, key::Symbol )
    nms = filter( x->x!=key, fieldnames( t ) )
    ty = create_namedtuple_type( nms )
    vals = [ getindex( t, nm ) for nm in nms ]
    return ty(vals...)
end

Base.Broadcast._containertype(::Type{<:NamedTuple}) = NamedTuple
Base.Broadcast.promote_containertype(::Type{NamedTuple}, ::Type{NamedTuple}) = NamedTuple
Base.Broadcast.promote_containertype(::Type{NamedTuple}, _) = error()
Base.Broadcast.promote_containertype(_, ::Type{NamedTuple}) = error()

@inline function Base.Broadcast.broadcast_c(f, ::Type{NamedTuple}, nts...)
    _map(f, nts...)
end

moduleof(t::DataType) = t.name.module

if VERSION < v"0.6.0-dev"
    uniontypes(u) = u.types
else
    const uniontypes = Base.uniontypes
    moduleof(t::UnionAll) = moduleof(Base.unwrap_unionall(t))
end

struct NTType end
struct NTVal end

function Base.serialize{NT<:NamedTuple}(io::AbstractSerializer, ::Type{NT})
    if NT === Union{}
        Base.Serializer.write_as_tag(io, Base.Serializer.BOTTOM_TAG)
    elseif isa(NT, Union)
        Base.serialize_type(io, NTType)
        serialize(io, Union)
        serialize(io, [uniontypes(NT)...])
    elseif isleaftype(NT)
        Base.serialize_type(io, NTType)
        serialize(io, fieldnames(NT))
        serialize(io, moduleof(NT))
        write(io.io, UInt8(0))
        serialize(io, NT.parameters)
    else
        u = Base.unwrap_unionall(NT)
        if isa(u, DataType) && NT === u.name.wrapper
            Base.serialize_type(io, NTType)
            serialize(io, fieldnames(NT))
            serialize(io, moduleof(NT))
            write(io.io, UInt8(1))
        else
            error("cannot serialize type $NT")
        end
    end
end

function Base.deserialize(io::AbstractSerializer, ::Type{NTType})
    fnames = deserialize(io)
    if fnames == Union
        types = deserialize(io)
        return Union{types...}
    else
        mod = deserialize(io)
        NT = create_namedtuple_type(fnames, mod)
        if read(io.io, UInt8) == 0
            params = deserialize(io)
            return NT{params...}
        else
            return NT
        end
    end
end

function Base.serialize(io::AbstractSerializer, x::NamedTuple)
    Base.serialize_type(io, NTVal)
    serialize(io, typeof(x))
    for i in 1:nfields(x)
        serialize(io, getfield(x, i))
    end
end

function Base.deserialize(io::AbstractSerializer, ::Type{NTVal})
    NT = deserialize(io)
    nf = nfields(NT)
    if nf == 0
        return NT()
    elseif nf == 1
        return NT(deserialize(io))
    elseif nf == 2
        return NT(deserialize(io), deserialize(io))
    elseif nf == 3
        return NT(deserialize(io), deserialize(io), deserialize(io))
    else
        return NT(Any[ deserialize(io) for i = 1:nf ]...)
    end
end

end # module


text:
using ElectromagneticFields
using Test

include("test_analytic.jl")


text:
using Printf
@userplot Specificationplot

"""
    specificationplot([S,CS,T], [WS,WU,WT])

This function visualizes the control synthesis using the hInf_synthesize with
the three weighting functions {WS(jω), WU(jω), WT(jω)} inverted and scaled by γ,
against the corresponding transfer fucntions {S(jω), C(jω)S(jω), T(jω)}, to
verify visually that the specifications are met. This may be run using both MIMO
and SISO systems.
"""
specificationplot
@recipe function specificationplot(
    p::Specificationplot;
    wint = (-3, 5),
    wnum = 101,
    hz = true,
    s_labels = [
        "σ(S(jω))",
        "σ(C(jω)S(jω))",
        "σ(T(jω))",
    ],
    w_labels = [
        "γ σ(Wₛ(jω)⁻¹)",
        "γ σ(Wᵤ(jω)⁻¹)",
        "γ σ(Wₜ(jω)⁻¹)",
    ],
    colors = [:red, :blue, :green],
)

    sensitivityfunctions, weightfunctions, γ = p.args[1:3]

    title --> "Specification sigma plot"
    xguide --> "Frequency ($(hz ? "Hz" : "rad/s"))", yguide --> "Singular Values $_PlotScaleStr"

    w = [10^i for i in range(wint[1], stop = wint[2], length = wnum)]

    # Plot the sensitivity functions
    for (index, G) in enumerate(sensitivityfunctions)
        if G isa Number || G isa LTISystem
            singval = sigma(ss(G), w)[1]
            if ControlSystems._PlotScale == "dB"
                singval = 20 * log10.(singval)
            end
            for i = 1:size(singval, 2)
                @series begin
                    xscale --> :log10
                    yscale --> ControlSystems._PlotScaleFunc
                    linestyle --> :solid
                    color --> colors[mod(index - 1, 3)+1]
                    label --> (i == 1 ? s_labels[mod(index - 1, 3)+1] : "")
                    w ./ (hz ? 2pi : 1), singval[:, i]
                end
            end
        end
    end

    ## Plot the weight functions
    for (index, W) in enumerate(weightfunctions)
        if W isa Number
            W = ss(float(W))
        elseif W == []
            continue
        end
        singval = sigma(γ / W, w)[1]
        if ControlSystems._PlotScale == "dB"
            singval = 20 * log10.(singval)
        end
        for i = 1:size(singval, 2)
            weightlabel = (i == 1 ? w_labels[mod(index - 1, 3)+1] : "")
            @series begin
                xscale --> :log10
                yscale --> ControlSystems._PlotScaleFunc
                linestyle --> :dash
                color --> colors[mod(index - 1, 3)+1]
                linewidth --> 2
                label --> (i == 1 ? w_labels[mod(index - 1, 3)+1] : "")
                w ./ (hz ? 2pi : 1), singval[:, i]
            end
        end
    end
end

# Case where a single sensitivity function (for instance the closed loop TF from
# disturbance to output) and the gain γ
specificationplot(sens::T, γ::Number; kwargs...) where {T<:LTISystem} =
    specificationplot([sens], [1], γ; kwargs...)

@userplot MvNyquistplot
"""
    fig = nyquistplot(sys, w;  unit_circle=false, hz = false, kwargs...)

Create a Nyquist plot of the `LTISystem`. A frequency vector `w` must be
provided.

- `unit_circle`: if the unit circle should be displayed
If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.

`kwargs` is sent as argument to plot.
"""
nyquistplot
@recipe function nyquistplot(p::MvNyquistplot;  unit_circle=false, hz=false)
    sys, w = p.args
    nw = length(w)
    framestyle --> :zerolines
    θ = range(0, stop=2π, length=100)
    S, C = sin.(θ), cos.(θ)
    L = freqresp(sys, w)
    dets = map(axes(L, 1)) do i
        det(I + L[i,:,:])
    end
    dets = vec(dets)
    redata = real.(dets)
    imdata = imag.(dets)
    if eltype(imdata) <: AbstractFloat
        ylims --> (clamp(minimum(imdata), -10, -1), clamp(maximum(imdata), 1, 10))
        xlims --> (clamp(minimum(redata), -10, -1), clamp(maximum(redata), 1, 10))
    end
    title --> "Mutivariable Nyquist plot"

    @series begin
        hover --> [hz ? Printf.@sprintf("f = %.3f", w/2π) : Printf.@sprintf("ω = %.3f", w) for w in w]
        (redata, imdata)
    end
    @series begin # Mark the critical point
        primary := false
        markershape := :xcross
        seriescolor := :red
        markersize := 5
        seriesstyle := :scatter
        [0], [0]
    end
    if unit_circle
        @series begin
            primary := false
            linestyle := :dash
            linecolor := :gray
            seriestype := :path
            markershape := :none
            (C, S)
        end
    end
end


text:
## ------------------------------------------------------------------
let
    @info("Testing FileContentEvent")
    e = FileContentEvent()
    file = joinpath(tempdir(), "FileContentEvent.test")
    write(file, "0")
    update!(e, file)

    events_count = 0
    trigger_at = [3, 4, 5]
    for it in 1:10
        touch(file)
        if has_event!(e, file)
            println(basename(file), " changed!!!")
            events_count += 1
        end

        if it in trigger_at
            write(file, it)
        end
        sleep(0.1)
    end
    @test events_count == length(trigger_at)
end


text:
abstract type PreAllocArrays end
abstract type AbstractDerivatives{T} <: PreAllocArrays end

struct Jacobian{T<:AbstractFloat} <: AbstractDerivatives{T}

    jac_phi::Matrix{T}
    jac_kick::Matrix{T}
    jac_copy::Matrix{T}
    jac_ij::Matrix{T}
    jac_tmp1::Matrix{T}
    jac_tmp2::Matrix{T}
    jac_err1::Matrix{T}
    dqdt_ij::Vector{T}
    dqdt_phi::Vector{T}
    dqdt_kick::Vector{T}

    function Jacobian(::Type{T},n::Integer) where T<:AbstractFloat
        sevn::Int64 = 7*n
        jac_phi = zeros(T,sevn,sevn)
        jac_kick = zeros(T,sevn,sevn)
        jac_copy = zeros(T,sevn,sevn)
        jac_ij = zeros(T,14,14)
        jac_tmp1 = zeros(T,14,sevn)
        jac_tmp2 = zeros(T,14,sevn)
        jac_err1 = zeros(T,14,sevn)
        dqdt_ij = zeros(T,14)
        dqdt_phi = zeros(T,sevn)
        dqdt_kick = zeros(T,sevn)

        return new{T}(jac_phi,jac_kick,jac_copy,jac_ij,jac_tmp1,jac_tmp2,jac_err1,
                    dqdt_ij,dqdt_phi,dqdt_kick)
    end
end

struct dTime{T<:AbstractFloat} <: AbstractDerivatives{T}

    jac_phi::Matrix{T}
    jac_kick::Matrix{T}
    jac_ij::Matrix{T}
    dqdt_phi::Vector{T}
    dqdt_kick::Vector{T}
    dqdt_ij::Vector{T}
    dqdt_tmp1::Vector{T}
    dqdt_tmp2::Vector{T}
    jac_kepler::Matrix{T}
    jac_mass::Vector{T}

    function dTime(::Type{T},n::Integer) where T<:AbstractFloat
        sevn::Int64 = 7*n
        jac_phi = zeros(T,sevn,sevn)
        jac_kick = zeros(T,sevn,sevn)
        jac_ij = zeros(T,14,14)
        dqdt_phi = zeros(T,sevn)
        dqdt_kick = zeros(T,sevn)
        dqdt_ij = zeros(T,14)
        dqdt_tmp1 = zeros(T,14)
        dqdt_tmp2 = zeros(T,14)
        jac_kepler = zeros(T,6,8)
        jac_mass = zeros(T,6)
        return new{T}(jac_phi,jac_kick,jac_ij,dqdt_phi,dqdt_kick,dqdt_ij,
            dqdt_tmp1,dqdt_tmp2,jac_kepler,jac_mass)
    end
end

struct Derivatives{T<:AbstractFloat} <: AbstractDerivatives{T}
    jac_phi::Matrix{T}
    jac_kick::Matrix{T}
    jac_copy::Matrix{T}
    jac_ij::Matrix{T}
    jac_tmp1::Matrix{T}
    jac_tmp2::Matrix{T}
    jac_err1::Matrix{T}
    dqdt_phi::Vector{T}
    dqdt_kick::Vector{T}
    dqdt_ij::Vector{T}
    dqdt_tmp1::Vector{T}
    dqdt_tmp2::Vector{T}
    jac_kepler::Matrix{T}
    jac_mass::Vector{T}
    dadq::Array{T,4}
    dotdadq::Matrix{T}
    tmp7n::Vector{T}
    tmp14::Vector{T}
    ctime::Vector{T}

    function Derivatives(::Type{T},n::Integer) where T<:AbstractFloat
        sevn::Int64 = 7*n
        jac_phi = zeros(T,sevn,sevn)
        jac_kick = zeros(T,sevn,sevn)
        jac_copy = zeros(T,sevn,sevn)
        jac_ij = zeros(T,14,14)
        jac_tmp1 = zeros(T,14,sevn)
        jac_tmp2 = zeros(T,14,sevn)
        jac_err1 = zeros(T,14,sevn)
        dqdt_phi = zeros(T,sevn)
        dqdt_kick = zeros(T,sevn)
        dqdt_ij = zeros(T,14)
        dqdt_tmp1 = zeros(T,14)
        dqdt_tmp2 = zeros(T,14)
        jac_kepler = zeros(T,6,8)
        jac_mass = zeros(T,6)
        dadq = zeros(T,3,n,4,n)
        dotdadq = zeros(T,4,n)

        tmp7n = zeros(T,sevn)
        tmp14 = zeros(T,14)
        ctime = zeros(T,1)

        return new{T}(jac_phi,jac_kick,jac_copy,jac_ij,jac_tmp1,jac_tmp2,jac_err1,
                      dqdt_phi,dqdt_kick,dqdt_ij,dqdt_tmp1,dqdt_tmp2,jac_kepler,jac_mass,
                      dadq,dotdadq,tmp7n,tmp14,ctime)
    end
end

"""Zero out each array in Derivatives."""
@generated function zero_out!(d::Derivatives{T}) where T<:AbstractFloat
    exprs = Vector{Expr}()
    for f in fieldnames(Derivatives)
        expr = :(d.$f .= 0.0)
        push!(exprs,expr)
    end
    return Expr(:block, exprs..., nothing)
end


text:
module TestReferences

using Test
using Mimi

import Mimi:
    reset_compdefs, getproperty, @defmodel

reset_compdefs()

@defcomp Foo begin
    input = Parameter()
    intermed = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.intermed[t] = p.input
    end
end

@defcomp Bar begin
    intermed = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.output[t] = p.intermed[t]
    end
end

@defmodel m begin
    index[time] = [1]
    component(Foo)
    component(Bar)

    Foo.input = 3.14
    Foo.intermed => Bar.intermed
end

run(m)

@test m[:Bar, :output][1] == 3.14

end


text:
export Odb, OdbObject, OdbWrite, OdbRead,
       exists, data, open_wstream, read_header

#TODO: rename to GitOdb
type Odb
    ptr::Ptr{Void}

    function Odb(ptr::Ptr{Void})
        @assert ptr != C_NULL
        o = new(ptr)
        finalizer(o, free!)
        return o
    end
end

free!(o::Odb) = begin
    if o.ptr != C_NULL
        ccall((:git_odb_free, libgit2), Void, (Ptr{Void},), o.ptr)
        o.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, o::Odb) = o.ptr

exists(o::Odb, id::Oid)  = bool(ccall((:git_odb_exists, libgit2), Cint, (Ptr{Void}, Ptr{Oid}), o, &id))
Base.in(id::Oid, o::Odb) = exists(o, id)

function read_header(odb::Odb, id::Oid)
    nbytes = Csize_t[0]
    otype  = Cint[0]
    @check ccall((:git_odb_read_header, libgit2), Cint,
                 (Ptr{Csize_t}, Ptr{Cint}, Ptr{Void}, Ptr{Oid}), nbytes, otype, odb, &id)
    return {:type => gitobj_const_type(otype[1]),
            :nbytes => nbytes[1]}
end

type OdbObject{T<:GitObject}
    ptr::Ptr{Void}
end

OdbObject(ptr::Ptr{Void}) = begin
    @assert ptr != C_NULL
    obj_type = ccall((:git_odb_object_type, libgit2), Cint, (Ptr{Void},), ptr)
    T = gitobj_const_type(obj_type)
    o = OdbObject{T}(ptr)
    finalizer(o, free!)
    return o
end

free!(o::OdbObject) = begin
    if o.ptr != C_NULL
        ccall((:git_odb_object_free, libgit2), Void, (Ptr{Void},), o.ptr)
        o.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, o::OdbObject) = o.ptr

Base.sizeof(o::OdbObject) = int(ccall((:git_odb_object_size, libgit2), Csize_t, (Ptr{Void},), o))
Base.length(o::OdbObject) = div(sizeof(o), sizeof(Uint8))

function data(o::OdbObject)
    data_ptr = ccall((:git_odb_object_data, libgit2), Ptr{Uint8}, (Ptr{Void},), o)
    @assert data_ptr != C_NULL
    return bytestring(data_ptr)
end

abstract OdbIO

free!(os::OdbIO) = begin
    if os.ptr != C_NULL
        #TODO: close before gc?
        ccall((:git_odb_stream_free, libgit2), Void, (Ptr{Void},), os.ptr)
        os.ptr = C_NULL
    end
end

type OdbWrite <: OdbIO
    ptr::Ptr{Void}
    id::Oid

    function OdbWrite(ptr::Ptr{Void})
        @assert ptr != C_NULL
        oid_ptr = ccall((:git_odb_object_id, libgit2), Ptr{Oid}, (Ptr{Void},), ptr)
        @assert oid_ptr != C_NULL
        s = new(ptr, Oid(oid_ptr))
        finalizer(s, free!)
        return s
    end
end

Base.convert(::Type{Ptr{Void}}, o::OdbWrite) = o.ptr

Oid(odbw::OdbWrite) = odbw.id

function open_wstream{T<:GitObject}(::Type{T}, odb::Odb, len::Int)
    @assert odb.ptr != C_NULL
    @assert len > 0
    gtype = git_otype(T)
    stream_ptr = Ptr{Void}[0]
    @check ccall((:git_odb_open_wstream, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Csize_t, Cint), stream_ptr, odb, len, gtype)
    return OdbWrite(stream_ptr[1])
end

#=
Base.isreadable(io::OdbWrite)  = false
Base.iswriteable(io::OdbWrite) = true
=#

#TODO: this needs to be reworked
Base.write(io::OdbWrite, buffer::ByteString) = begin
    len = length(buffer)
    @check ccall((:git_odb_stream_write, libgit2), Cint,
                 (Ptr{Void}, Ptr{Uint8}, Csize_t), io, buffer, len)
    return len
end

#TODO: this is broken...
Base.write{T}(io::OdbWrite, buffer::Array{T}) = begin
    @assert isbits(T)
    ptr = convert(Ptr{Uint8}, b)
    len = convert(Csize_t, div(length(b) * sizeof(T), sizeof(Uint8)))
    @check ccall((:git_odb_stream_write, libgit2), Cint,
                 (Ptr{Void}, Ptr{Uint8}, Csize_t), io, buffer, len)
    return io
end

Base.close(io::OdbWrite) = begin
    ioid = Oid(io)
    @check ccall((:git_odb_stream_finalize_write, libgit2), Cint, (Ptr{Oid}, Ptr{Void}), &ioid, io)
    return io
end

type OdbRead <: OdbIO
    ptr::Ptr{Void}
end

Base.convert(::Type{Ptr{Void}}, o::OdbRead) = o.ptr
#=
Base.isreadable(io::OdbRead)  = true
Base.iswriteable(io::OdbRead) = false
=#
Base.readbytes!(io::OdbRead, buffer::Vector{Uint8}, nb=length(b)) = begin
    @assert io.ptr != C_NULL
    len = convert(Csize_t, length(b))
    ret = @check ccall((:git_odb_stream_read, libgit2), Cint,
                       (Ptr{Void}, Ptr{Uint8}, Csize_t), io, buffer, len)
    @assert len > 0
    return len
end

Base.close(os::OdbRead) = return nothing # no op


text:
## --------------------------------------------------------------------
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Human1")

import Chemostat_Human1
const H1 = Chemostat_Human1

## --------------------------------------------------------------------
# Downloading
if isdir(H1.HUMAN1_PUBLICATION_RAW_DATA_DIR)
    println(relpath(H1.HUMAN1_PUBLICATION_RAW_DATA_DIR), " already exist, to force a re-unzip delete the folder")
elseif isfile(H1.HUMAN1_PUBLICATION_ZIP_FILE)
    println(relpath(H1.HUMAN1_PUBLICATION_ZIP_FILE), " already exist, to force a re-download delete it")
    run(`unzip $(H1.HUMAN1_PUBLICATION_ZIP_FILE) -d $(H1.RAW_DATA_DIR)`)
else
    println("downloading from: ", H1.HUMAN1_PUBLICATION_ZENODO_LINK)
    run(`curl $(H1.HUMAN1_PUBLICATION_ZENODO_LINK) --output $(H1.HUMAN1_PUBLICATION_ZIP_FILE)`)
    run(`unzip $(H1.HUMAN1_PUBLICATION_ZIP_FILE) -d $(H1.RAW_DATA_DIR)`)
end
!isdir(H1.HUMAN1_PUBLICATION_RAW_DATA_DIR) &&
    error("$(H1.HUMAN1_PUBLICATION_RAW_DATA_DIR) not found after download && unzip!!!")
println(relpath(H1.HUMAN1_PUBLICATION_RAW_DATA_DIR), " ready!!!")
flush(stdout)
flush(stderr)


text:
struct BackupConfig
    target::String
    dst::String
    id::Int
    backups::Vector{Int}
end

BackupConfig(target::AbstractString, dst::AbstractString, backups::AbstractVector{Int}) = BackupConfig(target, dst, isempty(backups) ? 1 : last(backups) + 1, backups)
BackupConfig(target::AbstractString, dst::AbstractString) = BackupConfig(target, dst, existing_backups(dst))

backup_path(dst, id) = joinpath(dst, string(id))

function flags(config::BackupConfig)
    _flags = ["-a", "--delete", "--one-file-system", "--info=progress2", "--no-inc-recursive"]
    !isempty(config.backups) && push!(_flags, string("--link-dest=", backup_path(config.dst, last(config.backups))))
    _flags
end

function make_backup(config::BackupConfig; confirm::Bool = true)
    dst = backup_path(config.dst, string(config.id))
    src = config.target

    includes = joinpath.(src, ["home", "etc"])
    excludes = joinpath.(src, ["home/Games", "home/.local/share/Steam"])

    cmd = `sudo rsync $(flags(config)) --exclude=$excludes $includes $dst/`
    println("The following command will be run: ", crayon"#aacc22", string(cmd), RESET)

    if confirm
        await_confirmation(true; msg=string("Ready to back up ", crayon"#aa77bb", src, RESET, " to ", crayon"#aabbff", dst, RESET, ". Continue?")) || return
    end

    @info "Creating backup of $src at $dst"
    mkpath(dst)
    origin_saves = joinpath(src, "home/belmant/Games/origin/drive_c/users/belmant/My Documents/Bioware")
    if isdir(origin_saves)
        # copy save files from the games
        cp(origin_saves, joinpath(src, "home/belmant/.origin_saves"), force=true)
    end
    try
        run(cmd)
    catch e
        # exit code 24 indicates that some files vanished before transfer; it is likely to happen if there is user activity during transfer and is considered fine
        (e isa ProcessFailedException && first(e.procs).exitcode == 24) || rethrow(e)
    end
end

function existing_backups(dst::AbstractString)
    !isdir(dst) && return Int[]
    parse.(Int, filter(x -> !isnothing(match(r"^\d+$", x)), readdir(dst)))
end