text:
export ParameterHomotopy, nparameters, set_parameters!

"""
    ParameterHomotopy(F, parameters;
        variables=setdiff(MP.variables(F), parameters),
        p₁=randn(ComplexF64, length(parameters)),
        p₀=randn(ComplexF64, length(parameters)),
        γ₁=nothing, γ₀=nothing)

Construct the homotopy
```math
H(x, t) = F(x, (tγ₁p₁+(1-t)γ₀p₀) / (tγ₁+(1-t)γ₀)),
```
where `p₁` and `p₀` are a vector of parameter values for ``F`` and
`γ₁` and `γ₀` are complex numbers. If `γ₁` or `γ₀` is `nothing`, it is assumed
that `γ₁` and `γ₀` are ``1``.
The input ``parameters`` specifies the parameter variables of ``F``.
Neccessarily, `length(parameters) == length(p₁) == length(p₀)`.

Note that `p₁` and `p₀` are stored as a tuple `p` of `SVectors` and `γ₁` and `γ₀`
are stored as a tuple `γ` or as `γ=nothing`

    ParameterHomotopy(F, parameters;
        variables=setdiff(MP.variables(F), parameters),
        start_parameters=randn(ComplexF64, length(parameters)),
        target_parameters=randn(ComplexF64, length(parameters)),
        start_gamma=nothing, target_gamma=nothing)

This is a non-unicode variant where `γ₁=start_parameters`, `γ₀=target_parameters`,
`γ₁=start_gamma`, `γ₀=target_gamma`.
"""
mutable struct ParameterHomotopy{Sys<:AbstractSystem, T<:Number} <: AbstractHomotopy
    F::Sys
    p::NTuple{2, Vector{T}}
    γ::Union{Nothing, NTuple{2, ComplexF64}}
end

function ParameterHomotopy(F::AbstractSystem, p₁::AbstractVector, p₀::AbstractVector;
    start_gamma=nothing, γ₀ = start_gamma,
    target_gamma=nothing, γ₁ = target_gamma)

    length(p₁) == length(p₀) || error("Length of parameters provided doesn't match.")

    γ = (γ₁ === nothing || γ₀ === nothing) ? nothing : (γ₁, γ₀)
    p = Vector.(promote(p₁, p₀))

    ParameterHomotopy(F, p, γ)
end

function ParameterHomotopy(F::AbstractSystem;
    start_parameters=nothing, p₁ = start_parameters,
    target_parameters=nothing, p₀ = target_parameters,
    kwargs...)
    (p₁ !== nothing && p₀ !== nothing) || error("Parameters not provided for ParameterHomotopy")
    ParameterHomotopy(F, p₁, p₀; kwargs...)
end

function ParameterHomotopy(F::Vector{T},
    parameters::AbstractVector{V};
    variables=setdiff(MP.variables(F), parameters),
    start_parameters=randn(ComplexF64, length(parameters)),
    target_parameters=randn(ComplexF64, length(parameters)),
    kwargs...) where {T<:MP.AbstractPolynomialLike, V<:MP.AbstractVariable}
    G = SPSystem(F; variables=variables, parameters=parameters)
    ParameterHomotopy(G; start_parameters=start_parameters,
        target_parameters=target_parameters, kwargs...)
end

struct ParameterHomotopyCache{C<:AbstractSystemCache, T1<:Number, T2<:Number} <: AbstractHomotopyCache
    F_cache::C
    pt::Vector{T1}
    ∂p∂t::Vector{T1}
    J_p::Matrix{T2}
end

(H::ParameterHomotopy)(x, t, c=cache(H, x, t)) = evaluate(H, x, t, c)

function cache(H::ParameterHomotopy, x, t)
    if H.γ === nothing
        pt = Vector{typeof(t * H.p[1][1])}(undef, length(H.p[1]))
    else
        pt = Vector{typeof(H.γ[1] * t * H.p[1][1])}(undef, length(H.p[1]))
    end
    p!(pt, H, t)
    F_cache = cache(H.F, x, pt)
    ∂p∂t = Vector{eltype(pt)}(undef, length(pt))
    J_p = Matrix(differentiate_parameters(H.F, x, pt, F_cache))

    ParameterHomotopyCache(F_cache, pt, ∂p∂t, J_p)
end

Base.size(H::ParameterHomotopy) = size(H.F)

"""
    nparameters(H::ParameterHomotopy)

Returns the number of parameters of `H`.
"""
nparameters(H::ParameterHomotopy) = length(H.p[1])

"""
    set_parameters!(H::ParameterHomotopy, p::Tuple, γ)

Update the parameters `p` and `γ` of `H`.
"""
function set_parameters!(H::ParameterHomotopy, p::Tuple, γ=nothing)
    H.p[1] .= p[1]
    H.p[2] .= p[2]
    H.γ = γ
    H
end

"""
    set_start_parameters!(H::ParameterHomotopy, p)

Update the start parameters of `H`.
"""
function set_start_parameters!(H::ParameterHomotopy, p)
    H.p[1] .= p
    H
end

"""
    set_target_parameters!(H::ParameterHomotopy, p)

Update the target parameters of `H`.
"""
function set_target_parameters!(H::ParameterHomotopy, p)
    H.p[2] .= p
    H
end

"""
    set_parameters!(H::ParameterHomotopy, p₁, p₀, γ)

Update the parameters `p` and `γ` of `H`.
"""
function set_parameters!(H::ParameterHomotopy, p₁::AbstractVector, p₀::AbstractVector, γ=nothing)
    set_parameters!(H, (p₁, p₀), γ)
end

@inline function p!(pt, H::ParameterHomotopy, t)
    p₁, p₀ = H.p
    if H.γ === nothing
        for i in eachindex(pt)
            @inbounds pt[i] = t * p₁[i] + (1 - t) * p₀[i]
        end
    else
        # compute (tγ₁p₁+(1-t)γ₀p₀) / (tγ₁+(1-t)γ₀)
        γ₁, γ₀ = H.γ
        tγ₁, γ₀_₁₋t = t * γ₁, (1 - t) * γ₀
        γ = (tγ₁ + γ₀_₁₋t)
        a = (@fastmath tγ₁ / γ)
        b = (@fastmath γ₀_₁₋t / γ)
        for i in eachindex(pt)
            @inbounds pt[i] = a * p₁[i] + b * p₀[i]
        end
    end
    pt
end

@inline function ∂p∂t!(u, H::ParameterHomotopy, t, c::ParameterHomotopyCache)
    p₁, p₀ = H.p
    if H.γ === nothing
        for i in eachindex(p₁)
            u[i] = p₁[i] - p₀[i]
        end
    else
        γ₁, γ₀ = H.γ
        tγ₁, γ₀_₁₋t = t * γ₁, (1 - t) * γ₀
        γ = (tγ₁ + γ₀_₁₋t)
        λ = @fastmath γ₁ * γ₀ / (γ * γ)
        for i in eachindex(p₁)
            u[i] = λ * (p₁[i] - p₀[i])
        end
    end
end

function evaluate!(u, H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    evaluate!(u, H.F, x, p!(c.pt, H, t), c.F_cache)
end
function evaluate(H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    evaluate(H.F, x, p!(c.pt, H, t), c.F_cache)
end

function jacobian!(u, H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    jacobian!(u, H.F, x, p!(c.pt, H, t), c.F_cache)
end
function jacobian(H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    jacobian(H.F, x, p!(c.pt, H, t), c.F_cache)
end

function evaluate_and_jacobian!(u, U, H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    evaluate_and_jacobian!(u, U, H.F, x, p!(c.pt, H, t), c.F_cache)
end

function dt!(u, H::ParameterHomotopy, x, t, c::ParameterHomotopyCache)
    # apply chain rule to H(x, p(t))
    p!(c.pt, H, t)
    ∂p∂t!(c.∂p∂t, H, t, c)
    differentiate_parameters!(c.J_p, H.F, x, c.pt, c.F_cache)
    LinearAlgebra.mul!(u, c.J_p, c.∂p∂t)
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

using Base.Meta

issum(s::Symbol) = (s == :sum) || (s == :∑) || (s == :Σ)
isprod(s::Symbol) = (s == :prod) || (s == :∏)

# for backward compatibility with 0.4
function comparison_to_call(ex)
    if isexpr(ex,:comparison) && length(ex.args) == 3
        return Expr(:call,ex.args[2],ex.args[1],ex.args[3])
    else
        return ex
    end
end

function flatten_error(ex)
    isexpr(ex, :flatten) || return
    error("The generator syntax \"$ex\" with multiple \"for\" statements is not yet supported. For JuMP, you should contract statements like \"for i in 1:N for j in 1:i\" into \"for i in 1:N, j in 1:i\"")
end

include("parseExpr_staged.jl")

###############################################################################
# buildrefsets
# Unexported. Takes as input an object representing a name, associated index
# sets, and conditions on those sets, for example
# buildrefsets(:(x[i=1:3,[:red,:blue]],k=S; i+k <= 6))
# Used internally in macros to build JuMPContainers and constraints. Returns
#       refcall:  Expr to reference a particular element, e.g. :(x[i,j,k])
#       idxvars:  Index names used in referencing, e.g.g {:i,:j,:k}
#       idxsets:  Index sets for indexing, e.g. {1:3, [:red,:blue], S}
#       idxpairs: Vector of IndexPair
#       condition: Expr containing any condition present for indexing
# Note in particular that it does not actually evaluate the condition, and so
# it returns just the cartesian product of possible indices.
function buildrefsets(expr::Expr, cname)
    c = copy(expr)
    idxvars = Any[]
    idxsets = Any[]
    idxpairs = IndexPair[]
    # Creating an indexed set of refs
    refcall = Expr(:ref, cname)
    if isexpr(c, :typed_vcat) || isexpr(c, :ref)
        shift!(c.args)
    end
    condition = :()
    if isexpr(c, :vcat) || isexpr(c, :typed_vcat)
        if isexpr(c.args[1], :parameters)
            @assert length(c.args[1].args) == 1
            condition = shift!(c.args).args[1]
        else
            condition = pop!(c.args)
        end
    end

    for s in c.args
        parse_done = false
        if isa(s, Expr)
            parse_done, idxvar, _idxset = tryParseIdxSet(s::Expr)
            if parse_done
                idxset = esc(_idxset)
                push!(idxpairs, IndexPair(idxvar, _idxset))
            end
        end
        if !parse_done # No index variable specified
            idxvar = gensym()
            idxset = esc(s)
            push!(idxpairs, IndexPair(nothing,s))
        end
        push!(idxvars, idxvar)
        push!(idxsets, idxset)
        push!(refcall.args, esc(idxvar))
    end
    return refcall, idxvars, idxsets, idxpairs, condition
end

buildrefsets(c, cname)  = (cname, Any[], Any[], IndexPair[], :())
buildrefsets(c) = buildrefsets(c, getname(c))

###############################################################################
# getloopedcode
# Unexported. Takes a bit of code and corresponding looping information and
# returns that code nested in corresponding loops, along with preceding code
# to construct an appropriate container. Input is:
#       c: symbolic representation of name and appropriate indexing sets, if
#          any. E.g. :(myvar) or :(x[i=1:3,[:red,:blue]])
#       code: inner loop code kernel to be nested in the loops
#       condition: a boolean expression to be evaluated before each kernel.
#                  If none, pass :().
#       idxvars: As defined for buildrefsets
#       idxsets: As defined for buildrefsets
#       idxpairs: As defined for buildrefsets
#       sym: A symbol or expression containing the element type of the
#            resulting container, e.g. :AffExpr or :Variable
function getloopedcode(varname, code, condition, idxvars, idxsets, idxpairs, sym; lowertri=false)

    # if we don't have indexing, just return to avoid allocating stuff
    if isempty(idxsets)
        return code
    end

    hascond = (condition != :())

    if lowertri
        @assert !hascond
        @assert length(idxvars)  == 2
        @assert length(idxpairs) == 2
        @assert !hasdependentsets(idxvars, idxsets)

        i, j = esc(idxvars[1]), esc(idxvars[2])
        expr = copy(code)
        vname = expr.args[1].args[1]
        expr.args[1] = :tmp
        code = quote
            let
                $(localvar(i))
                $(localvar(j))
                for $i in $(idxsets[1]), $j in $(idxsets[2])
                    $i <= $j || continue
                    $expr
                    $vname[$i,$j] = tmp
                    $vname[$j,$i] = tmp
                end
            end
        end
    else
        if hascond
            code = quote
                $(esc(condition)) || continue
                $code
            end
        end
        for (idxvar, idxset) in zip(reverse(idxvars),reverse(idxsets))
            code = quote
                let
                    $(localvar(esc(idxvar)))
                    for $(esc(idxvar)) in $idxset
                        $code
                    end
                end
            end
        end
    end
    if hascond || hasdependentsets(idxvars,idxsets)
        # force a JuMPDict
        N = length(idxsets)
        mac = :($varname = JuMPDict{$(sym),$N}())
    else
        mac = gendict(varname, sym, idxsets...)
    end
    return quote
        $mac
        $code
        nothing
    end
end

localvar(x::Symbol) = _localvar(x)
localvar(x::Expr) = Expr(:block, _localvar(x)...)
_localvar(x::Symbol) = :(local $(esc(x)))
function _localvar(x::Expr)
    @assert x.head == :escape
    args = Any[]
    for t in x.args
        if isa(t, Symbol)
            push!(args, :(local $(esc(t))))
        else
            @assert isa(t, Expr)
            if t.head == :tuple
                append!(args, map(_localvar, t.args))
            else
                error("Internal error defining local variables in macros; please file an issue at https://github.com/JuliaOpt/JuMP.jl/issues/new")
            end
        end
    end
    args
end

getname(c::Symbol) = c
getname(c::Void) = ()
getname(c::AbstractString) = c
function getname(c::Expr)
    if c.head == :string
        return c
    else
        return c.args[1]
    end
end

validmodel(m::AbstractModel, name) = nothing
validmodel(m::MathProgBase.MathProgCallbackData, name) = error("Expected $name to be a JuMP model, but it is a callback object. Use of this macro is not supported within callbacks.")
validmodel(m, name) = error("Expected $name to be a JuMP model, but it has type ", typeof(m))

function assert_validmodel(m, macrocode)
    # assumes m is already escaped
    quote
        validmodel($m, $(quot(m.args[1])))
        $macrocode
    end
end

function _canonicalize_sense(sns::Symbol)
    if sns == :(==)
        return (:(==),false)
    elseif sns == :(>=) || sns == :(≥)
        return (:(>=),false)
    elseif sns == :(<=) || sns == :(≤)
        return (:(<=),false)
    elseif sns == :(.==)
        return (:(==),true)
    elseif sns == :(.>=) || sns == :(.≥)
        return (:(>=),true)
    elseif sns == :(.<=) || sns == :(.≤)
        return (:(<=),true)
    else
        error("Unrecognized sense $sns")
    end
end

function _construct_constraint!(args...)
    warn("_construct_constraint! is deprecated. Use constructconstraint! instead")
    constructconstraint(args...)
end

# two-argument constructconstraint! is used for one-sided constraints.
# Right-hand side is zero.
constructconstraint!(v::Variable, sense::Symbol) = constructconstraint(convert(AffExpr,v), sense)
function constructconstraint!(aff::AffExpr, sense::Symbol)
    offset = aff.constant
    aff.constant = 0.0
    if sense == :(<=) || sense == :≤
        return LinearConstraint(aff, -Inf, -offset)
    elseif sense == :(>=) || sense == :≥
        return LinearConstraint(aff, -offset, Inf)
    elseif sense == :(==)
        return LinearConstraint(aff, -offset, -offset)
    else
        error("Cannot handle ranged constraint")
    end
end

function constructconstraint!(aff::AffExpr, lb, ub)
    offset = aff.constant
    aff.constant = 0.0
    LinearConstraint(aff, lb-offset, ub-offset)
end

constructconstraint!(quad::QuadExpr, sense::Symbol) = QuadConstraint(quad, sense)

function constructconstraint!(normexpr::SOCExpr, sense::Symbol)
    # check that the constraint is SOC representable
    if sense == :(<=)
        SOCConstraint(normexpr)
    elseif sense == :(>=)
        SOCConstraint(-normexpr)
    else
        error("Invalid sense $sense in SOC constraint")
    end
end

constructconstraint!(x::Array, sense::Symbol) = map(c->constructconstraint!(c,sense), x)

_vectorize_like(x::Number, y::Array{AffExpr}) = fill(x, size(y))
function _vectorize_like{R<:Number}(x::Array{R}, y::Array{AffExpr})
    for i in 1:max(ndims(x),ndims(y))
        size(x,i) == size(y,i) || error("Unequal sizes for ranged constraint")
    end
    x
end

function constructconstraint!(x::Array{AffExpr}, lb, ub)
    LB = _vectorize_like(lb,x)
    UB = _vectorize_like(ub,x)
    ret = similar(x, LinearConstraint)
    map!(ret, 1:length(ret)) do i
        constructconstraint!(x[i], LB[i], UB[i])
    end
end

# three-argument constructconstraint! is used for two-sided constraints.
function constructconstraint!(aff::AffExpr, lb::Real, ub::Real)
    offset = aff.constant
    aff.constant = 0.0
    LinearConstraint(aff,lb-offset,ub-offset)
end

constructconstraint!(aff::Variable, lb::Real, ub::Real) = constructconstraint!(convert(AffExpr,v),lb,ub)

constructconstraint!(q::QuadExpr, lb, ub) = error("Two-sided quadratic constraints not supported. (Try @NLconstraint instead.)")

constraint_error(args, str) = error("In @constraint($(join(args,","))): ", str)

macro constraint(args...)
    # Pick out keyword arguments
    if isexpr(args[1],:parameters) # these come if using a semicolon
        kwargs = args[1]
        args = args[2:end]
    else
        kwargs = Expr(:parameters)
    end
    append!(kwargs.args, collect(filter(x -> isexpr(x, :kw), args))) # comma separated
    args = collect(filter(x->!isexpr(x, :kw), args))

    if length(args) < 2
        if length(kwargs.args) > 0
            constraint_error(args, "Not enough positional arguments")
        else
            constraint_error(args, "Not enough arguments")
        end
    end
    m = args[1]
    x = args[2]
    extra = args[3:end]

    m = esc(m)
    # Two formats:
    # - @constraint(m, a*x <= 5)
    # - @constraint(m, myref[a=1:5], a*x <= 5)
    length(extra) > 1 && constraint_error(args, "Too many arguments.")
    # Canonicalize the arguments
    c = length(extra) == 1 ? x        : gensym()
    x = length(extra) == 1 ? extra[1] : x

    anonvar = isexpr(c, :vect) || isexpr(c, :vcat)
    variable = gensym()
    quotvarname = quot(getname(c))
    escvarname  = anonvar ? variable : esc(getname(c))

    if isa(x, Symbol)
        constraint_error(args, "Incomplete constraint specification $x. Are you missing a comparison (<=, >=, or ==)?")
    end

    (x.head == :block) &&
        constraint_error(args, "Code block passed as constraint. Perhaps you meant to use @constraints instead?")
    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end

    # Strategy: build up the code for non-macro addconstraint, and if needed
    # we will wrap in loops to assign to the ConstraintRefs
    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(c, variable)
    # Build the constraint
    if isexpr(x, :call)
        # Simple comparison - move everything to the LHS
        @assert length(x.args) == 3
        (sense,vectorized) = _canonicalize_sense(x.args[1])
        lhs = :($(x.args[2]) - $(x.args[3]))
        addconstr = (vectorized ? :addVectorizedConstraint : :addconstraint)
        newaff, parsecode = parseExprToplevel(lhs, :q)
        constraintcall = :($addconstr($m, constructconstraint!($newaff,$(quot(sense)))))
        for kw in kwargs.args
            @assert isexpr(kw, :kw)
            push!(constraintcall.args, esc(kw))
        end
        code = quote
            q = zero(AffExpr)
            $parsecode
            $(refcall) = $constraintcall
        end
    elseif isexpr(x, :comparison)
        # Ranged row
        (lsign,lvectorized) = _canonicalize_sense(x.args[2])
        (rsign,rvectorized) = _canonicalize_sense(x.args[4])
        if (lsign != :(<=)) || (rsign != :(<=))
            constraint_error(args, "Only ranged rows of the form lb <= expr <= ub are supported.")
        end
        ((vectorized = lvectorized) == rvectorized) || constraint_error("Signs are inconsistently vectorized")
        addconstr = (lvectorized ? :addVectorizedConstraint : :addconstraint)
        x_str = string(x)
        lb_str = string(x.args[1])
        ub_str = string(x.args[5])
        newaff, parsecode = parseExprToplevel(x.args[3],:aff)

        newlb, parselb = parseExprToplevel(x.args[1],:lb)
        newub, parseub = parseExprToplevel(x.args[5],:ub)

        constraintcall = :($addconstr($m, constructconstraint!($newaff,$newlb,$newub)))
        for kw in kwargs.args
            @assert isexpr(kw, :kw)
            push!(constraintcall.args, esc(kw))
        end
        code = quote
            aff = zero(AffExpr)
            $parsecode
            lb = 0.0
            $parselb
            ub = 0.0
            $parseub
        end
        if vectorized
            code = quote
                $code
                lbval, ubval = $newlb, $newub
            end
        else
            code = quote
                $code
                CoefType = coeftype($newaff)
                try
                    lbval = convert(CoefType, $newlb)
                catch
                    constraint_error(args, string("Expected ",$lb_str," to be a ", CoefType, "."))
                end
                try
                    ubval = convert(CoefType, $newub)
                catch
                    constraint_error(args, string("Expected ",$ub_str," to be a ", CoefType, "."))
                end
            end
        end
        code = quote
            $code
            $(refcall) = $constraintcall
        end
    else
        # Unknown
        constraint_error(args, string("Constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2\n" *
              "       expr1 == expr2\n" * "       lb <= expr <= ub"))
    end
    return assert_validmodel(m, quote
        $(getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :ConstraintRef))
        registercon($m, $quotvarname, $variable)
        $(anonvar ? variable : :($escvarname = $variable))
    end)
end

macro SDconstraint(m, x)
    m = esc(m)

    if isa(x, Symbol)
        error("in @SDConstraint: Incomplete constraint specification $x. Are you missing a comparison (<= or >=)?")
    end

    (x.head == :block) &&
        error("Code block passed as constraint.")
    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    isexpr(x,:call) && length(x.args) == 3 || error("in @SDconstraint ($(string(x))): constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2")
    # Build the constraint
    # Simple comparison - move everything to the LHS
    sense = x.args[1]
    if sense == :⪰
        sense = :(>=)
    elseif sense == :⪯
        sense = :(<=)
    end
    sense,_ = _canonicalize_sense(sense)
    lhs = :()
    if sense == :(>=)
        lhs = :($(x.args[2]) - $(x.args[3]))
    elseif sense == :(<=)
        lhs = :($(x.args[3]) - $(x.args[2]))
    else
        error("Invalid sense $sense in SDP constraint")
    end
    newaff, parsecode = parseExprToplevel(lhs, :q)
    assert_validmodel(m, quote
        q = zero(AffExpr)
        $parsecode
        c = SDConstraint($newaff)
        push!($(m).sdpconstr, c)
        c
    end)
end

macro LinearConstraint(x)
    (x.head == :block) &&
        error("Code block passed as constraint. Perhaps you meant to use @LinearConstraints instead?")

    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    if isexpr(x, :call) && length(x.args) == 3
        (sense,vectorized) = _canonicalize_sense(x.args[1])
        # Simple comparison - move everything to the LHS
        vectorized &&
            error("in @LinearConstraint ($(string(x))): Cannot add vectorized constraints")
        lhs = :($(x.args[2]) - $(x.args[3]))
        return quote
            newaff = @Expression($(esc(lhs)))
            c = constructconstraint!(newaff,$(quot(sense)))
            isa(c, LinearConstraint) ||
                error("Constraint in @LinearConstraint is really a $(typeof(c))")
            c
        end
    elseif isexpr(x, :comparison)
        # Ranged row
        (lsense,lvectorized) = _canonicalize_sense(x.args[2])
        (rsense,rvectorized) = _canonicalize_sense(x.args[4])
        if (lsense != :<=) || (rsense != :<=)
            error("in @constraint ($(string(x))): only ranged rows of the form lb <= expr <= ub are supported.")
        end
        (lvectorized || rvectorized) &&
            error("in @LinearConstraint ($(string(x))): Cannot add vectorized constraints")
        lb = x.args[1]
        ub = x.args[5]
        return quote
            if !isa($(esc(lb)),Number)
                error(string("in @LinearConstraint (",$(string(x)),"): expected ",$(string(lb))," to be a number."))
            elseif !isa($(esc(ub)),Number)
                error(string("in @LinearConstraint (",$(string(x)),"): expected ",$(string(ub))," to be a number."))
            end
            newaff = @Expression($(esc(x.args[3])))
            offset = newaff.constant
            newaff.constant = 0.0
            isa(newaff,AffExpr) || error("Ranged quadratic constraints are not allowed")
            LinearConstraint(newaff,$(esc(lb))-offset,$(esc(ub))-offset)
        end
    else
        # Unknown
        error("in @LinearConstraint ($(string(x))): constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2\n" *
              "       expr1 == expr2\n" * "       lb <= expr <= ub")
    end
end

macro QuadConstraint(x)
    (x.head == :block) &&
        error("Code block passed as constraint. Perhaps you meant to use @QuadConstraints instead?")

    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    if isexpr(x, :call) && length(x.args) == 3
        (sense,vectorized) = _canonicalize_sense(x.args[1])
        # Simple comparison - move everything to the LHS
        vectorized &&
            error("in @QuadConstraint ($(string(x))): Cannot add vectorized constraints")
        lhs = :($(x.args[2]) - $(x.args[3]))
        return quote
            newaff = @Expression($(esc(lhs)))
            q = constructconstraint!(newaff,$(quot(sense)))
            isa(q, QuadConstraint) || error("Constraint in @QuadConstraint is really a $(typeof(q))")
            q
        end
    elseif isexpr(x, :comparison)
        error("Ranged quadratic constraints are not allowed")
    else
        # Unknown
        error("in @QuadConstraint ($(string(x))): constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2\n" *
              "       expr1 == expr2")
    end
end

macro SOCConstraint(x)
    (x.head == :block) &&
        error("Code block passed as constraint. Perhaps you meant to use @SOCConstraints instead?")

    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    if isexpr(x, :call) && length(x.args) == 3
        (sense,vectorized) = _canonicalize_sense(x.args[1])
        # Simple comparison - move everything to the LHS
        vectorized &&
            error("in @SOCConstraint ($(string(x))): Cannot add vectorized constraints")
        lhs = :($(x.args[2]) - $(x.args[3]))
        return quote
            newaff = @Expression($(esc(lhs)))
            q = constructconstraint!(newaff,$(quot(sense)))
            isa(q, SOCConstraint) || error("Constraint in @SOCConstraint is really a $(typeof(q))")
            q
        end
    elseif isexpr(x, :comparison)
        error("Ranged second-order cone constraints are not allowed")
    else
        # Unknown
        error("in @SOCConstraint ($(string(x))): constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2")
    end
end

for (mac,sym) in [(:LinearConstraints, Symbol("@LinearConstraint")),
                  (:QuadConstraints,   Symbol("@QuadConstraint")),
                  (:SOCConstraints,    Symbol("@SOCConstraint"))]
    @eval begin
        macro $mac(x)
            x.head == :block || error(string("Invalid syntax for ", $(quot(sym))))
            @assert x.args[1].head == :line
            code = Expr(:vect)
            for it in x.args
                if it.head == :line
                    # do nothing
                elseif it.head == :comparison # regular constraint
                    push!(code.args, Expr(:macrocall, $sym, esc(it)))
                elseif it.head == :tuple # constraint ref
                    if all([isexpr(arg,:comparison) for arg in it.args]...)
                        # the user probably had trailing commas at end of lines, e.g.
                        # @LinearConstraints(m, begin
                        #     x <= 1,
                        #     x >= 1
                        # end)
                        error(string("Invalid syntax in ", $(quot(sym)), ". Do you have commas at the end of a line specifying a constraint?"))
                    end
                    error(string($(quot(sym)), " does not currently support the two argument syntax for specifying groups of constraints in one line."))
                else
                    error("Unexpected constraint expression $it")
                end
            end
            return code
        end
    end
end

for (mac,sym) in [(:constraints,  Symbol("@constraint")),
                  (:NLconstraints,Symbol("@NLconstraint")),
                  (:SDconstraints,Symbol("@SDconstraint")),
                  (:variables,Symbol("@variable")),
                  (:expressions, Symbol("@expression")),
                  (:NLexpressions, Symbol("@NLexpression"))]
    @eval begin
        macro $mac(m, x)
            x.head == :block || error("Invalid syntax for @$mac")
            @assert x.args[1].head == :line
            code = quote end
            for it in x.args
                if isexpr(it, :line)
                    # do nothing
                elseif isexpr(it, :tuple) # line with commas
                    args = []
                    for ex in it.args
                        if isexpr(ex, :tuple) # embedded tuple
                            append!(args, ex.args)
                        else
                            push!(args, ex)
                        end
                    end
                    args_esc = []
                    for ex in args
                        if isexpr(ex, :(=))
                            push!(args_esc,Expr(:kw, ex.args[1], esc(ex.args[2])))
                        else
                            push!(args_esc, esc(ex))
                        end
                    end
                    mac = Expr(:macrocall,$(quot(sym)), esc(m), args_esc...)
                    push!(code.args, mac)
                else # stand-alone symbol or expression
                    push!(code.args,Expr(:macrocall,$(quot(sym)), esc(m), esc(it)))
                end
            end
            push!(code.args, :(nothing))
            return code
        end
    end
end

macro objective(m, args...)
    m = esc(m)
    if length(args) != 2
        # Either just an objective sene, or just an expression.
        error("in @objective: needs three arguments: model, objective sense (Max or Min) and expression.")
    end
    sense, x = args
    if sense == :Min || sense == :Max
        sense = Expr(:quote,sense)
    end
    newaff, parsecode = parseExprToplevel(x, :q)
    code = quote
        q = zero(AffExpr)
        $parsecode
        setobjective($m, $(esc(sense)), $newaff)
    end
    return assert_validmodel(m, code)
end

# Return a standalone, unnamed expression
# ex = @Expression(2x + 3y)
# Currently for internal use only.
macro Expression(x)
    newaff, parsecode = parseExprToplevel(x, :q)
    return quote
        q = 0.0
        $parsecode
        $newaff
    end
end

macro expression(args...)
    if length(args) == 3
        m = esc(args[1])
        c = args[2]
        x = args[3]
    elseif length(args) == 1
        m = nothing
        c = gensym()
        x = args[1]
        Base.warn_once("The one-argument version of @defExpr is deprecated. The corresponding JuMP model is now required as the first argument, and a name for the expression or collection of expressions is required as the second argument. The new syntax is @expression(<JuMP model>, <name of expression(s)>, <expression>)")
    elseif length(args) == 2
        m = nothing
        c = args[1]
        x = args[2]
        Base.warn_once("The two-argument version of @defExpr is deprecated. The corresponding JuMP model is now required as the first argument. The new syntax is @expression(<JuMP model>, <name of expression(s)>, <expression>)")
    else
        error("@expression: needs three arguments.")
    end

    anonvar = isexpr(c, :vect) || isexpr(c, :vcat)
    variable = gensym()
    escvarname  = anonvar ? variable : esc(getname(c))

    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(c, variable)
    newaff, parsecode = parseExprToplevel(x, :q)
    code = quote
        q = 0.0
        $parsecode
    end
    if isa(c,Expr)
        code = quote
            $code
            (isa($newaff,AffExpr) || isa($newaff,Number) || isa($newaff,Variable)) || error("Collection of expressions with @expression must be linear. For quadratic expressions, use your own array.")
        end
    end
    code = quote
        $code
        $(refcall) = $newaff
    end
    code = getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :AffExpr)
    if m === nothing # deprecated usage
        return quote
            $code
            $(anonvar ? variable : :($escvarname = $variable))
        end
    else
        # don't do anything with the model, but check that it's valid anyway
        return assert_validmodel(m, quote
            $code
            $(anonvar ? variable : :($escvarname = $variable))
        end)
    end
end

function hasdependentsets(idxvars, idxsets)
    # check if any index set depends on a previous index var
    for i in 2:length(idxsets)
        for v in idxvars[1:(i-1)]
            if dependson(idxsets[i],v)
                return true
            end
        end
    end
    return false
end

dependson(ex::Expr,s::Symbol) = any(a->dependson(a,s), ex.args)
dependson(ex::Symbol,s::Symbol) = (ex == s)
dependson(ex,s::Symbol) = false
function dependson(ex1,ex2)
    @assert isa(ex2, Expr)
    @assert ex2.head == :tuple
    any(s->dependson(ex1,s), ex2.args)
end

function isdependent(idxvars,idxset,i)
    for (it,idx) in enumerate(idxvars)
        it == i && continue
        dependson(idxset, idx) && return true
    end
    return false
end

esc_nonconstant(x::Number) = x
esc_nonconstant(x) = esc(x)

const EMPTYSTRING = UTF8String("")

variable_error(args, str) = error("In @variable($(join(args,","))): ", str)

macro variable(args...)
    m = esc(args[1])

    extra = vcat(args[2:end]...)
    # separate out keyword arguments
    kwargs = filter(ex->isexpr(ex,:kw), extra)
    extra = filter(ex->!isexpr(ex,:kw), extra)

    # if there is only a single non-keyword argument, this is an anonymous
    # variable spec and the one non-kwarg is the model
    if length(kwargs) == length(args)-1
        x = gensym()
        anon_singleton = true
    else
        x = shift!(extra)
        anon_singleton = false
    end

    t = :Cont
    gottype = false
    haslb = false
    hasub = false
    # Identify the variable bounds. Five (legal) possibilities are "x >= lb",
    # "x <= ub", "lb <= x <= ub", "x == val", or just plain "x"
    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    explicit_comparison = false
    if isexpr(x,:comparison) # two-sided
        explicit_comparison = true
        haslb = true
        hasub = true
        if x.args[2] == :>= || x.args[2] == :≥
            # ub >= x >= lb
            x.args[4] == :>= || x.args[4] == :≥ || variable_error(args, "Invalid variable bounds")
            var = x.args[3]
            lb = esc_nonconstant(x.args[5])
            ub = esc_nonconstant(x.args[1])
        elseif x.args[2] == :<= || x.args[2] == :≤
            # lb <= x <= u
            var = x.args[3]
            (x.args[4] != :<= && x.args[4] != :≤) &&
                variable_error(args, "Expected <= operator after variable name.")
            lb = esc_nonconstant(x.args[1])
            ub = esc_nonconstant(x.args[5])
        else
            variable_error(args, "Use the form lb <= ... <= ub.")
        end
    elseif isexpr(x,:call)
        explicit_comparison = true
        if x.args[1] == :>= || x.args[1] == :≥
            # x >= lb
            var = x.args[2]
            @assert length(x.args) == 3
            lb = esc_nonconstant(x.args[3])
            haslb = true
            ub = Inf
        elseif x.args[1] == :<= || x.args[1] == :≤
            # x <= ub
            var = x.args[2]
            # NB: May also be lb <= x, which we do not support
            #     We handle this later in the macro
            @assert length(x.args) == 3
            ub = esc_nonconstant(x.args[3])
            hasub = true
            lb = -Inf
        elseif x.args[1] == :(==)
            # fixed variable
            var = x.args[2]
            @assert length(x.args) == 3
            lb = esc(x.args[3])
            haslb = true
            ub = esc(x.args[3])
            hasub = true
            gottype = true
            t = :Fixed
        else
            # Its a comparsion, but not using <= ... <=
            variable_error(args, "Unexpected syntax $(string(x)).")
        end
    else
        # No bounds provided - free variable
        # If it isn't, e.g. something odd like f(x), we'll handle later
        var = x
        lb = -Inf
        ub = Inf
    end

    anonvar = isexpr(var, :vect) || isexpr(var, :vcat) || anon_singleton
    anonvar && explicit_comparison && error("Cannot use explicit bounds via >=, <= with an anonymous variable")
    variable = gensym()
    quotvarname = anonvar ? :(:__anon__) : quot(getname(var))
    escvarname  = anonvar ? variable     : esc(getname(var))

    # process keyword arguments
    value = NaN
    obj = nothing
    inconstraints = nothing
    coefficients = nothing
    for ex in kwargs
        kwarg = ex.args[1]
        if kwarg == :start
            value = esc(ex.args[2])
        elseif kwarg == :objective
            obj = esc(ex.args[2])
        elseif kwarg == :inconstraints
            inconstraints = esc(ex.args[2])
        elseif kwarg == :coefficients
            coefficients = esc(ex.args[2])
        elseif kwarg == :basename
            quotvarname = esc(ex.args[2])
        elseif kwarg == :lowerbound
            haslb && variable_error(args, "Cannot specify variable lowerbound twice")
            lb = esc_nonconstant(ex.args[2])
            haslb = true
        elseif kwarg == :upperbound
            hasub && variable_error(args, "Cannot specify variable upperbound twice")
            ub = esc_nonconstant(ex.args[2])
            hasub = true
        else
            variable_error(args, "Unrecognized keyword argument $kwarg")
        end
    end

    if (obj !== nothing || inconstraints !== nothing || coefficients !== nothing) &&
       (obj === nothing || inconstraints === nothing || coefficients === nothing)
        variable_error(args, "Must provide 'objective', 'inconstraints', and 'coefficients' arguments all together for column-wise modeling")
    end

    sdp = any(t -> (t == :SDP), extra)
    symmetric = (sdp || any(t -> (t == :Symmetric), extra))
    extra = filter(x -> (x != :SDP && x != :Symmetric), extra) # filter out SDP and sym tag

    # Determine variable type (if present).
    # Types: default is continuous (reals)
    if length(extra) > 0
        if t == :Fixed
            variable_error(args, "Unexpected extra arguments when declaring a fixed variable")
        end
        if extra[1] in [:Bin, :Int, :SemiCont, :SemiInt]
            gottype = true
            t = extra[1]
        end

        if t == :Bin
            if (lb != -Inf || ub != Inf) && !(lb == 0.0 && ub == 1.0)
            variable_error(args, "Bounds other than [0, 1] may not be specified for binary variables.\nThese are always taken to have a lower bound of 0 and upper bound of 1.")
            else
                lb = 0.0
                ub = 1.0
            end
        end

        !gottype && variable_error(args, "Syntax error")
    end

    # Handle the column generation functionality
    if coefficients !== nothing
        !isa(var,Symbol) &&
        variable_error(args, "Can only create one variable at a time when adding to existing constraints.")

        return assert_validmodel(m, quote
            $variable = Variable($m,$lb,$ub,$(quot(t)),$obj,$inconstraints,$coefficients,UTF8String(string($quotvarname)),$value)
            $(anonvar ? variable : :($escvarname = $variable))
        end)
    end

    if isa(var,Symbol)
        # Easy case - a single variable
        sdp && variable_error(args, "Cannot add a semidefinite scalar variable")
        code = :($variable = Variable($m,$lb,$ub,$(quot(t)),UTF8String(string($quotvarname)),$value))
        if !anonvar
            code = quote
                $code
                registervar($m, $quotvarname, $variable)
                $escvarname = $variable
            end
        end
        return assert_validmodel(m, code)
    end
    isa(var,Expr) || variable_error(args, "Expected $var to be a variable name")

    # We now build the code to generate the variables (and possibly the JuMPDict
    # to contain them)
    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(var, variable)
    clear_dependencies(i) = (isdependent(idxvars,idxsets[i],i) ? () : idxsets[i])

    code = :( $(refcall) = Variable($m, $lb, $ub, $(quot(t)), EMPTYSTRING, $value) )
    if symmetric
        # Sanity checks on SDP input stuff
        condition == :() ||
            variable_error(args, "Cannot have conditional indexing for SDP variables")
        length(idxvars) == length(idxsets) == 2 ||
            variable_error(args, "SDP variables must be 2-dimensional")
        !symmetric || (length(idxvars) == length(idxsets) == 2) ||
            variable_error(args, "Symmetric variables must be 2-dimensional")
        hasdependentsets(idxvars, idxsets) &&
            variable_error(args, "Cannot have index dependencies in symmetric variables")
        for _rng in idxsets
            isexpr(_rng, :escape) ||
                variable_error(args, "Internal error 1")
            rng = _rng.args[1] # undo escaping
            (isexpr(rng,:(:)) && rng.args[1] == 1 && length(rng.args) == 2) ||
                variable_error(args, "Index sets for SDP variables must be ranges of the form 1:N")
        end

        if !(lb == -Inf && ub == Inf)
            variable_error(args, "Semidefinite or symmetric variables cannot be provided bounds")
        end
        return assert_validmodel(m, quote
            $(esc(idxsets[1].args[1].args[2])) == $(esc(idxsets[2].args[1].args[2])) || error("Cannot construct symmetric variables with nonsquare dimensions")
            (Compat.issymmetric($lb) && Compat.issymmetric($ub)) || error("Bounds on symmetric  variables must be symmetric")
            $(getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :Variable; lowertri=symmetric))
            $(if sdp
                quote
                    push!($(m).varCones, (:SDP, first($variable).col : last($variable).col))
                end
            end)
            push!($(m).dictList, $variable)
            registervar($m, $quotvarname, $variable)
            storecontainerdata($m, $variable, $quotvarname,
                               $(Expr(:tuple,idxsets...)),
                               $idxpairs, $(quot(condition)))
            $(anonvar ? variable : :($escvarname = $variable))
        end)
    else
        coloncheckcode = Expr(:call,:coloncheck,refcall.args[2:end]...)
        code = :($coloncheckcode; $code)
        return assert_validmodel(m, quote
            $(getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :Variable))
            isa($variable, JuMPContainer) && pushmeta!($variable, :model, $m)
            push!($(m).dictList, $variable)
            registervar($m, $quotvarname, $variable)
            storecontainerdata($m, $variable, $quotvarname,
                               $(Expr(:tuple,map(clear_dependencies,1:length(idxsets))...)),
                               $idxpairs, $(quot(condition)))
            $(anonvar ? variable : :($escvarname = $variable))
        end)
    end
end

storecontainerdata(m::Model, variable, varname, idxsets, idxpairs, condition) =
    m.varData[variable] = JuMPContainerData(varname, map(collect,idxsets), idxpairs, condition)

macro constraintref(var)
    if isa(var,Symbol)
        # easy case
        return esc(:(local $var))
    else
        if !isexpr(var,:ref)
            error("Syntax error: Expected $var to be of form var[...]")
        end

        varname = var.args[1]
        idxsets = var.args[2:end]
        idxpairs = IndexPair[]

        code = quote
            $(esc(gendict(varname, :ConstraintRef, idxsets...)))
            nothing
        end
        return code
    end
end

macro NLobjective(m, sense, x)
    m = esc(m)
    if sense == :Min || sense == :Max
        sense = Expr(:quote,sense)
    end
    return assert_validmodel(m, quote
        initNLP($m)
        setobjectivesense($m, $(esc(sense)))
        ex = @processNLExpr($m, $(esc(x)))
        $m.nlpdata.nlobj = ex
        $m.obj = zero(QuadExpr)
        $m.internalModelLoaded = false
        nothing
    end)
end

macro NLconstraint(m, x, extra...)
    m = esc(m)
    # Two formats:
    # - @NLconstraint(m, a*x <= 5)
    # - @NLconstraint(m, myref[a=1:5], sin(x^a) <= 5)
    length(extra) > 1 && error("in @NLconstraint: too many arguments.")
    # Canonicalize the arguments
    c = length(extra) == 1 ? x        : gensym()
    x = length(extra) == 1 ? extra[1] : x

    anonvar = isexpr(c, :vect) || isexpr(c, :vcat)
    variable = gensym()
    quotvarname = anonvar ? :(:__anon__) : quot(getname(c))
    escvarname  = anonvar ? variable : esc(getname(c))

    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end

    # Strategy: build up the code for non-macro addconstraint, and if needed
    # we will wrap in loops to assign to the ConstraintRefs
    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(c, variable)
    # Build the constraint
    if isexpr(x, :call) # one-sided constraint
        # Simple comparison - move everything to the LHS
        op = x.args[1]
        if op == :(==)
            lb = 0.0
            ub = 0.0
        elseif op == :(<=) || op == :(≤)
            lb = -Inf
            ub = 0.0
        elseif op == :(>=) || op == :(≥)
            lb = 0.0
            ub = Inf
        else
            error("in @NLconstraint ($(string(x))): expected comparison operator (<=, >=, or ==).")
        end
        lhs = :($(x.args[2]) - $(x.args[3]))
        code = quote
            c = NonlinearConstraint(@processNLExpr($m, $(esc(lhs))), $lb, $ub)
            push!($m.nlpdata.nlconstr, c)
            $(refcall) = ConstraintRef{Model,NonlinearConstraint}($m, length($m.nlpdata.nlconstr))
        end
    elseif isexpr(x, :comparison)
        # ranged row
        if (x.args[2] != :<= && x.args[2] != :≤) || (x.args[4] != :<= && x.args[4] != :≤)
            error("in @NLconstraint ($(string(x))): only ranged rows of the form lb <= expr <= ub are supported.")
        end
        lb = x.args[1]
        ub = x.args[5]
        code = quote
            if !isa($(esc(lb)),Number)
                error(string("in @NLconstraint (",$(string(x)),"): expected ",$(string(lb))," to be a number."))
            elseif !isa($(esc(ub)),Number)
                error(string("in @NLconstraint (",$(string(x)),"): expected ",$(string(ub))," to be a number."))
            end
            c = NonlinearConstraint(@processNLExpr($m, $(esc(x.args[3]))), $(esc(lb)), $(esc(ub)))
            push!($m.nlpdata.nlconstr, c)
            $(refcall) = ConstraintRef{Model,NonlinearConstraint}($m, length($m.nlpdata.nlconstr))
        end
    else
        # Unknown
        error("in @NLconstraint ($(string(x))): constraints must be in one of the following forms:\n" *
              "       expr1 <= expr2\n" * "       expr1 >= expr2\n" *
              "       expr1 == expr2")
    end
    looped = getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :(ConstraintRef{Model,NonlinearConstraint}))
    return assert_validmodel(m, quote
        initNLP($m)
        $m.internalModelLoaded = false
        $looped
        registercon($m, $quotvarname, $variable)
        $(anonvar ? variable : :($escvarname = $variable))
    end)
end

macro NLexpression(args...)
    if length(args) <= 2
        s = IOBuffer()
        print(s,args[1])
        if length(args) == 2
            print(s,",")
            print(s,args[2])
        end
        msg = """
        in @NLexpression($(takebuf_string(s))): three arguments are required.
        Note that the syntax of @NLexpression has recently changed:
        The first argument should be the model to which the expression is attached.
        The second is the name of the expression (or collection of expressions).
        The third is the expression itself.
        Example:
        @NLexpression(m, my_expr, x^2/y)
        @NLexpression(m, my_expr_collection[i=1:2], sin(z[i])^2)
        Support for the old syntax (with the model omitted) will be removed in an upcoming release.
        """
        Base.warn(msg)
        m = :(__last_model[1])
        if length(args) == 2
            c = args[1]
            x = args[2]
        else
            c = gensym()
            x = args[1]
        end
    else
        @assert length(args) == 3
        m, c, x = args
        m = esc(m)
    end

    anonvar = isexpr(c, :vect) || isexpr(c, :vcat)
    variable = gensym()
    escvarname  = anonvar ? variable : esc(getname(c))

    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(c, variable)
    code = quote
        $(refcall) = NonlinearExpression($m, @processNLExpr($m, $(esc(x))))
    end
    return assert_validmodel(m, quote
        $(getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :NonlinearExpression))
        $(anonvar ? variable : :($escvarname = $variable))
    end)
end

# syntax is @NLparameter(m, p[i=1] == 2i)
macro NLparameter(m, ex)
    m = esc(m)
    if VERSION < v"0.5.0-dev+3231"
        ex = comparison_to_call(ex)
    end
    @assert isexpr(ex, :call)
    @assert length(ex.args) == 3
    @assert ex.args[1] == :(==)
    c = ex.args[2]
    x = ex.args[3]

    anonvar = isexpr(c, :vect) || isexpr(c, :vcat)
    if anonvar
        error("In @NLparameter($m, $ex): Anonymous nonlinear parameter syntax is not currently supported")
    end
    variable = gensym()
    escvarname  = anonvar ? variable : esc(getname(c))

    refcall, idxvars, idxsets, idxpairs, condition = buildrefsets(c, variable)
    code = quote
        $(refcall) = newparameter($m, $(esc(x)))
    end
    return assert_validmodel(m, quote
        $(getloopedcode(variable, code, condition, idxvars, idxsets, idxpairs, :NonlinearParameter))
        $(anonvar ? variable : :($escvarname = $variable))
    end)
end

ext: jl


text:
# The references to special cases are to Table of Integrals, Series, and Products, § 9.121, followed by NIST's DLMF.

"""
Compute the Gauss hypergeometric function `₂F₁(a,b;c;z)`.
"""
function _₂F₁(a::Number,b::Number,c::Number,z::Number)
    if real(a) > real(b)
        return _₂F₁(b,a,c,z) # ensure a ≤ b
    elseif isequal(a,c) # 1. 15.4.6
        return exp(-b*log1p(-z))
    elseif isequal(b,c) # 1. 15.4.6
        return exp(-a*log1p(-z))
    elseif isequal(c,0.5)
        if isequal(a+b,0) # 31. 15.4.11 & 15.4.12
            return cosnasinsqrt(2b,z)
        elseif isequal(a+b,1) # 32. 15.4.13 & 15.4.14
            return cosnasinsqrt(1-2b,z)*exp(-0.5log1p(-z))
        elseif isequal(b-a,0.5) # 15.4.7 & 15.4.8
            return expnlog1pcoshatanhsqrt(-2a,z)
        end
    elseif isequal(c,1.5)
        if abeqcd(a,b,0.5) # 13. 15.4.4 & 15.4.5
            return sqrtasinsqrt(z)
        elseif abeqcd(a,b,1) # 14.
            return sqrtasinsqrt(z)*exp(-0.5log1p(-z))
        elseif abeqcd(a,b,0.5,1) # 15. 15.4.2 & 15.4.3
            return sqrtatanhsqrt(z)
        elseif isequal(a+b,1) # 29. 15.4.15 & 15.4.16
            return sinnasinsqrt(1-2b,z)
        elseif isequal(a+b,2) # 30.
            return sinnasinsqrt(2-2b,z)*exp(-0.5log1p(-z))
        elseif isequal(b-a,0.5) # 4. 15.4.9 & 15.4.10
            return expnlog1psinhatanhsqrt(1-2a,z)
        end
    elseif isequal(c,2)
        if abeqcd(a,b,1) # 6. 15.4.1
            return (s = -z; log1p(s)/s)
        elseif a ∈ ℤ && b == 1 # 5.
            return expm1nlog1p(1-a,-z)
        elseif a == 1 && b ∈ ℤ # 5.
            return expm1nlog1p(1-b,-z)
        end
    elseif isequal(c,2.5) && abeqcd(a,b,1,1.5)
         return speciallog(z)
    end
    _₂F₁general(a,b,c,z) # catch-all
end
_₂F₁(a::Number,b::Number,c::Number,z::AbstractArray) = reshape(promote_type(typeof(a),typeof(b),typeof(c),eltype(z))[ _₂F₁(a,b,c,z[i]) for i in eachindex(z) ], size(z))

"""
Compute the Gauss hypergeometric function `₂F₁(a,b;c;z)` with general parameters a,b, and c.
This polyalgorithm is designed based on the paper

    N. Michel and M. V. Stoitsov, Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl–Teller–Ginocchio potential wave functions, Comp. Phys. Commun., 178:535–551, 2008.
"""
function _₂F₁general(a::Number,b::Number,c::Number,z::Number)
    T = promote_type(typeof(a),typeof(b),typeof(c),typeof(z))

    real(b) < real(a) && (return _₂F₁general(b,a,c,z))
    real(c) < real(a)+real(b) && (return exp((c-a-b)*log1p(-z))*_₂F₁general(c-a,c-b,c,z))

    if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
        _₂F₁maclaurin(a,b,c,z)
    elseif abs(z/(z-1)) ≤ ρ && absarg(1-z) < convert(real(T),π)
        exp(-a*log1p(-z))_₂F₁maclaurin(a,c-b,c,z/(z-1))
    elseif abs(inv(z)) ≤ ρ && absarg(-z) < convert(real(T),π)
        _₂F₁Inf(a,b,c,z)
    elseif abs(1-inv(z)) ≤ ρ && absarg(z) < convert(real(T),π) && absarg(1-z) < convert(real(T),π)
        exp(-a*log1p(-z))*_₂F₁Inf(a,c-b,c,z/(z-1))
    elseif abs(1-z) ≤ ρ && absarg(z) < convert(real(T),π) && absarg(1-z) < convert(real(T),π)
        _₂F₁one(a,b,c,z)
    elseif abs(inv(1-z)) ≤ ρ && absarg(-z) < convert(real(T),π)
        exp(-a*log1p(-z))*_₂F₁one(a,c-b,c,z/(z-1))
    elseif absarg(1-z) < convert(real(T),π)
        _₂F₁taylor(a,b,c,z)
    else
        zero(T)#throw(DomainError())
    end
end

"""
Compute the Gauss hypergeometric function `₂F₁(a,b;c;z)` with general parameters a,b, and c.
This polyalgorithm is designed based on the review

    J. W. Pearson, S. Olver and M. A. Porter, Numerical methos for the computation of the confluent and Gauss hypergeometric functions, arXiv:1407.7786, 2014.
"""
function _₂F₁general2(a::Number,b::Number,c::Number,z::Number)
    T = promote_type(typeof(a),typeof(b),typeof(c),typeof(z))
    if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
        _₂F₁maclaurin(a,b,c,z)
    elseif abs(z/(z-1)) ≤ ρ && absarg(1-z) < convert(real(T),π) # 15.8.1
        w = z/(z-1)
        _₂F₁maclaurin(a,c-b,c,w)*exp(-a*log1p(-z))
    elseif abs(inv(z)) ≤ ρ && absarg(-z) < convert(real(T),π)
        w = inv(z)
        if isapprox(a,b) # 15.8.8
            gamma(c)/gamma(a)/gamma(c-a)*(-w)^a*_₂F₁logsumalt(a,c-a,z,w)
        elseif a-b ∉ ℤ # 15.8.2
            gamma(c)*((-w)^a*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a,a-c+1,a-b+1,w)+(-w)^b*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b,b-c+1,b-a+1,w))
        else
            zero(T) # TODO: full 15.8.8
        end
    elseif abs(inv(1-z)) ≤ ρ && absarg(-z) < convert(real(T),π)
        w = inv(1-z)
        if isapprox(a,b) # 15.8.9
            gamma(c)*exp(-a*log1p(-z))/gamma(a)/gamma(c-b)*_₂F₁logsum(a,c-a,z,w,1)
        elseif a-b ∉ ℤ # 15.8.3
            gamma(c)*(exp(-a*log1p(-z))*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a,c-b,a-b+1,w)+exp(-b*log1p(-z))*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b,c-a,b-a+1,w))
        else
            zero(T) # TODO: full 15.8.9
        end
    elseif abs(1-z) ≤ ρ && absarg(z) < convert(real(T),π) && absarg(1-z) < convert(real(T),π)
        w = 1-z
        if isapprox(c,a+b) # 15.8.10
            gamma(c)/gamma(a)/gamma(b)*_₂F₁logsum(a,b,z,w,-1)
        elseif c-a-b ∉ ℤ # 15.8.4
            gamma(c)*(gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a,b,a+b-c+1,w)+exp((c-a-b)*log1p(-z))*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a,c-b,c-a-b+1,w))
        else
            zero(T) # TODO: full 15.8.10
        end
    elseif abs(1-inv(z)) ≤ ρ && absarg(z) < convert(real(T),π) && absarg(1-z) < convert(real(T),π)
        w = 1-inv(z)
        if isapprox(c,a+b) # 15.8.11
            gamma(c)*z^(-a)/gamma(a)*_₂F₁logsumalt(a,b,z,w)
        elseif c-a-b ∉ ℤ # 15.8.5
            gamma(c)*(z^(-a)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a,a-c+1,a+b-c+1,w)+z^(a-c)*(1-z)^(c-a-b)*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a,1-a,c-a-b+1,w))
        else
            zero(T) # TODO: full 15.8.11
        end
    elseif abs(z-0.5) > 0.5
        if isapprox(a,b) && !isapprox(c,a+0.5)
            gamma(c)/gamma(a)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuationalt(a,c,0.5,z)
        elseif a-b ∉ ℤ
            gamma(c)*(gamma(b-a)/gamma(b)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuation(a,a+b,c,0.5,z) + gamma(a-b)/gamma(a)/gamma(c-b)*(0.5-z)^(-b)*_₂F₁continuation(b,a+b,c,0.5,z))
        else
            zero(T)
        end
    else
        zero(T)#throw(DomainError())
    end
end

"""
Compute the generalized hypergeometric function `₃F₂(a₁,a₂,a₃;b₁,b₂;z)` with the Maclaurin series.
"""
function _₃F₂(a₁::Number,a₂::Number,a₃::Number,b₁::Number,b₂::Number,z::Number)
    if abs(z) ≤ ρ
        _₃F₂maclaurin(a₁,a₂,a₃,b₁,b₂,z)
    else
        zero(z)
    end
end
_₃F₂(a₁::Number,a₂::Number,a₃::Number,b₁::Number,b₂::Number,z::AbstractArray) = reshape(promote_type(typeof(a₁),typeof(a₂),typeof(a₃),typeof(b₁),typeof(b₂),eltype(z))[ _₃F₂(a₁,a₂,a₃,b₁,b₂,z[i]) for i in eachindex(z) ], size(z))
_₃F₂(a₁::Number,b₁::Number,z) = _₃F₂(1,1,a₁,2,b₁,z)

"""
Compute the generalized hypergeometric function `mFn(a;b;z)` with the Maclaurin series.
"""
function mFn{S<:Number,V<:Number}(a::AbstractVector{S},b::AbstractVector{V},z::Number)
    if abs(z) ≤ ρ || length(a) ≤ length(b)
        mFnmaclaurin(a,b,z)
    else
        zero(promote_type(S,V,typeof(z)))
    end
end
mFn{S<:Number,V<:Number}(a::AbstractVector{S},b::AbstractVector{V},z::AbstractArray) = reshape(promote_type(S,V,eltype(z))[ mFn(a,b,z[i]) for i in eachindex(z) ], size(z))


text:
#############################################################################
#Copyright (c) 2020 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

"""
	gcvspl_(x,y,ese,SmoothSpline;SplineOrder = Int32(2),SplineMode = Int32(3))

This function allows you to calculate the spline coefficients. It calls gcvspline subroutine of the program GCVSPL.f.

INPUTS:

    x: Array{Float64}, the independent variables;

    y: Array{Float64}, the observations (we assume here that you want to use this spline only on 1 dataset... see gcvspl.f if not);

    ese: Array{Float64}, the errors on y;

    SplineSmooth: Float64, the smoothing factor;

    SplineOrder (M parameter in gcvspl.f): Int32, the half order of the required B-splines. default: splorder = 2 (cubic). SplineOrder = 1,2,3,4 correspond to linear, cubic, quintic, and heptic splines, respectively.

    SplineMode (Int32, MD parameter in gcvspl.f) is the Optimization mode switch:
		default:   SplineMode = 2 (General Cross Validated)

	               SplineMode = 1: Prior given value for p in VAL
	                         (VAL.ge.ZERO). This is the fastest
	                         use of GCVSPL, since no iteration
	                         is performed in p.

	               SplineMode = 2: Generalized cross validation.

	               SplineMode = 3: True predicted mean-squared error,
	                         with prior given variance in VAL.

	               SplineMode = 4: Prior given number of degrees of
	                         freedom in VAL (ZERO.le.VAL.le.N-M).

	               SplineMode  < 0: It is assumed that the contents of
	                         X, W, M, N, and WK have not been
	                         modified since the previous invoca-
	                         tion of GCVSPL. If MD < -1, WK(4)
	                         is used as an initial estimate for
	                         the smoothing parameter p.  At the
	                         first call to GCVSPL, MD must be > 0.

	               Other values for SplineMode, and inappropriate values
	               for VAL will result in an error condition, or
	               cause a default value for VAL to be selected.
	               After return from MD.ne.1, the same number of
	               degrees of freedom can be obtained, for identical
	               weight factors and knot positions, by selecting
	               SplineMode=1, and by copying the value of p from WK(4)
	               into VAL. In this way, no iterative optimization
	               is required when processing other data in Y.

OUPUTS:

	c: Array{Float64}, the spline coefficients;

	WK: Array{Float64}, working vector, see gcvspl.f;

	IER: error parameter

		IER = 0: Normal exit

		IER = 1: M.le.0 .or. N.lt.2*M

		IER = 2: Knot sequence is not strictly increasing, or some weight factor is not positive.

		IER = 3: Wrong mode parameter or value.

SEE GCVSPL.f and Woltring (1986) for even more information.

"""

function gcv_spl(x::Array{Float64,1},y::Array{Float64,1},ese::Array{Float64,1},SmoothSpline::Float64;SplineOrder::Int32 = Int32(2),SplineMode::Int32 = Int32(3))

	if size(y,2)>1
		error("This function does not accept multiple dataset. Please provide only one dataset at a time.")
		#elseif size(x,1) != size(y,1) || size(ese,1) != size(y,1) || size(x,1) != size(ese,1)
		#error("This function only accepts x, y and ese arrays of the same lengths.")
	end

	WX = 1. ./((SmoothSpline.*ese).^2) # relative variance of observations
	WY = zeros(1).+1. # systematic errors... not used so put them to 1
	VAL = (SmoothSpline.*ese).^2

	# M = SplineOrder # No need to redefine but just this comment to keep record
	N = Int32(size(x,1)) #same as length(y)
	K = Int32(1) # number of y columns
	# MD = SplineMode # No need to redefine but just this comment to keep record
	NC = Int32(length(y))

	c = ones(N,NC)
	WK = ones(6 .*(N.*SplineOrder.+1).+N,1)
	IER=Int32(1)

	ccall( (:gcvspl_, gcvlib), Cvoid,
   (Ptr{Float64},Ptr{Float64},Ref{Cint},Ptr{Float64},Ptr{Float64},Ref{Cint},Ref{Cint},Ref{Cint},Ref{Cint},Ptr{Float64},Ref{Float64},Ref{Cint},Ref{Float64},Ref{Cint}),
    x,y,N,WX,WY,SplineOrder,N,K,SplineMode,VAL,c,NC,WK,IER)
	return c, WK, IER
end

"""
	splderivative_julia(xfull::Array{Float64},xparse::Array{Float64},cparse::Array{Float64};SplineOrder::Int32 = Int32(2), L::Int32 = Int32(1), IDER::Int32 = Int32(0))

Wrapper to the SPLDER function of gcvspl.f, for interpolation purpose

INPUTS:

	xfull: Float64 Array, contains the entire x range where the spline has to be evaluated

	xparse: Float64 Array, contains the x values of interpolation regions

	WARNING!!! => xparse[0] <= xfull[0] <= xparse[n]

	cparse: Float64 Array, is the evaluated spline coefficients returned by gcvspl for xparse

OPTIONS:

	splineorder (integer): is the spline order, default: splineorder = 2 (cubic);

	L (integer): see gcvspl.f for details, default: L = 1;

    IDER: the Derivative order required, with 0.le.IDER and IDER.le.2*M. If IDER.eq.0, the function value is returned; otherwise, the IDER-th derivative of the spline is returned.

SEE gcvspl.f FOR MORE INFORMATION
"""

function spl_der(xfull::Array{Float64},xparse::Array{Float64},cparse::Array{Float64};SplineOrder::Int32 = Int32(2), L::Int32 = Int32(1), IDER::Int32 = Int32(0))
	# Note: In the following t, x and c in the fortran code are renamed  xfull, xparse and cparse for clarity
	N = Int32(size(xparse,1))
    q = zeros(2 .*SplineOrder,1)  # working array
    y_calc::Array{Float64} = zeros(size(xfull,1),1) # Output array
    # we loop other xfull to create the output values
	for i =1:size(y_calc,1)
		y_calc[i] = ccall( (:splder_, gcvlib),
        Float64, (Ref{Cint},Ref{Cint},Ref{Cint},Ref{Float64},Ptr{Float64},Ptr{Float64},Ref{Cint},Ptr{Float64}),
        IDER, SplineOrder, N, xfull[i], xparse, cparse, L, q)
	end
	return y_calc
end

max_stars_repo_name: charlesll/gcvspl

max_issues_repo_name: charlesll/gcvspl

max_forks_repo_name: charlesll/gcvspl


text:
# Base type for all predictors
import Base.size

abstract type BasePredictor{T}
end

"""
    add!( p::Predictor{SymbolType}, s::SymbolType)

Adds a symbols `s` to the prediction model.

## Example
```julia-repl
julia> p = LZ78{Char}()
DiscretePredictors.LZ78{Char}([*] (0)
, Char[], Char[])

julia> add!( p, 'd' )
```
"""
function add!( p::BasePredictor{T}, sym::T ) where {T}
    nothing
end

"""
    predict( p::Predictor{SymbolType} )

Returns a dictionary of Dict{SymbolType,Float64} with probabilities of next symbol.

## Examples
```julia-repl
julia> p = LZ78{Int}()
DiscretePredictors.LZ78{Int}([*] (0)
, Int[], Int[])

julia> add!( p, 2 )

julia> add!( p, 3 )

julia> predict( p )
Dict{Int,Float64} with 2 entries:
  2 => 0.5
  3 => 0.5
```
"""
function predict( p::BasePredictor{T} ) where {T}
    return Dict{T,Float64}()
end

"""
    info_string( p::Predictor{SymbolType} )

Returns a string with information about the predictor `p`.

## Example
```julia-repl
julia> p = DG{Char}(4)
DiscretePredictors.DG{Char}([*] (0)
, Char[], 4)

julia> info_string( p )
"DG(4)"
```
"""
function info_string( p::BasePredictor{T} ) where {T}
    return @sprintf( "BasePredictor" )
end

"""
    unique_string( p::Predictor{SymbolType} )

Returns an unique string to identify the predictor; Useful for naming log files.

## Example
```julia-repl
julia> p = DG{Char}(4)
DiscretePredictors.DG{Char}([*] (0)
, Char[], 4)

julia> unique_string( p )
"DG_04"
```
"""
function unique_string( p::BasePredictor{T} ) where {T}
    return @sprintf( "BASE" )
end

"""
    get_best_symbol( p::Predictor{SymbolType} [,default_sym = nothing] )

Returns the most probable symbol from predictor `p` on current context. If no symbol to
predict, then returns `default_sym`.

## Example
```julia-repl
julia> p = DG{Int}(4)
DiscretePredictors.DG{Int}([*] (0)
, Int[], 4)

julia> get_best_symbol( p )

julia> get_best_symbol( p, 2 )
2
```
"""
function get_best_symbol( p::BasePredictor{T}, default_sym = nothing ) where {T}
    symbols     = predict( p )
    return get_best_symbol( symbols, default_sym )
end

function get_best_symbol( symbols::Dict{T,Float64}, default_sym = nothing ) where {T}
    best_symbol = default_sym
    maxP        = 0.00
    for (symbol,prob) in symbols
        if prob > maxP
            best_symbol = symbol
            maxP        = prob
        end
    end
    return best_symbol
end

"""
    size( p::Predictor{SymbolType} )

Returns the number of nodes in the model of predictor `p` including root node.

## Example
```julia-repl
julia> p = KOM{Int}(3)
DiscretePredictors.KOM{Int}([*] (0)
, Int[], 3)

julia> add!(p,1); add!(p,3); add!(p,2); add!(p,2); add!(p,1)

julia> size( p )
13

julia> p.model
[*] (5)
+---[2] (2)
     +---[2] (1)
          +---[1] (1)
     +---[1] (1)
+---[3] (1)
     +---[2] (1)
          +---[2] (1)
               +---[1] (1)
+---[1] (2)
     +---[3] (1)
          +---[2] (1)
               +---[2] (1)
```
"""
function size( p::BasePredictor{T} ) where {T}
    return length(keys(p.model)) + 1
end


text:
using MieSeries
using Base.Test

include("test_harmonics.jl")

# write your own tests here
from_matlab = [
    8.098707839783351
    2.064213995903995
    1.251251398359278
    1.045470773461690
    0.958170336867781
    0.905146609154803
    0.854560436057232
    0.818937250719909
    0.862677950996668
    0.834386518680957]

radius = 1.2
eps0 = 1.5
mu0 = 1.0
sphere = MieSeries.PECSphere(radius, eps0, mu0)
expnsn = MieSeries.excite(sphere, 2π, :left)

Θ = linspace(0.0,   π, 10)
Φ = linspace(0.0, 0.0,  1)
pts = [ [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)] for θ in Θ for ϕ in Φ]

N = length(pts)
E = zeros(N)
H = zeros(N)
for (i,p) in enumerate(pts)
    Ex, Ey, Ez, Hx, Hy, Hz = MieSeries.farfield(sphere, expnsn, p)
    E[i] = norm([Ex,Ey,Ez])
    H[i] = norm([Ex,Ey,Ez])
end

for i in eachindex(E)
    @eval @test_approx_eq_eps E[$i] from_matlab[$i] sqrt(eps(Float64))
end

scattered_from_matlab = [
    0.334493148260097   0.448055290925219   0.474627521463338   0.321583691519947   0.217397836178558
    0.944332487407548                   0   0.982002821143344   0.397084866195122   0.239462161862970
    2.446818207854618                   0                   0   0.455389391557446   0.259887640573314
    0.944332487407548                   0   0.982002821143344   0.397084866195122   0.239462161862970
    0.334493148260097   0.448055290925219   0.474627521463338   0.321583691519947   0.217397836178558]

X = linspace(-2.0, 2.0, 5)
Z = linspace(-2.0, 4.0, 5)
P = [ [x, 0, z] for z in Z for x in X ]

A = zeros(Complex128, length(P))
for (i,p) in enumerate(P)
    norm(p) < 1.2 && continue
    ex, ey, ez, _ = scatteredfield(sphere, expnsn, p)
    A[i] = norm([ex,ey,ez])
end
A = reshape(A, length(X), length(Z))

for i in eachindex(A)
    @eval @test_approx_eq A[$i] scattered_from_matlab[$i]
end


text:
function make_x_and_w(n, p; rng)
    x1 = randn(rng, n, p)
    x2 = randn(rng, n, p) .+ 5
    x3 = randn(rng, n, p) .- 5
    x = table(vcat(x1, x2, x3))
    w = abs.(50 * randn(rng, nrows(x)))
    return x, w
end

function make_xtest(ntest, p; rng)
    xtest1 = randn(rng, ntest, p)
    xtest2 = randn(rng, ntest, p) .+ 5
    xtest3 = randn(rng, ntest, p) .- 5
    xtest = table(vcat(xtest1, xtest2, xtest3))
    return xtest
end


text:
#########################################
## check convex combination constraint ##
#########################################

const precision = 1e-9
function chkEqualsOne(x::Array{Float64, 1})
    if(abs(sum(x) - 1) > precision)
        error("weights must sum to one")
    end
end

function chkEqualsOne(x::Array{Float64, 2})
    if size(x, 2) == 1
        if(abs(sum(x) .- 1) > precision)
            error("weights must sum to one")
        end
    else
        if any(abs(sum(x, 2) .- 1) .> precision)
            error("weights must sum to one")
        end
    end
end

function chkEqualsOne(df::DataFrame)
    rowSums = sum(convert(Array, df), 2)

    if any(abs(rowSums .- 1) .> precision)
        error("weights must sum to one")
    end
end

function chkEqualsOne(tm::Timematr)
    rowSums = sum(core(tm), 2)

    if any(abs(rowSums .- 1) .> precision)
        error("weights must sum to one")
    end
end