text:
"""
    TestUtils

Some commom test utilities for portfolio optimization formulations.

Mainly:
    - `get_test_data` that returns a TimeArray of Prices for 6 assets.
    - `backtest_po` that provides a basic backtest using provided strategy and returns data.
"""
module TestUtils
    using Distributions
    using JuMP
    using MarketData
    using LinearAlgebra
    import Reexport

    include("auxilary_functions.jl")
    include("backtest.jl")

    export readjust_volumes,
        backtest_po,
        get_test_data,
        mean_variance,
        returns_montecarlo

    Reexport.@reexport using MarketData
end


text:
@with_kw mutable struct L2ballParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    r::T = 1
    c::Vector{T}
end

struct L2ball <: AbstractProx
    params::L2ballParameters{Float64}

    function (::Type{L2ball})(; kwargs...)
        return new(L2ballParameters(; kwargs...))
    end
end

function prox!(l2ball::L2ball,step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    @unpack r,c = l2ball.params
    radius = norm(xprev - step*gcurr,2)
    scaling = r/max(radius,r)
    xcurr .= c + scaling*(xprev-c)
end

ext: jl


text:
#=

Author: Benjamin Hemingway

This code returns the Gaussian quadrature nodes and weights associated
with a Truncated Normal distribution.

This code is based on the method outlined by John Burkardt in the following
presentation:
https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
=#

#---------------------------------------------------------------------#
# Calculate Nodes and Weights
#---------------------------------------------------------------------#

function momentGaussQuad(M::Array{T,2}) where T<:Real

  if size(M)[1]!=size(M)[2]
    error("Hankel Matrix M must be square")
  end

  # 1) Hankel matrix already constructed
  n = (size(M)[1]-1)

  # 2) Find Upper-trinagular Cholesky decomposition: M=R'*R
  # R = chol(M)
  R = (cholesky(Hermitian(M))).U

  # 3) Compute a vector α of length n and vector β of length (n-1).
  α = Array{eltype(M),1}(undef,n)
  β = Array{eltype(M),1}(undef,n-1)

  α[1] = R[1,2] / R[1,1]

  for i = 2:n
    α[i] = R[i,i+1] / R[i,i] - R[i-1,i] / R[i-1,i-1]
  end

  for i = 1:n-1
    β[i] = R[i+1,i+1] / R[i,i]
  end

  # 4) Use α and β to generate the Golub-Welsch matrix J.

  J = Array{eltype(M),2}(undef,n,n)
  for i=1:n, j=1:n
    if i==j
      J[i,j] = α[i]
    elseif i==j+1
      J[i,j] = β[j]
    elseif i==j-1
      J[i,j] = β[i]
    else
      J[i,j] = 0.0
    end
  end

  # 5) Determine the eigenvalues L and eigenvectors V of J
  #   (in matlab use: [V, L] = eig(J) )
  dcomp = eigen(J,permute=true,scale=true)

  # 6) Generate the nodes and weights using:
  #     x = diag(L)
  #     w = μ0*V(1,1:m).^2
  x = dcomp.values
  w = vec(dcomp.vectors[1,1:n].^2)

  #= Weights are unadjusted...
  # 7) adjust weights
  dist = Distributions.Normal(μ,σ)
  Gub = Distributions.cdf(dist,ub)
  Glb = Distributions.cdf(dist,lb)

  return x, (Gub-Glb)*w
  =#
  return x, w
end


text:
using Test
using Clang.wrap_c
include("strip_line_numbers.jl")

@testset "wrap_c" begin
    wc = wrap_c.init(;  headers = [joinpath(@__DIR__, "cxx/cbasic.h")])
    hdr_tunits = wrap_c.parse_c_headers(wc)

    # test func1
    tunit = hdr_tunits[joinpath(@__DIR__, "cxx/cbasic.h")]
    func1 = search(tunit, "func1")[1]
    buf = Any[]

    wrap_c.wrap(wc, buf, func1, "test")

    exc = :(
        function func1(a::Cint,b::Cdouble,c,d)
            ccall((:func1,test),Cint,(Cint,Cdouble,Ptr{Cdouble},Ptr{Cvoid}),a,b,c,d)
        end)

    strip_line_numbers!(exc)

    @test buf[1] == exc
end


text:
export BeijingRing, taylor_hessian, local_hessian, hdata
export rings_init!, nrings, beijingring!, collect_hessian

const rings = Vector{Float64}[]
rings_init!() = empty!(rings)
nrings() = length(rings)

function collect_hessian()
    N = nrings()
    hess = [hdata((i,j)) for i=1:N, j=1:N]
end

function beijingring!(x::Float64)
    nr = length(rings)+1
    r = zeros(Float64, nr*2-1)
    push!(rings, r)
    BeijingRing(x, 0.0, nr)
end

struct BeijingRing{T}
    x::T
    g::T
    index::Int
end

NiLang.AD.grad(hd::BeijingRing) = hd.g
NiLang.value(hd::BeijingRing) = hd.x

function hdata(h::Tuple{Int,Int})
    i1, i2 = h
    i1 > i2 ? rings[i1][i2] : rings[i2][end-i1+1]
end

hdata(h::Tuple{BeijingRing{T},BeijingRing{T}}) where T = hdata((h[1].index, h[2].index))

function NiLang.chfield(h::Tuple{Int, Int}, ::typeof(hdata), val)
    i1, i2 = h
    if i1 > i2
        rings[i1][i2] = val
    else
        rings[i2][end-i1+1] = val
    end
    h
end

function NiLang.chfield(h::Tuple{BeijingRing{T}, BeijingRing{T}}, ::typeof(hdata), val) where T
    NiLang.chfield((h[1].index, h[2].index), hdata, val)
    h
end

function NiLang.chfield(hd::BeijingRing, ::typeof(value), val)
    chfield(hd, Val(:x), val)
end
function NiLang.chfield(hd::BeijingRing, ::typeof(grad), val)
    chfield(hd, Val(:g), val)
end

function Base.zero(x::BeijingRing{T}) where T
    zero(BeijingRing{T})
end

function Base.zero(x::Type{BeijingRing{T}}) where T
    beijingring!(zero(T))
end

function NiLangCore.deanc(x::BeijingRing, val::BeijingRing)
    pop!(NiLang.AD.rings)
    pop!(NiLang.AD.rings)
    value(x) == value(val)
end

@i function NEG(x!::BeijingRing)
    NEG(x!.x)
    # hessian from hessian
    for i=1:nrings()
        NEG(hdata((x!.index, i)))
        NEG(hdata((i, x!.index)))
    end

    # update gradients
    NEG(grad(x!))
end

@i function CONJ(x!::BeijingRing)
    CONJ(x!.x)
    # hessian from hessian
    for i=1:nrings()
        CONJ(hdata((x!.index, i)))
        CONJ(hdata((i, x!.index)))
    end

    # update gradients
    CONJ(grad(x!))
end

@i function ⊖(identity)(out!::BeijingRing, x::BeijingRing)
    ⊖(identity)(out!.x, x.x)
    # hessian from hessian
    for i=1:nrings()
        hdata((x.index, i)) ⊕ hdata((out!.index, i))
    end
    for i=1:nrings()
        hdata((i, x.index)) ⊕ hdata((i, out!.index))
    end

    # update gradients
    grad(x) ⊕ grad(out!)
end

@i function ⊖(abs)(out!::BeijingRing, x::BeijingRing)
    ⊖(abs)(out!.x, x.x)
    # hessian from hessian
    if (x.x > 0, ~)
        for i=1:nrings()
            hdata((x.index, i)) ⊕ hdata((out!.index, i))
        end
        for i=1:nrings()
            hdata((i, x.index)) ⊕ hdata((i, out!.index))
        end

        # update gradients
        grad(x) ⊕ grad(out!)
    else
        for i=1:nrings()
            hdata((x.index, i)) ⊖ hdata((out!.index, i))
        end
        for i=1:nrings()
            hdata((i, x.index)) ⊖ hdata((i, out!.index))
        end

        # update gradients
        grad(x) ⊖ grad(out!)
    end
end

@i function ⊖(exp)(out!::BeijingRing{T}, x::BeijingRing) where T
    @anc expx = zero(T)
    expx += exp(x.x)
    out!.x ⊕ expx
    # hessian from hessian
    for i=1:nrings()
        hdata((x.index, i)) += expx * hdata((out!.index, i))
    end
    for i=1:nrings()
        hdata((i, x.index)) += expx * hdata((i, out!.index))
    end
    hdata((x, x)) += expx * grad(out!)

    # update gradients
    grad(x) += expx * grad(out!)
    expx -= exp(x.x)
end

@i function ⊖(log)(out!::BeijingRing{T}, x::BeijingRing) where T
    @anc g = zero(T)
    @anc h = zero(T)
    out!.x += log(x.x)

    @routine r1 begin
        g += 1.0/x.x
        h += g/x.x
        NEG(h)
    end
    # hessian from hessian
    for i=1:nrings()
        hdata((x.index, i)) += g * hdata((out!.index, i))
    end
    for i=1:nrings()
        hdata((i, x.index)) += g * hdata((i, out!.index))
    end
    hdata((x, x)) += h * grad(out!)

    # update gradients
    grad(x) += g * grad(out!)

    ~@routine r1
end

@i function ⊖(sin)(out!::BeijingRing{T}, x::BeijingRing) where T
    @anc sinx = zero(T)
    @anc cosx = zero(T)
    sinx += sin(x.x)
    cosx += cos(x.x)
    out!.x ⊕ sinx
    # hessian from hessian
    for i=1:nrings()
        hdata((x.index, i)) += cosx * hdata((out!.index, i))
    end
    for i=1:nrings()
        hdata((i, x.index)) += cosx * hdata((i, out!.index))
    end
    hdata((x, x)) -= sinx * grad(out!)

    # update gradients
    grad(x) += cosx * grad(out!)
    sinx -= sin(x.x)
    cosx -= cos(x.x)
end

@i function ⊖(cos)(out!::BeijingRing{T}, x::BeijingRing) where T
    @anc sinx = zero(T)
    @anc cosx = zero(T)
    sinx += sin(x.x)
    cosx += cos(x.x)
    out!.x ⊕ cosx
    # hessian from hessian
    for i=1:nrings()
        hdata((x.index, i)) -= sinx * hdata((out!.index, i))
    end
    for i=1:nrings()
        hdata((i, x.index)) -= sinx * hdata((i, out!.index))
    end
    hdata((x, x)) -= cosx * grad(out!)

    # update gradients
    grad(x) -= sinx * grad(out!)
    sinx -= sin(x.x)
    cosx -= cos(x.x)
end

@i function SWAP(x!::BeijingRing, y!::BeijingRing)
    SWAP(x!.x, y!.x)
    # hessian from hessian
    for i=1:nrings()
        SWAP(hdata((x!.index, i)), hdata((y!.index, i)))
    end
    for i=1:nrings()
        SWAP(hdata((i, x!.index)), hdata((i, y!.index)))
    end

    # update gradients
    SWAP(grad(x!), grad(y!))
end

# dL^2/dx/dy = ∑(dL^2/da/db)*da/dx*db/dy
# https://arxiv.org/abs/1206.6464
@i function ⊖(*)(out!::BeijingRing, x, y)
    ⊖(*)(out!.x, value(x), value(y))
    # hessian from hessian
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((x.index, i)) += value(y) * hdata((out!.index, i))
        end
        if (y isa BeijingRing, ~)
            hdata((y.index, i)) += value(x) * hdata((out!.index, i))
        end
    end
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((i, x.index)) += value(y) * hdata((i, out!.index))
        end
        if (y isa BeijingRing, ~)
            hdata((i, y.index)) += value(x) * hdata((i, out!.index))
        end
    end

    # hessian from jacobian
    if (x isa BeijingRing && y isa BeijingRing, ~)
        hdata((x, y)) ⊕ grad(out!)
        hdata((y, x)) ⊕ grad(out!)
    end

    # update gradients
    if (x isa BeijingRing, ~)
        grad(x) += grad(out!) * value(y)
    end
    if (y isa BeijingRing, ~)
        grad(y) += value(x) * grad(out!)
    end
end

@i function ⊖(/)(out!::BeijingRing{T}, x, y) where T
    ⊖(/)(out!.x, value(x), value(y))
    @anc binv = zero(T)
    @anc binv2 = zero(T)
    @anc binv3 = zero(T)
    @anc a3 = zero(T)
    @anc xjac = zero(T)
    @anc yjac = zero(T)
    @anc yyjac = zero(T)
    @anc xyjac = zero(T)

    @routine jacs begin
        # compute dout/dx and dout/dy
        xjac += 1.0/value(y)
        binv2 += xjac^2
        binv3 += xjac^3
        yjac -= value(x)*binv2
        a3 += value(x)*binv3
        yyjac += 2*a3
        xyjac ⊖ binv2
    end
    # hessian from hessian
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((x.index, i)) += xjac * hdata((out!.index, i))
        end
        if (y isa BeijingRing, ~)
            hdata((y.index, i)) += yjac * hdata((out!.index, i))
        end
    end
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((i, x.index)) += xjac * hdata((i, out!.index))
        end
        if (y isa BeijingRing, ~)
            hdata((i, y.index)) += yjac * hdata((i, out!.index))
        end
    end

    # hessian from jacobian
    if (x isa BeijingRing && y isa BeijingRing, ~)
        hdata((y.index, y.index)) += yyjac*grad(out!)
        hdata((x.index, y.index)) += xyjac*grad(out!)
        hdata((y.index, x.index)) += xyjac*grad(out!)
    end

    # update gradients
    if (x isa BeijingRing, ~)
        grad(x) += grad(out!) * xjac
    end
    if (y isa BeijingRing, ~)
        grad(y) += yjac * grad(out!)
    end

    ~@routine jacs
end

@i function ⊖(^)(out!::BeijingRing{T}, x, n) where T
    ⊖(^)(out!.x, value(x), value(n))
    @anc logx = zero(T)
    @anc logx2 = zero(T)
    @anc powerxn = zero(T)
    @anc anc1 = zero(T)
    @anc anc2 = zero(T)
    @anc xjac = zero(T)
    @anc njac = zero(T)
    @anc hxn = zero(T)
    @anc hxx = zero(T)
    @anc hnn = zero(T)
    @anc nminus1 = zero(T)

    # compute jacobians
    @routine getjac begin
        nminus1 ⊕ value(n)
        nminus1 ⊖ 1
        powerxn += value(x)^value(n)
        logx += log(value(x))
        out!.x ⊕ powerxn

        # dout!/dx = n*x^(n-1)
        anc1 += value(x)^nminus1
        xjac += anc1 * value(n)
        # dout!/dn = logx*x^n
        njac += logx*powerxn

        # for hessian
        logx2 += logx^2
        anc2 += xjac/value(x)
        hxn ⊕ anc1
        hxn += xjac * logx
        hxx += anc2 * nminus1
        hnn += logx2 * powerxn
    end

    # hessian from hessian
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((i, x.index)) += hdata((i, out!.index)) * xjac
        end
        if (n isa BeijingRing, ~)
            hdata((i, n.index)) += hdata((i, out!.index)) * njac
        end
    end
    for i=1:nrings()
        if (x isa BeijingRing, ~)
            hdata((x.index, i)) += hdata((out!.index, i)) * xjac
        end
        if (n isa BeijingRing, ~)
            hdata((n.index, i)) += hdata((out!.index, i)) * njac
        end
    end

    # hessian from jacobian
    # Dnn = x^n*log(x)^2
    # Dxx = (-1 + n)*n*x^(-2 + n)
    # Dxn = Dnx = x^(-1 + n) + n*x^(-1 + n)*log(x)
    if (x isa BeijingRing && n isa BeijingRing, ~)
        hdata((x, x)) += hxx * grad(out!)
        hdata((n, n)) += hnn * grad(out!)
        hdata((x, n)) += hxn * grad(out!)
        hdata((n, x)) += hxn * grad(out!)
    end

    # update gradients
    if (x isa BeijingRing, ~)
        grad(x) += grad(out!) * xjac
    end
    if (n isa BeijingRing, ~)
        grad(n) += grad(out!) * njac
    end

    ~@routine getjac
end

@i function IROT(a!::BeijingRing{T}, b!, θ) where T
    @anc s = zero(T)
    @anc c = zero(T)
    @anc ca = zero(T)
    @anc sb = zero(T)
    @anc sa = zero(T)
    @anc cb = zero(T)
    @anc θ2 = zero(T)
    IROT(value(a!), value(b!), value(θ))

    @routine temp begin
        θ2 ⊖ value(θ)
        θ2 ⊖ π/2
        s += sin(value(θ))
        c += cos(value(θ))
        ca += c * value(a!)
        sb += s * value(b!)
        sa += s * value(a!)
        cb += c * value(b!)
    end

    # update gradient, #1
    for i=1:nrings()
        ROT(hdata((i, a!.index)), hdata((i, b!.index)), θ2)
        hdata((i, θ.index)) += value(a!) * hdata((i, a!.index))
        hdata((i, θ.index)) += value(b!) * hdata((i, b!.index))
        ROT(hdata((i, a!.index)), hdata((i, b!.index)), π/2)
    end
    for i=1:nrings()
        ROT(hdata((a!.index, i)), hdata((b!.index, i)), θ2)
        hdata((θ.index, i)) += value(a!) * hdata((a!.index, i))
        hdata((θ.index, i)) += value(b!) * hdata((b!.index, i))
        ROT(hdata((a!.index, i)), hdata((b!.index, i)), π/2)
    end

    # update local hessian
    hdata((a!, θ)) -= s * grad(a!)
    hdata((b!, θ)) -= c * grad(a!)
    hdata((θ, a!)) -= s * grad(a!)
    hdata((θ, b!)) -= c * grad(a!)
    hdata((θ, θ)) -= ca * grad(a!)
    hdata((θ, θ)) += sb * grad(a!)

    hdata((a!, θ)) += c * grad(b!)
    hdata((b!, θ)) -= s * grad(b!)
    hdata((θ, a!)) += c * grad(b!)
    hdata((θ, b!)) -= s * grad(b!)
    hdata((θ, θ)) -= sa * grad(b!)
    hdata((θ, θ)) -= cb * grad(b!)

    # update gradients
    ROT(grad(a!), grad(b!), θ2)
    grad(θ) += value(a!) * grad(a!)
    grad(θ) += value(b!) * grad(b!)
    ROT(grad(a!), grad(b!), π/2)

    ~@routine temp
end

function local_hessian(f, args; kwargs=())
    nargs = length(args)
    hes = zeros(nargs,nargs,nargs)
    @instr f(args...)
    for j=1:nargs
        rings_init!()
        largs = [beijingring!(arg) for arg in args]
        @instr grad(largs[j]) ⊕ 1.0
        @instr (~f)(largs...)
        hes[:,:,j] .= collect_hessian()
    end
    hes
end

function (h::Hessian)(args...; kwargs...)
    @assert count(x -> x isa Loss, args) == 1
    N = length(args)

    iloss = 0
    for i=1:length(args)
        if tget(args,i) isa Loss
            iloss += identity(i)
        end
    end
    @instr (~Loss)(tget(args, iloss))

    @instr h.f(args...)
    rings_init!()
    args = [beijingring!(x) for x in args]
    @instr grad(args[iloss]) ⊕ 1.0
    @instr (~h.f)(args...)
    @instr Loss(tget(args, iloss))
    args
end

macro nohess(ex)
    @match ex begin
        :($f($(args...))) => begin
            newargs = []
            for arg in args
                push!(newargs, @match arg begin
                    :($x::BeijingRing) => :($x.x)
                    :($x::BeijingRing{$tp}) => :($x.x)
                    _ => arg
                end
                )
            end
            esc(quote
                @i function $f($(args...))
                    $f($(newargs...))
                end
            end)
        end
        _ => error("expect `f(args...)`, got $ex")
    end
end

@nohess ⊖(identity)(out!::BeijingRing, x)
@nohess ⊖(identity)(out!, x::BeijingRing)
for op in [:*, :/, :^]
    @eval @nohess ⊖($op)(out!::BeijingRing, x::Number, y::Number)
    @eval @nohess ⊖($op)(out!::Number, x, y::BeijingRing)
    @eval @nohess ⊖($op)(out!::Number, x::BeijingRing, y::BeijingRing)
    @eval @nohess ⊖($op)(out!::Number, x::BeijingRing, y)
end


text:
function test_orthogonality(chars)
    projs = [SymbolicWedderburn.matrix_projection(χ) for χ in chars]

    res = Matrix{Bool}(undef, length(projs), length(projs))

    for j in axes(res, 2)
        for i in axes(res, 1)
            p, q = projs[i], projs[j]
            if p == q
                # p^2 == p
                res[i,j] = all(x -> isapprox(0.0, x; atol = 1e-12), p*p - p)
            else
                res[i,j] = all(x -> isapprox(0.0, x; atol = 1e-12), p*q)
            end
        end
    end
    return all(res)
end

@testset "Orthogonality of projections" begin
    G = PermGroup([perm"(1,2,3,4)"])
    chars = let irr = SymbolicWedderburn.irreducible_characters(G)
        if !all(isreal, irr)
            irr, _ = SymbolicWedderburn.affordable_real(irr)
        end
        SymbolicWedderburn.Character{Float64}.(irr)
    end
    @test test_orthogonality(chars)
    @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) ==
          PermutationGroups.degree(G)

    G = PermGroup([perm"(1,2,3,4)", perm"(1,2)"])
    chars = SymbolicWedderburn.irreducible_characters(G)
    @test test_orthogonality(chars)
    @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) ==
          PermutationGroups.degree(G)

    G = PermGroup([perm"(1,2,3)", perm"(2,3,4)"])
    chars = SymbolicWedderburn.irreducible_characters(G)
    @test test_orthogonality(chars)
    @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) ==
          PermutationGroups.degree(G)

    @time for ord = 2:16
        # @testset "SmallGroup($ord, $n)"
        for (n, G) in enumerate(SmallPermGroups[ord])
            chars = SymbolicWedderburn.irreducible_characters(Rational{Int}, G)

            @test test_orthogonality(chars)
            @test sum(
                first ∘ size,
                SymbolicWedderburn.image_basis.(chars),
            ) == PermutationGroups.degree(G)
        end
    end
end


text:
using InteractiveUtils

function fn(n)
    a = Array{Float32}(undef, 5)
    a .= n
    b = Array{UInt32}([n, n, n, n, n])
    c = Array{UInt32}([n, n, n, n, n])
    d = Array{UInt32}([n, n, n, n, n])
    return a, b, c, d
end

code_warntype(fn, (UInt32,))

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
using Revise
using ZonalFlow
using JLD2

tspan   = (0.0,50.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            dense=false,
            save_noise=false,
            saveat=5
           );

domain  = Domain(extent=(2π,2π),res=(8,8));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

dn      = "data/8x8/"

eqs = [GQL(0),CE2(eigmax=true)];
labels = ["ql_test","ce2_kf_eigmax_qlic_cb"];

ql      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;u0=ql.u[end],tsargs...)
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])


text:
# Helper functions
function d_subproblem{T}(s::Array{T, 2}, mu_s::T, lambda_s::Array{T, 2},
                         ZtZ::Array{T, 2}, Ztx::Array{T, 2})
  filtLen = size(s, 2);
  dictSize = size(ZtZ, 1);
  # d = ( ZtZ + mu_s.*I ) \ ( Ztx + mu_s.*s - lambda_s )
  A = copy(ZtZ);
  for j=1:dictSize; A[j, j] += mu_s; end
  b = copy(s);
  scale!(b, mu_s);
  b += Ztx;
  b -= lambda_s;
  d = A\b;
  return d::Array{T, 2};
end

function s_subproblem{T}(d::Array{T, 2}, mu_s::T, lambda_s::Array{T, 2})
  # s = d + 1/mu_s.*lambda_s;
  s = copy(lambda_s);
  scale!(s, 1/mu_s);
  s += d;
  normalize!(s);
  return s::Array{T, 2};
end

function fsc_learning_multi{T}(data::LearningData{T}, optParams::LearningOpts{T})
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const ZtZ       = data.ZtZ;
  const Ztx       = data.Ztx;
  const d         = optParams.Dinit;
        s         = optParams.s;
        lambda_s  = optParams.lambda_s;
        mu_s::T   = optParams.mu_s;
  const max_iter  = optParams.max_iter;
  const tol       = optParams.tol;
  const DEBUG     = optParams.DEBUG;
        debugData = optParams.debugData;

  filtLen = size(Ztx, 2);
  dictSize = size(ZtZ, 1);

  assert(size(Ztx, 1)==dictSize, "Ztx does not have the correct dimensions.")
  assert(isempty(d) || size(d, 1)==dictSize,
    "Number of filters in Dinit does not match the number of channels in Z.")

  # Initialization
  if isempty(lambda_s);
    lambda_s = zeros(T, dictSize, filtLen);
  end
  if isempty(d);
    d = zeros(T, dictSize, filtLen);
  end
  if isempty(s) || vecnorm(s)==0;
    s = s_subproblem(d, mu_s, lambda_s);
  end

  if DEBUG;
    debugData["d"] = d;
    debugData["s"] = s;
    debugData["lambda_s"] = lambda_s;
    (prev_obj, rerr, spar) = compute_lagrangian(debugData);
  end

  iter = 0;
  converged = false;
  while ~converged;
    iter += 1;

    if max_iter>1; @printf("\b\b\b\b\b%05d", iter); end

    # Solve for d* in the frequency domain
    d = d_subproblem(s, mu_s, lambda_s, ZtZ, Ztx);
    if DEBUG;
      debugData["d"] = d;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("            D: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Solve for s* in the spatial domain
    s = s_subproblem(d, mu_s, lambda_s);
    if DEBUG;
      debugData["s"] = s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("            S: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    if max_iter==1; return (d, s); end # Don't do unnecessary work

    # Update the Lagrangian multipliers
    for i=1:filtLen*dictSize; lambda_s[i] += mu_s*( d[i]-s[i] ); end
    if DEBUG;
      debugData["lambda_s"] = lambda_s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("            L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Update penalty term
    mu_s = min(tau*mu_s, mu_max);
    if DEBUG;
      debugData["mu_s"] = mu_s;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("            M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    if iter>=max_iter || mean((d-s).^2)<=tol;
      converged = true;
    end
  end

  return (d, s);

end

max_stars_repo_name: bkong/FSC

max_issues_repo_name: bkong/FSC

max_forks_repo_name: bkong/FSC


text:
abstract type ExportFormat <: WeaveFormat end

function Base.getproperty(sf::T, s::Symbol) where {T<:ExportFormat}
    hasfield(T, s) && return getfield(sf, s)
    return getproperty(sf.primaryformat, s)
end
function Base.setproperty!(sf::T, s::Symbol, v) where {T<:ExportFormat}
    if hasfield(T, s)
        setfield!(sf, s, v)
    else
        setproperty!(sf.primaryformat, s, v)
    end
end
function Base.hasproperty(sf::T, s::Symbol) where {T<:ExportFormat}
    hasfield(T, s) || hasproperty(sf.primaryformat, s)
end

render_doc(df::ExportFormat, body, doc) = render_doc(df.primaryformat, body, doc)

render_chunk(df::ExportFormat, chunk) = render_chunk(df.primaryformat, chunk)
# Need to define these to avoid ambiguities
render_chunk(df::ExportFormat, chunk::DocChunk) = render_chunk(df.primaryformat, chunk)
render_chunk(df::ExportFormat, chunk::CodeChunk) = render_chunk(df.primaryformat, chunk)
render_output(df::ExportFormat, output) = render_output(df.primaryformat, output)


text:
using Zygote
using Zygote: gradient, @adjoint
using Zygote, LuxurySparse, Yao, YaoBase, SparseArrays, BitBasis
import Zygote: Context

using Yao
using YaoBlocks: ConstGate
import Yao: apply!, ArrayReg, statevec, RotationGate

using LuxurySparse, SparseArrays, LinearAlgebra
using BitBasis: controller, controldo
using TupleTools

Zygote.@adjoint (::Type{T})(perms, vals) where T <: PermMatrix = T(perms, vals), Δ -> nothing
Zygote.@adjoint (::Type{T})() where T <: IMatrix = T(), Δ -> nothing
Zygote.@adjoint Base.:(*)(A::Number, B::PermMatrix) = A * B, Δ->(sum(Δ .* B), A * Δ)
Zygote.@adjoint Base.:(*)(A::PermMatrix, B::Number) = A * B, Δ->(A * Δ, sum(Δ .* B))

Zygote.@adjoint SparseArrays.SparseMatrixCSC(A::PermMatrix) = SparseMatrixCSC(A), Δ->(Δ, )

Zygote.@adjoint BitBasis.onehot(::Type{T}, nbits::Int, x::Integer, nbatch::Int) where T = onehot(T, nbits, x, nbatch), Δ -> nothing

# upstreams
Zygote.@adjoint Base.:(-)(a, b) = a-b, Δ -> (Δ, -Δ)
Zygote.@adjoint Base.:(+)(a, b) = a+b, Δ -> (Δ, Δ)

# require mutate
Zygote.@adjoint! function copyto!(xs::AbstractVector, ys::Tuple)
    xs_ = copy(xs)
    copyto!(xs, ys), function (dxs)
        copyto!(xs_, xs)
        return (nothing, Tuple(dxs))
    end
end

@adjoint function Base.Iterators.Zip(tp)
    Base.Iterators.Zip(tp), adjy-> ((@show adjy;zip(adjy...)),)
end

@adjoint function reduce(func, xs; kwargs...)
    backs = Any[]
    ys = Any[]
    function nfunc(x, x2)
        y, back = forward(func, x, x2)
        push!(backs, back)
        push!(ys, y)
        return y
    end
    reduce(nfunc, xs; kwargs...),
    function (adjy)
        res = Vector{Any}(undef, length(ys))
        for i=length(ys):-1:1
            back, y = backs[i], ys[i]
            adjy, res[i] = back(adjy)
        end
        if !haskey(kwargs, :init)
            insert!(res, 1, adjy)
        end
        return (nothing, res)
    end
end

@adjoint function mapreduce(op, func, xs; kwargs...)
    opbacks = Any[]
    backs = Any[]
    ys = Any[]
    function nop(x)
        y, back = forward(op,x)
        push!(opbacks, back)
        y
    end
    function nfunc(x, x2)
        y, back = forward(func, x, x2)
        push!(backs, back)
        push!(ys, y)
        return y
    end
    mapreduce(nop, nfunc, xs; kwargs...),
    function (adjy)
        offset = haskey(kwargs, :init) ? 0 : 1
        res = Vector{Any}(undef, length(ys)+offset)
        for i=length(ys):-1:1
            opback, back, y = opbacks[i+offset], backs[i], ys[i]
            adjy, adjthis = back(adjy)
            res[i+offset], = opback(adjthis)
        end
        if offset==1
            res[1], = opbacks[1](adjy)
        end
        return (nothing, nothing, res)
    end
end

@adjoint function collect(::Type{T}, source::TS) where {T, TS}
    collect(T, source),
    adjy -> (adjy,)   # adjy -> convert(TS, adjy)
end

@adjoint Iterators.reverse(x::T) where T = Iterators.reverse(x), adjy->(collect(Iterators.reverse(adjy)),)  # convert is better

# data projection
@adjoint function *(sp::Union{SDSparseMatrixCSC, SDDiagonal, SDPermMatrix}, v::AbstractVector)
    sp*v, adjy -> (outer_projection(sp, adjy, v'), sp'*adjy)
end

@adjoint YaoBlocks.decode_sign(args...) = YaoBlocks.decode_sign(args...), adjy->nothing

@adjoint function *(v::LinearAlgebra.Adjoint{T, V}, sp::Union{SDSparseMatrixCSC, SDDiagonal, SDPermMatrix}) where {T, V<:AbstractVector}
    v*sp, adjy -> (adjy*sp', outer_projection(sp, v', adjy))
end

@adjoint function *(v::LinearAlgebra.Adjoint{T, V}, sp::SDDiagonal, v2::AbstractVector) where {T, V<:AbstractVector}
    v*sp, adjy -> (adjy*(sp*v2)', adjy*projection(sp, v', v2'), adjy*(v*sp)')
end

function outer_projection(y::SDSparseMatrixCSC, adjy, v)
    # adjy*v^T
    out = zero(y)
    is, js, vs = findnz(y)
    for (k,(i,j)) in enumerate(zip(is, js))
        @inbounds out.nzval[k] = adjy[i]*v[j]
    end
    out
end

outer_projection(y::SDDiagonal, adjy, v) = Diagonal(adjy.*v)

"""
Project a dense matrix to a sparse matrix
"""
function projection(y::AbstractSparseMatrix, m::AbstractMatrix)
    out = zero(y)
    is, js, vs = findnz(y)
    for (k,(i,j)) in enumerate(zip(is, js))
        @inbounds out.nzval[k] = m[i,j]
    end
    out
end

Base.zero(pm::PermMatrix) = PermMatrix(pm.perm, zero(pm.vals))
projection(y::SDDiagonal, m::AbstractMatrix) = Diagonal(diag(m))
function projection(y::PermMatrix, m::AbstractMatrix)
    res = zero(y)
    for i=1:size(res, 1)
        @inbounds res.vals[i] = m[i,res.perm[i]]
    end
    res
end
projection(x::RT, adjx::Complex) where RT<:Real = RT(real(adjx))

function rotgrad(::Type{T}, rb::RotationGate{N}) where {N, T}
    -sin(rb.theta / 2)/2 * IMatrix{1<<N}() + im/2 * cos(rb.theta / 2) * mat(T, rb.block)
end

@adjoint function mat(::Type{T}, rb::RotationGate{N, RT}) where {T, N, RT}
    mat(T, rb), adjy -> (nothing, projection(rb.theta, sum(adjy .* rotgrad(T, rb))),)
end

@adjoint function mat(::Type{T}, rb::Union{PutBlock{N, C, RT}, RT}) where {T, N, C, RT<:ConstGate.ConstantGate}
    mat(T, rb), adjy -> (nothing, nothing)
end

@adjoint function RotationGate(G, θ)
    RotationGate(G, θ), adjy->(nothing, adjy)
end

@adjoint function PutBlock{N}(block::GT, locs::NTuple{C, Int}) where {N, M, C, GT <: AbstractBlock{M}}
    PutBlock{N}(block, locs), adjy->(adjy.content, adjy.locs)
end

@adjoint function ChainBlock(blocks) where N
    ChainBlock(blocks),
    adjy -> ((@show adjy.blocks; adjy.blocks),)
end

@adjoint function chain(blocks) where N
    chain(blocks),
    adjy -> ((@show adjy.blocks; adjy.blocks),)
end

@adjoint function KronBlock{N, MT}(slots::Vector{Int}, locs::Vector{Int}, blocks::Vector{MT}) where {N, MT<:AbstractBlock}
    KronBlock{N, MT}(slots, locs, blocks),
    adjy -> (adjy.slots, adjy.locs, adjy.blocks)
end

@adjoint function ControlBlock{N, BT, C, M}(ctrl_locs, ctrl_config, block, locs) where {N, C, M, BT<:AbstractBlock}
    ControlBlock{N, BT, C, M}(ctrl_locs, ctrl_config, block, locs),
    adjy->(adjy.ctrl_locs, adjy.ctrl_config, adjy.content, adjy.locs)
end

@adjoint function YaoBlocks.cunmat(nbit::Int, cbits::NTuple{C, Int}, cvals::NTuple{C, Int}, U0::AbstractMatrix, locs::NTuple{M, Int}) where {C, M}
    y = YaoBlocks.cunmat(nbit, cbits, cvals, U0, locs)
    y, adjy-> (nothing, nothing, nothing, adjcunmat(y, adjy, nbit, cbits, cvals, U0, locs), nothing)
end

@inline function adjsetcol!(csc::SparseMatrixCSC, icol::Int, rowval::AbstractVector, nzval::SubArray)
    @inbounds begin
        S = csc.colptr[icol]
        E = csc.colptr[icol+1]-1
        nzval .+= view(csc.nzval, S:E)
    end
    csc
end

@inline function adjunij!(mat::SparseMatrixCSC, locs, U::Matrix)
    for j = 1:size(U, 2)
        @inbounds adjsetcol!(mat, locs[j], locs, view(U,:,j))
    end
    return U
end

@inline function adjunij!(mat::SparseMatrixCSC, locs, U::SparseMatrixCSC)
    for j = 1:size(U, 2)
        S = U.colptr[j]
        E = U.colptr[j+1]-1
        @inbounds adjsetcol!(mat, locs[j], locs, view(U.nzval,S:E))
    end
    return U
end

@inline function adjunij!(mat::SDDiagonal, locs, U::Diagonal)
    @inbounds U.diag .+= mat.diag[locs]
    return U
end

@inline function adjunij!(mat::SDPermMatrix, locs, U::PermMatrix)
    @inbounds U.vals .+= mat.vals[locs]
    return U
end

function adjcunmat(y::AbstractMatrix, adjy::AbstractMatrix, nbit::Int, cbits::NTuple{C, Int}, cvals::NTuple{C, Int}, U0::AbstractMatrix{T}, locs::NTuple{M, Int}) where {C, M, T}
    U, ic, locs_raw = YaoBlocks.reorder_unitary(nbit, cbits, cvals, U0, locs)
    adjy = _render_adjy(adjy, y)
    adjU = _render_adjU(U)

    ctest = controller(cbits, cvals)

    controldo(ic) do i
        adjunij!(adjy, locs_raw+i, adjU)
    end

    adjU = all(TupleTools.diff(locs).>0) ? adjU : YaoBase.reorder(adjU, collect(locs)|>sortperm|>sortperm)
    adjU
end
_render_adjy(adjy, y) = projection(y, adjy)

_render_adjU(U0::AbstractMatrix{T}) where T = zeros(T, size(U0)...)
_render_adjU(U0::SDSparseMatrixCSC{T}) where T = SparseMatrixCSC(size(U0)..., dynamicize(U0.colptr), dynamicize(U0.rowval), zeros(T, U0.nzval|>length))
_render_adjU(U0::SDDiagonal{T}) where T = Diagonal(zeros(T, size(U0, 1)))
_render_adjU(U0::SDPermMatrix{T}) where T = PermMatrix(U0.perm, zeros(T, length(U0.vals)))

function collect_gradients(st, out=Any[])
    for blk in st
        collect_gradients(blk, out)
    end
    out
end

collect_gradients(st::Number, out=any[]) = push!(out, st)
collect_gradients(st::Nothing, out=any[]) = out

############# Build a circuit (Copied from `QuAlgorithmZoo`) ###############
"""
    pair_ring(n::Int) -> Vector

Pair ring.
"""
pair_ring(n::Int) = [i=>mod(i, n)+1 for i=1:n]

"""
    cnot_entangler([n::Int, ] pairs::Vector{Pair}) = ChainBlock

Arbitrary rotation unit, support lazy construction.
"""
cnot_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>X) for (ctrl, target) in pairs)
cnot_entangler(pairs) = n->cnot_entangler(n, pairs)

"""
    rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{nbit, ComplexF64}

Arbitrary rotation unit (put in `nbit` space), set parameters notrailing, noleading true to remove trailing and leading Z gates.
"""
function rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false)
    rt = chain(nbit, [put(nbit, ibit=>Rz(0.0)), put(nbit, ibit=>Rx(0.0)), put(nbit, ibit=>Rz(0.0))])
    noleading && popfirst!(rt)
    notrailing && pop!(rt)
    rt
end

rotorset(nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [rotor(nbit, j, noleading, notrailing) for j=1:nbit])

"""
    random_diff_circuit(nbit, nlayer, pairs; do_cache=false)

A kind of widely used differentiable quantum circuit, angles in the circuit is randomely initialized.

ref:
    1. Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Chow, J. M., & Gambetta, J. M. (2017).
       Hardware-efficient Quantum Optimizer for Small Molecules and Quantum Magnets. Nature Publishing Group, 549(7671), 242–246.
       https://doi.org/10.1038/nature23879.
"""
function random_diff_circuit(nbit, nlayer, pairs; do_cache=false)
    circuit = chain(nbit)

    ent = cnot_entangler(pairs)
    if do_cache
        ent = ent |> cache
    end
    for i = 1:(nlayer + 1)
        i!=1 && push!(circuit, ent)
        push!(circuit, rotorset(nbit, i==1, i==nlayer+1))
    end
    circuit
end

"""
    heisenberg(nbit::Int; periodic::Bool=true)

heisenberg hamiltonian, for its ground state, refer `PRB 48, 6141`.
"""
function heisenberg(nbit::Int; periodic::Bool=true)
    sx = i->put(nbit, i=>X)
    sy = i->put(nbit, i=>Y)
    sz = i->put(nbit, i=>Z)
    mapreduce(i->(j=i%nbit+1; sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j)), +, 1:(periodic ? nbit : nbit-1))
end

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
# -------------------------------------------------------------------------------
# Prepare function
# -------------------------------------------------------------------------------
function prepare(::PatMat, X::AbstractMatrix, y::AbstractVector{Bool})

    ind_pos  = findall(y)
    ind_neg  = findall(.~y)
    nα       = length(ind_pos)
    nβ       = length(y)
    inv_perm = nα .+ (1:nβ)

    XX  = vcat(X[ind_pos, :], X)

    return XX, nα, nβ, ind_pos, ind_neg, inv_perm
end

function prepare(
    ::Union{AbstractTopPushK, PatMatNP},
    X::AbstractMatrix,
    y::AbstractVector{Bool}
)

    ind_pos  = findall(y)
    ind_neg  = findall(.~y)
    nα       = length(ind_pos)
    nβ       = length(ind_neg)
    inv_perm = invperm(vcat(ind_pos, ind_neg))

    XX  = vcat(X[ind_pos, :], X[ind_neg, :])

    return XX, nα, nβ, ind_pos, ind_neg, inv_perm
end

# -------------------------------------------------------------------------------
# DTrain kernel matrix
# -------------------------------------------------------------------------------
function kernelmatrix(M::AbstractModel,
                      Xtrain::AbstractArray,
                      ytrain::AbstractVector{Bool};
                      kernel::Kernel = LinearKernel(1, 0),
                      T::DataType = Float32,)

    X, nα, nβ, ind_pos, ind_neg, inv_perm = prepare(M, Xtrain, ytrain)
    n = length(ytrain)

    K = T.(MLKernels.kernelmatrix(Val(:row), kernel, X, true))
    K[1:nα, nα+1:end]   .*= -1
    K[(nα+1):end, 1:nα] .*= -1
    return K, n, nα, nβ, ind_pos, ind_neg, inv_perm
end

function save_kernelmatrix(M::AbstractModel,
                           file::AbstractString,
                           Xtrain::AbstractMatrix,
                           ytrain::AbstractVector{Bool};
                           kernel::Kernel = LinearKernel(1, 0),
                           T::DataType    = Float32,
                           partial_fill::Bool = false)

    X, nα, nβ, ind_pos, ind_neg, inv_perm = prepare(M, Xtrain, ytrain)

    n, npos, nneg = length(inv_perm), length(ind_pos), length(ind_neg)
    N             = nα + nβ

    # auxiliary
    io = open(file, "w+");
    write(io, 0)
    write(io, [nα, nβ, n, npos, nneg, N, N])
    write(io, ind_pos)
    write(io, ind_neg)
    write(io, inv_perm)

    # kernel matrix
    K = Mmap.mmap(io, Matrix{T}, (N, N))
    if partial_fill
        fill_kernelmatrix!(K, kernel, X, X)
    else
        K .= MLKernels.kernelmatrix(Val(:row), kernel, X, true)
    end
    K[1:nα, nα+1:end]   .*= -1
    K[(nα+1):end, 1:nα] .*= -1
    Mmap.sync!(K)
    close(io)

    return
end

# -------------------------------------------------------------------------------
# DValidation kernel matrix
# -------------------------------------------------------------------------------
function kernelmatrix(M::AbstractModel,
                      Xtrain::AbstractArray,
                      ytrain::AbstractVector{Bool},
                      Xvalid::AbstractArray,
                      yvalid::AbstractVector{Bool};
                      kernel::Kernel = LinearKernel(1, 0),
                      T::DataType = Float32,)

    X, nα, nβ, = prepare(M, Xtrain, ytrain)
    n          = length(yvalid)
    ind_pos    = findall(yvalid)
    ind_neg    = findall(.~yvalid)
    inv_perm   = 1:length(yvalid)

    K = T.(MLKernels.kernelmatrix(Val(:row), kernel, X, Xvalid))
    K[(nα+1):end, :] .*= -1
    return K, n, nα, nβ, ind_pos, ind_neg, inv_perm
end

function save_kernelmatrix(M::AbstractModel,
                           file::AbstractString,
                           Xtrain::AbstractMatrix,
                           ytrain::AbstractVector{Bool},
                           Xvalid::AbstractMatrix,
                           yvalid::AbstractVector{Bool};
                           kernel::Kernel = LinearKernel(1, 0),
                           T::DataType    = Float32,
                           partial_fill::Bool = false)

    X, nα, nβ, = prepare(M, Xtrain, ytrain)
    ind_pos    = findall(yvalid)
    ind_neg    = findall(.~yvalid)
    inv_perm   = 1:length(yvalid)

    n, npos, nneg = length(inv_perm), length(ind_pos), length(ind_neg)
    M, N          = size(X,1), size(Xvalid,1)

    # auxiliary
    io = open(file, "w+");
    write(io, 1)
    write(io, [nα, nβ, n, npos, nneg, M, N])
    write(io, ind_pos)
    write(io, ind_neg)
    write(io, inv_perm)

    # kernel matrix
    K = Mmap.mmap(io, Matrix{T}, (M, N))
    if partial_fill
        fill_kernelmatrix!(K, kernel, X, Xvalid)
    else
        K .= MLKernels.kernelmatrix(Val(:row), kernel, X, Xvalid)
    end
    K[nα+1:end, :] .*= -1
    Mmap.sync!(K)
    close(io)

    return
end

# -------------------------------------------------------------------------------
# DTest kernel matrix
# -------------------------------------------------------------------------------
function kernelmatrix(M::AbstractModel,
                      Xtrain::AbstractArray,
                      ytrain::AbstractVector{Bool},
                      Xtest::AbstractArray;
                      kernel::Kernel = LinearKernel(1, 0),
                      T::DataType = Float32,)

    X, nα, nβ, = prepare(M, Xtrain, ytrain)
    n = size(Xtest, 1)

    K = T.(MLKernels.kernelmatrix(Val(:row), kernel, X, Xtest))
    K[(nα+1):end, :] .*= -1

    return K, n, nα, nβ
end

function save_kernelmatrix(M::AbstractModel,
                           file::AbstractString,
                           Xtrain::AbstractMatrix,
                           ytrain::AbstractVector{Bool},
                           Xtest::AbstractMatrix;
                           kernel::Kernel = LinearKernel(1, 0),
                           T::DataType = Float32,
                           partial_fill::Bool = false)

    X, nα, nβ, = prepare(M, Xtrain, ytrain)
    n          = size(Xtest, 1)
    M, N       = size(X,1), size(Xtest,1)

    # auxiliary
    io = open(file, "w+");
    write(io, 2)
    write(io, [nα, nβ, n, M, N])

    # kernel matrix
    K = Mmap.mmap(io, Matrix{T}, (M, N))
    if partial_fill
        fill_kernelmatrix!(K, kernel, X, Xtest)
    else
        K .= MLKernels.kernelmatrix(Val(:row), kernel, X, Xtest)
    end
    K[nα+1:end, :] .*= -1
    Mmap.sync!(K)
    close(io)

    return
end

# -------------------------------------------------------------------------------
# Fill function for large matrices
# -------------------------------------------------------------------------------
function fill_kernelmatrix!(K::AbstractMatrix,
                            kernel::Kernel,
                            X::AbstractMatrix,
                            Y::AbstractMatrix;
                            max_chunk_size::Real = 1e6)

    n    = size(X,1)
    Rows = Iterators.partition(1:n, floor(Int, max(max_chunk_size, n)/n))

    ProgressMeter.@showprogress "Kernel matrix calculation in progress: " for rows in Rows
        K[rows, :] .= MLKernels.kernelmatrix(Val(:row), kernel, X[rows,:], Y)
        Mmap.sync!(K)
    end
end

# -------------------------------------------------------------------------------
# Load function
# -------------------------------------------------------------------------------
function load_kernelmatrix(file::AbstractString; T::DataType = Float32)
    io   = open(file, "r");
    type = read(io, Int)

    type in [0,1,2] || @error "Unknown kernel matrix type $(type)"

    if type == 0 || type == 1
        nα, nβ, n, npos, nneg, M, N = [read(io, Int) for k in 1:7]

        ind_pos  = [read(io, Int) for k in 1:npos]
        ind_neg  = [read(io, Int) for k in 1:nneg]
        inv_perm = [read(io, Int) for k in 1:n]
        K        = Mmap.mmap(io, Matrix{T}, (M, N))

        return type, io, (K, n, nα, nβ, ind_pos, ind_neg, inv_perm)
    elseif type == 2
        nα, nβ, n, M, N = [read(io, Int) for k in 1:5]

        K = Mmap.mmap(io, Matrix{T}, (M, N))
        return type, io, (K, n, nα, nβ)
    end
end


text:
using CoordinateTransformations

@testset "operations" begin
    @testset "dimension error test" begin
        a = GeoArray(rand(5, 5))
        b = GeoArray(rand(5, 6))
        @test_throws DimensionMismatch a - b
        @test_throws DimensionMismatch a + b
        @test_throws DimensionMismatch a * b
        @test_throws DimensionMismatch a / b
    end

    @testset "affine error tests" begin
        c = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 -1.0], [375000.03, 380000.03]), "")
        d = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]), "")
        @test_throws DimensionMismatch c - d
        @test_throws DimensionMismatch c + d
        @test_throws DimensionMismatch c * d
        @test_throws DimensionMismatch c / d
    end

    @testset "crs error test" begin
        e = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 -1.0], [375000.03, 380000.03]), "")
        f = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 -1.0], [375000.03, 380000.03]), "PROJCS[\"NAD27 / UTM zone 11N\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.9786982138982,AUTHORITY[\"EPSG\",\"7008\"]],AUTHORITY[\"EPSG\",\"6267\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4267\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-117],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"26711\"]]")
        @test_throws DimensionMismatch e - f
        @test_throws DimensionMismatch e + f
        @test_throws DimensionMismatch e * f
        @test_throws DimensionMismatch e / f
    end

    @testset "operations" begin
        g = GeoArray(rand(5, 5))
        h = GeoArray(rand(5, 5))
        @test typeof(g - h).name.name == :GeoArray
        @test typeof(g + h).name.name == :GeoArray
        @test typeof(g * h).name.name == :GeoArray
        @test typeof(g / h).name.name == :GeoArray
    end
end


text:
# Julia wrapper for header: ../deps/openjpeg/src/lib/openjp3d/openjp3d.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const LIB_OPEN_JP3D_PATH = "../deps/openjpeg/build/bin/libopenjp3d"

function opj_version()
    ccall((:opj_version, LIB_OPEN_JP3D_PATH), Cstring, ())
end

function opj_volume_create(numcmpts::Cint, cmptparms, clrspc::OPJ_COLOR_SPACE)
    ccall((:opj_volume_create, LIB_OPEN_JP3D_PATH), Ptr{opj_volume_t}, (Cint, Ptr{opj_volume_cmptparm_t}, OPJ_COLOR_SPACE), numcmpts, cmptparms, clrspc)
end

function opj_volume_destroy(volume)
    ccall((:opj_volume_destroy, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_volume_t},), volume)
end

function opj_cio_open(cinfo::opj_common_ptr, buffer, length::Cint)
    ccall((:opj_cio_open, LIB_OPEN_JP3D_PATH), Ptr{opj_cio_t}, (opj_common_ptr, Ptr{Cuchar}, Cint), cinfo, buffer, length)
end

function opj_cio_close(cio)
    ccall((:opj_cio_close, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_cio_t},), cio)
end

function cio_tell(cio)
    ccall((:cio_tell, LIB_OPEN_JP3D_PATH), Cint, (Ptr{opj_cio_t},), cio)
end

function cio_seek(cio, pos::Cint)
    ccall((:cio_seek, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_cio_t}, Cint), cio, pos)
end

function opj_set_event_mgr(cinfo::opj_common_ptr, event_mgr, context)
    ccall((:opj_set_event_mgr, LIB_OPEN_JP3D_PATH), Ptr{opj_event_mgr_t}, (opj_common_ptr, Ptr{opj_event_mgr_t}, Ptr{Void}), cinfo, event_mgr, context)
end

function opj_create_decompress(format::OPJ_CODEC_FORMAT)
    ccall((:opj_create_decompress, LIB_OPEN_JP3D_PATH), Ptr{opj_dinfo_t}, (OPJ_CODEC_FORMAT,), format)
end

function opj_destroy_decompress(dinfo)
    ccall((:opj_destroy_decompress, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_dinfo_t},), dinfo)
end

function opj_set_default_decoder_parameters(parameters)
    ccall((:opj_set_default_decoder_parameters, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_dparameters_t},), parameters)
end

function opj_setup_decoder(dinfo, parameters)
    ccall((:opj_setup_decoder, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_dinfo_t}, Ptr{opj_dparameters_t}), dinfo, parameters)
end

function opj_decode(dinfo, cio)
    ccall((:opj_decode, LIB_OPEN_JP3D_PATH), Ptr{opj_volume_t}, (Ptr{opj_dinfo_t}, Ptr{opj_cio_t}), dinfo, cio)
end

function opj_create_compress(format::OPJ_CODEC_FORMAT)
    ccall((:opj_create_compress, LIB_OPEN_JP3D_PATH), Ptr{opj_cinfo_t}, (OPJ_CODEC_FORMAT,), format)
end

function opj_destroy_compress(cinfo)
    ccall((:opj_destroy_compress, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_cinfo_t},), cinfo)
end

function opj_set_default_encoder_parameters(parameters)
    ccall((:opj_set_default_encoder_parameters, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_cparameters_t},), parameters)
end

function opj_setup_encoder(cinfo, parameters, volume)
    ccall((:opj_setup_encoder, LIB_OPEN_JP3D_PATH), Void, (Ptr{opj_cinfo_t}, Ptr{opj_cparameters_t}, Ptr{opj_volume_t}), cinfo, parameters, volume)
end

function opj_encode(cinfo, cio, volume, index)
    ccall((:opj_encode, LIB_OPEN_JP3D_PATH), Cint, (Ptr{opj_cinfo_t}, Ptr{opj_cio_t}, Ptr{opj_volume_t}, Cstring), cinfo, cio, volume, index)
end


text:
# Wisdom & Hernandez version of Kepler solver, but with quartic
# convergence.

function calc_ds_opt(y::T,yp::T,ypp::T,yppp::T) where {T <: Real}
# Computes quartic Newton's update to equation y=0 using first through 3rd derivatives.
# Uses techniques outlined in Murray & Dermott for Kepler solver.
# Rearrange to reduce number of divisions:
num = y*yp
den1 = yp*yp-y*ypp*.5
den12 = den1*den1
den2 = yp*den12-num*.5*(ypp*den1-third*num*yppp)
return -y*den12/den2
end

function solve_kepler!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,
   r0::T,dr0dt::T,s0::T,state::Array{T,1}) where {T <: Real}
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)
# Solves elliptic Kepler's equation for both elliptic and hyperbolic cases.
# Initial guess (if s0 = 0):
r0inv = inv(r0)
beta0inv = inv(beta0)
signb = sign(beta0)
if s0 == zero
  s = h*r0inv
else
  s = copy(s0)
end
s0 = copy(s)
sqb = sqrt(signb*beta0)
y = zero; yp = one
iter = 0
ds = Inf
fac1 = k-r0*beta0
fac2 = r0*dr0dt
while iter == 0 || (abs(ds) > KEPLER_TOL && iter < 10)
  xx = sqb*s
#  sx = sqb*sin(xx)
#  cx = cos(xx)
#  eval(Expr(:call,trig,xx,cx,sx))
  if beta0 > 0
    sx = sin(xx); cx = cos(xx)
  else
    cx = cosh(xx); sx = exp(xx)-cx
  end
  sx *= sqb
# Third derivative:
  yppp = fac1*cx - signb*fac2*sx
# Take derivative:
  yp = (-yppp+ k)*beta0inv
# Second derivative:
  ypp = signb*fac1*beta0inv*sx + fac2*cx
  y  = (-ypp + fac2 +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
  ds = calc_ds_opt(y,yp,ypp,yppp)
  s += ds
  iter +=1
end
#println("iter: ",iter," s0: ",s0," s: ",s," ds: ",ds)
# Since we updated s, need to recompute:
xx = 0.5*sqb*s #; sx = sin(xx) ; cx = cos(xx)
#xx = 0.5*sqb*s; cx = cosh(xx); sx = exp(xx)-cx
#eval(Expr(:call,trig,xx,cx,sx))
if beta0 > 0
  sx = sin(xx); cx = cos(xx)
else
  cx = cosh(xx); sx = exp(xx)-cx
end
#if beta0 > 0
#  cos_sin!(xx,cx,sx)
#else
#  cosh_sinh!(xx,cx,sx)
#end
#println("beta0: ",beta0," xx: ",xx," cx: ",cx," sx: ",sx)
# Now, compute final values:
g1bs = 2.*sx*cx/sqb
g2bs = 2.*signb*sx^2*beta0inv
f = one - k*r0inv*g2bs # eqn (25)
g = r0*g1bs + fac2*g2bs # eqn (27)
#println("f: ",f," g: ",g)
#println("k: ",k," r0inv: ",r0inv," g2bs: ",g2bs," r0: ",r0," g1bs: ",g1bs," fac2: ",fac2)
for j=1:3
# Position is components 2-4 of state:
  state[1+j] = x0[j]*f+v0[j]*g
end
r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
rinv = inv(r)
dfdt = -k*g1bs*rinv*r0inv
dgdt = r0*(one-beta0*g2bs+dr0dt*g1bs)*rinv
for j=1:3
# Velocity is components 5-7 of state:
  state[4+j] = x0[j]*dfdt+v0[j]*dgdt
end
return s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter
end

#function kep_elliptic!(x0::Array{Float64,1},v0::Array{Float64,1},r0::Float64,dr0dt::Float64,k::Float64,h::Float64,beta0::Float64,s0::Float64,state::Array{Float64,1})
function kep_ell_hyp!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,
  beta0::T,s0::T,state::Array{T,1}) where {T <: Real}
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)
# Now, solve for s in elliptical Kepler case:
f = zero; g=zero; dfdt=zero; dgdt=zero; cx=zero;sx=zero;g1bs=zero;g2bs=zero
s=zero; ds = zero; r = zero;rinv=zero; iter=0
if beta0 > zero || beta0 < zero
   s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter = solve_kepler!(h,k,x0,v0,beta0,r0,dr0dt,
    s0,state)
else
  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero; fill!(state,zero); rinv=zero; s=zero; ds=zero; iter = 0
end
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# recompute beta:
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

#function kep_elliptic!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,beta0::T,s0::T,state::Array{T,1},jacobian::Array{T,2})
function kep_ell_hyp!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,beta0::T,s0::T,state::Array{T,1},jacobian::Array{T,2}) where {T <: Real}
# Computes the Jacobian as well
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)
r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in elliptical Kepler case:
f = zero; g=zero; dfdt=zero; dgdt=zero; cx=zero;sx=zero;g1bs=zero;g2bs=zero
s=zero; ds=zero; r = zero;rinv=zero; iter=0
if beta0 > zero || beta0 < zero
#  solve_kepler!(h,k,x0,v0,beta0,s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r0,dr0dt,r,
#    rinv,ds,s0,state,iter)
   s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r,rinv,ds,iter = solve_kepler!(h,k,x0,v0,beta0,r0,dr0dt,
    s0,state)
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v} & q0 = {x0,v0}.
  fill!(jacobian,zero)
  compute_jacobian!(h,k,x0,v0,beta0,s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r0,dr0dt,r,jacobian)
else
  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero; fill!(state,zero); rinv=zero; s=zero; ds=zero; iter = 0
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

#function kep_hyperbolic!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,beta0::T,s0::T,state::Array{T,1})
function kep_hyperbolic!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,beta0::T,s0::T,state::Array{T,1}) where {T <: Real}
# Solves equation (35) from Wisdom & Hernandez for the hyperbolic case.
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)

r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in hyperbolic Kepler case:
if beta0 < zero
#  solve_kepler!()
else
  println("Not hyperbolic",beta0," x0 ",x0)
  r= zero; fill!(state,zero); rinv=zero; s=zero; ds=zero; iter = 0
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

#function kep_hyperbolic!(x0::Array{Float64,1},v0::Array{Float64,1},r0::Float64,dr0dt::Float64,k::Float64,h::Float64,beta0::Float64,s0::Float64,state::Array{Float64,1},jacobian::Array{Float64,2})
function kep_hyperbolic!(x0::Array{T,1},v0::Array{T,1},r0::T,dr0dt::T,k::T,h::T,beta0::T,s0::T,state::Array{T,1},jacobian::Array{T,2}) where {T <: Real}
# Solves equation (35) from Wisdom & Hernandez for the hyperbolic case.
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)

r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in hyperbolic Kepler case:
if beta0 < zero
# Initial guess (if s0 = 0):
  if s0 == zero
    s = h*r0inv
  else
    s = copy(s0)
  end
  s0 = copy(s)
  sqb = sqrt(-beta0)
  y = zero; yp = one
  iter = 0
  ds = Inf
  fac1 = k-r0*beta0
  fac2 = r0*dr0dt
  while iter == 0 || (abs(ds) > KEPLER_TOL && iter < 10)
    xx = sqb*s; cx = cosh(xx); sx = sqb*(exp(xx)-cx)
# Third derivative:
    yppp = fac1*cx + fac2*sx
# Take derivative:
    yp = (-yppp+ k)*beta0inv
# Second derivative:
    ypp = -fac1*beta0inv*sx  + fac2*cx
    y  = (-ypp +fac2 +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
    ds = calc_ds_opt(y,yp,ypp,yppp)
    s += ds
    iter +=1
#    println(iter," ds: ",ds)
  end
#  if iter > 1
#    #println("iter: ",iter," ds/s: ",ds/s0)
#    println(iter," ",s," ",s/s0-1," ds: ",ds)
#  end
  xx = 0.5*sqb*s; cx = cosh(xx); sx = exp(xx)-cx
# Now, compute final values:
  g1bs = 2.0*sx*cx/sqb
  g2bs = -2.0*sx^2*beta0inv
  f = one - k*r0inv*g2bs # eqn (25)
  g = r0*g1bs + fac2*g2bs # eqn (27)
  for j=1:3
    state[1+j] = x0[j]*f+v0[j]*g
  end
  # r = norm(x)
  r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv
  dgdt = r0*(one-beta0*g2bs+dr0dt*g1bs)*rinv
  for j=1:3
# Velocity is components 5-7 of state:
    state[4+j] = x0[j]*dfdt+v0[j]*dgdt
  end
# Now, compute the jacobian:
  fill!(jacobian,zero)
  compute_jacobian!(h,k,x0,v0,beta0,s,f,g,dfdt,dgdt,cx,sx,g1bs,g2bs,r0,dr0dt,r,jacobian)
else
  println("Not hyperbolic",beta0," x0 ",x0)
  r= zero; fill!(state,zero); rinv=zero; s=zero; ds=zero; iter = 0
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

# function compute_jacobian!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,s::T,f::T,g::T,dfdt::T,dgdt::T,cx::T,sx::T,g1::T,g2::T,r0::T,dr0dt::T,r::T,jacobian::Array{T,2})
function compute_jacobian!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,s::T,f::T,g::T,dfdt::T,dgdt::T,cx::T,sx::T,g1::T,g2::T,r0::T,dr0dt::T,r::T,jacobian::Array{T,2}) where {T <: Real}
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v,k} & q0 = {x0,v0,k}.
# Now, compute the Jacobian: (9/18/2017 notes)
#g0 = cx^2-sx^2
zero = convert(typeof(h),0.0); one = convert(typeof(h),1.0)
g0 = one-beta0*g2
g3 = (s-g1)/beta0
dotalpha0 = r0*dr0dt  # unnecessary to divide by r0 for dr0dt & multiply for \dot\alpha_0
absv0 = sqrt(dot(v0,v0))
dsdbeta = (2h-r0*(s*g0+g1)+k/beta0*(s*g0-g1)-dotalpha0*s*g1)/(2beta0*r)
dsdr0 = -(2k/r0^2*dsdbeta+g1/r)
dsda0 = -g2/r
dsdv0 = -2absv0*dsdbeta
dsdk = 2/r0*dsdbeta-g3/r
dbetadr0 = -2k/r0^2
dbetadv0 = -2absv0
dbetadk  = 2/r0
# "p" for partial derivative:
#pxpr0 = zeros(Float64,3); pxpa0=zeros(Float64,3); pxpk=zeros(Float64,3); pxps=zeros(Float64,3); pxpbeta=zeros(Float64,3)
#dxdr0 = zeros(Float64,3); dxda0=zeros(Float64,3); dxdk=zeros(Float64,3); dxdv0 =zeros(Float64,3)
#prvpr0 = zeros(Float64,3); prvpa0=zeros(Float64,3); prvpk=zeros(Float64,3); prvps=zeros(Float64,3); prvpbeta=zeros(Float64,3)
#drvdr0 = zeros(Float64,3); drvda0=zeros(Float64,3); drvdk=zeros(Float64,3); drvdv0=zeros(Float64,3)
#vtmp = zeros(Float64,3); dvdr0 = zeros(Float64,3); dvda0=zeros(Float64,3); dvdv0=zeros(Float64,3); dvdk=zeros(Float64,3)
for i=1:3
  pxpr0[i] = k/r0^2*g2*x0[i]+g1*v0[i]
  pxpa0[i] = g2*v0[i]
  pxpk[i]  = -g2/r0*x0[i]
  pxps[i]  = -k/r0*g1*x0[i]+(r0*g0+dotalpha0*g1)*v0[i]
  pxpbeta[i] = -k/(2beta0*r0)*(s*g1-2g2)*x0[i]+1/(2beta0)*(s*r0*g0-r0*g1+s*dotalpha0*g1-2*dotalpha0*g2)*v0[i]
  prvpr0[i] = k*g1/r0^2*x0[i]+g0*v0[i]
  prvpa0[i] = g1*v0[i]
  prvpk[i] = -g1/r0*x0[i]
  prvps[i] = -k*g0/r0*x0[i]+(-beta0*r0*g1+dotalpha0*g0)*v0[i]
  prvpbeta[i] = -k/(2beta0*r0)*(s*g0-g1)*x0[i]+1/(2beta0)*(-s*r0*beta0*g1+dotalpha0*s*g0-dotalpha0*g1)*v0[i]
end
prpr0 = g0
prpa0 = g1
prpk  = g2
prps = (k-beta0*r0)*g1+dotalpha0*g0
prpbeta = 1/(2beta0)*(s*(k-beta0*r0)*g1+dotalpha0*s*g0-dotalpha0*g1-2k*g2)
for i=1:3
  dxdr0[i] = pxps[i]*dsdr0 + pxpbeta[i]*dbetadr0 + pxpr0[i]
  dxda0[i] = pxps[i]*dsda0 + pxpa0[i]
  dxdv0[i] = pxps[i]*dsdv0 + pxpbeta[i]*dbetadv0
  dxdk[i]  = pxps[i]*dsdk  + pxpbeta[i]*dbetadk + pxpk[i]
  drvdr0[i] = prvps[i]*dsdr0 + prvpbeta[i]*dbetadr0 + prvpr0[i]
  drvda0[i] = prvps[i]*dsda0 + prvpa0[i]
  drvdv0[i] = prvps[i]*dsdv0 + prvpbeta[i]*dbetadv0
  drvdk[i]  = prvps[i]*dsdk  + prvpbeta[i]*dbetadk +prvpk[i]
end
drdr0 = prpr0 + prps*dsdr0 + prpbeta*dbetadr0
drda0 = prpa0 + prps*dsda0
drdv0 = prps*dsdv0 + prpbeta*dbetadv0
drdk  = prpk + prps*dsdk + prpbeta*dbetadk
for i=1:3
  vtmp[i] = dfdt*x0[i]+dgdt*v0[i]
  dvdr0[i] = (drvdr0[i]-drdr0*vtmp[i])/r
  dvda0[i] = (drvda0[i]-drda0*vtmp[i])/r
  dvdv0[i] = (drvdv0[i]-drdv0*vtmp[i])/r
  dvdk[i]  = (drvdk[i] -drdk *vtmp[i])/r
end
# Now, compute Jacobian:
for i=1:3
  jacobian[  i,  i] = f
  jacobian[  i,3+i] = g
  jacobian[3+i,  i] = dfdt
  jacobian[3+i,3+i] = dgdt
  jacobian[  i,7] = dxdk[i]
  jacobian[3+i,7] = dvdk[i]
end
for j=1:3
  for i=1:3
    jacobian[  i,  j] += dxdr0[i]*x0[j]/r0 + dxda0[i]*v0[j]
    jacobian[  i,3+j] += dxdv0[i]*v0[j]/absv0 + dxda0[i]*x0[j]
    jacobian[3+i,  j] += dvdr0[i]*x0[j]/r0 + dvda0[i]*v0[j]
    jacobian[3+i,3+j] += dvdv0[i]*v0[j]/absv0 + dvda0[i]*x0[j]
  end
end
jacobian[7,7]=one
return
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test !Sys.isunix(:Windows)
@test !Sys.islinux(:Windows)
@test Sys.islinux(:Linux)
@test Sys.iswindows(:Windows)
@test Sys.iswindows(:NT)
@test !Sys.iswindows(:Darwin)
@test Sys.isapple(:Darwin)
@test Sys.isapple(:Apple)
@test !Sys.isapple(:Windows)
@test Sys.isunix(:Darwin)
@test Sys.isunix(:FreeBSD)
@test_throws ArgumentError Sys.isunix(:BeOS)
if !Sys.iswindows()
    @test Sys.windows_version() == v"0.0.0"
else
    @test Sys.windows_version() >= v"1.0.0-"
end

@test (@static true ? 1 : 2) === 1
@test (@static false ? 1 : 2) === 2
@test (@static if true 1 end) === 1
@test (@static if false 1 end) === nothing
@test (@static true && 1) === 1
@test (@static false && 1) === false
@test (@static true || 1) === true
@test (@static false || 1) === 1

# test that path variables use correct path delimiters on windows
if Sys.iswindows()
    @test !contains(Base.SYSCONFDIR, "/")
    @test !contains(Base.DATAROOTDIR, "/")
    @test !contains(Base.DOCDIR, "/")
    @test !contains(Base.LIBDIR, "/")
    @test !contains(Base.PRIVATE_LIBDIR, "/")
    @test !contains(Base.INCLUDEDIR, "/")
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
"""
    entanglementplot(state; site=0[, kwargs])

Plot the entanglement spectrum of a given InfiniteMPS.

# Arguments
- `site::Int=0`: mps index for multisite unit cells.
- `expand_symmetry::Logical=false`: add quantum dimension degeneracies.
- `sortby=maximum`: the method of sorting the sectors.
- `sector_margin=1//10`: the amount of whitespace between sectors.
- `sector_formatter=string`: how to convert sectors to strings.
- `kwargs...: other kwargs are passed on to the plotting backend.
"""
function entanglementplot end
@userplot EntanglementPlot

@recipe function f(h::EntanglementPlot; site=0, expand_symmetry=false, sortby=maximum, sector_margin=1//10, sector_formatter=string)

    mps = h.args[1]
    site <= length(mps) ||
        throw(ArgumentError("Not a valid site for the given mps."))

    (_, s, _) = tsvd(mps.CR[site])
    sectors = blocksectors(s)
    spectrum = []
    for sector in sectors
        partial_spectrum = diag(block(s, sector))

        # Duplicate entries according to the quantum dimension.
        if expand_symmetry
            partial_spectrum = repeat(partial_spectrum, dim(sector))
            sort!(partial_spectrum, rev = true)
        end

        push!(spectrum, diag(block(s, sector)))
    end

    if length(spectrum) > 1
        order = sortperm(spectrum, by=sortby, rev=true)
        spectrum = spectrum[order]
        sectors = sectors[order]
    end

    for (i, (partial_spectrum, sector)) in enumerate(zip(spectrum, sectors))
        @series begin
            seriestype := :scatter
            label := sector_formatter(sector)
            n_spectrum = length(partial_spectrum)

            # Put single dot in the middle, or a linear range with padding.
            if n_spectrum == 1
                x = [i+1//2]
            else
                x = range(i + sector_margin, i + 1 - sector_margin,
                    length=n_spectrum)
            end
            return x, partial_spectrum
        end
    end

    title --> "Entanglement Spectrum"
    legend --> false
    grid --> :xy
    widen --> true

    xguide --> "χ = $(dim(domain(s)))"
    xticks --> (1:length(sectors), sector_formatter.(sectors))
    xtickfonthalign --> :center
    xtick_direction --> :out
    xrotation --> 45
    xlims --> (1, length(sectors)+1)

    ylims --> (-Inf, 1+1e-1)
    yscale --> :log10

    return ([])
end

"""
    transferplot(above, below[, sectors[, transferkwargs[, kwargs]]])

Plot the partial transfer matrix spectrum of two InfiniteMPS's.

# Arguments
- `above::InfiniteMPS`: above mps for ``transfer_spectrum``.
- `below::InfiniteMPS`: below mps for ``transfer_spectrum``.
- `sectors=[]`: vector of sectors for which to compute the spectrum.
- `transferkwargs`: kwargs for call to ``transfer_spectrum``.
- `kwargs`: other kwargs are passed on to the plotting backend.
- `thetaorigin=0`: origin of the angle range.
- `sector_formatter=string`: how to convert sectors to strings.
"""
function transferplot end
@userplot TransferPlot
@recipe function f(h::TransferPlot; sectors=nothing, transferkwargs=(;), thetaorigin=0, sector_formatter=string)

    if sectors === nothing
        sectors = [one(TensorKit.sectortype(space(h.args[1], 1)))]
    end

    for sector in sectors
        below = length(h.args) == 1 ? h.args[1] : h.args[2];
        spectrum = transfer_spectrum(h.args[1]; below=below, sector=sector, transferkwargs...)

        @series begin
            yguide --> "r"
            ylims --> (-Inf,1.05)

            xguide --> "θ"
            xlims --> (thetaorigin,thetaorigin+2pi)
            xticks --> range(0, 2pi,length=7)
            xformatter --> x-> "$(rationalize(x/π, tol=0.05))π"
            xwiden --> true
            seriestype := :scatter
            markershape --> :auto
            label := sector_formatter(sector)
            return mod2pi.(angle.(spectrum).+thetaorigin).-thetaorigin, abs.(spectrum)
        end
    end

    title --> "Transfer Spectrum"
    legend --> false
    grid --> :xy
    framestyle --> :zerolines

    return nothing
end


text:
function _slice(f::MOI.VectorOfVariables, slice)
    f_sliced = MOIU.eachscalar(f)[slice]

    if MOI.output_dimension(f_sliced) == 1
        return f_sliced.variables[1]
    else
        return f_sliced
    end
end

function _slice(f::MOI.VectorAffineFunction{T}, slice) where {T}
    f_sliced = MOIU.eachscalar(f)[slice]

    if MOI.output_dimension(f_sliced) == 1
        return MOIU.eachscalar(f)[slice[1]]
    else
        return f_sliced
    end
end


text:
###
# helpers
###

# make sure a passed distance matrix is a square
function check_square(m, msg)
   	n = size(m, 1)
   	if n != size(m, 2)
      		error(msg)
   	end
   	return n
end

"""
    legal_circuit(circuit::AbstractArray{<:Integer})

Check that an array of integers is a valid circuit. A valid circuit over `n` locations has
length `n+1`. The first `n` entries are a permutation of `1, ..., n`, and the `(n+1)`-st entry
is equal to the first entry.
"""
function legal_circuit(circuit)
   	n = length(circuit) - 1
   	return circuit[1] == circuit[end] && sort(circuit[1:(end - 1)]) == 1:n
end

"""
    repetitive_heuristic(distmat::Matrix, heuristic::Function, repetitive_kw::Symbol; keywords...)

For each `i` in `1:n`, run `heuristic` with keyword `repetitive_kw` set to `i`. Return the tuple
`(bestpath, bestcost)`. By default, `repetitive_kw` is `:firstcity`, which is appropriate for the
`farthest_insertion`, `cheapest_insertion`, and `nearest_neighbor` heuristics.

Any keywords passed to `repetitive_heuristic` are passed along to each call of `heuristic`. For
example, `repetitive_heuristic(dm, nearest_neighbor; do2opt = true)` will perform 2-opt for
each of the `n` nearest neighbor paths.

!!! note
    The repetitive heuristic calls are parallelized with threads. For optimum speed make sure
    Julia is running with multiple threads.
"""
function repetitive_heuristic(dm::Matrix{T},
                              heuristic::Function,
                              repetitive_kw = :firstcity;
							  kwargs...) where {T <: Real}
	# call the heuristic with varying starting cities
   	n = size(dm, 1)
   	results_list = Vector{Tuple{Vector{Int},T}}(undef, n)
   	Threads.@threads for i in 1:n
      		results_list[i] = heuristic(dm; kwargs..., repetitive_kw => i)
   	end

   	bestind, bestcost = 1, results_list[1][2]
   	for i in 2:n
      		if results_list[i][2] < bestcost
         			bestind, bestcost = i, results_list[i][2]
      		end
   	end
   	return results_list[bestind]
end

# helper for readable one-line path costs
# optionally specify the bounds for the subpath we want the cost of
# defaults to the whole path
# but when calculating reversed path costs can help to have subpath costs
function pathcost(distmat::Matrix{T}, path::AbstractArray{S},
                  lb::Int = 1, ub::Int = length(path)) where {T <: Real,S <: Integer}
   	cost = zero(T)
   	for i in lb:(ub - 1)
      	@inbounds cost += distmat[path[i], path[i + 1]]
   	end
   	return cost
end

# calculate the cost of reversing part of a path
# cost of walking along the entire path specified but reversing the sequence from revLow to revHigh, inclusive
function pathcost_rev(distmat::Matrix{T}, path::AbstractArray{S},
                      revLow::Int, revHigh::Int) where {T <: Real,S <: Integer}
   	cost = zero(T)
	# if there's an initial unreversed section
   	if revLow > 1
		for i in 1:(revLow - 2)
			@inbounds cost += distmat[path[i], path[i + 1]]
		end
		# from end of unreversed section to beginning of reversed section
      	@inbounds cost += distmat[path[revLow - 1], path[revHigh]]
   	end
	# main reverse section
   	for i in revHigh:-1:(revLow + 1)
      	@inbounds cost += distmat[path[i], path[i - 1]]
   	end
	# if there's an unreversed section after the reversed bit
   	n = length(path)
   	if revHigh < length(path)
		# from end of reversed section back to regular
		@inbounds cost += distmat[path[revLow], path[revHigh + 1]]
		for i in (revHigh + 1):(n - 1)
			@inbounds cost += distmat[path[i], path[i + 1]]
		end
   	end
   	return cost
end

# calculate the cost of two opt reverse
# cost of walking the reversed two edges of two opt (only works for symmetric cases)
function pathcost_two_opt_rev(distmat::Matrix{T}, path::AbstractArray{S},
	revLow::Int, revHigh::Int) where {T <: Real,S <: Integer}
    cost = zero(T)
    @inbounds cost += distmat[path[revLow - 1], path[revHigh]]
    @inbounds cost += distmat[path[revLow], path[revHigh + 1]]
    return cost
end

# calculate the cost of two opt
# cost of walking the normal two edges of two opt (only works for symmetric cases)
function pathcost_two_opt(distmat::Matrix{T}, path::AbstractArray{S},
	revLow::Int, revHigh::Int) where {T <: Real,S <: Integer}
    cost = zero(T)
    @inbounds cost += distmat[path[revLow - 1], path[revLow]]
    @inbounds cost += distmat[path[revHigh], path[revHigh + 1]]
    return cost
end

#Cost of inserting city `k` after index `after` in path `path` with costs `distmat`.
function inscost(k::Int, after::Int, path::AbstractArray{S}, distmat::Matrix{T}) where {T <: Real,S <: Integer}
   	return distmat[path[after], k] +
		  distmat[k, path[after + 1]] -
		  distmat[path[after], path[after + 1]]
end


text:
abstract AbstractTuner

###Factory function, currently implemented are :monitor and :scale
tuner(s::Symbol,p::Int,args...;keyargs...) = _tuner(Val{s},p,args...;keyargs...)

period(tuner::AbstractTuner) = tuner.period
verbose(tuner::AbstractTuner) = tuner.verbose
needstuning(tuner::AbstractTuner,i::Int) = mod(i,period(tuner))==0

############################################################################

abstract AbstractTunerState

###Factory function for tunerstates
tunerstate{T<:AbstractFloat}(t::AbstractTuner,nburnin::Int,::Type{T}) = _tunerstate(t,nburnin,T)

accepted(state::AbstractTunerState) = state.accepted
proposed(state::AbstractTunerState) = state.proposed
rate(state::AbstractTunerState) = accepted(state)./proposed(state)
total(state::AbstractTunerState) = sum(proposed(state))
index(state::AbstractTunerState) = state.index
numtunesteps(state::AbstractTunerState) = length(state.accepted)

current(state::AbstractTunerState) = (a = state.accepted[state.index] ; p = state.proposed[state.index] ; tuple(a,p))

tune(tuner::AbstractTuner,state::AbstractTunerState) = ()

@inline function accepted!(state::AbstractTunerState,indicator::AbstractIndicatorMatrix)
    if state.index <= numtunesteps(state)
        state.accepted[state.index] += accepted(indicator)
        state.proposed[state.index] += numsamples(indicator)
    else
        warn("TunerState is full. Additional tuning steps cannot be stored")
    end
end

@inline nextindex!(state::AbstractTunerState) = (state.index+=1)


text:
module Compute

include("module_common.jl")

const _libcloud_compute_types = PyCall.PyNULL()
const _libcloud_compute_providers = PyCall.PyNULL()
const _libcloud_compute_base = PyCall.PyNULL()
const _libcloud_compute_deployment = PyCall.PyNULL()

immutable NodeDriver
    o::PyObject

    function NodeDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_compute_providers[:get_driver](provider)
        NodeDriver(cls(args...; kwargs...))
    end
    NodeDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::NodeDriver) = print(io, c.o[:__str__]())
features(driver::NodeDriver) = driver.o[:features]

immutable NodeSize
    o::PyObject

    id::Compat.String
    name::Compat.String
    ram::Nullable{Int}
    disk::Nullable{Int}
    bandwidth::Nullable{Int}
    price::Nullable{Float64}
    extra::Dict
    driver::NodeDriver

    function NodeSize(o::PyObject)
        new(o, o[:id], o[:name], o[:ram], o[:disk], o[:bandwidth], o[:price], o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::NodeSize) = o.o
convert(::Type{NodeSize}, o::PyObject) = NodeSize(o)
show(io::IO, o::NodeSize) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "NodeSize") NodeSize

immutable NodeImage
    o::PyObject

    id::Compat.String
    name::Compat.String
    extra::Dict
    driver::NodeDriver

    function NodeImage(o::PyObject)
        new(o, o[:id], o[:name], o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::NodeImage) = o.o
convert(::Type{NodeImage}, o::PyObject) = NodeImage(o)
show(io::IO, o::NodeImage) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "NodeImage") NodeImage

immutable Node
    o::PyObject

    id::Compat.String
    name::Compat.String
    public_ips::Vector
    private_ips::Vector
    size::Nullable{Any}    # can be NodeSize, String, Nothing
    image::Nullable{Any}   # can be NodeImage, String, Nothing
    created_at::Nullable{DateTime}
    extra::Dict
    driver::NodeDriver

    function Node(o::PyObject)
        created_at = (o[:created_at] == nothing) ? Nullable{DateTime}() : Nullable(DateTime(o[:created_at]))
        size = (o[:size] == nothing) ? Nullable{Any}() : Nullable{Any}(o[:size])
        image = (o[:image] == nothing) ? Nullable{Any}() : Nullable{Any}(o[:image])
        new(o, o[:id], o[:name], o[:public_ips], o[:private_ips], size, image, created_at, o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::Node) = o.o
convert(::Type{Node}, o::PyObject) = Node(o)
show(io::IO, o::Node) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "Node") Node
state(n::Node) = n.o[:state]

immutable NodeLocation
    o::PyObject

    id::Compat.String
    name::Compat.String
    country::Compat.String
    driver::NodeDriver

    function NodeLocation(o::PyObject)
        new(o, o[:id], o[:name], o[:country], NodeDriver(o[:driver]))
    end
end
PyObject(o::NodeLocation) = o.o
convert(::Type{NodeLocation}, o::PyObject) = NodeLocation(o)
show(io::IO, o::NodeLocation) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "NodeLocation") NodeLocation

immutable NodeAuthSSHKey
    o::PyObject

    pubkey::Compat.String

    function NodeAuthSSHKey(o::PyObject)
        new(o, o[:pubkey])
    end
    function NodeAuthSSHKey(pubkey::Compat.String)
        new(_libcloud_compute_base["NodeAuthSSHKey"](pubkey), pubkey)
    end
end
PyObject(o::NodeAuthSSHKey) = o.o
convert(::Type{NodeAuthSSHKey}, o::PyObject) = NodeAuthSSHKey(o)
show(io::IO, o::NodeAuthSSHKey) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "NodeAuthSSHKey") NodeAuthSSHKey

immutable NodeAuthPassword
    o::PyObject

    password::Compat.String
    generated::Bool

    function NodeAuthPassword(o::PyObject)
        new(o, o[:password], o[:generated])
    end
end
PyObject(o::NodeAuthPassword) = o.o
convert(::Type{NodeAuthPassword}, o::PyObject) = NodeAuthPassword(o)
show(io::IO, o::NodeAuthPassword) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "NodeAuthPassword") NodeAuthPassword

immutable StorageVolume
    o::PyObject

    id::Compat.String
    name::Compat.String
    size::Int
    state::Compat.String
    extra::Dict
    driver::NodeDriver

    function StorageVolume(o::PyObject)
        new(o, o[:id], o[:name], o[:size], o[:state], o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::StorageVolume) = o.o
convert(::Type{StorageVolume}, o::PyObject) = StorageVolume(o)
show(io::IO, o::StorageVolume) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "StorageVolume") StorageVolume

immutable VolumeSnapshot
    o::PyObject

    id::Compat.String
    size::Int
    created::DateTime
    state::Compat.String
    extra::Dict
    driver::NodeDriver

    function VolumeSnapshot(o::PyObject)
        created = o[:created]
        new(o, o[:id], o[:size], created, o[:state], o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::VolumeSnapshot) = o.o
convert(::Type{VolumeSnapshot}, o::PyObject) = VolumeSnapshot(o)
show(io::IO, o::VolumeSnapshot) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "VolumeSnapshot") VolumeSnapshot

immutable KeyPair
    o::PyObject

    name::Compat.String
    fingerprint::Compat.String
    public_key::Nullable{Compat.String}
    private_key::Nullable{Compat.String}
    extra::Dict
    driver::NodeDriver

    function KeyPair(o::PyObject)
        new(o, o[:name], o[:fingerprint], o[:public_key], o[:private_key], o[:extra], NodeDriver(o[:driver]))
    end
end
PyObject(o::KeyPair) = o.o
convert(::Type{KeyPair}, o::PyObject) = KeyPair(o)
show(io::IO, o::KeyPair) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_base, "KeyPair") KeyPair

immutable SSHKeyDeployment
    o::PyObject

    key::Compat.String

    function SSHKeyDeployment(o::PyObject)
        new(o, o[:key])
    end
    function SSHKeyDeployment(key::Compat.String)
        new(_libcloud_compute_deployment["SSHKeyDeployment"](key), key)
    end
end
PyObject(o::SSHKeyDeployment) = o.o
convert(::Type{SSHKeyDeployment}, o::PyObject) = SSHKeyDeployment(o)
show(io::IO, o::SSHKeyDeployment) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_deployment, "SSHKeyDeployment") SSHKeyDeployment

immutable FileDeployment
    o::PyObject

    source::Compat.String
    target::Compat.String

    function FileDeployment(o::PyObject)
        new(o, o[:source], o[:target])
    end
    function FileDeployment(source::Compat.String, target::Compat.String)
        new(_libcloud_compute_deployment["FileDeployment"](source, target), source, target)
    end
end
PyObject(o::FileDeployment) = o.o
convert(::Type{FileDeployment}, o::PyObject) = FileDeployment(o)
show(io::IO, o::FileDeployment) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_deployment, "FileDeployment") FileDeployment

immutable ScriptDeployment
    o::PyObject

    script::Compat.String
    args::Vector{Compat.String}
    name::Compat.String
    delete::Bool

    function ScriptDeployment(o::PyObject)
        new(o, o[:script], o[:args], o[:name], o[:delete])
    end
    function ScriptDeployment(script::Compat.String, args...; kwargs...)
        ScriptDeployment(_libcloud_compute_deployment["ScriptDeployment"](script, args...; kwargs...))
    end
end
PyObject(o::ScriptDeployment) = o.o
convert(::Type{ScriptDeployment}, o::PyObject) = ScriptDeployment(o)
show(io::IO, o::ScriptDeployment) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_deployment, "ScriptDeployment") ScriptDeployment

immutable ScriptFileDeployment
    o::PyObject

    script::Compat.String
    args::Vector{Compat.String}
    name::Compat.String
    delete::Bool

    function ScriptFileDeployment(o::PyObject)
        new(o, o[:script], o[:args], o[:name], o[:delete])
    end
    function ScriptFileDeployment(script::Compat.String, args...; kwargs...)
        ScriptDeployment(_libcloud_compute_deployment["ScriptFileDeployment"](script, args...; kwargs...))
    end
end
PyObject(o::ScriptFileDeployment) = o.o
convert(::Type{ScriptFileDeployment}, o::PyObject) = ScriptFileDeployment(o)
show(io::IO, o::ScriptFileDeployment) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_deployment, "ScriptFileDeployment") ScriptFileDeployment

immutable MultiStepDeployment
    o::PyObject

    function MultiStepDeployment(o::PyObject)
        new(o)
    end
    function MultiStepDeployment(add::Vector=[])
        if isempty(add)
            MultiStepDeployment(_libcloud_compute_deployment["MultiStepDeployment"]())
        else
            MultiStepDeployment(_libcloud_compute_deployment["MultiStepDeployment"](add=add))
        end
    end
end
PyObject(o::MultiStepDeployment) = o.o
convert(::Type{MultiStepDeployment}, o::PyObject) = MultiStepDeployment(o)
show(io::IO, o::MultiStepDeployment) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_compute_deployment, "MultiStepDeployment") MultiStepDeployment

const _nodedriver_fns = [
    # Node management methods
    :list_nodes, :list_sizes, :list_locations, :create_node, :deploy_node, :reboot_node, :destroy_node, :wait_until_running,
    # Volume and snapshot management methods
    :list_volumes, :list_volume_snapshots, :create_volume, :create_volume_snapshot, :attach_volume, :detach_volume, :destroy_volume, :destroy_volume_snapshot,
    # Image management methods
    :list_images, :create_image, :delete_image, :get_image, :copy_image,
    # SSH key pair management methods
    :list_key_pairs, :get_key_pair, :create_key_pair, :import_key_pair_from_string, :import_key_pair_from_file, :delete_key_pair
]

const _uuid_fns = [ :get_uuid ]
const _node_fns = [ :reboot, :destroy ]
const _storage_volume_fns = [ :list_snapshots, :attach, :detach, :snapshot, :destroy ]
const _volume_snapshot_fns = [ :destroy ]
const _deployment_fns = [ :run ]
const _multi_step_deployment_fns = [ :add ]
const _script_deployment_accessor_fns = [ :stdout, :stdin, :exit_status ]
const _multi_step_deployment_accessor_fns = [ :steps ]

for f in _uuid_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_base, "UuidMixin", $sf) $(f)(n::Node, args...; kwargs...) = n.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_base, "UuidMixin", $sf) $(f)(n::NodeSize, args...; kwargs...) = n.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_base, "UuidMixin", $sf) $(f)(n::NodeImage, args...; kwargs...) = n.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_base, "UuidMixin", $sf) $(f)(n::StorageVolume, args...; kwargs...) = n.o[$(sf)](args...; kwargs...)
end

for f in union(Set(_nodedriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_base, "NodeDriver", $sf) $(f)(driver::NodeDriver, args...; kwargs...) = driver.o[$(sf)](args...; kwargs...)
end

for f in _node_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_base, "Node", $sf) $(f)(node::Node, args...; kwargs...) = node.o[$(sf)](args...; kwargs...)
end

for f in _storage_volume_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_base, "StorageVolume", $sf) $(f)(vol::StorageVolume, args...; kwargs...) = vol.o[$(sf)](args...; kwargs...)
end

for f in _volume_snapshot_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_base, "VolumeSnapshot", $sf) $(f)(vol::VolumeSnapshot, args...; kwargs...) = vol.o[$(sf)](args...; kwargs...)
end

for f in _deployment_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "SSHKeyDeployment", $sf) $(f)(d::SSHKeyDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "FileDeployment", $sf) $(f)(d::FileDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "ScriptDeployment", $sf) $(f)(d::ScriptDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "ScriptDeployment", $sf) $(f)(d::ScriptFileDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "MultiStepDeployment", $sf) $(f)(d::MultiStepDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
end

for f in _multi_step_deployment_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_compute_deployment, "MultiStepDeployment", $sf) $(f)(d::MultiStepDeployment, args...; kwargs...) = d.o[$(sf)](args...; kwargs...)
end

for f in _script_deployment_accessor_fns
    sf = string(f)
    @eval $(f)(d::ScriptDeployment) = d.o[sf]
    @eval $(f)(d::ScriptFileDeployment) = d.o[sf]
end

for f in _multi_step_deployment_accessor_fns
    sf = string(f)
    @eval $(f)(d::MultiStepDeployment) = d.o[sf]
end

# Initialize cloud compute
function __init__()
    copy!(_libcloud_compute_types, pyimport("libcloud.compute.types"))
    copy!(_libcloud_compute_providers, pyimport("libcloud.compute.providers"))
    copy!(_libcloud_compute_base, pyimport("libcloud.compute.base"))
    copy!(_libcloud_compute_deployment, pyimport("libcloud.compute.deployment"))

    global const ComputeProvider = pywrap(_libcloud_compute_types[:Provider])
    global const NodeState = pywrap(_libcloud_compute_types[:NodeState])
    global const Architecture = pywrap(_libcloud_compute_types[:Architecture])
    global const StorageVolumeState = pywrap(_libcloud_compute_types[:StorageVolumeState])
    global const VolumeSnapshotState = pywrap(_libcloud_compute_types[:VolumeSnapshotState])

    _map_types(_libcloud_compute_base, (Node, NodeSize, NodeImage, NodeLocation, NodeAuthSSHKey, NodeAuthPassword, StorageVolume, VolumeSnapshot, KeyPair))
    _map_types(_libcloud_compute_deployment, (SSHKeyDeployment, FileDeployment, ScriptDeployment, ScriptFileDeployment, MultiStepDeployment))
end

# types
export ComputeProvider, NodeDriver, NodeSize, NodeImage, Node, NodeLocation, NodeAuthSSHKey, NodeAuthPassword, StorageVolume, VolumeSnapshot, KeyPair
export NodeState, Architecture, StorageVolumeState, VolumeSnapshotState
export SSHKeyDeployment, FileDeployment, ScriptDeployment, ScriptFileDeployment, MultiStepDeployment, stdout, stderr, exit_status, steps
# Node management methods
export list_nodes, list_sizes, list_locations, create_node, deploy_node, reboot_node, destroy_node, wait_until_running
export features
# Volume and snapshot management methods
export list_volumes, list_volume_snapshots, create_volume, create_volume_snapshot, attach_volume, detach_volume, destroy_volume, destroy_volume_snapshot
# Image management methods
export list_images, create_image, delete_image, get_image, copy_image
# SSH key pair management methods
export list_key_pairs, get_key_pair, create_key_pair, import_key_pair_from_string, import_key_pair_from_file, delete_key_pair

export get_uuid
# node methods
export reboot, destroy, state
export list_snapshots, attach, detach, snapshot, destroy
export destroy

end # module Compute


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined formats.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export unicode, ascii_dots, ascii_rounded, borderless, compact, markdown, mysql, simple, unicode_rounded, matrix

const unicode = TextFormat()

const ascii_dots = TextFormat(
    up_right_corner     = '.',
    up_left_corner      = '.',
    bottom_left_corner  = ':',
    bottom_right_corner = ':',
    up_intersection     = '.',
    left_intersection   = ':',
    right_intersection  = ':',
    middle_intersection = ':',
    bottom_intersection = ':',
    column              = ':',
    left_border         = ':',
    right_border        = ':',
    row                 = '.'
)

const ascii_rounded = TextFormat(
    up_right_corner     = '.',
    up_left_corner      = '.',
    bottom_left_corner  = '\'',
    bottom_right_corner = ''',
    up_intersection     = '.',
    left_intersection   = ':',
    right_intersection  = ':',
    middle_intersection = '+',
    bottom_intersection = ''',
    column              = '|',
    left_border         = '|',
    right_border        = '|',
    row                 = '-'
)

const borderless = TextFormat(
    up_right_corner     = ' ',
    up_left_corner      = ' ',
    bottom_left_corner  = ' ',
    bottom_right_corner = ' ',
    up_intersection     = ' ',
    left_intersection   = ' ',
    right_intersection  = ' ',
    middle_intersection = ' ',
    bottom_intersection = ' ',
    column              = ' ',
    left_border         = ' ',
    right_border        = ' ',
    row                 = ' ',
    top_line            = false,
    bottom_line         = false
)

const compact = TextFormat(
    up_right_corner     = ' ',
    up_left_corner      = ' ',
    bottom_left_corner  = ' ',
    bottom_right_corner = ' ',
    up_intersection     = ' ',
    left_intersection   = ' ',
    right_intersection  = ' ',
    middle_intersection = ' ',
    bottom_intersection  = ' ',
    column              = ' ',
    left_border         = ' ',
    right_border        = ' ',
    row                 = '-'
   )

const markdown = TextFormat(
    left_intersection   = '|',
    right_intersection  = '|',
    middle_intersection = '|',
    column              = '|',
    left_border         = '|',
    right_border        = '|',
    row                 = '-',
    top_line            = false,
    bottom_line         = false,
)

const mysql = TextFormat(
    up_right_corner     = '+',
    up_left_corner      = '+',
    bottom_left_corner  = '+',
    bottom_right_corner = '+',
    up_intersection     = '+',
    left_intersection   = '+',
    right_intersection  = '+',
    middle_intersection = '+',
    bottom_intersection = '+',
    column              = '|',
    left_border         = '|',
    right_border        = '|',
    row                 = '-'
   )

const simple = TextFormat(
    up_right_corner     = '=',
    up_left_corner      = '=',
    bottom_left_corner  = '=',
    bottom_right_corner = '=',
    up_intersection     = ' ',
    left_intersection   = '=',
    right_intersection  = '=',
    middle_intersection = ' ',
    bottom_intersection  = ' ',
    column              = ' ',
    left_border         = ' ',
    right_border        = ' ',
    row                 = '='
)

const unicode_rounded = TextFormat(
    up_right_corner     = '╮',
    up_left_corner      = '╭',
    bottom_left_corner  = '╰',
    bottom_right_corner = '╯',
    up_intersection     = '┬',
    left_intersection   = '├',
    right_intersection  = '┤',
    middle_intersection = '┼',
    bottom_intersection = '┴',
    column              = '│',
    left_border         = '│',
    right_border        = '│',
    row                 = '─'
)

const matrix = TextFormat(
    up_intersection     = ' ',
    left_intersection   = '│',
    right_intersection  = '│',
    middle_intersection = ' ',
    bottom_intersection = ' ',
    left_border         = '│',
    right_border        = '│',
    column              = ' ',
    row                 = ' '
)


text:
# PowellSingular function
#
# source: Nielsen, Hans Bruun and others, Damping parameter in Marquardt's method, 1999, IMM
#
function PowellSingular()
NLS_ForwardDiff(4,4) do θ
        T = eltype(θ)

        T[θ[1]+10*θ[2],
          sqrt(5)*(θ[3]-θ[4]),
          (θ[2]-2*θ[3])^2,
          sqrt(10)*(θ[1]-θ[4])^2]

    end
end


text:
## Test higher criticism ##
module Test_HigherCriticism

using MultipleTesting
using Test

@testset "Higher criticism" begin

    pval1 = [0.00001, 0.0001, 0.002, 0.002, 0.1, 0.1, 0.2]
    hcs1_ref = [1.080048, 1.672734, 3.044358, 3.044358, 5.724654, 5.724654, 6.048691]
    hcv1_ref = 0.2

    pval2 = [0, 1e-10, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0]
    hcs2_ref = [1.032796, 1.511858, 1.921538, 2.309400, 2.696791, 3.098304, 3.526862, 3.992000,
                3.729249, 3.511505, 3.344031, 3.233162, 3.202563, 3.326087, 2.272150, 0.0]
    hcv2_ref = 0.001

    @testset "Higher criticism scores" begin

        hcs1 = estimate(PValues(pval1), HigherCriticismScores())
        @test isapprox(hcs1, hcs1_ref, atol = 1e-5)

        hcs2 = estimate(PValues(pval2), HigherCriticismScores())
        @test isapprox(hcs2, hcs2_ref, atol = 1e-5)

        @test_throws MethodError estimate(PValues(pval1))
        @test_throws MethodError estimate(HigherCriticismScores())

        # currently not supported
        @test_throws MethodError estimate(pval1, HigherCriticismScores())
        @test_throws MethodError estimate(PValues(pval1), HigherCriticismScores)

    end

    @testset "Higher criticism threshold" begin

        hcv1 = estimate(PValues(pval1), HigherCriticismThreshold())
        @test hcv1 == hcv1_ref

        hcv2 = estimate(PValues(pval2), HigherCriticismThreshold())
        @test hcv2 == hcv2_ref

        @test_throws MethodError estimate(PValues(pval1))
        @test_throws MethodError estimate(HigherCriticismThreshold())

        # currently not supported
        @test_throws MethodError estimate(pval1, HigherCriticismThreshold())
        @test_throws MethodError estimate(PValues(pval1), HigherCriticismThreshold)

    end

end

end


text:
@testset "Collections" begin
  @testset "PointSet" begin
    @test !isgrid(PointSet{1,T})
    @test !isgrid(PointSet{2,T})
    @test !isgrid(PointSet{3,T})

    pset = PointSet(rand(P1, 100))
    @test embeddim(pset) == 1
    @test coordtype(pset) == T
    @test nelements(pset) == 100
    @test eltype(pset) <: P1

    pset = PointSet(rand(P2, 100))
    @test embeddim(pset) == 2
    @test coordtype(pset) == T
    @test nelements(pset) == 100
    @test eltype(pset) <: P2

    pset = PointSet(rand(P3, 100))
    @test embeddim(pset) == 3
    @test coordtype(pset) == T
    @test nelements(pset) == 100
    @test eltype(pset) <: P3

    pset1 = PointSet([P3(1,2,3), P3(4,5,6)])
    pset2 = PointSet(P3(1,2,3), P3(4,5,6))
    pset3 = PointSet([T.((1,2,3)), T.((4,5,6))])
    pset4 = PointSet(T.((1,2,3)), T.((4,5,6)))
    pset5 = PointSet([T[1,2,3], T[4,5,6]])
    pset6 = PointSet(T[1,2,3], T[4,5,6])
    pset7 = PointSet(T[1 4; 2 5; 3 6])
    @test pset1 == pset2 == pset3 == pset4 ==
          pset5 == pset6 == pset7
    for pset in [pset1, pset2, pset3, pset4,
                pset5, pset6, pset7]
      @test embeddim(pset) == 3
      @test coordtype(pset) == T
      @test nelements(pset) == 2
      @test pset[1] == P3(1,2,3)
      @test pset[2] == P3(4,5,6)
    end

    pset = PointSet(P2[(0,0), (1,0), (0,1)])
    @test centroid(pset) == P2(1/3, 1/3)

    pset = PointSet(P2[(1,0), (0,1)])
    @test nelements(pset) == 2
    @test centroid(pset, 1) == P2(1, 0)
    @test centroid(pset, 2) == P2(0, 1)

    pset = PointSet(P2[(1,0), (0,1)])
    @test sprint(show, pset) == "2 PointSet{2,$T}"
    if T == Float32
      @test sprint(show, MIME"text/plain"(), pset) == "2 PointSet{2,Float32}\n  └─Point(1.0f0, 0.0f0)\n  └─Point(0.0f0, 1.0f0)"
    elseif T == Float64
      @test sprint(show, MIME"text/plain"(), pset) == "2 PointSet{2,Float64}\n  └─Point(1.0, 0.0)\n  └─Point(0.0, 1.0)"
    end

    if visualtests
      ps1 = P1[0,1,1,2]
      ps2 = P2[(0,1),(1,1),(1,2),(2,3)]
      ps3 = P3[(0,1,1),(1,1,2),(1,2,3),(2,3,5)]
      @test_reference "data/pset-1D-$T.png" plot(PointSet(ps1))
      @test_reference "data/pset-2D-$T.png" plot(PointSet(ps2))
      @test_reference "data/pset-3D-$T.png" plot(PointSet(ps3))
      @test_reference "data/pset-1D-$T-data.png" plot(PointSet(ps1),1:10)
      @test_reference "data/pset-2D-$T-data.png" plot(PointSet(ps2),1:10)
      @test_reference "data/pset-3D-$T-data.png" plot(PointSet(ps3),1:10)
    end
  end

  @testset "GeometrySet" begin
    @test !isgrid(GeometrySet)

    s = Segment(P2(0,0), P2(1,1))
    t = Triangle(P2(0,0), P2(1,0), P2(0,1))
    p = PolyArea(P2[(0,0), (1,0), (1,1), (0,1), (0,0)])
    gset = GeometrySet([s, t, p])
    @test [centroid(gset, i) for i in 1:3] == P2[(1/2,1/2), (1/3,1/3), (1/2,1/2)]
  end
end