text:
################################################################################
# PHeap.jl
#
# The Persistent heap (priority queue) type.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020-2021 Noah C. Benson

import Random

# ==============================================================================
# PHeap definition.
"""
    PHeap{T,W,F}

A `PHeap` object is a priority queue to which objects can be `push`ed and given
a particular weight. `PHeap`s can then be `iterate`d in weighted order, `pop`ped
in weighted order, and the `first`  element may be obtained in weighted order.
Additionally, a random sample may be drawn from a heap, which uses the
distribution implied by the relative weights of the items in the heap.

`PHeap` objects may be initialized with a single argument of type
`F <: Function`, in which case this function is used to compare the weights for
ordering. By default this is `>` such that high weights are retrieved first by
the `first()` function and removed first by the `pop()` function. Note, however,
that the relative weights are still used by the `rand` function regardless of
the ordering function, so negative values cannot be used.

See also: [`PWSet`](@ref), [`PWDict`](@ref), [`pop`](@ref), `Base.first`.
"""
struct PHeap{T,W<:Number,F<:Function,D<:AbstractPDict{T,Int}}
    _heap::PVector{Tuple{T,W,W}} # value, weight, total subtree weight
    _index::D
    _compare::F
end
PHeap{T,W,F,D}(f::F) where {T, W<:Number, F<:Function, D<:AbstractPDict} = PHeap{T,W,F,D}(
    PVector{Tuple{T,W,W}}(),
    D{T,Int}(),
    f)
PHeap{T,W,F,D}(f::F) where {T, W<:Number, F<:Function, D<:AbstractPDict{T,Int}} = PHeap{T,W,F,D}(
    PVector{Tuple{T,W,W}}(),
    D(),
    f)
PHeap{T,W,D}(f::F) where {T, W<:Number, F<:Function, D<:AbstractPDict} = PHeap{T,W,F,D}(
    PVector{Tuple{T,W,W}}(),
    D{T,Int}(),
    f)
PHeap{T,W,D}(f::F) where {T, W<:Number, F<:Function, D<:AbstractPDict{T,Int}} = PHeap{T,W,F,D}(
    PVector{Tuple{T,W,W}}(),
    D(),
    f)
PHeap{T,W,D}()   where {T, W<:Number, D<:AbstractPDict} = PHeap{T,W,D}(>)
PHeap{T,W}(f::F) where {T, W<:Number, F<:Function} = PHeap{T,W,PDict}(f)
PHeap{T,W}()     where {T, W<:Number} = PHeap{T,W,PDict}(>)
PHeap{T,D}(f::F) where {T, F<:Function, D<:AbstractPDict} = PHeap{T,Float64,D}(f)
PHeap{T,D}()     where {T, D<:AbstractPDict} = PHeap{T,Float64,D}(>)
PHeap{T}(f::F)   where {T, F<:Function} = PHeap{T,Float64,F,PDict}(f)
PHeap{T}()       where {T} = PHeap{T,Float64,PDict}(>)
PHeap(f::F)      where {F <: Function} = PHeap{Any,Float64,F,PDict}(f)
PHeap() = PHeap{Any,Float64,PDict}(>)
_pheap_swap(heap::PVector{Tuple{T,W,W}}, index::D, ii::Int, jj::Int, subii::Bool, subjj::Bool) where {T,W,D<:AbstractPDict{T,Int}} = begin
    (tii,wii,totii) = heap[ii]
    (tjj,wjj,totjj) = heap[jj]
    newtotii = totjj + wii - (subjj ? wjj : 0)
    newtotjj = totii + wjj - (subii ? wii : 0)
    heap = setindex(heap, (tjj, wjj, newtotjj), ii)
    heap = setindex(heap, (tii, wii, newtotii), jj)
    index = setindex(index, ii, tjj)
    index = setindex(index, jj, tii)
    return (heap, index)
end
_pheap_fix_tot(heap::PVector{Tuple{T,W,W}}, ii::Int, dw::W) where {T,W} = begin
    dw == 0 && return heap
    while ii > 0
        node = heap[ii]
        heap = setindex(heap, (node[1], node[2], node[3] + dw), ii)
        ii = div(ii, 2)
    end
    return heap
end
_pheap_fix_up(heap::PVector{Tuple{T,W,W}}, index::D, cmp::F, ii::Int, oldw::W) where {T,W,F<:Function,D<:AbstractPDict{T,Int}} = begin
    # This change affects the weight totals of nodes above us by this abount:
    newnode = heap[ii]
    neww = newnode[2]
    dw = neww - oldw
    zW = zero(W)
    # if ii is 1, then we don't need to worry about the parent nodes and
    # swapping on the way up in order to fix things.
    while ii > 1
        pp = div(ii, 2)
        parent = heap[pp]
        # Break when we don't need to percolate up anymore.
        cmp(neww, parent[2]) || break
        # We need to swap with the parent
        (heap, index) = _pheap_swap(heap, index, ii, pp, true, false)
        ii = pp
    end
    # At this point we may have gone all the way up the heap, but we might have
    # only traversed up some of it. In the latter case, if the change in the
    # weight is not 0, then we need to fix weight totals up the heap.
    heap = _pheap_fix_tot(heap, div(ii, 2), dw)
    return (heap, index)
end
_pheap_fix_down(heap::PVector{Tuple{T,W,W}}, index::D, cmp::F, ii::Int) where {T,W,F<:Function,D<:AbstractPDict{T,Int}} = begin
    # We can assume when this function is called that the totals are correct for
    # the entire heap at this point.
    n = length(heap)
    zW = zero(W)
    while true
        lch = ii*2
        # If we've reached the end, we can return the (updated?) heap and index.
        (lch > n) && break
        # Otherwise, we need to check the children
        rch = lch + 1
        node = heap[ii]
        lnode = heap[lch]
        if rch > n
            # Only check the left...
            if cmp(lnode[2], node[2])
                # Swap these two
                (heap, index) = _pheap_swap(heap, index, ii, lch, true, true)
            end
            # There can't be children beyond this point.
            break
        else
            rnode = heap[rch]
            if cmp(lnode[2], rnode[2])
                cmp(node[2], lnode[2]) && break
                # lnode < node and lnode < rnode
                (heap, index) = _pheap_swap(heap, index, ii, lch, true, true)
                ii = lch
            else
                cmp(node[2], rnode[2]) && break
                # rnode < node and rnode < lnode
                (heap, index) = _pheap_swap(heap, index, ii, rch, true, true)
                ii = rch
            end
        end
    end
    return (heap, index)
end
_pheap_fixw(heap::PVector{Tuple{T,W,W}}, index::D, cmp::F, ii::Int, neww::W) where {T,W,F<:Function,D<:AbstractPDict{T,Int}} = begin
    node = heap[ii]
    oldw = node[2]
    dw = neww - oldw
    (dw == 0) && return (heap, index)
    # Fix this in the current node.
    heap = setindex(heap, (node[1], neww, node[3] + dw), ii)
    # If dw is positive then there is no chance that we will need to
    # swap with the children; on the flip side, if it is negative,
    # we still need to fix totals up the heap. So we can do that
    # now.
    (heap, index) = _pheap_fix_up(heap, index, cmp, ii, oldw)
    # Now, we only need to check the lower heap nodes if the weight
    # is now lower. If it is now lower, then ii will still point
    # to the same node after the _pheap_fix_up() call.
    if cmp(oldw, neww)
        (heap, index) = _pheap_fix_down(heap, index, cmp, ii)
    end
    # That is all.
    return (heap, index)
end
push(p::PHeap{T,W,F,D}, tw::Tuple{S,X}) where {T,W,F,D,S,X<:Number} = begin
    (t,w) = tw
    (w <= 0) && throw(ArgumentError("Weights must be positive"))
    cmp = p._compare
    heap = p._heap
    index = p._index
    # If t is already in the heap, we are just editing the weight.
    ii = get(index, t, 0)
    if ii == 0
        # Since t isn't already in the heap, we need to add it to the end
        # with a weight of 0 (which keeps the tree valid and all of the
        # weight totals in tact) then fix its weight.
        zW = zero(W)
        heap = push(heap, (t,zW,zW))
        ii = length(heap)
        index = setindex(index, ii, t)
    end
    (heap, index) = _pheap_fixw(heap, index, cmp, ii, W(w))
    (p._heap === heap && p._index === index) && return p
    return PHeap{T,W,F,D}(heap, index, cmp)
end
push(p::PHeap{T,W,F,D}, tw::Pair{S,X}) where {T,W,F,D,S,X<:Number} = begin
    return push(p, (tw[1],tw[2]))
end
"""
    setweight(h, x, w)

Sets the weight of the value x in the given persistent heap or persistent
weighted collection object h to be w and yields the new updated version of h.
If x is not already in the object h, then an error is thrown.

See also: [`getweight`](@ref), [`PWDict`](@ref), [`PWSet`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWSet{Symbol, ?Float64} with 3 elements:"
julia> d = PWSet{Symbol}(:a => 1.0, :b => 2.0, :c => 3.0)
PWSet{Symbol,Float64} with 3 elements:
  :c
  :b
  :a

julia> setweight(d, :a, 4.0)
PWSet{Symbol,Float64} with 3 elements:
  :a
  :c
  :b
```
"""
setweight(p::PHeap{T,W,F,D}, t::S, w::X) where {T,W,F,D,S,X<:Number} = begin
    (w <= 0) && throw(ArgumentError("Weights must be positive"))
    ii = get(p._index, t, 0)
    (ii == 0) && error("cannot setweight on item that is not in the PHeap")
    (heap, index) = _pheap_fixw(p._heap, p._index, p._compare, ii, w)
    (heap === p._heap && index === p._index) && return p
    return PHeap{T,W,F,D}(heap, index, p._compare)
end
"""
    getweight(h, x)

Yields the weight for the given value x in the given persistent heap or
persistent weighted collection h. If h is a weighted dictionary, then x
should be the key.

If the key or object x is not found in the collection h, then 0 is returned.

See also: [`setweight`](@ref), [`PWDict`](@ref), [`PWSet`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWSet{Symbol, ?Float64} with 3 elements:"
julia> d = PWSet{Symbol}(:a => 1.0, :b => 2.0, :c => 3.0)
PWSet{Symbol,Float64} with 3 elements:
  :c
  :b
  :a

julia> getweight(d, :a)
1.0
```
"""
getweight(h::PHeap{T,W,F,D}, t::S) where {T,W,F,D,S} = begin
    ii = get(h._index, t, 0)
    (ii == 0) && return 0
    return h._heap[ii][2]
end
_pheap_delete(p::PHeap{T,W,F,D}, ii::Int) where {T,W,F,D} = begin
    n = length(p)
    cmp = p._compare
    (n == 1) && return PHeap{T,W,F,D}(pop(p._heap), empty(p._index), cmp)
    # Grab the very last node, then we can prep the removal by updating its
    # weight to be zero (this subtracts its weight from the tree's totals).
    (tn,wn,totn) = p._heap[n]
    (heap, index) = _pheap_fixw(p._heap, p._index, cmp, n, zero(W))
    # Now we replace the node in slot ii with the last node, but keep the
    # old weight so that the tree's totals are still valid.
    (tii, wii, totii) = heap[ii]
    heap = setindex(pop(heap), (tn, wii, totii), ii)
    index = setindex(delete(index, tii), ii, tn)
    # Then set the weight at that node to be that of its old value!
    (heap, index) = _pheap_fixw(heap, index, cmp, ii, wn)
    # Thwt will leave us in a valid state; return the new heap.
    (heap === p._heap && index === p._index) && return p
    return PHeap{T,W,F,D}(heap, index, cmp)
end
pop(p::PHeap{T,W,F,D}) where {T,W,F,D} = begin
    n = length(p)
    (n == 0) && throw(ArgumentError("PHeap must be non-empty"))
    return _pheap_delete(p, 1)
end
delete(p::PHeap{T,W,F,D}, t::S) where {T,W,F,D,S} = begin
    # Find this particular node we need to remove
    ii = get(p._index, t, 0)
    return ii == 0 ? p : _pheap_delete(p, ii)
end
Base.length(p::PHeap{T,W,F,D}) where {T,W,F,D} = length(p._index)
Base.iterate(p::PHeap{T,W,F,D}) where {T,W,F,D} = iterate(p, p)
Base.iterate(::PHeap{T,W,F,D}, p::PHeap{T,W,F,D}) where {T,W,F,D} = begin
    return (length(p) == 0 ? nothing : (first(p), pop(p)))
end
Base.first(p::PHeap{T,W,F,D}) where {T,W,F,D} = begin
    (length(p) == 0) && throw(ArgumentError("PHeap must be non-empty"))
    return p._heap[1][1]
end
Base.in(x::S, p::PHeap{T,W,F,D}) where {S,T,W,F,D} = begin
    return (get(p._index, x, 0) != 0)
end
Random.rand(p::PHeap{T,W,F,D}) where {T,W,F,D} = begin
    n = length(p)
    (n == 0) && throw(ArgumentError("PHeap must be non-empty"))
    x = Random.rand(Float64) * p._heap[1][3]
    ii = 1
    (t,w,tot) = p._heap[ii]
    while true
        (x <= w) && return t
        x -= w
        ii *= 2
        (t,w,tot) = p._heap[ii]
        (x <= tot) && continue
        x -= tot
        ii += 1
        (t,w,tot) = p._heap[ii]
    end
    error("invalid state reached")
end

# ==============================================================================
# Core Air API methods and related Base methods
Base.isequal(a::PHeap, b::PHeap) = begin
    (length(a) == length(b)) || return false
    while length(a) > 0
        f = first(a)
        isequal(f, first(b)) || return false
        (getweight(a, f) == getweight(b, f)) || return false
        a = pop(a)
        b = pop(b)
    end
    return true
end
Base.hash(a::PHeap{T,W,F,D}) where {T,W,F,D}  = begin
    h = 0x11 * UInt(length(a))
    for (t,w,tot) in a._heap
        h += hash(t) ⊻ hash(w)
    end
    return h
end

# Export the relevant symbols (PHeap is considered internal).
export getweight, setweight


text:
function startofbags(k::Vector{Int})
  s = zeros(Int,length(k))
  for i in 2:length(k)
    s[i] = s[i-1] + k[i - 1]
  end
  s
end

"""
    scatter(x::Matrix,k::Int)
    scatter(x::Matrix,k::Vector{Int})

    repeat each column of `x` k-times
    repeat j-th column of `x` k[j]-times

    ```juliadoc
    x = [1 3; 2 4]
    julia> FluxExtensions.scatter(x, 3)

    2×6 Array{Int64,2}:
     1  1  1  3  3  3
     2  2  2  4  4  4
     ```
"""
function scatter(x::Matrix,k::Int)
  xx = similar(x,size(x,1),k*size(x,2))
  @inbounds for j in 1:size(x,2)
    for b in 1:k
      for i in 1:size(x,1)
        xx[i,(j-1)*k+b] = x[i,j]
      end
    end
  end
  xx
end

function scatter(x::Matrix,k::Vector{Int})
  @assert size(x,2) >= length(k)
  xx = similar(x,size(x,1),sum(k))
  starts = startofbags(k)
  for j in 1:length(k)
    for b in 1:k[j]
      for i in 1:size(x,1)
        xx[i,starts[j] + b] = x[i,j]
      end
    end
  end
  xx
end

"""
    gather(x,k::Int)
    gather(f,x,k::Vector{Int})
    gather(f,x,k::Int)

    sum blocks of size k-columns
    reduce with function `f`

    ```juliadoc
    julia>  x = [ 1 1 1 3 3 3; 2 2 2 4 4 4]
    julia> FluxExtensions.gather(x,3)
    2×2 Array{Int64,2}:
     3   9
     6  12
    ```
"""
function gather(x,k::Int)
  d, l = size(x, 1), div(size(x,2), k)
  reshape(sum(reshape(x, d, k, l), dims=2), d, l)
end

function gather(x,k::Vector{Int})
  xx = similar(x,size(x,1),length(k)) .= 0
  starts = startofbags(k)
  for j in 1:length(k)
    for b in 1:k[j]
      for i in 1:size(x,1)
        xx[i,j] += x[i,starts[j] + b]
      end
    end
  end
  xx
end

function gather(f::Function, x, k::Int)
  d, l = size(x,1), div(size(x,2),k)
  reshape(sum(reshape(x,d,k,l),2),d,l)
end

function logsumexp(x,k::Int)
  xm = maximum(x,dims = k)
  xm .+ log.( sum( exp.(x .- xm),dims = k))
end

function logsumexp(x)
  xm = maximum(x)
  xm .+ log.( sum( exp.(x .- xm)))
end

scatter(x::Flux.Tracker.TrackedMatrix, k) = Flux.Tracker.track(scatter,x, k )
Flux.Tracker.@grad function scatter(x,k)
  return(scatter(Flux.data(x), k), Δ -> (gather(Δ, k), nothing))
end

ext: jl

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
using ProtoStructs, Test

@proto struct TestMe{T, V <: Real}
    A::Int
    B
    C::T
    D::V
end
test_me = @test_nowarn TestMe(1, "2", complex(1), 5)
test_me_kw = @test_nowarn TestMe(A=1, B="2", C=complex(1), D=5)

@testset "Construction" begin
    @test TestMe((A=1,)) isa TestMe
    @test_throws UndefKeywordError TestMe(A=1)
end # testset

@testset "Access" begin
    @test test_me.A == 1
    @test test_me.B == "2"
    @test test_me.C == complex(1)
end # testset

@testset "Properties" begin
    @test propertynames( test_me ) == (:A, :B, :C, :D)
end # testset

@proto struct TestMe{T, V <: Real}
    A::Int
    B
    C::T
    D::V
    E::String
end
test_me2 = @test_nowarn TestMe(1, "2", complex(1), 5, "tadaa")
test_me_kw2 = @test_nowarn TestMe(A=1, B="2", C=complex(1), D=5, E="tadaa")

@testset "Redefinition" begin
    @test length(methods(TestMe)) == 3
end # testset

ProtoStructs.@proto struct TestKw{T, V <: Real}
    A::Int = 1
    B = :no
    C::T = nothing
    D::V
    E::String
end

@testset "kwdef" begin
    tw = TestKw(D = 1.2, E = "yepp")
    @test tw isa TestKw
    @test tw.A == 1
    @test tw.B == :no
    @test tw.C === nothing
    @test tw.D == 1.2
    @test tw.E == "yepp"
end


text:
using DataFrames, CSV, DataFramesMeta, Hose, Statistics, StatsBase
using Plots, PlotThemes, StatPlots

theme(:dark)

optionally_convert_to_df(df::DataFrame) = df
optionally_convert_to_df(df) = df |> DataFrame

log_df = CSV.File("logs/log.csv") |> optionally_convert_to_df

bs = @hose log_df |>
@where(:Dataset .== "Bitstamp", :TaskTypes .== "Execution") |>
@select(:Language, :TaskNames, :TaskTypes, :TaskExecTimes) |>
@by([:Language, :TaskNames], TotalTime = sum(:TaskExecTimes))

@df bs groupedbar(:Language, :TotalTime, group = :TaskNames,
                  ylabel = "Time (s)", xticks = :all,
                  bar_position = :stack, legend=:topright, size=(400,600))
savefig("img/bitstamp_exec.pdf")
run(`convert -density 150 img/bitstamp_exec.pdf -quality 90 img/bitstamp_exec.png`)

hfl = @hose log_df |>
@where(:Dataset .== "HFlights", :TaskTypes .== "Execution") |>
@select(:Language, :TaskNames, :TaskTypes, :TaskExecTimes) |>
@by([:Language, :TaskNames], TotalTime = sum(:TaskExecTimes))

@df hfl groupedbar(:Language, :TotalTime, group = :TaskNames,
                  ylabel = "Time (s)", xticks = :all,
                  bar_position = :stack, legend=:topleft, size=(400,600))
savefig("img/hflights_exec.pdf")
run(`convert -density 150 img/hflights_exec.pdf -quality 90 img/hflights_exec.png`)

overhead = @hose log_df |>
@where(:TaskTypes .!= "Execution") |>
@select(:Dataset, :Language, :TaskNames, :TaskTypes, :TaskExecTimes) |>
@by([:Language, :TaskTypes], TotalTime = sum(:TaskExecTimes)) |>
@transform(TaskTypeLanguage = :TaskTypes .* " in " .* :Language)

@df overhead groupedbar(:Language, :TotalTime, group = :TaskTypes,
                  ylabel = "Time (s)", xticks = :all,
                  bar_position = :stack, legend=:topright, size=(400,600))
savefig("img/overhead.pdf")
run(`convert -density 150 img/overhead.pdf -quality 90 img/overhead.png`)
# execs = @hose log_df |>
# @where(:TaskTypes .== "Execution") |>
# @select(:Dataset, :Language, :TaskNames, :TaskTypes, :TaskExecTimes) |>
# @by([:Dataset, :Language, :TaskNames], TotalTime = sum(:TaskExecTimes)) |>
# @transform(DatasetLanguage = :Dataset .* " by " .* :Language)
#
# Plots.reset_defaults()
# @df execs groupedbar(:DatasetLanguage, :TotalTime, group = :TaskNames,
#                   ylabel = "Time (s)", xticks = :all,
#                   bar_position = :stack, legend=:topright, size=(2000,1200))
# png("img/exec.png")

max_stars_repo_name: RevoltBI/Julia-Meetup

max_issues_repo_name: RevoltBI/Julia-Meetup

max_forks_repo_name: RevoltBI/Julia-Meetup


text:
using PyPlot,
using Random, ProgressMeter
#using Plots

Base.@kwdef mutable struct TwoLayer
    d::Int  # input dim
    m::Int  # output dim
    W::Array{Float64, 2}  # the input and output weighst
    TwoLayer(d::Int, m::Int) = new(d, m, randn(m, d+1))
end

"""
Initialise a model
"""
function init!(l::TwoLayer, alpha1::Real=1.0, alpha2::Real=1.0, sym::Bool=false)
    d, m = l.d, l.m
    # input weights initialized uniformly on the unit sphere
    l.W[:,1:d] ./= sqrt.(sum(l.W[:,1:d].^2, dims=2))
    l.W[:, 1:d] .*= alpha1
    # output weights initialized by -1 or 1 with equal proba
    l.W[1:div(m,2),end] .= alpha2
    l.W[div(m,2)+1:end,end] .= -alpha2
    if sym
        @assert iseven(m)
        W[div(m,2)+1:end,1:d]  .=  W[1:div(m,2),1:d]
    end
end

##
#
"""
Train 2-layer relu NN with GD on square loss
INPUT
 m is the number of neurons,
 alpha1 (resp. alpha2): scale of the initialization of the input (resp. output) weights
 ouly_output: train just the output layer (instead of both)
 sym_init: trick to make the output = 0 exactly at initialization
OUTPUT
 Ws: the whole training path
 loss: evolution of loss during training
"""
function twonet_square(model::Layer, X, Y, m, stepsize, niter;alpha1=1.0, alpha2=1.0, weight_decay=0.0, sym_init=false, only_output = false)
    (n,d) = size(X) # n samples in R^d
    # initialize
    #
    # input weights initialized uniformly on the unit sphere
    init!(model)

    Ws    = zeros(m, d+1, niter) # store optimization path
    loss  = zeros(niter) # loss is -log of the empirical risk

    @showprogress 1 "Training neural network..." for iter = 1:niter
        Ws[:,:,iter] = W
        act  =  max.( W[:,1:end-1] * X', 0.0) # activations
        out  =  (1/m) * sum( W[:,end] .* act , dims=1)[:] # predictions of the network
        gradR   = (out .- Y)/n  # size n
        grad_w1 = (W[:,end] .* float.(act .> 0) * ( X .* gradR  )) + weight_decay * W[:,1:end-1] # gradient for input weights
        grad_w2 = act * gradR + weight_decay * W[:,end]  # gradient for output weights
        if only_output
            grad = cat(zeros(m,d), grad_w2, dims=2) # size (m × d+1)
        else
            grad = cat(grad_w1, grad_w2, dims=2) # size (m × d+1)
        end
        loss[iter] = (1/2)*sum( ( out - Y).^2 )/n + (weight_decay/2)*sum(W.^2)/m
        W = W - stepsize * grad
    end

    return Ws, loss
end

##

Random.seed!(5)
d=3  # dim(X) + bias
n = 10
m = 1000
X = cat(2*(rand(n,2).-0.5), ones(n), dims=2)
Y = rand([-1,1],n)
stepsize = 0.1
niter = 20000
Ws, loss = twonet_square(X, Y, m, stepsize, niter, weight_decay=5e-5, only_output =false)

loglog(loss,label="both")

stepsize = 8
niter = 40000
Ws2, loss2 = twonet_square(X, Y, m, stepsize, niter, weight_decay=5e-6, only_output =true)

loglog!(loss2)

#legend()

colm = :magma
# cm = colormap(colm)

#figure(figsize=[8,4])
f(x1,x2) = (1/m) * sum( Ws[:,end,end] .* max.( Ws[:,1:3,end] * [x1;x2;1], 0.0)) # prediction function
f2(x1,x2) = (1/m) * sum( Ws2[:,end,end] .* max.( Ws2[:,1:3,end] * [x1;x2;1], 0.0)) # prediction function
reso = 100
xs = range(-1.0,1.0,length=reso)
tab = [f(xs[i],xs[j]) for i=eachindex(xs), j=eachindex(xs)]
tab2 = [f2(xs[i],xs[j]) for i=eachindex(xs), j=eachindex(xs)]

# contour(xs, xs, f, fill=true, c=:magma)

#figure(figsize=[6.75,3])
#subplot(121)
#title("Training both layers")

#pcolormesh(xs', xs, tab', cmap=colm, shading="gouraud",vmin=-4.5,vmax = 1.5,edgecolor="face")

# cm = plt.get_cmap(colm)

scatter!(X[:,1],X[:,2],50,c=cm.(Y/6 .+0.7),edgecolors="k")

#axis("equal")
xticks([-1,0,1]);xlabel(L"(a)")
yticks([-1,0,1])#;ylabel(L"y")
subplot(122)
pcolormesh(xs', xs, tab2', cmap=colm, shading="gouraud",vmin=-4.5,vmax = 1.5,edgecolor="face")
cm = plt.get_cmap(colm)
scatter(X[:,1],X[:,2],50,c=cm.(Y/6 .+0.7),edgecolors="k")
#axis("equal")
#title("Training output layer")
xticks([-1,0,1]);xlabel(L"(b)")
yticks([-1,0,1])#;ylabel(L"y")
#savefig("regularized.png",dpi=150,bbox_inches="tight")
###

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


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

"""
    Geometry{Dim,T}

A geometry embedded in a `Dim`-dimensional space with coordinates of type `T`.
"""
abstract type Geometry{Dim,T} end

"""
    embeddim(geometry)

Return the number of dimensions of the space where the `geometry` is embedded.
"""
embeddim(::Type{<:Geometry{Dim,T}}) where {Dim,T} = Dim
embeddim(g::Geometry) = embeddim(typeof(g))

"""
    paramdim(geometry)

Return the number of parametric dimensions of the `geometry`. For example, a
sphere embedded in 3D has 2 parametric dimension (polar and azimuthal angles).
"""
paramdim(g::Geometry) = paramdim(typeof(g))

"""
    coordtype(geometry)

Return the machine type of each coordinate used to describe the `geometry`.
"""
coordtype(::Type{<:Geometry{Dim,T}}) where {Dim,T} = T
coordtype(g::Geometry) = coordtype(typeof(g))

"""
    centroid(geometry)

Return the centroid of the `geometry`.
"""
centroid(g::Geometry) = center(g)

"""
    measure(geometry)

Return the measure of the `geometry`, i.e. the length, area, or volume.
"""
function measure end

"""
    extrema(geometry)

Return the top left and bottom right corners of the
bounding box of the `geometry`.
"""
Base.extrema(g::Geometry) = extrema(boundingbox(g))

"""
    boundary(geometry)

Return the boundary of the `geometry`.
"""
function boundary end

"""
    g₁ ⊆ g₂

Tells whether or not the geometry `g₁` is a subset of geometry `g₂`.
"""
Base.issubset(g₁::Geometry, g₂::Geometry) = all(p ∈ g₂ for p in vertices(g₁))

"""
    p ∈ g

Tells whether or not the point `p` is in the geometry `g`.
"""
Base.in(p::Point, g::Geometry)

"""
    isconvex(geometry)

Tells whether or not the `geometry` is convex.
"""
isconvex(g::Geometry) = isconvex(typeof(g))

"""
    issimplex(geometry)

Tells whether or not the `geometry` is simplex.
"""
issimplex(::Type{<:Geometry}) = false
issimplex(g::Geometry) = issimplex(typeof(g))

# ----------------
# IMPLEMENTATIONS
# ----------------

include("polytopes.jl")
include("primitives.jl")

"""
    PointOrGeometry{Dim,T}

A union type that can represent either a point or a geometry.
"""
const PointOrGeometry{Dim,T} = Union{Point{Dim,T},Geometry{Dim,T}}

"""
    Multi(entities)

A collection of points or geometries seen as a single entity.

In geographic information systems (GIS) it is common to represent
multiple polygons as a single entity. In this case the polygons
refer to the same object in the real world (e.g. country with islands).
"""
struct Multi{Dim,T,I<:PointOrGeometry{Dim,T}} <: Geometry{Dim,T}
  items::Vector{I}
end

Base.getindex(multi::Multi, ind) = getindex(multi.items, ind)
Base.length(multi::Multi) = length(multi.items)
Base.eltype(multi::Multi) = eltype(multi.items)
Base.firstindex(multi::Multi) = firstindex(multi.items)
Base.lastindex(multi::Multi) = lastindex(multi.items)
Base.iterate(multi::Multi, state=1) =
  state > length(multi) ? nothing : (multi[state], state+1)

paramdim(multi::Multi) = maximum(paramdim, multi.items)

==(multi₁::Multi, multi₂::Multi) =
  length(multi₁) == length(multi₂) &&
  all(g -> g[1] == g[2], zip(multi₁, multi₂))

vertices(multi::Multi) =
  [vertex for geom in multi for vertex in vertices(geom)]

function centroid(multi::Multi)
  cs = coordinates.(centroid.(multi.items))
  Point(sum(cs) / length(cs))
end

measure(multi::Multi) = sum(measure, multi.items)

area(multi::Multi{Dim,T,<:Polygon}) where{Dim,T} = measure(multi)

boundary(multi::Multi) = Multi([boundary(geom) for geom in multi])

chains(multi::Multi{Dim,T,<:Polygon}) where {Dim,T} =
  [chain for geom in multi for chain in chains(geom)]

Base.in(point::Point, multi::Multi) = any(geom -> point ∈ geom, multi)

function Base.show(io::IO, multi::Multi{Dim,T}) where {Dim,T}
  n = length(multi.items)
  G = eltype(multi.items)
  print(io, "$n Multi-$(nameof(G)){$Dim,$T}")
end

function Base.show(io::IO, ::MIME"text/plain", multi::Multi)
  println(io, multi)
  print(io, io_lines(multi, "  "))
end


text:
using HDF5
using Test
using DSGE

m = PoolModel()
Φ1, F_ϵ1, F_λ1 = transition(m)
Φ2, F_ϵ2, F_λ2 = solve(m)

@testset "Check state-space system matches reference" begin
    @test Φ1 == Φ2
    @test F_ϵ1 == F_ϵ2
    @test F_λ1 == F_λ2
end


text:
using Pkg
Pkg.activate(".")

using BenchmarkTools
using Flux.Optimise:update!
using Revise
using Flux.Optimise
using Random
using Plots

includet("../models/REN.jl")
includet("../models/ffREN.jl")
includet("../models/solvers.jl")
includet("../models/models.jl")
includet("../models/robust_rnn.jl")

includet("../utils.jl")
includet("../data_handler.jl")
includet("../train.jl")

Random.seed!(123);

# Make a model
nu = 1
ny = 1
nx = 75
nv = 150

# nl(x) = (x + cos(x)*sin(x)) / 2
nl(x) = Flux.relu(x)

# model = stable_ff_rnn(Float64, nu, nx, nv, ny; ϵ=0.01)
# model = robust_rnn_stable(nu, nx, nv, ny; nl=nl)

model = bounded_ren(Float64, nu, nx, nv, ny, 0.25f0; ϵ=0.01)
solver = PeacemanRachford(tol=1E-4, α=1.0, maxIter=2000, verbose=false, cg=false)
(model::dissipative_ren)(x0, u) = model(x0, u, solver)  # capture solver

# download and wing flutter example data
seq_len = 1024
download_and_extract("wing_flutter")
train, val, stats = load_wing_flutter()

# Make training and validation datasets
ut = Iterators.partition(train[1], seq_len)
yt = Iterators.partition(train[2], seq_len)
train_data = zip(Iterators.partition(train[2], seq_len))

train_data = zip(ut, yt)
val_data = [val]

opt = Flux.Optimise.ADAM(1E-3)

# Test model on sample input
tloss = my_train!(train_data, val_data, model, opt; Epochs=200)
val_res = validate(model, val_data, stats)

# Try model on example input
x0 = init_state(model, 1)
yest, xend = model(x0, train[1])
plot(extract(yest, 1, 1))
plot!(extract(train[2], 1, 1))

norm(extract(yest, 1, 1)) / norm(extract(train[1], 1, 1))
gamma = estimate_lipschitz_lower(model; seq_len=200, maxIter=150,step_size=1E-1)

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
function ℱ(dmodel, μ0, x0, y0, z0)

end
using StaticArrays

"""
Equalizes offer and demand
# Arguments
* `dmodel::`:
* `μ0`:
* `x0`:
* `y0`:
* `z0`:
# Optional Argument
* `diff::boolean`: Indicates whether we want to compute differentials
# Returns
* `A`:
"""
function 𝒜(dmodel, μ0, x0, y0::SVector{d, Float64}, z0; diff=false) where d

    N = length(μ0)

    μ = μ0[:]

    s0 = repeat(dmodel.F.s0, length(x0.views))

    zvec = [z0 for n=1:N]
    yvec = [y0 for n=1:N]

    parms = SVector(dmodel.hmodel.calibration[:parameters]...)

    res = Dolark.equilibrium(dmodel.hmodel, s0, x0.data, yvec, zvec, yvec, zvec, parms)

    @assert length(res[1])==1 # only one aggregate condition for now.

    res = reinterpret(Float64, res)

    if diff==false

        res = sum(μ .* res) # this is 1d only
        return [res]

    else

        A = sum(μ .* res)

        e_x0, e_y0, e_z0, e_y1, e_z1 = Dolark.equilibrium(dmodel.hmodel, Val{(2,3,4,5,6)}, s0, x0.data, yvec, zvec, yvec, zvec, parms)

        n_x = size(e_x0[1],2)
        n_y = size(e_y0[1],2)
        n_z = size(e_z0[1],2)

        A_mu = d_mu -> sum(d_mu .* res)
        A_x = d_x -> sum(μ .* (e_x0 .* reinterpret(SVector{n_x, Float64}, d_x)))
        A_y = d_y -> sum(μ .* [(e_y0 .+ e_y1)[n]*(SVector{n_y, Float64}(d_y...)) for n=1:N] )
        A_z = d_z -> sum(μ .* [(e_z0 .+ e_z1)[n]*(SVector{n_z, Float64}(d_z...)) for n=1:N] )

        # A_mu = d_mu -> sum(d_mu .* res)
        # A_x = d_x -> sum(μ .* e_x0 .* d_x.data)
        # A_y = d_y -> sum(μ .* [(e_y0 .+ e_y1)[n]*d_y for n=1:N] )
        # A_z = d_z -> sum(μ .* [(e_z0 .+ e_z1)[n]*d_z for n=1:N] )

        n_y = length(e_x0[1])
        N_μ = length(μ)
        N_x = length(x0.data)
        n_z = size(e_z1[1],2)

        R_A_mu = LinearMaps.LinearMap(A_mu, n_y, N_μ)
        R_A_x = LinearMaps.LinearMap(A_x, n_y, N_x)
        R_A_y = LinearMaps.LinearMap(A_y, n_y, n_y)
        R_A_z = LinearMaps.LinearMap(A_z, n_y, n_z)

        A = [A] # this is 1d only
        return A, R_A_mu, R_A_x, R_A_y, R_A_z

    end

end

using LinearAlgebra: I

using LinearMaps

using LinearAlgebra: I

function Residual(dmodel, u::Unknown)

    z = SVector(dmodel.hmodel.calibration[:exogenous]...)
    p_ = SVector(dmodel.hmodel.calibration[:parameters]...)

    μ = u.μ[:]
    x = u.x
    y = u.y
    p = u.p

    N_μ = length(μ)
    N_x = length(x.data)*length(x.data[1])
    n_p = length(p)
    n_y = length(y)
    n_z = length(z)

    F = dmodel.F
    G = dmodel.G

    # quick computations if we don't compute derivatives
    # r_μ = μ - G(u.μ, u.x)
    r_F = F(x, x; exo=(p,p), set_future=true) # this is not equal to the normalized version
    # r_A = 𝒜(dmodel, μ, x, y, z)
    # r_p = Dolark.projection(dmodel.hmodel, y, z, p_)

    μ1, G_mu, G_x = G(u.μ, u.x, diff=true)

    r_μ = μ - μ1

    r_F = F(x, x; exo=(p,p), set_future=true) # this is not equal to the normalized version
    JJ = Dolo.df_A(F, x,x; exo=(p,p))
    L = Dolo.df_B(F, x,x; exo=(p,p))
    F_p1, F_p2 = Dolo.df_e(F, x, x, p,p)
    F_p = F_p1 + F_p2 # that one looks wrong

    # no preconditioning
    # Dolo.prediv!(L, JJ)
    # r_F = JJ\r_F
    # F_p = JJ\F_p

    Ft_x = LinearMaps.LinearMap(z->JJ*z, N_x, N_x)
    # Ft_x = LinearMaps.LinearMap(z->z, N_x, N_x)
    Ft_xx = LinearMaps.LinearMap(z->L*z, N_x, N_x)
    Ft_p = LinearMaps.LinearMap(z->F_p*z, N_x, n_p)

    r_p, r_p_y, r_p_z = Dolark.projection(dmodel.hmodel, Val{(0,1,2)}, y, z, p_)

    R_p_y = LinearMaps.LinearMap(z->r_p_y*z, n_p, n_y)
    R_p_z = LinearMaps.LinearMap(z->r_p_z*z, n_p, n_z)

    r_A, r_A_mu, r_A_x, r_A_y, r_A_z = 𝒜(dmodel, μ, x, y, z; diff=true)

    u = Dolark.Unknown(μ - μ1, p - r_p, r_F, r_A)

    J = [ I-G_mu                 zeros(N_μ, n_p)        -G_x               zeros(N_μ, n_y)         ; # μ
          zeros(n_p, N_μ)         I                     zeros(n_p, N_x)     -R_p_y        # p;
          zeros(N_x, N_μ)        Ft_p                  Ft_x+Ft_xx          zeros(N_x, n_y)         ; # x
          r_A_mu                 zeros(n_y, n_p)       r_A_x               r_A_y          # y
        ]

    # Y = [ -G_mu                 zeros(N_μ, n_p)        -G_x               zeros(N_μ, n_y)         ; # μ
    #     zeros(n_p, N_μ)           zeros(n_p, n_p)                   zeros(n_p, N_x)     -R_p_y        # p;
    #     zeros(N_x, N_μ)        Ft_p                  Ft_xx          zeros(N_x, n_y)         ; # x
    #     r_A_mu                 zeros(n_y, n_p)       r_A_x               r_A_y          # y
    #   ]

    return u, J

end

function Residual(dmodel, v::AbstractVector; diff=false)

    n1 = N_μ = length(dmodel.G.grid)
    n2 = n_p = length(dmodel.hmodel.factories[:projection].equations)
    n3 = N_x = length(dmodel.F.grid)
    n4 = n_y = length(dmodel.hmodel.symbols[:aggregate])

    n_z = length(dmodel.hmodel.symbols[:exogenous])
    n_x = length(dmodel.F.x0.data[1])

    μ = v[1:N_μ]
    p = SVector(v[n1+1:n1+n2]...)
    data = copy( reinterpret(SVector{n_x, Float64}, v[n1+n2+1:n1+n2+n3]))
    x = MSM(data, dmodel.F.x0.sizes)
    y = SVector(v[n1+n2+n3+1:end]...)

    u = Dolark.Unknown(μ, p, x, y)

    res, J = Residual(dmodel, u)
    r = flatten(res)

    if diff
        return r, J
    else
        return r
    end

end

import Dolo:MSM
import Base: *

function flatten(u::Dolark.Unknown)
    cat(
        u.μ[:],
        u.p,
        cat(u.x.data...; dims=1),
        u.y
    ;dims=1
    )
end

function unflatten(u::Dolark.Unknown, v::AbstractVector)

    n1 = length(u.μ)
    μ = reshape(v[1:n1], size(u.μ)...)
    n2 = length(u.p)
    p = SVector(v[n1+1:n1+n2]...)
    n_x = size(u.x.data[1],1)
    n3 = length(u.x.data)*n_x
    data = copy( reinterpret(SVector{n_x, Float64}, v[n1+n2+1:n1+n2+n3]))
    x = MSM(data, u.x.sizes)
    y = SVector(v[n1+n2+n3+1:end]...)

    u = Dolark.Unknown(μ, p, x, y)

    return u

end

using IterativeSolvers

using Plots

import Base: +

function +(a::Dolark.Unknown, b::Dolark.Unknown)
    p = a.p + b.p
    x = a.x + b.x
    y = a.y + b.y
    μ = a.μ + b.μ
    return Dolark.Unknown(μ,p,x,y)
end

using FiniteDiff

function proto_solve_steady_state(dmodel, u0; numdiff=true, use_blas=true, maxit=10)

    backsteps = [2.0^(-i) for i=0:10]

    local guess
    u0_ = flatten(u0)

    for i=1:maxit

        println(u0_[end])
        r_, J = Residual(dmodel, u0_; diff=true)

        ε = maximum(abs, r_)

        println("r: $(ε)")
        if numdiff
            J = FiniteDiff.finite_difference_jacobian(u->Dolark.Residual(dmodel, u; diff=false), u0_) # ;
        end

        if use_blas
            J = convert(Matrix, J)
            δ_ = -J \ r_
        else
            δ_ = -gmres(J,r_)
        end

        for λ in backsteps
            guess = u0_ + δ_*λ  # (J,λ*δ_)
            εg_ = [1.0]
            try
                rg = Residual(dmodel, guess)
                εg_[1] = maximum(abs,rg)
            catch
                εg_[1] = 1000.0
            end
            εg = εg_[1]
            if εg<ε
                println(εg)
                u0_ = guess
                break
            end
        end
        # δ = unflatten(J, δ_)

        # u0 = u0 + δ

    end

    return unflatten(u0, guess)
end


text:
# Autogenerated wrapper script for Libepoxy_jll for powerpc64le-linux-gnu
export libepoxy

using Libglvnd_jll
using Xorg_libX11_jll
JLLWrappers.@generate_wrapper_header("Libepoxy")
JLLWrappers.@declare_library_product(libepoxy, "libepoxy.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Xorg_libX11_jll)
    JLLWrappers.@init_library_product(
        libepoxy,
        "lib/libepoxy.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module NonlinearRHS
import FFTW
import Hankel
import Cubature
import Base: show
import LinearAlgebra: mul!, ldiv!
import NumericalIntegration: integrate, SimpsonEven
import Luna: PhysData, Modes, Maths, Grid
import Luna.PhysData: wlfreq

"""
    to_time!(Ato, Aω, Aωo, IFTplan)

Transform ``A(ω)`` on normal grid to ``A(t)`` on oversampled time grid.
"""
function to_time!(Ato::Array{<:Real, D}, Aω, Aωo, IFTplan) where D
    N = size(Aω, 1)
    No = size(Aωo, 1)
    scale = (No-1)/(N-1) # Scale factor makes up for difference in FFT array length
    fill!(Aωo, 0)
    copy_scale!(Aωo, Aω, N, scale)
    mul!(Ato, IFTplan, Aωo)
end

function to_time!(Ato::Array{<:Complex, D}, Aω, Aωo, IFTplan) where D
    N = size(Aω, 1)
    No = size(Aωo, 1)
    scale = No/N # Scale factor makes up for difference in FFT array length
    fill!(Aωo, 0)
    copy_scale_both!(Aωo, Aω, N÷2, scale)
    mul!(Ato, IFTplan, Aωo)
end

"""
    to_freq!(Aω, Aωo, Ato, FTplan)

Transform oversampled A(t) to A(ω) on normal grid
"""
function to_freq!(Aω, Aωo, Ato::Array{<:Real, D}, FTplan) where D
    N = size(Aω, 1)
    No = size(Aωo, 1)
    scale = (N-1)/(No-1) # Scale factor makes up for difference in FFT array length
    mul!(Aωo, FTplan, Ato)
    copy_scale!(Aω, Aωo, N, scale)
end

function to_freq!(Aω, Aωo, Ato::Array{<:Complex, D}, FTplan) where D
    N = size(Aω, 1)
    No = size(Aωo, 1)
    scale = N/No # Scale factor makes up for difference in FFT array length
    mul!(Aωo, FTplan, Ato)
    copy_scale_both!(Aω, Aωo, N÷2, scale)
end

"""
    copy_scale!(dest, source, N, scale)

Copy first N elements from source to dest and simultaneously multiply by scale factor.
For multi-dimensional `dest` and `source`, work along first axis.
"""
function copy_scale!(dest::Vector, source::Vector, N, scale)
    for i = 1:N
        dest[i] = scale * source[i]
    end
end

"""
    copy_scale_both!(dest::Vector, source::Vector, N, scale)

Copy first and last N elements from source to first and last N elements in dest
and simultaneously multiply by scale factor.
For multi-dimensional `dest` and `source`, work along first axis.
"""
function copy_scale_both!(dest::Vector, source::Vector, N, scale)
    for i = 1:N
        dest[i] = scale * source[i]
    end
    for i = 1:N
        dest[end-i+1] = scale * source[end-i+1]
    end
end

function copy_scale!(dest, source, N, scale)
    (size(dest)[2:end] == size(source)[2:end]
     || error("dest and source must be same size except along first dimension"))
    idcs = CartesianIndices(size(dest)[2:end])
    _cpsc_core(dest, source, N, scale, idcs)
end

function _cpsc_core(dest, source, N, scale, idcs)
    for i in idcs
        for j = 1:N
            dest[j, i] = scale * source[j, i]
        end
    end
end

function copy_scale_both!(dest, source, N, scale)
    (size(dest)[2:end] == size(source)[2:end]
     || error("dest and source must be same size except along first dimension"))
    idcs = CartesianIndices(size(dest)[2:end])
    _cpscb_core(dest, source, N, scale, idcs)
end

function _cpscb_core(dest, source, N, scale, idcs)
    for i in idcs
        for j = 1:N
            dest[j, i] = scale * source[j, i]
        end
        for j = 1:N
            dest[end-j+1, i] = scale * source[end-j+1, i]
        end
    end
end

"""
    Et_to_Pt!(Pt, Et, responses, density)

Accumulate responses induced by Et in Pt.
"""
function Et_to_Pt!(Pt, Et, responses, density::Number)
    fill!(Pt, 0)
    for resp! in responses
        resp!(Pt, Et, density)
    end
end

function Et_to_Pt!(Pt, Et, responses, density::AbstractVector)
    fill!(Pt, 0)
    for ii in eachindex(density)
        for resp! in responses[ii]
            resp!(Pt, Et, density[ii])
        end
    end
end

function Et_to_Pt!(Pt, Et, responses, density, idcs)
    for i in idcs
        Et_to_Pt!(view(Pt, :, i), view(Et, :, i), responses, density)
    end
end

"""
    TransModal

Transform E(ω) -> Pₙₗ(ω) for modal fields.
"""
mutable struct TransModal{tsT, lT, TT, FTT, rT, gT, dT, ddT, nT}
    ts::tsT
    full::Bool
    dimlimits::lT
    Emω::Array{ComplexF64,2}
    Erω::Array{ComplexF64,2}
    Erωo::Array{ComplexF64,2}
    Er::Array{TT,2}
    Pr::Array{TT,2}
    Prω::Array{ComplexF64,2}
    Prωo::Array{ComplexF64,2}
    Prmω::Array{ComplexF64,2}
    FT::FTT
    resp::rT
    grid::gT
    densityfun::dT
    density::ddT
    norm!::nT
    ncalls::Int
    z::Float64
    rtol::Float64
    atol::Float64
    mfcn::Int
end

function show(io::IO, t::TransModal)
    grid = "grid type: $(typeof(t.grid))"
    modes = "modes: $(t.ts.nmodes)\n"*" "^4*join([string(mi) for mi in t.ts.ms], "\n    ")
    p = t.ts.indices == 1:2 ? "x,y" : t.ts.indices == 1 ? "x" : "y"
    pol = "polarisation: $p"
    samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))"
    resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n    ")
    full = "full: $(t.full)"
    out = join(["TransModal", modes, pol, grid, samples, full, resp], "\n  ")
    print(io, out)
end

"""
    TransModal(grid, ts, FT, resp, densityfun, norm!; rtol=1e-3, atol=0.0, mfcn=300, full=false)

Construct a `TransModal`, transform E(ω) -> Pₙₗ(ω) for modal fields.

# Arguments
- `grid::AbstractGrid` : the grid used in the simulation
- `ts::Modes.ToSpace` : pre-created `ToSpace` for conversion from modal fields to space
- `FT::FFTW.Plan` : the time-frequency Fourier transform for the oversampled time grid
- `resp` : `Tuple` of response functions
- `densityfun` : callable which returns the gas density as a function of `z`
- `norm!` : normalisation function as fctn of `z`, can be created via [`norm_modal`](@ref)
- `rtol::Float=1e-3` : relative tolerance on the `HCubature` integration
- `atol::Float=0.0` : absolute tolerance on the `HCubature` integration
- `mfcn::Int=300` : maximum number of function evaluations for one modal integration
- `full::Bool=false` : if `true`, use full 2-D mode integral, if `false`, only do radial integral
"""
function TransModal(tT, grid, ts::Modes.ToSpace, FT, resp, densityfun, norm!;
                    rtol=1e-3, atol=0.0, mfcn=300, full=false)
    Emω = Array{ComplexF64,2}(undef, length(grid.ω), ts.nmodes)
    Erω = Array{ComplexF64,2}(undef, length(grid.ω), ts.npol)
    Erωo = Array{ComplexF64,2}(undef, length(grid.ωo), ts.npol)
    Er = Array{tT,2}(undef, length(grid.to), ts.npol)
    Pr = Array{tT,2}(undef, length(grid.to), ts.npol)
    Prω = Array{ComplexF64,2}(undef, length(grid.ω), ts.npol)
    Prωo = Array{ComplexF64,2}(undef, length(grid.ωo), ts.npol)
    Prmω = Array{ComplexF64,2}(undef, length(grid.ω), ts.nmodes)
    IFT = inv(FT)
    TransModal(ts, full, Modes.dimlimits(ts.ms[1]), Emω, Erω, Erωo, Er, Pr, Prω, Prωo, Prmω,
               FT, resp, grid, densityfun, densityfun(0.0), norm!, 0, 0.0, rtol, atol, mfcn)
end

function TransModal(grid::Grid.RealGrid, args...; kwargs...)
    TransModal(Float64, grid, args...; kwargs...)
end

function TransModal(grid::Grid.EnvGrid, args...; kwargs...)
    TransModal(ComplexF64, grid, args...; kwargs...)
end

function reset!(t::TransModal, Emω::Array{ComplexF64,2}, z::Float64)
    t.Emω .= Emω
    t.ncalls = 0
    t.z = z
    t.dimlimits = Modes.dimlimits(t.ts.ms[1], z=z)
    t.density = t.densityfun(z)
end

function pointcalc!(fval, xs, t::TransModal)
    # TODO: parallelize this in Julia 1.3
    for i in 1:size(xs, 2)
        x1 = xs[1, i]
        # on or outside boundaries are zero
        if x1 <= t.dimlimits[2][1] || x1 >= t.dimlimits[3][1]
            fval[:, i] .= 0.0
            continue
        end
        if size(xs, 1) > 1 # full 2-D mode integral
            x2 = xs[2, i]
            if t.dimlimits[1] == :polar
                pre = x1
            else
                if x2 <= t.dimlimits[2][2] || x1 >= t.dimlimits[3][2]
                    fval[:, i] .= 0.0
                    continue
                end
                pre = 1.0
            end
        else
            if t.dimlimits[1] == :polar
                x2 = 0.0
                pre = 2π*x1
            else
                x2 = 0.0
                pre = 1.0
            end
        end
        x = (x1,x2)
        Modes.to_space!(t.Erω, t.Emω, x, t.ts, z=t.z)
        to_time!(t.Er, t.Erω, t.Erωo, inv(t.FT))
        # get nonlinear pol at r,θ
        Et_to_Pt!(t.Pr, t.Er, t.resp, t.density)
        t.ncalls += 1
        @. t.Pr *= t.grid.towin
        to_freq!(t.Prω, t.Prωo, t.Pr, t.FT)
        @. t.Prω *= t.grid.ωwin
        t.norm!(t.Prω)
        # now project back to each mode
        # matrix product (nω x npol) * (npol x nmodes) -> (nω x nmodes)
        mul!(t.Prmω, t.Prω, transpose(t.ts.Ems))
        fval[:, i] .= pre.*reshape(reinterpret(Float64, t.Prmω), length(t.Emω)*2)
    end
end

function (t::TransModal)(nl, Eω, z)
    reset!(t, Eω, z)
    _, ll, ul = t.dimlimits
    if t.full
        val, err = Cubature.hcubature_v(
            length(Eω)*2,
            (x, fval) -> pointcalc!(fval, x, t),
            ll, ul,
            reltol=t.rtol, abstol=t.atol, maxevals=t.mfcn, error_norm=Cubature.L2)
    else
        val, err = Cubature.pcubature_v(
            length(Eω)*2,
            (x, fval) -> pointcalc!(fval, x, t),
            (ll[1],), (ul[1],),
            reltol=t.rtol, abstol=t.atol, maxevals=t.mfcn, error_norm=Cubature.L2)
    end
    nl .= reshape(reinterpret(ComplexF64, val), size(nl))
end

"""
    norm_modal(grid; shock=true)

Normalisation function for modal propagation. If `shock` is `false`, the intrinsic frequency
dependence of the nonlinear response is ignored, which turns off optical shock formation/
self-steepening.
"""
function norm_modal(grid; shock=true)
    ω0 = PhysData.wlfreq(grid.referenceλ)
    withshock!(nl) = @. nl *= (-im * grid.ω/4)
    withoutshock!(nl) = @. nl *= (-im * ω0/4)
    shock ? withshock! : withoutshock!
end

"""
    TransModeAvg

Transform E(ω) -> Pₙₗ(ω) for mode-averaged single-mode propagation.
"""
struct TransModeAvg{TT, FTT, rT, gT, dT, nT, aT}
    Pto::Vector{TT}
    Eto::Vector{TT}
    Eωo::Vector{ComplexF64}
    Pωo::Vector{ComplexF64}
    FT::FTT
    resp::rT
    grid::gT
    densityfun::dT
    norm!::nT
    aeff::aT # function which returns effective area
end

function show(io::IO, t::TransModeAvg)
    grid = "grid type: $(typeof(t.grid))"
    samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))"
    resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n    ")
    out = join(["TransModeAvg", grid, samples, resp], "\n  ")
    print(io, out)
end

function TransModeAvg(TT, grid, FT, resp, densityfun, norm!, aeff)
    Eωo = zeros(ComplexF64, length(grid.ωo))
    Eto = zeros(TT, length(grid.to))
    Pto = similar(Eto)
    Pωo = similar(Eωo)
    TransModeAvg(Pto, Eto, Eωo, Pωo, FT, resp, grid, densityfun, norm!, aeff)
end

function TransModeAvg(grid::Grid.RealGrid, FT, resp, densityfun, norm!, aeff)
    TransModeAvg(Float64, grid, FT, resp, densityfun, norm!, aeff)
end

function TransModeAvg(grid::Grid.EnvGrid, FT, resp, densityfun, norm!, aeff)
    TransModeAvg(ComplexF64, grid, FT, resp, densityfun, norm!, aeff)
end

const nlscale = sqrt(PhysData.ε_0*PhysData.c/2)

function (t::TransModeAvg)(nl, Eω, z)
    to_time!(t.Eto, Eω, t.Eωo, inv(t.FT))
    @. t.Eto /= nlscale*sqrt(t.aeff(z))
    Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z))
    @. t.Pto *= t.grid.towin
    to_freq!(nl, t.Pωo, t.Pto, t.FT)
    t.norm!(nl, z)
    for i in eachindex(nl)
        !t.grid.sidx[i] && continue
        nl[i] *= t.grid.ωwin[i]
    end
end

function norm_mode_average(grid, βfun!, aeff; shock=true)
    β = zeros(Float64, length(grid.ω))
    shockterm = shock ? grid.ω.^2 : grid.ω .* PhysData.wlfreq(grid.referenceλ)
    pre = @. -im*shockterm/(2*PhysData.c^(3/2)*sqrt(2*PhysData.ε_0))
    function norm!(nl, z)
        βfun!(β, z)
        sqrtaeff = sqrt(aeff(z))
        for i in eachindex(nl)
            !grid.sidx[i] && continue
            nl[i] *= pre[i]*sqrtaeff/β[i]
        end
    end
end

"""
    TransRadial

Transform E(ω) -> Pₙₗ(ω) for radially symetric free-space propagation
"""
struct TransRadial{TT, HTT, FTT, nT, rT, gT, dT, iT}
    QDHT::HTT # Hankel transform (space to k-space)
    FT::FTT # Fourier transform (time to frequency)
    normfun::nT # Function which returns normalisation factor
    resp::rT # nonlinear responses (tuple of callables)
    grid::gT # time grid
    densityfun::dT # callable which returns density
    Pto::Array{TT,2} # Buffer array for NL polarisation on oversampled time grid
    Eto::Array{TT,2} # Buffer array for field on oversampled time grid
    Eωo::Array{ComplexF64,2} # Buffer array for field on oversampled frequency grid
    Pωo::Array{ComplexF64,2} # Buffer array for NL polarisation on oversampled frequency grid
    idcs::iT # CartesianIndices for Et_to_Pt! to iterate over
end

function show(io::IO, t::TransRadial)
    grid = "grid type: $(typeof(t.grid))"
    samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))"
    resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n    ")
    nr = "radial points: $(t.QDHT.N)"
    R = "aperture: $(t.QDHT.R)"
    out = join(["TransRadial", grid, samples, nr, R, resp], "\n  ")
    print(io, out)
end

function TransRadial(TT, grid, HT, FT, responses, densityfun, normfun)
    Eωo = zeros(ComplexF64, (length(grid.ωo), HT.N))
    Eto = zeros(TT, (length(grid.to), HT.N))
    Pto = similar(Eto)
    Pωo = similar(Eωo)
    idcs = CartesianIndices(size(Pto)[2:end])
    TransRadial(HT, FT, normfun, responses, grid, densityfun, Pto, Eto, Eωo, Pωo, idcs)
end

"""
    TransRadial(grid, HT, FT, responses, densityfun, normfun)

Construct a `TransRadial` to calculate the reciprocal-domain nonlinear polarisation.

# Arguments
- `grid::AbstractGrid` : the grid used in the simulation
- `HT::QDHT` : the Hankel transform which defines the spatial grid
- `FT::FFTW.Plan` : the time-frequency Fourier transform for the oversampled time grid
- `responses` : `Tuple` of response functions
- `densityfun` : callable which returns the gas density as a function of `z`
- `normfun` : normalisation factor as fctn of `z`, can be created via [`norm_radial`](@ref)
"""
function TransRadial(grid::Grid.RealGrid, args...)
    TransRadial(Float64, grid, args...)
end

function TransRadial(grid::Grid.EnvGrid, args...)
    TransRadial(ComplexF64, grid, args...)
end

"""
    (t::TransRadial)(nl, Eω, z)

Calculate the reciprocal-domain (ω-k-space) nonlinear response due to the field `Eω` and
place the result in `nl`
"""
function (t::TransRadial)(nl, Eω, z)
    to_time!(t.Eto, Eω, t.Eωo, inv(t.FT)) # transform ω -> t
    ldiv!(t.Eto, t.QDHT, t.Eto) # transform k -> r
    Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs) # add up responses
    @. t.Pto *= t.grid.towin # apodisation
    mul!(t.Pto, t.QDHT, t.Pto) # transform r -> k
    to_freq!(nl, t.Pωo, t.Pto, t.FT) # transform t -> ω
    nl .*= t.grid.ωwin .* (-im.*t.grid.ω)./(2 .* t.normfun(z))
end

"""
    const_norm_radial(ω, q, nfun)

Make function to return normalisation factor for radial symmetry without re-calculating at
every step.
"""
function const_norm_radial(grid, q, nfun)
    nfunω = (ω; z) -> nfun(wlfreq(ω))
    normfun = norm_radial(grid, q, nfunω)
    out = copy(normfun(0.0))
    function norm(z)
        return out
    end
    return norm
end

"""
    norm_radial(ω, q, nfun)

Make function to return normalisation factor for radial symmetry.

!!! note
    Here, `nfun(ω; z)` needs to take frequency `ω` and a keyword argument `z`.
"""
function norm_radial(grid, q, nfun)
    ω = grid.ω
    out = zeros(Float64, (length(ω), q.N))
    kr2 = q.k.^2
    k2 = zeros(Float64, length(ω))
    function norm(z)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z).*grid.ω[grid.sidx]./PhysData.c).^2
        for ir = 1:q.N
            for iω in eachindex(ω)
                if ω[iω] == 0
                    out[iω, ir] = 1.0
                    continue
                end
                βsq = k2[iω] - kr2[ir]
                if βsq <= 0
                    out[iω, ir] = 1.0
                    continue
                end
                out[iω, ir] = sqrt(βsq)/(PhysData.μ_0*ω[iω])
            end
        end
        return out
    end
    return norm
end

mutable struct TransFree{TT, FTT, nT, rT, gT, xygT, dT, iT}
    FT::FTT # 3D Fourier transform (space to k-space and time to frequency)
    normfun::nT # Function which returns normalisation factor
    resp::rT # nonlinear responses (tuple of callables)
    grid::gT # time grid
    xygrid::xygT
    densityfun::dT # callable which returns density
    Pto::Array{TT, 3} # buffer for oversampled time-domain NL polarisation
    Eto::Array{TT, 3} # buffer for oversampled time-domain field
    Eωo::Array{ComplexF64, 3} # buffer for oversampled frequency-domain field
    Pωo::Array{ComplexF64, 3} # buffer for oversampled frequency-domain NL polarisation
    scale::Float64 # scale factor to be applied during oversampling
    idcs::iT # iterating over these slices Eto/Pto into Vectors, one at each position
end

function show(io::IO, t::TransFree)
    grid = "grid type: $(typeof(t.grid))"
    samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))"
    resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n    ")
    y = "y grid: $(minimum(t.xygrid.y)) to $(maximum(t.xygrid.y)), N=$(length(t.xygrid.y))"
    x = "x grid: $(minimum(t.xygrid.x)) to $(maximum(t.xygrid.x)), N=$(length(t.xygrid.x))"
    out = join(["TransFree", grid, samples, y, x, resp], "\n  ")
    print(io, out)
end

function TransFree(TT, scale, grid, xygrid, FT, responses, densityfun, normfun)
    Ny = length(xygrid.y)
    Nx = length(xygrid.x)
    Eωo = zeros(ComplexF64, (length(grid.ωo), Ny, Nx))
    Eto = zeros(TT, (length(grid.to), Ny, Nx))
    Pto = similar(Eto)
    Pωo = similar(Eωo)
    idcs = CartesianIndices((Ny, Nx))
    TransFree(FT, normfun, responses, grid, xygrid, densityfun,
              Pto, Eto, Eωo, Pωo, scale, idcs)
end

"""
    TransFree(grid, xygrid, FT, responses, densityfun, normfun)

Construct a `TransFree` to calculate the reciprocal-domain nonlinear polarisation.

# Arguments
- `grid::AbstractGrid` : the grid used in the simulation
- `xygrid` : the spatial grid (instances of [`Grid.FreeGrid`](@ref))
- `FT::FFTW.Plan` : the full 3D (t-y-x) Fourier transform for the oversampled time grid
- `responses` : `Tuple` of response functions
- `densityfun` : callable which returns the gas density as a function of `z`
- `normfun` : normalisation factor as fctn of `z`, can be created via [`norm_free`](@ref)
"""
function TransFree(grid::Grid.RealGrid, args...)
    N = length(grid.ω)
    No = length(grid.ωo)
    scale = (No-1)/(N-1)
    TransFree(Float64, scale, grid, args...)
end

function TransFree(grid::Grid.EnvGrid, args...)
    N = length(grid.ω)
    No = length(grid.ωo)
    scale = No/N
    TransFree(ComplexF64, scale, grid, args...)
end

"""
    (t::TransFree)(nl, Eω, z)

Calculate the reciprocal-domain (ω-kx-ky-space) nonlinear response due to the field `Eω`
and place the result in `nl`.
"""
function (t::TransFree)(nl, Eωk, z)
    fill!(t.Eωo, 0)
    copy_scale!(t.Eωo, Eωk, length(t.grid.ω), t.scale)
    ldiv!(t.Eto, t.FT, t.Eωo) # transform (ω, ky, kx) -> (t, y, x)
    Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs) # add up responses
    @. t.Pto *= t.grid.towin # apodisation
    mul!(t.Pωo, t.FT, t.Pto) # transform (t, y, x) -> (ω, ky, kx)
    copy_scale!(nl, t.Pωo, length(t.grid.ω), 1/t.scale)
    nl .*= t.grid.ωwin .* (-im.*t.grid.ω)./(2 .* t.normfun(z))
end

"""
    const_norm_free(grid, xygrid, nfun)

Make function to return normalisation factor for 3D propagation without re-calculating at
every step.
"""
function const_norm_free(grid, xygrid, nfun)
    nfunω = (ω; z) -> nfun(wlfreq(ω))
    normfun = norm_free(grid, xygrid, nfunω)
    out = copy(normfun(0.0))
    function norm(z)
        return out
    end
    return norm
end

"""
    norm_free(grid, xygrid, nfun)

Make function to return normalisation factor for 3D propagation.

!!! note
    Here, `nfun(ω; z)` needs to take frequency `ω` and a keyword argument `z`.
"""
function norm_free(grid, xygrid, nfun)
    ω = grid.ω
    kperp2 = @. (xygrid.kx^2)' + xygrid.ky^2
    idcs = CartesianIndices((length(xygrid.ky), length(xygrid.kx)))
    k2 = zero(grid.ω)
    out = zeros(Float64, (length(grid.ω), length(xygrid.ky), length(xygrid.kx)))
    function norm(z)
        k2[grid.sidx] = (nfun.(grid.ω[grid.sidx]; z=z).*grid.ω[grid.sidx]./PhysData.c).^2
        for ii in idcs
            for iω in eachindex(ω)
                if ω[iω] == 0
                    out[iω, ii] = 1.0
                    continue
                end
                βsq = k2[iω] - kperp2[ii]
                if βsq <= 0
                    out[iω, ii] = 1.0
                    continue
                end
                out[iω, ii] = sqrt(βsq)/(PhysData.μ_0*ω[iω])
            end
        end
        return out
    end
end

end


text:
using TiledIteration

### Construction of loop expressions in type domain
### This is the fallback implementation for AbstractArrays
function arrayop_body{A<:AbstractArray, L,R,F,E}(name, ::Type{A},
                                                 op::Type{Assign{L,R,F,E}})

    acc = kernel_expr(:($name.lhs), L) # :() will be ignored
    rhs_inner = kernel_expr(:($name.rhs), R)

    if hasreduceddims(op)
        # we need to wrap the expression in a call to the reducer
        expr = :($acc = $name.reducefn($acc, $rhs_inner))
    else
        expr = :($acc = $rhs_inner)
    end

    checks = :()
    input_ranges = Any[]

    simded = false

    for (sym, spaces) in index_spaces(name, op) # sort this based on # of potential cache misses
        if length(spaces) > 1
            # check dimensions for equality
            equal_dims = [:(indices($(d[3]), $(d[2]))) for d in spaces]
            checks = :($checks; @assert allequal($(equal_dims...),))
        end
        T,dim,nm = first(spaces)
        sym_range = Symbol("$(sym)_range")
        push!(input_ranges, sym_range => :(indices($nm, $dim)))
        expr = :(for $sym in $sym_range
                    $expr
                end)
        if !simded
            # innermost loop
            expr = :(@simd $expr)
            simded = true
        end
    end

    expr = quote
        ranges = ($(map(last, input_ranges)...),)
        @inbounds for tile in TileIterator(ranges, tilesize(ranges))
            ($(map(first, input_ranges)...),) = tile
            $expr
        end
    end

    :($checks; $expr; $name.lhs.array)
end

# constructs Dict
# Keys:  indexing symbol
# Value: a list of (array type, dimension, expr) tuples which correspond to that
#        indexing symbol
# Generate the expression corresponding to a type
function kernel_expr(name, itr)
    kernel_expr(name, arraytype(itr), itr)
end

function kernel_expr{X, Idx, A<:AbstractArray}(name, ::Type{A},
                                               itr::Type{Indexing{X, Idx}})
    idx = get_subscripts(name, itr)
    :($name.array[$(idx...)])
end

function kernel_expr{A<:AbstractArray, F, Ts}(name, ::Type{A},
                                              itr::Type{Map{F, Ts}})

    inner_kernels = [kernel_expr(:($name.arrays[$i]), arraytype(T), T)
                        for (i, T) in enumerate(Ts.parameters)]

    :($name.f($(inner_kernels...)))
end

allequal(x) = true
function allequal(x, xs...)
    x == xs[1] && allequal(xs...)
end

function tilesize(ranges)
    map(x->16, ranges)
end

@inline @generated function arrayop!{L,R,A<:AbstractArray}(::Type{A}, t::Assign{L,R})
    arrayop_body(:t, arraytype(L), t)
end


text:
"Fix to ensure positive definiteness by dividing each off-diagonal element by sum of absolute values of off-diagonal elements in its row"
function fixMatrix(A::Array{Float64}, denom_factor::Float64)
	p = size(A, 1);
	for cur_row in 1:p
		cur_sum = sum(abs(A[cur_row, :])) - 1
		if cur_sum != 1
			A[cur_row, :] = A[cur_row, :] / (denom_factor * cur_sum);
		end
		# Make sure diagonal entries are still 1
		A[cur_row, cur_row] = 1;
	end
	# Final matrix is average of matrix with its transpose
	A = (A + A')/2
    return A
end

"Calculate Covariance Matrix using Gaussian Kernel"
function calcSigma(X, inv_Sigma)
	Sigma = zeros(size(X, 1), size(X, 1))
	for i in 1:size(Sigma, 1)
	  for j in 1:size(Sigma, 2)
          Sigma[i, j] = exp(- 1/2 * (X[i, :] - X[j, :])' * inv_Sigma * (X[i, :] - X[j, :]))[1]
	  end
	end
	return Sigma
end


text:
module TupleVectors

using Requires: @require

include("tuplevector.jl")
include("summarize.jl")
include("chainvec.jl")
include("with.jl")
include("hypercube.jl")

function __init__()
    @require GoldenSequences = "e13314f6-e9b3-11e9-0763-f13dff015e8a" begin
        include("optional/GoldenSequences.jl")
    end

    @require MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" begin
        include("optional/MonteCarloMeasurements.jl")
    end

    @require Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" begin
        include("optional/Sobol.jl")
    end

    @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" begin
        include("optional/Distributions.jl")
    end

    @require MeasureTheory = "eadaa1a4-d27c-401d-8699-e962e1bbc33b" begin
        include("optional/MeasureTheory.jl")
    end

end

end


text:
struct DSA3200Status
    ptype::Int16
    pad::NTuple{78,UInt8}
    status::NTuple{20,UInt8}
    pad2::NTuple{80,UInt8}
end

struct DSA3200ScanRaw
    ptype::Int16
    pad::Int16
    frame::Int32
    press::NTuple{16,Int16}
    temp::NTuple{16,Int16}
end

struct DSA3200ScanEU
    ptype::Int16
    pad::Int16
    frame::Int32
    press::NTuple{16,Float32}
    temp::NTuple{16,Int16}
end

struct DSA3200ScanRawTime
    ptype::Int16
    pad::Int16
    frame::Int32
    press::NTuple{16,Int16}
    temp::NTuple{16,Int16}
    time::Int32
    timeunit::Int32
end

struct DSA3200ScanEUTime
    ptype::Int16
    pad::Int16
    frame::Int32
    press::NTuple{16,Float32}
    temp::NTuple{16,Int16}
    time::Int32
    timeunit::Int32
end


text:
using TropicalTensors
using TropicalTensors: getbonds, getvertices, rand_Js, rand_hs

function run(L)
    lt = CubicLattice(L, L, L)
    sg = Spinglass(lt, rand_Js(Float16, Randpm(), lt), rand_hs(Float16, Zero(), lt))
    res = @time solve(Tropical{Float16}, sg; usecuda=false)
end

const L = parse(Int, ARGS[1])
res = run(L)
println("The maximum energy is $(res.n).")


text: # put in processing styled copies of the GLFW keycodes


text:
function getdist(v1::Array{Float64,1},v2::Array{Float64,1},geometry="planar")
	if geometry=="planar"
		return norm(v1-v2)
	elseif geometry=="sphere"
		if (v1[1]-v2[1])*(v1[1]-v2[1])+
			(v1[2]-v2[2])*(v1[2]-v2[2])+
			(v1[3]-v2[3])*(v1[3]-v2[3])<1e-5
			#if norm(v1-v2)<1e-15
			return norm(v1-v2)
		end
		#https://en.wikipedia.org/wiki/Great-circle_distance
		#return atan(norm(cross(v1,v2),2)/dot(v1,v2))
		return acos(max(-1.0,v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3]))
	else
		error("Type has to be planar or sphere.")
	end
end

function getangle(B::Array{Float64,1},A::Array{Float64,1},C::Array{Float64,1}, geometry="planar")
	if geometry=="planar"
		u,v = A-B, C-A;
		return acos(dot(u,v)/norm(u,2)/norm(v,2));
	elseif geometry=="sphere"
		c= getdist(B,A,"sphere");
		b =getdist(A,C,"sphere");
		a= getdist(C,B,"sphere");
		if c<1e-10 || b<1e-10 || a<1e-10
			return 0;
		end
		#https://en.wikipedia.org/wiki/Spherical_trigonometry#Cosine_rules_and_sine_rules
		tmp = (cos(a)-cos(b)*cos(c))/(sin(b)*sin(c));
		# It might happen, that due to computations, tmp is not inside [-1,1]
		# but inside [-1-eps,1+eps]
		tmp = max(-1.0,min(1.0,tmp))

		angle = acos(tmp);
		return angle;
	else
		error("Type has to be planar or sphere.")
	end
end

function getarea(A::Array{Float64,1},B::Array{Float64,1},C::Array{Float64,1}, geometry="planar")
	if geometry=="planar"
		alpha = getangle(B,A,C)
		lb = norm(B-A)
		la = norm(C-A)
		return 0.5*sin(alpha)*lb*la
	elseif geometry=="sphere"
		alpha = getangle(B,A,C,"sphere")
		beta =  getangle(C,B,A,"sphere")
		gamma = getangle(A,C,B,"sphere")
		#https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess
		#Excess = max(0.0,alpha+beta+gamma-pi)
		#R = 1; # We always consider the uni sphere
		#return Excess*R*R;
		return alpha+beta+gamma-pi
	else
		error("Type has to be planar or sphere.")
	end
end

function muphi2xyz!(muphi::Array{Float64,2},xyz::Array{Float64,2})
	n = size(xyz)[1]
	for i=1:n
		xyz[i,1] = sqrt(1-muphi[i,1]^2)*cos(muphi[i,2])
		xyz[i,2] = sqrt(1-muphi[i,1]^2)*sin(muphi[i,2])
		xyz[i,3] =        muphi[i,1]
	end
end

function xyz2muphi!(xyz::Array{Float64,2},muphi::Array{Float64,2})
	n = size(xyz)[1]
	for i=1:n
		muphi[i,1] = xyz[i,3]
		muphi[i,2] =  atan(xyz[i,2],xyz[i,1])
	end
end

#Check if the projection of p onto the plane defined by p1 p2 p3 is inside the triangle p1 p2 p3
#https://math.stackexchange.com/questions/544946/determine-if-projection-of-3d-point-onto-plane-is-within-a-triangle
function isprojectedpointinsidetriangle(p::Array{Float64,1},p1::Array{Float64,1},
	p2::Array{Float64,1} ,p3::Array{Float64,1})
	n1 = cross(p1,p2);
	if n1'p3<0
		if p'n1>0
			return false
		end
	elseif n1'p3>0
		if p'n1<0
			return false
		end
	end

	n2 = cross(p2,p3);
	if n2'p1<0
		if p'n2>0
			return false
		end
	elseif n2'p1>0
		if p'n2<0
			return false
		end
	end

	n3 = cross(p3,p1);
	if n3'p2<0
		if p'n3>0
			return false
		end
	elseif n3'p2>0
		if p'n3<0
			return false
		end
	end
	return true
end

#Compute the rotation that rotates theta around
#the vector u
function getrotationmatrix!(Rot::Array{Float64,2},u::Array{Float64,1},theta::Float64,order="forward")
	#https://en.wikipedia.org/wiki/Rotation_matrix
	ux = u[1];
	uy = u[2];
	uz = u[3];
	costheta = cos(theta);
	sintheta = sin(theta);

	Rot[1,1] = ux*ux*(1-costheta)+   costheta;
	Rot[1,2] = ux*uy*(1-costheta)-uz*sintheta;
	Rot[1,3] = ux*uz*(1-costheta)+uy*sintheta;

	Rot[2,1] = uy*ux*(1-costheta)+uz*sintheta;
	Rot[2,2] = uy*uy*(1-costheta)+   costheta;
	Rot[2,3] = uy*uz*(1-costheta)-ux*sintheta;

	Rot[3,1] = uz*ux*(1-costheta)-uy*sintheta;
	Rot[3,2] = uz*uy*(1-costheta)+ux*sintheta;
	Rot[3,3] = uz*uz*(1-costheta)+   costheta;
	#niceprint(Rot)
	#if (abs(norm(Rot)-1)>1e-6)
	#	error("No rotation matrix.")
	#end
end

#
#Compute the rotation that rotas alpha around the x axis, beta around the y axis
#and gamma around the z axis.
#
function getrotationmatrix!(Rot::Array{Float64,2},alpha::Float64,beta::Float64,gamma::Float64,order="forward")
	RotX = zeros(3,3)
	RotY = zeros(3,3)
	RotZ = zeros(3,3)
	# set nonzero entries

	RotX[2,2] =  cos(alpha);
	RotX[2,3] = -sin(alpha);
	RotX[3,2] =  sin(alpha);
	RotX[3,3] =  cos(alpha);

	RotY[1,1] =  cos(beta);
	RotY[1,3] =  sin(beta);
	RotY[3,1] = -sin(beta);
	RotY[3,3] =  cos(beta);

	RotZ[1,1] =  cos(gamma);
	RotZ[1,2] = -sin(gamma);
	RotZ[2,1] =  sin(gamma);
	RotZ[2,2] =  cos(gamma);

	if order=="forward"
		Rot = deepcopy(RotZ*RotY*RotX);
	elseif order=="reverse"
		Rot = deepcopy(RotX*RotY*RotY);
	else
		error("Order has to be forward or reverse!")
	end
end

function computeinterpolationweights(v::Array{Float64,1},va::Array{Float64,1},vb::Array{Float64,1},
									 vc::Array{Float64,1})
	#https://classes.soe.ucsc.edu/cmps160/Fall10/resources/barycentricInterpolation.pdf

	weights = zeros(3)

	areac = getarea(v,va,vb,"sphere");
	areab = getarea(v,va,vc,"sphere");
	areaa = getarea(v,vb,vc,"sphere");

	area = areaa+areab+areac
	weights[1] = areaa/area;
	weights[2] = areab/area;
	weights[3] = areac/area;

	return weights

end

function interpolate(v::Array{Float64,1},va::Array{Float64,1},vb::Array{Float64,1},
					 vc::Array{Float64,1},phia::Float64,phib::Float64,phic::Float64)
	w = computeinterpolationweights(v,va,vb,vc)
	return phia*w[1]+phib*w[2]+phic*w[3]
end

function slerp(pt1::Array{Float64,1},pt2::Array{Float64,1},n::Int64,unevenfix=false::Bool)
	if norm(pt1-pt2)<1e-10 # same points
		return repeat(pt1,1,n) # return n copies of that point
	end

	omega =acos(dot(pt1,pt2));
	t = collect(range(0,stop=1,length=n))

	return (sin.((1 .-t)*omega)/sin.(omega))'.*pt1+
	(sin.(t*omega)/sin.(omega))'.*pt2
end

function unique(Points::Array{Float64,2},Triangles::Array{Int64,2})
	# Given points P and a triangulation T, this function removes
	# all duplicate (non unique) entries from P and adjusts the entries of T accordingly.
	# The final points and triangulation are stored in
	# XYZ and Triangulation (i.e.) the class attributes

	nPoints = size(Points)[2]
	nTriangles = size(Triangles)[2]

	map    = collect(1:nPoints)
	unique = ones(Int64,nPoints)

	#Okay so what is map?
	#map is a vector, where each entry contains the mapping to the ID in the unique vector we want to have.
	#We start with map(i) = i and then iterate through all the points
	#if it turns out, that point 10 == point 4, we will set map(10) = 4 (and map(4)=4)
	#if then also point 20 == point 10 == point 4, we will set map(10) == 4 as well.
	#unique simply stores a boolean as for whether or not the value is unique or a duplicate.
	#This does however keep the first occurrence as unique.
	#Meaning for a vector
	#[1 2 3 3 4 5 3]
	#unique would be
	#[1 1 1 0 1 1 0]
	#
	for i=1:nPoints
		for j=i+1:nPoints
			# apparently calling norm(p1-p2) allocates a lot of
			dist = sqrt(
						(Points[1,i]-Points[1,j])^2+
						(Points[2,i]-Points[2,j])^2+
						(Points[3,i]-Points[3,j])^2)
			if (dist<1e-6) # Points are equal
				map[j] = min(map[j],i); # Take the minimal ID
				unique[j] = 0; # This point is no longer unique
			end
		end
	end

	#We are still not done
	#Consider again the vector from above
	#[1 2 3 3 4 5 3]
	#map was at first
	#[0 1 2 3 4 5 6]
	#and is now
	#[0 1 2 2 4 5 2]
	#meaning that we have to renumber certain entries because
	#we want to get
	#[0 1 2 2 3 4 2] , i.e. the "real count
	#
	uniqueCounter = 1 # The "real" count
	for i=1:nPoints
		if unique[i]==1
			map[i] = uniqueCounter;
			uniqueCounter +=1;
		else
			map[i] = map[map[i]]
		end
	end

	# Now we go through the triangles and update the IDs
	triangulation = zeros(Int64,3,nTriangles)
	for i=1:nTriangles
		for j=1:3
			idx = Triangles[j,i];
			triangulation[j,i] = min(idx,map[idx])
		end
	end

	#Count unique points
	nQuadPoints=0;
	for i=1:nPoints
		nQuadPoints += (unique[i]==1);
	end

	#Write unique points to XYZ
	xyz = zeros(Float64,3,nQuadPoints);
	counter = 1;
	for i=1:nPoints
		if unique[i]==1
			xyz[:,counter] = Points[:,i]
			counter +=1;
		end
	end
	return xyz,triangulation
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using ITensors,
      Test,
      Random

include("util.jl")

function isingMPO(sites)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N+1)
  for n=1:N+1
    link[n] = Index(3,"Link,Ising,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(dag(ll),dag(s),s',rl)
    H[n] += setelt(ll[1])*setelt(rl[1])*op(sites,"Id",n)
    H[n] += setelt(ll[3])*setelt(rl[3])*op(sites,"Id",n)
    H[n] += setelt(ll[2])*setelt(rl[1])*op(sites,"Sz",n)
    H[n] += setelt(ll[3])*setelt(rl[2])*op(sites,"Sz",n)
  end
  LE = ITensor(link[1]); LE[3] = 1.0;
  RE = ITensor(dag(link[N+1])); RE[1] = 1.0;
  H[1] *= LE
  H[N] *= RE
  return H
end

function heisenbergMPO(sites,
                       h::Vector{Float64},
                       onsite::String="Sz")::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N+1)
  for n=1:N+1
    link[n] = Index(5,"Link,Heis,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setelt(ll[1])*setelt(rl[1])*op(sites,"Id",n)
    H[n] += setelt(ll[5])*setelt(rl[5])*op(sites,"Id",n)
    H[n] += setelt(ll[2])*setelt(rl[1])*op(sites,"S+",n)
    H[n] += setelt(ll[3])*setelt(rl[1])*op(sites,"S-",n)
    H[n] += setelt(ll[4])*setelt(rl[1])*op(sites,"Sz",n)
    H[n] += setelt(ll[5])*setelt(rl[2])*op(sites,"S-",n)*0.5
    H[n] += setelt(ll[5])*setelt(rl[3])*op(sites,"S+",n)*0.5
    H[n] += setelt(ll[5])*setelt(rl[4])*op(sites,"Sz",n)
    H[n] += setelt(ll[5])*setelt(rl[1])*op(sites,onsite,n)*h[n]
  end
  H[1] *= setelt(link[1][5])
  H[N] *= setelt(link[N+1][1])
  return H
end

function NNheisenbergMPO(sites,
                         J1::Float64,
                         J2::Float64)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N+1)
  for n=1:N+1
    link[n] = Index(8,"Link,H,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setelt(ll[1])*setelt(rl[1])*op(sites,"Id",n)
    H[n] += setelt(ll[8])*setelt(rl[8])*op(sites,"Id",n)

    H[n] += setelt(ll[2])*setelt(rl[1])*op(sites,"S-",n)
    H[n] += setelt(ll[5])*setelt(rl[2])*op(sites,"Id",n)
    H[n] += setelt(ll[8])*setelt(rl[2])*op(sites,"S+",n)*J1/2
    H[n] += setelt(ll[8])*setelt(rl[5])*op(sites,"S+",n)*J2/2

    H[n] += setelt(ll[3])*setelt(rl[1])*op(sites,"S+",n)
    H[n] += setelt(ll[6])*setelt(rl[3])*op(sites,"Id",n)
    H[n] += setelt(ll[8])*setelt(rl[3])*op(sites,"S-",n)*J1/2
    H[n] += setelt(ll[8])*setelt(rl[6])*op(sites,"S-",n)*J2/2

    H[n] += setelt(ll[4])*setelt(rl[1])*op(sites,"Sz",n)
    H[n] += setelt(ll[7])*setelt(rl[4])*op(sites,"Id",n)
    H[n] += setelt(ll[8])*setelt(rl[4])*op(sites,"Sz",n)*J1
    H[n] += setelt(ll[8])*setelt(rl[7])*op(sites,"Sz",n)*J2
  end
  H[1] *= setelt(link[1][8])
  H[N] *= setelt(link[N+1][1])
  return H
end

function threeSiteIsingMPO(sites,
                           h::Vector{Float64})::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N+1)
  for n=1:N+1
    link[n] = Index(4,"Link,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setelt(ll[1])*setelt(rl[1])*op(sites,"Id",n)
    H[n] += setelt(ll[4])*setelt(rl[4])*op(sites,"Id",n)
    H[n] += setelt(ll[2])*setelt(rl[1])*op(sites,"Sz",n)
    H[n] += setelt(ll[3])*setelt(rl[2])*op(sites,"Sz",n)
    H[n] += setelt(ll[4])*setelt(rl[3])*op(sites,"Sz",n)
    H[n] += setelt(ll[4])*setelt(rl[1])*op(sites,"Sx",n)*h[n]
  end
  H[1] *= setelt(link[1][4])
  H[N] *= setelt(link[N+1][1])
  return H
end

function fourSiteIsingMPO(sites)::MPO
  H = MPO(sites)
  N = length(H)
  link = Vector{Index}(undef, N+1)
  for n=1:N+1
    link[n] = Index(5,"Link,l=$(n-1)")
  end
  for n=1:N
    s = sites[n]
    ll = link[n]
    rl = link[n+1]
    H[n] = ITensor(ll,s,s',rl)
    H[n] += setelt(ll[1])*setelt(rl[1])*op(sites,"Id",n)
    H[n] += setelt(ll[5])*setelt(rl[5])*op(sites,"Id",n)
    H[n] += setelt(ll[2])*setelt(rl[1])*op(sites,"Sz",n)
    H[n] += setelt(ll[3])*setelt(rl[2])*op(sites,"Sz",n)
    H[n] += setelt(ll[4])*setelt(rl[3])*op(sites,"Sz",n)
    H[n] += setelt(ll[5])*setelt(rl[4])*op(sites,"Sz",n)
  end
  H[1] *= setelt(link[1][5])
  H[N] *= setelt(link[N+1][1])
  return H
end

@testset "AutoMPO" begin

  N = 10

  @testset "Show MPOTerm" begin
    ampo = AutoMPO()
    add!(ampo,"Sz",1,"Sz",2)
    @test length(sprint(show,ITensors.data(ampo)[1])) > 1
  end

  @testset "Show AutoMPO" begin
    ampo = AutoMPO()
    add!(ampo,"Sz",1,"Sz",2)
    add!(ampo,"Sz",2,"Sz",3)
    @test length(sprint(show,ampo)) > 1
  end

  @testset "Single creation op" begin
    ampo = AutoMPO()
    add!(ampo,"Adagup",3)
    sites = siteinds("Electron",N)
    W = MPO(ampo,sites)
    psi = makeRandomMPS(sites)
    cdu_psi = copy(psi)
    cdu_psi[3] = noprime(cdu_psi[3]*op(sites,"Adagup",3))
    @test inner(psi,W,psi) ≈ inner(cdu_psi,psi)
  end

  @testset "Ising" begin
    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
    end
    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Ising-Different Order" begin
    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j+1,"Sz",j)
    end
    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)
    He = isingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Heisenberg" begin
    ampo = AutoMPO()
    h = rand(N) #random magnetic fields
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo,h[j],"Sz",j)
    end

    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)
    He = heisenbergMPO(sites,h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Multiple Onsite Ops" begin
    sites = siteinds("S=1",N)
    ampo1 = AutoMPO()
    for j=1:N-1
      add!(ampo1,"Sz",j,"Sz",j+1)
      add!(ampo1,0.5,"S+",j,"S-",j+1)
      add!(ampo1,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo1,"Sz*Sz",j)
    end
    Ha1 = MPO(ampo1,sites)

    ampo2 = AutoMPO()
    for j=1:N-1
      add!(ampo2,"Sz",j,"Sz",j+1)
      add!(ampo2,0.5,"S+",j,"S-",j+1)
      add!(ampo2,0.5,"S-",j,"S+",j+1)
    end
    for j=1:N
      add!(ampo2,"Sz",j,"Sz",j)
    end
    Ha2 = MPO(ampo2,sites)

    He = heisenbergMPO(sites,ones(N),"Sz*Sz")
    psi = makeRandomMPS(sites)
    Oe = inner(psi,He,psi)
    Oa1 = inner(psi,Ha1,psi)
    @test Oa1 ≈ Oe
    Oa2 = inner(psi,Ha2,psi)
    @test Oa2 ≈ Oe
  end

  @testset "Three-site ops" begin
    ampo = AutoMPO()
    # To test version of add! taking a coefficient
    add!(ampo,1.0,"Sz",1,"Sz",2,"Sz",3)
    @test length(ITensors.data(ampo)) == 1
    for j=2:N-2
      add!(ampo,"Sz",j,"Sz",j+1,"Sz",j+2)
    end
    h = ones(N)
    for j=1:N
      add!(ampo,h[j],"Sx",j)
    end
    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)
    He = threeSiteIsingMPO(sites,h)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Four-site ops" begin
    ampo = AutoMPO()
    for j=1:N-3
      add!(ampo,"Sz",j,"Sz",j+1,"Sz",j+2,"Sz",j+3)
    end
    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)
    He = fourSiteIsingMPO(sites)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
  end

  @testset "Next-neighbor Heisenberg" begin
    ampo = AutoMPO()
    J1 = 1.0
    J2 = 0.5
    for j=1:N-1
      add!(ampo,J1,  "Sz",j,"Sz",j+1)
      add!(ampo,J1*0.5,"S+",j,"S-",j+1)
      add!(ampo,J1*0.5,"S-",j,"S+",j+1)
    end
    for j=1:N-2
      add!(ampo,J2,  "Sz",j,"Sz",j+2)
      add!(ampo,J2*0.5,"S+",j,"S-",j+2)
      add!(ampo,J2*0.5,"S-",j,"S+",j+2)
    end
    sites = siteinds("S=1/2",N)
    Ha = MPO(ampo,sites)

    He = NNheisenbergMPO(sites,J1,J2)
    psi = makeRandomMPS(sites)
    Oa = inner(psi,Ha,psi)
    Oe = inner(psi,He,psi)
    @test Oa ≈ Oe
    #@test maxlinkdim(Ha) == 8
  end

  @testset "Onsite Regression Test" begin
    sites = siteinds("S=1",4)
    ampo = AutoMPO()
    add!(ampo, 0.5, "Sx",1)
    add!(ampo, 0.5, "Sy",1)
    H = MPO(ampo, sites)
    l = commonind(H[1],H[2])
    T = setelt(l[1])*H[1]
    O = op(sites[1],"Sx")+op(sites[1],"Sy")
    @test norm(T-0.5*O) < 1E-8

    sites = siteinds("S=1",2)
    ampo = AutoMPO()
    add!(ampo, 0.5im, "Sx",1)
    add!(ampo, 0.5, "Sy",1)
    H = MPO(ampo, sites)
    T = H[1]*H[2]
    O = im*op(sites[1],"Sx")*op(sites[2],"Id")+op(sites[1],"Sy")*op(sites[2],"Id")
    @test norm(T-0.5*O) < 1E-8
  end

  @testset "+ syntax" begin

    @testset "Single creation op" begin
      ampo = AutoMPO()
      ampo += ("Adagup",3)
      sites = siteinds("Electron",N)
      W = MPO(ampo,sites)
      psi = makeRandomMPS(sites)
      cdu_psi = copy(psi)
      cdu_psi[3] = noprime(cdu_psi[3]*op(sites,"Adagup",3))
      @test inner(psi,W,psi) ≈ inner(cdu_psi,psi)
    end

    @testset "Ising" begin
      ampo = AutoMPO()
      for j=1:N-1
        ampo += ("Sz",j,"Sz",j+1)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Ising-Different Order" begin
      ampo = AutoMPO()
      for j=1:N-1
        ampo += ("Sz",j+1,"Sz",j)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Heisenberg" begin
      ampo = AutoMPO()
      h = rand(N) #random magnetic fields
      for j=1:N-1
        ampo += ("Sz",j,"Sz",j+1)
        ampo += (0.5,"S+",j,"S-",j+1)
        ampo += (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo += (h[j],"Sz",j)
      end

      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = heisenbergMPO(sites,h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Multiple Onsite Ops" begin
      sites = siteinds("S=1",N)
      ampo1 = AutoMPO()
      for j=1:N-1
        ampo1 += ("Sz",j,"Sz",j+1)
        ampo1 += (0.5,"S+",j,"S-",j+1)
        ampo1 += (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo1 += ("Sz*Sz",j)
      end
      Ha1 = MPO(ampo1,sites)

      ampo2 = AutoMPO()
      for j=1:N-1
        ampo2 += ("Sz",j,"Sz",j+1)
        ampo2 += (0.5,"S+",j,"S-",j+1)
        ampo2 += (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo2 += ("Sz",j,"Sz",j)
      end
      Ha2 = MPO(ampo2,sites)

      He = heisenbergMPO(sites,ones(N),"Sz*Sz")
      psi = makeRandomMPS(sites)
      Oe = inner(psi,He,psi)
      Oa1 = inner(psi,Ha1,psi)
      @test Oa1 ≈ Oe
      Oa2 = inner(psi,Ha2,psi)
      @test Oa2 ≈ Oe
    end

    @testset "Three-site ops" begin
      ampo = AutoMPO()
      # To test version of add! taking a coefficient
      ampo += (1.0,"Sz",1,"Sz",2,"Sz",3)
      @test length(ITensors.data(ampo)) == 1
      for j=2:N-2
        ampo += ("Sz",j,"Sz",j+1,"Sz",j+2)
      end
      h = ones(N)
      for j=1:N
        ampo += (h[j],"Sx",j)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = threeSiteIsingMPO(sites,h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Four-site ops" begin
      ampo = AutoMPO()
      for j=1:N-3
        ampo += ("Sz",j,"Sz",j+1,"Sz",j+2,"Sz",j+3)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = fourSiteIsingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Next-neighbor Heisenberg" begin
      ampo = AutoMPO()
      J1 = 1.0
      J2 = 0.5
      for j=1:N-1
        ampo += (J1,  "Sz",j,"Sz",j+1)
        ampo += (J1*0.5,"S+",j,"S-",j+1)
        ampo += (J1*0.5,"S-",j,"S+",j+1)
      end
      for j=1:N-2
        ampo += (J2,  "Sz",j,"Sz",j+2)
        ampo += (J2*0.5,"S+",j,"S-",j+2)
        ampo += (J2*0.5,"S-",j,"S+",j+2)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)

      He = NNheisenbergMPO(sites,J1,J2)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
      #@test maxlinkdim(Ha) == 8
    end

    #@testset "-= syntax" begin
    #  ampo = AutoMPO()
    #  ampo += (-1,"Sz",1,"Sz",2)
    #  ampo2 = AutoMPO()
    #  ampo2 -= ("Sz",1,"Sz",2)
    #  @test ampo == ampo2
    #end

    @testset "Onsite Regression Test" begin
      sites = siteinds("S=1",4)
      ampo = AutoMPO()
      ampo += (0.5, "Sx",1)
      ampo += (0.5, "Sy",1)
      H = MPO(ampo, sites)
      l = commonind(H[1],H[2])
      T = setelt(l[1])*H[1]
      O = op(sites[1],"Sx")+op(sites[1],"Sy")
      @test norm(T-0.5*O) < 1E-8

      sites = siteinds("S=1",2)
      ampo = AutoMPO()
      ampo += (0.5im, "Sx",1)
      ampo += (0.5, "Sy",1)
      H = MPO(ampo, sites)
      T = H[1]*H[2]
      O = im*op(sites[1],"Sx")*op(sites[2],"Id")+op(sites[1],"Sy")*op(sites[2],"Id")
      @test norm(T-0.5*O) < 1E-8
    end

  end

  @testset ".+= and .-= syntax" begin

    #@testset ".-= syntax" begin
    #  ampo = AutoMPO()
    #  ampo .+= (-1,"Sz",1,"Sz",2)
    #  ampo2 = AutoMPO()
    #  ampo2 .-= ("Sz",1,"Sz",2)
    #  @test ampo == ampo2
    #end

    @testset "Single creation op" begin
      ampo = AutoMPO()
      ampo .+= ("Adagup",3)
      sites = siteinds("Electron",N)
      W = MPO(ampo,sites)
      psi = makeRandomMPS(sites)
      cdu_psi = copy(psi)
      cdu_psi[3] = noprime(cdu_psi[3]*op(sites,"Adagup",3))
      @test inner(psi,W,psi) ≈ inner(cdu_psi,psi)
    end

    @testset "Ising" begin
      ampo = AutoMPO()
      for j=1:N-1
        ampo .+= ("Sz",j,"Sz",j+1)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Ising-Different Order" begin
      ampo = AutoMPO()
      for j=1:N-1
        ampo .+= ("Sz",j+1,"Sz",j)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = isingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Heisenberg" begin
      ampo = AutoMPO()
      h = rand(N) #random magnetic fields
      for j=1:N-1
        ampo .+= ("Sz",j,"Sz",j+1)
        ampo .+= (0.5,"S+",j,"S-",j+1)
        ampo .+= (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo .+= (h[j],"Sz",j)
      end

      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = heisenbergMPO(sites,h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Multiple Onsite Ops" begin
      sites = siteinds("S=1",N)
      ampo1 = AutoMPO()
      for j=1:N-1
        ampo1 .+= ("Sz",j,"Sz",j+1)
        ampo1 .+= (0.5,"S+",j,"S-",j+1)
        ampo1 .+= (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo1 .+= ("Sz*Sz",j)
      end
      Ha1 = MPO(ampo1,sites)

      ampo2 = AutoMPO()
      for j=1:N-1
        ampo2 .+= ("Sz",j,"Sz",j+1)
        ampo2 .+= (0.5,"S+",j,"S-",j+1)
        ampo2 .+= (0.5,"S-",j,"S+",j+1)
      end
      for j=1:N
        ampo2 .+= ("Sz",j,"Sz",j)
      end
      Ha2 = MPO(ampo2,sites)

      He = heisenbergMPO(sites,ones(N),"Sz*Sz")
      psi = makeRandomMPS(sites)
      Oe = inner(psi,He,psi)
      Oa1 = inner(psi,Ha1,psi)
      @test Oa1 ≈ Oe
      Oa2 = inner(psi,Ha2,psi)
      @test Oa2 ≈ Oe
    end

    @testset "Three-site ops" begin
      ampo = AutoMPO()
      # To test version of add! taking a coefficient
      ampo .+= (1.0,"Sz",1,"Sz",2,"Sz",3)
      @test length(ITensors.data(ampo)) == 1
      for j=2:N-2
        ampo .+= ("Sz",j,"Sz",j+1,"Sz",j+2)
      end
      h = ones(N)
      for j=1:N
        ampo .+= (h[j],"Sx",j)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = threeSiteIsingMPO(sites,h)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Four-site ops" begin
      ampo = AutoMPO()
      for j=1:N-3
        ampo .+= ("Sz",j,"Sz",j+1,"Sz",j+2,"Sz",j+3)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)
      He = fourSiteIsingMPO(sites)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
    end

    @testset "Next-neighbor Heisenberg" begin
      ampo = AutoMPO()
      J1 = 1.0
      J2 = 0.5
      for j=1:N-1
        ampo .+= (J1,  "Sz",j,"Sz",j+1)
        ampo .+= (J1*0.5,"S+",j,"S-",j+1)
        ampo .+= (J1*0.5,"S-",j,"S+",j+1)
      end
      for j=1:N-2
        ampo .+= (J2,  "Sz",j,"Sz",j+2)
        ampo .+= (J2*0.5,"S+",j,"S-",j+2)
        ampo .+= (J2*0.5,"S-",j,"S+",j+2)
      end
      sites = siteinds("S=1/2",N)
      Ha = MPO(ampo,sites)

      He = NNheisenbergMPO(sites,J1,J2)
      psi = makeRandomMPS(sites)
      Oa = inner(psi,Ha,psi)
      Oe = inner(psi,He,psi)
      @test Oa ≈ Oe
      #@test maxlinkdim(Ha) == 8
    end

    @testset "Onsite Regression Test" begin
      sites = siteinds("S=1",4)
      ampo = AutoMPO()
      ampo .+= (0.5, "Sx",1)
      ampo .+= (0.5, "Sy",1)
      H = MPO(ampo, sites)
      l = commonind(H[1],H[2])
      T = setelt(l[1])*H[1]
      O = op(sites[1],"Sx")+op(sites[1],"Sy")
      @test norm(T-0.5*O) < 1E-8

      sites = siteinds("S=1",2)
      ampo = AutoMPO()
      ampo .+= (0.5im, "Sx",1)
      ampo .+= (0.5, "Sy",1)
      H = MPO(ampo, sites)
      T = H[1]*H[2]
      O = im*op(sites[1],"Sx")*op(sites[2],"Id")+op(sites[1],"Sy")*op(sites[2],"Id")
      @test norm(T-0.5*O) < 1E-8
    end

  end

  @testset "Fermionic Operators" begin
    N = 5
    s = siteinds("Fermion",N)

    a1 = AutoMPO()
    a1 += ("Cdag",1,"C",3)
    M1 = MPO(a1,s)

    a2 = AutoMPO()
    a2 += (-1,"C",3,"Cdag",1)
    M2 = MPO(a2,s)

    a3 = AutoMPO()
    a3 += ("Cdag",1,"N",2,"C",3)
    M3 = MPO(a3,s)

    p011 = productMPS(s,[1,2,2,1,1])
    p110 = productMPS(s,[2,2,1,1,1])

    @test inner(p110,M1,p011) ≈ -1.0
    @test inner(p110,M2,p011) ≈ -1.0
    @test inner(p110,M3,p011) ≈ -1.0

    p001 = productMPS(s,[1,1,2,1,1])
    p100 = productMPS(s,[2,1,1,1,1])

    @test inner(p100,M1,p001) ≈ +1.0
    @test inner(p100,M2,p001) ≈ +1.0
    @test inner(p100,M3,p001) ≈  0.0

    #
    # Repeat similar test but
    # with Electron sites
    #

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

    a1 = AutoMPO()
    a1 += ("Cdagup",1,"Cup",3)
    M1 = MPO(a1,s)

    a2 = AutoMPO()
    a2 += (-1,"Cdn",3,"Cdagdn",1)
    M2 = MPO(a2,s)

    p0uu = productMPS(s,[1,2,2,1,1])
    puu0 = productMPS(s,[2,2,1,1,1])
    p0ud = productMPS(s,[1,2,3,1,1])
    pdu0 = productMPS(s,[3,2,1,1,1])
    p00u = productMPS(s,[1,1,2,1,1])
    pu00 = productMPS(s,[2,1,1,1,1])
    p00d = productMPS(s,[1,1,3,1,1])
    pd00 = productMPS(s,[3,1,1,1,1])

    @test inner(puu0,M1,p0uu) ≈ -1.0
    @test inner(pdu0,M2,p0ud) ≈ -1.0
    @test inner(pu00,M1,p00u) ≈ +1.0
    @test inner(pd00,M2,p00d) ≈ +1.0

  end

  @testset "Complex AutoMPO Coefs" begin
    N = 4

    for use_qn in [false,true]
      sites = siteinds("S=1/2",N;conserve_qns=use_qn)
      ampo = AutoMPO()
      for i=1:N-1
        ampo += (+1im,"S+",i,"S-",i+1)
        ampo += (-1im,"S-",i,"S+",i+1)
      end
      H = MPO(ampo,sites)
      psiud = productMPS(sites,[1,2,1,2])
      psidu = productMPS(sites,[2,1,1,2])
      @test inner(psiud,H,psidu) ≈ +1im
      @test inner(psidu,H,psiud) ≈ -1im
    end
  end

end

nothing


text:
module HessianTest

import Calculus

using Test
using ForwardDiff
using ForwardDiff: Dual, Tag
using StaticArrays
using DiffTests

include(joinpath(dirname(@__FILE__), "utils.jl"))

#############################
# rosenbrock hardcoded test #
#############################

f = DiffTests.rosenbrock_1
x = [0.1, 0.2, 0.3]
v = f(x)
g = [-9.4, 15.6, 52.0]
h = [-66.0  -40.0    0.0;
     -40.0  130.0  -80.0;
       0.0  -80.0  200.0]

for c in (1, 2, 3), tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x)))
    println("  ...running hardcoded test with chunk size = $c and tag = $(repr(tag))")
    cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}(), tag)
    resultcfg = ForwardDiff.HessianConfig(f, DiffResults.HessianResult(x), x, ForwardDiff.Chunk{c}(), tag)

    @test eltype(resultcfg) == eltype(cfg)

    @test isapprox(h, ForwardDiff.hessian(f, x))
    @test isapprox(h, ForwardDiff.hessian(f, x, cfg))

    out = similar(x, 3, 3)
    ForwardDiff.hessian!(out, f, x)
    @test isapprox(out, h)

    out = similar(x, 3, 3)
    ForwardDiff.hessian!(out, f, x, cfg)
    @test isapprox(out, h)

    out = DiffResults.HessianResult(x)
    ForwardDiff.hessian!(out, f, x)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.gradient(out), g)
    @test isapprox(DiffResults.hessian(out), h)

    out = DiffResults.HessianResult(x)
    ForwardDiff.hessian!(out, f, x, resultcfg)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.gradient(out), g)
    @test isapprox(DiffResults.hessian(out), h)
end

cfgx = ForwardDiff.HessianConfig(sin, x)
@test_throws ForwardDiff.InvalidTagException ForwardDiff.hessian(f, x, cfgx)
@test ForwardDiff.hessian(f, x, cfgx, Val{false}()) == ForwardDiff.hessian(f,x)

########################
# test vs. Calculus.jl #
########################

for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    v = f(X)
    g = ForwardDiff.gradient(f, X)
    h = ForwardDiff.hessian(f, X)
    # finite difference approximation error is really bad for Hessians...
    @test isapprox(h, Calculus.hessian(f, X), atol=0.02)
    for c in CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x)))
        println("  ...testing $f with chunk size = $c and tag = $(repr(tag))")
        cfg = ForwardDiff.HessianConfig(f, X, ForwardDiff.Chunk{c}(), tag)
        resultcfg = ForwardDiff.HessianConfig(f, DiffResults.HessianResult(X), X, ForwardDiff.Chunk{c}(), tag)

        out = ForwardDiff.hessian(f, X, cfg)
        @test isapprox(out, h)

        out = similar(X, length(X), length(X))
        ForwardDiff.hessian!(out, f, X, cfg)
        @test isapprox(out, h)

        out = DiffResults.HessianResult(X)
        ForwardDiff.hessian!(out, f, X, resultcfg)
        @test isapprox(DiffResults.value(out), v)
        @test isapprox(DiffResults.gradient(out), g)
        @test isapprox(DiffResults.hessian(out), h)
    end
end

##########################################
# test specialized StaticArray codepaths #
##########################################

println("  ...testing specialized StaticArray codepaths")

x = rand(3, 3)
sx = StaticArrays.SArray{Tuple{3,3}}(x)

cfg = ForwardDiff.HessianConfig(nothing, x)
scfg = ForwardDiff.HessianConfig(nothing, sx)

actual = ForwardDiff.hessian(prod, x)
@test ForwardDiff.hessian(prod, sx) == actual
@test ForwardDiff.hessian(prod, sx, cfg) == actual
@test ForwardDiff.hessian(prod, sx, scfg) == actual
@test ForwardDiff.hessian(prod, sx, scfg) isa StaticArray
@test ForwardDiff.hessian(prod, sx, scfg, Val{false}()) == actual
@test ForwardDiff.hessian(prod, sx, scfg, Val{false}()) isa StaticArray

out = similar(x, 9, 9)
ForwardDiff.hessian!(out, prod, sx)
@test out == actual

out = similar(x, 9, 9)
ForwardDiff.hessian!(out, prod, sx, cfg)
@test out == actual

out = similar(x, 9, 9)
ForwardDiff.hessian!(out, prod, sx, scfg)
@test out == actual

result = DiffResults.HessianResult(x)
result = ForwardDiff.hessian!(result, prod, x)

result1 = DiffResults.HessianResult(x)
result2 = DiffResults.HessianResult(x)
result3 = DiffResults.HessianResult(x)
result1 = ForwardDiff.hessian!(result1, prod, sx)
result2 = ForwardDiff.hessian!(result2, prod, sx, ForwardDiff.HessianConfig(prod, result2, x, ForwardDiff.Chunk(x), nothing))
result3 = ForwardDiff.hessian!(result3, prod, sx, ForwardDiff.HessianConfig(prod, result3, x, ForwardDiff.Chunk(x), nothing))
@test DiffResults.value(result1) == DiffResults.value(result)
@test DiffResults.value(result2) == DiffResults.value(result)
@test DiffResults.value(result3) == DiffResults.value(result)
@test DiffResults.gradient(result1) == DiffResults.gradient(result)
@test DiffResults.gradient(result2) == DiffResults.gradient(result)
@test DiffResults.gradient(result3) == DiffResults.gradient(result)
@test DiffResults.hessian(result1) == DiffResults.hessian(result)
@test DiffResults.hessian(result2) == DiffResults.hessian(result)
@test DiffResults.hessian(result3) == DiffResults.hessian(result)

sresult1 = DiffResults.HessianResult(sx)
sresult2 = DiffResults.HessianResult(sx)
sresult3 = DiffResults.HessianResult(sx)
sresult1 = ForwardDiff.hessian!(sresult1, prod, sx)
sresult2 = ForwardDiff.hessian!(sresult2, prod, sx, ForwardDiff.HessianConfig(prod, sresult2, x, ForwardDiff.Chunk(x), nothing))
sresult3 = ForwardDiff.hessian!(sresult3, prod, sx, ForwardDiff.HessianConfig(prod, sresult3, x, ForwardDiff.Chunk(x), nothing))
@test DiffResults.value(sresult1) == DiffResults.value(result)
@test DiffResults.value(sresult2) == DiffResults.value(result)
@test DiffResults.value(sresult3) == DiffResults.value(result)
@test DiffResults.gradient(sresult1) == DiffResults.gradient(result)
@test DiffResults.gradient(sresult2) == DiffResults.gradient(result)
@test DiffResults.gradient(sresult3) == DiffResults.gradient(result)
@test DiffResults.hessian(sresult1) == DiffResults.hessian(result)
@test DiffResults.hessian(sresult2) == DiffResults.hessian(result)
@test DiffResults.hessian(sresult3) == DiffResults.hessian(result)

end # module