text:
mutable struct RelaxedStateSequencingNode{T<:Real,U<:Int,V<:Int}<:Node
    uuid                    ::UUID
    length                  ::U
    state                   ::BitVector
    allUp                   ::BitVector
    allDown                 ::BitVector
    someUp                  ::BitVector
    someDown                ::BitVector
    lengthToRoot            ::T
    lengthToTerminal        ::T
    bestSeq                 ::Vector{V}
    exact                   ::Bool
    arcsIn                  ::Vector{SimpleArc}
    arcsOut                 ::Vector{SimpleArc}
end

#change hashing of nodes to use uuids
Base.hash(node::RelaxedStateSequencingNode) = hash(node.uuid)
Base.:(==)(node1::RelaxedStateSequencingNode, node2::RelaxedStateSequencingNode) = node1.uuid == node2.uuid
Base.length(node::RelaxedStateSequencingNode) = node.length
#do UUIDs automatically
function RelaxedStateSequencingNode(pathLength::U,state::BitVector,allUp::BitVector,allDown::BitVector,someUp::BitVector,someDown::BitVector,lengthToRoot::T,lengthToTerminal::T, bestSeq::Vector{V}, exact::Bool,arcsIn::Vector{SimpleArc},arcsOut::Vector{SimpleArc})where{T<:Real,U<:Int, V<:Int}
    return RelaxedStateSequencingNode(uuid4(), pathLength,state,allUp,allDown,someUp,someDown,lengthToRoot,lengthToTerminal, bestSeq, exact,arcsIn,arcsOut)
end

"""
    deepCopy(node::RelaxedStateSequencingNode)

Deep copy a node

# Arguments
- `node::RelaxedStateSequencingNode`: the node
"""
function deepCopy(node::RelaxedStateSequencingNode)
    return RelaxedStateSequencingNode(
        node.uuid, #uuid
        node.length, #path length
        copy(node.state), #current state
        copy(node.allUp), #allUp
        copy(node.allDown), #allDown
        copy(node.someUp), #someUp
        copy(node.someDown), #someDown
        node.lengthToRoot, #to root
        node.lengthToTerminal, #to terminal
        copy(node.bestSeq),
        node.exact, #exact
        deepCopy(node.arcsIn), #arcsin
        deepCopy(node.arcsOut) #arcsout
    )
end

Base.show(io::IO,node::RelaxedStateSequencingNode) = Base.print(
    io,
    "\n   pathLength: ", node.length, "\n",
    "   state: ", showBitVector(node.state), "\n",
    "   allUp: ", node.allUp, "\n",
    "   allDown: ", node.allDown, "\n",
    "   someUp: ", node.someUp, "\n",
    "   someDown: ", node.someDown, "\n",
    "   lengthToRoot: ", node.lengthToRoot, "\n",
    "   lengthToTerminal: ", node.lengthToTerminal, "\n",
    "   pathToNode: ", node.bestSeq, "\n",
    "   isExact: ", node.exact, "\n",
    "   arcsIn: ", showArcAssignments(node.arcsIn), "\n",
    "   arcsOut: ", showArcAssignments(node.arcsOut), "\n",
)

function showBitVector(b::BitVector)
    str = ""
    for (i,bit) in enumerate(b)
        if bit
            str = string(str, "$i ")
        end
    end
    return str
end

function showArcAssignments(arcs::Vector{SimpleArc})
    str = ""
    for arc in arcs
        str = string(str, arc.assignment, " ")
    end
    return str
end

"""
    createRootRelaxedStateSequencingNodeFromModel(model::SequencingModel)

Create a root node for a relaxed DD from the sequencing model

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
"""
function createRootRelaxedStateSequencingNodeFromModel(model::SequencingModel)
    restrictedNode = createRestrictedRootNode(model)

    state = falses(length(model))
    state[last(restrictedNode.path)] = true
    relaxedNode = RelaxedStateSequencingNode(
        length(restrictedNode.path), #length
        state, #current state
        falses(length(model)), #allUp
        copy(restrictedNode.visited), #allDown
        falses(length(model)), #someUp
        copy(restrictedNode.visited), #someDown
        0, #to root
        0, #to terminal
        copy(restrictedNode.path),
        true, #exact
        Vector{SimpleArc}(), #arcsin
        Vector{SimpleArc}() #arcsout
    )

    return relaxedNode, restrictedNode.domain
end

"""
createPartialRelaxedStateNode(;length::T=0, state::BitVector=falses(1), allUp::BitVector=BitVector(), allDown::BitVector=BitVector(), someUp::BitVector=BitVector(), someDown::BitVector=BitVector(), bestSeq::Vector{U}=Vector{Int}(),lengthToRoot::V=0, lengthToTerminal::V=0, exact::Bool=false, arcsIn::Vector{SimpleArc}= Vector{SimpleArc}(), arcsOut::Vector{SimpleArc}= Vector{SimpleArc}())where{T<:Int, U<:Int, V<:Real}

Create a node using default values for any combination of fields

# Arguments
- `length::T`: the length field
- `state::BitVector`: the state field
- `allUp::BitVector`: the allUp field
- `allDown::BitVector`: the allDown field
- `someUp::BitVector`: the someUp field
- `someDown::BitVector`: the someDown field
- `bestSeq::Vector{Int}`: the bestSeq field
- `lengthToRoot::V`: the lengthToRoot field
- `lengthToTerminal::V`: the lengthToTerminal field
- `exact::Bool`: the exact field
- `arcsIn::Vector{SimpleArc}`: the arcsIn field
- `arcsOut::Vector{SimpleArc}`: the arcsOut field
"""
function createPartialRelaxedStateNode(;length::T=0, state::BitVector=falses(1), allUp::BitVector=BitVector(), allDown::BitVector=BitVector(), someUp::BitVector=BitVector(), someDown::BitVector=BitVector(), bestSeq::Vector{U}=Vector{Int}(),lengthToRoot::V=0, lengthToTerminal::V=0, exact::Bool=false, arcsIn::Vector{SimpleArc}= Vector{SimpleArc}(), arcsOut::Vector{SimpleArc}= Vector{SimpleArc}())where{T<:Int, U<:Int, V<:Real}
    return RelaxedStateSequencingNode(
        length, #length
        state, #current state
        allUp, #allUp
        allDown, #allDown
        someUp, #someUp
        someDown, #someDown
        lengthToRoot, #to root
        lengthToTerminal, #to terminal
        bestSeq,
        exact, #exact
        arcsIn, #arcsin
        arcsOut #arcsout
    )
end

"""
    getParents(node::RelaxedStateSequencingNode)

Find and return the parents of a node

# Arguments
- `node::RelaxedStateSequencingNode`: the node whose parents to find
"""
function getParents(node::RelaxedStateSequencingNode)
    parents = Set{RelaxedStateSequencingNode}()
    for arc in node.arcsIn
        push!(parents, arc.origin)
    end
    return parents
end

"""
    getChildren(node::RelaxedStateSequencingNode)

Find and return the children
# Arguments
- `node::RelaxedSequencingNode`: the node whose children to find
"""
function getChildren(node::RelaxedStateSequencingNode)
    children = Set{RelaxedStateSequencingNode}()
    for arc in node.arcsOut
        push!(children, arc.destination)
    end
    return children
end

"""
    updateNodeDownVariables!(node::RelaxedStateSequencingNode,model::SequencingModel)

Update the state variable for allDown and someDown and lengthToRoot and exact

# Arguments
- `node::RelaxedStateSequencingNode`: the node to update
- `model::SequencingModel`: the sequencing problem being evaluated
"""
function updateNodeDownVariables!(node::RelaxedStateSequencingNode,model::SequencingModel)
    node.allDown = falses(length(model))
    node.someDown = falses(length(model))
    node.state = falses(length(model))
    node.exact = true

    for arc in node.arcsIn
        node.state[arc.assignment] = true
        node.someDown[arc.assignment] = true
        node.allDown[arc.assignment] = true
    end
    #update the best path variables
    bestArc, node.lengthToRoot = findBestArc(model, node.arcsIn)
    node.bestSeq = push!(copy(bestArc.origin.bestSeq), bestArc.assignment)
    #union domains with parents and check for exactness
    for parent in getParents(node)
        node.someDown = node.someDown .| parent.someDown
        node.allDown = node.allDown .| parent.allDown
        if !parent.exact
            node.exact = false
        end
    end
    #check each all down to make sure its valid
    for i in 1:length(node.allDown)
        if node.allDown[i]
            for arc in node.arcsIn
                if arc.assignment != i
                    if !arc.origin.allDown[i]
                        node.allDown[i]=false
                        break
                    end
                end
            end
        end
    end
    #final check for exactness (first part is just a check to see if the second part should run)
    if node.exact && (sum(node.state)!=1 || node.someDown != node.allDown || node.length!=sum(node.allDown))
        node.exact = false
    end
end

"""
    updateOutgoingArcValues!(model::SequencingModel, node::RelaxedStateSequencingNode)

Update the arcs leaving a node
# Arguments
- `model::SequencingModel`: the sequencing problem being evaluated
- `node::RelaxedStateSequencingNode`: the node to update
"""
function updateOutgoingArcValues!(model::SequencingModel, node::RelaxedStateSequencingNode)
    toFilter = Vector{SimpleArc}()
    for arc in node.arcsOut
        tempVal = getArcValue(model, arc)
        if isnothing(tempVal)
            push!(toFilter, arc)
        else
            arc.value = tempVal
        end
    end
    deleteArcs!(toFilter)
end

"""
    updateNodeUpVariables!(node::RelaxedStateSequencingNode,model::SequencingModel)

Update the state variable for someUp, allUp, and lengthToTerminal

# Arguments
- `node::RelaxedStateSequencingNode`: the node to update
- `model::SequencingModel`: the sequencing problem being evaluated
"""
function updateNodeUpVariables!(node::RelaxedStateSequencingNode,model::SequencingModel)
    node.allUp = falses(length(model))
    node.someUp = falses(length(model))
    if isempty(node.arcsOut)
        node.lengthToTerminal = 0
    else
        node.lengthToTerminal = node.arcsOut[1].value + node.arcsOut[1].destination.lengthToTerminal
    end
    for arc in node.arcsOut
        node.someUp[arc.assignment] = true
        node.allUp[arc.assignment] = true
        #update the length to terminal
        if isBetterSolutionValue(model, arc.value + arc.destination.lengthToTerminal, node.lengthToTerminal)
            node.lengthToTerminal = arc.value + arc.destination.lengthToTerminal
        end
    end

    #union domains with children
    for child in getChildren(node)
        node.someUp = node.someUp .| child.someUp
        node.allUp = node.allUp .| child.allUp
    end
    #check each all up to make sure its valid
    for i in 1:length(node.allUp)
        if node.allUp[i]
            for arc in node.arcsOut
                if arc.assignment != i
                    if !arc.destination.allUp[i]
                        node.allUp[i]=false
                        break
                    end
                end
            end
        end
    end

end

"""
    createWidthOneDDFromNode(model::SequencingModel; rootNode::Union{RelaxedStateSequencingNode,Nothing}=nothing)

Create a width one decision diagram for a relaxed DD from the sequencing model

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
- `rootNode::Union{RelaxedStateSequencingNode,Nothing}`: optional root node, default is creating one from the model
"""
function createWidthOneDDFromNode(model::SequencingModel; rootNode::Union{RelaxedStateSequencingNode,Nothing}=nothing,masterDomain::BitVector=BitVector())
    if isnothing(rootNode)
        rootNode, masterDomain = createRootRelaxedStateSequencingNodeFromModel(model::SequencingModel)
    end

    relaxedDD = Vector{Vector{RelaxedStateSequencingNode}}()
    push!(relaxedDD, [rootNode])

    #for each layer going down
    for layer in 1:length(model)-length(rootNode)
        #new node starts as a copy of the old node
        currentNode = relaxedDD[layer][1]
        newNode = deepCopy(currentNode)
        push!(relaxedDD,[newNode])
        #update length
        newNode.length = currentNode.length+1
        #reset arcs
        newNode.arcsIn = Vector{SimpleArc}()
        newNode.arcsOut = Vector{SimpleArc}()
        #iterate over domain
        for i in 1:length(masterDomain)
            if masterDomain[i]
                #add arcs connecting the two
                newArc = SimpleArc(currentNode, newNode, i, 0)
                tempVal = getArcValue(model, newArc)
                if !isnothing(tempVal)
                    newArc.value = tempVal
                    push!(currentNode.arcsOut, newArc)
                    push!(newNode.arcsIn, newArc)
                end
            end
        end
        updateNodeDownVariables!(newNode,model)
        updateOutgoingArcValues!(model, newNode)

        if model.precedenceConstraints
            for i in 1:length(masterDomain)
                if !masterDomain[i] && !newNode.allDown[i]
                    if haskey(model.precedenceF2P,i)
                        addToDomain = true
                        for j in model.precedenceF2P[i]
                            if !newNode.someDown[j]
                                addToDomain = false
                                break
                            end
                        end
                        if addToDomain
                            masterDomain[i] = true
                        end
                    else
                        masterDomain[i] = true
                    end
                end
            end
        end

        #check for penultimate node
        if hasEnd(model) && layer == length(model)-length(rootNode)-1
            masterDomain = falses(length(masterDomain))
            masterDomain[model.seqEnd] = true
        end
    end

    if model.precedenceConstraints
        for layer in relaxedDD
            for arc in layer[1].arcsOut
                #handle precedence constraints
                if haskey(model.precedenceF2P,arc.assignment)
                    for i in model.precedenceF2P[arc.assignment]
                        #check each of the relevant precedence constraints
                        #if satisfied, put in the domain
                        if !layer[1].someDown[i]
                            deleteArc!(arc)
                            break
                        end
                    end
                end
            end
        end
    end

    #for each layer going up
    for layer in Iterators.reverse(1:length(relaxedDD))
        updateNodeUpVariables!(relaxedDD[layer][1],model)
    end

    return relaxedDD
end


text:
module TaylorTests

using ForwardDiff
using Morbit
using Test

using LinearAlgebra: norm

const RFD = Morbit.RFD
const Trees = RFD.Trees

ensure_vec( x ) = vec(x)
ensure_vec( x :: Number ) = [x,]

const test_funcs = [
	x -> [x[1],],
	x -> [π + ℯ * sum(x.^2),],
	x -> x.^2,
	x -> [ exp(sum(x)), x[1] + x[end]^3 ]
]

const test_stamps64 = [
	RFD.CFDStamp(1,2),
	RFD.CFDStamp(1,4),
	RFD.CFDStamp(1,6),
	RFD.FFDStamp(1,1),
	RFD.FFDStamp(1,2),
	RFD.FFDStamp(1,3),
	RFD.BFDStamp(1,1),
	RFD.BFDStamp(1,2),
	RFD.BFDStamp(1,3),
]

h32 = RFD.stepsize(RFD.stepsize(Float32))
const test_stamps32 = [
	RFD.CFDStamp(1,2, h32),
	RFD.CFDStamp(1,4, h32),
	RFD.CFDStamp(1,6, h32),
	RFD.FFDStamp(1,1, h32),
	RFD.FFDStamp(1,2, h32),
	RFD.FFDStamp(1,3, h32),
	RFD.BFDStamp(1,1, h32),
	RFD.BFDStamp(1,2, h32),
	RFD.BFDStamp(1,3, h32),
]

h16 = RFD.stepsize(RFD.stepsize(Float16))
const test_stamps16 = [
	RFD.CFDStamp(1,2, h16),
	RFD.CFDStamp(1,4, h16),
	RFD.CFDStamp(1,6, h16),
	RFD.FFDStamp(1,1, h16),
	RFD.FFDStamp(1,2, h16),
	RFD.FFDStamp(1,3, h16),
	RFD.BFDStamp(1,1, h16),
	RFD.BFDStamp(1,2, h16),
	RFD.BFDStamp(1,3, h16),
]

@testset "Recursive Finite Difference Trees" begin
	for ET in [Float16, Float32, Float64]
		stamps = ET == Float16 ? test_stamps16 : ET == Float32 ? test_stamps32 : test_stamps64

		for (j,func) in enumerate(test_funcs)
			n_in = rand(1:5)
			x0 = ensure_vec(rand(ET, n_in))
			fx0 = func(x0)
			for stamp in stamps
				for order in [1,2]
					# first, use convenience function to fill tree
					dw_1 = RFD.DiffWrapper(; x0, fx0, stamp, order)
					RFD.prepare_tree!(dw_1, func)

					# secondly, use the Morbit routine
					dw_2 = RFD.DiffWrapper(; x0, fx0, stamp, order)
					RFD.substitute_leaves!(dw_2)
					all_sites = RFD.collect_leave_sites( dw_2 )
					all_vals = func.(all_sites)
					RFD.set_leave_values!(dw_2, all_vals)

					if order > 1
						for i = 1 : length(fx0)
							Hi = RFD.hessian(dw_1; output_index = i)
							@test Hi ≈ RFD.hessian( dw_2; output_index = i)
							# I don't test against automatic differencing because
							# finite diff hessians become somewhat instable without careful tuning
							# of stepsize …
							# H_ad = ForwardDiff.hessian( x -> func(x)[i], x0 )
							# ET == Float16 || @test norm( Hi .- H_ad ) <= .1
						end
					end

					J1 = RFD.jacobian(dw_1)

					for i = 1 : length(fx0)
						Gi = RFD.gradient( dw_1; output_index = i )
						@test Gi ≈ RFD.gradient( dw_2; output_index = i)
						@test Gi ≈ J1[i,:]
					end

					# does it roughly compare to a jacobian determined by automatic differencing
					J_ad = ForwardDiff.jacobian( func, x0 )
					ET == Float16 || @test norm( J1 .- J_ad ) <= 0.1

				end
			end
		end
	end

end#testset

@testset "Accurate Linear Models, Unconstrained " begin
	mop = Morbit.MixedMOP(2)

	# The gradients for a linear objective will be exact and a linear model
	# should then equal the linear objective globally
	Morbit.add_objective!( mop, x -> sum(x), Morbit.TaylorConfig(;degree=1) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
	Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ [1, 1]
	@test Morbit.get_jacobian(sc, x0) ≈ [ 1 1 ]

end#testset

@testset "Accurate Linear Models, Constrained " begin
	mop = Morbit.MixedMOP( fill(-10.0, 2), fill(10.0, 2) )

	Morbit.add_objective!( mop, x -> sum(x), Morbit.TaylorConfig(;degree=1) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
	#Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	x̂0 = Morbit.scale( x0, smop)

	@test Morbit.eval_models( mod, x̂0 ) ≈ Morbit.eval_objf(objf, x0 )
	@test Morbit.get_gradient( mod, x̂0, 1 ) ≈ ForwardDiff.gradient( ξ -> Morbit.eval_objf(objf, Morbit.unscale(ξ, smop) )[end], x̂0 )
	@test Morbit.get_jacobian( mod, x̂0 ) ≈ ForwardDiff.jacobian( ξ -> Morbit.eval_all_objectives(smop, ξ), x̂0 )

end#testset

@testset "Quadratic Models, Unconstrained " begin
	mop = Morbit.MixedMOP(2)

	Morbit.add_objective!( mop, x -> sum(x.^2), Morbit.TaylorConfig(;degree=2) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
	#Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ 2 .* x0
	@test Morbit.get_jacobian(sc, x0) ≈ 2 .* x0'
end#testset

@testset "Two Parabolas, Taylor Models, 3D2D" begin
	mop = MixedMOP(3)

	add_objective!( mop, x -> sum( ( x .- 1 ).^2 ), Morbit.TaylorApproximateConfig(;degree=2,mode=:fdm) )
	add_objective!( mop, x -> sum( ( x .+ 1 ).^2 ), Morbit.TaylorApproximateConfig(;degree=2,mode=:autodiff) )

	x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])

	@test isapprox(x_fin[1],x_fin[2], atol = 1e-3) && isapprox(x_fin[2], x_fin[3]; atol = 1e-3)

	mop = MixedMOP(3)

	add_objective!( mop, x -> sum( ( x .- 1 ).^2 ),
		TaylorCallbackConfig(; degree = 1, gradients = [ x -> 2 .* ( x .- 1 ), ] ))
	add_objective!( mop, x -> sum( ( x .+ 1 ).^2 ),
		TaylorCallbackConfig(; degree=2, gradients = [ x -> 2 .* ( x .+ 1 ), ],
		hessians = [ x -> [ 2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 2.0] ] ) )

	x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])

	@test isapprox(x_fin[1],x_fin[2], atol = 1e-3) && isapprox(x_fin[2], x_fin[3]; atol = 1e-3)

end#testset

@testset "Two Parabolas, Derivative Accuracy" begin
	# deg 1 autodiff
	mop = Morbit.MixedMOP( 2 )

	Morbit.add_objective!( mop, x -> sum(x), Morbit.TaylorApproximateConfig(;degree=1, mode=:autodiff) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
#	Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	x̂0 = Morbit.scale( x0, smop)

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ [1, 1]
	@test Morbit.get_jacobian(sc, x0) ≈ [ 1 1 ]

	#===================================================================#
	# deg 2 autodiff
	mop = Morbit.MixedMOP( 2 )

	Morbit.add_objective!( mop, x -> sum(x.^2), Morbit.TaylorApproximateConfig(;degree=2, mode=:autodiff) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
	# Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	x̂0 = Morbit.scale( x0, smop)

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ 2 .* x0
	@test Morbit.get_jacobian(sc, x0) ≈ 2 .* x0'
	@test mod.H[1] ≈ [ 2.0 0.0; 0.0 2.0 ]	# hessian of Taylor Model

	#===================================================================#
	# deg 1 finite diff
	mop = Morbit.MixedMOP( 2 )

	Morbit.add_objective!( mop, x -> sum(x), Morbit.TaylorApproximateConfig(;degree=1, mode=:fdm) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
	#Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	x̂0 = Morbit.scale( x0, smop)

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ [1, 1]
	@test Morbit.get_jacobian(sc, x0) ≈ [ 1 1 ]

	#===================================================================#
	# deg 2 autodiff
	mop = Morbit.MixedMOP( 2 )

	Morbit.add_objective!( mop, x -> sum(x.^2), Morbit.TaylorApproximateConfig(;degree=2, mode=:fdm) )

	x0 = [π, -ℯ]

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )
#	Morbit.update_surrogates!( sc, mop, iter_data, data_base, ac )

	mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]

	x̂0 = Morbit.scale( x0, smop)

	@test Morbit.eval_models( mod, x0 ) ≈ Morbit.eval_objf(objf,x0)
	@test Morbit.get_gradient( mod, x0, 1 ) ≈ 2 .* x0
	@test Morbit.get_jacobian(sc, x0) ≈ 2 .* x0'
	@test norm(mod.H[1] .- [ 2.0 0.0; 0.0 2.0 ]) <= 0.1	# hessian of Taylor Model

end

end#module

using .TaylorTests

ext: jl


text:
"Exception type for error thrown by the NetCDF library"
mutable struct NetCDFError <: Exception
    code::Cint
    msg::String
end

"Construct a NetCDFError from the error code"
NetCDFError(code::Cint) = NetCDFError(code, unsafe_string(nc_strerror(code)))

function Base.showerror(io::IO, err::NetCDFError)
    print(io, "NetCDF error code $(err.code):\n\t$(err.msg)")
end

"Check the NetCDF error code, raising an error if nonzero"
function check(code::Cint)
    # zero means success, return
    if code == Cint(0)
        return
    # otherwise throw an error message
    else
        throw(NetCDFError(code))
    end
end

const funext = [
    (Float64, "double", "float64a"),
    (Float32, "float", "float32a"),
    (Int32, "int", "int32a"),
    (UInt32, "uint", "uint32a"),
    (UInt8, "uchar", "uint8a"),
    (Int8, "schar", "int8a"),
    (Int16, "short", "int16a"),
    (UInt16, "ushort", "uint16a"),
    (Int64, "longlong", "int64a"),
    (UInt64, "ulonglong", "uint64a"),
]

const ida = zeros(Int32, 1)
const namea = zeros(UInt8, NC_MAX_NAME + 1)
const lengtha = zeros(Csize_t, 1)
const dimida = zeros(Int32, NC_MAX_VAR_DIMS)
const chunk_sizea = zeros(Csize_t, NC_MAX_VAR_DIMS)
const storagep = zeros(Int32, 1)
const ndima = zeros(Int32, 1)
const nvara = zeros(Int32, 1)
const ngatta = zeros(Int32, 1)
const nunlimdimida = zeros(Int32, 1)
const typea = zeros(Int32, 1)
const natta = zeros(Int32, 1)
const nvals = zeros(Csize_t, 1)
const dima = zeros(Int32, 1)
const varida = zeros(Int32, 1)
const vara = zeros(Int32, 1)
const dumids = zeros(Int32, NC_MAX_DIMS)
const gstart = zeros(UInt, NC_MAX_DIMS)
const gcount = zeros(UInt, NC_MAX_DIMS)

for (t, ending, aname) in funext
    @eval const $(Symbol(aname)) = zeros($t, 1)
end

"""
    nc_char2string(x::Matrix{ASCIIChar})

Convert a `ASCIIChar` array read from a NetCDF variable of type `NC_CHAR` to a Julia array of Strings
"""
function nc_char2string(x::Array{ASCIIChar})
    mapslices(x, dims = 1) do r
        i = findfirst(isequal(ASCIIChar(0)), r)
        if i === nothing
            String(r)
        else
            String(r[1:i-1])
        end
    end |> vec
end

"""
    nc_string2char(x::Array{String})

Convert a Julia `String` array to an `Array{ASCIIChar}` so that it can be written to a
NetCDF variable of type `NC_CHAR`.
"""
function nc_string2char(s::Array{String})
    maxlen = maximum(length, s)
    c = fill(ASCIIChar(0), maxlen + 1, size(s)...)
    for ii in eachindex(s)
        for j = 1:length(s[ii])
            cp = s[ii][j]
            !isascii(cp) && error("Only ascii string can be written to NC_CHAR arrays.")
            cp2 = ASCIIChar(UInt8(cp))
            c[j, ii] = cp2
        end
    end
    c
end

function nc_open(fname::AbstractString, omode::UInt16)
    # Function to open file fname, returns a NetCDF file ID
    # the library doesn't handle this case
    isempty(fname) && throw(ArgumentError("Cannot open empty filename."))
    nc_open(fname, omode, ida)
    return ida[1]
end

function nc_create(name, mode)
    # the library doesn't handle this case
    isempty(name) && throw(ArgumentError("Cannot create empty filename."))
    nc_create(name, mode, ida)
    return ida[1]
end

function nc_inq_dim(id::Integer, idim::Integer)
    # File to inquire dimension idimm returns dimension name and length
    nc_inq_dim(id, idim, namea, lengtha)
    namea[end] = 0
    name = unsafe_string(pointer(namea))
    return (name, lengtha[1])
end

function nc_inq_dimid(id::Integer, name::AbstractString)
    # Function to read dimension id for a given function name
    NetCDF.nc_inq_dimid(id, name, dimida)
    return dimida[1]
end

function nc_redef(nc::NcFile)
    if (!nc.in_def_mode)
        nc_redef(nc.ncid)
        nc.in_def_mode = true
    end
end

function nc_enddef(nc::NcFile)
    if nc.in_def_mode
        nc_enddef(nc.ncid)
        nc.in_def_mode = false
    end
end

function nc_inq(id::Integer)
    # Inquire NetCDF file, return number of dims, number of variables, number of global attributes and number of unlimited dimensions
    nc_inq(id, ndima, nvara, ngatta, nunlimdimida)
    return (ndima[1], nvara[1], ngatta[1], nunlimdimida[1])
end

function nc_inq_attname(ncid::Integer, varid::Integer, attnum::Integer)
    # Get attribute name from attribute number
    nc_inq_attname(ncid, varid, attnum, namea)
    namea[end] = 0
    name = unsafe_string(pointer(namea))
    return name
end

function nc_get_att(ncid::Integer, varid::Integer, att::Union{Integer,AbstractString})
    #Reads attribute name, type and number of values
    name = typeof(att) <: Integer ? nc_inq_attname(ncid, varid, att) : string(att)
    nc_inq_att(ncid, varid, name, typea, nvals)
    text = nc_get_att(ncid, varid, string(name), typea[1], nvals[1])
    return name, text
end

#Define methods for writing attributes
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Int8}) =
    nc_put_att_schar(ncid, varid, name, NC_BYTE, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{UInt8}) =
    nc_put_att_ubyte(ncid, varid, name, NC_UBYTE, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Int16}) =
    nc_put_att_short(ncid, varid, name, NC_SHORT, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{UInt16}) =
    nc_put_att_ushort(ncid, varid, name, NC_USHORT, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Int32}) =
    nc_put_att_int(ncid, varid, name, NC_INT, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{UInt32}) =
    nc_put_att_uint(ncid, varid, name, NC_UINT, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Int64}) =
    nc_put_att_longlong(ncid, varid, name, NC_INT64, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{UInt64}) =
    nc_put_att_ulonglong(ncid, varid, name, NC_UINT64, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Float32}) =
    nc_put_att_float(ncid, varid, name, NC_FLOAT, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Float64}) =
    nc_put_att_double(ncid, varid, name, NC_DOUBLE, length(val), val)
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Array{Any}) =
    error("Writing attribute array of type Any is not possible")

nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Int8) = begin
    int8a[1] = val
    nc_put_att(ncid, varid, name, int8a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::UInt8) = begin
    uint8a[1] = val
    nc_put_att(ncid, varid, name, uint8a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Int16) = begin
    int16a[1] = val
    nc_put_att(ncid, varid, name, int16a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::UInt16) = begin
    uint16a[1] = val
    nc_put_att(ncid, varid, name, uint16a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Int32) = begin
    int32a[1] = val
    nc_put_att(ncid, varid, name, int32a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::UInt32) = begin
    uint32a[1] = val
    nc_put_att(ncid, varid, name, uint32a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Int64) = begin
    int64a[1] = val
    nc_put_att(ncid, varid, name, int64a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::UInt64) = begin
    uint64a[1] = val
    nc_put_att(ncid, varid, name, uint64a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Float32) = begin
    float32a[1] = val
    nc_put_att(ncid, varid, name, float32a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val::Float64) = begin
    float64a[1] = val
    nc_put_att(ncid, varid, name, float64a)
end
nc_put_att(ncid::Integer, varid::Integer, name::AbstractString, val) =
    error("Writing attribute of type $(typeof(val)) is currently not possible.")

function nc_put_att(
    ncid::Integer,
    varid::Integer,
    name::AbstractString,
    val::AbstractString,
)
    val = string(val)
    len = sizeof(val)
    nc_put_att_text(ncid, varid, name, len, val)
end

function nc_put_att(
    ncid::Integer,
    varid::Integer,
    name::AbstractString,
    val::AbstractArray{String},
)
    vals_p = map(x -> pointer(x), val)
    nc_put_att_string(ncid, varid, name, length(val), vals_p)
end

function nc_get_att(
    ncid::Integer,
    varid::Integer,
    name::AbstractString,
    attype::Integer,
    attlen::Integer,
)
    if attype == NC_CHAR
        valsa = zeros(ASCIIChar, attlen + 1)
        nc_get_att_text(ncid, varid, name, valsa)
        valsa[end] = ASCIIChar(0)
        return String(valsa[1:findfirst(isequal(ASCIIChar(0)), valsa)-1])
    else
        valsa = Array{nctype2jltype[attype]}(undef, attlen)
        return nc_get_att!(ncid, varid, name, valsa)
    end
end

nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{UInt8}) =
    begin
        nc_get_att_ubyte(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Int8}) = begin
    nc_get_att_schar(ncid, varid, name, valsa)
    valsa
end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Int16}) =
    begin
        nc_get_att_short(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{UInt16}) =
    begin
        nc_get_att_ushort(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Int32}) =
    begin
        nc_get_att_int(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{UInt32}) =
    begin
        nc_get_att_uint(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Int64}) =
    begin
        nc_get_att_longlong(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{UInt64}) =
    begin
        nc_get_att_ulonglong(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Float32}) =
    begin
        nc_get_att_float(ncid, varid, name, valsa)
        valsa
    end
nc_get_att!(ncid::Integer, varid::Integer, name::AbstractString, valsa::Array{Float64}) =
    begin
        nc_get_att_double(ncid, varid, name, valsa)
        valsa
    end

function nc_get_att!(
    ncid::Integer,
    varid::Integer,
    name::AbstractString,
    valsa::Array{T,N},
) where {T<:AbstractString,N}
  #Assert that length of the attribute matches length of output
    nc_inq_attlen(ncid, varid, name, lengtha)
    @assert lengtha[1] == length(valsa)
    valsa_c = Array{Ptr{UInt8}}(undef, length(valsa))
    nc_get_att_string(ncid, varid, name, valsa_c)
    for i = 1:length(valsa)
        valsa[i] = unsafe_string(valsa_c[i])
    end
    nc_free_string(length(valsa_c), valsa_c)
    valsa
end

function nc_inq_var(nc::NcFile, varid::Integer)
    # Inquire variables in the file
    nc_inq_var(nc.ncid, varid, namea, typea, ndima, dimida, natta)
    dimids = ndima[1] > 0 ? dimida[1:ndima[1]] : Int32[]
    namea[end] = 0
    name = unsafe_string(pointer(namea))
    #Find out chunks
    nc_inq_var_chunking(nc.ncid, varid, storagep, chunk_sizea)
    chunksize = storagep[1] == NC_CONTIGUOUS ? ntuple(i -> 0, ndima[1]) :
                ntuple(i -> chunk_sizea[i], ndima[1])
    return name, typea[1], dimids, natta[1], ndima[1], isdimvar(nc, name), chunksize
end

#Test if a variable name is also a dimension name
isdimvar(v::NcVar) = v.name == v.dim[1].name ? true : false

function isdimvar(nc::NcFile, name::AbstractString)
    for n in nc.dim
        if n[2].name == name
            return true
        end
    end
    return false
end

function getdimnamebyid(nc::NcFile, dimid::Integer)
    da = ""
    for d in nc.dim
        da = d[2].dimid == dimid ? d[2].name : da
    end
    return da
end

function getatts_all(ncid::Integer, varid::Integer, natts::Integer)
    atts = Dict{Any,Any}()
    for attnum = 0:natts-1
        attname, attval = nc_get_att(ncid, varid, attnum)
        if ((length(attval) == 1) && !(typeof(attval) <: AbstractString))
            attval = attval[1]
        end
        atts[attname] = attval
    end
    NC_VERBOSE && println(atts)
    return atts
end

function getatts(
    ncid::Integer,
    varid::Integer,
    attnames::Union{AbstractString,AbstractArray},
)
    atts = Dict{Any,Any}()
    for aattname in (typeof(attnames) <: AbstractArray ? attnames : [attnames])
        attname, attval = nc_get_att(ncid, varid, aattname)
        if ((length(attval) == 1) && !(typeof(attval) <: AbstractString))
            attval = attval[1]
        end
        atts[attname] = attval
    end
    NC_VERBOSE && println(atts)
    return atts
end

function readdimvars(nc::NcFile)
    for v in nc.vars
        if isdimvar(v[2])
            v[2].dim[1].vals = readvar(nc, v[1])
            d[2].atts = v[2].atts
        end
    end
end

function checkboundsNC(v)
    nd = ndims(v)
    for i = 1:nd
        ci = nd + 1 - i
        gstart[ci] < 0 && error("Start index must not be smaller than 1")
        if gstart[ci] + gcount[ci] > v.dim[i].dimlen
            if v.dim[i].unlim
                #Reset length of unlimited dimension
                v.dim[i].dimlen = gstart[ci] + gcount[ci]
            else
                error("Start + Count exceeds dimension length in dimension $(v.dim[i].name)")
            end
        end
        nothing
    end
end

function preparestartcount(start, count, v::NcVar)

    length(start) == v.ndim || error("Length of start ($(length(start))) must equal the number of variable dimensions ($(v.ndim))")
    length(count) == v.ndim || error("Length of start ($(length(count))) must equal the number of variable dimensions ($(v.ndim))")

    p = one(eltype(gcount))
    nd = length(start)

    for i = 1:nd
        ci = nd + 1 - i
        gstart[ci] = start[i] - 1
        gcount[ci] = count[i] < 0 ? v.dim[i].dimlen - gstart[ci] : count[i]

        p *= gcount[ci]
    end

    checkboundsNC(v)

    return p
end

function _readdimvars(nc::NcFile)
    for d in nc.dim
        for v in nc.vars
            if d[2].name == v[2].name
                NC_VERBOSE && println(d[2].name, " ", v[2].name)
                d[2].vals = readvar(nc, v[2].name)
                d[2].atts = v[2].atts
            end
        end
    end
end

ischunked(v::NcVar) = v.chunksize[1] > 0

defaultstart(v::NcVar) = ones(Int, v.ndim)
defaultcount(v::NcVar) = Int[i for i in size(v)]

function parsedimargs(dim)
    isempty(dim) && return NcDim[]
    idim = 0
    dimlen = nothing
    dimvals = nothing
    dimatts = nothing
    name = nothing
    d = NcDim[]
    for a in dim
        NC_VERBOSE && println(a, idim)
        if isa(a, AbstractString)
            #Assume a name is given
            #first create an NcDim object from the last dim
            if (name != nothing)
                push!(d, finalizedim(dimlen, dimvals, dimatts, name))
                idim = idim + 1
                dimlen = nothing
                dimvals = nothing
                dimatts = nothing
                name = nothing
            end
            name = a
        elseif isa(a, Integer)
            #Assume a dimension length is given
            if a > 0
                dimlen = a
            else
                #Assume unlimited dimension
                dimlen = 0
            end
        elseif isa(a, AbstractFloat) && isinf(a)
            #Generate unlimited dimension
            dimlen = 0
        elseif isa(a, AbstractArray)
            #Assume dimension values are given
            if dimvals == nothing
                dimvals = a
                dimlen = length(dimvals)
            else
                error("Dimension values of $name defined more than once")
            end
        elseif isa(a, Dict)
            #Assume attributes are given
            dimatts = dimatts == nothing ? a :
                      error("Dimension attributes of $name defined more than once")
        else
            error("Could not parse argument $a in nccreate")
        end
    end
    push!(d, finalizedim(dimlen, dimvals, dimatts, name))
    return d
end

function finalizedim(dimlen, dimvals, dimatts, name)
    if (dimlen == nothing) && (dimvals == nothing)
        dimlen = 1
    end
    if (dimlen != nothing) && (dimvals == nothing)
        dimvals = Array{Float64}(undef, 0)
    end
    if dimatts == nothing
        dimatts = Dict("missval" => -9999)
    end
    isunlimited = dimlen == 0 ? true : false
    return NcDim(name, dimlen, atts = dimatts, values = dimvals, unlimited = isunlimited)
end


text:
const empty_nb = "{\n \"cells\": [\n  {\n   \"outputs\": [],\n   \"cell_type\": \"markdown\",\n   \"source\": [\n    \"---\\n\",\n    \"\\n\",\n    \"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*\"\n   ],\n   \"metadata\": {}\n  }\n ],\n \"nbformat_minor\": 3,\n \"metadata\": {\n  \"language_info\": {\n   \"file_extension\": \".jl\",\n   \"mimetype\": \"application/julia\",\n   \"name\": \"julia\",\n   \"version\": \"1.3.1\"\n  },\n  \"kernelspec\": {\n   \"name\": \"julia-1.3\",\n   \"display_name\": \"Julia 1.3.1\",\n   \"language\": \"julia\"\n  }\n },\n \"nbformat\": 4\n}\n"
const empty_script = "# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl\n\n"

function rm_if_empty(filepath::String)
    s = open(filepath) do file
        read(file, String)
    end
    if s == empty_nb || s == empty_script
        @info "Removing empty file: $filepath"
        rm(filepath)
    end
end

function postprocess_notebook(nb)
    txt = read(nb, String)
    open(nb, "w") do f
        write(
            f,
            replace(
                txt,
                r"\"outputs\":\ \[\]\,\n\ \ \ \"cell_type\":\ \"markdown\"" =>
                    "\"cell_type\": \"markdown\"",
            ),
        )
    end
end

"""
`literate_file(folder::AbstractString, file::AbstractString)`

Checks if the file has been modified since the last Weave and updates the notebook and tests accordingly.

# Arguments
* `folder`: Name of the folder the tutorial is in
* `file`: Name of the tutorial file
# Key word arguments
* `force`: force literate irrespective of file changes
* `notebook_target_dir`: folder to write notebook
* `test::Bool`: build test file
"""
function literate_file(folder, file; force = false, kwargs...)
    filename = splitext(file)[1]
    srcpath = joinpath(SCRIPT_DIR, folder, file)
    testpath = joinpath(TEST_DIR, folder)
    notebook_path = get(kwargs, :notebook_target_dir, NB_DIR)
    notebookpath = joinpath(notebook_path, folder)
    notebookfilepath = joinpath(notebookpath, join([filename, ".ipynb"]))
    configpath = joinpath(SCRIPT_DIR, folder, filename * "_config.json")

    config = get(kwargs, :config, Dict())
    if isfile(configpath)
        @info "found config file for $filename"
        config = read_json(configpath)
    end

    literate = get(config, "literate", true)

    if literate
        make_test = get(kwargs, :test, get(config, "test", true))
        make_notebook = get(config, "notebook", true)

        if make_test &&
           (mtime(srcpath) > mtime(testpath) || mtime(testpath) == 0.0 || force)
            @warn "Updating tests for $filename."
            fn = Literate.script(srcpath, testpath; config = config, kwargs...)
            rm_if_empty(fn)
        else
            @warn "Skipping tests for $filename."
        end
        if make_notebook && (
            mtime(srcpath) > mtime(notebookfilepath) ||
            mtime(notebookfilepath) == 0.0 ||
            force
        )
            @warn "Converting $filename to Jupyter Notebook."
            fn = Literate.notebook(srcpath, notebookpath; config = config, kwargs...)
            postprocess_notebook(fn)
            rm_if_empty(fn)
        else
            @warn "Skipping Jupyter Notebook for $filename."
        end
    else
        @warn "Skipping literate for $filename per config"
    end
end

"""
`literate_folder(folder::AbstractString)`

Checks the files present in the specified folder for modifications and updates the notebook and tests accordingly.

* `folder` = Name of the folder to check
"""
function literate_folder(folder; force = false, kwargs...)
    for file in readdir(joinpath(SCRIPT_DIR, folder))
        if splitext(file)[end] == ".jl"
            println("")
            println("Building $(joinpath(folder, file))")
            try
                literate_file(folder, file, force = force; kwargs...)
            catch
                @error "failed to build $(joinpath(folder, file))"
            end
            println("")
        end
    end
end

function literate_folder(example::Type{<:Examples}; kwargs...)
    folder = get_dir(example)
    literate_folder(folder; kwargs...)
end

"""
`literate_all()`

Checks every tutorial for modifications and updates the notebook accordingly.
"""
function literate_all(; force = false, kwargs...)
    for folder in readdir(SCRIPT_DIR)
        literate_folder(folder; force = force, kwargs...)
    end
end

"""
`notebook(example::Type{<:Examples}; kwargs...)`

Launches a notebook server for the specified `Examples` folder.

# Arguments
- `example::Type{<:Examples}`: The example category `JuliaExamples`, `PSYExamples`, `PSIExamples`, or `PSDExamples`
- `notebook_target_dir = mktempdir()`: directory to create notebooks and launch server

# Example
`notebook(PSYExamples,  ".")`
"""
function notebook(example::Type{<:Examples}, notebook_target_dir = nothing)
    if isnothing(notebook_target_dir)
        in_pkg_path = startswith(pathof(SIIPExamples), pwd())
        notebook_target_dir = in_pkg_path ? NB_DIR : mktempdir()
    end
    literate_folder(
        example;
        execute = false,
        test = false,
        notebook_target_dir = notebook_target_dir,
    )
    IJulia.notebook(dir = joinpath(notebook_target_dir, get_dir(example)))
end


text:
# test/runtests.jl
using EE, Unitful, PyPlot, Compat.Test

function tests()
  @testset "Subset of tests" begin
    @test ∥(2,3) ≈ 1.2 atol=1E-6
    @test abs(pol(1,45°)) ≈ 1 atol=1E-6
    @test angle(pol(1,45°)) ≈ π/4 atol=1E-6
    @test angle(∠(45°)) ≈ π/4 atol=1E-6
  end
end

tests()