text:
export inferType, inferModule, inferExpr, inferDecls

function inferType(globalTC::GlobalTC, localTC::LocalTC, exp::Surf.TyExpr)::HMT
    ln = localTC.ln
    typetype::HMT = @switch exp begin
    @case Surf.TNew(tn)
        T(Nom(Symbol(tn, '@', ln.line, "_", symgen(globalTC))))
    @case Surf.TImplicit(t)
        T(Implicit(inferType(globalTC, localTC, t)))
    @case Surf.TQuery(label, exp)
        let ret = T(inferType(globalTC, localTC, exp))
            push!(globalTC.queries, label => ret)
            ret
        end

    @case Surf.TSym(a)
        T(Nom(a))

    @case Surf.TVar(:_)
        T(globalTC.tcstate.new_tvar())

    @case Surf.TVar(a)
        t = get(localTC.typeEnv, a, nothing)
        t === nothing && throw(MLError(
            localTC.ln,
            UnboundTypeVar(a)))
        t::HMT

    @case Surf.TApp(f, b)
        f = inferType(globalTC, localTC, f)
        b = inferType(globalTC, localTC, b)
        T(App(f, b))

    @case Surf.TTuple(args)
        T(Tup(HMT[inferType(globalTC, localTC, arg) for arg in args]))

    @case Surf.TForall(tvars, p)
        let typeEnv = localTC.typeEnv

            uniqueNames = UN[UN(k) for k in tvars]
            typeEnv = typeEnv[[un.name => T(Bound(un)) for un in uniqueNames]]
            localTC = @set localTC.typeEnv = typeEnv
            p = inferType(globalTC, localTC, p)
            T(Forall(uniqueNames, p))
        end

    @case Surf.TArrow(a, b)
        a = inferType(globalTC, localTC, a)
        b = inferType(globalTC, localTC, b)
        T(Arrow(a, b))
    end

    tcstate = globalTC.tcstate
    unify = tcstate.unify
    new_tvar = tcstate.new_tvar
    prune = tcstate.prune

    @match typetype begin
        T(x) => x
        _ => begin
            targ = new_tvar()
            unify(T(targ), typetype) || throw(MLError(ln, UnificationFail))
            prune(targ)
        end
    end
end

const _InferPostponed = CFunc{IR.Decl, Tuple{}}

function inferModule(globalTC::GlobalTC, localTC::LocalTC, decls::Vector{Surf.Decl})
    props, localTC = inferDecls(globalTC, localTC, decls, Val(true))
    IR.Decl[f() for f in props], localTC
end

function inferDecls(globalTC::GlobalTC, localTC::LocalTC, decls::Vector{Surf.Decl})
    inferDecls(globalTC, localTC, decls, Val(false))
end

function inferDecls(globalTC::GlobalTC, localTC::LocalTC, decls::Vector{Surf.Decl}, ::Val{isGlobal}) where isGlobal
    annotated :: Dict{Symbol, Vector{HM.UN}} = Dict{Symbol, Tuple}()
    loweredDeclFs = _InferPostponed[]
    @inline addDecl!(f::Function) = push!(loweredDeclFs, _InferPostponed(f))

    for decl in decls
        @switch decl begin
        @case Surf.DQuery(label, sym)
            annTy = if sym in localTC.typeEnv
                globalTC.tcstate.prune(localTC.typeEnv[sym])
            else
                throw(MLError(localTC.ln, UnboundVar(sym)))
            end
            push!(globalTC.queries, label=>annTy)
            nothing
        @case Surf.DLoc(ln)
            localTC = @set localTC.ln = ln
            nothing
        @case Surf.DOpen(exp)
            let exprTyProp = inferExpr(globalTC, localTC, exp),
                ex = exprTyProp(NoProp, localTC.localImplicits),
                tcstate = globalTC.tcstate,
                namespaceTy = tcstate.new_tvar(),
                target = App(Nom(:namespace), Tup(HMT[ex.ty, namespaceTy])),
                (fromTy, _) = instanceResolveAndType(
                    globalTC, localTC.localImplicits, target, localTC.ln)

                tcstate.unifyImplicits!(target, fromTy, HMT[]) || throw(MLError(ln, UnificationFail))

                namespaceTy = tcstate.prune(namespaceTy)
                @match namespaceTy begin
                    Tup(fields) =>
                        for field in fields
                            @match field begin
                                Tup([Nom(name), Nom(gensym), t]) =>
                                    begin
                                        localTC = @set localTC.symmap = localTC.symmap[name => gensym]
                                        localTC = @set localTC.typeEnv = localTC.typeEnv[name => t]
                                    end
                                _ => error(MLError(localTC.ln, InvalidNamespaceType(namespaceTy)))
                            end
                        end
                    _ => error(MLError(localTC.ln, InvalidNamespaceType(namespaceTy)))
                end
                IR.Perform(ex)
            end

        @case Surf.DAnn(sym, tyExpr)
            haskey(annotated, sym) && throw(MLError(localTC.ln, UnusedAnnotation(sym)))
            gensym = Symbol(symgen(globalTC), :_, sym)
            localTC = @set localTC.symmap = localTC.symmap[sym => gensym]
            annTy = inferType(globalTC, localTC, tyExpr)
            localTC = @set localTC.typeEnv = localTC.typeEnv[sym => annTy]
            @match annTy begin
                Forall(ns, _) => begin
                        annotated[sym] = ns
                    end
                Implicit(t) =>
                    if isGlobal
                        t = globalTC.tcstate.prune(t)
                        tConsHead = getTConsHead(t)
                        globalImplicits = get!(globalTC.globalImplicits, tConsHead) do
                            InstRec[]
                        end
                        push!(globalImplicits, InstRec(t, localTC.ln, gensym, false))
                        globalTC.globalImplicitDeltas[tConsHead] =
                            get(globalTC.globalImplicitDeltas, tConsHead, 0) + 1

                        annotated[sym] = HM.UN[]
                    else
                        inst = InstRec(t, localTC.ln, gensym, false)
                        localImplicits = cons(inst, localTC.localImplicits)
                        localTC = @set localTC.localImplicits = localImplicits
                        annotated[sym] = HM.UN[]
                    end
                _ => begin annotated[sym] = HM.UN[] end
            end
            nothing
        @case Surf.DBind(sym, expr)
            if sym === :_

                let exprTyProp = inferExpr(globalTC, localTC, expr),
                    localImplicits = localTC.localImplicits

                    addDecl!(() -> IR.Perform(exprTyProp(NoProp, localImplicits)))
                end
                continue
            end
            (gensym, annTy, typevars) = if haskey(annotated, sym)
                typevars = pop!(annotated, sym)
                localTC.symmap[sym], globalTC.tcstate.prune(localTC.typeEnv[sym]), typevars
            else
                gensym = Symbol(symgen(globalTC), :_, sym)
                localTC = @set localTC.symmap = localTC.symmap[sym => gensym]
                tvar = globalTC.tcstate.new_tvar()
                localTC = @set localTC.typeEnv = localTC.typeEnv[sym => tvar]
                gensym, tvar, ()
            end
            # @info :annotation sym annTy
            let localTC = if isempty(typevars)
                            localTC
                        else
                           @set localTC.typeEnv =
                                localTC.typeEnv[[un.name => T(Bound(un)) for un in typevars]]
                        end,
                exprTyProp = inferExpr(globalTC, localTC, expr),
                gensym = gensym,
                annTy = annTy,
                prune = globalTC.tcstate.prune,
                localImplicits = localTC.localImplicits

                addDecl!(
                    () ->
                    IR.Assign(
                        gensym,
                        annTy,
                        exprTyProp(InstTo(annTy), localImplicits)))
            end
        end
    end
    loweredDeclFs, localTC
end

intTypes = Dict{Int, Nom}()
floatTypes = Dict{Int, Nom}()
for bit in (8, 16, 32, 64)
    intTypes[bit] = Nom(Symbol(:i, bit))
    floatTypes[bit] = Nom(Symbol(:f, bit))
end
strT = Nom(:str)
charT = Nom(:char)
boolT = Nom(:bool)

function inferExpr(globalTC::GlobalTC, localTC::LocalTC, expr::Surf.Expr)
    tcstate = globalTC.tcstate
    unifyINST = tcstate.unifyINST
    unifyImplicits! = tcstate.unifyImplicits!
    unify = tcstate.unify
    prune = tcstate.prune
    new_tvar = tcstate.new_tvar
    localImplicits = localTC.localImplicits
    ln = localTC.ln

    @inline function applyTI(ti::TypeInfo, me::HMT)
        local implicits
        implicits = HMT[]
        @match ti begin
            NoProp => (me, implicits)
            InstTo(instTy) => begin
                unifyImplicits!(instTy, me, implicits) || begin
                    # @info :UNIFY_FAIL prune(instTy) prune(me)
                    throw(MLError(ln, UnificationFailOf(prune(instTy), prune(me))))
                end
                (instTy, implicits)
            end
        end
    end

    @inline function prop(varTy::HMT, eImpl::IR.ExprImpl)
        function propInner(ti::TypeInfo, localImplicits::InstResolCtx)
            local implicits
            varTy = prune(varTy)
            varTy, implicits = applyTI(ti, varTy)
            IR.applyImplicits(eImpl, implicits, varTy, ln, localImplicits)
        end
    end

    @switch expr begin
    @case Surf.EExt(jlex)

        function (ti::TypeInfo, localImplicits::InstResolCtx)
            local t, implicits
            t = new_tvar()
            t, implicits = applyTI(ti, t)
            IR.applyImplicits(IR.EExt(jlex), implicits, t, ln, localImplicits)
        end

    @case Surf.EQuery(label, e)
        exprTyProp = inferExpr(globalTC, localTC, e)
        queries = globalTC.queries
        function propQuery(ti::TypeInfo, localImplicits::InstResolCtx)
            e = exprTyProp(ti, localImplicits)
            push!(queries, label=>e.ty)
            e
        end
    @case Surf.ELoc(ln, expr)
        inferExpr(globalTC, (@set localTC.ln = ln), expr)

    @case Surf.EField(subject, field)

        exprTyProp = inferExpr(globalTC, localTC, subject)
        function propField(ti::TypeInfo, localImplicits::InstResolCtx)
            local eSub, retTV, retT, eRet, op_Getter
            eSub = exprTyProp(NoProp, localImplicits)
            retTV = new_tvar()
            target = App(Nom(:field), Tup(HMT[eSub.ty, Nom(field), retTV]))
            retT, eAccess = instanceResolveAndType(globalTC, localImplicits, target, ln)
            unifyImplicits!(target, retT, HMT[]) || throw(MLError(ln, UnificationFail))
            retTV, implicits = applyTI(ti, retTV)

            op_Getter = IR.Expr(ln, nothing, IR.EVar(:op_Getter))
            eRet = IR.EApp(IR.Expr(ln, nothing, IR.EApp(op_Getter, eAccess)), eSub)
            IR.applyImplicits(eRet, implicits, retTV, ln, localImplicits)
        end

    @case Surf.ETup(xs)

        let elty = CFunc{IR.Expr, Tuple{TypeInfo, InstResolCtx}},
            props = elty[elty(inferExpr(globalTC, localTC, x)) for x in xs],
            n_xs = length(xs)

            function propTup(ti::TypeInfo, localImplicits::InstResolCtx)
                local ts, tupT, elts, implicits
                ts = HMT[new_tvar() for i = 1:n_xs]
                tupT = Tup(ts)
                tupT, implicits = applyTI(ti, tupT)
                elts = IR.Expr[(prop(InstTo(t), localImplicits)) for (prop, t) in zip(props, ts)]
                IR.applyImplicits(IR.ETup(elts), implicits, tupT, ln, localImplicits)
            end
        end
    @case Surf.EApp(f, arg)
        let fProp = inferExpr(globalTC, localTC, f),
            argProp = inferExpr(globalTC, localTC, arg)

            function propApp(ti::TypeInfo, localImplicits::InstResolCtx)
                local argT, retT, eF, eArg, implicits
                argT = new_tvar()
                retT = new_tvar()
                retT, implicits = applyTI(ti, retT)
                eF = fProp(InstTo(Arrow(argT, retT)), localImplicits)
                eArg = argProp(InstTo(argT), localImplicits)
                IR.applyImplicits(IR.EApp(eF, eArg), implicits, retT, ln, localImplicits)
            end
        end

    @case Surf.EFun(n, expr)
        gensym = Symbol(symgen(globalTC), :_, n)
        localTC = @set localTC.symmap = localTC.symmap[n => gensym]
        argT = new_tvar()
        localTC = @set localTC.typeEnv = localTC.typeEnv[n => argT]
        retT = new_tvar()
        let exprTyProp = inferExpr(globalTC, localTC, expr),
            gensym = gensym

            function propFun(ti::TypeInfo, localImplicits::InstResolCtx)
                local arrowT, eBody, implicits
                arrowT, implicits = applyTI(ti, Arrow(argT, retT))
                localImplicits = @match prune(argT) begin
                    Implicit(instance) =>
                        cons(InstRec(instance, ln, gensym, false), localImplicits)
                    _ => localImplicits
                end
                eBody = exprTyProp(InstTo(retT), localImplicits)
                IR.applyImplicits(
                    IR.EFun(gensym, eBody),
                    implicits, arrowT, ln, localImplicits)
            end
        end

    @case Surf.ELet(decls, expr)
        let (loweredDeclFs, localTC) = inferDecls(globalTC, localTC, decls),
            exprTyProp = inferExpr(globalTC, localTC, expr)

            function propLet(t::TypeInfo, localImplicits::InstResolCtx)
                local eBody, loweredDecls, eI
                eBody = exprTyProp(t, localImplicits)
                loweredDecls = [f() for f in loweredDeclFs]
                eI = IR.ELet(loweredDecls, eBody)
                IR.Expr(ln, eBody.ty, eI)
            end
        end

    @case Surf.EITE(cond, arm1, arm2)
        let condProp = inferExpr(globalTC, localTC, cond),
            arm1Prop = inferExpr(globalTC, localTC, arm1),
            arm2Prop = inferExpr(globalTC, localTC, arm2),
            new_tvar = globalTC.tcstate.new_tvar

            function propITE(ti::TypeInfo, localImplicits::InstResolCtx)
                local t, eArm1, eArm2, eCond, eI, implicits
                t, implicits = applyTI(ti, new_tvar())
                ti = InstTo(t)
                eArm1 = arm1Prop(ti, localImplicits)
                eArm2 = arm2Prop(ti, localImplicits)
                eCond = condProp(InstTo(boolT), localImplicits)
                eI = IR.EITE(eCond, eArm1, eArm2)
                IR.applyImplicits(eI, implicits, t, ln, localImplicits)
            end
        end

    @case Surf.EVal(val)
        (eValI, eT) = @switch val begin
            @case ::Integer
                bit = 8 * sizeof(val)
                IR.EInt(val, bit), intTypes[bit]
            @case ::AbstractFloat
                bit = 8 * sizeof(val)
                IR.EFloat(val, bit), floatTypes[bit]
            @case ::String
                IR.EStr(val), strT
            @case ::Char
                IR.EChar(val), charT
            @case ::Bool
                IR.EBool(val), boolT
        end
        prop(eT, eValI)

    @case Surf.EVar(n)
        n in localTC.typeEnv || throw(MLError(localTC.ln, UnboundVar(n)))
        varT = prune(localTC.typeEnv[n])
        @match varT begin
            T(t) =>
                function propTVar(ti::TypeInfo, _::InstResolCtx)
                    @match ti begin
                        InstTo(t′) =>
                            begin
                                unify(t′, varT) || throw(MLError(ln, UnificationFail))
                            end
                        _ => nothing
                    end

                    IR.Expr(ln, varT, IR.ETypeVal(t))
                end
            _ =>
                let gensym = localTC.symmap[n]
                    prop(varT, IR.EVar(gensym))
                end
        end
    end
end


text:
module OrthPolys

using SparseArrays
using LinearAlgebra: dot

import ACE

import ACE: evaluate!, evaluate_d!, read_dict, write_dict,
            transform, transform_d, transform_dd, inv_transform,
            ACEBasis, ScalarACEBasis,
            valtype, gradtype,
            acquire!, release!, acquire_B!, release_B!,
            acquire_dB!, release_dB!

using ACE.Transforms: DistanceTransform

import ACE: VectorPool

using ForwardDiff: derivative

import Base: ==

import ChainRules: rrule, NoTangent
import ACE: evaluate, evaluate_d, evaluate_dd,
            _rrule_evaluate, _rrule_evaluate_d

export transformed_jacobi

# these inner functions have been timed to run at
#    6.7ns, 9.2ns, 11.7ns => no need to hand-optimise
_fcut_inner(pl, tl, pr, tr, t) = (t - tl)^pl * (t - tr)^pr

_fcut_d_inner(pl, tl, pr, tr, t) =
      derivative( t -> _fcut_inner(pl, tl, pr, tr, t),  t )

_fcut_dd_inner(pl, tl, pr, tr, t) =
      derivative( t -> _fcut_d_inner(pl, tl, pr, tr, t),  t )

function _fcut_(pl, tl, pr, tr, t)
   if (pl > 0 && t < tl) || (pr > 0 && t > tr)
      return zero(t)
   end
   return _fcut_inner(pl, tl, pr, tr, t)
end

function _fcut_d_(pl, tl, pr, tr, t)
   if (pl > 0 && t < tl) || (pr > 0 && t > tr)
      return zero(t)
   end
   return _fcut_d_inner(pl, tl, pr, tr, t)
end

function _fcut_dd_(pl, tl, pr, tr, t)
   if (pl > 0 && t < tl) || (pr > 0 && t > tr)
      return zero(t)
   end
   return _fcut_dd_inner(pl, tl, pr, tr, t)
end

@doc raw"""
`OrthPolyBasis:` defined a basis of orthonormal polynomials in terms of the
recursion coefficients. What is slightly unusual is that the polynomials have
an "envelope". This results in the recursion
```math
\begin{aligned}
   J_1(x) &= A_1 (x - x_l)^{p_l} (x - x_r)^{p_r} \\
   J_2 &= (A_2 x + B_2) J_1(x) \\
   J_{n} &= (A_n x + B_n) J_{n-1}(x) + C_n J_{n-2}(x)
\end{aligned}
```
Orthogonality is achieved with respect to a user-specified distribution, which
can be either continuous or discrete.

TODO: say more on the distribution! Maybe generalize to non-diagonal
inner products?
"""
struct OrthPolyBasis{T} <: ScalarACEBasis
   # ----------------- the parameters for the cutoff function
   pl::Int        # cutoff power left
   tl::T          # cutoff left (transformed variable)
   pr::Int        # cutoff power right
   tr::T          # cutoff right (transformed variable)
   # ----------------- the recursion coefficients
   A::Vector{T}
   B::Vector{T}
   C::Vector{T}
   # ----------------- used only for construction ...
   #                   but useful to have since it defines the notion of orth.
   tdf::Vector{T}
   ww::Vector{T}
   # -------------
   B_pool::VectorPool{T}
   dB_pool::VectorPool{T}
end

OrthPolyBasis(pl, tl::T, pr, tr::T, A::Vector{T}, B::Vector{T}, C::Vector{T},
              tdf, ww) where {T} =
   OrthPolyBasis(pl, tl, pr, tr, A, B, C, tdf, ww,
                 VectorPool{T}(), VectorPool{T}())

valtype(P::OrthPolyBasis{T}, x::TX = one(T)) where {T, TX <: Number} =
      promote_type(T, TX)

gradtype(P::OrthPolyBasis{T}, x::TX = one(T)) where {T, TX <: Number} =
      promote_type(T, TX)

Base.length(P::OrthPolyBasis) = length(P.A)

==(J1::OrthPolyBasis, J2::OrthPolyBasis) =
      all( getfield(J1, sym) == getfield(J2, sym)
           for sym in (:pr, :tr, :pl, :tl, :A, :B, :C) )

write_dict(J::OrthPolyBasis{T}) where {T} = Dict(
      "__id__" => "ACE_OrthPolyBasis",
      "T" => write_dict(T),
      "pr" => J.pr,
      "tr" => J.tr,
      "pl" => J.pl,
      "tl" => J.tl,
      "A" => J.A,
      "B" => J.B,
      "C" => J.C
   )

OrthPolyBasis(D::Dict, T=read_dict(D["T"])) =
   OrthPolyBasis(
      D["pl"], D["tl"], D["pr"], D["tr"],
      Vector{T}(D["A"]), Vector{T}(D["B"]), Vector{T}(D["C"]),
      T[], T[]
   )

read_dict(::Val{:ACE_OrthPolyBasis}, D::Dict) = OrthPolyBasis(D)

# rand applied to a J will return a random transformed distance drawn from
# the measure w.r.t. which the polynomials were constructed.
# TODO: allow non-constant weights!
function ACE.rand_radial(J::OrthPolyBasis)
   @assert maximum(abs, diff(J.ww)) == 0
   return rand(J.tdf)
end

function OrthPolyBasis(N::Integer,
                       pcut::Integer,
                       tcut::T,
                       pin::Integer,
                       tin::T,
                       tdf::AbstractVector{T},
                       ww::AbstractVector{T} = ones(T, length(tdf))
                       ) where {T <: AbstractFloat}
   @assert pcut >= 0  && pin >= 0
   @assert N > 0

   if tcut < tin
      tl, tr = tcut, tin
      pl, pr = pcut, pin
   else
      tl, tr = tin, tcut
      pl, pr = pin, pcut
   end

   if minimum(tdf) < tl || maximum(tdf) > tr
      @warn("OrthoPolyBasis: t range outside [tl, tr]")
   end

   A = zeros(T, N)
   B = zeros(T, N)
   C = zeros(T, N)

   # normalise the weights s.t. <1, 1> = 1
   ww = ww ./ sum(ww)
   # define inner products
   dotw = (f1, f2) -> dot(f1, ww .* f2)

   # start the iteration
   _J1 = _fcut_.(pl, tl, pr, tr, tdf)
   a = sqrt( dotw(_J1, _J1) )
   A[1] = 1/a
   J1 = A[1] * _J1

   if N > 1
      # a J2 = (t - b) J1
      b = dotw(tdf .* J1, J1)
      _J2 = (tdf .- b) .* J1
      a = sqrt( dotw(_J2, _J2) )
      A[2] = 1/a
      B[2] = -b / a
      J2 = (A[2] * tdf .+ B[2]) .* J1

      # keep the last two for the 3-term recursion
      Jprev = J2
      Jpprev = J1
   end

   for n = 3:N
      # a Jn = (t - b) J_{n-1} - c J_{n-2}
      b = dotw(tdf .* Jprev, Jprev)
      c = dotw(tdf .* Jprev, Jpprev)
      _J = (tdf .- b) .* Jprev -c * Jpprev
      a = sqrt( dotw(_J, _J) )
      A[n] = 1/a
      B[n] = - b / a
      C[n] = - c / a
      Jprev, Jpprev = _J / a, Jprev
   end

   return OrthPolyBasis(pl, tl, pr, tr, A, B, C, collect(tdf), collect(ww))
end

evaluate_P1(J::OrthPolyBasis, t) =
   J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)

function evaluate!(P, J::OrthPolyBasis, t; maxn=length(J))
   @assert length(P) >= maxn
   P[1] = evaluate_P1(J, t)
   if maxn == 1; return P; end
   P[2] = (J.A[2] * t + J.B[2]) * P[1]
   if maxn == 2; return P; end
   @inbounds for n = 3:maxn
      P[n] = (J.A[n] * t + J.B[n]) * P[n-1] + J.C[n] * P[n-2]
   end
   return P
end

function evaluate_d!(dP, J::OrthPolyBasis, t; maxn=length(J))
   @assert maxn <= length(dP)

   P1 = evaluate_P1(J, t)
   dP[1] = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
   if maxn == 1; return dP; end

   α = J.A[2] * t + J.B[2]
   P2 = α * P1
   dP[2] = α * dP[1] + J.A[2] * P1
   if maxn == 2; return dP; end

   @inbounds for n = 3:maxn
      α = J.A[n] * t + J.B[n]
      P3 = α * P2 + J.C[n] * P1
      P2, P1 = P3, P2
      dP[n] = α * dP[n-1] + J.C[n] * dP[n-2] + J.A[n] * P1
   end
   return dP
end

# INCORRECT???
# function evaluate_dd!(ddP, J::OrthPolyBasis, t; maxn=length(J))
#    @assert maxn <= length(dP)

#    P1 = J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)
#    dP1 = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
#    ddP[1] = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
#    if maxn == 1; return dP; end

#    α = J.A[2] * t + J.B[2]
#    P2 = α * P1
#    dP[2] = α * dP[1] + J.A[2] * P1
#    if maxn == 2; return dP; end

#    @inbounds for n = 3:maxn
#       α = J.A[n] * t + J.B[n]
#       P3 = α * P2 + J.C[n] * P1
#       P2, P1 = P3, P2
#       dP[n] = α * dP[n-1] + J.C[n] * dP[n-2] + J.A[n] * P1
#    end
#    return dP
# end

"""
`discrete_jacobi(N; pcut=0, tcut=1.0, pin=0, tin=-1.0, Nquad = 1000)`

A utility function to generate a jacobi-type basis
"""
function discrete_jacobi(N; pcut=0, tcut=1.0, pin=0, tin=-1.0, Nquad = 1000)
   tl, tr = minmax(tin, tcut)
   dt = (tr - tl) / Nquad
   tdf = range(tl + dt/2, tr - dt/2, length=Nquad)
   return OrthPolyBasis(N, pcut, tcut, pin, tin, tdf)
end

# ----------------------------------------------------------------
#   Transformed Polynomials Basis
# ----------------------------------------------------------------

struct TransformedPolys{T, TT, TJ} <: ScalarACEBasis
   J::TJ          # the actual basis
   trans::TT      # coordinate transform
   rl::T          # lower bound r
   ru::T          # upper bound r = rcut
   B_pool::VectorPool{T}
   dB_pool::VectorPool{T}
end

function TransformedPolys(J::OrthPolyBasis{T}, trans, rl, ru)  where {T}
   B_pool = VectorPool{T}()
   return TransformedPolys(J, trans, T(rl), T(ru), B_pool, B_pool)
end

==(J1::TransformedPolys, J2::TransformedPolys) = (
   (J1.J == J2.J) &&
   (J1.trans == J2.trans) &&
   (J1.rl == J2.rl) &&
   (J1.ru == J2.ru) )

TransformedPolys(J, trans, rl, ru) =
   TransformedPolys(J, trans, rl, ru)

write_dict(J::TransformedPolys) = Dict(
      "__id__" => "ACE_TransformedPolys",
      "J" => write_dict(J.J),
      "rl" => J.rl,
      "ru" => J.ru,
      "trans" => write_dict(J.trans)
   )

TransformedPolys(D::Dict) =
   TransformedPolys(
      read_dict(D["J"]),
      read_dict(D["trans"]),
      D["rl"],
      D["ru"]
   )

read_dict(::Val{:ACE_TransformedPolys}, D::Dict) = TransformedPolys(D)

Base.length(J::TransformedPolys) = length(J.J)

valtype(P::TransformedPolys, args...) = valtype(P.J, args...)

gradtype(P::TransformedPolys, args...) = gradtype(P.J, args...)

function ACE.rand_radial(J::TransformedPolys)
   t = ACE.rand_radial(J.J)
   return inv_transform(J.trans, t)
end

# This is NOT the JuLIP or ACEatoms cutoff, but an internal function
_cutoff(J::TransformedPolys) = J.ru

function evaluate!(P, J::TransformedPolys, r; maxn=length(J))
   # transform coordinates
   t = transform(J.trans, r)
   # evaluate the actual polynomials
   evaluate!(P, J.J, t; maxn=maxn)
   return P
end

function evaluate_d!(dP, J::TransformedPolys, r; maxn=length(J))
   # transform coordinates
   t = transform(J.trans, r)
   dt = transform_d(J.trans, r)
   # evaluate the actual Jacobi polynomials + derivatives w.r.t. x
   evaluate_d!(dP, J.J, t, maxn=maxn)
   @. dP *= dt
   return dP
end

evaluate_dd(J::TransformedPolys, r) = derivative(r -> evaluate_d(J, r), r)

"""
`transformed_jacobi(maxdeg, trans, rcut, rin = 0.0; kwargs...)` : construct
a `TransformPolys` basis with an inner polynomial basis of `OrthPolys` type.

* `maxdeg` : maximum degree
* `trans` : distance transform; normally `PolyTransform(...)`
* `rin, rcut` : inner and outer cutoff

**Keyword arguments:**

* `pcut = 2` : cutoff parameter
* `pin = 0` : inner cutoff parameter
* `Nquad = 1000` : number of quadrature points
"""
function transformed_jacobi(maxdeg::Integer,
                            trans::DistanceTransform,
                            rcut::Real, rin::Real = 0.0;
                            kwargs...)
   J =  discrete_jacobi(maxdeg; tcut = transform(trans, rcut),
                                tin = transform(trans, rin),
                                pcut = 2,
                                kwargs...)
   return TransformedPolys(J, trans, rin, rcut)
end

# ------------- AD

function _rrule_evaluate(J::OrthPolyBasis, t::Number,
                         w::AbstractVector{<: Number})
   maxn = length(w)
   @assert maxn <= length(J)

   P1 = J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)
   dP1 = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
   a = dP1 * w[1]
   if maxn == 1
      return a
   end

   α = J.A[2] * t + J.B[2]
   P2 = α * P1
   dP2 = α * dP1 + J.A[2] * P1
   a += dP2 * w[2]
   if maxn == 2
      return a
   end

   @inbounds for n = 3:maxn
      α = J.A[n] * t + J.B[n]
      P3 = α * P2 + J.C[n] * P1
      dP3 = α * dP2 + J.C[n] * dP1 + J.A[n] * P2
      a += dP3 * w[n]
      P2, P1 = P3, P2
      dP2, dP1 = dP3, dP2
   end

   return a
end

function _rrule_evaluate(P::TransformedPolys, x::Number,
                         dx::AbstractVector{<: Number})
   t = transform(P.trans, x)
   dt = transform_d(P.trans, x)
   a = _rrule_evaluate(P.J, t, dx)
   return a * dt
end

function rrule(::typeof(evaluate), P::TransformedPolys, x::Number)
   B = evaluate(P, x)
   return B, dx -> (NoTangent(), NoTangent(), _rrule_evaluate(P, x, dx))
end

function _rrule_evaluate_d(J::OrthPolyBasis, t::Number,
                           w::AbstractVector{<: Number},
                           dt = 1.0, ddt = 0.0)
   maxn = length(w)
   @assert maxn <= length(J)

   P1 = J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)
   dP1 = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
   ddP1 = J.A[1] * _fcut_dd_(J.pl, J.tl, J.pr, J.tr, t)
   a = (ddP1 * dt^2 + dP1 * ddt) * w[1]
   if maxn == 1
      return a
   end

   α = J.A[2] * t + J.B[2]
   P2 = α * P1
   dP2 = α * dP1 + J.A[2] * P1
   ddP2 = α * ddP1 + 2 * J.A[2] * dP1
   a += (ddP2 * dt^2 + dP2 * ddt) * w[2]
   if maxn == 2
      return a
   end

   @inbounds for n = 3:maxn
      α = J.A[n] * t + J.B[n]
      P3 = α * P2 + J.C[n] * P1
      dP3 = α * dP2 + J.C[n] * dP1 + J.A[n] * P2
      ddP3 = α * ddP2 + J.C[n] * ddP1 + 2 * J.A[n] * dP2
      a += (ddP3 * dt^2 + dP3 * ddt) * w[n]
      P2, P1 = P3, P2
      dP2, dP1 = dP3, dP2
      ddP2, ddP1 = ddP3, ddP2
   end
   return a
end

function _rrule_evaluate_d(P::TransformedPolys, x::Number,
                           w::AbstractVector{<: Number})
   t = transform(P.trans, x)
   dt = transform_d(P.trans, x)
   ddt = transform_dd(P.trans, x)
   return _rrule_evaluate_d(P.J, t, w, dt, ddt)
end

function rrule(::typeof(evaluate_d), P::TransformedPolys, x::Number)
   dB = evaluate_d(P, x)
   return dB, dx -> (NoTangent(), NoTangent(), _rrule_evaluate_d(P, x, dx))
end

end

ext: jl


text:
using CSV
using DataFrames

main_fn = ARGS[1]
aft_fn = ARGS[2]

main_csv = CSV.read(main_fn)
w = open(aft_fn, "w")
write(w, main_csv)
close(w)


text:
# new_LP.jl
# New LP formulation for travel time estimation
# Authored by Arthur J Delarue on 7/2/15

MODEL = "relaxed_smooth"		# relaxed/strict _ smooth/nothing/random
LAST_SMOOTH = false 			# true if last iteration is smooth, false otherwise. Irrelevant if model is smooth already.
ERROR_COMPUTATION = "average" 		# average/single-ride/both

MAX_ROUNDS = 2					# max number of iterations

MIN_RIDES = 1					# min number of rides
RADIUS = 140					# radius (4D Euclidian)
YEAR = 2013
START_TIME = 12
END_TIME = 14
START_MONTH = 1
END_MONTH = 1
WEEKDAYS = true 				# true for weekdays, false for weekends
METHOD = "average"				# average/nearest [neighbor]

RANDOM_CONSTRAINTS = false 		# true for purely random constraints, false otw
DYNAMIC_CONSTRAINTS = true 		# true for dynamic constraints, false otw
GLOBAL_CONSTRAINT_UPDATE = false # true if we remove worst paths globally, false otw
SAMPLE_SIZE = 1000 				# starting number of constraints
NUM_OD_ADDED = 1000 			# number of (O,D) pairs to add
NUM_OD_REMOVED = 1000			# number of (O,D) pairs to remove
UPDATE_EVERY_N_ITERATIONS = 1 	# add number of (O,D) above every $N iterations
ADD_IF_ABOVE_PERCENTILE = 0.7		# only add rides that have an error above this percentile
REMOVE_IF_BELOW_PERCENTILE = 0.3	# only remove rides that have an error below this percentile

TURN_COST = 0.	 				# turning cost initial value
TURN_COST_AS_VARIABLE = false 	# true if LP updates turning cost, false otw

DELTA_BOUND = [0.06, 0.05, 0.04]	# relaxation base (sketchy)

MAX_NUM_PATHS_PER_OD = 3 		# at least 1 please

COMPUTE_FINAL_SHORTEST_PATHS = true 	# Should always be true unless trying to save time somehow

START_SIMPLE = false 			# true if initial simple LP is used, false otw

METROPOLIS = false 				# always set this to false unless you wish to crash everything

WARMSTART = false

LOCALIZED_TESTING = false

function fast_LP(
	manhattan::Manhattan,
	travelTimes::Array{Float64,2},
	numRides::Array{Int,2},
	testingData::Union{DataFrame, Array{Float64,2}},
	startTimes::AbstractArray{Float64,2};
	errorComputation::AbstractString=ERROR_COMPUTATION,
	model_type::AbstractString=MODEL,
	last_smooth::Bool=LAST_SMOOTH,
	max_rounds::Int=MAX_ROUNDS,
	min_rides::Int=MIN_RIDES,
	radius::Int=RADIUS,
	year::Int=YEAR,
	startTime::Int=START_TIME,
	endTime::Int=END_TIME,
	startMonth::Int=START_MONTH,
	endMonth::Int=END_MONTH,
	method::AbstractString=METHOD,
	sample_size::Int=SAMPLE_SIZE,
	turnCost::Float64=TURN_COST,
	turnCostAsVariable::Bool=TURN_COST_AS_VARIABLE,
	delta_bound::Union{Array{Float64,1},Array{Any,1}}=DELTA_BOUND,
	maxNumPathsPerOD::Int=MAX_NUM_PATHS_PER_OD,
	computeFinalSP::Bool = COMPUTE_FINAL_SHORTEST_PATHS,
	startWithSimpleLP::Bool=START_SIMPLE,
	randomConstraints::Bool=RANDOM_CONSTRAINTS,
	dynamicConstraints::Bool=DYNAMIC_CONSTRAINTS,
	numPairsToAdd::Int = NUM_OD_ADDED,
	numPairsToRemove::Int = NUM_OD_REMOVED,
	iterationMultiple::Int = UPDATE_EVERY_N_ITERATIONS,
	addOnlyIfAbovePercentile::Float64 = ADD_IF_ABOVE_PERCENTILE,
	removeOnlyIfBelowPercentile::Float64 = REMOVE_IF_BELOW_PERCENTILE,
	globalConstraintUpdate::Bool=GLOBAL_CONSTRAINT_UPDATE,
	warmstart::Bool=WARMSTART,
	metropolis::Bool=METROPOLIS,
	localized_testing::Bool=LOCALIZED_TESTING,
	real_TOD_metropolis::AbstractArray{Float64}=zeros(1,1),
	real_tij_metropolis::AbstractArray{Float64}=zeros(1,1),
	prob::Float64=0.0,
	testingData2::Array{Float64,2}=zeros(1,1))

	graph = manhattan.network
	roadTimes = deepcopy(manhattan.roadTime)
	distances = manhattan.distances
	positions = manhattan.positions
	roads = edges(graph)
	nodes = vertices(graph)
	out = [copy(out_neighbors(graph,i)) for i in nodes]
	inn = [copy(in_neighbors(graph,i)) for i in nodes]

	# Find nonzero travel times
	if DYNAMIC_CONSTRAINTS
		numDataPoints = sample_size
		pairs = [copy(Int[]) for i = 1:length(nodes)]
	else
		pairs = [find(travelTimes[i,:]) for i=nodes]
		numDataPoints = sum([length(pairs[i]) for i=nodes])
	end

	# Create output directory name (sorry this is so complicated)
	if ((!last_smooth) && (!metropolis))
		TESTDIR = "fast_r$(radius)_minr$(min_rides)_i$(max_rounds)_wd_$(startTime)$(endTime)_m$(startMonth)$(endMonth)_y$(year)_$(model_type)_ppc$(maxNumPathsPerOD)_data$(method)_error$errorComputation"
	elseif last_smooth && (!metropolis)
		TESTDIR = "fast_r$(radius)_minr$(min_rides)_i$(max_rounds)_wd_$(startTime)$(endTime)_m$(startMonth)$(endMonth)_y$(year)_$(model_type)_lsmooth_ppc$(maxNumPathsPerOD)_data$(method)_error$errorComputation"
	elseif (!last_smooth) && metropolis
		TESTDIR = "fast_metropolis$(prob)_minr$(min_rides)_i$(max_rounds)_$(model_type)_ppc$(maxNumPathsPerOD)"
	elseif last_smooth && metropolis
		TESTDIR = "fast_metropolis$(prob)_minr$(min_rides)_i$(max_rounds)_$(model_type)_lsmooth_ppc$(maxNumPathsPerOD)"
	end
	if randomConstraints
		TESTDIR=string(TESTDIR, "_rnd_rides$(sample_size)")
	elseif dynamicConstraints
		TESTDIR=string(TESTDIR, "_dynConstr_st$(sample_size)_add$(numPairsToAdd)_thr$(addOnlyIfAbovePercentile)_rem$(numPairsToRemove)_thr$(removeOnlyIfBelowPercentile)_every$(iterationMultiple)")
	end
	if globalConstraintUpdate
		TESTDIR=string(TESTDIR, "_gcu")
	end
	if turnCostAsVariable
		TESTDIR=string(TESTDIR, "_tcvar_start$(turnCost)")
	else
		TESTDIR=string(TESTDIR, "_tc$(turnCost)")
	end
	if startWithSimpleLP
		TESTDIR=string(TESTDIR, "_ss")
	end
	if warmstart
		TESTDIR=string(TESTDIR, "_ws")
	end
	if localized_testing
		TESTDIR=string(TESTDIR, "_loctest")
	end

	# Create directory if necessary:
	if !isdir("Outputs/$TESTDIR")
		mkdir("Outputs/$TESTDIR")
	end
	# Add vital information to file
	writeDataToFile("Outputs/$TESTDIR", model_type, max_rounds, min_rides, radius, false, 0, sample_size, turnCost)
	# Turn cost file
	if turnCostAsVariable
		turnCostFile = open("Outputs/$(TESTDIR)/tc.csv", "w")
	end
	# Metropolis error output
	if metropolis
		# specific error output
	else
		errorFile = open("Outputs/$(TESTDIR)/errorstats.csv", "w")
	end
	costFunctionFile = open("Outputs/$(TESTDIR)/cf.csv", "w")
	# Runtime info output
	timeFile = open("Outputs/$(TESTDIR)/timestats.csv", "w")
	if dynamicConstraints
		numConstraintsFile = open("Outputs/$(TESTDIR)/numODpairs.csv", "w")
	end

	# Tell us where output is going
	println("-- Saving outputs to directory Outputs/$(TESTDIR)/")

	# Create path storing array and initialize. By convention totalPaths[i,1] is the equality constraint
	numPaths = zeros(Int, numDataPoints)
	totalPaths = Array(Any, numDataPoints)
	totalNumExpensiveTurns = Array(Any, numDataPoints)
	for i = 1:numDataPoints
		totalPaths[i] = Array{Int}[]
		totalNumExpensiveTurns[i] = Int[]
		sizehint!(totalPaths[i], maxNumPathsPerOD)
		sizehint!(totalNumExpensiveTurns[i], maxNumPathsPerOD)
	end

	if errorComputation == "single-ride" || errorComputation == "both"
		println("**** Building KDTree of nodes for testing ****")
		nodeTree, nodePairs = buildNodeKDTree(manhattan.positions)
	end

	println("**** Initializing dataset ****")
	if dynamicConstraints
		srcs, dsts, newTravelTimes, newNumRides, pairs = chooseStartingConstraints(travelTimes, numRides, pairs, sample_size = sample_size, num_nodes = length(nodes))
	else
		# Run over all pairs of nodes that have data
		srcs = Int[]
		dsts = Int[]
		sizehint!(srcs, sample_size)
		sizehint!(dsts, sample_size)
		for i in nodes, j in nodes
			if travelTimes[i,j] > 0
				# Load sources and destinations: will be useful when adding constraints
				push!(srcs, i)
				push!(dsts, j)
			end
		end
	end

	status = 0
	newTimes = startTimes
	old_objective = 0

	# Compute shortest paths (with turn cost)
	if startWithSimpleLP
		println("**** Finding initial set of paths ****")
		status, newTimes, turnCost, totalPaths, totalNumExpensiveTurns = simple_LP(manhattan, travelTimes, numRides, testingData, TESTDIR, model_type=split(model_type, "_")[1], max_rounds=max_rounds, turnCost=turnCost, turnCostAsVariable=turnCostAsVariable, delta_bound=delta_bound, maxNumPathsPerOD=maxNumPathsPerOD, computeFinalSP=false)
	else
		println("**** Computing shortest paths ****")
		new_graph, new_edge_dists, new_nodes, new_edge_isExpensive = modifyGraphForDijkstra(graph, newTimes, positions, turn_cost=turnCost)
		old_nodes = getInverseMapping(new_nodes, nv(new_graph))
		@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
	end

	# l is round number
	l = 1
	while l <= max_rounds
		println("###### ROUND $l ######")
		tic()
		actual_l = l
		# Update bound on relaxation variables
		if split(model_type, "_")[1] == "relaxed"
			delta = (1/2)^(div(actual_l-1,length(delta_bound)))*delta_bound[((actual_l-1)%length(delta_bound)) + 1]
			# Avoid numerical issues
			if delta < 1e-4 || l == max_rounds
				delta = 0
			end
		end

		# Create JuMP model for LP and QP
		println("**** Creating LP instance ****")
		if l == 1 && !startWithSimpleLP && !dynamicConstraints
			tolerance = 1e-4
		else
			tolerance = 1e-6
		end
		m = Model(solver=GurobiSolver(TimeLimit=10000, OutputFlag=1, Method=3, BarConvTol=tolerance))

		# Add one variable for each road, for each model
		if !metropolis
			@defVar(m, t[i=nodes,j=out[i]] >= 0.038 * manhattan.distances[i,j])
		else
			@defVar(m, t[i=nodes,j=out[i]] >= startTimes[i,j])
		end

		# Add decision variable for turn cost if necessary
		if turnCostAsVariable
			@defVar(m, tc >= 0)
		else
			tc = turnCost
		end

		# Epsilon is the variable used for the absolute value difference between T and \hat{T}
		@defVar(m, epsilon[i=nodes,j=pairs[i]] >= 0)

		# Create handles for inequality constraints
		if maxNumPathsPerOD > 1
			@defConstrRef inequalityPath[1:numDataPoints, 1:(max_rounds-1)]
		end

		# Set first LP objective
		@setObjective(m, Min, sum{ sqrt(numRides[i,j]/travelTimes[i,j]) * epsilon[i,j], i=nodes, j=pairs[i]})

		# Add path constraints
		println("**** Adding constraints ****")
		if !(startWithSimpleLP) || l > 1
			println("---Fetching paths")
			paths, numExpensiveTurns = reconstructMultiplePathsWithExpensiveTurns(new_sp.previous, srcs, dsts, old_nodes, new_sp.real_destinations, new_edge_isExpensive)
			println("---Updating path array")
		end
		# Update paths
		if !(startWithSimpleLP) || l > 1
			totalPaths, totalNumExpensiveTurns, numPaths = updatePaths(paths, numExpensiveTurns, totalPaths, totalNumExpensiveTurns, numPaths, maxNumPathsPerOD=maxNumPathsPerOD, times=newTimes, turnCost=turnCost, travelTimes=travelTimes, numRides=numRides, srcs=srcs, dsts=dsts, dynamicConstraints=dynamicConstraints, globalConstraintUpdate=globalConstraintUpdate)
		else
			for i = 1:numDataPoints
				numPaths[i] = length(totalPaths[i])
			end
		end
		for i=1:numDataPoints
			# Add inequality constraints
			if numPaths[i] > 1
				for j = 1:(numPaths[i]-1)
					if split(model_type, "_")[1] == "strict"
						inequalityPath[i,j] = @addConstraint(m, sum{t[totalPaths[i][j+1][a],totalPaths[i][j+1][a+1]], a=1:(length(totalPaths[i][j+1])-1)} - sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} >= - tc * totalNumExpensiveTurns[i][j+1] + tc * totalNumExpensiveTurns[i][1])
					elseif split(model_type, "_")[1] == "relaxed"
						inequalityPath[i,j] = @addConstraint(m, sum{t[totalPaths[i][j+1][a],totalPaths[i][j+1][a+1]], a=1:(length(totalPaths[i][j+1])-1)} - sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} >= - tc * totalNumExpensiveTurns[i][j+1] + tc * totalNumExpensiveTurns[i][1] - delta * travelTimes[srcs[i],dsts[i]])
					end
				end
			end
		end
		# Equality constraints (shortest path close to travel time data)
		@addConstraint(m, TLowerBound[i=1:numDataPoints], sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} + tc * totalNumExpensiveTurns[i][1] - travelTimes[srcs[i],dsts[i]] >= - epsilon[srcs[i],dsts[i]])
		@addConstraint(m, TUpperBound[i=1:numDataPoints], sum{t[totalPaths[i][1][a],totalPaths[i][1][a+1]], a=1:(length(totalPaths[i][1])-1)} + tc * totalNumExpensiveTurns[i][1] - travelTimes[srcs[i],dsts[i]] <= epsilon[srcs[i],dsts[i]])

		# Define objective variables and constraints for second part of m
		if split(model_type, "_")[2] == "smooth" || (last_smooth && l == max_rounds)
			@defVar(m, delta2[i=nodes,j=out[i]] >= 0)
			@addConstraint(m, objConstrLower[i=nodes,j=out[i]], -1 * t[i,j]/distances[i,j] + 1/(length(inn[i]) + length(out[j])) * (sum{1/distances[j,k] * t[j,k], k = out[j]} + sum{1/distances[h,i] * t[h,i], h=inn[i]}) <= delta2[i,j])
			@addConstraint(m, objConstrUpper[i=nodes,j=out[i]], t[i,j]/distances[i,j] - 1/(length(inn[i]) + length(out[j])) * (sum{1/distances[j,k] * t[j,k], k = out[j]} + sum{1/distances[h,i] * t[h,i], h=inn[i]}) <= delta2[i,j])
		elseif split(model_type, "_")[2] == "random"
			coeffVec = 2 * rand(length(roads)) - 1
			coeffs = spzeros(length(nodes), length(nodes))
			k = 1
			for i = nodes, j = out[i]
				coeffs[i,j] = coeffVec[k]
				k += 1
			end
			if turnCostAsVariable
				coeff_tc = 2 * rand() - 1
			end
		end
		# Solve LP
		println("**** Solving LP ****")
		status = solve(m)
		if status == :Infeasible
			buildInternalModel(m)
			print_iis_gurobi(m)
			break
		end
		result = getValue(t)
		objective = getObjectiveValue(m)
		# Update cost function log
		write(costFunctionFile, string(l, ",", objective, "\n"))
		flush(costFunctionFile)
		if turnCostAsVariable
			println("Left turn cost: ", getValue(tc))
		end

		println("**** Setting up second LP ****")
		if split(model_type, "_")[2] == "smooth" || (last_smooth && l == max_rounds)
			if turnCostAsVariable
				@setObjective(m, Min, sum{delta2[i,j], i=nodes, j=out[i]} + tc)
			else
				@setObjective(m, Min, sum{delta2[i,j], i=nodes, j=out[i]})
			end
		elseif split(model_type, "_")[2] == "random"
			# Need to minimize absolute value (otherwise unbounded) so define new variable
			@defVar(m, zAbsVal >= 0)
			# Depending on whether turn cost is a variable, add it (or not) to the objective
			if turnCostAsVariable
				@addConstraint(m, randomObjectiveConstraint1, sum{coeffs[i,j] * t[i,j]/distances[i,j], i=nodes, j=out[i]} + coeff_tc * tc <= zAbsVal)
				@addConstraint(m, randomObjectiveConstraint2, sum{coeffs[i,j] * t[i,j]/distances[i,j], i=nodes, j=out[i]} + coeff_tc * tc >= - zAbsVal)
			else
				@addConstraint(m, randomObjectiveConstraint1, sum{coeffs[i,j] * t[i,j]/distances[i,j], i=nodes, j=out[i]} <= zAbsVal)
				@addConstraint(m, randomObjectiveConstraint2, sum{coeffs[i,j] * t[i,j]/distances[i,j], i=nodes, j=out[i]} >= - zAbsVal)
			end
			@setObjective(m, Min, zAbsVal)
		end

		if split(model_type, "_")[2] == "nothing" && !(last_smooth && l == max_rounds)
			# do nothing
			if turnCostAsVariable
				turnCost = getValue(tc)
				write(turnCostFile, string(l,",", turnCost, "\n"))
			end
		else
			# Solve second LP
			println("**** Solving second LP ****")
			@addConstraint(m, fixObjective, sum{sqrt(numRides[i,j]/travelTimes[i,j]) * epsilon[i,j], i=nodes, j=pairs[i]} <=objective + 1e-1)
			status = solve(m)
			# Get new turning cost
			if turnCostAsVariable
				turnCost = getValue(tc)
				println(turnCost)
				write(turnCostFile, string(l,",", turnCost, "\n"))
			end
		end

		# Debug if infeasible
		if status == :Infeasible
			println("!!!! Computing IIS !!!!")
			buildInternalModel(m)
			print_iis_gurobi(m)
			break
		# Prepare output
		elseif status == :Optimal || status == :Suboptimal || status == :UserLimit
			if status != :UserLimit
				result = getValue(t)
			end
			newTimes = spzeros(length(nodes), length(nodes))
			for key in keys(result)
				newTimes[key[1], key[2]] = result[key[1], key[2]]
			end
			# Save updated Manhattan road times to file
			if metropolis
				saveRoadTimes(newTimes, "$TESTDIR/metropolis-times-$l")
			else
				saveRoadTimes(newTimes, "$TESTDIR/manhattan-times-$l")
			end
			if abs(objective - old_objective)/old_objective < 1e-6
				save("Outputs/$TESTDIR/end.jld", "num_iter", l)
				l = max_rounds
			else
				old_objective = objective
			end
		end
		if l < max_rounds || computeFinalSP
			println("**** Computing shortest paths ****")
			new_graph, new_edge_dists, new_nodes, new_edge_isExpensive = modifyGraphForDijkstra(graph, newTimes, manhattan.positions, turn_cost=turnCost)
			old_nodes = getInverseMapping(new_nodes, nv(new_graph))
			@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
			if metropolis
				printErrorStatsToFile("Outputs/$TESTDIR/errorstats-$(actual_l).txt", real_TOD_metropolis, travelTimes, new_sp.traveltime, numRides, real_tij_metropolis, newTimes, num_nodes = nv(graph))
			else
				if errorComputation == "single-ride"
					@time avg_sq_err, avg_rel_err, avg_bias = computeTestingError(testingData, nodeTree, nodePairs, new_sp.traveltime, method)
					write(errorFile, string(l, ",", avg_sq_err, ",", avg_rel_err, ",", avg_bias,"\n"))
					flush(errorFile)
				elseif errorComputation == "average"
					@time avg_sq_err, avg_rel_err, avg_bias = computeAverageTestingError(testingData, new_sp.traveltime, num_nodes=nv(graph))
					write(errorFile, string(l, ",", avg_sq_err, ",", avg_rel_err, ",", avg_bias,"\n"))
					flush(errorFile)
				elseif errorComputation == "both"
					@time avg_sq_err, avg_rel_err, avg_bias = computeTestingError(testingData, nodeTree, nodePairs, new_sp.traveltime, method)
					@time avg_sq_err_2, avg_rel_err_2, avg_bias_2 = computeAverageTestingError(testingData2, new_sp.traveltime, num_nodes=nv(graph))
					write(errorFile, string(l, ",", avg_sq_err, ",", avg_rel_err, ",", avg_bias, ",", avg_sq_err_2, ",", avg_rel_err_2, ",", avg_bias_2,"\n"))
					flush(errorFile)
				end
			end
			if dynamicConstraints
				write(numConstraintsFile, string(l,",",length(srcs),",", count_number_edges_used(totalPaths, manhattan),"\n"))
				flush(numConstraintsFile)
				if l < max_rounds && l % iterationMultiple == 0
					println("**** Updating constraint set ****")
					srcs, dsts, totalPaths, totalNumExpensiveTurns, numPaths, pairs = updateConstraints(travelTimes, numRides, new_sp.traveltime, totalPaths, totalNumExpensiveTurns, numPaths, srcs, dsts, pairs, numNodePairsToAdd = numPairsToAdd, numNodePairsToRemove = numPairsToRemove, addOnlyIfAbovePercentile = addOnlyIfAbovePercentile, removeOnlyIfBelowPercentile = removeOnlyIfBelowPercentile)
					numDataPoints = length(srcs)
				end
			end
		end
		write(timeFile, string(l, ",", toc(), "\n"))
		flush(timeFile)
		l += 1
	end
	println("-- Saved outputs to directory Outputs/$(TESTDIR)/")
	# Close log files
	if turnCostAsVariable
		close(turnCostFile)
	end
	if !metropolis
		close(errorFile)
	end
	close(costFunctionFile)
	close(timeFile)
	if dynamicConstraints
		close(numConstraintsFile)
	end
	return status, newTimes, "Outputs/$TESTDIR/"
end

# manhattan = loadCityGraph()
# travel_times, num_rides, trainTestDf, testing_travel_times, testing_num_rides, testDf = loadInputTravelTimes(manhattan.positions, METHOD, year=YEAR, startTime=START_TIME, endTime=END_TIME, startMonth=START_MONTH, endMonth=END_MONTH, radius = RADIUS)
# if RANDOM_CONSTRAINTS
# 	travel_times, num_rides = chooseConstraints(travel_times, num_rides, sample_size=SAMPLE_SIZE);
# end
# if ERROR_COMPUTATION == "single-ride"
# 	@time status, new_times = fast_LP(manhattan, travel_times, num_rides, trainTestDf, manhattan.roadTime)
# elseif ERROR_COMPUTATION == "average"
# 	@time status, new_times = fast_LP(manhattan, travel_times, num_rides, testing_travel_times, manhattan.roadTime)
# elseif ERROR_COMPUTATION == "both"
# 	@time status, new_times = fast_LP(manhattan, travel_times, num_rides, trainTestDf, manhattan.roadTime, testingData2 = testing_travel_times)
# end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
# functions that populate the initial conditions
# List of functions:

"""
  This macro facilitates make initial condition functions.  Specifically,
  IC functions that use a function in common_funcs.jl to populate the initial
  condition can be created with this macro

  **Inputs**

   * namestub: the name stub.  If the function from common_funcs.jl is named
               `calcFoo`, then `namestub` should be `Foo` and the new initial
               condition function will be `ICFoo`.
   * docstring: docstring for the newly created function `ICFoo`.  This argument
                is optional, but good programming practices dictate that it
                should be provided in all cases.

  **Example Usage**

   ```
     @makeIC Foo \"\"\"
Docstring for ICFoo
\"\"\"
  ```
"""
macro makeIC(namestub, docstring="""""")

  # names of IC function, common_funcs.jl function
  fname = Symbol(string("IC", namestub))
  calcname = Symbol(string("calc", namestub))

  # introducing an empty docstring is not the same as not supplying a docstring
  # figure out which to do
  if docstring != ""
    docex = quote
              @doc $docstring $fname
            end
  else
    docex = :()
  end

  return esc(
    quote

      function $fname(mesh::AbstractMesh{Tmsh}, sbp, eqn::EulerData{Tsol}, opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsol,}

        q = zeros(Tsol, mesh.numDofPerNode)
        for i=1:mesh.numEl
          for j=1:mesh.numNodesPerElement
            dofs = sview(mesh.dofs, :, j, i)
            coords = sview(mesh.coords, :, j, i)
            $calcname(eqn.params, coords, q)
            for k=1:mesh.numDofPerNode
              u0[dofs[k]] = q[k]
            end
          end
        end

        return nothing
      end

    end  # end quote
  )

end

@doc """
### EulerEquationMod.ICZero

  Sets all components of the solution to zero
"""->
function ICZero(mesh::AbstractMesh{Tmsh},
operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol}, opts,
u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
# populate u0 with initial values
# this is a template for all other initial conditions

numEl = mesh.numEl
nnodes = operator.numnodes
dofpernode = mesh.numDofPerNode
for i=1:numEl

  for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      # get dof numbers for each variable
      dofnum_rho = dofnums_j[1]
      dofnum_rhou = dofnums_j[2]
      dofnum_rhov = dofnums_j[3]
      dofnum_e = dofnums_j[4]

      x = coords_j[1]
      y = coords_j[2]

      # apply initial conditions here
      u0[dofnum_rho] = 0.0
      u0[dofnum_rhou] = 0.0
      u0[dofnum_rhov] = 0.0
      u0[dofnum_e] = 0.0
  end
end

return nothing

end  # end function

@doc """
### EulerEquationMod.ICOnes

  Sets all components of the solution to 1.0
"""->

function ICOnes(mesh::AbstractMesh{Tmsh},
operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol}, opts,
u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}

  numEl = mesh.numEl
  nnodes = operator.numnodes
  dofpernode = mesh.numDofPerNode
  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      # get dof numbers for each variable
      dofnum_rho = dofnums_j[1]
      dofnum_rhou = dofnums_j[2]
      dofnum_rhov = dofnums_j[3]
      dofnum_e = dofnums_j[4]

      x = coords_j[1]
      y = coords_j[2]

      # apply initial conditions here
      u0[dofnum_rho] = 1.0
      u0[dofnum_rhou] = 1.0
      u0[dofnum_rhov] = 1.0
      u0[dofnum_e] = 1.0

      # u0 = 2*u0
    end
  end

  return nothing
end # end function ICOnes

@doc """
### EulerEquationMod.ICRho1E2

  Sets all density values to 1.0 and energy values to 2.0
"""
function ICRho1E2(mesh::AbstractMesh{Tmsh},
operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol}, opts,
u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
# populate u0 with initial values
# this is a template for all other initial conditions

numEl = mesh.numEl
nnodes = operator.numnodes
dofpernode = mesh.numDofPerNode
for i=1:numEl
  for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      # get dof numbers for each variable
      dofnum_rho = dofnums_j[1]
      dofnum_rhou = dofnums_j[2]
      dofnum_rhov = dofnums_j[3]
      dofnum_e = dofnums_j[4]

      x = coords_j[1]
      y = coords_j[2]

      # apply initial conditions here
      u0[dofnum_rho] = 1.0
      u0[dofnum_rhou] = 0.0
      u0[dofnum_rhov] = 0.0
      u0[dofnum_e] = 2.0
  end
end

return nothing

end  # end function

@makeIC Rho1E2U1VW0 """
### EulerEquationMod.ICRho1E2U1VW0

  Sets the density values 1.0, x momentum to 1.0,
  v & w momenta to 0.0, and energy to 2.0 at a node.

  It should work for 2D and 3D meshes.
"""

@makeIC Rho1E2U3 """
### EulerEquationMod.ICRho1E2U3

  Sets all components density values to 1.0, x and y momenta to 0.35355, and
  energy to 2.0
"""

@makeIC FreeStream """
### EulerEquationMod.ICFreeStream

  Sets all components of the solution to the free stream condition according
  to the angle of attack and and Mach number.

"""

@makeIC FreeStream0 """
  Like [`calcFreeStream`](@ref), but uses [`calcFreeStream0`](@ref) instead of
  [`calcFreeStream`](@ref).
"""

@makeIC Vortex """
 what is this? how is it different than ICIsentropic Vortex?
"""

@doc """
### EulerEquationMod.ICsmoothHeavisideder

  Sets the density to the derivative of the smooth Heaviside function, all
  other components to zero.
"""->
function ICsmoothHeavisideder(mesh::AbstractMesh{Tmsh},
            operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol},
            opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
# calculate the value of the smooth heaviside function derivative at a location x
# x0 is specified within this function

# smooth heaviside  parameters
  x0 = 0
  L = 5
  k = 5

  numEl = mesh.numEl
  nnodes = operator.numnodes
  dofpernode = mesh.numDofPerNode
  for i=1:numEl
  #  dofnums_i = sview(mesh, i)  # get dof nums for this element
  #  coords = sview(mesh, [i])

    for j=1:nnodes
        coords_j = sview(mesh.coords, :, j, i)
        dofnums_j = sview(mesh.dofs, :, j, i)

        # get dof numbers for each variable
        dofnum_rho = dofnums_j[1]
        dofnum_rhou = dofnums_j[2]
        dofnum_rhov = dofnums_j[3]
        dofnum_e = dofnums_j[4]

        x = coords_j[1]
        y = coords_j[2]

        # apply initial conditions here
        u0[dofnum_rho] = L*(2*k*e^(-2*k*x))/(e^(-2*k*x) +1 )^2
        u0[dofnum_rhou] = 0.0
        u0[dofnum_rhov] = 0.0
        u0[dofnum_e] = 0.0
    end
  end

  return nothing
end

@doc """
### EulerEquationMod.ICsmoothHeaviside

  Sets the density to the smooth Heaviside function, all other components to
  zero.

"""->
function ICsmoothHeaviside(mesh::AbstractMesh{Tmsh},
                           operator::AbstractOperator{Tsbp},
                           eqn::EulerData{Tsol},
                           opts, u0::AbstractArray{Tsol, 1}) where {Tmsh, Tsbp, Tsol}
  # calculate the value of the smooth heaviside function at a location x
  # x0 is specified within this function

  # smooth heaviside  parameters
  x0 = 0
  L = 5
  k = 5

  numEl = mesh.numEl
  nnodes = operator.numnodes
  dofpernode = mesh.numDofPerNode
  for i=1:numEl
    for j=1:nnodes
      coords = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)

      # get dof numbers for each variable
      dofnum_rho = dofnums_j[1]
      dofnum_rhou = dofnums_j[2]
      dofnum_rhov = dofnums_j[3]
      dofnum_e = dofnums_j[4]

      x = coords_j[1]
      y = coords_j[2]

      # apply initial conditions here
      u0[dofnum_rho] = L/(1 + e^(-k*(x-x0)))
      u0[dofnum_rhou] = 0.0
      u0[dofnum_rhov] = 0.0
      u0[dofnum_e] = 0.0
    end
  end

  return nothing
end

@makeIC IsentropicVortex """

### EulerEquationMod.ICIsentropicVortex

Sets the solution to the steady isentropic vortex solution.

"""

@doc """
### EulerEquationMod.ICIsentropicVortexWithNoise

  Sets the solution to the steady isentropic vortex solution plus
  a small random noise component.

"""->
function ICIsentropicVortexWithNoise(mesh::AbstractMesh{Tmsh},
                   operator::AbstractOperator{Tsbp},
                   eqn::EulerData{Tsol},
                   opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values
  # this is a template for all other initial conditions

  numEl = mesh.numEl
  nnodes = operator.numnodes
  dofpernode = mesh.numDofPerNode
  sol = zeros(Tsol, 4)
  for i=1:numEl
    for j=1:nnodes
      coords_j = sview(mesh.coords, :, j, i)
      dofnums_j = sview(mesh.dofs, :, j, i)
      calcIsentropicVortex(eqn.params, coords_j, sol)

      # apply initial conditions here
      for k=1:dofpernode
        u0[dofnums_j[k]] = sol[k] + 0.1*rand()
      end

    end
  end

  return nothing

end  # end function

@makeIC UnsteadyVortex """
### EulerEquationMod.ICUnsteadyVortex

  Sets the solution to the unsteady vortex problem.  eqn.params.t is used to
  determine what time to use for the solution.

"""
@makeIC UnsteadyVortex2 """
  Vortex travelling at an angle
"""

@makeIC InvChannelIC """
  Initial condition for SU2 bump in inviscid channel case.  See also
  The subsonic inflow and subsonic outflow boundary conditions.
"""
function ICInvChannel(mesh::AbstractMesh{Tmsh},
        operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol},
        opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol}

  sol = zeros(Tsol, mesh.numDofPerNode)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
      coords_j = sview(mesh.coords, :, j, i)
      calcInvChannelIC(eqn.params, coords_j, sol)

      for k=1:mesh.numDofPerNode
        u0[dofnums_j[k]] = sol[k]
      end
    end
  end

  return nothing
end

@doc """
### EulerEquationMod.ICFile

  This function reads a vector from a file on disk and set the solution to it.
  The vector must contain the same number of entries as there are degrees of
  freedom in the mesh.

  This function is useful for things like warm-starting the solver.
  The file should be the output of `writeSolutionFiles`, and must be loaded on
  an *identical* mesh to the one used when saving.

"""->
function ICFile(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator{Tsbp}, eqn::EulerData{Tsol}, opts,
                u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}
  # populate u0 with initial values from a disk file
  # the file name comes from opts["ICfname"]

  readSolutionFiles(mesh, sbp, eqn, opts, opts["ICfname"], u0)
end

"""
  This initial condition uses the current `eqn.q_vec` as the initial condition.
  This is useful for repeatedly solving a PDE with, using the previous solution
  as the initial guess.
"""
function ICPassThrough(mesh::AbstractMesh{Tmsh},
                       operator::AbstractOperator{Tsbp}, eqn::EulerData{Tsol},
                       opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsbp, Tsol}

  # use current eqn.q_vec as IC
  copy!(u0, eqn.q_vec)

  return nothing
end

@makeIC Exp """
  Assigns exp(k*x*y*z) as the initial condition, of each node, where k is
  the index of the degree of freedom of the node
"""

@makeIC PeriodicMMS """
  Writes calcPeriodicMMS to the initial condition vector u0
"""

"""
  This function applies the initial condition for the Taylor Green vortex,
  using the constants in Gassner, Winters, and Kopriva's Split form Nodal
  DG paper
"""
function ICTaylorGreen(mesh::AbstractMesh{Tmsh}, sbp,
          eqn::EulerData{Tsol}, opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsol,}

  # parameters
  M = 1.0  # Mach number
  gamma_1 = eqn.params.gamma_1
  gamma = eqn.params.gamma

  q = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      dofs = sview(mesh.dofs, :, j, i)
      coords = sview(mesh.coords, :, j, i)
      x = coords[1]
      y = coords[2]
      z = coords[3]

      p = 100/gamma + (1/16)*( cos(2*x)*cos(2*z) + 2*cos(2*y) + 2*cos(2*x) + cos(2*y)*cos(2*z) )
      q[1] = 1
      q[2] = q[1]*M*sin(x)*cos(y)*cos(z)
      q[3] = -q[1]*M*cos(x)*sin(y)*cos(z)
      q[4] = 0
      q[5] = p/gamma_1 + (q[2]*q[2] + q[3]*q[3] + q[4]*q[4])/(2*q[1])

      for k=1:mesh.numDofPerNode
        u0[dofs[k]] = q[k]
      end
    end
  end

  return nothing
end

@makeIC ChannelMMS """
  Initial condition of channel MMS
"""

@makeIC Square1D """
  Initial for square wave in 1D
"""

@makeIC Square2D """
  Initial for square wave in 2D
"""

@makeIC SedovExplosion """
  Initial for Sedov Explosion
"""

@makeIC LaplaceSolution """
  IC for testing shock capturing diffusion terms
"""

@makeIC Wedge20

@makeIC Wedge20Potential

# declare a const dictionary here that maps strings to function (used for input arguments)
"""
  Map IC names to functions.  Generally the name is the same as the function
  name

  Every function in this dictionary must have the signature:

  ```
    Foo(mesh::AbstractMesh{Tmsh}, sbp,
          eqn::EulerData{Tsol}, opts, u0::AbstractVector{Tsol}) where {Tmsh, Tsol,}
  ```

  where `u0` is the vector to be overwritten with the initial condition,
  `mesh` is the mesh object`, `sbp` is the SBP operator, `eqn` is the equation
  object, and `opts` is the optios dictionary.

  Many of the initial condition functions call a function in common_funcs.jl
  to compute the state.  See the macro [`makeIC`](@ref) for an easy way
  to make such functions.  If these functions depend on time, they use
  `eqn.params.t` to get the current time value.  This means IC functions can
  be used to get analytical solutions at times other than `t = 0`.

"""
global const ICDict = Dict{Any, Function}(
"ICZero" => ICZero,
"ICOnes" => ICOnes,
"ICRho1E2" => ICRho1E2,
"ICRho1E2U1VW0" => ICRho1E2U1VW0,
"ICRho1E2U3" => ICRho1E2U3,
"ICFreeStream" => ICFreeStream,
"ICFreeStream0" => ICFreeStream0,
"ICVortex" => ICVortex,
#"ICLinear" => ICLinear,
"ICsmoothHeavisideder" => ICsmoothHeavisideder,
"ICsmoothHeaviside" => ICsmoothHeaviside,
"ICIsentropicVortex" => ICIsentropicVortex,
"ICUnsteadyVortex" => ICUnsteadyVortex,
"ICUnsteadyVortex2" => ICUnsteadyVortex2,
"ICIsentropicVortexWithNoise" => ICIsentropicVortexWithNoise,
"ICInvChannel" => ICInvChannelIC,
"ICFile" => ICFile,
"ICPassThrough" => ICPassThrough,
"ICExp" => ICExp,
"ICPeriodicMMS" => ICPeriodicMMS,
"ICTaylorGreen" => ICTaylorGreen,
"ICChannelMMS" => ICChannelMMS,
"ICSquare1D" => ICSquare1D,
"ICSquare2D" => ICSquare2D,
"ICSedovExplosion" => ICSedovExplosion,
"ICLaplace" => ICLaplaceSolution,
"ICWedge20" => ICWedge20,
"ICWedge20Potential" => ICWedge20Potential,
)


text:
using FIX
using Test

import DataStructures: OrderedDict

@test checksum("123ioн") == 171

sample = OrderedDict(1 => "g", 2=> "h")

@test join(sample, Char(1)) == "1=g" * Char(1) * "2=h" * Char(1)

@test convert(TAGS_INT_STRING, sample) == Dict("AdvId" => "h", "Account" => "g")


text:
doc"""
    InverseGamma(α, θ)

The *inverse gamma distribution* with shape parameter `α` and scale `θ` has probability
density function

$f(x; \alpha, \theta) = \frac{\theta^\alpha x^{-(\alpha + 1)}}{\Gamma(\alpha)}
e^{-\frac{\theta}{x}}, \quad x > 0$

It is related to the [`Gamma`](:func:`Gamma`) distribution: if $X \sim \operatorname{Gamma}(\alpha, \beta)$, then $1 / X \sim \operatorname{InverseGamma}(\alpha, \beta^{-1})$.

.. code-block:: julia

    InverseGamma()        # Inverse Gamma distribution with unit shape and unit scale, i.e. InverseGamma(1, 1)
    InverseGamma(a)       # Inverse Gamma distribution with shape a and unit scale, i.e. InverseGamma(a, 1)
    InverseGamma(a, b)    # Inverse Gamma distribution with shape a and scale b

    params(d)        # Get the parameters, i.e. (a, b)
    shape(d)         # Get the shape parameter, i.e. a
    scale(d)         # Get the scale parameter, i.e. b

External links

* [Inverse gamma distribution on Wikipedia](http://en.wikipedia.org/wiki/Inverse-gamma_distribution)

"""
immutable InverseGamma{T<:Real} <: ContinuousUnivariateDistribution
    invd::Gamma{T}
    θ::T

    function InverseGamma(α, θ)
        @check_args(InverseGamma, α > zero(α) && θ > zero(θ))
        new(Gamma(α, 1 / θ), θ)
    end
end

InverseGamma{T<:Real}(α::T, θ::T) = InverseGamma{T}(α, θ)
InverseGamma(α::Real, θ::Real) = InverseGamma(promote(α, θ)...)
InverseGamma(α::Integer, θ::Integer) = InverseGamma(Float64(α), Float64(θ))
InverseGamma(α::Real) = InverseGamma(α, 1.0)
InverseGamma() = InverseGamma(1.0, 1.0)

@distr_support InverseGamma 0.0 Inf

#### Conversions
convert{T <: Real, S <: Real}(::Type{InverseGamma{T}}, α::S, θ::S) = InverseGamma(T(α), T(θ))
convert{T <: Real, S <: Real}(::Type{InverseGamma{T}}, d::InverseGamma{S}) = InverseGamma(T(shape(d.invd)), T(d.θ))

#### Parameters

shape(d::InverseGamma) = shape(d.invd)
scale(d::InverseGamma) = d.θ
rate(d::InverseGamma) = scale(d.invd)

params(d::InverseGamma) = (shape(d), scale(d))
@inline partype{T<:Real}(d::InverseGamma{T}) = T

#### Parameters

mean{T<:Real}(d::InverseGamma{T}) = ((α, θ) = params(d); α  > 1 ? θ / (α - 1) : T(Inf))

mode(d::InverseGamma) = scale(d) / (shape(d) + 1)

function var{T<:Real}(d::InverseGamma{T})
    (α, θ) = params(d)
    α > 2 ? θ^2 / ((α - 1)^2 * (α - 2)) : T(Inf)
end

function skewness{T<:Real}(d::InverseGamma{T})
    α = shape(d)
    α > 3 ? 4sqrt(α - 2) / (α - 3) : T(NaN)
end

function kurtosis{T<:Real}(d::InverseGamma{T})
    α = shape(d)
    α > 4 ? (30α - 66) / ((α - 3) * (α - 4)) : T(NaN)
end

function entropy(d::InverseGamma)
    (α, θ) = params(d)
    α + lgamma(α) - (1 + α) * digamma(α) + log(θ)
end

#### Evaluation

pdf(d::InverseGamma, x::Real) = exp(logpdf(d, x))

function logpdf(d::InverseGamma, x::Real)
    (α, θ) = params(d)
    α * log(θ) - lgamma(α) - (α + 1) * log(x) - θ / x
end

cdf(d::InverseGamma, x::Real) = ccdf(d.invd, 1 / x)
ccdf(d::InverseGamma, x::Real) = cdf(d.invd, 1 / x)
logcdf(d::InverseGamma, x::Real) = logccdf(d.invd, 1 / x)
logccdf(d::InverseGamma, x::Real) = logcdf(d.invd, 1 / x)

quantile(d::InverseGamma, p::Real) = 1 / cquantile(d.invd, p)
cquantile(d::InverseGamma, p::Real) = 1 / quantile(d.invd, p)
invlogcdf(d::InverseGamma, p::Real) = 1 / invlogccdf(d.invd, p)
invlogccdf(d::InverseGamma, p::Real) = 1 / invlogcdf(d.invd, p)

function mgf{T<:Real}(d::InverseGamma{T}, t::Real)
    (a, b) = params(d)
    t == zero(t) ? one(T) : 2(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*b*t))
end

function cf{T<:Real}(d::InverseGamma{T}, t::Real)
    (a, b) = params(d)
    t == zero(t) ? one(T)+zero(T)*im : 2(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*im*b*t))
end

#### Evaluation

rand(d::InverseGamma) = 1 / rand(d.invd)

function _rand!(d::InverseGamma, A::AbstractArray)
    s = sampler(d.invd)
    for i = 1:length(A)
    	v = 1 / rand(s)
        @inbounds A[i] = v
    end
    A
end


text:
module InterVal

import Base: convert, promote_type, show, showcompact,
    read, write, string, hash, sizeof,
    issubnormal, isfinite, isreal, isinf, isnan,
    zero, one, Inf, NaN,
    (==),(!=),(<),(<=),(>=),(>), isless, isequal,
    (+),(-),(*),(/),sqrt

export ClCl, ClOp, OpCl, OpOp,
       ClCl!, ClOp!, OpCl!, OpOp!,
       lowerbound, upperbound

include("jas/ErrorFree.jl/src/ErrorFree.jl")
using .ErrorFree
include("jas/Adjacent.jl/src/Adjacent.jl")
using .Adjacent
include("jas/RoundFast.jl/src/RoundFast.jl")
using .RoundFast

include("type/interval.jl")
include("type/easily.jl")
include("type/bounding.jl")
include("type/io.jl")

include("number/compares.jl")
include("number/givenvals.jl")
include("number/predicates.jl")

include("arith/add.jl")
include("arith/sub.jl")
include("arith/mul.jl")
include("arith/div.jl")
include("arith/sqrt.jl")
include("arith/shared.jl")

end # module


text:
module Poisson2D

using LinearAlgebra
using SparseArrays

include("structures.jl")
export FDPoisson2DProblem, DirichletBC, NeumannBC, PeriodicBC
include("constructors_sparse.jl")
export FDPoisson2DProblem
include("assembly_fd.jl")
export assemble_system!
include("solve_fd.jl")
export direct_solve_poisson2d

end # module


text:
NumTestRepetitions = 100

@testset "Adaptive differential evolution optimizer" begin

ss = symmetric_search_space(1, (0.0, 10.0))
fake_problem = FunctionBasedProblem(x -> 0.0, "test_problem", MinimizingFitnessScheme, ss)

ade = adaptive_de_rand_1_bin(fake_problem, ParamsDict(
    :Population => rand(1, 100)))

@testset "parameters adjust!()" begin
    for i in 1:NumTestRepetitions
        cur_cr, cur_f = BlackBoxOptim.crossover_parameters(ade.modify, ade.population, rand(1:popsize(ade.population)))
        @test (0.0 <= cur_cr <= 1.0)
        @test (0.0 <= cur_f <= 1.0)
        BlackBoxOptim.adjust!(ade.modify, 0, 1, 0.0, 0.0, false)
        # FIXME this fails too often; needs adjusting the distribution?
        #new_cr, new_f = BlackBoxOptim.crossover_parameters(ade.params, 1)
        #@fact new_cr != cur_cr --> true
        #@fact new_f != cur_f --> true
    end
end

# FIXME actually this test is not required as the standard DE already tests for that
@testset "ask()" begin
    for i in 1:NumTestRepetitions
        res = BlackBoxOptim.ask(ade)
        @test length(res) == 2

        trial, target = res

        @test ndims(trial.params) == 1
        @test (1 <= trial.index <= popsize(ade))
        @test in(trial.params, ade.embed.searchSpace)

        @test ndims(target.params) == 1
        @test (1 <= target.index <= popsize(ade))
        @test in(target.params, ade.embed.searchSpace)

        @test trial.index == target.index
    end
end

end


text:
using DemoTDD
using Test

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

    @test 3 == 3
	# Super cool y nuevo proceso:
	#@test T == T broken=(T = NaN)
    #
    @test_broken 1 == 2

    @testset "Nivel 1" begin
        @test 1 == 1
        @testset "Nivel 2" begin
            @test 2 == 2
            @testset "Nivel 3" begin
                @test 3 == 2
            end
        end
    end
    for i in (1,2,3)
        @test i == i
    end

    @testset "Nivel 1 otra vez" begin
        @test true
        @test false
        @test true
    end

    @test all(true for i in 1:100000)

    @test_throws BoundsError [1, 2][3]

end


text:
function getAnalyticalConstGrad2D(n::Array{Int64,1},h::Array{Float64,1})
src = [1,div(n[2],2)];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];

(X1,X2) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

a  = -0.4;
s0 = 2.0;

kappaSquared = s0^2 + 2*a*(X1-source1);

# matshow(kappaSquared); colorbar();
SBarSquared = s0^2 + a*(X1-source1);
rSquared = ((X1-source1).^2 + (X2-source2).^2);
sigmaSquarred = 2*rSquared./(SBarSquared + sqrt.(SBarSquared.^2 - (a^2).*rSquared));
sigma = sqrt.(sigmaSquarred);
T_exact =  SBarSquared.*sigma - (a^2).*(sigma.^3)./6;

return kappaSquared,src,T_exact;
end

function getAnalyticalConstGrad3D(n::Array{Int64,1},h::Array{Float64,1})
src = [div(n[1],2),div(n[2],2),1];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];
source3 = (src[3]-1)*h[3];

(X1,X2,X3) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

# a  = -1.75;
a  = -1.65;
s0 = 2.0;

kappaSquared = s0^2 + 2*a*(X3-source3);
SBarSquared = s0^2 + a*(X3-source3);
rSquared = (X1-source1).^2 + (X2-source2).^2 + (X3-source3).^2;
sigmaSquarred = 2*rSquared./(SBarSquared + sqrt.(SBarSquared.^2 - (a^2).*rSquared));
sigma = sqrt.(sigmaSquarred);
T_exact =  SBarSquared.*sigma - (a^2).*(sigma.^3)./6;

return kappaSquared,src,T_exact;
end

function runAnalyticalConstGradExperiment3D(h::Array{Float64,1},WU::Float64)
I = [1.6,1.6,0.8];

n = zeros(Int64,3);
n[1] = round(Int64,(I[1]/h[1])+1);
n[2] = round(Int64,(I[2]/h[2])+1);
n[3] = round(Int64,(I[3]/h[3])+1);

# src = [div(n[1],2),div(n[2],2),div(n[3],2)];
src = [div(n[1],2),div(n[2],2),1];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];
source3 = (src[3]-1)*h[3];

(X1,X2,X3) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

a  = -1.75;
s0 = 2.0;

kappaSquared = s0^2 + 2*a*(X3-source3);
SBarSquared = s0^2 + a*(X3-source3);
rSquared = (X1-source1).^2 + (X2-source2).^2 + (X3-source3).^2;
sigmaSquarred = 2*rSquared./(SBarSquared + sqrt.(SBarSquared.^2 - (a^2).*rSquared));
sigma = sqrt.(sigmaSquarred);
T_exact =  SBarSquared.*sigma - (a^2).*(sigma.^3)./6;

sigmaSquarred = sigma = SBarSquared = X1 = X2 = X3 = rSquared = 0.0;

return T_exact,src,kappaSquared;
end

function getAnalyticalConstGradInv2D(n::Array{Int64,1},h::Array{Float64,1})

src = [1,div(n[2],2)];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];

(X1,X2) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

a  = 1.0;
s0 = 2.0;

Kappa = 1./(1./s0 + a*(X1-source1));
T_exact = (1/a)*acosh.(1+(0.5*s0*a*a).*Kappa.*((X1-source1).^2 + (X2-source2).^2));
G2_exact = (1/a)*(1./sqrt.((1+(0.5*s0*a*a).*Kappa.*((X1-source1).^2 + (X2-source2).^2)).^2 - 1)).*(0.5*s0*a*a).*Kappa.*2.*(X2-source2);
G1_exact = (1/a)*(1./sqrt.((1+(0.5*s0*a*a).*Kappa.*((X1-source1).^2 + (X2-source2).^2)).^2 - 1)).*((0.5*s0*a*a).*Kappa.*2.*(X1-source1)
			+ (0.5*s0*a*a).*((X1-source1).^2 + (X2-source2).^2).*(-(Kappa.^2)).*a);
kappaSquared = Kappa.^2;

return kappaSquared,src,T_exact;

end

function getAnalyticalConstGradInv3D(n::Array{Int64,1},h::Array{Float64,1})

src = [div(n[1],2),div(n[1],2),1];

source1 = (src[1]-1)*h[1];
source2 = (src[2]-1)*h[2];
source3 = (src[3]-1)*h[3];

(X1,X2,X3) = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

a  = 1.0;
s0 = 2.0;
Kappa = 1./(1./s0 + a*(X3-source3));
Rsquared = (X1-source1).^2 + (X2-source2).^2 + (X3-source3).^2;
T_exact = (1/a)*acosh.(1+(0.5*s0*a*a).*Kappa.*(Rsquared));
G1_exact = (1/a)*(1./sqrt.((1+(0.5*s0*a*a).*Kappa.*Rsquared).^2 - 1)).*(0.5*s0*a*a).*Kappa.*2.*(X1-source1);
G2_exact = (1/a)*(1./sqrt.((1+(0.5*s0*a*a).*Kappa.*Rsquared).^2 - 1)).*(0.5*s0*a*a).*Kappa.*2.*(X2-source2);
G3_exact = (1/a)*(1./sqrt.((1+(0.5*s0*a*a).*Kappa.*Rsquared).^2 - 1)).*((0.5*s0*a*a).*Kappa.*2.*(X3-source3)
			+ (0.5*s0*a*a).*(Rsquared).*(-(Kappa.^2)).*a);
kappaSquared = Kappa.^2;

return kappaSquared,src,T_exact;

end

function getSmoothGaussianMedium(n::Array{Int64,1},h::Array{Float64,1})
src = div.(n,4);
kappaSquared = [];
T_exact = [];
if length(n)==2
	(T1_exact,G11_exact,G12_exact) = getSmoothFactoredModel(n,h);
	(T0,G01,G02,L0) = getAnalytic2DeikonalSolutionAll(n,h,src);
	G1_exact = T0.*G11_exact + G01.*T1_exact;
	G2_exact = T0.*G12_exact + G02.*T1_exact;
	kappaSquared = G1_exact.*G1_exact + G2_exact.*G2_exact;
	# matshow(sqrt.(kappaSquared)); colorbar();
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],0:1:4);

	#figure()
	#CS = contour(kappaSquared,50);colorbar();title("kappaSquared");
	#clabel(CS, inline=1, fontsize=10)
	#figure()
	#contour(G1_exact);colorbar();title("G1_exact");
	T_exact = T0.*T1_exact;
	#figure()
	#contour(T_exact,100);colorbar();title("T exact");

else

	(T1_exact,G11_exact,G12_exact,G13_exact) = getSmoothFactoredModel3D(n,h);

	(T0,G01,G02,G03) = getAnalytic3DeikonalSolutionAll(n,h,src);

	G1_exact = T0.*G11_exact + G01.*T1_exact;
	G2_exact = T0.*G12_exact + G02.*T1_exact;
	G3_exact = T0.*G13_exact + G03.*T1_exact;
	kappaSquared = G1_exact.*G1_exact + G2_exact.*G2_exact + G3_exact.*G3_exact;

	# Er = kappaSquared;
	# for k=1:n[3]
		# matshow(reshape(Er[:,:,k],n[1],n[2]));colorbar();title("kappaSQ");
	# end

	T_exact = T0.*T1_exact;

end

return kappaSquared,src,T_exact;
end

function getSmoothFactoredModel(n::Array{Int64,1},h::Array{Float64,1})
xsrc = zeros(2);
src_kappa = zeros(Int64,2);
src_kappa[1] = div(n[1],3);
src_kappa[2] = div(n[2],4);

xsrc[1] = (src_kappa[1]-1)*h[1];
xsrc[2] = (src_kappa[2]-1)*h[2];

X1,X2 = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2]);

sigma = 0.1;
T1_exact = (exp.( - (sigma*((X1 - xsrc[1]).^2) + 4*sigma*((X2-xsrc[2]).^2))) + 1)/2;
G11_exact = -2*sigma*(X1 - xsrc[1]).*exp.( - (sigma*((X1 - xsrc[1]).^2) + 4*sigma*((X2-xsrc[2]).^2)))/2;
G12_exact = -8*sigma*(X2 - xsrc[2]).*exp.( - (sigma*((X1 - xsrc[1]).^2) + 4*sigma*((X2-xsrc[2]).^2)))/2;

return T1_exact,G11_exact,G12_exact;
end

function getSmoothFactoredModel3D(n::Array{Int64,1},h::Array{Float64,1})
xsrc = zeros(3);
src_kappa = zeros(Int64,3);
src_kappa[1] = div(n[1],3);
src_kappa[2] = div(n[2],4);
src_kappa[3] = div(n[3],2);

xsrc[1] = (src_kappa[1]-1)*h[1];
xsrc[2] = (src_kappa[2]-1)*h[2];
xsrc[3] = (src_kappa[3]-1)*h[3];

X1,X2,X3 = ndgrid((0:(n[1]-1))*h[1],(0:(n[2]-1))*h[2],(0:(n[3]-1))*h[3]);

Sigma = [0.1,0.4,0.2];
EXPRSkewedSquared = exp.(-Sigma[1]*((X1 - xsrc[1]).^2) - Sigma[2]*((X2-xsrc[2]).^2) - Sigma[3]*((X3-xsrc[3]).^2));

T1_exact = EXPRSkewedSquared/2 + 1/2;
G11_exact = -Sigma[1]*(X1 - xsrc[1]).*EXPRSkewedSquared;
G12_exact = -Sigma[2]*(X2 - xsrc[2]).*EXPRSkewedSquared;
G13_exact = -Sigma[3]*(X3 - xsrc[3]).*EXPRSkewedSquared;

return T1_exact,G11_exact,G12_exact,G13_exact;
end