text:
using Observables2
using Test

@testset "Basics" begin
    xx = Observable([1, 2, 3])

    yy = observe!(xx, 2) do xs, factor
        xs .* factor
    end

    zz = observe!(yy) do y
        y ./ 10
    end

    xx[!] = [2, 3, 4]
    @test xx[] == [2, 3, 4]
    @test yy[] == [4, 6, 8]
    @test zz[] == [0.4, 0.6, 0.8]

    # check that all three observables are disabled through xx
    @test !is_disabled(xx)
    @test !is_disabled(yy)
    @test !is_disabled(zz)

    @test disable!(xx) == 3

    @test is_disabled(xx)
    @test is_disabled(yy)
    @test is_disabled(zz)
end

@testset "disabling middleman" begin

    objectcounter = Ref(0)

    mutable struct ObjectInMemory
        function ObjectInMemory()
            objectcounter[] += 1
            o = new()
            finalizer(o) do o
                objectcounter[] -= 1
                o
            end
            o
        end
    end

    qq = Observable([1, 2, 3])

    xx = observe!(identity, qq)

    yy = observe!(xx, 2) do xs, factor
        xs .* factor
    end

    @test objectcounter[] == 0

    zz = let
        memoryobject = ObjectInMemory()

        observe!(yy) do y
            # close over memoryobject
            if length(y) > 5
                println(memoryobject)
            end
            y ./ 10
        end
    end

    @test objectcounter[] == 1
    GC.gc()
    # the object shouldn't be garbage collected because it is referenced
    # in a closure in zz
    @test objectcounter[] == 1

    @test !is_disabled(qq)
    @test !is_disabled(xx)
    @test !is_disabled(yy)
    @test !is_disabled(zz)

    @test disable!(yy) == 2

    # yy and zz should be disabled, qq and xx shouldn't
    @test !is_disabled(qq)
    @test !is_disabled(xx)
    @test is_disabled(yy)
    @test is_disabled(zz)

    GC.gc()
    @test isnothing(zz.f)
    # this should work, but doesn't: @test objectcounter[] == 0
    # maybe that's due to garbage collector implementation, as I think the closure reference should be deleted correctly
end

@testset "Printing" begin
    xx = Observable([1, 2, 3])
    yy = observe!(identity, xx)

    @test string(xx) == """
        Observable{Array{Int64,1}} with 0 observable, 0 ordinary inputs, and 1 listeners.
        Value: [1, 2, 3]"""

    @test string(yy) == """
        Observable{Array{Int64,1}} with 1 observable, 0 ordinary inputs, and 0 listeners.
        Value: [1, 2, 3]"""

    stop_observing!(yy)

    @test string(xx) == """
        Observable{Array{Int64,1}} with 0 observable, 0 ordinary inputs, and 0 listeners.
        Value: [1, 2, 3]"""

    @test string(yy) == """
        Observable{Array{Int64,1}} with 0 observable, 1 ordinary inputs, and 0 listeners.
        Value: [1, 2, 3]"""
end

@testset "Typing" begin
    xx = Observable([1, 2, 3])
    @test typeof(xx) == Observable{Array{Int,1}}
    xx = Observable([1, 2, 3], type = Any)
    @test typeof(xx) == Observable{Any}
end

@testset "on onany off" begin
    x = Observable(1)
    y = Observable(2)

    testref = Ref(0)
    z = onany(x, y) do x, y
        testref[] = x + y
    end

    @test z isa ObservingFunction
    @test testref[] == 0
    x[] = 3
    @test testref[] == 5

    @test n_observable_inputs(z) == 2
    off(x, z)
    @test n_observable_inputs(z) == 1

    x[] = 4
    @test testref[] == 5
    y[] = 5
    @test testref[] == 8

end

@testset "ObservingFunction extras" begin
    x = Observable(1)

    r = Ref(0)

    # create an ObservingFunction
    y = on(x) do x
        r[] += 1
    end

    notify!(x)
    @test r[] == 1

    stop_observing!(y)
    notify!(x)
    # r should not have incremented further
    @test r[] == 1
end


text:
"""
    ParentLinks(::Type{T})
    ParentLinks(tree)

A trait which indicates whether a tree node stores references to its parents (`StoredParents()`) or
if the parents must be inferred from the tree structure (`ImplicitParents()`).

Trees for which `parentlinks` returns `StoredParents()` *MUST* implement [`parent`](@ref).

If `StoredParents()`, all nodes in the tree must also have `StoredParents()`, otherwise use
`ImplicitParents()`.

**OPTIONAL**: This should be implemented for a tree if parents of nodes are stored
```julia
AbstractTrees.parentlinks(::Type{<:TreeType}) = AbstractTrees.StoredParents()
parent(t::TreeType) = get_parent(t)
```
"""
abstract type ParentLinks end

"""
    StoredParents <: ParentLinks

Indicates that this node stores parent links explicitly. The implementation
is responsible for defining the parentind function to expose this
information.

If a node in a tree has this trait, so must all connected nodes.  If this is not the case,
use [`ImplicitParents`](@ref) instead.

## Required Methods
- [`parent`](@ref)
"""
struct StoredParents <: ParentLinks; end

"""
    ImplicitParents <: ParentLinks

Indicates that the tree does not store parents.  In these cases parents must be inferred from
the tree structure and cannot be inferred through a single node.
"""
struct ImplicitParents <: ParentLinks; end

ParentLinks(::Type) = ImplicitParents()
ParentLinks(tree) = ParentLinks(typeof(tree))

"""
    SiblingLinks(::Type{T})
    SiblingLinks(tree)

A trait which indicates whether a tree node stores references to its siblings (`StoredSiblings()`) or
must be inferred from the tree structure (`ImplicitSiblings()`).

If a node has the trait `StoredSiblings()`, so must all connected nodes in the tree.  Otherwise,
use `ImplicitSiblings()` instead.

**OPTIONAL**: This should be implemented for a tree if siblings of nodes are stored
```julia
AbstractTrees.SiblingLinks(::Type{<:TreeType}) = AbstractTrees.StoredSiblings()
```
"""
abstract type SiblingLinks end

"""
    StoredSiblings <: SiblingLinks

Indicates that this tree node stores sibling links explicitly, or can compute them
quickly (e.g. because the tree has a (small) fixed branching ratio, so the
current index of a node can be determined by quick linear search).

If a node has this trait, so must all connected nodes in the tree.  Otherwise, use `ImplicitSiblings()` instead.

## Required Methods
- [`nextsibling`](@ref)
"""
struct StoredSiblings <: SiblingLinks; end

"""
    ImplicitSiblings <: SiblingLinks

Indicates that tree nodes do not store references to siblings so that they must be inferred
from the tree structure.
"""
struct ImplicitSiblings <: SiblingLinks; end

SiblingLinks(::Type) = ImplicitSiblings()
SiblingLinks(tree) = SiblingLinks(typeof(tree))

"""
    ChildIndexing(::Type{N})
    ChildIndexing(node)

A trait indicating whether the tree node `n` has children (as returned by [`children`](@ref)) which can be
indexed using 1-based indexing.

If a node has the `IndexedChildren()` so must all connected nodes in the tree.  Otherwise, use
`NonIndexedChildren()` instead.

Options are either [`NonIndexedChildren`](@ref) (default) or [`IndexedChildren`](@ref).
"""
abstract type ChildIndexing end

"""
    IndexedChildren <: ChildIndexing

Indicates that the object returned by `children(n)` where `n` is a tree node is indexable (1-based).

If a node has this trait so must all connected nodes in the tree.  Otherwise, use `NonIndexedChildren()` instead.

## Required Methods
- A node `node` with this trait must return an indexable object from `children(node)`.
"""
struct IndexedChildren <: ChildIndexing end

"""
    NonIndexedChildren <: ChildIndexing

Indicates that the object returned by `children(n)` where `n` is a tree node is not necessarily indexable.
This trait applies to any tree which cannot guarantee indexable children in all cases, regardless of whether
the tree is indexable in special cases.  For example, `Array` has this trait even though there is a large
class of indexable trees consisting of arrays.
"""
struct NonIndexedChildren <: ChildIndexing end

ChildIndexing(::Type) = NonIndexedChildren()
ChildIndexing(node) = ChildIndexing(typeof(node))

"""
    childrentype(::Type{T})
    childrentype(n)

Indicates the type of the children (the *collection* of children, not individual children) of the tree node `n`
or its type `T`.  [`children`](@ref) should return an object of this type.

If the `childrentype` can be inferred from the type of the node alone, the type `::Type{T}` definition is preferred
(the latter will fall back to it).

**OPTIONAL**: In most cases, [`childtype`](@ref) is used instead.  If `childtype` is not defined it will fall back
to `eltype ∘ childrentype`.
"""
childrentype(nodetype::Type) = Base._return_type(children, Tuple{nodetype})
childrentype(node) = typeof(children(node))

"""
    childtype(::Type{T})
    childtype(n)

Indicates the type of children of the tree node `n` or its type `T`.

If `childtype` can be inferred from the type of the node alone, the type `::Type{T}` definition is preferred
(the latter will fall back to it).

**OPTIONAL**: It is strongly recommended to define this wherever possible, as without it almost no tree algorithms
can be type-stable.  If `childrentype` is defined and can be known from the node type alone, this function will
fall back to `eltype(childrentype(T))`.  If this gives a correct result it's not necessary to define `childtype`.
"""
childtype(::Type{T}) where {T} = eltype(childrentype(T))
childtype(node) = eltype(childrentype(node))

"""
    childstatetype(::Type{T})
    childstatetype(n)

Indicates the type of the iteration state of the tree node `n` or its type `T`.  This is used by tree traversal
algorithms which must retain this state.  It therefore is necessary to define this to ensure that most tree
traversal is type stable.

**OPTIONAL**: Type inference is used to attempt to
"""
childstatetype(::Type{T}) where {T} = Iterators.approx_iter_type(childrentype(T))
childstatetype(node) = childstatetype(typeof(node))

"""
    NodeType(::Type)
    NodeType(node)

A trait which specifiees whether a tree has a predictable node type (`HasNodeType()`) or not (`NodeTypeUnknown()`).

This is analogous to `Base.IteratorEltype`.  In particular the `IteratorEltype` of [`TreeIterator`](@ref) is dictated
by this trait.

The default value is `NodeTypeUnknown()`.
"""
abstract type NodeType end

"""
    HasNodeType <: NodeType

Indicates that this node is connected to a tree for which *all* nodes have types descended from `eltype(node)`.
"""
struct HasNodeType <: NodeType end

"""
    NodeTypeUnknown <: NodeType

Indicates that this node is connected to a tree for which it cannot be guaranteed that all nodes have the same
type.
"""
struct NodeTypeUnknown <: NodeType end

NodeType(::Type) = NodeTypeUnknown()
NodeType(node) = NodeType(typeof(node))

"""
    nodetype(::Type{T})
    nodetype(node) = nodetype(typeof(node))

Returns a type which must be a parent type of all nodes in the tree connected to `node`.  This can be used to,
for example, specify the `eltype` of any `TreeIterator` on `node`.
"""
nodetype(::Type) = Any
nodetype(node) = nodetype(typeof(node))

ext: jl


text:
# Cubic interpolation
a = 4.3
b = 1.7
c = 6.2
d = 5.9
cfunc = x -> a*x.^3 .+ b*x.^2 .+ c*x .+ d
cgrad = x -> 3a*x.^2 .+ 2b*x .+ c
chess = x -> 6a*x .+ 2b
xg = Float64[1:8;]
y = cfunc(xg)

# Low-level interface
yc = copy(y)
interp_invert!(yc,BCnan,InterpCubic)
ic = InterpGridCoefs(yc,InterpCubic)

x = 3.7
set_position(ic, BCnil, true, true, [x])
yi = interp(ic,yc)
@test_approx_eq_eps cfunc(x) yi 0.01*(abs(cfunc(x))+abs(yi))
set_gradient_coordinate(ic, 1)
gi = interp(ic, yc)
@test_approx_eq_eps cgrad(x) gi 0.1*(abs(cgrad(x))+abs(yi))
set_hessian_coordinate(ic, 1, 1)
hi = interp(ic, yc)
@test_approx_eq_eps chess(x) hi 0.1*(abs(chess(x))+abs(yi))

# High-level interface
yi = InterpGrid(y, BCnil, InterpCubic)
@test_approx_eq_eps yi[x] cfunc(x) 0.01*(abs(cfunc(x))+abs(yi[x]))
v,g = valgrad(yi,x)
@test_approx_eq_eps v cfunc(x) 0.01*(abs(cfunc(x))+abs(yi[x]))
@test_approx_eq_eps g cgrad(x) 0.01*(abs(cgrad(x))+abs(yi[x]))
v,g,h = valgradhess(yi,x)
@test_approx_eq_eps v cfunc(x) 0.01*(abs(cfunc(x))+abs(yi[x]))
@test_approx_eq_eps g cgrad(x) 0.01*(abs(cgrad(x))+abs(yi[x]))
@test_approx_eq_eps h chess(x) 0.01*(abs(chess(x))+abs(yi[x]))

# Test derivatives with other boundary conditions
func = x -> yi[x]
gnum = derivative_numer(func,x)
hnum = derivative2_numer(func,x)

for BC in (BCnan, BCna)
    yi = InterpGrid(y, BC, InterpCubic)
    @test_throws BoundsError yi[1.1,2.8]  # wrong dimensionality
    func = x -> yi[x]
    v,g = valgrad(yi,x)
    gnum = derivative_numer(func,x)
    @test_approx_eq_eps g gnum 10*cbrt(eps())*(abs(g)+abs(gnum))
    g = zeros(2)
    @test_throws ErrorException valgrad(g, yi, x)
    v,g,h = valgradhess(yi,x)
    gnum = derivative_numer(func,x)
    hnum = derivative2_numer(func,x)
    @test_approx_eq_eps g gnum 10*cbrt(eps())*(abs(g)+abs(gnum))
    @test_approx_eq_eps h hnum 10*cbrt(eps())*(abs(h)+abs(hnum))
    g = zeros(1)
    h = zeros(2,2)
    @test_throws ErrorException valgradhess(g,h,yi,x)
    g = zeros(2)
    h = zeros(1,1)
    @test_throws ErrorException valgradhess(g,h,yi,x)
end

# Test derivatives of 2D interpolation

y = rand(7,8)
yi = InterpGrid(y,BCnan,InterpCubic)
x = [3.4,6.1]
v,g = valgrad(yi,x...)
v2,g2,h2 = valgradhess(yi,x...)
func = (x1,x2) -> yi[x1,x2]

gnum = derivative_numer(func, tuple(x...), 1)
@test_approx_eq_eps g[1] gnum cbrt(eps())*(abs(g[1])+abs(gnum))
@test_approx_eq_eps g2[1] gnum cbrt(eps())*(abs(g2[1])+abs(gnum))

gnum = derivative_numer(func, tuple(x...), 2)
@test_approx_eq_eps g[2] gnum cbrt(eps())*(abs(g[2])+abs(gnum))
@test_approx_eq_eps g2[2] gnum cbrt(eps())*(abs(g2[2])+abs(gnum))

hnum = derivative2_numer(func, tuple(x...), (1, 1))
@test_approx_eq_eps h2[1,1] hnum 10*sqrt(sqrt(eps()))*(abs(h2[1,1])+abs(hnum))
hnum = derivative2_numer(func, tuple(x...), (1, 2))
@test_approx_eq_eps h2[1,2] hnum 10*sqrt(sqrt(eps()))*(abs(h2[1,2])+abs(hnum))
@test_approx_eq_eps h2[2,1] hnum 10*sqrt(sqrt(eps()))*(abs(h2[2,1])+abs(hnum))
hnum = derivative2_numer(func, tuple(x...), (2, 2))
@test_approx_eq_eps h2[2,2] hnum 10*sqrt(sqrt(eps()))*(abs(h2[2,2])+abs(hnum))


text:
module  TestPopData

using PopGenCore
using DataFrames
using Test

x = @nancycats

@testset "PopData" begin
    @testset "Struct types" begin
        @test x isa PopData
        @test x.metadata isa PopDataInfo
        @test x.genodata isa DataFrame
        @test x.sampleinfo isa DataFrame
        @test x.locusinfo isa DataFrame
    end

    @testset "PopDataInfo" begin
        @test x.metadata == x.info
        @test x.sampleinfo == x.metadata.sampleinfo
        @test PopData(x.genodata) isa PopData
        @test x.metadata.locusinfo == x.locusinfo
        @test x.metadata.ploidy == 2
        @test x.metadata.loci == 9
        @test x.metadata.samples == 237
        @test x.metadata.populations == 17
        @test x.metadata.biallelic == false
    end

    @testset "Indexing DataFrame" begin
        @test x[x.genodata.locus .== "fca8", :] isa DataFrame
        @test x[x.genodata.name .== "N100", :] isa DataFrame
        @test x[x.genodata.name .∈ Ref(["N100", "N217"]), :] isa DataFrame
        @test x[x.genodata.locus .∈ Ref(["fca8", "fca37"]), :] isa DataFrame
        @test x[x.genodata.locus .== "fca8", :name] isa AbstractArray
        @test x[(x.genodata.locus .== "fca8") .& (ishom(x.genodata.genotype, 139)), :name] isa AbstractArray

    end
    @testset "Indexing PopData" begin
        @test x[x.genodata.locus .== "fca8"] isa PopData
        @test x[x.genodata.name .== "N100"] isa PopData
        @test x[x.genodata.name .∈ Ref(["N100", "N217"])] isa PopData
        @test x[x.genodata.locus .∈ Ref(["fca8", "fca37"])] isa PopData
        @test x[(x.genodata.locus .== "fca8") .& (ishom(x.genodata.genotype, 139))] isa PopData
    end
end

end #module


text:
module CameraModels

using DocStringExtensions
using StaticArrays

# Exports
export AbstractCameraModel
export PinholeCamera
export RadialDistortion

# Abstract types
abstract type AbstractCameraModel end

"""
    $TYPEDEF

Pinhole Camera model is the most simplistic.

Notes
- https://en.wikipedia.org/wiki/Pinhole_camera
- Standard Julia *[Images.jl](https://juliaimages.org/latest/)-frame* convention is, `size(img) <==> (i,j) <==> (height,width) <==> (y,x)`,
  - Common *camera-frame* in computer vision and robotics, `(x,y) <==> (width,height) <==> (j,i)`,
  - Using top left corner of image as `(0,0)` in all cases.
  - Direct comparison with [OpenCV convention](https://docs.opencv.org/3.4/d9/d0c/group__calib3d.html) is:
    - `(x,y) [CamXYZ] <==> (j,i) [Images.jl] <==> (u,v) [OpenCV]` -- look very carfully at `(u,v) <==> (j,i)` in *image-frame*
- Always follow right hand rule for everything.

DevNotes
- https://en.wikipedia.org/wiki/Distortion_(optics)
"""
struct PinholeCamera{R <: Real} <: AbstractCameraModel
  K::SMatrix{3,3,R}
end

f_w(pc::PinholeCamera) = pc.K[1,1]
f_h(pc::PinholeCamera) = pc.K[2,2]
shear(pc::PinholeCamera) = pc.K[1,2]
c_w(pc::PinholeCamera) = pc.K[1,3]
c_h(pc::PinholeCamera) = pc.K[2,3]

set_f_w!(pc::PinholeCamera, val::Real) = (pc.K[1,1] = val)
set_f_h!(pc::PinholeCamera, val::Real) = (pc.K[2,2] = val)
set_shear!(pc::PinholeCamera, val::Real) = (pc.K[1,2] = val)
set_c_w!(pc::PinholeCamera, val::Real) = (pc.K[1,3] = val)
set_c_h!(pc::PinholeCamera, val::Real) = (pc.K[2,3] = val)

PinholeCamera(; f_w::Real=300,
                f_h::Real=f_w,
                c_w::R=320.0,
                c_h::Real=240,
                shear::Real=0) where {R <: Real} = SMatrix{3,3,R}([f_w shear c_w;
                                                                     0  f_h  c_h;
                                                                     0    0    1])
"""
    $TYPEDSIGNATURES

Constructor helper assuming you just have a camera image and need to start somewhere for a basic camera model.

Notes:
- Calibration will incorrect but hopefully in a distant ballpark to get the calibration process started.
- See [AprilTags.jl Calibration section](https://juliarobotics.org/AprilTags.jl/latest/#Camera-Calibration-1) for code and help.
"""
function PinholeCamera(img::AbstractArray{T,2}) where T
  f_w, c_w, c_h = size(img, 1), size(img, 2)/2, size(img, 1)/2
  f_h = f_w
  @info "Assuming default PinholeCamera from image $(size(img)):" f_w f_h c_w c_h
  PinholeCamera(f_w=f_w, f_h=f_h, c_w=c_w, c_h=c_h)
end

# Radial distortion model

"""
    $TYPEDEF

Slightly general Radial Distortion type, currently limited to StaticArrays.jl on CPU, but can later be extended to utilize GPUs -- see notes.

Notes
- Make sure `dest` image is large enough to encapsulate the resulting image after un-distortion

Example
```julia
using Images, FileIO, CameraModels

# load the image
img = load("myimg.jpg")

# genereate a radial distortion object
radialdistortion = RadialDistortion()
```

Reference:
From Wikipedia: https://en.wikipedia.org/wiki/Distortion_(optics)
  ( xd ,   yd ) = distorted image point as projected on image plane using specified lens,
  ( xu ,   yu ) = undistorted image point as projected by an ideal pinhole camera,
  ( xc ,   yc ) = distortion center,
              r = sqrt( (xd - xc)^2 + (yd - yc)^2 )

```math
xu = xc + (xd + xc) / (1 + K1*(r^2) + K2*(r^4) + ...)
yu = yc + (yd + yc) / (1 + K1*(r^2) + K2*(r^4) + ...)
```

DevNotes (Contributions welcome):
- TODO manage image clamping if `dest` is too small and data should be cropped out.
- TODO buffer radii matrix for better reuse on repeat image size sequences
- TODO dispatch with either CUDA.jl or AMDGPU.jl <:AbstractArray objects.
- TODO use Tullio.jl with multithreading and GPU
- TODO make sure LoopVectorization.jl tools like `@avx` is working
"""
struct RadialDistortion{N, R <: Real, K <: AbstractMatrix, C <: AbstractVector}
  Ki::K # SVector{N,R}
  center::C # SVector{2,R} # [h,w]
  # _radius2::Matrix{R} # perhaps SizedArray{R,2} or StaticArray{R,2} or GPUArray{R,2} depending on performance
end

# construction helper on CPU
# TODO better dispatch on Ki conversion to SVector
RadialDistortion(;Ki::AbstractVector{R}=SVector(0.0),
                  center::AbstractVector{R}=SVector{2,R}(0.0,0.0)) where {R <: Real} = RadialDistortion{length(Ki),R,SVector{length(Ki),R},SVector{2,R}}(Ki isa SVector ? Ki : (@SVector Ki), center)
#

function (rd::RadialDistortion{N,R,K,C})(dest::AbstractMatrix, src::AbstractMatrix) where {N,R,THR}
  # loop over entire image
  for h_d in size(src,1), w_d in size(src,2)
    # temporary coordinates
    @inbounds h_ = h_d - rd.center[1]
    @inbounds w_ = w_d - rd.center[2]
    # calculate the radius from distortion center
    _radius2 = h_^2 + w_^2
    # calculate the denominator
    _denomin = 1
    @inbounds @simd @fastmath for k in 1:N
      _denomin += rd.Ki[k]*(_radius2^k)
    end
    # calculate the new 'undistorted' coordinates and set equal to incoming image
    @inbounds @fastmath h_u = rd.center[1] + h_/_denomin
    @inbounds @fastmath w_u = rd.center[2] + h_/_denomin
    dest[h_u,w_u] = src[h_d,w_d]
  end
  nothing
end

end # module


text:
module TestParameterTypes

using Mimi
using Test

import Mimi:
    external_params, external_param, TimestepMatrix, TimestepVector,
    ArrayModelParameter, ScalarModelParameter, FixedTimestep, import_params!,
    set_first_last!, _get_param_times

#
# Test that parameter type mismatches are caught
#
expr = :(
    @defcomp BadComp1 begin
        a = Parameter(index=[time, regions], default=[10, 11, 12])  # should be 2D default
        function run_timestep(p, v, d, t)
        end
    end
)
@test_throws ErrorException eval(expr)

expr = :(
    @defcomp BadComp2 begin
        a = Parameter(default=[10, 11, 12])  # should be scalar default
        function run_timestep(p, v, d, t)
        end
    end
)
@test_throws ErrorException eval(expr)

#
# Test that the old type parameterization syntax errors
#
expr = :(
    @defcomp BadComp3 begin
        a::Int = Parameter()
        function run_timestep(p, v, d, t)
        end
    end
)
@test_throws LoadError eval(expr)

@defcomp MyComp begin
    a = Parameter(index=[time, regions], default=ones(101,3))
    b = Parameter(index=[time], default=1:101)
    c = Parameter(index=[regions])
    d = Parameter()
    e = Parameter(index=[four])
    f = Parameter{Array{Float64, 2}}()
    g = Parameter{Int}(default=10.0)    # value should be Int despite Float64 default
    h = Parameter(default=10)           # should be "numtype", despite Int default
    j = Parameter{Int}(index = [regions])

    function run_timestep(p, v, d, t)
    end
end

# Check that explicit number type for model works as expected
numtype = Float32
arrtype = Union{Missing, numtype}

m = Model(numtype)

set_dimension!(m, :time, 2000:2100)
set_dimension!(m, :regions, 3)
set_dimension!(m, :four, 4)

add_comp!(m, MyComp)
set_param!(m, :MyComp, :c, [4,5,6])
set_param!(m, :MyComp, :d, 0.5)   # 32-bit float constant
set_param!(m, :MyComp, :e, [1,2,3,4])
set_param!(m, :MyComp, :f, reshape(1:16, 4, 4))
set_param!(m, :MyComp, :j, [1,2,3])

Mimi.build!(m)    # applies defaults, creating external params in the model instance's copied definition
extpars = external_params(m.mi.md)

@test isa(extpars[:a], ArrayModelParameter)
@test isa(extpars[:b], ArrayModelParameter)
@test _get_param_times(extpars[:a]) == _get_param_times(extpars[:b]) == 2000:2100

@test isa(extpars[:c], ArrayModelParameter)
@test isa(extpars[:d], ScalarModelParameter)
@test isa(extpars[:e], ArrayModelParameter)
@test isa(extpars[:f], ScalarModelParameter) # note that :f is stored as a scalar parameter even though its values are an array

@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1, 2100}, arrtype, 1, Array{arrtype, 2}}
@test typeof(extpars[:b].values) == TimestepVector{FixedTimestep{2000, 1, 2100}, arrtype, Array{arrtype, 1}}

@test typeof(extpars[:c].values) == Array{arrtype, 1}
@test typeof(extpars[:d].value) == numtype
@test typeof(extpars[:e].values) == Array{arrtype, 1}
@test typeof(extpars[:f].value) == Array{Float64, 2}
@test typeof(extpars[:g].value) <: Int
@test typeof(extpars[:h].value) == numtype

# test updating parameters
@test_throws ErrorException update_param!(m, :a, 5) # expects an array
@test_throws ErrorException update_param!(m, :a, ones(101)) # wrong size
@test_throws ErrorException update_param!(m, :a, fill("hi", 101, 3)) # wrong type

set_param!(m, :a, Array{Int,2}(zeros(101, 3))) # should be able to convert from Int to Float
@test_throws ErrorException update_param!(m, :d, ones(5)) # wrong type; should be scalar
update_param!(m, :d, 5) # should work, will convert to float
new_extpars = external_params(m)    # Since there are changes since the last build, need to access the updated dictionary in the model definition
@test extpars[:d].value == 0.5      # The original dictionary still has the old value
@test new_extpars[:d].value == 5.   # The new dictionary has the updated value
@test_throws ErrorException update_param!(m, :e, 5) # wrong type; should be array
@test_throws ErrorException update_param!(m, :e, ones(10)) # wrong size
update_param!(m, :e, [4,5,6,7])

@test length(extpars) == 9          # The old dictionary has the default values that were added during build, so it has more entries
@test length(new_extpars) == 6
@test typeof(new_extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1, 2100}, arrtype, 1, Array{arrtype, 2}}

@test typeof(new_extpars[:d].value) == numtype
@test typeof(new_extpars[:e].values) == Array{arrtype, 1}

#------------------------------------------------------------------------------
# Test updating TimestepArrays with update_param!
#------------------------------------------------------------------------------

@defcomp MyComp2 begin
    x=Parameter(index=[time])
    y=Variable(index=[time])
    function run_timestep(p,v,d,t)
        v.y[t]=p.x[t]
    end
end

# 1. update_param! with Fixed Timesteps

m = Model()
set_dimension!(m, :time, 2000:2004)
add_comp!(m, MyComp2, first=2001, last=2003)
set_param!(m, :MyComp2, :x, [1, 2, 3, 4, 5])
# Year      x       Model   MyComp2
# 2000      1       first
# 2001      2               first
# 2002      3
# 2003      4               last
# 2004      5      last

update_param!(m, :x, [2.,3.,4.,5.,6.])
update_param!(m, :x, zeros(5))
update_param!(m, :x, [1,2,3,4,5])

set_dimension!(m, :time, 1999:2001)
# Year      x       Model   MyComp2
# 1999      missing first
# 2000      1
# 2001      2       last    first, last

x = external_param(m.md, :x)
@test ismissing(x.values.data[1])
@test x.values.data[2:3] == [1.0, 2.0]
@test _get_param_times(x) == 1999:2001
run(m) # should be runnable

update_param!(m, :x, [2, 3, 4]) # change x to match
# Year      x       Model   MyComp2
# 1999      2       first
# 2000      3
# 2001      4       last    first, last

x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.FixedTimestep{1999, 1, 2001}, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4.]
run(m)
@test ismissing(m[:MyComp2, :y][1])  # 1999
@test ismissing(m[:MyComp2, :y][2])  # 2000
@test m[:MyComp2, :y][3] == 4   # 2001

set_first_last!(m, :MyComp2, first = 1999, last = 2001)
# Year      x       Model   MyComp2
# 1999      2       first   first
# 2000      3
# 2001      4       last    last

run(m)
@test m[:MyComp2, :y] == [2, 3, 4]

# 2. Test with Variable Timesteps

m = Model()
set_dimension!(m, :time, [2000, 2005, 2020])
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])
# Year      x       Model   MyComp2
# 2000      1       first   first
# 2005      2
# 2010      3       last    last

set_dimension!(m, :time, [2000, 2005, 2020, 2100])
# Year      x       Model   MyComp2
# 2000      1       first   first
# 2005      2
# 2020      3               last
# 2100      missing last

x = external_param(m.md, :x)
@test ismissing(x.values.data[4])
@test x.values.data[1:3] == [1.0, 2.0, 3.0]

update_param!(m, :x, [2, 3, 4, 5]) # change x to match
# Year      x       Model   MyComp2
# 2000      2       first   first
# 2005      3
# 2020      4               last
# 2100      5        last

x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.VariableTimestep{(2000, 2005, 2020, 2100)}, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4., 5.]
run(m)
@test m[:MyComp2, :y][1] == 2   # 2000
@test m[:MyComp2, :y][2] == 3   # 2005
@test m[:MyComp2, :y][3] == 4   # 2020
@test ismissing(m[:MyComp2, :y][4]) # 2100 - past last attribute for component

set_first_last!(m, :MyComp2, first = 2000, last = 2020)
# Year      x       Model   MyComp2
# 2000      1       first   first
# 2005      2
# 2020      3       last    last

run(m)
@test m[:MyComp2, :y][1:3] == [2., 3., 4.]
@test ismissing(m[:MyComp2, :y][4])

# 3. Test updating from a dictionary

m = Model()
set_dimension!(m, :time, [2000, 2005, 2020])
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])

set_dimension!(m, :time, [2000, 2005, 2020, 2100])

update_params!(m, Dict(:x=>[2, 3, 4, 5]))
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.VariableTimestep{(2000, 2005, 2020, 2100)}, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4., 5.]
run(m)

@test m[:MyComp2, :y][1] == 2   # 2000
@test m[:MyComp2, :y][2] == 3   # 2005
@test m[:MyComp2, :y][3] == 4   # 2020
@test ismissing(m[:MyComp2, :y][4])   # 2100

# 4. Test updating the time index to a different length

m = Model()
set_dimension!(m, :time, 2000:2002)     # length 3
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])
# Year      x       Model   MyComp2
# 2000      1       first   first
# 2001      2
# 2002      3       last    last

set_dimension!(m, :time, 1999:2003)     # length 5
update_param!(m, :x, [2, 3, 4, 5, 6])
# Year      x       Model   MyComp2
# 1999      2       first
# 2000      3               first
# 2001      4
# 2002      5               last
# 2003      6       last

x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.FixedTimestep{1999, 1, 2003}, Union{Missing, Float64}, 1, 1}
@test x.values.data == [2., 3., 4., 5., 6.]

run(m)
@test ismissing(m[:MyComp2, :y][1])
@test m[:MyComp2, :y][2:4] == [3., 4., 5.]
@test ismissing(m[:MyComp2, :y][5])

set_first_last!(m, :MyComp2, first = 1999, last = 2001)
# Year      x       Model   MyComp2
# 1999      2       first   first
# 2000      3
# 2001      4               last
# 2002      5
# 2003      6       last

run(m)
@test ismissing(m[:MyComp2, :y][4])
@test ismissing(m[:MyComp2, :y][5])
@test m[:MyComp2, :y][1:3] == [2., 3., 4.]

# 5. Test all the warning and error cases

@defcomp MyComp3 begin
    regions=Index()
    x=Parameter(index=[time])       # One timestep array parameter
    y=Parameter(index=[regions])    # One non-timestep array parameter
    z=Parameter()                   # One scalar parameter
end

m = Model()                             # Build the model
set_dimension!(m, :time, 2000:2002)     # Set the time dimension
set_dimension!(m, :regions, [:A, :B])
add_comp!(m, MyComp3)
set_param!(m, :MyComp3, :x, [1, 2, 3])
set_param!(m, :MyComp3, :y, [10, 20])
set_param!(m, :MyComp3, :z, 0)

@test_throws ErrorException update_param!(m, :x, [1, 2, 3, 4]) # Will throw an error because size
update_param!(m, :y, [10, 15])
@test external_param(m.md, :y).values == [10., 15.]
update_param!(m, :z, 1)
@test external_param(m.md, :z).value == 1

# Reset the time dimensions
set_dimension!(m, :time, 1999:2001)

update_params!(m, Dict(:x=>[3,4,5], :y=>[10,20], :z=>0)) # Won't error when updating from a dictionary

@test external_param(m.md, :x).values isa Mimi.TimestepArray{Mimi.FixedTimestep{1999,1, 2001},Union{Missing,Float64},1}
@test external_param(m.md, :x).values.data == [3.,4.,5.]
@test external_param(m.md, :y).values == [10.,20.]
@test external_param(m.md, :z).value == 0

#------------------------------------------------------------------------------
# Test the three different set_param! methods for a Symbol type parameter
#------------------------------------------------------------------------------

@defcomp A begin
    p1 = Parameter{Symbol}()
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 10)
    add_comp!(m, A)
    return m
end

# Test the 3-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :p1, 3)  # Can't set it with an Int

set_param!(m, :p1, :foo)    # Set it with a Symbol
run(m)
@test m[:A, :p1] == :foo

# Test the 4-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :A, :p1, 3)

set_param!(m, :A, :p1, :foo)
run(m)
@test m[:A, :p1] == :foo

# Test the 5-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :A, :p1, :A_p1, 3)

set_param!(m, :A, :p1, :A_p1, :foo)
run(m)
@test m[:A, :p1] == :foo

#------------------------------------------------------------------------------
# Test that if set_param! errors in the connection step,
#       the created param doesn't remain in the model's list of params
#------------------------------------------------------------------------------

@defcomp A begin
    p1 = Parameter(index = [time])
end

@defcomp B begin
    p1 = Parameter(index = [time])
end

m = Model()
set_dimension!(m, :time, 10)
add_comp!(m, A)
add_comp!(m, B)

@test_throws ErrorException set_param!(m, :p1, 1:5)     # this will error because the provided data is the wrong size
@test isempty(m.md.external_params)                     # But it should not be added to the model's dictionary

end #module


text:
include("opfdata_structjump.jl")
using StructJuMP, JuMP
using StructJuMPSolverInterface

type SCOPFData
  raw::RawData
  lines_off::Array
  #Float64::gener_ramp #generator ramp limit for contingency (percentage)
end

function scopf_solve(scopfmodel, scopfdata, solver, prof::Bool)

  status = solve(scopfmodel;solver=solver, with_prof=prof)

  # getVarValue(scopfmodel)
  return scopfmodel,status
end

function scopf_model(raw::RawData)
  opfdata = opf_loaddata(raw) #load root node
  lines_off=Array{Line}(length(ARGS)-1)
  for l in 1:length(lines_off)
    lines_off[l] = opfdata.lines[parse(Int,ARGS[l+1])]
  end
  sd = SCOPFData(raw,lines_off)

  nscen = length(sd.lines_off)

  opfmodel = StructuredModel(num_scenarios=nscen)

  #shortcuts for compactness
  lines = opfdata.lines; buses = opfdata.buses; generators = opfdata.generators; baseMVA = opfdata.baseMVA;
  busIdx = opfdata.BusIdx; FromLines = opfdata.FromLines; ToLines = opfdata.ToLines; BusGeners = opfdata.BusGenerators;
  nbus  = length(buses); nline = length(lines); ngen  = length(generators);

  #branch admitances
  YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(lines, buses, baseMVA)

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax)
  @variable(opfmodel, 0<=extra[i=1:ngen]<=0.05*generators[i].Pmax)
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax)

  @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax)
  @variable(opfmodel,Va[1:nbus])
  #fix the voltage angle at the reference bus
  setlowerbound(Va[opfdata.bus_ref], buses[opfdata.bus_ref].Va)
  setupperbound(Va[opfdata.bus_ref], buses[opfdata.bus_ref].Va)
  # setlowerbound(Va[opfdata.bus_ref], -.1)
  # setupperbound(Va[opfdata.bus_ref], .1)

  #objective function
  @NLobjective(opfmodel, Min, (1/(nscen+1))*sum( generators[i].coeff[generators[i].n-2]*(baseMVA*(Pg[i] + extra[i]))^2
                         +generators[i].coeff[generators[i].n-1]*(baseMVA*(Pg[i]+extra[i]))
                     +generators[i].coeff[generators[i].n  ] for i=1:ngen))

  # @show "Root - Buses: %d  Lines: %d  Generators", nbus, nline, ngen

  # power flow balance
  for b in 1:nbus
    #real part
    @NLconstraint(
      opfmodel,
      ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b]^2
      + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
      + sum( Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
      - ( sum(baseMVA*(Pg[g]+extra[g]) for g in BusGeners[b]) - buses[b].Pd ) / baseMVA      # Sbus part
      ==0)

    #imaginary part
    @NLconstraint(
      opfmodel,
      ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b]^2
      + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
      + sum( Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
      - ( sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd ) / baseMVA      #Sbus part
      ==0)
  end

  # generator max
  @NLconstraint(opfmodel, ex[i=1:ngen], generators[i].Pmin <= Pg[i] + extra[i] <= generators[i].Pmax)

  #
  # branch/lines flow limits
  #
  nlinelim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nlinelim += 1
      flowmax=(lines[l].rateA/baseMVA)^2

      #branch apparent power limits (from bus)
      Yff_abs2=YffR[l]^2+YffI[l]^2; Yft_abs2=YftR[l]^2+YftI[l]^2
      Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
      @NLconstraint(opfmodel,
        Vm[busIdx[lines[l].from]]^2 *
        ( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
          + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
        )
        - flowmax <=0)

      #branch apparent power limits (to bus)
      Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
      Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
      @NLconstraint(
        opfmodel,
        Vm[busIdx[lines[l].to]]^2 *
        ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
          + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
        )
        - flowmax <=0)
    end
  end

  # @show "Lines with limits  ", nlinelim

  for c in getLocalChildrenIds(opfmodel)
  # @declare_second_stage(opfmodel, c, begin
    opfmodel_c = StructuredModel(parent=opfmodel,id=c)
    opfdata_c=opf_loaddata(raw, sd.lines_off[c])

    #shortcuts for compactness
    lines = opfdata_c.lines
    buses = opfdata_c.buses
    generators = opfdata_c.generators
    baseMVA = opfdata_c.baseMVA

    busIdx = opfdata_c.BusIdx
    FromLines = opfdata_c.FromLines
    ToLines = opfdata_c.ToLines
    BusGeners = opfdata_c.BusGenerators

    nbus  = length(buses)
    nline = length(lines)
    ngen  = length(generators)

    # @printf("Contingency %d -> Buses: %d  Lines: %d  Generators: %d\n", c, nbus, nline, ngen)

    YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(opfdata_c.lines, opfdata_c.buses, opfdata_c.baseMVA)

    @variable(opfmodel_c, 0<=extra[i=1:ngen]<=0.05*generators[i].Pmax)
    @variable(opfmodel_c, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax)
    @variable(opfmodel_c, Va[1:nbus])

    #fix the voltage angle at the reference bus
    setlowerbound(Va[opfdata_c.bus_ref], buses[opfdata_c.bus_ref].Va)
    setupperbound(Va[opfdata_c.bus_ref], buses[opfdata_c.bus_ref].Va)

    @NLobjective(opfmodel_c, Min, (1/(nscen+1))*sum( generators[i].coeff[generators[i].n-2]*(baseMVA*(Pg[i] + extra[i]))^2
                         +generators[i].coeff[generators[i].n-1]*(baseMVA*(Pg[i]+extra[i]))
                     +generators[i].coeff[generators[i].n  ] for i=1:ngen))

    # power flow balance
    for b in 1:nbus
      #real part
      @NLconstraint(
        opfmodel_c,
        ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b]^2
        + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
        + sum( Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
        - ( sum(baseMVA*(Pg[g] + extra[g]) for g in BusGeners[b]) - buses[b].Pd ) / baseMVA      # Sbus part
        ==0)

      #imaginary part
      @NLconstraint(
        opfmodel_c,
        ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b]^2
        + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
        + sum( Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
        - ( sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd ) / baseMVA      #Sbus part
        ==0)
    end

    #generator max
    @NLconstraint(opfmodel_c, ex[i=1:ngen], generators[i].Pmin <= Pg[i] + extra[i] <= generators[i].Pmax)

    #
    # branch/lines flow limits
    #
    nlinelim=0
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        nlinelim += 1
        flowmax=(lines[l].rateA/baseMVA)^2

        #branch apparent power limits (from bus)
        Yff_abs2=YffR[l]^2+YffI[l]^2; Yft_abs2=YftR[l]^2+YftI[l]^2
        Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
        @NLconstraint(
          opfmodel_c,
          Vm[busIdx[lines[l].from]]^2 *
          ( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
            + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
          )
          - flowmax <=0)

        #branch apparent power limits (to bus)
        Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
        Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
        @NLconstraint(
          opfmodel_c,
          Vm[busIdx[lines[l].to]]^2 *
          ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
            + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]*(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])-Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]]))
          )
          - flowmax <=0)
      end
    end

    # @printf("Contingency %d -> Buses: %d  Lines: %d  Generators: %d\n", c, nbus, nline, ngen)
    # println("     lines with limits:  ", nlinelim)
  end
  # ) #end second stage
  # @show "acopf_model_sc  - done"
  return opfmodel, sd
end

function scopf_structjump_outputAll(opfmodel, scopf_data)
  #shortcuts for compactness
  sd=scopf_data; opf_data = opf_loaddata(sd.raw);
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", getObjectiveVal(opfmodel), "USD/hr")
  VM=getvalue(getindex(opfmodel,:Vm)); VA=getvalue(getindex(opfmodel,:Va))
  PG=getvalue(getindex(opfmodel,:Pg)); QG=getvalue(getindex(opfmodel,:Qg))

  EX=getvalue(getindex(opfmodel,:extra));

  # printing the first stage variables
  println("============================= BUSES ==================================")
  println("  Generator  |  extra ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")
  println("----------------------------------------------------------------------")
  for i in 1:ngen
      @printf("  %10d | %6.2f \n",generators[i].bus, EX[i])
  end
  println("\n")

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
        buses[i].bus_i, VM[i], VA[i]*180/pi,
        length(BusGeners[i])==0?"   --          --  ":@sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end
  return

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MathProgBase.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MathProgBase.eval_g(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "\%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "\%" )
          #@printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  return
end

function scopf_init_x(scopfmodel,scopfdata)
  raw = scopfdata.raw
  lines_off = scopfdata.lines_off
  for i in getLocalBlocksIds(scopfmodel)
    if(i==0)
      opfdata = opf_loaddata(raw)
      Pg0,Qg0,Vm0,Va0 = scopf_compute_x0(opfdata)
      extra0 = 0.025*Pg0
      setvalue(getindex(scopfmodel, :Pg), Pg0)
      setvalue(getindex(scopfmodel, :extra), extra0)
      setvalue(getindex(scopfmodel, :Qg), Qg0)
      setvalue(getindex(scopfmodel, :Vm), Vm0)
      setvalue(getindex(scopfmodel, :Va), Va0)
    else
      mm = getchildren(scopfmodel)[i]
      opfdata_c=opf_loaddata(raw,lines_off[i])
      Pg0,Qg0,Vm0,Va0 = scopf_compute_x0(opfdata_c)
      extra0 = 0.025*Pg0
      setvalue(getindex(mm, :extra), extra0)
      setvalue(getindex(mm, :Vm), Vm0)
      setvalue(getindex(mm, :Va), Va0)
    end
  end
end

# Compute initial point for IPOPT based on the values provided in the case data
function scopf_compute_x0(opfdata)
  Pg=zeros(length(opfdata.generators));
  Qg=zeros(length(opfdata.generators));
  i=1
  for g in opfdata.generators
    # set the power levels in in between the bounds as suggested by matpower
    # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
    Pg[i]=0.5*(g.Pmax+g.Pmin)
    Qg[i]=0.5*(g.Qmax+g.Qmin)
    i=i+1
  end
  @assert i-1==length(opfdata.generators)

  Vm=zeros(length(opfdata.buses)); i=1;
  for b in opfdata.buses
    # set the ini val for voltage magnitude in between the bounds
    # (case data contains initials values in Vm - not used with IPOPT)
    Vm[i]=0.5*(b.Vmax+b.Vmin);
    i=i+1
  end
  @assert i-1==length(opfdata.buses)

  # set all angles to the angle of the reference bus
  Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

  return Pg,Qg,Vm,Va
end


text:
# These 3 methods are defined on concrete formulations of the branches to avoid ambiguity
construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::DeviceModel{<:PSY.ACBranch, StaticBranch},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchBounds},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::DeviceModel{<:PSY.DCBranch, <:AbstractDCLineFormulation},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::Type{<:PM.AbstractPowerModel},
) = nothing

# For DC Power only. Implements constraints
function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranch},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end

# For DC Power only
function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranch},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end

    add_variables!(optimization_container, S, devices, StaticBranch())
    branch_flow_values!(optimization_container, devices, model, S)
    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranchBounds},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end

    add_variables!(optimization_container, S, devices, StaticBranchBounds())
    branch_flow_values!(optimization_container, devices, model, S)
    branch_rate_bounds!(optimization_container, devices, model, S)
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranchUnbounded},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end

    add_variables!(optimization_container, S, devices, StaticBranchUnbounded())
    branch_flow_values!(optimization_container, devices, model, S)
    return
end

# For AC Power only. Implements Bounds on the active power and rating constraints on the aparent power
function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranch},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_bounds!(optimization_container, devices, model, S)
    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, StaticBranchBounds},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_bounds!(optimization_container, devices, model, S)
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, <:AbstractDCLineFormulation},
    ::Type{S},
) where {B <: PSY.DCBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, U},
    ::Type{S},
) where {
    B <: PSY.DCBranch,
    U <: Union{HVDCLossless, HVDCUnbounded},
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end

    add_variables!(optimization_container, FlowActivePowerVariable, devices, U())
    add_variable_to_expression!(optimization_container, devices, model, S)
    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{B, U},
    ::Type{S},
) where {
    B <: PSY.DCBranch,
    U <: AbstractDCLineFormulation,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end

    add_variables!(optimization_container, FlowActivePowerVariable, devices, U())

    branch_rate_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    return
end


text:
export codemirror

immutable CodeMirror <: Widget
    code::String
    language::String
    theme::String
    linenumbers::Bool
    tabsize::Int
end

codemirror(;
            name=:_code,
            code="",
            language="julia",
            theme="elegant",
            linenumbers=true,
            tabsize=4) =
    CodeMirror(code, language, theme, linenumbers, tabsize) |>
        c -> hasstate(c, name=name, attr="currentValue", trigger="change")

# Render to virtual DOM
render(c::CodeMirror) =
    Elem("code-mirror",
        value=c.code,
        mode=c.language,
        theme=c.theme,
        lineNumbers=c.linenumbers)


text:
# ---------------------------------------------------------------------------
# Violin Plot

const _violin_warned = [false]

function violin_coords(y; trim::Bool=false)
    kd = KernelDensity.kde(y, npoints = 200)
    if trim
        xmin, xmax = Plots.ignorenan_extrema(y)
        inside = Bool[ xmin <= x <= xmax for x in kd.x]
        return(kd.density[inside], kd.x[inside])
    end
    kd.density, kd.x
end

@recipe function f(::Type{Val{:violin}}, x, y, z; trim=true, side=:both)
    xsegs, ysegs = Segments(), Segments()
    glabels = sort(collect(unique(x)))
    bw = d[:bar_width]
    bw == nothing && (bw = 0.8)
    for (i,glabel) in enumerate(glabels)
        widths, centers = violin_coords(y[filter(i -> _cycle(x,i) == glabel, 1:length(y))], trim=trim)
        isempty(widths) && continue

        # normalize
        hw = 0.5_cycle(bw, i)
        widths = hw * widths / Plots.ignorenan_maximum(widths)

        # make the violin
        xcenter = Plots.discrete_value!(d[:subplot][:xaxis], glabel)[1]
        if (side==:right)
          xcoords = vcat(widths, zeros(length(widths))) + xcenter
        elseif (side==:left)
          xcoords = vcat(zeros(length(widths)), -reverse(widths)) + xcenter
        else
          xcoords = vcat(widths, -reverse(widths)) + xcenter
        end
        ycoords = vcat(centers, reverse(centers))

        push!(xsegs, xcoords)
        push!(ysegs, ycoords)
    end

    seriestype := :shape
    x := xsegs.pts
    y := ysegs.pts
    ()
end
Plots.@deps violin shape


text:
module FunctionBarrier

using Base.Meta

export @barrier

struct Var
    name::Symbol
    num_esc::Int
end

function makelet(v::Var)
    ex = v.name
    for i=1:v.num_esc
        ex = esc(ex)
    end
    return ex
end

# Wrap `barrier_expression` in a `function` block to improve efficiency.
function wrap_barrier(module_, barrier_expression)
    ex = macroexpand(module_, barrier_expression)
    bound_vars, captured_vars = [Var(:end, 0), Var(:begin, 0)], Var[]
    find_var_uses!(captured_vars, bound_vars, ex, 0)
    fname = gensym()
    lastex = ex.args[end]
    args = map(makelet, captured_vars)
    if lastex isa Symbol || last isa Expr &&
        lastex.head == :tuple &&
        all(a -> a isa Symbol, lastex.args)
        retvar = lastex
    else
        retvar = gensym()
    end
    quote
        function $fname($(args...))
            $barrier_expression
        end
        $retvar = $fname($(args...))
    end
end

macro barrier(ex)
    esc(wrap_barrier(__module__, ex))
end

function find_var_uses!(varlist, bound_vars, ex, num_esc)
    if isa(ex, Symbol)
        var = Var(ex,num_esc)
        if !(var in bound_vars)
            var ∈ varlist || push!(varlist, var)
        end
        return varlist
    elseif isa(ex, Expr)
        if ex.head == :quote || ex.head == :line || ex.head == :inbounds
            return varlist
        end
        if ex.head == :(=)
            find_var_uses_lhs!(varlist, bound_vars, ex.args[1], num_esc)
            find_var_uses!(varlist, bound_vars, ex.args[2], num_esc)
        elseif ex.head == :kw
            find_var_uses!(varlist, bound_vars, ex.args[2], num_esc)
        elseif ex.head == :for || ex.head == :while || ex.head == :comprehension || ex.head == :let
            # New scopes
            inner_bindings = copy(bound_vars)
            find_var_uses!(varlist, inner_bindings, ex.args, num_esc)
        elseif ex.head == :try
            # New scope + ex.args[2] is a new binding
            find_var_uses!(varlist, copy(bound_vars), ex.args[1], num_esc)
            catch_bindings = copy(bound_vars)
            !isa(ex.args[2], Symbol) || push!(catch_bindings, Var(ex.args[2],num_esc))
            find_var_uses!(varlist,catch_bindings,ex.args[3], num_esc)
            if length(ex.args) > 3
                finally_bindings = copy(bound_vars)
                find_var_uses!(varlist,finally_bindings,ex.args[4], num_esc)
            end
        elseif ex.head == :call
            find_var_uses!(varlist, bound_vars, ex.args[2:end], num_esc)
        elseif ex.head == :local
            foreach(ex.args) do e
                if !isa(e, Symbol)
                    find_var_uses!(varlist, bound_vars, e, num_esc)
                end
            end
        elseif ex.head == :(::)
            find_var_uses_lhs!(varlist, bound_vars, ex, num_esc)
        elseif ex.head == :escape
            find_var_uses!(varlist, bound_vars, ex.args[1], num_esc+1)
        else
            find_var_uses!(varlist, bound_vars, ex.args, num_esc)
        end
    end
    varlist
end

find_var_uses!(varlist, bound_vars, exs::Vector, num_esc) =
    foreach(e->find_var_uses!(varlist, bound_vars, e, num_esc), exs)

function find_var_uses_lhs!(varlist, bound_vars, ex, num_esc)
    if isa(ex, Symbol)
        var = Var(ex,num_esc)
        var ∈ bound_vars || push!(bound_vars, var)
    elseif isa(ex, Expr)
        if ex.head == :tuple
            find_var_uses_lhs!(varlist, bound_vars, ex.args, num_esc)
        elseif ex.head == :(::)
            find_var_uses!(varlist, bound_vars, ex.args[2], num_esc)
            find_var_uses_lhs!(varlist, bound_vars, ex.args[1], num_esc)
        else
            find_var_uses!(varlist, bound_vars, ex.args, num_esc)
        end
    end
end

find_var_uses_lhs!(varlist, bound_vars, exs::Vector, num_esc) = foreach(e->find_var_uses_lhs!(varlist, bound_vars, e, num_esc), exs)

end


text:
#!/usr/bin/env julia

using Luxor, Colors, Test, Random
Random.seed!(42)

Drawing(800, 800, "morepolytests.png")
origin()
background("white")

@testset "collinear points tests" begin
    pgon = box(BoundingBox(), vertices=true)
    @test length(pgon) == 4

    pgon = polysample(pgon, 12)
    @test length(pgon) == 12

    pgon1 = polyremovecollinearpoints(pgon)
    @test length(pgon1) == 4

    pgon = [Point(x, 0) for x in -250:1.4:250]
    pgon1 = polyremovecollinearpoints(pgon)
    # they've all gone!
    @test length(pgon1) == 0
    pgon = [Point(0, y) for y in -250:1.4:250]
    pgon1 = polyremovecollinearpoints(pgon)
    # they've all gone!
    @test length(pgon1) == 0

    badpgon = Point[O, O + (10, 0)]
    l = length(polyremovecollinearpoints(badpgon))
    @test l == 0

    badpgon = Point[O]
    l = length(polyremovecollinearpoints(badpgon))
    @test l == 0

end

@testset "clockwise tests" begin
    pgon = box(BoundingBox(), vertices=true)
    @test ispolyclockwise(pgon) == true
    pgon1 = reverse(pgon)
    @test ispolyclockwise(pgon1) == false

    pgon = star(O, 100, 21, 0.1, vertices=true)
    @test ispolyclockwise(pgon) == true
    pgon1 = reverse(pgon)
    @test ispolyclockwise(pgon1) == false

    pgon = [polar(100, θ) for θ in 0:π/12:2π]
    @test ispolyclockwise(pgon) == true
    pgon1 = reverse(pgon)
    @test ispolyclockwise(pgon1) == false

    pgon = [Point(x, 0) for x in -250:20:250]
    # move last point lower, makes it clockwise
    pgon[end] += (0, 50)
    @test ispolyclockwise(pgon) == true
    # move it back
    pgon[end] -= (0, 50)
    #  what should a straight line really be, clockwise or not?
    # looks like it's not clockwise
    @test ispolyclockwise(pgon) == false
    @test ispolyclockwise([O, O + (10, 0)]) == false
end

@testset "polyorientation tests" begin
    # clockwise
    orientations = Float64[]
    for i in 1.0:-0.1:0.1
        pgon = box(BoundingBox() * i, vertices=true)
        # polyarea and polyorientation are basically the same, but
        # the latter is doubled and carries sign
        @test polyarea(pgon) == polyorientation(pgon)/2
        push!(orientations, polyorientation(pgon))
    end
    # each polyorientation is area, and should all be clockwise...
    @test all(isequal(25600.00), diff(diff(orientations)))
    # not clockwise
    orientations = Float64[]
    for i in 1.0:-0.1:0.1
        pgon = reverse(box(BoundingBox() * i, vertices=true))
        @test polyarea(pgon) == -polyorientation(pgon)/2
        push!(orientations, polyorientation(pgon))
    end
    @test all(isequal(-25600.00), diff(diff(orientations)))

    # I don't like these, but we've never defined what a polygon _should_ be:
    badpgon = Point[O, O + (10, 0)]
    @test iszero(polyorientation(badpgon))
    @test !ispolyclockwise(badpgon)

    badpgon = Point[O, O + (0, 10)]
    @test iszero(polyorientation(badpgon))
    @test !ispolyclockwise(badpgon)

    badpgon = Point[O, O + (0, 0)]
    @test iszero(polyorientation(badpgon))
    @test !ispolyclockwise(badpgon)

    badpgon = Point[O, O]
    @test iszero(polyorientation(badpgon))
    @test !ispolyclockwise(badpgon)

    badpgon = Point[O]
    @test iszero(polyorientation(badpgon))
    @test !ispolyclockwise(badpgon)
end

@testset "point inside triangle tests" begin
    goodtriangle = Point[O, Point(100, 100), Point(-100, 100)]
    # ispointinsidetriangle() accepts points on the borders
    @test ispointinsidetriangle(O, goodtriangle)
    @test ispointinsidetriangle(Point(100, 100), goodtriangle)
    @test ispointinsidetriangle(Point(50, 100), goodtriangle)
    @test ispointinsidetriangle(Point(0, 100), goodtriangle)
    @test ispointinsidetriangle(Point(-100, 100), goodtriangle)
    @test ispointinsidetriangle(Point(0.00000000001, 0), goodtriangle)
    @test !ispointinsidetriangle(Point(-0.000000001, 0), goodtriangle)

    @test ispointinsidetriangle(Point(50, 50), goodtriangle)
    @test ispointinsidetriangle(Point(10e9, 50), goodtriangle)

    bigtriangle = Point[O, Point(10e9, 0), Point(-10e9, 10e9)]
    @test ispointinsidetriangle(Point(10e4, 10e4), bigtriangle)
    #shouldn't be able to fool this with an anticlockwise triangle
    @test ispointinsidetriangle(Point(10e4, 10e4), reverse(bigtriangle))

    badtriangle = Point[O, Point(100, 0), Point(200, 0)]
    # apparently this is OK... ;|
    @test ispointinsidetriangle(O, badtriangle)

    # perhaps we should catch these not-triangles better than throwing BoundsError
    badtriangle = Point[O, Point(100, 0)]
    @test_throws BoundsError ispointinsidetriangle(O, badtriangle)

    badtriangle = Point[O]
    @test_throws BoundsError ispointinsidetriangle(O, badtriangle)
end

@testset "polytriangulate tests" begin
    pgon = star(O, 200, 12, 0.7, vertices=true)

    sethue("black")
    poly(pgon, :stroke, close=true)

    triangles = polytriangulate!(copy(pgon))

    for t in triangles
        randomhue()
        poly(t, :fill)
    end

    # area of all the little triangles should add up to area of big pgon
    @test isapprox(sum(polyarea.(triangles)), polyarea(pgon))

    @test length(triangles) == 22

    pgon = box(BoundingBox(), vertices=true)
    triangles = polytriangulate!(copy(pgon))

    # areas should add up to area of big pgon
    @test isapprox(sum(polyarea.(triangles)), polyarea(pgon))

    @test length(triangles) == 2
end

@testset "insertvertices tests" begin
    pgon = star(Point(20, 20), 200, 12, 0.7, π/7, vertices=true)
    originallength = length(pgon)
    @test length(insertvertices!(pgon)) == 2originallength
end

@test finish() == true


text:
function surf_props(ebm::EBM, cn::Constants, Sf)

    # Snow albedo

    if ebm.am == 0  # Diagnosed snow albedo
        ebm.albs = ebm.asmn + (ebm.asmx - ebm.asmn) * (ebm.Tsurf - cn.Tm) / ebm.Talb
    elseif ebm.am == 1  # Prognostic snow albedo
        tau = 3600.0 * ebm.tcld
        if (ebm.Tsurf >= cn.Tm)
            tau = 3600.0 * ebm.tmlt
        end
        rt = 1.0 / tau + Sf / ebm.Salb
        alim = (ebm.asmn / tau + Sf * ebm.asmx / ebm.Salb) / rt
        ebm.albs = alim + (ebm.albs - alim) * exp(-rt * ebm.dt)
    end

    if (ebm.albs < min(ebm.asmx, ebm.asmn))
        ebm.albs = min(ebm.asmx, ebm.asmn)
    end

    if (ebm.albs > max(ebm.asmx, ebm.asmn))
        ebm.albs = max(ebm.asmx, ebm.asmn)
    end

    # Density of fresh snow

    if (ebm.dm == 0)
        ebm.rfs = ebm.rho0
    elseif (ebm.dm == 1)
        ebm.rfs = ebm.rhof
    end

    # Thermal conductivity of snow

    ebm.ksnow[:] .= ebm.kfix
    if (ebm.cm == 0) # Fixed
        ebm.ksnow[:] .= ebm.kfix
    elseif (ebm.cm == 1)  # Density function
        for k in 1:ebm.Nsnow
            rhos = ebm.rfs
            if (ebm.dm == 1 && ebm.Ds[k] > eps(Float64))
                rhos = (ebm.Sice[k] + ebm.Sliq[k]) / ebm.Ds[k]
            end
            ebm.ksnow[k] = cn.hcon_ice * (rhos / cn.rho_ice)^ebm.bthr
        end
    end

    # Partial snow cover

    snowdepth = sum(ebm.Ds)
    ebm.fsnow = tanh(snowdepth / ebm.hfsn)
    ebm.alb = ebm.fsnow * ebm.albs + (1.0 - ebm.fsnow) * ebm.alb0
    ebm.z0 = (ebm.z0sn^ebm.fsnow) * (ebm.z0sf^(1.0 - ebm.fsnow))

    # Soil

    dPsidT = - cn.rho_ice * cn.Lf / (cn.rho_wat * cn.g * cn.Tm)

    for k in 1:ebm.Nsoil
        ebm.csoil[k] = ebm.hcap_soil * ebm.Dzsoil[k]
        ebm.ksoil[k] = ebm.hcon_soil
        if (ebm.theta[k] > eps(Float64))
            dthudT = 0.0
            sthu = ebm.theta[k]
            sthf = 0.0
            Tc = ebm.Tsoil[k] - cn.Tm
            Tmax = cn.Tm + (ebm.sathh / dPsidT) * (ebm.Vsat / ebm.theta[k])^ebm.b
            if (ebm.Tsoil[k] < Tmax)
                dthudT = (-dPsidT * ebm.Vsat / (ebm.b * ebm.sathh)) * (dPsidT * Tc / ebm.sathh)^(-1 / ebm.b - 1)
                sthu = ebm.Vsat * (dPsidT * Tc / ebm.sathh)^(-1.0 / ebm.b)
                sthu = min(sthu, ebm.theta[k])
                sthf = (ebm.theta[k] - sthu) * cn.rho_wat / cn.rho_ice
            end
            Mf = cn.rho_ice * ebm.Dzsoil[k] * sthf
            Mu = cn.rho_wat * ebm.Dzsoil[k] * sthu
            ebm.csoil[k] = ebm.hcap_soil * ebm.Dzsoil[k] + cn.hcap_ice * Mf + cn.hcap_wat * Mu + cn.rho_wat * ebm.Dzsoil[k] * ((cn.hcap_wat - cn.hcap_ice) * Tc + cn.Lf) * dthudT
            Smf = cn.rho_ice * sthf / (cn.rho_wat * ebm.Vsat)
            Smu = sthu / ebm.Vsat
            thice = 0.0
            if (Smf > 0.0)
                thice = ebm.Vsat * Smf / (Smu + Smf)
            end
            thwat = 0.0
            if (Smu > 0.0)
                thwat = ebm.Vsat * Smu / (Smu + Smf)
            end
            hcon_sat = ebm.hcon_soil * (cn.hcon_wat^thwat) * (cn.hcon_ice^thice) / (cn.hcon_air^ebm.Vsat)
            ebm.ksoil[k] = (hcon_sat - ebm.hcon_soil) * (Smf + Smu) + ebm.hcon_soil
            if (k == 1)
                ebm.gs = ebm.gsat * max((Smu * ebm.Vsat / ebm.Vcrit)^2.0, 1.0)
            end
        end
    end

    # Surface layer

    ebm.Dz1 = max(ebm.Dzsoil[1], ebm.Ds[1])

    ebm.Ts1 = ebm.Tsoil[1] + (ebm.Tsnow[1] - ebm.Tsoil[1]) * ebm.Ds[1] / ebm.Dzsoil[1]

    ebm.ksurf = ebm.Dzsoil[1] / (2.0 * ebm.Ds[1] / ebm.ksnow[1] + (ebm.Dzsoil[1] - 2 * ebm.Ds[1]) / ebm.ksoil[1])

    if (ebm.Ds[1] > 0.5 * ebm.Dzsoil[1])
        ebm.ksurf = ebm.ksnow[1]
    end

    if (ebm.Ds[1] > ebm.Dzsoil[1])
        ebm.Ts1 = ebm.Tsnow[1]
    end

    return nothing

end


text:
const GD = GeoData

const LayerItr = Union{AbstractArray,Tuple}

"""
    geoarray(T::Type{<:RasterDataSource}, args...; kw...) => AbstractArray

Load a `RasterDataSource` as an `AbstractGeoStack`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoStack`.
"""
function geoarray end
"""
    stack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractGeoStack

Load a `RasterDataSource` as an `AbstractGeoStack`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoStack`.
"""
function stack end
stack(T, layer::Symbol; kw...) = stack(T, (layer,); kw...)

"""
    series(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractGeoSeries

Load a `RasterDataSource` as an `AbstractGeoSeries`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoSeries`.
"""
function series end
series(T, layer::Symbol; kw...) = series(T, (layer,); kw...)

### WorldClim ###

# Weather
function geoarray(T::Type{<:WorldClim{Weather}}, layer::Symbol; date, kw...)
    GDALarray(getraster(T, layer; date=date); kw...)
end
function stack(T::Type{WorldClim{Weather}}, layers::LayerItr=layers(T); date, kw...)
    GDALstack(map(l -> getraster(T, l; date=date), layers); keys=layers, kw...)
end
function series(T::Type{WorldClim{Weather}}, layers::LayerItr=layers(T); date, window=(), kw...)
    step = Month(1)
    dates = _date_sequence(date, step)
    timedim = Ti(dates; mode=Sampled(Ordered(), Regular(step), Intervals(Start())))
    stacks = [stack(T, layers; date=d, window=window) for d in dates]
    GeoData.GeoSeries(stacks, timedim; kw...)
end

# Climate
function geoarray(T::Type{<:WorldClim{Climate}}, layer::Symbol; month, res=defres(T), kw...)
    GDALarray(getraster(T, layer; res=res, month=month); kw...)
end
function stack(T::Type{WorldClim{Climate}}, layers::LayerItr=layers(T); month, res=defres(T), kw...)
    filenames = map(l -> getraster(T, l; res=res, month=month), layers)
    GDALstack(filenames; keys=layers, kw...)
end
function series(T::Type{WorldClim{Climate}}, layers::LayerItr=layers(T);
    res=defres(T), month=1:12, window=(), kw...
)
    timedim = Ti(month; mode=Sampled(span=Regular(1), sampling=Intervals(Start())))
    stacks = [stack(T, layers; res=res, month=m, window=window) for m in month]
    GeoData.GeoSeries(stacks, timedim; kw...)
end

# BioClim
function geoarray(T::Type{<:WorldClim{BioClim}}, layer::Int; res=defres(T), kw...)
    GDALarray(getraster(T, layer; res=res); kw...)
end
function stack(T::Type{WorldClim{BioClim}}, layers::LayerItr=layers(T); res=defres(T), kw...)
    filenames = [getraster(T, l; res=res) for l in layers]
    GDALstack(filenames; keys=_asbioclim(layers), kw...)
end

_asbioclim(ns) = map(_asbioclim, ns)
_asbioclim(n::Int) = string("BIO", n)

#### CHELSA ####

# BioClim
function geoarray(T::Type{<:CHELSA{BioClim}}, layer::Int; kw...)
    GDALarray(getraster(T, layer); kw...)
end
function stack(T::Type{CHELSA{BioClim}}, layers::LayerItr=layers(T); kw...)
    filenames = [getraster(T, l) for l in layers]
    GDALstack(filenames; keys=_asbioclim(layers), kw...)
end

#### EarthEnv ####

# HabitatHeterogeneity
function geoarray(T::Type{<:EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res=defres(T), kw...)
    GDALarray(getraster(T, layer; res=res); kw...)
end
function stack(T::Type{EarthEnv{HabitatHeterogeneity}}, layers::LayerItr=layers(T); res=defres(T), kw...)
    filenames = [getraster(T, l; res=res) for l in layers]
    GDALstack(filenames; keys=layers, kw...)
end

# LandCover
function geoarray(T::Type{<:EarthEnv{LandCover}}, layer::Int; discover=false, kw...)
    GDALarray(getraster(T, layer; discover=discover); kw...)
end
function stack(T::Type{EarthEnv{LandCover}}, layers::LayerItr=layers(T); discover=false, kw...)
    filenames = [getraster(T, l; discover=discover) for l in layers]
    GDALstack(filenames; keys=layers, kw...)
end

#### ALWB ####

function geoarray(T::Type{<:ALWB}, layer::Symbol; date::TimeType, kw...)
    NCDarray(getraster(T, layer; date=date), layer; kw...)
end
function stack(T::Type{<:ALWB}, layers::LayerItr=layers(T); date, kw...)
    filenames = [getraster(T, l; date=date) for l in layers]
    NCDstack(filenames; keys=layers, kw...)
end
function series(T::Type{<:ALWB{M,P}}, layers::LayerItr=layers(T); date, window=(), kw...) where {M,P}
    step = P(1)
    dates = _date_sequence(date, step)
    timedim = Ti(dates; mode=Sampled(span=Regular(step), sampling=Intervals(Start())))
    stacks = [stack(T, layers; date=d, window=window) for d in dates]
    GeoData.GeoSeries(stacks, timedim; kw...)
end

#### AWAP ####

# AWAP asciii files don't have crs and BOM don't specify what it is besides
# being evenly spaced lat/lon. So we assume it's EPSG(4326) and set it manually
function geoarray(T::Type{AWAP}, layer; date, kw...)
    GDALarray(getraster(T, layer; date=date); crs=EPSG(4326), kw...)
end
function stack(T::Type{AWAP}, layers::LayerItr=layers(T); date, kw...)
    GDALstack(map(l -> getraster(T, l; date=date), layers);
        childkwargs=(crs=EPSG(4326),), keys=layers, kw...
    )
end
function series(T::Type{AWAP}, layers::LayerItr=layers(T); date, window=(), kw...)
    step = Day(1)
    dates = _date_sequence(date, step)
    timedim = Ti(dates; mode=Sampled(span=Regular(step), sampling=Intervals(Start())))
    stacks = [stack(T, layers; date=d, window=window) for d in dates]
    GeoData.GeoSeries(stacks, timedim; kw...)
end


text:
function FiniteDifferencesDifferentiator(arg0::jint, arg1::jdouble)
    return FiniteDifferencesDifferentiator((jint, jdouble), arg0, arg1)
end

function FiniteDifferencesDifferentiator(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return FiniteDifferencesDifferentiator((jint, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function differentiate(obj::FiniteDifferencesDifferentiator, arg0::UnivariateFunction)
    return jcall(obj, "differentiate", UnivariateDifferentiableFunction, (UnivariateFunction,), arg0)
end

function differentiate(obj::FiniteDifferencesDifferentiator, arg0::UnivariateMatrixFunction)
    return jcall(obj, "differentiate", UnivariateDifferentiableMatrixFunction, (UnivariateMatrixFunction,), arg0)
end

function differentiate(obj::FiniteDifferencesDifferentiator, arg0::UnivariateVectorFunction)
    return jcall(obj, "differentiate", UnivariateDifferentiableVectorFunction, (UnivariateVectorFunction,), arg0)
end

function get_nb_points(obj::FiniteDifferencesDifferentiator)
    return jcall(obj, "getNbPoints", jint, ())
end

function get_step_size(obj::FiniteDifferencesDifferentiator)
    return jcall(obj, "getStepSize", jdouble, ())
end


text:
# sim_multithread_della.jl: main simulation code with multithreading
# on della
using Distributed
using Dates
using Pkg; Pkg.activate(".")
using CooperationPolarization
addprocs(15; exeflags="--project")

@everywhere begin
    ###################
    ## load packages ##
    ###################
    using Distributed
    using Dates
    using Pkg; Pkg.activate(".")
    using CooperationPolarization
    using BenchmarkTools
    using PyPlot: plt, savefig
    using Statistics: mean, var
    using SharedArrays
    using DelimitedFiles
    println("packages loaded")

end

###########################
## initialize parameters ##
###########################
# model parameters (fixed)
const b = 1.0 # benefit to cooperating
const c = 0.2 # cost to cooperating

const N = 40 # population size
const P = 2 # number of parties

const β = 0.001    # selection strength
# const u = 0.001    # strategy mutation rate
# const v = 0.001    # set mutation rate
# const p = 0.       # party bias (0 = low, 1 = high)
# const ϵ = 1.       # damping factor (0 = no bias in mutation, 1 = max bias (= p) in mutation)
# note: α = 0 automatically

# model parameters (varied)
const M = 5 # number of sets
const K = 3 # number of issues to care about

const pvals = repeat([0.0, 0.25, 0.5, 0.75, 1.0], 20) # uniform populations for now
const us    = [0.001]
const vs    = [0.005]
const ϵ     = 1.
const parameters = collect(Iterators.product(pvals,us,vs))

# simulation parameters
const verbose       = false
const verbose_track = false
const generations   = 20000000
const testtype      = "della"  # specify location

# plotting and data processing parameters
const gen_skip   = 1 # max( div(generations, 10000), 1)
const gen_start  = max( 1, generations-500000 )
const gen_end    = generations
const gen_cutoff = 2000000

println("parameters defined")

const start_time = Dates.format(now(),"yymmdd-HHMMSS.sss")

################
## simulation ##
################
# initialize global tracker
total_summary = SharedArray{Float64, 2}(length(pvals)*length(us)*length(vs),543) # SharedArray{Real}[]

@time @sync @distributed for i in 1:length(pvals)*length(us)*length(vs)

    # set parameter(s)
    p::Float64 = parameters[i][1]
    u::Float64 = parameters[i][2]
    v::Float64 = parameters[i][3]

    # set file name and type
    filename = "b_$(b)_c_$(c)_N_$(N)_M_$(M)_K_$(K)_P_$(P)_beta_$(β)_u_$(u)_v_$(v)_p_$(p)_ϵ_$(ϵ)_gens_$(generations)_$(start_time)"

    ########################################
    ## initialize population and trackers ##
    ########################################
    sets    = random_sets(N, M, K, P)
    game    = Game(b, c, β, u, v, p, ϵ)
    pop     = Population(sets, game, verbose)
    tracker = Tracker(pop, generations, verbose_track)

    println("\nstarting $(filename)")

    ##########################
    ## main simulation loop ##
    ##########################
    for g::Int64 in 1:generations
        # if generations > 10 && (g-1) % (generations ÷ 10) == 0 println("\tinitiating generation $g") end

        evolve!(tracker.pop, 1)
        track!(tracker, g)

        if g == generations
            println("replicate complete, ran $(generations) generations:\n")
            println("\tnum_imitating:\t\t\t$(tracker.pop.num_imitating)")
            println("\tnum_not_imitating:\t\t$(tracker.pop.num_not_imitating)")
            println("\tnum_strategies_changed:  \t$(tracker.pop.num_strategies_changed)")
            println("\t\tnum_random_strategies: \t$(tracker.pop.num_random_strategies)")
            println("\tnum_sets_changed: \t\t$(tracker.pop.num_sets_changed)")
            println("\t\tnum_random_opinions: \t$(tracker.pop.num_random_opinions)")
            println("\t\tnum_biased_opinions: \t$(tracker.pop.num_biased_opinions)")
            println("\tnum_both_changed: \t\t$(tracker.pop.num_both_changed)\n")
        end
    end

    ############################
    ## prep data for plotting ##
    ############################
    # make a summary table
    ens_summary, colnames = make_summary_table( tracker, gen_cutoff, generations )
    total_summary[i,:]    = ens_summary

    # overwrite file after every parameter set (in case code stops due to memory issues)
    # writedlm("data/$(testtype)_run1.csv"; [colnames; total_summary], ',', header=true)
    open("data/gens_$(generations)/$(testtype)_run_multi_M$(M)K$(K)_$(start_time).csv", "w") do f
        write(f, "$(colnames)\n")
        writedlm(f, total_summary, ',')
    end
    println("data processed and exported")

end

#################################
## plotting simulation results ##
#################################
# # plotting
# fig = plotter( tracker, gen_start, gen_skip, gen_end );

# plt.suptitle("Parameters: b=$b, c=$c, N=$N, M=$M, K=$K, β=$β, u=$u, v=$v, q=$q, gamma=$gamma", fontsize=12);
# fig.tight_layout(rect=[0, 0.03, 1, 0.97]);

# # save file
# savefig("plots/gens_$(generations)/$(testtype)_$(filename).pdf")
# plt.close(fig);
# println("plots saved")


text:
mutable struct Location
    x::Int
    y::Int
    heading::Char
end

struct Action
    instruction::Char
    amount::Int
end

function getRoute(file)
    actions::Vector{Action} = []
    for line in readlines(file)
        push!(actions, Action(line[1], parse(Int,line[2:length(line)])))
    end
    #println(actions)
    return actions
end

function turnRight(direction, degrees)
    directions = ['N','E','S','W','N','E','S','W']
    dirIndex = findfirst(isequal(direction), directions)[1]
    quadrants = degrees ÷ 90
    result = directions[dirIndex + quadrants]
    println("Heading ", direction, " turning right for ", degrees, " degrees -> ", result )
    return result
end

function turnLeft(direction, degrees)
    directions = ['E','S','W','N','E','S','W','N']
    dirIndex = findlast(isequal(direction), directions)
    quadrants = degrees ÷ 90
    result = directions[dirIndex - quadrants]
    println("Heading ", direction, " turning left for ", degrees, " degrees -> ", result )
    return result
end

function newLocation(currentLocation, newDirection, amountToMove, changeDirection)

    newLocation = currentLocation
    if newDirection == 'N'
        newLocation.y += amountToMove
    end
    if newDirection == 'W'
        newLocation.x -= amountToMove
    end
    if newDirection == 'S'
        newLocation.y -= amountToMove
    end
    if newDirection == 'E'
        newLocation.x += amountToMove
    end

    newLocation.heading = changeDirection ? newDirection : currentLocation.heading
    #println("newLocation ", newLocation)
    return newLocation
end

function move(currentLocation, action)
    println(action)
    currentDirection = currentLocation.heading
    directionToMove = '_'
    changeDirecton = false
    if action.instruction == 'F'
        directionToMove = currentDirection
        amountToMove = action.amount
        changeDirection = false
    end
    if findfirst(action.instruction, "NWES") != nothing
        directionToMove = action.instruction
        amountToMove = action.amount
        changeDirection = false
    end
    if action.instruction == 'R'
        directionToMove = turnRight(currentDirection, action.amount)
        amountToMove = 0
        changeDirection = true
    end
    if action.instruction == 'L'
        directionToMove = turnLeft(currentDirection, action.amount)
        amountToMove = 0
        changeDirection = true
    end
    println("Instruction=", action.instruction, " directionToMove=", directionToMove, " changeDir=", changeDirection)
    println("Moving ", amountToMove, " to the ", directionToMove)
    return (directionToMove, amountToMove, changeDirection)
end

function navigate(route)
    startingPoint = Location(0,0,'E')
    currentLocation = Location(0,0,'E')
    for action in route
        println(currentLocation)
        directionToMove, amountToMove, changeDirection = move(currentLocation, action)
        currentLocation = newLocation(currentLocation, directionToMove, amountToMove, changeDirection)
        println("")
    end
    traveledNS = abs(currentLocation.y - startingPoint.y)
    traveledEW = abs(currentLocation.x - startingPoint.x)
    manhattanDistance = traveledNS + traveledEW
    println("Manhattan Distance = ", manhattanDistance)
end

route = getRoute("data.txt")
navigate(route)

# function testNewLocation()
#     currentLocation = Location(0,0,'E')
#     currentLocation = newLocation(currentLocation, 'N', 2)
#     println("")
#     currentLocation = newLocation(currentLocation, 'E', 2)
#     println("")
#     currentLocation = newLocation(currentLocation, 'S', 2)
#     println("")
#     currentLocation = newLocation(currentLocation, 'W', 2)
#     println("")
# end
# testNewLocation()

# function testLR()
#     println("N")
#     println("turnLeft('N')=", turnLeft('N') )
#     println("turnRight('N')=", turnRight('N') )
#     println("")
#     println("E")
#     println("turnLeft('E')=", turnLeft('E') )
#     println("turnRight('E')=", turnRight('E') )
#     println("")
#     println("S")
#     println("turnLeft('S')=", turnLeft('S') )
#     println("turnRight('S')=", turnRight('S') )
#     println("")
#     println("W")
#     println("turnLeft('W')=", turnLeft('W') )
#     println("turnRight('W')=", turnRight('W') )
#     println("")

# end
# testLR()


text:
import LineSearches

@testset "Extrapolation" begin
    methods = [LBFGS(),
               ConjugateGradient(),
               LBFGS(alphaguess = LineSearches.InitialQuadratic(),
                     linesearch = LineSearches.BackTracking(order=2))]
    msgs = ["LBFGS Default Options: ",
            "CG Default Options: ",
            "LBFGS + Backtracking + Extrapolation: "]

    if debug_printing
        println("--------------------")
        println("Rosenbrock Example: ")
        println("--------------------")
    end
    rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    for (method, msg) in zip(methods, msgs)
        results = Optim.optimize(rosenbrock, zeros(2), method)
        debug_printing && println(msg, "g_calls = ", results.g_calls, ", f_calls = ", results.f_calls)
    end

    if debug_printing
        println("--------------------------------------")
        println("p-Laplacian Example (preconditioned): ")
        println("--------------------------------------")
    end
    plap(U; n=length(U)) = (n-1) * sum((0.1 .+ diff(U).^2).^2) - sum(U) / (n-1)
    plap1(U; n=length(U), dU = diff(U), dW = 4 .* (0.1 .+ dU.^2) .* dU) =
                            (n - 1) .* ([0.0; dW] .- [dW; 0.0]) .- ones(n) / (n - 1)
    precond(x::Vector) = precond(length(x))
    precond(n::Number) = spdiagm(-1 => -ones(n-1), 0 => 2*ones(n),  1 => -ones(n-1)) * (n+1)
    f(X) = plap([0;X;0])
    g!(g, X) = copyto!(g, (plap1([0;X;0]))[2:end-1])
    N = 100
    initial_x = zeros(N)
    P = precond(initial_x)
    methods = [LBFGS(P=P),
               ConjugateGradient(P=P),
               LBFGS(alphaguess = LineSearches.InitialQuadratic(),
                     linesearch = LineSearches.BackTracking(order=2), P=P)]

    for (method, msg) in zip(methods, msgs)
        results = Optim.optimize(f, g!, copy(initial_x), method)
        debug_printing && println(msg, "g_calls = ", Optim.g_calls(results), ", f_calls = ", Optim.f_calls(results))
    end
end