text:
import BEAST; BE = BEAST;
using CompScienceMeshes
using Test

T = Float64
P = SVector{3,T}
tol = eps(T)*10^3

struct IntegrationDomain{T,P}
    signed_height::T
    dist_to_plane::T
    center::P
    plane_center::P
    unit_normal::P
    segments::Vector{SVector{2,P}}
end

function IntegrationDomain(p1, p2, p3, center)

    t1 = p2 - p1
    t2 = p3 - p1

    unit_normal = t1 × t2
    unit_normal /= norm(unit_normal)
    signed_height = dot(center-p1, unit_normal)
    dist_to_plane = abs(signed_height)
    plane_center = center - signed_height * unit_normal

    s1 = SVector((p1, p2,))
    s2 = SVector((p2, p3,))
    s3 = SVector((p3, p1,))
    segments = [s1, s2, s3]

    #T = eltype(p1)
    IntegrationDomain(
        signed_height,
        dist_to_plane,
        center,
        plane_center,
        unit_normal,
        segments)

end

wiltonints(a,b,c,center) = wiltonints(IntegrationDomain(a,b,c,center))

# type WiltonIntsWS{T}
#     segments_through_c::Vector{T}
# end

# type WiltonIntsRV{T,N,M}
#     scalars::SVector{N,T}
#     gradients::SVector{M,T}
# end

# This function implements the Wilton algorithm for Integration of
# linear potentials on triangular domains
function wiltonints(domain::IntegrationDomain{T,Q}, tol=eps(T)*10^3) where {T,Q}

    dim = 3

    d = domain.signed_height
    D = domain.dist_to_plane
    c = domain.plane_center
    n = domain.unit_normal

    num_segments = length(domain.segments)

    # outside the domain of integration, or whether it is on the boundary.
    # keep track of whether the projected singularity is inside or
    # If it is on the boundary, make a note of on which part of the boundary
    # it is located.
    c_inside = true
    num_segments_through_c = 0
    segments_through_c = zeros(Int, num_segments)

    s = zeros(T,6) # TODO alloc
    v = zeros(Q, 6) # TODO alloc

    for i in 1 : num_segments

        a = domain.segments[i][1]
        b = domain.segments[i][2]

        ab = b - a     # segment tangent
        l = norm(ab)   # segment length
        t = ab / l     # unit tangent
        m = cross(t,n) # unit binormal

        ca = a - c
        cb = b - c

        # coordinates of a and b w.r.t. the projection
        # of c on the line [a,b]
        xa = dot(ca,t)
        xb = dot(cb,t)

        # signed and absolute distance of plane center to line
        p = dot(ca,m)
        P = abs(p)

        if p < 0
          c_inside = false
        end

        w2 = p^2 + d^2
        ra = sqrt(xa^2 + w2)
        rb = sqrt(xb^2 + w2)
        w = sqrt(w2)

        # if the plane center is NOT on the line through this segment
        # use the general formula; no singularities are encountered in
        # the computation
        if P > tol
            lg = log((xb+rb) / (xa+ra))
            at = atan(d*xb/(p*rb)) - atan(d*xa/(p*ra))

            s[1] += -at
            s[2] += p*lg + d*at
            s[3] += 1/2 * p * (xb-xa)
            pt1 = 1/6*p*(xb*rb-xa*ra)
            pt2 = (1/2*d*d+1/6*p*p)*p*lg
            pt3 = 1/3*d*d*d*at
            s[4] += pt1 + pt2 + pt3
            pt1 = 1/2*xb*p*(1/6*xb*xb+1/2*p*p+d*d)
            pt2 = 1/2*xa*p*(1/6*xa*xa+1/2*p*p+d*d)
            s[5] += pt1 - pt2
            pt1b = 1/20*p*xb*rb*(xb*xb+1/2*(9*d*d+5*p*p))
            pt1a = 1/20*p*xa*ra*(xa*xa+1/2*(9*d*d+5*p*p))
            pt2 = 1/40*p*(15*d*d*d*d+10*d*d*p*p+3*p*p*p*p)*lg + 0.2*d*d*d*d*d*at;
            s[6] += (pt1b - pt1a) + pt2

            v[1] += -lg * m
            v[2] += 0.5 *((xb*rb-xa*ra) + w*w*lg) * m
            pt1b = 0.5 * xb * (1/3*xb*xb + p*p)
            pt1a = 0.5 * xa * (1/3*xa*xa + p*p)
            v[3] += (pt1b - pt1a) * m
            pt1 = 1/12*(xb*rb*rb*rb-xa*ra*ra*ra) + 1/8*w*w*(xb*rb-xa*ra)
            pt2 = 1/8*w*w*w*w*lg
            v[4] += (pt1 + pt2) * m
            pt1b = 1/20*xb*xb*xb*xb*xb + 1/6*p*p*xb*xb*xb + 0.25*p*p*p*p*xb + 0.5*d*d*(1/3*xb*xb*xb + p*p*xb)
            pt1a = 1/20*xa*xa*xa*xa*xa + 1/6*p*p*xa*xa*xa + 0.25*p*p*p*p*xa + 0.5*d*d*(1/3*xa*xa*xa + p*p*xa)
            v[5] += (pt1b - pt1a) * m

        else # plane center lies on the line ab

            if xa*xb < tol^2
                # projection c on ab lies in between a and b
                num_segments_through_c += 1
                segments_through_c[num_segments_through_c] = i
            end

            if D > tol
                # singularity not in the plane of the triangle and thus not on ab

                lg = log((xb+rb)/(xa+ra))

                v[1] += -lg * m
                pt1 = 0.5*(xb*rb-xa*ra)
                pt2 = 0.5*w*w*lg
                v[2] += (pt1 + pt2) * m
                pt1b = 1/6*xb*xb*xb
                pt1a = 1/6*xa*xa*xa
                v[3] += (pt1b - pt1a) * m
                pt1 = 1/12*(xb*rb*rb*rb-xa*ra*ra*ra) + 1/8*w*w*(xb*rb-xa*ra)
                pt2 = 1/8*w*w*w*w*lg
                v[4] += (pt1 + pt2) * m
                pt1b = 1/20*xb*xb*xb*xb*xb + 1/6*p*p*xb*xb*xb + 0.25*p*p*p*p*xb + 0.5*d*d*(1/3*xb*xb*xb + p*p*xb)
                pt1a = 1/20*xa*xa*xa*xa*xa + 1/6*p*p*xa*xa*xa + 0.25*p*p*p*p*xa + 0.5*d*d*(1/3*xa*xa*xa + p*p*xa)
                v[5] += (pt1b - pt1a) * m

            else # D <= tol

                # the projected singularity lies on ab. It is assumed the singularity
                # itself does not lie on [a,b] however, since this would lead to some
                # of the integrals to become infinite.

                lg = (xb > 0) ? log(xb/xa) : log(xa/xb)

                v[1] += -lg * m
                v[2] += 0.5*(xb*rb-xa*ra) * m
                pt1b = 1/6*xb*xb*xb
                pt1a = 1/6*xa*xa*xa
                v[3] += (pt1b - pt1a) * m
                pt1 = 1/12*(xb*rb*rb*rb-xa*ra*ra*ra) + 1/8*w*w*(xb*rb-xa*ra)
                pt2 = 1/8*w*w*w*w*lg
                v[4] += (pt1 + pt2) * m
                pt1b = 1/20*xb*xb*xb*xb*xb + 1/6*p*p*xb*xb*xb + 0.25*p*p*p*p*xb + 0.5*d*d*(1/3*xb*xb*xb + p*p*xb)
                pt1a = 1/20*xa*xa*xa*xa*xa + 1/6*p*p*xa*xa*xa + 0.25*p*p*p*p*xa + 0.5*d*d*(1/3*xa*xa*xa + p*p*xa)
                v[5] += (pt1b - pt1a) * m

            end # if D < tol

        end # if p < tol

    end # next segment

    # contributions from intersection of epsilon cirlce around c with triangle
    sgn = D < tol ? zero(d) : sign(d)
    alpha = zero(d)
    if num_segments_through_c == 0 && c_inside
        alpha = 2pi
    elseif num_segments_through_c == 1
        alpha = pi
    elseif num_segments_through_c == 2
        i = segments_through_c[1]
        t1 = domain.segments[i][2] - domain.segments[i][1]
        i = segments_through_c[2]
        t2 = domain.segments[i][2] - domain.segments[i][1]
        cosalpha = dot(-t1,t2)/(norm(t1)*norm(t2))
        alpha = acos(clamp(cosalpha,-1,+1))
    elseif num_segments_through_c > 2
        error("num_intersections must be either 0, 1, or 2.")
    end

    s[1] += alpha * sgn
    s[2] += -alpha * D
    s[4] += -1/3 * alpha * D*D*D
    s[6] += -1/5 * alpha * D*D*D*D*D

    # build gradients
    g = Vector{Q}(undef,4)
    g[1] = -1( v[1] - s[1] * n )
    g[2] = +1( v[2] - d*s[2]*n )
    g[3] = +2( v[3] - d*s[3]*n )
    g[4] = +3( v[4] - d*s[4]*n )
    return s, v, g

end

# s1 = int d/R^3
# s2 = int 1/R
# s3 = int 1
# s4 = int R
# s5 = int R^2
# s6 = int R^3
# v1 = int (rho'-rho)/R^3
# v2 = int (rho'-rho)/R
# v3 = int (rho'-rho)
# v4 = int (rho'-rho)*R
# v5 = int (rho'-rho)*R
# g1 = int grad' (1/R)
# g2 = int grad' R
# g3 = int grad' R^2
# g4 = int grad' R^3
function withquadrules(triangle, r, n)

	u, w = BE.trgauss(n)

    n = (triangle[1] - triangle[3]) × (triangle[2] - triangle[3])
    n /= norm(n)
    d = (r - triangle[1]) ⋅ n
    ρ = r - d*n

    scalar = zeros(T,6)
    vector = zeros(P,6)
    gradgr = zeros(P,4)

	m = simplex(triangle[1], triangle[2], triangle[3])

    for i in 1 : length(w)
        q = neighborhood(m, u[:,i])
        s = cartesian(q)
        dq = w[i] * jacobian(q)

        R = norm(s - r)
        scalar[1] += d / R^3 * dq
        scalar[2] += 1 / R * dq
        scalar[3] += 1 * dq
        scalar[4] += R * dq
        scalar[5] += R^2 * dq
        scalar[6] += R^3 * dq
        vector[1] += (s - ρ) / R^3 * dq
        vector[2] += (s - ρ) / R * dq
        vector[3] += (s - ρ) * dq
        vector[4] += (s - ρ) * R * dq
        vector[5] += (s - ρ) * R^2 * dq
    end

    gradgr[1] = -1( vector[1] - scalar[1] * n )
    gradgr[2] = +1( vector[2] - d*scalar[2]*n )
    gradgr[3] = +2( vector[3] - d*scalar[3]*n )
    gradgr[4] = +3( vector[4] - d*scalar[4]*n )

    return scalar, vector, gradgr
end

tr = [
	point(0.0, 0.0, 0.0),
	point(2.0, 0.0, 0.0),
	point(0.0, 3.0, 0.0)]

pts = [
    (tr[1]+2*tr[2]-1*tr[3]),
    (tr[1]+tr[2])/2 + point(0.0, 0.0, 0.3),
    point(3.0, 0.0, 0.0),
    (tr[1]+tr[2]+tr[3])/3 + point(0.0, 0.0, 0.3),
    tr[2] + point(0.0, 0.0, 0.3),
    (tr[1]+tr[2]+tr[3])/3 + point(0.0, 0.0, 100.0)
]

for _p in pts
    s1, v1, g1 = wiltonints(tr[1],tr[2],tr[3],_p)
    s2, v2, g2 = withquadrules(tr,_p,13)

    for i in length(s1)
        @test norm(s1[i]-s2[i]) < 1.0e-5
    end

    for i in length(v1)
        @test norm(v1[i]-v2[i]) < 1.0e-5
    end

    for i in 1:length(g1)
        for _j in 1:3
            @test norm(g1[i][_j]-g2[i][_j]) < 1.0e-5
        end
    end
end

Γ = readmesh(joinpath(dirname(@__FILE__),"assets","sphere2.in"))
nc = numcells(Γ)
t = chart(Γ, first(cells(Γ)))
s = chart(Γ, last(cells(Γ)))

X = BEAST.raviartthomas(Γ)
x = BEAST.refspace(X)

κ = 1.0
op = BEAST.MWSingleLayer3D(κ)

n = BE.numfunctions(x)
z1 = zeros(ComplexF64, n, n)
z2 = zeros(ComplexF64, n, n)

BE = BEAST

tqd = BE.quadpoints(x, [t], (12,))
bqd = BE.quadpoints(x, [s], (13,))

DQ_strategy = BE.DoubleQuadRule(tqd[1,1], bqd[1,1])
BEAST.momintegrals!(op, x, x, t, s, z1, DQ_strategy)

SE_strategy = BE.WiltonSERule(
  tqd[1,1],
  BE.DoubleQuadRule(
	tqd[1,1],
	bqd[1,1],
  ),
)
BEAST.momintegrals!(op, x, x, t, s, z2, SE_strategy)

@test norm(z1-z2) < 1.0e-7


text:
module DirectConvolution

const RootDir = pathof(DirectConvolution)

using StaticArrays
using LinearAlgebra

include("utils.jl")
include("linearFilter.jl")
include("directConvolution.jl")
include("SG_Filter.jl")
include("udwt.jl")

end # module

ext: jl


text:
module ODDS

    using OutlierDetection: CLASS_NORMAL, CLASS_OUTLIER, to_categorical
    using DataFrames: DataFrame, convert
    using DataDeps
    using MAT

    export list, load

    const ODDS_DATASETS = [
        (
            name = "annthyroid",
            md5 = "17a1012d005f2179eca3b87c0028eb5ca807794b359c33fa1f8ca5ac846c514f",
            url = "https://www.dropbox.com/s/aifk51owxbogwav/annthyroid.mat?dl=1"
        ),
        (
            name = "arrhythmia",
            md5 = "e4d47e499591ae31e5d7f0cd983367092353eac240f48e14c9b36457edb3460d",
            url = "https://www.dropbox.com/s/lmlwuspn1sey48r/arrhythmia.mat?dl=1"
        ),
        (
            name = "breastw",
            md5 = "a2e3f7c7f61ddde87c4373f0c0539de996227b9659c70b3c227948ebf0c6062b",
            url = "https://www.dropbox.com/s/g3hlnucj71kfvq4/breastw.mat?dl=1"
        ),
        (
            name = "cardio",
            md5 = "5e568b6a73711012c481fa6cc984ab9c68e9a4fa44d6aeb3b43efa381c3f695b",
            url = "https://www.dropbox.com/s/galg3ihvxklf0qi/cardio.mat?dl=1"
        ),
        (
            name = "cover",
            md5 = "39acf6da4d31f08c2c28dc047249d8cd8ccd9d6fc299d9b5a08874ac92290c98",
            url = "https://www.dropbox.com/s/awx8iuzbu8dkxf1/cover.mat?dl=1"
        ),
        # the MAT import for ecoli is currently not working
        # (
        #     name = "ecoli",
        #     md5 =  "76886bf14f549680450a0720c58b461db965f8b762a5ab4d165c03bbdcb355fc",
        #     url =  "https://www.dropbox.com/s/f38mquv0g08xxgp/ecoli.mat?dl=1"
        # ),
        (
            name = "glass",
            md5 = "2c506e6e1e2ca1e82a905d10a90773151ff812c4bd969c71a12db5b2bb6e7921",
            url = "https://www.dropbox.com/s/iq3hjxw77gpbl7u/glass.mat?dl=1"
        ),
        (
            name = "http",
            md5 = "81d689cfec5407305ab98cd26794e2cefb3aef5a521bf9c59744107aabcaed7e",
            url = "https://www.dropbox.com/s/iy9ucsifal754tp/http.mat?dl=1"
        ),
        (
            name = "ionosphere",
            md5 = "8ac5dfa6f47db4acdfe10615489b7f7dd015d79d5b381ce641e72baca3cfefb0",
            url = "https://www.dropbox.com/s/lpn4z73fico4uup/ionosphere.mat?dl=1"
        ),
        (
            name = "letter",
            md5 = "f4a229b0e2a6b0c9157ce27fa760bcbac0f948cd690b2f5d51c120e5450f2195",
            url = "https://www.dropbox.com/s/rt9i95h9jywrtiy/letter.mat?dl=1"
        ),
        (
            name = "lympho",
            md5 = "031b907ecb50347720b08b38f8139d05635adeca9335a064d89d5569c27a5834",
            url = "https://www.dropbox.com/s/ag469ssk0lmctco/lympho.mat?dl=1"
        ),
        (
            name = "mammography",
            md5 = "271ebb568314a856666d3504b4882e21b0ea6e1ba9e648ad256d572a36df597e",
            url = "https://www.dropbox.com/s/tq2v4hhwyv17hlk/mammography.mat?dl=1"
        ),
        (
            name = "mnist",
            md5 = "9479e221f236684c09f3e11c7ad124725d34ecd301b9450425e8b113368413a6",
            url = "https://www.dropbox.com/s/n3wurjt8v9qi6nc/mnist.mat?dl=1"
        ),
        (
            name = "mulcross",
            md5 = "f64a10efb641229356e7ae06e8805f9025e1ce29a2b0abb1207447423746276b",
            url = "https://www.dropbox.com/s/tmfan6f7qf46bkx/mulcross.mat?dl=1"
        ),
        (
            name = "musk",
            md5 = "50bbf414caceb66cfe6f6b1eb7c1ff7c6eb1e6d87173f7bb597ba55cd3dc5493",
            url = "https://www.dropbox.com/s/we6aqhb0m38i60t/musk.mat?dl=1"
        ),
        (
            name = "optdigits",
            md5 = "30c664cd6f3740a173277a27cd1028d471b06f0835a9902c9712749cec58c1c4",
            url = "https://www.dropbox.com/s/w52ndgz5k75s514/optdigits.mat?dl=1"
        ),
        (
            name = "pendigits",
            md5 = "7bc32ac83b3d583443acea8b1342a54aaa59eaedb9996af405322f473f3c22bc",
            url = "https://www.dropbox.com/s/1x8rzb4a0lia6t1/pendigits.mat?dl=1"
        ),
        (
            name = "pima",
            md5 = "b6a5be8d548b9293d6958beceb1c77055591714d25d40d957eba87808456468f",
            url = "https://www.dropbox.com/s/mvlwu7p0nyk2a2r/pima.mat?dl=1"
        ),
        (
            name = "satellite",
            md5 = "6feac3112b9c14e1c3e60afc437f2f3d29dc1000119c9f950c628078778d6aa0",
            url = "https://www.dropbox.com/s/dpzxp8jyr9h93k5/satellite.mat?dl=1"
        ),
        (
            name = "satimage-2",
            md5 = "95c37a07cea2cb201f6900d14f521510e10b6b2f0074c06ac46d3157440ec67b",
            url = "https://www.dropbox.com/s/hckgvu9m6fs441p/satimage-2.mat?dl=1"
        ),
        (
            name = "shuttle",
            md5 = "b6e396c92354d5dcf5036901e36dc58e2913f2d1a5e138c35a2bedf9c69db899",
            url = "https://www.dropbox.com/s/mk8ozgisimfn3dw/shuttle.mat?dl=1"
        ),
        (
            name = "smtp",
            md5 = "5b6d13dd23e8c5e4235081ea58b606df216a2b9b86e7f5403cde26f62dafc0d8",
            url = "https://www.dropbox.com/s/dbv2u4830xri7og/smtp.mat?dl=1"
        ),
        (
            name = "speech",
            md5 = "193d85404c5aff8b33e7e0540a8be211a7c94948d1a33162ce500b7e526ea937",
            url = "https://www.dropbox.com/s/w6xv51ctea6uauc/speech.mat?dl=1"
        ),
        (
            name = "thyroid",
            md5 = "49155d7861b977e90808faf68343ff67a67f63090c2744d514dbb67ac68e7e33",
            url = "https://www.dropbox.com/s/bih0e15a0fukftb/thyroid.mat?dl=1"
        ),
        (
            name = "vertebral",
            md5 = "84b50b055f118f54b5c976650ccd284ab7c1adee822f8b3eb517c7e2c890d15e",
            url = "https://www.dropbox.com/s/5kuqb387sgvwmrb/vertebral.mat?dl=1"
        ),
        (
            name = "vowels",
            md5 = "51d5397bad1ceeda22b7cd44fd58125c5d603c9c2edb53aa59933330e259ac7b",
            url = "https://www.dropbox.com/s/pa26odoq6atq9vx/vowels.mat?dl=1"
        ),
        (
            name = "wbc",
            md5 = "afe6b26953cc196b9a480318c903e0edef22529c2f61bd02082438c6e70630a7",
            url = "https://www.dropbox.com/s/ebz9v9kdnvykzcb/wbc.mat?dl=1"
        ),
        (
            name = "wine",
            md5 = "d1b7cc643379e14397ee25bc1125ba08e4c929f2ec44d8782b1340d7045d4cf7",
            url = "https://www.dropbox.com/s/uvjaudt2uto7zal/wine.mat?dl=1"
        )
    ]

    list() = [name for (name, _, _) in ODDS_DATASETS]
    to_name(dataset::AbstractString) = "odds-$dataset"

    function load(dataset::AbstractString)
        @assert dataset in list()

        dep = @datadep_str to_name(dataset)
        data = MAT.matread("$dep/$dataset.mat")
        X = DataFrame(data["X"], :auto)
        y = ifelse.(data["y"][:, 1] .== 0, CLASS_NORMAL, CLASS_OUTLIER)
        X, to_categorical(y)
    end

    function __init__()
        for (name, md5, url) in ODDS_DATASETS
            register(DataDep(
                to_name(name),
                """Dataset: $name
                Collection: Outlier Detection DataSets (ODDS)
                Authors: Shebuti Rayana
                Website: http://odds.cs.stonybrook.edu""",
                url,
                md5
            ))
        end
    end

end


text:
using MultivariateAnomalies
#using LinearAlgebra
using Statistics
using Test

include("testFeatureExtr.jl")
include("testAlgorithms.jl")
include("testPostprocessing.jl")


text:
struct BlackBox()
end

struct BellScenario()
end

struct Polytope()
end


text:
include(normpath(joinpath(@__DIR__, "IOM", "src", "models", "EMOM", "EMOM.jl")))
using NCDatasets
using MPI
using Formatting
using ArgParse

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--config-file"
            help = "Config file used to to produce a snapshot file."
            arg_type = String
            required = true

    end

    return parse_args(s)
end

parsed = parse_commandline()

MPI.Init()

println("Processing data...")

include(parsed["config-file"])

init_POP_file = config[:MODEL_CORE][:cdata_file]
domain_file = config[:MODEL_CORE][:domain_file]
init_file  = config[:MODEL_MISC][:init_file]

Nz = length(config[:MODEL_CORE][:z_w]) - 1

Dataset(init_POP_file, "r") do ds
    global TEMP  = permutedims(nomissing(ds["TEMP"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
    global SALT  = permutedims(nomissing(ds["SALT"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

valid_idx = isfinite.(TEMP)

println("Create a model to save initial condition.")

ev = EMOM.Env(config[:MODEL_CORE])
mb = EMOM.ModelBlock(ev; init_core=false)

mb.fi.sv[:TEMP][valid_idx] .= TEMP[valid_idx]
mb.fi.sv[:SALT][valid_idx] .= SALT[valid_idx]

output_dir = dirname(init_file)
mkpath(dirname(init_file))

println(format("Output file: {}.", init_file))
EMOM.takeSnapshot(DateTimeNoLeap(1,1,1), mb, init_file)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module TestShow
    using DataFrames
    df = DataFrame(A = 1:3, B = ["x", "y", "z"])

    io = IOBuffer()
    show(io, df)
    show(io, df, true)
    showall(io, df)
    showall(io, df, true)

    subdf = df[df[:A] .> 1.0, :]
    show(io, subdf)
    show(io, subdf, true)
    showall(io, subdf)
    showall(io, subdf, true)

    dfr = DataFrameRow(df, 1)
    show(io, dfr)

    df = DataFrame(A = Array(UTF8String, 3))
end


text:
""
function constraint_thermal_limit_from_soft(pm::AbstractPowerModel, n::Int, f_idx, rate_a)
    l,i,j = f_idx
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_fr^2 + q_fr^2 <= (rate_a + sm_slack)^2)
end

""
function constraint_thermal_limit_to_soft(pm::AbstractPowerModel, n::Int, t_idx, rate_a)
    l,i,j = t_idx
    p_to = var(pm, n, :p, t_idx)
    q_to = var(pm, n, :q, t_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_to^2 + q_to^2 <= (rate_a + sm_slack)^2)
end

""
function constraint_gen_active_deviation(pm::AbstractPowerModel, n::Int, i, pg)
    pg_var = var(pm, n, :pg, i)
    pg_delta = var(pm, n, :pg_delta, i)

    JuMP.@constraint(pm.model, pg_var == pg + pg_delta)
end

""
function constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm::AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate + var(pm, :branch_cont_flow_vio, i))
end

""
function constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm::AbstractPowerModel, n::Int, i::Int, cut_map, rate)
    bus_injection = var(pm, :bus_pg)
    bus_withdrawal = var(pm, :bus_wdp)

    @constraint(pm.model, -sum(weight*(bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut_map) <= rate + var(pm, :branch_cont_flow_vio, i))
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

## original example, following here
move_to(cr, 16.0, 32.0);
curve_to(cr, 16.0, 16.0, 16.0, 16.0, 32.0, 16.0);

opath = Cairo.convert_cairo_path_data(Cairo.copy_path(cr));
dx,dy,ex,ey = path_extents(cr);

fpath = Cairo.convert_cairo_path_data(Cairo.copy_path_flat(cr));

stroke(cr);

restore(cr);

l = 2 # something like a line counter

# loop over all entries and just print out command and available points
for x in opath
    if x.element_type == Cairo.CAIRO_PATH_MOVE_TO
        s0 = "moveto"
    elseif x.element_type == Cairo.CAIRO_PATH_LINE_TO
        s0 = "lineto"
    elseif x.element_type == Cairo.CAIRO_PATH_CURVE_TO
        s0 = "curveto"
    elseif x.element_type == Cairo.CAIRO_PATH_CLOSE_PATH
        s0 = "closepath"
    end
    s1 = repr(x.points)
    move_to(cr,10.0,16.0+(14.0*l))
    l += 1
    show_text(cr,s0*s1); #
end

l = 2 # something like a line counter

# loop over all entries and just print out command and available points
for x in fpath
    if x.element_type == Cairo.CAIRO_PATH_MOVE_TO
        s0 = "moveto"
    elseif x.element_type == Cairo.CAIRO_PATH_LINE_TO
        s0 = "lineto"
    elseif x.element_type == Cairo.CAIRO_PATH_CURVE_TO
        s0 = "curveto"
    elseif x.element_type == Cairo.CAIRO_PATH_CLOSE_PATH
        s0 = "closepath"
    end
    s1 = repr(x.points)
    move_to(cr,10.0,50.0+(14.0*l))
    l += 1
    show_text(cr,s0*s1); #
end

## mark picture with current date
restore(cr);

move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));

write_to_png(c,"sample_copy_path_flat.png");


text:
module TestMultiplePackages

const SCRIPT = joinpath(@__DIR__, "../../samples/multiplepackages/runtests.jl")

test() = include(SCRIPT)

end  # module


text:
using LinearAlgebra: lu, I, ldiv!, rmul!, norm, qr, rank, nullspace
using SparseArrays: sparse, spzeros, SparseMatrixCSC, findnz
using Combinatorics: combinations
using MatrixNetworks: bfs, scomponents

export sparse_nullbasis, nullbasis

"""
    closest_permutation(mat::AbstractMatrix)

Given an `m × n` matrix `mat` where `m ≥ n`, find an assignment of rows to columns such
that each column is assigned a distinct row and the absolute values of the entries
corresponding to the assignment is as large as possible under lexicographical ordering.
Return a matrix of Bools where each column has a single true in it corresponding to its
assignment.
"""
function closest_permutation(mat::AbstractMatrix)
    nrows, ncols = size(mat)
    @assert nrows >= ncols
    ans = similar(BitArray, axes(mat)) # initialize answer as all false
    fill!(ans, false)

    ax2 = axes(mat, 2)

    # Do a first pass, ignoring collisions.
    colidxs = similar(ax2)
    for col in ax2
        @inbounds begin
            _, i = findmax(mat[:, col])
            ans[i, col] = true
            colidxs[col] = i
        end
    end

    # Prepare a BitArray `collisions` and mark collided columns in it.
    sp = sortperm(colidxs)
    sortedidxs = colidxs[sp]
    collisions = similar(BitArray, ax2)
    fill!(collisions, false)

    previndex = firstindex(sortedidxs)
    @inbounds for i in eachindex(sortedidxs)
        isequal(i, firstindex(sortedidxs)) && continue
        if isequal(sortedidxs[i], sortedidxs[previndex])
            collisions[sp[i]] = true
            collisions[sp[previndex]] = true
        end
        previndex = i
    end

    @inbounds stack = sp[collisions]
    isempty(stack) && return ans

    # reset columns in `ans` where collisions happened
    for col in stack
        @inbounds ans[:, col] .= false
    end

    # fall-back to single-threaded collision resolution on collided columns
    return _closest_permutation(mat, ans, stack)
end

function _closest_permutation(mat::AbstractMatrix, ans, stack)
    ax1 = axes(mat, 1)
    ax1l = length(ax1)
    possible_assignments = Vector{eltype(ax1)}() # row indices

    @inbounds while length(stack) > 0
        col = pop!(stack)
        resize!(possible_assignments, ax1l)
        possible_assignments .= axes(mat, 1)
        assigned = false
        while !assigned
            value, i = findmax(mat[possible_assignments, col])
            assignment = possible_assignments[i]
            collision = findfirst(ans[assignment, :])
            if collision == nothing
                ans[assignment, col] = true
                assigned = true
            elseif mat[assignment, collision] < value
                ans[assignment, col] = true
                assigned = true
                ans[assignment, collision] = false
                push!(stack, collision) # need to reassign collision
            else
                deleteat!(possible_assignments, i) # try somewhere else
            end
        end
    end
    return ans
end

"""
    match_vectors(port_vectors::AbstractMatrix{<:Number},
        eigenvectors::AbstractMatrix{<:Number})

For each column of `port_vectors` find a corresponding column of `eigenvectors` that has a
similar "shape". Match each `port_vector` to a different `eigenvector`.
"""
function match_vectors(port_vectors::AbstractMatrix{<:Number},
    eigenvectors::AbstractMatrix{<:Number})
    @assert size(port_vectors, 1) == size(eigenvectors, 1)
    @assert size(port_vectors, 2) <= size(eigenvectors, 2)
    mat = abs2.(eigenvectors' * port_vectors) # tall skinny matrix
    _, indices = findmax(closest_permutation(mat), dims=1)
    return [index[1] for index in vcat(indices...)]
end

"""
    inv_power_eigen(M::AbstractMatrix{<:Number};
        v0::AbstractVector{<:Number}=rand(Complex{real(eltype(M))}, size(M, 1)),
        λ0::Union{Number, Nothing}=nothing,
        pred::Function=(x,y) -> norm(x-y) < eps())

Use the inverse power method to find an eigenvalue and eigenvector.

# Arguments
- `v0::AbstractVector{<:Number}`: an initial value for the eigenvector. A random vector is
    used if `v0` is not given.
- `λ0::Union{Number, Nothing}`: an initial value for the eigenvalue. `v' * M * v` is used
    if `λ0` is not given.
- `pred::Function`: a function taking two successive eigenvalue estimates. The algorithm
    halts when `pred` returns true.
"""
function inv_power_eigen(M::AbstractMatrix{<:Number};
    v0::AbstractVector{<:Number}=rand(Complex{real(eltype(M))}, size(M, 1)),
    λ0::Union{Number, Nothing}=nothing,
    pred::Function=(x,y) -> norm(x-y) < eps())
    v = v0/norm(v0)
    l = length(v)
    @assert size(M) == (l, l)
    λ = (λ0 == nothing) ? v' * M * v : λ0
    done = false
    μ = λ
    A = lu(M - λ * I)
    while !done
        ldiv!(A, v)
        rmul!(v, 1/norm(v))
        μ = v' * M * v
        done = pred(λ, μ)::Bool
        λ = μ
    end
    return μ, v
end

#######################################################
# sparse nullbases
#######################################################

"""
    row_column_graph(mat::SparseMatrixCSC{<:Number, Int})

Let `mat` be a matrix with the following properties
    - for any two rows of `mat` there is at most one column on which both rows are nonzero
    - any column of `mat` has at most 2 nonzero values
The rows of `mat` form a simple graph where two rows have an edge if and only if they share
a nonzero column in `mat`. We add one additional vertex to this graph and for each column
with exactly one nonzero value, we add an edge between its nonzero row and the additional
vertex. In this graph every column with at least one nonzero value is represented by an
edge. Produce a matrix `A` where `A[i,j]` is the column shared by rows `i` and `j` if such a
column exists, and `0` otherwise. If `i == size(mat, 1) + 1`, then `A[i,j]` is chosen
arbitrarily from the columns whose unique nonzero value is in row `j`, or 0 if no such
columns exist. The case where `j == size(mat, 1) + 1` is analogous. `A` is an adjacency
matrix for the graph with extra information relating the graph to `mat`.
"""
function row_column_graph(mat::SparseMatrixCSC{<:Number, Int})
    num_vertices = size(mat, 1) + 1
    row_inds, col_inds, vals = Int[], Int[], Int[]
    matT = sparse(transpose(mat)) # SparseMatrixCSC row access is slow
    for (r1, r2) in combinations(1:(num_vertices-1), 2)
        inds = findall((matT[:, r1] .!= 0) .& (matT[:, r2] .!= 0))
        @assert length(inds) <= 1 # any two rows share at most 1 common nonzero column
        if length(inds) == 1
            c = inds[1]
            push!(row_inds, r1); push!(col_inds, r2); push!(vals, c)
            push!(row_inds, r2); push!(col_inds, r1); push!(vals, c)
        end
    end
    root = num_vertices
    connected_to_root = Set{Int}()
    for c in 1:size(mat, 2)
        inds = findall(mat[:, c] .!= 0)
        @assert length(inds) <= 2 # any column has at most 2 nonzero values
        if length(inds) == 1
            r = inds[1]
            if !(r in connected_to_root)
                push!(connected_to_root, r)
                push!(row_inds, root); push!(col_inds, r); push!(vals, c)
                push!(row_inds, r); push!(col_inds, root); push!(vals, c)
            end
        end
    end
    adj = sparse(row_inds, col_inds, vals, num_vertices, num_vertices)
    return adj
end

"""
    connected_row_column_graph(mat::SparseMatrixCSC{<:Number, Int})

Let `mat` be a matrix with the following properties
    - for any two rows of `mat` there is at most one column on which both rows are nonzero
    - any column of `mat` has at most 2 nonzero values
Produce a matrix `mat_connected` with the same nullspace as `mat` such that
`row_column_graph(mat_connected)` encodes a connected graph. Return `mat_connected` and
`row_column_graph(mat_connected)`.
"""
function connected_row_column_graph(mat::SparseMatrixCSC{<:Number, Int})
    adj = row_column_graph(mat)
    components = scomponents(adj)
    if components.number > 1 # adj is not connected
        root_component = components.map[end]
        remove_rows = Int[]
        zero_columns = Set{Int}()
        matT = sparse(transpose(mat)) # SparseMatrixCSC row access is slow
        for component in 1:components.number
            if component != root_component
                inds = findall(components.map .== component)
                num_dependent_rows = length(inds) - find_rank(matT[:, inds])
                if num_dependent_rows > 0 # some rows are redundant, remove them
                    append!(remove_rows, inds[1:num_dependent_rows])
                else # submat has full rank and trivial nullspace
                    append!(remove_rows, inds) # remove these rows
                    for i in inds # set columns in these rows to 0
                        push!(zero_columns, findall(matT[:, i] .!= 0)...)
                    end
                end
            end
        end
        num_zeros = length(zero_columns)
        zero_rows = sparse(1:num_zeros, collect(zero_columns), ones(num_zeros),
            num_zeros, size(mat, 2))
        remaining_inds = setdiff(1:size(mat, 1), remove_rows)
        mat = vcat(transpose(matT[:, remaining_inds]), zero_rows)
        adj = row_column_graph(mat) # adj is now connected
    end
    return mat, adj
end

find_rank(m::Matrix) = rank(m)
find_rank(m::SparseMatrixCSC) = rank(qr(1.0 * m))

"""
    bfs_tree(adj::AbstractMatrix{Int}, root::Int)

Produce the vertices and edges of the tree found by performing breadth first on the graph
with adjacency matrix `adj` starting from the vertex `root`. The vertices includes all
vertices of the tree besides the root.
"""
function bfs_tree(adj::AbstractMatrix{Int}, root::Int)
    dists, _, predecessors = bfs(adj, root)
    tree_vertices, tree_edges = Int[], Tuple{Int, Int}[]
    for i1 in reverse(sortperm(dists))
        i2 = predecessors[i1]
        if dists[i1] > 0 # exclude unreachables and the root
            push!(tree_vertices, i1)
            push!(tree_edges, tuple(sort([i1, i2])...))
        end
    end
    return tree_vertices, tree_edges
end

# Find the nullspace of mat by solving the equations corresponding to each of its rows. We
# solve row `tree_rows[i]`` by eliminating variable `tree_columns[i]`.
function sparse_nullbasis(mat::SparseMatrixCSC{<:Number, Int},
    tree_rows::AbstractVector{Int}, tree_columns::AbstractVector{Int})
    @assert length(tree_rows) == length(tree_columns) == size(mat, 1)
    non_tree_columns = setdiff(1:size(mat, 2), tree_columns)
    # first create mapping that preserves the non_tree_columns
    row_inds, col_inds, vals = Int[], Int[], Float64[]
    for (col_ind, row_ind) in enumerate(non_tree_columns)
        push!(row_inds, row_ind); push!(col_inds, col_ind); push!(vals, 1)
    end
    # work with transposes because SparseMatrixCSC row access is slow
    nullbasisT = sparse(col_inds, row_inds, vals, length(non_tree_columns), size(mat, 2))
    matT = sparse(transpose(mat))
    # now solve row constraints by walking the tree from the leaves to the root
    for (mat_row, mat_col) in zip(tree_rows, tree_columns)
        null_row = spzeros(size(nullbasisT, 1))
        for other_mat_col in findall(matT[:, mat_row] .!= 0)
            if other_mat_col != mat_col
                null_row += (-mat[mat_row, other_mat_col]/mat[mat_row, mat_col]) *
                    nullbasisT[:, other_mat_col]
            end
        end
        nullbasisT[:, mat_col] = null_row
    end
    return sparse(transpose(nullbasisT))
end

"""
    sparse_nullbasis(mat::SparseMatrixCSC{<:Number, Int})

Let `mat` be a matrix with the following properties
    - for any two rows of `mat` there is at most one column on which both rows are nonzero
    - any column of `mat` has at most 2 nonzero values
Produce a sparse matrix whose columns form a basis for the nullspace of `mat`.
"""
function sparse_nullbasis(mat::SparseMatrixCSC{<:Number, Int})
    mat, adj = connected_row_column_graph(mat)
    tree_rows, tree_edges = bfs_tree(adj, size(adj, 1))
    tree_columns = [adj[r1, r2] for (r1, r2) in tree_edges]
    return sparse_nullbasis(mat, tree_rows, tree_columns)
end

"""
    nullbasis(mat::AbstractMatrix{<:Number}; warn::Bool=true)

Return a `SparseMatrixCSC` whose columns form a basis for the null space of mat. If `mat`
satisfies the conditions of `sparse_nullbasis`, a sparse matrix is returned. Otherwise a
warning is issued and the `LinearAlgebra.nullspace` function is used.
"""
function nullbasis(mat::AbstractMatrix{<:Number}; warn::Bool=true)
    try
        return sparse_nullbasis(sparse(mat))
    catch AssertionError
        if warn
            @warn "sparse nullbasis not found"
        end
        return sparse(nullspace(collect(mat)))
    end
end


text:
using Drawing

p = 1.02
q = 0.9375

with(
     PNG("ndy.png", 300, 300),
#     TK(300, 300),
     Axes(negative=true)) do
         paint(Ink("darkgreen")) do
             square(100, align=5)
         end
         paint(Ink("red")) do
             square(2, align=5)
         end
         paint(Ink("black"), Font("Arial Black, Bold"; size=0.58)) do
             move(-q, p)
             text("NOT", align=1)
             move(0.0, 0.0)
             text("DEAD", align=5)
             move(q, -p)
             text("YET", align=9)
         end
     end


text:
using FactCheck
using Lint

max_stars_repo_name: jbnicolai/jlbox

max_issues_repo_name: jbnicolai/jlbox

max_forks_repo_name: jbnicolai/jlbox


text:
@doc "STANFORD UNIVERSITY INTERIM Propagation Model"
type SuiModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  terrainKind::Int
  shadowFading ::Float64
  checkFreqRange::Bool
  SuiModel() = new(2500,50,1.5,TerrainKind.A,8.2,true)
end

@doc "Calculate the path loss with the SUI Model"
function pathloss_op(model::SuiModel,dist::AbstractFloat)

  	d::Float64 = dist * 1000	#km to m
    f = model.freq
    if (model.checkFreqRange)
      if (f < 1900 || f > 11000)
    		error("The SUI model frequency range is  1.9-11GHz\n");
      end
    end
    txH::Float64 = model.txH
    rxH::Float64 = model.rxH
    # Terrain mode A
    mode = model.terrainKind
  	 a::Float64 = 4.6
  	 b::Float64 = 0.0075
  	 c::Float64 = 12.6
  	 s::Float64 = model.shadowFading  # shadow fading 8.2 dB and 10.6
  	 XhCF::Float64 = -10.8 #
     # Terrain mode B
  	if (mode == TerrainKind.B)
  		a = 4.0
  		b = 0.0065
  		c = 17.1
  	end
     # Terrain mode C
  	if (mode == TerrainKind.C)
  		a = 3.6
  		b = 0.005
  		c = 20
  			#// shadow fading 8.2 dB and 10.6c
  		XhCF = -20
  	end

  	d0::Float64 = 100;
  	A::Float64 = 20 * log10((4 * pi * d0) / (300 / f))
  	y::Float64 = (a - b * txH) + c / txH
  	Xf::Float64 = 6 * log10(f / 2000)
  	Xh::Float64 = XhCF * log10(rxH / 2)
  	return A + (10 * y * log10(d / d0)) + Xf + Xh + s
end


text:
#=
ir_dump.jl
2019-03-05 Jeff Fessler
=#

export ir_dump

"""
`ir_dump(x::Any ; io::IO = stdout)`

Show all the fields of a structure or `NamedTuple` more nicely than dump() does
"""
function ir_dump(y::Any ; io::IO = stdout)
	x = typeof(y)
	print(io, x)
	fields = fieldnames(x)
	fieldtypes = x.types
	for (idx,fd) in enumerate(fields)
		println(io)
		print(io, " ", fd, "::")
		ft = fieldtypes[idx]
		print(io, ft)
		if ft <: Number
			print(io, " ", getfield(y, fd))
		end
		if ft == String
			print(io, " '", getfield(y, fd), "'")
		end
	end
	println(io)
	nothing
end

"""
`ir_dump(:test; io::IO = IOBuffer())`

self test
"""
function ir_dump(test::Symbol ; io=IOBuffer())
	test != :test && throw(ArgumentError("test $test"))
	x = (a=1, b=2)
	ir_dump(x, io=io)
	true
end


text:
##################
# JacobianResult #
##################

type JacobianResult{V,J} <: ForwardDiffResult
    value::V
    jacobian::J
end

JacobianResult(x) = JacobianResult(similar(x), similar(x, length(x), length(x)))

value(result::JacobianResult) = result.value
jacobian(result::JacobianResult) = result.jacobian

###############
# API methods #
###############

function jacobian{N}(f, x, chunk::Chunk{N} = pickchunk(x);
                     multithread::Bool = false,
                     usecache::Bool = true)
    if N == length(x)
        return vector_mode_jacobian(f, x, chunk, usecache)
    elseif multithread
        return multithread_chunk_mode_jacobian(f, x, chunk, usecache)
    else
        return chunk_mode_jacobian(f, x, chunk, usecache)
    end
end

function jacobian{N}(f!, y, x, chunk::Chunk{N} = pickchunk(x);
                      multithread::Bool = false,
                      usecache::Bool = true)
    if N == length(x)
        return vector_mode_jacobian(f!, y, x, chunk, usecache)
    elseif multithread
        return multithread_chunk_mode_jacobian(f!, y, x, chunk, usecache)
    else
        return chunk_mode_jacobian(f!, y, x, chunk, usecache)
    end
end

function jacobian!{N}(out, f, x, chunk::Chunk{N} = pickchunk(x);
                      multithread::Bool = false,
                      usecache::Bool = true)
    if N == length(x)
        vector_mode_jacobian!(out, f, x, chunk, usecache)
    elseif multithread
        multithread_chunk_mode_jacobian!(out, f, x, chunk, usecache)
    else
        chunk_mode_jacobian!(out, f, x, chunk, usecache)
    end
    return out
end

function jacobian!{N}(out, f!, y, x, chunk::Chunk{N} = pickchunk(x);
                      multithread::Bool = false,
                      usecache::Bool = true)
    @assert !(isa(out, JacobianResult)) "use jacobian!(JacobianResult(out, y), f!, x, ...) instead of jacobian!(out, f!, y, x, ...)"
    if N == length(x)
        vector_mode_jacobian!(out, f!, y, x, chunk, usecache)
    elseif multithread
        multithread_chunk_mode_jacobian!(out, f!, y, x, chunk, usecache)
    else
        chunk_mode_jacobian!(out, f!, y, x, chunk, usecache)
    end
    return out
end

function jacobian!{N}(out::JacobianResult, f!, x, chunk::Chunk{N} = pickchunk(x); kwargs...)
    jacobian!(out.jacobian, f!, out.value, x, chunk; kwargs...)
    return out
end

#######################
# workhorse functions #
#######################

# result extraction #
#-------------------#

function load_jacobian_value!(y, ydual)
    for i in eachindex(y)
        y[i] = value(ydual[i])
    end
    return y
end

function load_jacobian_value!(out::JacobianResult, ydual)
    load_jacobian_value!(out.value, ydual)
    return out
end

function load_jacobian!(out, ydual)
    for col in 1:size(out, 2), row in 1:size(out, 1)
        out[row, col] = partials(ydual[row], col)
    end
    return out
end

@inline function load_jacobian!(out::JacobianResult, ydual)
    load_jacobian!(out.jacobian, ydual)
    return out
end

function load_jacobian_chunk!(out, ydual, index, chunksize)
    offset = index - 1
    for i in 1:chunksize
        col = i + offset
        for row in eachindex(ydual)
            out[row, col] = partials(ydual[row], i)
        end
    end
    return out
end

@inline function load_jacobian_chunk!(out::JacobianResult, ydual, index, chunksize)
    load_jacobian_chunk!(out.jacobian, ydual, index, chunksize)
    return out
end

# vector mode #
#-------------#

function compute_vector_mode_jacobian(f, x, chunk, usecache)
    cache = jacobian_cachefetch!(x, chunk, usecache)
    xdual = cache.duals
    seed!(xdual, x, cache.seeds, 1)
    return f(xdual)
end

function compute_vector_mode_jacobian(f!, y, x, chunk, usecache)
    cache = jacobian_cachefetch!(x, chunk, usecache)
    ycache = jacobian_cachefetch!(y, chunk, usecache, true)
    xdual = cache.duals
    ydual = ycache.duals
    seed!(xdual, x, cache.seeds, 1)
    seedall!(ydual, y, zero(eltype(ycache.seeds)))
    f!(ydual, xdual)
    return ydual
end

function vector_mode_jacobian{N}(f, x, chunk::Chunk{N}, usecache)
    ydual = compute_vector_mode_jacobian(f, x, chunk, usecache)
    out = similar(ydual, numtype(eltype(ydual)), length(ydual), N)
    return load_jacobian!(out, ydual)
end

function vector_mode_jacobian{N}(f!, y, x, chunk::Chunk{N}, usecache)
    ydual = compute_vector_mode_jacobian(f!, y, x, chunk, usecache)
    load_jacobian_value!(y, ydual)
    out = similar(y, length(y), N)
    return load_jacobian!(out, ydual)
end

function vector_mode_jacobian!(out, f, x, chunk, usecache)
    ydual = compute_vector_mode_jacobian(f, x, chunk, usecache)
    load_jacobian!(reshape(out, length(ydual), length(x)), ydual)
    return out
end

function vector_mode_jacobian!(out, f!, y, x, chunk, usecache)
    ydual = compute_vector_mode_jacobian(f!, y, x, chunk, usecache)
    load_jacobian_value!(y, ydual)
    load_jacobian!(reshape(out, length(y), length(x)), ydual)
    return out
end

# chunk mode #
#------------#

function jacobian_chunk_mode_expr(out_definition::Expr, cache_definition::Expr,
                                  ydual_compute::Expr, y_definition::Expr)
    return quote
        @assert length(x) >= N "chunk size cannot be greater than length(x) ($(N) > $(length(x)))"

        # precalculate loop bounds
        xlen = length(x)
        remainder = xlen % N
        lastchunksize = ifelse(remainder == 0, N, remainder)
        lastchunkindex = xlen - lastchunksize + 1
        middlechunks = 2:div(xlen - lastchunksize, N)

        # fetch and seed work arrays
        $(cache_definition)
        xdual = cache.duals
        seeds = cache.seeds
        zeroseed = zero(eltype(seeds))
        seedall!(xdual, x, zeroseed)

        # do first chunk manually to calculate output type
        seed!(xdual, x, seeds, 1)
        $(ydual_compute)
        seed!(xdual, x, zeroseed, 1)
        $(out_definition)
        out_reshaped = reshape(out, length(ydual), length(xdual))
        load_jacobian_chunk!(out_reshaped, ydual, 1, N)

        # do middle chunks
        for c in middlechunks
            i = ((c - 1) * N + 1)
            seed!(xdual, x, seeds, i)
            $(ydual_compute)
            seed!(xdual, x, zeroseed, i)
            load_jacobian_chunk!(out_reshaped, ydual, i, N)
        end

        # do final chunk
        seed!(xdual, x, seeds, lastchunkindex, lastchunksize)
        $(ydual_compute)
        load_jacobian_chunk!(out_reshaped, ydual, lastchunkindex, lastchunksize)

        $(y_definition)

        return reshape(out_reshaped, size(out))
    end
end

@eval function chunk_mode_jacobian{N}(f, x, chunk::Chunk{N}, usecache)
    $(jacobian_chunk_mode_expr(:(out = similar(x, numtype(eltype(ydual)), length(ydual), xlen)),
                               :(cache = jacobian_cachefetch!(x, chunk, usecache)),
                               :(ydual = f(xdual)),
                               :()))
end

@eval function chunk_mode_jacobian{N}(f!, y, x, chunk::Chunk{N}, usecache)
    $(jacobian_chunk_mode_expr(:(out = similar(y, numtype(eltype(ydual)), length(ydual), xlen)),
                               quote
                                   cache = jacobian_cachefetch!(x, chunk, usecache)
                                   ycache = jacobian_cachefetch!(y, chunk, usecache, true)
                                   ydual = ycache.duals
                                   yzeroseed = zero(eltype(ycache.seeds))
                               end,
                               :(f!(seedall!(ydual, y, yzeroseed), xdual)),
                               :(load_jacobian_value!(y, ydual))))
end

@eval function chunk_mode_jacobian!{N}(out, f, x, chunk::Chunk{N}, usecache)
    $(jacobian_chunk_mode_expr(:(),
                               :(cache = jacobian_cachefetch!(x, chunk, usecache)),
                               :(ydual = f(xdual)),
                               :()))
end

@eval function chunk_mode_jacobian!{N}(out, f!, y, x, chunk::Chunk{N}, usecache)
    $(jacobian_chunk_mode_expr(:(),
                               quote
                                   cache = jacobian_cachefetch!(x, chunk, usecache)
                                   ycache = jacobian_cachefetch!(y, chunk, usecache, true)
                                   ydual = ycache.duals
                                   yzeroseed = zero(eltype(ycache.seeds))
                               end,
                               :(f!(seedall!(ydual, y, yzeroseed), xdual)),
                               :(load_jacobian_value!(y, ydual))))
end