text:
##

tbl = DataFrame([1 4]', [:latency])
X = ones(size(tbl))
shouldBeNeg = zeros(4, 4)
shouldBeNeg[1, :] = [1, 0, 0, 1]
shouldBeNeg[2, :] = [0, 1, 0, 0]
shouldBeNeg[3, :] = [0, 0, 1, 0]
shouldBeNeg[4, :] = [0, 0, 0, 1]

shouldBePos = zeros(4, 4)
shouldBePos[1, :] = [0, 0, 0, 0]
shouldBePos[2, :] = [1, 0, 0, 0]
shouldBePos[3, :] = [0, 1, 0, 0]
shouldBePos[4, :] = [0, 0, 1, 0]

## test negative
basisfunction = firbasis(τ = (-3, 0), sfreq = 1, name = "testing")
timeexpandterm = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency);
Xdc = Unfold.time_expand(X, timeexpandterm, tbl)
@test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBeNeg, atol = 1e-15))

## Test Positive only
basisfunction = firbasis(τ = (1, 4), sfreq = 1, name = "testing")
timeexpandterm = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency);
Xdc = Unfold.time_expand(X, timeexpandterm, tbl)
println(Matrix(Xdc))

@test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBePos, atol = 1e-15))

# customized eventfields
tbl2 = tbl = DataFrame([1 4]', [:onset])

timeexpandterm_latency = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction);
timeexpandterm_onset =
    Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, eventfields = [:onset]);
Xdc = Unfold.time_expand(X, timeexpandterm_onset, tbl)
@test_throws ArgumentError Unfold.time_expand(X, timeexpandterm_latency, tbl)

## combining designmatrices
tbl = DataFrame([1 4]', [:latency])
X = ones(size(tbl))
basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1")
basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2")
f = @formula 0 ~ 1
Xdc1 = designmatrix(UnfoldLinearModel, f, tbl, basisfunction1)
Xdc2 = designmatrix(UnfoldLinearModel, f, tbl .+ 1, basisfunction2)

Xdc = Xdc1 + Xdc2
@test size(Xdc.Xs, 2) == size(Xdc1.Xs, 2) + size(Xdc2.Xs, 2)
@test length(Xdc.events) == 2

Xdc_3 = Xdc1 + Xdc2 + Xdc2

@test size(Xdc_3.Xs, 2) == size(Xdc1.Xs, 2) + 2*size(Xdc2.Xs, 2)
@test length(Xdc_3.events) == 3

basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1")
basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2")

tbl = DataFrame(
    [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]',
    [:latency, :subject, :item],
)
tbl2 = DataFrame(
    [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]',
    [:latency, :subject, :itemB],
)
y = Float64.([collect(range(1, stop = 100))...])'
transform!(tbl, :subject => categorical => :subject)
transform!(tbl2, :itemB => categorical => :itemB)
transform!(tbl, :item => categorical => :item)
#tbl.itemB = tbl.item
f3 = @formula 0 ~ 1 + (1 | subject) + (1 | item)
f4 = @formula 0 ~ 1 + (1 | itemB)
f4_wrong = @formula 0 ~ 1 + (1 | item)
Xdc3 = designmatrix(UnfoldLinearMixedModel, f3, tbl, basisfunction1)
Xdc4 = designmatrix(UnfoldLinearMixedModel, f4, tbl2, basisfunction2)
Xdc4_wrong = designmatrix(UnfoldLinearMixedModel, f4_wrong, tbl, basisfunction2)

Xdc = Xdc3 + Xdc4;
@test typeof(Xdc.Xs[1]) <: SparseArrays.SparseMatrixCSC
@test size(Xdc.Xs[1], 2) == size(Xdc3.Xs[1], 2) + size(Xdc4.Xs[1], 2)
@test length(Xdc.Xs) == 4 # one FeMat  + 3 ReMat
@test_throws String Xdc3 + Xdc4_wrong
uf = UnfoldLinearMixedModelContinuousTime(Dict(), Xdc, [])
Unfold.fit!(uf, y);

Xdc = Xdc3 + Xdc4;
df = fit(
    UnfoldLinearMixedModelContinuousTime,
    Xdc,
    rand(1, size(Unfold.modelmatrix(Xdc)[1], 1)),
)
@test size(Unfold.coef(df), 2) == 17

## Speedtest
if 1 == 0
    x = collect(range(10, stop = 10000, step = 10))
    tbl = DataFrame(reshape(x, 1000, 1), [:latency])
    X = ones(size(tbl, 1), 3) .* [1, 2, 3]'

    basisfunction = firbasis(τ = (0, 1), sfreq = 100, name = "test")
    term = Unfold.TimeExpandedTerm(Term, basisfunction, :latency)
    @time Xdc = Matrix(Unfold.time_expand(X, term, tbl))

end

## Test equalizeReMatLengths
bf1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1")
bf2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2")

tbl1 = DataFrame(
    [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]',
    [:latency, :subject, :item],
)
tbl2 = DataFrame(
    [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]',
    [:latency, :subject, :itemB],
)

transform!(tbl1, :subject => categorical => :subject)
transform!(tbl1, :item => categorical => :item)
transform!(tbl2, :itemB => categorical => :itemB)
#tbl.itemB = tbl.item
f1 = @formula 0 ~ 1 + (1 | subject) + (1 | item)
f2 = @formula 0 ~ 1 + (1 | itemB)

form = apply_schema(f1, schema(f1, tbl1), MixedModels.LinearMixedModel)
form = Unfold.apply_basisfunction(form, bf1, nothing)
X1 = modelcols.(form.rhs, Ref(tbl1))

form = apply_schema(f2, schema(f2, tbl2), MixedModels.LinearMixedModel)
form = Unfold.apply_basisfunction(form, bf2, nothing)
X2 = modelcols.(form.rhs, Ref(tbl2))

# no missmatch, shouldnt change anything then
X = deepcopy(X1[2:end])
Unfold.equalizeReMatLengths!(X)
@test all([x[1] for x in size.(X)] .== 48)

X = (deepcopy(X1[2:end])..., deepcopy(X2[2:end])...)
@test !all([x[1] for x in size.(X)] .== 48) # not alllenghts the same
Unfold.equalizeReMatLengths!(X)
@test all([x[1] for x in size.(X)] .== 49) # now all lengths the same :-)

# Test  changeReMatSize & changeMatSize

X = deepcopy(X2[2])
@test size(X)[1] == 49
Unfold.changeReMatSize!(X, 52)
@test size(X)[1] == 52

X = deepcopy(X2[2])
@test size(X)[1] == 49
Unfold.changeReMatSize!(X, 40)
@test size(X)[1] == 40

X = (deepcopy(X1)..., deepcopy(X2[2:end])...)
@test size(X[1])[1] == 48
@test size(X[2])[1] == 48
@test size(X[3])[1] == 48
@test size(X[4])[1] == 49
XA, XB = Unfold.changeMatSize!(52, X[1], X[2:end])
@test size(XA)[1] == 52
@test size(XB)[1] == 52

XA, XB = Unfold.changeMatSize!(40, X[1], X[2:end])
@test size(XA)[1] == 40
@test size(XB)[1] == 40

XA, XB = Unfold.changeMatSize!(30, Matrix(X[1]), X[2:end])
@test size(XA)[1] == 30
@test size(XB)[1] == 30

#----- Some LinearMixedModel tests

data, evts = loadtestdata("testCase3", dataPath = (@__DIR__) * "/data") #
evts.subject = categorical(evts.subject)

f_zc = @formula 0 ~ 1 + condA + condB + zerocorr(1 + condA + condB | subject)
basisfunction = firbasis(τ = (-0.1, 0.1), sfreq = 10, name = "ABC")
Xdc_zc = designmatrix(UnfoldLinearMixedModel, f_zc, evts, basisfunction)

@test length(Xdc_zc.Xs[2].inds) == 9
f = @formula 0 ~ 1 + condA + condB + (1 + condA + condB | subject)
Xdc = designmatrix(UnfoldLinearMixedModel, f, evts, basisfunction)
@test length(Xdc.Xs[2].inds) == (9 * 9 + 9) / 2

# test bug with not sequential subjects
evts_nonseq = copy(evts)
evts_nonseq = evts_nonseq[.!(evts_nonseq.subject .== 2), :]
Xdc_nonseq = designmatrix(UnfoldLinearMixedModel, f_zc, evts_nonseq, basisfunction)
# This used to lead to problems here:
fit(UnfoldLinearMixedModel, Xdc_nonseq, data');

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
using QznbTestPackage
using Test

@testset "QznbTestPackage.jl" begin
    # Write your tests here.

    @test hello("qznb") == "Hello, qznb"

end

ext: jl

max_stars_repo_name: qz-michael/QznbTestPackage

max_issues_repo_name: qz-michael/QznbTestPackage

max_forks_repo_name: qz-michael/QznbTestPackage


text:
#=

JuliaBolt is ported from Neobolt: Neo4j Bolt connector for Python
https://github.com/neo4j-drivers/neobolt

Original notices:

# Copyright (c) 2002-2018 "Neo4j,"
# Neo4j Sweden AB [http://neo4j.com]
#
# This file is part of Neo4j.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

=#

#####################################
#              Errors               #
#####################################

abstract type JuliaBoltError <: Exception end

# Errors -------------------------- #

abstract type AbstractSecurityError <: JuliaBoltError end

# Errors -------------------------- #

abstract type AbstractCypherError <: JuliaBoltError end

# Errors -------------------------- #

abstract type AbstractTransientError <: AbstractCypherError end

# Errors -------------------------- #

abstract type AbstractClientError <: AbstractCypherError end

# Errors -------------------------- #

struct ProtocolError <: JuliaBoltError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct ServiceUnavailable <: JuliaBoltError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct ConnectionExpired <: JuliaBoltError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct SecurityError <: AbstractSecurityError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct CypherError <: AbstractCypherError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct ClientError <: AbstractClientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct DatabaseError <: AbstractCypherError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct TransientError <: AbstractTransientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct DatabaseUnavailableError <: AbstractTransientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct ConstraintError <: AbstractClientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct CypherSyntaxError <: AbstractClientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct CypherTypeError <: AbstractClientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct NotALeaderError <: AbstractClientError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct Forbidden <: AbstractSecurityError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct ForbiddenOnReadOnlyDatabaseError <: AbstractSecurityError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

struct AuthError <: AbstractSecurityError
    message
    code
    classification
    category
    title
    metadata
end

# Errors -------------------------- #

const client_errors = Dict(
    # ConstraintError
    "Neo.ClientError.Schema.ConstraintValidationFailed" => ConstraintError,
    "Neo.ClientError.Schema.ConstraintViolation" => ConstraintError,
    "Neo.ClientError.Statement.ConstraintVerificationFailed" => ConstraintError,
    "Neo.ClientError.Statement.ConstraintViolation" => ConstraintError,

    # CypherSyntaxError
    "Neo.ClientError.Statement.InvalidSyntax" => CypherSyntaxError,
    "Neo.ClientError.Statement.SyntaxError" => CypherSyntaxError,

    # CypherTypeError
    "Neo.ClientError.Procedure.TypeError" => CypherTypeError,
    "Neo.ClientError.Statement.InvalidType" => CypherTypeError,
    "Neo.ClientError.Statement.TypeError" => CypherTypeError,

    # Forbidden
    "Neo.ClientError.General.ForbiddenOnReadOnlyDatabase" => ForbiddenOnReadOnlyDatabaseError,
    "Neo.ClientError.General.ReadOnly" => Forbidden,
    "Neo.ClientError.Schema.ForbiddenOnConstraintIndex" => Forbidden,
    "Neo.ClientError.Schema.IndexBelongsToConstraint" => Forbidden,
    "Neo.ClientError.Security.Forbidden" => Forbidden,
    "Neo.ClientError.Transaction.ForbiddenDueToTransactionType" => Forbidden,

    # AuthError
    "Neo.ClientError.Security.AuthorizationFailed" => AuthError,
    "Neo.ClientError.Security.Unauthorized" => AuthError,

    # NotALeaderError
    "Neo.ClientError.Cluster.NotALeader" => NotALeaderError
)

# Errors -------------------------- #

const transient_errors = Dict(
    # DatabaseUnavailableError
    "Neo.TransientError.General.DatabaseUnavailable" => DatabaseUnavailableError
)

# Errors -------------------------- #

function hydrate_cypher_error(metadata::Dict)
    if haskey(metadata, "message")
        message = metadata["message"]
    else
        message = "An unknown error occurred."
    end

    if haskey(metadata, "code")
        code = metadata["code"]
    else
        code = "Neo.DatabaseError.General.UnknownError"
    end

    classification, category, title = "DatabaseError", "General", "UnknownError"
    try _, classification, category, title = split(code, ".") catch end

    error_type = extract_error_type(classification, code)

    return error_type(message, code, classification, category, title, metadata)
end

# Errors -------------------------- #

function extract_error_type(classification::AbstractString, code::AbstractString)
    if classification == "ClientError"
        if haskey(client_errors, code)
            client_errors[code]
        else
            return ClientError
        end
    elseif classification == "TransientError"
        if haskey(transient_errors, code)
            transient_errors[code]
        else
            return TransientError
        end
    elseif classification == "DatabaseError"
        return DatabaseError
    else
        return CypherError
    end
end

# Errors -------------------------- #
# # # # # # # # # # # # # # # # # # #

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
using DataDeps

register(DataDep(
    "Twitter Sentiment Dataset",
    """
    Twitter Sentiment Corpus by Niek Sanders
    This is downloaded from:
    http://help.sentiment140.com/for-students

    This dataset contains Twitter tweets:
    1600000 training examples
    498 test examples

    If you are using this data, cite Sentiment140
    https://cs.stanford.edu/people/alecmgo/papers/TwitterDistantSupervision09.pdf
    """,
    "https://cs.stanford.edu/people/alecmgo/trainingandtestdata.zip",
    "004a3772c8a7ff9bbfeb875880f47f0679d93fc63e5cf9cff72d54a8a6162e57",
    post_fetch_method=unpack
))


text:
struct QuilAddress{VT<:Vector}
    address_map::VT
end

Base.copy(info::QuilAddress) = QuilAddress(info.address_map)

function Base.show(io::IO, ::MIME"quil", blk::AbstractBlock)
    qcode = quil(blk)
    println(io, qcode)
end
quil(blk::AbstractBlock, N::Union{Nothing, Int}=nothing) = quil(blk, QuilAddress(collect(1:(N isa Nothing ? nqubits(blk) : N))))
quil(blk::CachedBlock, info::QuilAddress) = quil(content(blk), info)

######################## Basic Buiding Blocks ########################
const BASIC_GATES = [:X, :Y, :Z, :H, :S, :Sdag, :T, :Tdag, :SWAP, :CNOT]

locs2str(locs) = join(["$(i-1)" for i in locs], " ")
floats2str(params) = join(["$p" for p in params], ", ")
args2str(pcounts) = join(["p$i" for i=pcounts], ", ")
# place an operator on some location
import Base: >, /
>(gate::String, locs) = gate * " " * locs2str(locs)
/(locs, info::QuilAddress) = [info.address_map[loc] for loc in locs]

function quil(blk::Measure{N,K,ComputationalBasis}, info::QuilAddress) where {N, K}
    if blk.postprocess isa ResetTo
        return "MEASURE_AND_RESET($(blk.postprocess.x))" > info.address_map[blk.locations]
    elseif blk.postprocess isa NoPostProcess
        if N == K
            locs = 1:N
        else
            locs = blk.locations
        end
        return "MEASURE" > info.address_map[locs]
    else
        error("unable to parse $blk")
    end
end
quil(blk::PutBlock, info::QuilAddress) = quil(content(blk), info) > blk.locs/info
function quil(blk::Union{KronBlock, PauliString}, info::QuilAddress)
    locs = blk isa PauliString ? (1:nqubits(blk)) : blk.locs
    qs = String[]
    for (g, loc) in zip(subblocks(blk), locs)
        if !(g isa I2Gate)
            push!(qs, quil(g, info) > loc/info)
        end
    end
    join(qs, "\n")
end
quil(blk::RepeatedBlock, info::QuilAddress) = join([quil(content(blk), info) > loc/info for loc in blk.locs], "\n")

for G in BASIC_GATES
    GT = Symbol(G, :Gate)
    g = G |> String
    @eval quil(blk::$GT, info::QuilAddress) = "$($g)"
end
quil(blk::I2Gate, info::QuilAddress) = "I"

quil(blk::RotationGate{N,T}, info::QuilAddress) where {N,T} = "R$(quil(content(blk), info))($(blk.theta))"
quil(blk::RotationGate{N,T,G}, info::QuilAddress) where {N,T,G<:SWAPGate} = "P$(quil(content(blk), info))($(blk.theta))"  # @
quil(blk::ShiftGate, info::QuilAddress) = "PHASE($(blk.theta))"
_cstring(blk::ControlBlock) = prod(v==1 ? "C" : "¬C" for v in blk.ctrl_config)
function quil(blk::ControlBlock, info::QuilAddress)
    _cstring(blk) * quil(content(blk), info) > (blk.ctrl_locs..., blk.locs...)/info
end
quil(blk::ControlBlock{<:Any, <:XGate}, info::QuilAddress) = _cstring(blk) * "NOT" > (blk.ctrl_locs..., blk.locs...)/info  # @
function quil(blk::Subroutine, info::QuilAddress)
    info = copy(info)
    subinfo = QuilAddress(info.address_map[[blk.locs...]])
    quil(content(blk), subinfo)
end

function quil(blk::ChainBlock, info::QuilAddress)
    join([quil(b, info) for b in subblocks(blk)], "\n")
end

# TODO:
#
# QUIL: https://github.com/rigetti/pyquil/blob/master/pyquil/gates.py
# Quantum ->
# CPHASE00 -> diag([exp(1j*phi), 1, 1, 1])
# CPHASE01 -> diag([1, exp(1j*phi), 1, 1])
# CPHASE10 -> diag([1, 1, exp(1j*phi), 1])
# ISWAP -> [1000; 00i0; 0i00; 0001]
#
# CLASSICAL ->
# WAIT, NOP, HALT, RESET
# TRUE, FALSE, NEG, NOT, AND, OR, IOR, XOR, MOVE, EXCHANGE, LOAD, STORE, CONVERT, ADD, SUB, MUL, DIV, EQ, LT, LE, GT, GE
#
# if statement

# NOTE: give up support to roller


text:
#using Crystalline, KdotP

sgnum, klab, index = 213, "B", 1 #4, "G", 1#80, "P", 2#218, "R", 3
lgirsd = lgirreps(sgnum, Val(3))
lgirs = realify(lgirsd[klab])
lgir = lgirs[end]
H = kdotp(lgir)

## Visualize dispersion with a kz-slider
using GLMakie

f = Figure()
ax = f[1,1] = Axis3(f, aspect=(.5,.5,1), perspectiveness=1)
sl = f[2,1] = Slider(f, startvalue=0, range=-0.5:.025:.5) # k_z slider
N = 100

kxys = range(-.5,.5, N)
coefs = rand(length(H.cs))
tmp = lift(sl.value) do kz
    [eigvals(H((kx,ky,kz), coefs)) for kx in kxys, ky in kxys]
end
Es = [ lift(val -> getindex.(val,i),tmp) for i in 1:irdim(H)]
surface!.(ax, Ref(kxys), Ref(kxys), Es, colormap=:balance, colorrange=(-1,1).*.75, alpha=.5, transparency=true)
contourplots = contour3d!.(ax, Ref(kxys), Ref(kxys), Es, levels=range(-.75,.75,25), linewidth=1.5, transparency=true, colormap=:balance)

for elem in contourplots # fix until https://github.com/JuliaPlots/Makie.jl/pull/1727
    elem.plots[1].transparency=true
end

f

## Density of states near degeneracy
function dos(H, es, ks, γ = 0.01)
    d = zeros(Float64, length(es))
    for k in Iterators.product(ks,ks,ks)
        Eₖ = eigvals!(H(k))
        for Eₙₖ in Eₖ
            ωₙₖ = ComplexF64(Eₙₖ, γ)
            d .+= imag.(1.0 ./ (es .- ωₙₖ))
        end
    end
    return d ./ length(ks)^3
end

es = range(-0.25,0.25,50)
plot(es, dos(H, es, kxys))


text:
"""
Compute velocity `u` and angular velocity `ω` from position.

Assumes θ in radians.
Assumes values at equally spaced time intervals.
"""
function kinematics_from_position(
    x::Vector{Float64},
    y::Vector{Float64},
    θ::Vector{Float64};
    fps::Int=60,
    smooth::Bool=false,
    smooth_wnd::Float64=0.1,  # smoothing window size in seconds
)
    # compute speed
    speed = sqrt.(Δ(x) .^ 2 .+ Δ(y) .^ 2) .* fps

    # compute longitudinal and lateral velocities
    β = atan.(Δ(y), (Δ(x) .+ eps())) .- θ
    u = @. speed * cos(β)
    v = @. speed * sin(β)

    # compute angular velocity
    time = (collect(0:length(speed)) ./ fps)[2:end]
    ω = ∂(time, unwrap(θ))

    # smooth
    if smooth
        wnd = int(fps * smooth_wnd)
        speed = movingaverage(speed, wnd)
        u = movingaverage(u, wnd)
        v = movingaverage(v, wnd)
        ω = movingaverage(ω, wnd)
    end

    return speed, u, v, ω
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
using StochasticDiffEq, OrdinaryDiffEq, Zygote
using DiffEqSensitivity, Test, ForwardDiff

abstol=1e-12
reltol=1e-12
savingtimes=0.5

function test_discrete_callback(cb, tstops, g, dg!, cboop=nothing)
  function fiip(du,u,p,t)
    du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
    du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  end
  function foop(u,p,t)
    dx = p[1]*u[1] - p[2]*u[1]*u[2]
    dy = -p[3]*u[2] + p[4]*u[1]*u[2]
    [dx,dy]
  end

  p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]

  prob = ODEProblem(fiip,u0,(0.0,10.0),p)
  proboop = ODEProblem(foop,u0,(0.0,10.0),p)

  sol1 = solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes)
  sol2 = solve(prob,Tsit5(),u0=u0,p=p,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes)

  if cb.save_positions == [1,1]
    @test length(sol1.t) != length(sol2.t)
  else
    @test length(sol1.t) == length(sol2.t)
  end

  du01,dp1 = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint())),
    u0,p)

  du01b,dp1b = Zygote.gradient(
    (u0,p)->g(solve(proboop,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint())),
    u0,p)

  du01c,dp1c = Zygote.gradient(
    (u0,p)->g(solve(proboop,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint(checkpointing=false))),
    u0,p)

  if cboop === nothing
    du02,dp2 = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=ReverseDiffAdjoint()))
    ,u0,p)
  else
    du02,dp2 = Zygote.gradient(
      (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cboop,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=ReverseDiffAdjoint()))
      ,u0,p)
  end

  du03,dp3 = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=InterpolatingAdjoint(checkpointing=true))),
    u0,p)

  du03c,dp3c = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=InterpolatingAdjoint(checkpointing=false))),
    u0,p)

  du04,dp4 = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=QuadratureAdjoint())),
    u0,p)

  dstuff = ForwardDiff.gradient(
    (θ)->g(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:6],callback=cb,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes)),
    [u0;p])

  @info dstuff

  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:6]
  @test du01b ≈ dstuff[1:2]
  @test dp1b ≈ dstuff[3:6]
  @test du01c ≈ dstuff[1:2]
  @test dp1c ≈ dstuff[3:6]
  @test du01 ≈ du02
  @test du01 ≈ du03 rtol=1e-7
  @test du01 ≈ du03c rtol=1e-7
  @test du03 ≈ du03c
  @test du01 ≈ du04
  @test dp1 ≈ dp2
  @test dp1 ≈ dp3
  @test dp1 ≈ dp3c
  @test dp1 ≈ dp4 rtol=1e-7

  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:6]

  cb2 = DiffEqSensitivity.track_callbacks(CallbackSet(cb),prob.tspan[1],prob.u0,prob.p)
  sol_track = solve(prob,Tsit5(),u0=u0,p=p,callback=cb2,tstops=tstops,abstol=abstol,reltol=reltol,saveat=savingtimes)
  #cb_adj = DiffEqSensitivity.setup_reverse_callbacks(cb2,BacksolveAdjoint())

  adj_prob = ODEAdjointProblem(sol_track,BacksolveAdjoint(),dg!,sol_track.t,nothing,
  						 callback = cb2,
  						 abstol=abstol,reltol=reltol)
  adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  @test du01 ≈ -adj_sol[1:2,end]
  @test dp1 ≈ adj_sol[3:6,end]

  # adj_prob = ODEAdjointProblem(sol_track,InterpolatingAdjoint(),dg!,sol_track.t,nothing,
  # 						 callback = cb2,
  # 						 abstol=abstol,reltol=reltol)
  # adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  #
  # @test du01 ≈ -adj_sol[1:2,end]
  # @test dp1 ≈ adj_sol[3:6,end]
end

function test_continuous_wrt_discrete_callback()
  # test the continuous callbacks wrt to the equivalent discrete callback
  function f(du,u,p,t)
    #Bouncing Ball
    du[1] = u[2]
    du[2] = -p[1]
  end

  # no saving in Callbacks; prescribed vafter and vbefore; loss on the endpoint

  tstop = 3.1943828249997
  vbefore = -31.30495168499705
  vafter = 25.04396134799764

  u0 = [50.0,0.0]
  tspan = (0.0,5.0)
  p = [9.8, 0.8]

  prob = ODEProblem(f,u0,tspan,p)

  function condition(u,t,integrator) # Event when event_f(u,t) == 0
     t - tstop
  end
  function affect!(integrator)
    integrator.u[2] += vafter-vbefore
  end
  cb = ContinuousCallback(condition,affect!,save_positions=(false,false))

  condition2(u,t,integrator) = t == tstop
  cb2 = DiscreteCallback(condition2,affect!,save_positions=(false,false))

  du01,dp1 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb2,tstops=[tstop],
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  du02,dp2 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb,
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  dstuff = ForwardDiff.gradient((θ)-> sum(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],
     	callback=cb,saveat=tspan[2], save_start=false)),[u0;p])

  @info dstuff
  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:4]

  # no saving in Callbacks; prescribed vafter and vbefore; loss on the endpoint by slicing
  du01,dp1 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb2,tstops=[tstop],
     sensealg=BacksolveAdjoint())[end]),u0,p)

  du02,dp2 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb,
     sensealg=BacksolveAdjoint())[end]),u0,p)

  dstuff = ForwardDiff.gradient((θ)-> sum(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],
     	callback=cb)[end]),[u0;p])

  @info dstuff
  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:4]

  # with saving in Callbacks; prescribed vafter and vbefore; loss on the endpoint
  cb = ContinuousCallback(condition,affect!,save_positions=(true,true))
  cb2 = DiscreteCallback(condition2,affect!,save_positions=(true,true))

  du01,dp1 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb2,tstops=[tstop],
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  du02,dp2 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb,
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  dstuff = ForwardDiff.gradient((θ)-> sum(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],
     	callback=cb,saveat=tspan[2], save_start=false)),[u0;p])

  @info dstuff
  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:4]

  # with saving in Callbacks; prescribed vafter and vbefore; loss on the endpoint by slicing
  du01,dp1 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb2,tstops=[tstop],
     sensealg=BacksolveAdjoint())[end]),u0,p)

  du02,dp2 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb,
     sensealg=BacksolveAdjoint())[end]),u0,p)

  dstuff = ForwardDiff.gradient((θ)-> sum(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],
     	callback=cb)[end]),[u0;p])

  @info dstuff
  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:4]

  # with saving in Callbacks;  different affect function
  function affect2!(integrator)
    integrator.u[2] = -integrator.p[2]*integrator.u[2]
  end
  cb = ContinuousCallback(condition,affect2!,save_positions=(true,true))

  cb2 = DiscreteCallback(condition2,affect2!,save_positions=(true,true))

  du01,dp1 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb2,tstops=[tstop],
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  du02,dp2 = Zygote.gradient(
    (u0,p)->sum(solve(prob,Tsit5(),u0=u0,p=p,
     callback=cb,
     sensealg=BacksolveAdjoint(),
     saveat=tspan[2], save_start=false)),u0,p)

  dstuff = ForwardDiff.gradient((θ)-> sum(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],
     	callback=cb,saveat=tspan[2], save_start=false)),[u0;p])

  @info dstuff
  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du02 ≈ dstuff[1:2]
  @test dp2 ≈ dstuff[3:4]
  @test du01 ≈ du02
  @test dp1 ≈ dp2
end

function test_continuous_callback(cb, g, dg!)
  function fiip(du,u,p,t)
    du[1] = u[2]
    du[2] = -p[1]
  end
  function foop(u,p,t)
    dx = u[2]
    dy = -p[1]
    [dx,dy]
  end

  u0 = [50.0,0.0]
  tspan = (0.0,5.0)
  p = [9.8, 0.8]

  prob = ODEProblem(fiip,u0,tspan,p)
  proboop = ODEProblem(fiip,u0,tspan,p)

  sol1 = solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes)
  sol2 = solve(prob,Tsit5(),u0=u0,p=p,abstol=abstol,reltol=reltol,saveat=savingtimes)

  if cb.save_positions == [1,1]
    @test length(sol1.t) != length(sol2.t)
  else
    @test length(sol1.t) == length(sol2.t)
  end

  du01,dp1 = @time Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint())),
    u0,p)

  du01b,dp1b = Zygote.gradient(
    (u0,p)->g(solve(proboop,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint())),
    u0,p)

  du01c,dp1c = Zygote.gradient(
    (u0,p)->g(solve(proboop,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=BacksolveAdjoint(checkpointing=false))),
    u0,p)

  @test_broken du02,dp2 = @time Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=ReverseDiffAdjoint()))
    ,u0,p)

  du03,dp3 = @time Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=InterpolatingAdjoint(checkpointing=true))),
    u0,p)

  du03c,dp3c = Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=InterpolatingAdjoint(checkpointing=false))),
    u0,p)

  du04,dp4 = @time Zygote.gradient(
    (u0,p)->g(solve(prob,Tsit5(),u0=u0,p=p,callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes,sensealg=QuadratureAdjoint())),
    u0,p)

  dstuff = @time ForwardDiff.gradient(
    (θ)->g(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:4],callback=cb,abstol=abstol,reltol=reltol,saveat=savingtimes)),
    [u0;p])

  @info dstuff

  @test du01 ≈ dstuff[1:2]
  @test dp1 ≈ dstuff[3:4]
  @test du01b ≈ dstuff[1:2]
  @test dp1b ≈ dstuff[3:4]
  @test du01c ≈ dstuff[1:2]
  @test dp1c ≈ dstuff[3:4]
  @test_broken du01 ≈ du02
  @test du01 ≈ du03 rtol=1e-7
  @test du01 ≈ du03c rtol=1e-7
  @test du03 ≈ du03c
  @test du01 ≈ du04
  @test_broken dp1 ≈ dp2
  @test dp1 ≈ dp3
  @test dp1 ≈ dp3c
  @test dp1 ≈ dp4 rtol=1e-7

  @test_broken du02 ≈ dstuff[1:2]
  @test_broken dp2 ≈ dstuff[3:4]

  cb2 = DiffEqSensitivity.track_callbacks(CallbackSet(cb),prob.tspan[1],prob.u0,prob.p)
  sol_track = solve(prob,Tsit5(),u0=u0,p=p,callback=cb2,abstol=abstol,reltol=reltol,saveat=savingtimes)

  adj_prob = ODEAdjointProblem(sol_track,BacksolveAdjoint(),dg!,sol_track.t,nothing,
  						 callback = cb2,
  						 abstol=abstol,reltol=reltol)
  adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  @test du01 ≈ -adj_sol[1:2,end]
  @test dp1 ≈ adj_sol[3:4,end]

  # adj_prob = ODEAdjointProblem(sol_track,InterpolatingAdjoint(),dg!,sol_track.t,nothing,
  # 						 callback = cb2,
  # 						 abstol=abstol,reltol=reltol)
  # adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  #
  # @test du01 ≈ -adj_sol[1:2,end]
  # @test dp1 ≈ adj_sol[3:6,end]
end

function test_SDE_callbacks()
  function dt!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α*x - β*x*y
    du[2] = dy = -δ*y + γ*x*y
  end

  function dW!(du, u, p, t)
    du[1] = 0.1u[1]
    du[2] = 0.1u[2]
  end

  u0 = [1.0,1.0]
  tspan = (0.0, 10.0)
  p = [2.2, 1.0, 2.0, 0.4]
  prob_sde = SDEProblem(dt!, dW!, u0, tspan,p)

  condition(u,t,integrator) = integrator.t>9.0 #some condition
  function affect!(integrator)
  	 #println("Callback")  #some callback
  end
  cb = DiscreteCallback(condition,affect!,save_positions=(false,false))

  function predict_sde(p)
    return Array(solve(prob_sde, EM(), p=p, saveat=savingtimes,sensealg = ForwardDiffSensitivity(), dt=0.001, callback=cb))
  end

  loss_sde(p)= sum(abs2, x-1 for x in predict_sde(p))

  loss_sde(p)
  @time dp = gradient(p) do p
  	loss_sde(p)
  end

  @test !iszero(dp[1])
end

@testset "Test callbacks" begin
  println("Discrete Callbacks")
  @testset "Discrete callbacks" begin
  	@testset "ODEs" begin
  	  println("ODEs")
      @testset "simple loss function" begin
        g(sol) = sum(sol)
        function dg!(out,u,p,t,i)
          (out.=-1)
        end
        @testset "callbacks with no effect" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 0.0
          cb = DiscreteCallback(condition,affect!,save_positions=(false,false))
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callbacks with no effect except saving the state" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 0.0
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 2.0
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callback at multiple time points" begin
          affecttimes = [2.03,4.0,8.0]
          condition(u,t,integrator) = t ∈ affecttimes
          affect!(integrator) = integrator.u[1] += 2.0
          cb = DiscreteCallback(condition,affect!)
          test_discrete_callback(cb,affecttimes,g,dg!)
        end
        @testset "state-dependent += callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = (integrator.u .+= integrator.p[2]/8*sin.(integrator.u))
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "other callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = (integrator.u[1] = 2.0; @show "triggered!")
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "parameter changing callback at single time point" begin
          condition(u,t,integrator) = t == 5.1
          affect!(integrator) = (integrator.p .= 2*integrator.p .- 0.5)
          affect(integrator) = (integrator.p = 2*integrator.p .- 0.5)
          cb = DiscreteCallback(condition,affect!)
          cboop = DiscreteCallback(condition,affect)
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.1]
          test_discrete_callback(cb,tstops,g,dg!,cboop)
        end
      end
      @testset "MSE loss function" begin
        g(u) = sum((1.0.-u).^2)./2
        dg!(out,u,p,t,i) = (out.=1.0.-u)
        @testset "callbacks with no effect" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 0.0
          cb = DiscreteCallback(condition,affect!,save_positions=(false,false))
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callbacks with no effect except saving the state" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 0.0
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = integrator.u[1] += 2.0
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "callback at multiple time points" begin
          affecttimes = [2.03,4.0,8.0]
          condition(u,t,integrator) = t ∈ affecttimes
          affect!(integrator) = integrator.u[1] += 2.0
          cb = DiscreteCallback(condition,affect!)
          test_discrete_callback(cb,affecttimes,g,dg!)
        end
        @testset "state-dependent += callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = (integrator.u .+= integrator.p[2]/8*sin.(integrator.u))
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "other callback at single time point" begin
          condition(u,t,integrator) = t == 5
          affect!(integrator) = (integrator.u[1] = 2.0; @show "triggered!")
          cb = DiscreteCallback(condition,affect!)
          tstops=[5.0]
          test_discrete_callback(cb,tstops,g,dg!)
        end
        @testset "parameter changing callback at single time point" begin
          condition(u,t,integrator) = t == 5.1
          affect!(integrator) = (integrator.p .= 2*integrator.p .- 0.5)
          affect(integrator) = (integrator.p = 2*integrator.p .- 0.5)
          cb = DiscreteCallback(condition,affect!)
          cboop = DiscreteCallback(condition,affect)
          tstops=[5.1]
          test_discrete_callback(cb,tstops,g,dg!,cboop)
        end
      end
  	end
  	@testset "SDEs" begin
  	  println("SDEs")
      test_SDE_callbacks()
  	end
  end

  println("Continuous Callbacks")
  @testset "Continuous callbacks" begin
    @testset "Compare with respect to discrete callback" begin
      test_continuous_wrt_discrete_callback()
    end
    @testset "simple loss function" begin
      g(sol) = sum(sol)
      function dg!(out,u,p,t,i)
        (out.=-1)
      end
      @testset "callbacks with no effect" begin
        condition(u,t,integrator) = u[1] # Event when event_f(u,t) == 0
        affect!(integrator) = (integrator.u[2] += 0)
        cb = ContinuousCallback(condition,affect!,save_positions=(false,false))
        test_continuous_callback(cb,g,dg!)
      end
    #   @testset "callbacks with no effect except saving the state" begin
    #     condition(u,t,integrator) = u[1]
    #     affect!(integrator) = (integrator.u[2] += 0)
    #     cb = ContinuousCallback(condition,affect!)
    #     test_continuous_callback(cb,g,dg!)
    #   end
    #   @testset "+= callback" begin
    #     condition(u,t,integrator) = u[1]
    #     affect!(integrator) = (integrator.u[2] += 50.0)
    #     cb = ContinuousCallback(condition,affect!)
    #     test_continuous_callback(cb,g,dg!)
    #   end
    #   @testset "= callback" begin
    #     condition(u,t,integrator) = u[1]
    #     affect!(integrator) = (integrator.u[2] = -integrator.p[2]*integrator.u[2])
    #     cb = ContinuousCallback(condition,affect!)
    #     test_continuous_callback(cb,g,dg!)
    #   end
    end
  end
end


text:
using AdversarialGridworld
using Test
using Random

# Generate a sample task
mdp, policy, adv_mdp = generate_task(N_rewards=1, N_walls=1, Nx=4, Ny=4,
                                     return_orig=true, n_generative_samples_train=1,
                                     max_iterations_train=10, verbose = false)

@test mdp.agent_gets_action == :ego
@test adv_mdp.agent_gets_action == :adversary


text:
matdescra(A::LowerTriangular) = "TLNF"
matdescra(A::UpperTriangular) = "TUNF"
matdescra(A::Diagonal) = "DUNF"
matdescra(A::UnitLowerTriangular) = "TLUF"
matdescra(A::UnitUpperTriangular) = "TUUF"
matdescra(A::Symmetric) = string('S', A.uplo, 'N', 'F')
matdescra(A::Hermitian) = string('H', A.uplo, 'N', 'F')
matdescra(A::SparseMatrixCSC) = "GUUF"


text:
using NearlyNewton
using Test
using StaticArrays
#using Optim
#using LineSearches
using Printf
using LinearAlgebra: norm, I

include("noalloc.jl")

@testset "NearlyNewton.jl" begin

    function theta(x)
       if x[1] > 0
           return atan(x[2] / x[1]) / (2.0 * pi)
       else
           return (pi + atan(x[2] / x[1])) / (2.0 * pi)
       end
    end
    f(x) = 100.0 * ((x[3] - 10.0 * theta(x))^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

    function f∇f!(∇f, x)
        if !(∇f==nothing)
            if ( x[1]^2 + x[2]^2 == 0 )
                dtdx1 = 0;
                dtdx2 = 0;
            else
                dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
                dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
            end
            ∇f[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
            ∇f[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
            ∇f[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
        end

        fx = f(x)
        return ∇f==nothing ? fx : (fx, ∇f)
    end

    function f∇f(∇f, x)
        if !(∇f == nothing)
            gx = similar(x)
            return f∇f!(gx, x)
        else
            return f∇f!(∇f, x)
        end
    end
    function f∇fs(∇f, x)
        if !(∇f == nothing)
            if ( x[1]^2 + x[2]^2 == 0 )
                dtdx1 = 0;
                dtdx2 = 0;
            else
                dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) )
                dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) )
            end

            s1 = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 )
            s2 = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
                200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 )
            s3 = 200.0*(x[3]-10.0*theta(x)) + 2.0*x[3]
            ∇f = @SVector [s1, s2, s3]
            return f(x), ∇f
        else
            return f(x)
        end
    end

    x0 = [-1.0, 0.0, 0.0]
    xopt = [1.0, 0.0, 0.0]
    x0s = @SVector [-1.0, 0.0, 0.0]

    println("Starting  from: ", x0)
    println("Targeting     : ", xopt)

    shortname(::GradientDescent) = "GD  "
    shortname(::InverseApprox) = "(inverse)"
    shortname(::NearlyNewton.DirectApprox) = " (direct)"

    function printed_minimize(f∇f, x0, method, B0=nothing)
        res = minimize(f∇f, x0, method, B0, NearlyNewton.OptOptions())
        print("NN  $(shortname(method)) $(shortname(method.approx)): ")
        @printf("%2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    end
    function printed_minimize!(f∇f!, x0, method, B0=nothing)
        res = minimize!(f∇f!, x0, method, B0, NearlyNewton.OptOptions())
        print("NN! $(shortname(method)) $(shortname(method.approx)): ")
        @printf("%2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    end

    printed_minimize(f∇f, x0, GradientDescent(InverseApprox()))
    printed_minimize!(f∇f!, copy(x0), GradientDescent(InverseApprox()))
    res = minimize(f∇fs, x0s, GradientDescent(InverseApprox()))
    @printf("NN  GD(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

    printed_minimize(f∇f, x0, GradientDescent(DirectApprox()))
    printed_minimize!(f∇f!, copy(x0), GradientDescent(DirectApprox()))
    res = minimize(f∇fs, x0s, GradientDescent(DirectApprox()))
    @printf("NN  GD(S)  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()

    res = minimize(f∇f, x0, BFGS(InverseApprox()), nothing, NearlyNewton.OptOptions())
    @test res[4] == 30
    @printf("NN  BFGS    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), (BFGS(InverseApprox()), BackTracking()), nothing, NearlyNewton.OptOptions())
    @test res[4] == 30
    @printf("NN! BFGS    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, BFGS(InverseApprox()))
    @test res[4] == 30
    @printf("NN  BFGS(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

    res = minimize(f∇f, x0, BFGS(DirectApprox()), nothing, NearlyNewton.OptOptions())
    @printf("NN  BFGS    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), BFGS(DirectApprox()), nothing, NearlyNewton.OptOptions())
    @printf("NN! BFGS    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, BFGS(DirectApprox()))
    @printf("NN  BFGS(S) (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()
    # res = optimize(f∇f!, copy(x0), Optim.BFGS(linesearch=BackTracking()))
    # @time optimize(f∇f!, copy(x0), Optim.BFGS(linesearch=BackTracking()))
    # @printf("OT! BFGS (inverse): %2.2e  %2.2e\n", norm(Optim.minimizer(res)-xopt,Inf), Optim.g_residual(res))

    res = minimize(f∇f, x0, DFP(InverseApprox()))
    @printf("NN  DFP    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), DFP(InverseApprox()))
    @printf("NN! DFP    (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, DFP(InverseApprox()))
    @printf("NN  DFP(S) (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, x0, DFP(DirectApprox()))
    @printf("NN  DFP    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), DFP(DirectApprox()))
    @printf("NN! DFP    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, DFP(DirectApprox()))
    @printf("NN  DFP(S) (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()

    res = minimize(f∇f, x0, SR1(InverseApprox()))
    @printf("NN  SR1     (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), SR1(InverseApprox()))
    @printf("NN! SR1     (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, SR1(InverseApprox()))
    @printf("NN  SR1(S)  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, x0, SR1(DirectApprox()))
    @printf("NN  SR1     (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(x0), SR1(DirectApprox()))
    @printf("NN! SR1     (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, x0s, SR1(DirectApprox()))
    @printf("NN  SR1(S)  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

    xrand = rand(3)
    xrands = SVector{3}(xrand)
    println("\nFrom a random point: ", xrand)
    res = minimize(f∇f, xrand, GradientDescent(InverseApprox()))
    @printf("NN  GD   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), GradientDescent(InverseApprox()))
    @printf("NN! GD   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, xrand, GradientDescent(DirectApprox()))
    @printf("NN  GD    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), GradientDescent(DirectApprox()))
    @printf("NN! GD    (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()
    res = minimize(f∇f, xrand, BFGS(InverseApprox()))
    @printf("NN  BFGS (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), (BFGS(InverseApprox()), BackTracking()))
    @printf("NN! BFGS (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, xrand, BFGS(DirectApprox()))
    @printf("NN  BFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), BFGS(DirectApprox()))
    @printf("NN! BFGS  (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()

    res = minimize(f∇f, xrand, DFP(InverseApprox()))
    @printf("NN  DFP  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), DFP(InverseApprox()))
    @printf("NN! DFP  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, xrand, DFP(DirectApprox()))
    @printf("NN  DFP   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), DFP(DirectApprox()))
    @printf("NN! DFP   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    println()

    res = minimize(f∇f, xrand, SR1(InverseApprox()))
    @printf("NN  SR1  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), SR1(InverseApprox()))
    @printf("NN! SR1  (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, xrand, SR1(InverseApprox()))
    @printf("NN  SR1(S)   (inverse): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇f, xrand, SR1(DirectApprox()))
    @printf("NN  SR1   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize!(f∇f!, copy(xrand), SR1(DirectApprox()))
    @printf("NN! SR1   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])
    res = minimize(f∇fs, xrand, SR1(DirectApprox()))
    @printf("NN  SR1(S)   (direct): %2.2e  %2.2e  %d\n", norm(res[1]-xopt,Inf), norm(res[3], Inf), res[4])

    println()

    function himmelblau!(∇f, x)
        if !(∇f == nothing)
            ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
                44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
            ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
        end

        fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        return ∇f == nothing ? fx : (fx, ∇f)
    end

    function himmelblaus(∇f, x)
        f = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        if !(∇f == nothing)
            ∇f1 = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
                44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
            ∇f2 = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
            return f, @SVector([∇f1, ∇f2])
        else
            return f
        end
    end

    function himmelblau(∇f, x)
        g = ∇f == nothing ? ∇f : similar(x)

        return himmelblau!(g, x)
    end

    println("\nHimmelblau function")
    x0 = [3.0, 1.0]
    x0s = @SVector([3.0, 1.0])
    res = minimize(himmelblau, x0, GradientDescent(InverseApprox()))
    @printf("NN  GD    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), GradientDescent(InverseApprox()))
    @printf("NN! GD    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblaus, x0s, GradientDescent(InverseApprox()))
    @printf("NN  GD(S) (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()
    res = minimize(himmelblau, x0, GradientDescent(DirectApprox()))
    @printf("NN  GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), GradientDescent(DirectApprox()))
    @printf("NN! GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblaus, x0s, GradientDescent(DirectApprox()))
    @printf("NN  GD(S) (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()

    res = minimize(himmelblau, x0, BFGS(InverseApprox()))
    @printf("NN  BFGS    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), BFGS(InverseApprox()))
    @printf("NN! BFGS    (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, BFGS(InverseApprox()))
    @printf("NN  BFGS(S) (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()
    res = minimize(himmelblau, x0, BFGS(DirectApprox()))
    @printf("NN  BFGS    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), BFGS(DirectApprox()))
    @printf("NN! BFGS    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, BFGS(DirectApprox()))
    @printf("NN  BFGS(S) (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()

    res = minimize(himmelblau, x0, DFP(InverseApprox()))
    @printf("NN  DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), DFP(InverseApprox()))
    @printf("NN! DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, DFP(InverseApprox()))
    @printf("NN  DFP(S)  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0, DFP(DirectApprox()))
    println()
    @printf("NN  DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), DFP(DirectApprox()))
    @printf("NN! DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, DFP(DirectApprox()))
    @printf("NN  DFP(S)   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()

    res = minimize(himmelblau, x0, SR1(InverseApprox()))
    @printf("NN  SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), SR1(InverseApprox()))
    @printf("NN! SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, SR1(InverseApprox()))
    @printf("NN  SR1(S)  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0, SR1(DirectApprox()))
    @printf("NN  SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(x0), SR1(DirectApprox()))
    @printf("NN! SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, x0s, SR1(DirectApprox()))
    @printf("NN  SR1(S)   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    println()
    xrand = rand(2)
    println("\nFrom a random point: ", xrand)

    res = minimize(himmelblau, xrand, GradientDescent(InverseApprox()))
    @printf("NN  GD   (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(xrand), GradientDescent(InverseApprox()))
    @printf("NN! GD   (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, xrand, GradientDescent(DirectApprox()))
    @printf("NN  GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(xrand), GradientDescent(DirectApprox()))
    @printf("NN! GD    (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

    res = minimize(himmelblau, xrand, BFGS(InverseApprox()))
    @printf("NN   BFGS (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(xrand), BFGS(InverseApprox()))
    @printf("NN!  BFGS (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, xrand, BFGS(DirectApprox()))
    @printf("NN(S) BFGS (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, xrand, BFGS(DirectApprox()))
    @printf("NN  BFGS  (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!, copy(xrand), BFGS(DirectApprox()))
    @printf("NN! BFGS  (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

    res = minimize(himmelblau, xrand, DFP(InverseApprox()))
    @printf("NN  DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!,  copy(xrand), DFP(InverseApprox()))
    @printf("NN! DFP  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, xrand, DFP(DirectApprox()))
    @printf("NN  DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!,  copy(xrand), DFP(DirectApprox()))
    @printf("NN! DFP   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

    res = minimize(himmelblau, xrand, SR1(InverseApprox()))
    @printf("NN  SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!,  copy(xrand), SR1(InverseApprox()))
    @printf("NN! SR1  (inverse): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize(himmelblau, xrand, SR1(DirectApprox()))
    @printf("NN  SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
    res = minimize!(himmelblau!,  copy(xrand), SR1(DirectApprox()))
    @printf("NN! SR1   (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])

end

include("geometrytypes.jl")
include("newton.jl")


text:
using Distributions
using Iterators
reload("src/nsum.jl")

#@show closest_grid_loc([1,2,3], 30)
#@show closest_grid_loc([1,15,25], 30)
#@show closest_grid_loc([60,59,61,32], 30)

function gen_grid(mins, maxs, N=30)
    D = length(mins)
    stepsizes = ceil(float(maxs .- mins) ./N)
    ranges = [mins[i]:stepsizes[i]:maxs[i] for i=1:D]
    grid = hcat(map(collect, Iterators.product(ranges...))...)'
    return map(length,ranges), stepsizes, grid
end

function gen_unit_grid(mins, maxs)
    D = length(mins)
    ranges = [mins[i]:maxs[i] for i=1:D]
    hcat(map(collect, Iterators.product(ranges...))...)'
end

function f(x, vals)
    global iters,evals
    iters += 1
    evals += length(vals)
    pdfs = pdf(Poisson(rate), floor(vec(x)))
    vals[:] = vec(prod(reshape(pdfs,size(x)...), 1))
end

function ftree(x)
    #println("Computing fvals for $(x[1])")
    global iters,evals
    iters += 1
    evals += length(x)
    for eval in x
        pt = eval.location
        vals = eval.vals
        assert(length(vals) == 0)
        pdfs = pdf(Poisson(rate), pt)
        push!(vals, prod(pdfs))
    end
end

function ftree_skinny(x)
    #println("Computing fvals for $(x[1])")
    global iters,evals
    iters += 1
    evals += length(x)
    for eval in x
        pt = eval.location
        vals = eval.vals
        assert(length(vals) == 0)
        val = 1.0
        val *= pdf(Poisson(0.01), pt[1])
        for v in pt[2:end]
            val *= pdf(Poisson(rate), v)
        end
        push!(vals, val)
    end
end

feasy(x) = (v = zeros(size(x)...); f(x,v); sum(v))

rate = 90000
iters = 0
evals = 0

upp = 130000
D = 2

lens = [2,16]
r = NSum.RegionTree(lens)
NSum.cut!(r)

println("####################")
ndivs=1000
xs=linspace(0,upp,ndivs)
spacing=upp/ndivs
est = spacing*sum(map(x->pdf(Poisson(rate), ifloor(x)), xs))
println("Estimate using $ndivs ndivs is $est")

println("####################")
@time tot1,r = NSum.nsum(1, ftree, zeros(D).+upp, abstol=0.03, maxevals=50)
@show tot1
println("With $iters iters and $evals fun evals")

println("####################")
iters = 0
evals = 0
@time tot2,r = NSum.nsum(1, ftree_skinny, zeros(D).+upp, abstol=0.03, maxevals=50)
@show tot2
println("With $iters iters and $evals fun evals")

#println("####################")
#iters = 0
#evals = 0
#@time tot,r = IntSum.nsum(ftree, zeros(D).+upp)
#@show tot
#println("With $iters iters and $evals fun evals")

#allnodes = collect(r)

#using PyPlot
#close("all")
#figure()
#for n in allnodes
    #if length(n.vals)!=0
        #plot(n.mins[1], n.mins[2], "g.")
    #end
#end
#lens, steps, grid = gen_grid(zeros(D), zeros(D).+upp.-1, 50)

#vals = zeros(size(grid,1))
#f(grid',vals)

#imshow(reshape(vals,lens...), extent=[0,upp,0,upp], origin="lower")

## NOTE TO SELF: trying to get above code to work to verify that my new
# find uneven branch code with midpoint is working

#using Cubature
#println("####################")
#@time tot1 = hcubature_v(f, zeros(D), zeros(D).+upp, abstol=0.01)
#iters = 0
#evals = 0
#@time tot2 = hcubature_v(f, zeros(D), zeros(D).+upp, abstol=0.01)
#@show tot1, tot2
#println("With $iters iters and $evals fun evals")

#iters = 0
#evals = 0
#lens, steps, grid = gen_grid(zeros(D).+100, zeros(D).+140, 20)
#vals = Array(Float64, size(grid,1))
#f(grid',vals)
#@show steps
#@show sum(vals)


text:
using COSMO, SparseArrays, LinearAlgebra, FileIO, Test

if !@isdefined test_problem_path
  test_problem_path = "./testproblems/";
end

data = load(test_problem_path * "decomp_testProblem.jld2")

P = data["P"];
q = data["q"];
A = data["A"];
b = data["b"];
obj_true = data["objTrue"];
cs = COSMO.Constraint(-A, b, COSMO.PsdCone)

model = COSMO.Model()
settings = COSMO.Settings(decompose = false, obj_true = obj_true)
assemble!(model, P, q, [cs], settings)
# # solve problem with chordal decomposition turned off
res = COSMO.optimize!(model);

P = data["P"];
q = data["q"];
A = data["A"];
b = data["b"];
cs = COSMO.Constraint(-A, b, COSMO.PsdCone)
# # solve problem with sparsity exploitation
model = COSMO.Model()
settings_decomp = COSMO.Settings(decompose = true, obj_true = obj_true, verbose_timing = true)
assemble!(model, P, q, [cs], settings_decomp);
res_decomp = COSMO.optimize!(model);

#  compare results
@testset "Chordal Decomposition" begin
  @test abs(res_decomp.obj_val - res.obj_val) < 1e-3
  @test abs(res_decomp.obj_val - obj_true) < 1e-3
  @test norm(res_decomp.x - res.x,Inf) < 1e-3
end
nothing


text:
#---------------------------------------------------------
# # [A coupled PO₄–POP model](@id P-model)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/3_Pmodel.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/3_Pmodel.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# In this tutorial we will explicitly simulate 2 tracers whose distributions control and feed back on each other.

# We consider a simple model for the cycling of phosphorus with 2 state variables consisting of phosphate (PO₄) AKA dissolved inorganic phosphorus (DIP) and particulate organic phosphorus (POP).
# The dissolved phases are transported by advection and diffusion whereas the particulate phase sinks rapidly down the water column without any appreciable transport by the circulation.
#
# The governing equations that couple the 3D concentration fields of DIP and POP, denoted $x_\mathsf{DIP}$ and $x_\mathsf{POP}$, respectively, are:
#
# $$\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \mathbf{K}\nabla )\right] x_\mathsf{DIP} = -U(x_\mathsf{DIP}) + R(x_\mathsf{POP}),$$
#
# and
#
# $$\left[\frac{\partial}{\partial t} + \nabla \cdot \boldsymbol{w}\right] x_\mathsf{POP} = U(x_\mathsf{DIP}) - R(x_\mathsf{POP}).$$

# The $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ and $\nabla \cdot \boldsymbol{w}$ operators represent the ocean circulation and the sinking of particles, respectively.
# ([Tracer transport operators are described in the documentation](@ref tracer-transport-operators).)

# The function $U$ represents the biological uptake of DIP by phytoplankton, which we model here as
#
# $$U(x_\mathsf{DIP}) = \frac{x_\mathsf{DIP}}{\tau_\mathsf{DIP}} \, \frac{x_\mathsf{DIP}}{x_\mathsf{DIP} + k} \, (z < z_0),$$
#
# with the timescale, $\tau$, the half-saturation rate $k$, and the depth $z_0$ as parameters.

# The function $R$ defines the remineralization rate of POP, which converts POP back into DIP.
# For the remineralization, we simply use a linear rate constant, i.e.,
#
# $$R(x_\mathsf{POP}) = \frac{x_\mathsf{POP}}{\tau_\mathsf{POP}}.$$

# We start by telling Julia we want to use the AIBECS and the OCIM0.1 circulation for DIP.

using AIBECS
grd, T_OCIM = OCIM0.load()
T_DIP(p) = T_OCIM

# For the sinking of particles, we use the `transportoperator` function

T_POP(p) = transportoperator(grd, z -> w(z,p))

# for which we need to define the sinking speed `w(z,p)` as a function of depth `z` and of the parameters `p`.
# Following the assumption that $w(z) = w_0 + w' z$ increases linearly with depth, we write it as

function w(z,p)
    @unpack w₀, w′ = p
    return @. w₀ + w′ * z
end

# ##### Uptake (DIP → POP)

# For the uptake, $U$, we write

z = depthvec(grd)
function U(x,p)
    @unpack τ_DIP, k, z₀ = p
    return @. x/τ_DIP * x/(x+k) * (z≤z₀) * (x≥0)
end

# where we have "unpacked" the parameters to make the code clearer and as close to the mathematical equation as possible.
# (Note we have also added a constraint that `x` must be positive for uptake to happen.)

# ##### Remineralization (POP → DIP)

# For the remineralization, $R$, we write

function R(x,p)
    @unpack τ_POP = p
    return x / τ_POP
end

# ##### Net sources and sinks

# We lump the sources and sinks into `G` functions for DIP and POP.

function G_DIP(DIP, POP, p)
    @unpack DIP_geo, τ_geo = p
    return @. -$U(DIP,p) + $R(POP,p) + (DIP_geo - DIP) / τ_geo
end
function G_POP(DIP, POP, p)
    @unpack τ_geo = p
    return @. $U(DIP,p) - $R(POP,p) - POP / τ_geo
end

# where we have imposed a slow restoring of DIP to the global mean `DIP_geo` to prescribe the global mean concentration.
# (The `$` signs in front of `U` and `R` protect them from the broadcast macro `@.`)

# We now define and build the parameters.

# In this tutorial we will specify some initial values for the parameters
# and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
using Unitful: m, d, s, yr, Myr, mol, mmol, μmol, μM
@initial_value @units struct PmodelParameters{U} <: AbstractParameters{U}
    w₀::U       |  0.64 | m/d
    w′::U       |  0.13 | m/d/m
    τ_DIP::U    | 230.0 | d
    k::U        |  6.62 | μmol/m^3
    z₀::U       |  80.0 | m
    τ_POP::U    |   5.0 | d
    τ_geo::U    |   1.0 | Myr
    DIP_geo::U  |  2.12 | mmol/m^3
end

# Finally, thanks to the initial values we provided, we can instantiate the parameter vector succintly as

p = PmodelParameters()

# We generate the state function `F`,

nb = sum(iswet(grd))
F = AIBECSFunction((T_DIP, T_POP), (G_DIP, G_POP), nb)

# generate the steady-state problem,

@unpack DIP_geo = p
x = DIP_geo * ones(2nb) # initial guess
prob = SteadyStateProblem(F, x, p)

# and solve it

sol = solve(prob, CTKAlg()).u

# We can look at different the DIP and POP fields using the Plots.jl recipes.

DIP, POP = state_to_tracers(sol, grd) # unpack tracers

# First, let's look at the mean profile

using Plots
plothorizontalmean(DIP * (mol/m^3) .|> μM, grd)

# We can plot the concentration of DIP at a given depth via, e.g.,

plothorizontalslice(DIP * (mol/m^3) .|> μM, grd, depth=1000m, color=:viridis)

# Or have a look at a map of the uptake at the surface

plotverticalintegral(U(DIP,p) * (mol/m^3/s) .|> mmol/yr/m^3, grd, color=:algae)

# Or look at what is exported below 500 m

plothorizontalslice(POP .* w(z,p) * (mol/m^3*m/s) .|> mmol/yr/m^2, grd, depth=500m, color=:inferno, rev=true)

# Now let's make our model a little fancier and use a fine topographic map to refine the remineralization profile.
# For this, we will use the ETOPO dataset, which can be downloaded by AIBECS via

f_topo = ETOPO.fractiontopo(grd)

# `f_topo` is the fraction of sediments located in each wet box of the `grd` grid.
# We can use it to redefine the transport operator for sinking particles to take into consideration the subgrid topography,
# such that the fine-resolution sediments intercept settling POP.

T_POP2(p) = transportoperator(grd, z -> w(z,p); frac_seafloor=f_topo)

# With this new vertical transport for POP, we can recreate our problem, solve it again

F2 = AIBECSFunction((T_DIP, T_POP2), (G_DIP, G_POP), nb)
prob2 = SteadyStateProblem(F2, x, p)
sol2 = solve(prob2, CTKAlg()).u
DIP2, POP2 = state_to_tracers(sol2, grd) # unpack tracers

# and check the difference

plotzonalaverage((DIP2 - DIP) ./ DIP .|> u"percent", grd, color=:balance, clim=(-5,5))

# This zonal average shows how much DIP is redistributed as it is prevented from sinking out of the surface layers with the new subgrid parameterization.

# Let's look at the vertical average.

plotverticalaverage((DIP2 - DIP) ./ DIP .|> u"percent", grd, color=:balance, clim=(-10,10))

# This shows minor changes on the order of 1%, on the global scale,
# except along the coasts, which retain much more DIP with the subgrid topography.


text:
"""
Fetch the "natural" block size of this band.

GDAL contains a concept of the natural block size of rasters so that
applications can organized data access efficiently for some file formats.
The natural block size is the block size that is most efficient for accessing
the format. For many formats this is simple a whole scanline in which case
`*pnXSize` is set to `GetXSize()`, and *pnYSize is set to 1.

However, for tiled images this will typically be the tile size.

Note that the X and Y block sizes don't have to divide the image size evenly,
meaning that right and bottom edge blocks may be incomplete. See `ReadBlock()`
for an example of code dealing with these issues.

###Parameters
* `pnXSize`     integer to put the X block size into or NULL.
* `pnYSize`     integer to put the Y block size into or NULL.
"""
_getblocksize(hBand::GDALRasterBandH, pnXSize::Ptr{Cint}, pnYSize::Ptr{Cint}) =
    GDALGetBlockSize(hBand, pnXSize, pnYSize)

function getblocksize(rasterband::GDALRasterBandH)
    xy = Array(Cint, 2); x = pointer(xy); y = x + sizeof(Cint)
    _getblocksize(rasterband, x, y)
    xy
end

"Fetch the pixel data type for this band."
_getrasterdatatype(rasterband::GDALRasterBandH) = GDALGetRasterDataType(rasterband)

dtype(rasterband::GDALRasterBandH) = _jltype(_getrasterdatatype(rasterband))

"""
Read/write a region of image data for this band.

This method allows reading a region of a `GDALRasterBand` into a buffer, or
writing data from a buffer into a region of a `GDALRasterBand`. It
automatically takes care of data type translation if the data type (`eBufType`)
of the buffer is different than that of the `GDALRasterBand`. The method also
takes care of image decimation / replication if the buffer size
`(nBufXSize x nBufYSize)` is different than the size of the region being
accessed `(nXSize x nYSize)`.

The `nPixelSpace` and `nLineSpace` parameters allow reading into or writing
from unusually organized buffers. This is primarily used for buffers containing
more than one bands raster data in interleaved format.

Some formats may efficiently implement decimation into a buffer by reading from
lower resolution overview images.

For highest performance full resolution data access, read and write on "block
boundaries" returned by `GetBlockSize()`, or use the `ReadBlock()` and
`WriteBlock()` methods.

### Parameters
* `eRWFlag`     Either GF_Read to read a region of data, or GF_Write to write a
region of data.
* `nXOff`       The pixel offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the left side.
* `nYOff`       The line offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the top.
* `nXSize`      The width of the region of the band to be accessed in pixels.
* `nYSize`      The height of the region of the band to be accessed in lines.
* `pData`       The buffer into which the data should be read, or from which it
should be written. This buffer must contain at least `(nBufXSize * nBufYSize)`
words of type `eBufType`. It is organized in left to right, top to bottom pixel
order. Spacing is controlled by the `nPixelSpace`, and `nLineSpace` parameters.
* `nBXSize`     The width of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `nBYSize`     The height of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `eBufType`    The type of the pixel values in the `buffer`. The pixel values
will be auto-translated to/from the `GDALRasterBand` data type as needed.
* `nPixelSpace` The byte offset from the start of one pixel value in `buffer`
to the start of the next pixel value within a scanline. If defaulted (0) the
size of the datatype `eBufType` is used.
* `nLineSpace`  The byte offset from the start of one scanline in `buffer` to
the start of the next. If defaulted (0) the size of the datatype
`(eBufType * nBufXSize)` is used.
* `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg structure
with additional arguments to specify resampling and progress callback, or
`NULL` for default behaviour. The `GDAL_RASTERIO_RESAMPLING` configuration
option can also be defined to override the default resampling to one of
`BILINEAR`, `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.

### Returns
`CE_Failure` if the access fails, otherwise `CE_None`.
"""
_rasterio(hRBand::GDALRasterBandH,
          eRWFlag::GDALRWFlag,
          nXOff::Integer,
          nYOff::Integer,
          nXSize::Integer,
          nYSize::Integer,
          pData::Ptr{Void},
          nBXSize::Integer,
          nBYSize::Integer,
          eBufType::GDALDataType,
          nPixelSpace::Integer,
          nLineSpace::Integer) =
    GDALRasterIO(hRBand,eRWFlag,nXOff,nYOff,nXSize,nYSize,pData,nBXSize,
                 nBYSize,eBufType,nPixelSpace,nLineSpace)::CPLErr

function rasterio!{T <: Real}(rasterband::GDALRasterBandH,
                              buffer::Array{T,2},
                              width::Integer,
                              height::Integer,
                              xoffset::Integer = 0,
                              yoffset::Integer = 0,
                              access::GDALRWFlag = GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0)
    xsize, ysize = size(buffer)
    io_error = _rasterio(rasterband, access, xoffset, yoffset, width,
                         height, Ptr{Void}(pointer(buffer)), xsize,
                         ysize, _gdaltype(T), nPixelSpace, nLineSpace)
    (io_error == CE_Failure) && error("Failed to access raster band")
    buffer
end

function rasterio!{T <: Real}(rasterband::GDALRasterBandH,
                              buffer::Array{T,2},
                              access::GDALRWFlag = GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0)
    rasterio!(rasterband, buffer, _getrasterbandxsize(rasterband),
              _getrasterbandysize(rasterband), 0, 0, access,
              nPixelSpace, nLineSpace)
end

function rasterio!{T <: Real, U <: Integer}(rasterband::GDALRasterBandH,
                                            buffer::Array{T,2},
                                            rows::UnitRange{U},
                                            cols::UnitRange{U},
                                            access::GDALRWFlag = GF_Read,
                                            nPixelSpace::Integer = 0,
                                            nLineSpace::Integer = 0)
    width = cols[end] - cols[1] + 1
    width < 0 && error("invalid window width")
    height = rows[end] - rows[1] + 1
    height < 0 && error("invalid window height")
    rasterio!(rasterband, buffer, width, height, cols[1], rows[1], access,
              nPixelSpace, nLineSpace)
end

"Fetch the width in pixels of this band."
_getrasterbandxsize(rasterband::GDALRasterBandH) =
    GDALGetRasterBandXSize(rasterband)::Cint

"Fetch the height in pixels of this band."
_getrasterbandysize(rasterband::GDALRasterBandH) =
    GDALGetRasterBandYSize(rasterband)::Cint

"Find out if we have update permission for this band."
_getrasteraccess(band::GDALRasterBandH) = GDALGetRasterAccess(band)::GDALAccess

"""Fetch the band number (1+) within its dataset, or 0 if unknown.

This method may return a value of 0 to indicate overviews, or free-standing
`GDALRasterBand` objects without a relationship to a dataset.
"""
_getbandnumber(rasterband::GDALRasterBandH) =
    GDALGetBandNumber(rasterband)::Cint

"""
Fetch the handle to its dataset handle, or `NULL` if this cannot be determined.

Note that some `GDALRasterBands` are not considered to be a part of a dataset,
such as overviews or other "freestanding" bands.
"""

_getbanddataset(rasterband::GDALRasterBandH) = GDALGetBandDataset(rasterband)

"""
Return raster unit type.

Return a name for the units of this raster's values. For instance, it might be
"m" for an elevation model in meters, or "ft" for feet. If no units are
available, a value of "" will be returned. The returned string should not be
modified, nor freed by the calling application.

### Returns
unit name string.
"""
_getrasterunittype(hBand::GDALRasterBandH) =
    GDALGetRasterUnitType(hBand)::Ptr{UInt8}

getrasterunittype(hBand::GDALRasterBandH) =
    unsafe_string(_getrasterunittype(hBand))

"Set unit type."
_setrasterunittype(hBand::GDALRasterBandH, newvalue::Ptr{UInt8}) =
    GDALSetRasterUnitType(hBand, newvalue)::CPLErr

function setrasterunittype!(hBand::GDALRasterBandH, newvalue::String)
    result = _setrasterunittype(hBand, pointer(newvalue))
    (result == CE_Failure) && error("Failed to set raster offset")
end

"""
Fetch the raster value offset.

This (in combination with `GetScale()`) is used to transform raw pixel values
into the units returned by `GetUnits()`. For e.g. this might be used to store
elevations in `GUInt16` bands with a precision of 0.1, starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically, a value of 0 is returned.

### Parameters
* `pbSuccess`   pointer to a boolean to use to indicate if the returned value
is meaningful or not. May be `NULL` (default).

### Returns
the raster offset.
"""
_getrasteroffset(hBand::GDALRasterBandH, pbSuccess::Ptr{Cint}) =
    GDALGetRasterOffset(hBand, pbSuccess)::Cdouble

function getrasteroffset(rasterband::GDALRasterBandH)
    success = Ref{Cint}()
    offset = _getrasteroffset(rasterband, success)
    (offset, Bool(success[]))
end

"Set scaling offset."
_setrasteroffset(band::GDALRasterBandH, offset::Cdouble) =
    GDALSetRasterOffset(band, offset)::CPLErr

function setrasteroffset!(hBand::GDALRasterBandH, offset::Cdouble)
    result = _setrasteroffset(hBand, offset)
    (result == CE_Failure) && error("Failed to set raster offset")
end

"""
Fetch the raster value scale.

This value (in combination with the `GetOffset()` value) is used to transform
raw pixel values into the units returned by `GetUnits()`. For example this
might be used to store elevations in GUInt16 bands with a precision of 0.1,
and starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically a value of one is returned.

### Parameters
* `pbSuccess`   pointer to a boolean to use to indicate if the returned value
is meaningful or not. May be `NULL` (default).

### Returns
the raster scale.
"""
_getrasterscale(hBand::GDALRasterBandH, pbSuccess::Ptr{Cint}) =
    GDALGetRasterScale(hBand, pbSuccess)::Cdouble

function getrasterscale(rasterband::GDALRasterBandH)
    success = Ref{Cint}()
    scale = _getrasterscale(rasterband, success)
    (scale, Bool(success[]))
end

"Set scaling ratio."
_setrasterscale(band::GDALRasterBandH, scale::Cdouble) =
    GDALSetRasterScale(band, scale)::CPLErr

function setrasterscale!(hBand::GDALRasterBandH, scale::Cdouble)
    result = _setrasterscale(hBand, scale)
    (result == CE_Failure) && error("Failed to set raster scale")
end

"""
Copy all raster band raster data.

This function copies the complete raster contents of one band to another
similarly configured band. The source and destination bands must have the same
width and height. The bands do not have to have the same data type.

It implements efficient copying, in particular "chunking" the copy in
substantial blocks.

Currently the only `papszOptions` value supported is : "COMPRESSED=YES" to
force alignment on target dataset block sizes to achieve best compression.
More options may be supported in the future.

### Parameters
* `hSrcBand`        the source band
* `hDstBand`        the destination band
* `papszOptions`    transfer hints in "StringList" Name=Value format.
* `pfnProgress`     progress reporting function.
* `pProgressData`   callback data for progress function.

### Returns
`CE_None` on success, or `CE_Failure` on failure.
"""
_rasterbandcopywholeraster(hSrcBand::GDALRasterBandH,
                           hDstBand::GDALRasterBandH,
                           papszOptions::Ptr{Ptr{UInt8}},
                           pfnProgress::GDALProgressFunc,
                           pProgressData::Ptr{Void}) =
    GDALRasterBandCopyWholeRaster(hSrcBand, hDstBand, papszOptions,
                                  pfnProgress, pProgressData)::CPLErr

function rasterbandcopywholeraster(hSrcBand::GDALRasterBandH,
                                   hDstBand::GDALRasterBandH,
                                   options::Vector{String})
    result = _rasterbandcopywholeraster(hSrcBand, hDstBand,
                                        Ptr{Ptr{UInt8}}(pointer(options)),
                                        C_NULL, C_NULL)
    (result == CE_Failure) && error("Failed to copy raster band")
end

function rasterbandcopywholeraster(hSrcBand::GDALRasterBandH,
                                   hDstBand::GDALRasterBandH)
    result = _rasterbandcopywholeraster(hSrcBand, hDstBand,
                                        Ptr{Ptr{UInt8}}(C_NULL),
                                        C_NULL, C_NULL)
    (result == CE_Failure) && error("Failed to copy raster band")
end

"""
Generate downsampled overviews.

This function will generate one or more overview images from a base image using
the requested downsampling algorithm. Its primary use is for generating
overviews via BuildOverviews(), but it can also be used to generate downsampled
images in one file from another outside the overview architecture.

The output bands need to exist in advance.

This function will honour properly `NODATA_VALUES` tuples (special dataset
metadata) so that only a given RGB triplet (in case of a RGB image) will be
considered as the nodata value and not each value of the triplet independantly
per band.

### Parameters
* `hSrcBand`        the source (base level) band.
* `nOverviewCount`  the number of downsampled bands being generated.
* `pahOvrBands`     the list of downsampled bands to be generated.
* `pszResampling`   Resampling algorithm (eg. "AVERAGE").
* `pfnProgress`     progress report function.
* `pProgressData`   progress function callback data.

### Returns
`CE_None` on success or `CE_Failure` on failure.
"""
_regenerateoverviews(hSrcBand::GDALRasterBandH,
                     nOverviewCount::Integer,
                     pahOvrBands::Ptr{GDALRasterBandH},
                     pszResampling::Ptr{UInt8},
                     pfnProgress::GDALProgressFunc,
                     pProgressData::Ptr{Void}) =
    GDALRegenerateOverviews(hSrcBand, nOverviewCount, pahOverviewBands,
                            pszResampling, pfnProgress, pProgressData)::CPLErr

"Advise driver of upcoming read requests."
_rasteradviseread(hRB::GDALRasterBandH,
                  nDSXOff::Integer,
                  nDSYOff::Integer,
                  nDSXSize::Integer,
                  nDSYSize::Integer,
                  nBXSize::Integer,
                  nBYSize::Integer,
                  eBDataType::GDALDataType,
                  papszOptions::Ptr{Ptr{UInt8}}) =
    GDALRasterAdviseRead(hRB, nDSXOff, nDSYOff, nDSXSize, nDSYSize, nBXSize,
                         nBYSize, eBDataType, papszOptions)::CPLErr

"""
Read a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use RasterIO().

See the `GetLockedBlockRef()` method for a way of accessing internally cached
block oriented data without an extra copy into an application buffer.

### Parameters
* `nXBlockOff`  the horizontal block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `nYBlockOff`  the vertical block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `pImage`      the buffer into which the data will be read. The buffer must be
large enough to hold `GetBlockXSize()*GetBlockYSize()` words of type
`GetRasterDataType()`.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
_readblock(band::GDALRasterBandH,
           nXBlockOff::Integer,
           nYBlockOff::Integer,
           pImage::Ptr{Void}) =
    GDALReadBlock(band, nXBlockOff, nYBlockOff, pImage)::CPLErr

"""
Write a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use `RasterIO()`.

See `ReadBlock()` for an example of block oriented data access.

### Parameters
* `nXBlockOff`  the horizontal block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `nYBlockOff`  the vertical block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `pImage`      the buffer from which the data will be written. The buffer must
be large enough to hold `GetBlockXSize()*GetBlockYSize()` words of type
`GetRasterDataType()`.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
_writeblock(band::GDALRasterBandH,
            nXBlockOff::Integer,
            nYBlockOff::Integer,
            pImage::Ptr{Void}) =
    GDALWriteBlock(band, nXBlockOff, nYBlockOff, pImage)::CPLErr

"""
Check for arbitrary overviews.

This returns `TRUE` if the underlying datastore can compute arbitrary overviews
efficiently, e.g. with OGDI over a network. Datastores with arbitrary overviews
don't generally have any fixed overviews, but the `RasterIO()` method can be
used in downsampling mode to get overview data efficiently.

### Returns
`TRUE` if arbitrary overviews available (efficiently), otherwise `FALSE`.
"""
_hasarbitraryoverviews(band::GDALRasterBandH) =
    GDALHasArbitraryOverviews(band)::Cint

"""
Return the number of overview layers available.

### Returns
overview count, zero if none.
"""
_overviewcount(band::GDALRasterBandH) = GDALGetOverviewCount(band)::Cint

"""
Fetch overview raster band object.

### Parameters
* `i`   overview index between `0` and `GetOverviewCount()-1`.

### Returns
overview GDALRasterBand.
"""
_getoverview(band::GDALRasterBandH, i::Integer) =
    GDALGetOverview(band, i)::GDALRasterBandH

getoverview(band::GDALRasterBandH, i::Int) = _getoverview(band, i-1)

"""
Fetch best sampling overview.

Returns the most reduced overview of the given band that still satisfies the
desired number of samples. This function can be used with zero as the number
of desired samples to fetch the most reduced overview. The same band as was
passed in will be returned if it has not overviews, or if none of the overviews
have enough samples.

### Parameters
* `nDesiredSamples`     the returned band will have at least this many pixels.

### Returns
optimal overview or the band itself.
"""
_getrastersampleoverview(band::GDALRasterBandH, nDesiredSamples::Integer) =
    GDALGetRasterSampleOverview(band, nDesiredSamples)::GDALRasterBandH

"""
Fetch the no data value for this band.

If there is no out of data value, an out of range value will generally be
returned. The no data value for a band is generally a special marker value
used to mark pixels that are not valid data. Such pixels should generally
not be displayed, nor contribute to analysis operations.

### Parameters
* `pbSuccess`   pointer to a boolean to use to indicate if a value is actually
associated with this layer. May be `NULL` (default).

### Returns
the nodata value for this band.
"""
_getrasternodatavalue(band::GDALRasterBandH, pbSuccess::Ptr{Cint}) =
    GDALGetRasterNoDataValue(band, pbSuccess)::Cdouble

"Set the no data value for this band."
_setrasternodatavalue(band::GDALRasterBandH, value::Cdouble) =
    GDALSetRasterNoDataValue(band, value)::CPLErr

"""
Fetch the list of category names for this raster.

The return list is a "StringList" in the sense of the CPL functions. That is a
NULL terminated array of strings. Raster values without associated names will
have an empty string in the returned list. The first entry in the list is for
raster values of zero, and so on.

The returned stringlist should not be altered or freed by the application.
It may change on the next GDAL call, so please copy it if it is needed for any
period of time.

### Returns
list of names, or `NULL` if none.
"""
_getrastercategorynames(band::GDALRasterBandH) =
    GDALGetRasterCategoryNames(band)::Ptr{Ptr{UInt8}}

"Set the category names for this band."
_setrastercategorynames(band::GDALRasterBandH, names::Ptr{Ptr{UInt8}}) =
    GDALSetRasterCategoryNames(band, names)::CPLErr

"""
Flush raster data cache.

This call will recover memory used to cache data blocks for this raster band,
and ensure that new requests are referred to the underlying driver.

### Returns
`CE_None` on success.
"""
_flushrastercache(band::GDALRasterBandH) = GDALFlushRasterCache(band)::CPLErr

"""
Fill this band with a constant value.

GDAL makes no guarantees about what values pixels in newly created files are
set to, so this method can be used to clear a band to a specified "default"
value. The fill value is passed in as a double but this will be converted to
the underlying type before writing to the file. An optional second argument
allows the imaginary component of a complex constant value to be specified.

### Parameters
* `dfRealValue`         Real component of fill value
* `dfImaginaryValue`    Imaginary component of fill value, defaults to zero

### Returns
`CE_Failure` if the write fails, otherwise `CE_None`
"""
_fillraster(band::GDALRasterBandH,
            dfRealValue::Cdouble,
            dfImaginaryValue::Cdouble = 0) =
    GDALFillRaster(band, dfRealValue, dfImaginaryValue)::CPLErr

"""
Return the mask band associated with the band. (Since: GDAL 1.6.0)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
GDALNoDataValuesMaskBand class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will
return `GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new
`GDALAllValidRasterBand` class will be returned that has 255 values for all
pixels. The null flags will return `GMF_ALL_VALID`.

Note that the `GetMaskBand()` should always return a `GDALRasterBand` mask,
even if it is only an all 255 mask with the flags indicating `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
_getmaskband(band::GDALRasterBandH) = GDALGetMaskBand(band)::GDALRasterBandH

"""
Return the status flags of the mask band associated with the band.

The GetMaskFlags() method returns an bitwise OR-ed set of status flags with the
following available definitions that may be extended in the future:

* `GMF_ALL_VALID` (`0x01`):    There are no invalid pixels, all mask values
will be 255. When used this will normally be the only flag set.
* `GMF_PER_DATASET` (`0x02`):  The mask band is shared between all bands on
the dataset.
- `GMF_ALPHA` (`0x04`):        The mask band is actually an alpha band and may
have values other than 0 and 255.
- `GMF_NODATA` (`0x08`):       Indicates the mask is actually being generated
from nodata values. (mutually exclusive of `GMF_ALPHA`)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
`GDALNoDataValuesMaskBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags `GMF_PER_DATASET`
and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new `GDALAllValidRasterBand`
class will be returned that has 255 values for all pixels. The null flags will
return `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
_getmaskflags(band::GDALRasterBandH) = GDALGetMaskFlags(band)::Cint

"""
Adds a mask band to the current band.

The default implementation of the `CreateMaskBand()` method is implemented
based on similar rules to the `.ovr` handling implemented using the
`GDALDefaultOverviews` object. A `TIFF` file with the extension `.msk` will be
created with the same basename as the original file, and it will have as many
bands as the original image (or just one for `GMF_PER_DATASET`). The mask
images will be deflate compressed tiled images with the same block size as the
original image if possible.

If you got a mask band with a previous call to `GetMaskBand()`, it might be
invalidated by `CreateMaskBand()`. So you have to call `GetMaskBand()` again.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
_createmaskband(band::GDALRasterBandH,nFlags::Integer) =
    GDALCreateMaskBand(band, nFlags)::CPLErr


text:
function compute_parameter_jacobian(obj::NamedParameterJacobianProvider, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::JString)
    return jcall(obj, "computeParameterJacobian", Vector{jdouble}, (jdouble, Vector{jdouble}, Vector{jdouble}, JString), arg0, arg1, arg2, arg3)
end

function get_parameters_names(obj::Parameterizable)
    return jcall(obj, "getParametersNames", List, ())
end

function is_supported(obj::Parameterizable, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end


text:
@extension Cedrus begin
"""
[Extension Website](https://github.com/haberdashPI/WeberCedrus.jl)

    @Cedrus()

Creates an extension for Weber allowing experiments to respond to events from
Cedrus response-pad hardware. You can use [`iskeydown`](@ref) and
[`iskeyup`](@ref) to check for events. To find the keycodes of the buttons for
your response pad, run the following code, and press each of the buttons on the
response pad.

    run_keycode_helper(extensions=[@Cedrus()])

"""
end

@extension DAQmx begin
"""
[Extension Website](https://github.com/haberdashPI/WeberDAQmx.jl)

    @DAQmx(port;eeg_sample_rate,[codes])

Create a Weber extension that writes `record` events to a digital
out line via the DAQmx API. This can be used to send trigger
codes during eeg recording.

# Arguments

* port: should be `nothing`, to disable the extension, or
  the port name for the digital output line.
* eeg_sample_rate: should be set to the sampling rate for
  eeg recording. This calibrates the code length for triggers.
* codes: a Dict that maps record event codes (a string) to a number.
  This should be an Integer less than 256. Any codes not
  specified here will be automatically set, based on the order
  in which codes are recieved.

# Example

The following experiment sends the code 0x01 to port0 on TestDevice.

    port = "/TestDevice/port0/line0:7"
    experiment = Experiment(extensions=[
      @DAQmx(port;eeg_sample_rate=512,codes=Dict("test" => 0x01))])
    setup(experiment) do
      addtrial(moment(record,"test"))
    end
    run(experiment)
"""
end


text:
# PROX_QUAD    The proximal operator of a quadratic.
#
#   prox_quad(v,lambda,A,b)
function prox_quad(
	v,
	lambda::Real,
	A,
	b
)
    rho = 1 / lambda
    n = size(A)

    if issparse(A)
        x = (A + rho * speye(m)) \ (rho * v .- b);
    else
        x = (A + rho * eye(m)) \ (rho * v .- b);
    end

    return x
end


text:
using HMCExamples

# The main purpose is for tracing the compiler, not for execution speed or results.  Keep the number of samples small
HMCExamples.main_rbc_1_kalman(["--num_samples", "30", "--max_depth", "2"])
HMCExamples.main_rbc_1_kalman(["--num_samples", "30", "--max_depth", "2", "--print_level", "2"])
HMCExamples.main_rbc_1_joint(["--num_samples", "30", "--max_depth", "2"])
HMCExamples.main_rbc_2_joint(["--num_samples", "30", "--max_depth", "2"])
HMCExamples.main_FVGQ_1_kalman(["--num_samples", "20", "--print_level", "2"])
HMCExamples.main_FVGQ_1_joint(["--num_samples", "20", "--print_level", "3"])
HMCExamples.main_FVGQ_2_joint(["--num_samples", "10", "--max_depth", "2"])

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples